diff options
Diffstat (limited to 'vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib')
7 files changed, 1782 insertions, 0 deletions
diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.bundle b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.bundle Binary files differnew file mode 100755 index 0000000..5a56d47 --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.bundle diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.rb new file mode 100644 index 0000000..5b492ec --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal.rb @@ -0,0 +1,404 @@ +if RUBY_ENGINE == 'jruby' + JRuby::Util.load_ext("org.jruby.ext.bigdecimal.BigDecimalLibrary") + + class BigDecimal + def _decimal_shift(i) # :nodoc: + to_java.move_point_right(i).to_d + end + end +else + require 'bigdecimal.so' +end + +class BigDecimal + module Internal # :nodoc: + # Default extra precision for intermediate calculations + # This value is currently the same as BigDecimal.double_fig, but defined separately for future changes. + EXTRA_PREC = 16 + + # Coerce x to BigDecimal with the specified precision. + # TODO: some methods (example: BigMath.exp) require more precision than specified to coerce. + def self.coerce_to_bigdecimal(x, prec, method_name) # :nodoc: + case x + when BigDecimal + return x + when Integer, Float + return BigDecimal(x, 0) + when Rational + return BigDecimal(x, [prec, 2 * BigDecimal.double_fig].max) + end + raise ArgumentError, "#{x.inspect} can't be coerced into BigDecimal" + end + + def self.coerce_validate_prec(prec, method_name, accept_zero: false) # :nodoc: + unless Integer === prec + original = prec + # Emulate Integer.try_convert for ruby < 3.1 + if prec.respond_to?(:to_int) + prec = prec.to_int + else + raise TypeError, "no implicit conversion of #{original.class} into Integer" + end + raise TypeError, "can't convert #{original.class} to Integer" unless Integer === prec + end + + if accept_zero + raise ArgumentError, "Negative precision for #{method_name}" if prec < 0 + else + raise ArgumentError, "Zero or negative precision for #{method_name}" if prec <= 0 + end + prec + end + + def self.infinity_computation_result # :nodoc: + if BigDecimal.mode(BigDecimal::EXCEPTION_ALL).anybits?(BigDecimal::EXCEPTION_INFINITY) + raise FloatDomainError, "Computation results in 'Infinity'" + end + BigDecimal::INFINITY + end + + def self.underflow_computation_result # :nodoc: + if BigDecimal.mode(BigDecimal::EXCEPTION_ALL).anybits?(BigDecimal::EXCEPTION_UNDERFLOW) + raise FloatDomainError, 'Exponent underflow' + end + BigDecimal(0) + end + + def self.nan_computation_result # :nodoc: + if BigDecimal.mode(BigDecimal::EXCEPTION_ALL).anybits?(BigDecimal::EXCEPTION_NaN) + raise FloatDomainError, "Computation results to 'NaN'" + end + BigDecimal::NAN + end + + # Iteration for Newton's method with increasing precision + def self.newton_loop(prec, initial_precision: BigDecimal.double_fig / 2, safe_margin: 2) # :nodoc: + precs = [] + while prec > initial_precision + precs << prec + prec = (precs.last + 1) / 2 + safe_margin + end + precs.reverse_each do |p| + yield p + end + end + + # Calculates Math.log(x.to_f) considering large or small exponent + def self.float_log(x) # :nodoc: + Math.log(x._decimal_shift(-x.exponent).to_f) + x.exponent * Math.log(10) + end + + # Calculating Taylor series sum using binary splitting method + # Calculates f(x) = (x/d0)*(1+(x/d1)*(1+(x/d2)*(1+(x/d3)*(1+...)))) + # x.n_significant_digits or ds.size must be small to be performant. + def self.taylor_sum_binary_splitting(x, ds, prec) # :nodoc: + fs = ds.map {|d| [0, BigDecimal(d)] } + # fs = [[a0, a1], [b0, b1], [c0, c1], ...] + # f(x) = a0/a1+(x/a1)*(1+b0/b1+(x/b1)*(1+c0/c1+(x/c1)*(1+d0/d1+(x/d1)*(1+...)))) + while fs.size > 1 + # Merge two adjacent fractions + # from: (1 + a0/a1 + x/a1 * (1 + b0/b1 + x/b1 * rest)) + # to: (1 + (a0*b1+x*(b0+b1))/(a1*b1) + (x*x)/(a1*b1) * rest) + xn = xn ? xn.mult(xn, prec) : x + fs = fs.each_slice(2).map do |(a, b)| + b ||= [0, BigDecimal(1)._decimal_shift([xn.exponent, 0].max + 2)] + [ + (a[0] * b[1]).add(xn * (b[0] + b[1]), prec), + a[1].mult(b[1], prec) + ] + end + end + BigDecimal(fs[0][0]).div(fs[0][1], prec) + end + end + + # call-seq: + # self ** other -> bigdecimal + # + # Returns the \BigDecimal value of +self+ raised to power +other+: + # + # b = BigDecimal('3.14') + # b ** 2 # => 0.98596e1 + # b ** 2.0 # => 0.98596e1 + # b ** Rational(2, 1) # => 0.98596e1 + # + # Related: BigDecimal#power. + # + def **(y) + case y + when BigDecimal, Integer, Float, Rational + power(y) + when nil + raise TypeError, 'wrong argument type NilClass' + else + x, y = y.coerce(self) + x**y + end + end + + # call-seq: + # power(n) + # power(n, prec) + # + # Returns the value raised to the power of n. + # + # Also available as the operator **. + # + def power(y, prec = 0) + prec = Internal.coerce_validate_prec(prec, :power, accept_zero: true) + x = self + y = Internal.coerce_to_bigdecimal(y, prec.nonzero? || n_significant_digits, :power) + + return Internal.nan_computation_result if x.nan? || y.nan? + return BigDecimal(1) if y.zero? + + if y.infinite? + if x < 0 + return BigDecimal(0) if x < -1 && y.negative? + return BigDecimal(0) if x > -1 && y.positive? + raise Math::DomainError, 'Result undefined for negative base raised to infinite power' + elsif x < 1 + return y.positive? ? BigDecimal(0) : BigDecimal::Internal.infinity_computation_result + elsif x == 1 + return BigDecimal(1) + else + return y.positive? ? BigDecimal::Internal.infinity_computation_result : BigDecimal(0) + end + end + + if x.infinite? && y < 0 + # Computation result will be +0 or -0. Avoid overflow. + neg = x < 0 && y.frac.zero? && y % 2 == 1 + return neg ? -BigDecimal(0) : BigDecimal(0) + end + + if x.zero? + return BigDecimal(1) if y.zero? + return BigDecimal(0) if y > 0 + if y.frac.zero? && y % 2 == 1 && x.sign == -1 + return -BigDecimal::Internal.infinity_computation_result + else + return BigDecimal::Internal.infinity_computation_result + end + elsif x < 0 + if y.frac.zero? + if y % 2 == 0 + return (-x).power(y, prec) + else + return -(-x).power(y, prec) + end + else + raise Math::DomainError, 'Computation results in complex number' + end + elsif x == 1 + return BigDecimal(1) + end + + limit = BigDecimal.limit + frac_part = y.frac + + if frac_part.zero? && prec.zero? && limit.zero? + # Infinite precision calculation for `x ** int` and `x.power(int)` + int_part = y.fix.to_i + int_part = -int_part if (neg = int_part < 0) + ans = BigDecimal(1) + n = 1 + xn = x + while true + ans *= xn if int_part.allbits?(n) + n <<= 1 + break if n > int_part + xn *= xn + # Detect overflow/underflow before consuming infinite memory + if (xn.exponent.abs - 1) * int_part / n >= 0x7FFFFFFFFFFFFFFF + return ((xn.exponent > 0) ^ neg ? BigDecimal::Internal.infinity_computation_result : BigDecimal(0)) * (int_part.even? || x > 0 ? 1 : -1) + end + end + return neg ? BigDecimal(1) / ans : ans + end + + result_prec = prec.nonzero? || [x.n_significant_digits, y.n_significant_digits, BigDecimal.double_fig].max + BigDecimal.double_fig + result_prec = [result_prec, limit].min if prec.zero? && limit.nonzero? + + prec2 = result_prec + BigDecimal::Internal::EXTRA_PREC + + if y < 0 + inv = x.power(-y, prec2) + return BigDecimal(0) if inv.infinite? + return BigDecimal::Internal.infinity_computation_result if inv.zero? + return BigDecimal(1).div(inv, result_prec) + end + + if frac_part.zero? && y.exponent < Math.log(result_prec) * 5 + 20 + # Use exponentiation by squaring if y is an integer and not too large + pow_prec = prec2 + y.exponent + n = 1 + xn = x + ans = BigDecimal(1) + int_part = y.fix.to_i + while true + ans = ans.mult(xn, pow_prec) if int_part.allbits?(n) + n <<= 1 + break if n > int_part + xn = xn.mult(xn, pow_prec) + end + ans.mult(1, result_prec) + else + if x > 1 && x.finite? + # To calculate exp(z, prec), z needs prec+max(z.exponent, 0) precision if z > 0. + # Estimate (y*log(x)).exponent + logx_exponent = x < 2 ? (x - 1).exponent : Math.log10(x.exponent).round + ylogx_exponent = y.exponent + logx_exponent + prec2 += [ylogx_exponent, 0].max + end + BigMath.exp(BigMath.log(x, prec2).mult(y, prec2), result_prec) + end + end + + # Returns the square root of the value. + # + # Result has at least prec significant digits. + # + def sqrt(prec) + prec = Internal.coerce_validate_prec(prec, :sqrt, accept_zero: true) + return Internal.infinity_computation_result if infinite? == 1 + + raise FloatDomainError, 'sqrt of negative value' if self < 0 + raise FloatDomainError, "sqrt of 'NaN'(Not a Number)" if nan? + return self if zero? + + if prec == 0 + limit = BigDecimal.limit + prec = n_significant_digits + BigDecimal.double_fig + prec = [limit, prec].min if limit.nonzero? + end + + ex = exponent / 2 + x = _decimal_shift(-2 * ex) + y = BigDecimal(Math.sqrt(x.to_f), 0) + Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p| + y = y.add(x.div(y, p), p).div(2, p) + end + y._decimal_shift(ex).mult(1, prec) + end +end + +# Core BigMath methods for BigDecimal (log, exp) are defined here. +# Other methods (sin, cos, atan) are defined in 'bigdecimal/math.rb'. +module BigMath + module_function + + # call-seq: + # BigMath.log(decimal, numeric) -> BigDecimal + # + # Computes the natural logarithm of +decimal+ to the specified number of + # digits of precision, +numeric+. + # + # If +decimal+ is zero or negative, raises Math::DomainError. + # + # If +decimal+ is positive infinity, returns Infinity. + # + # If +decimal+ is NaN, returns NaN. + # + def log(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :log) + raise Math::DomainError, 'Complex argument for BigMath.log' if Complex === x + + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log) + return BigDecimal::Internal.nan_computation_result if x.nan? + raise Math::DomainError, 'Negative argument for log' if x < 0 + return -BigDecimal::Internal.infinity_computation_result if x.zero? + return BigDecimal::Internal.infinity_computation_result if x.infinite? + return BigDecimal(0) if x == 1 + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + + # Reduce x to near 1 + if x > 1.01 || x < 0.99 + # log(x) = log(x/exp(logx_approx)) + logx_approx + logx_approx = BigDecimal(BigDecimal::Internal.float_log(x), 0) + x = x.div(exp(logx_approx, prec2), prec2) + else + logx_approx = BigDecimal(0) + end + + # Solve exp(y) - x = 0 with Newton's method + # Repeat: y -= (exp(y) - x) / exp(y) + y = BigDecimal(BigDecimal::Internal.float_log(x), 0) + exp_additional_prec = [-(x - 1).exponent, 0].max + BigDecimal::Internal.newton_loop(prec2) do |p| + expy = exp(y, p + exp_additional_prec) + y = y.sub(expy.sub(x, p).div(expy, p), p) + end + y.add(logx_approx, prec) + end + + private_class_method def _exp_binary_splitting(x, prec) # :nodoc: + return BigDecimal(1) if x.zero? + # Find k that satisfies x**k / k! < 10**(-prec) + log10 = Math.log(10) + logx = BigDecimal::Internal.float_log(x.abs) + step = (1..).bsearch { |k| Math.lgamma(k + 1)[0] - k * logx > prec * log10 } + # exp(x)-1 = x*(1+x/2*(1+x/3*(1+x/4*(1+x/5*(1+...))))) + 1 + BigDecimal::Internal.taylor_sum_binary_splitting(x, [*1..step], prec) + end + + # call-seq: + # BigMath.exp(decimal, numeric) -> BigDecimal + # + # Computes the value of e (the base of natural logarithms) raised to the + # power of +decimal+, to the specified number of digits of precision. + # + # If +decimal+ is infinity, returns Infinity. + # + # If +decimal+ is NaN, returns NaN. + # + def exp(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :exp) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :exp) + return BigDecimal::Internal.nan_computation_result if x.nan? + if x.infinite? || x.exponent >= 21 # exp(10**20) and exp(-10**20) overflows/underflows 64-bit exponent + if x.positive? + return BigDecimal::Internal.infinity_computation_result + elsif x.infinite? + # exp(-Infinity) is +0 by definition, this is not an underflow. + return BigDecimal(0) + else + return BigDecimal::Internal.underflow_computation_result + end + end + + return BigDecimal(1) if x.zero? + + # exp(x * 10**cnt) = exp(x)**(10**cnt) + cnt = x < -1 || x > 1 ? x.exponent : 0 + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + cnt + x = x._decimal_shift(-cnt) + + # Decimal form of bit-burst algorithm + # Calculate exp(x.xxxxxxxxxxxxxxxx) as + # exp(x.xx) * exp(0.00xx) * exp(0.0000xxxx) * exp(0.00000000xxxxxxxx) + x = x.mult(1, prec2) + n = 2 + y = BigDecimal(1) + BigDecimal.save_limit do + BigDecimal.limit(0) + while x != 0 do + partial_x = x.truncate(n) + x -= partial_x + y = y.mult(_exp_binary_splitting(partial_x, prec2), prec2) + n *= 2 + end + end + + # calculate exp(x * 10**cnt) from exp(x) + # exp(x * 10**k) = exp(x * 10**(k - 1)) ** 10 + cnt.times do + y2 = y.mult(y, prec2) + y5 = y2.mult(y2, prec2).mult(y, prec2) + y = y5.mult(y5, prec2) + end + + y.mult(1, prec) + end +end diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/jacobian.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/jacobian.rb new file mode 100644 index 0000000..9a5e5be --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/jacobian.rb @@ -0,0 +1,92 @@ +# frozen_string_literal: false + +require 'bigdecimal' + +warn "'bigdecimal/jacobian' is deprecated and will be removed in a future release." + +# require 'bigdecimal/jacobian' +# +# Provides methods to compute the Jacobian matrix of a set of equations at a +# point x. In the methods below: +# +# f is an Object which is used to compute the Jacobian matrix of the equations. +# It must provide the following methods: +# +# f.values(x):: returns the values of all functions at x +# +# f.zero:: returns 0.0 +# f.one:: returns 1.0 +# f.two:: returns 2.0 +# f.ten:: returns 10.0 +# +# f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal. +# +# x is the point at which to compute the Jacobian. +# +# fx is f.values(x). +# +module Jacobian + module_function + + # Determines the equality of two numbers by comparing to zero, or using the epsilon value + def isEqual(a,b,zero=0.0,e=1.0e-8) + aa = a.abs + bb = b.abs + if aa == zero && bb == zero then + true + else + if ((a-b)/(aa+bb)).abs < e then + true + else + false + end + end + end + + + # Computes the derivative of +f[i]+ at +x[i]+. + # +fx+ is the value of +f+ at +x+. + def dfdxi(f,fx,x,i) + nRetry = 0 + n = x.size + xSave = x[i] + ok = 0 + ratio = f.ten*f.ten*f.ten + dx = x[i].abs/ratio + dx = fx[i].abs/ratio if isEqual(dx,f.zero,f.zero,f.eps) + dx = f.one/f.ten if isEqual(dx,f.zero,f.zero,f.eps) + until ok>0 do + deriv = [] + nRetry += 1 + if nRetry > 100 + raise "Singular Jacobian matrix. No change at x[" + i.to_s + "]" + end + dx = dx*f.two + x[i] += dx + fxNew = f.values(x) + for j in 0...n do + if !isEqual(fxNew[j],fx[j],f.zero,f.eps) then + ok += 1 + deriv <<= (fxNew[j]-fx[j])/dx + else + deriv <<= f.zero + end + end + x[i] = xSave + end + deriv + end + + # Computes the Jacobian of +f+ at +x+. +fx+ is the value of +f+ at +x+. + def jacobian(f,fx,x) + n = x.size + dfdx = Array.new(n*n) + for i in 0...n do + df = dfdxi(f,fx,x,i) + for j in 0...n do + dfdx[j*n+i] = df[j] + end + end + dfdx + end +end diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/ludcmp.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/ludcmp.rb new file mode 100644 index 0000000..d112fc1 --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/ludcmp.rb @@ -0,0 +1,91 @@ +# frozen_string_literal: false +require 'bigdecimal' + +warn "'bigdecimal/ludcmp' is deprecated and will be removed in a future release." + +# +# Solves a*x = b for x, using LU decomposition. +# +module LUSolve + module_function + + # Performs LU decomposition of the n by n matrix a. + def ludecomp(a,n,zero=0,one=1) + prec = BigDecimal.limit(nil) + ps = [] + scales = [] + for i in 0...n do # pick up largest(abs. val.) element in each row. + ps <<= i + nrmrow = zero + ixn = i*n + for j in 0...n do + biggst = a[ixn+j].abs + nrmrow = biggst if biggst>nrmrow + end + if nrmrow>zero then + scales <<= one.div(nrmrow,prec) + else + raise "Singular matrix" + end + end + n1 = n - 1 + for k in 0...n1 do # Gaussian elimination with partial pivoting. + biggst = zero; + for i in k...n do + size = a[ps[i]*n+k].abs*scales[ps[i]] + if size>biggst then + biggst = size + pividx = i + end + end + raise "Singular matrix" if biggst<=zero + if pividx!=k then + j = ps[k] + ps[k] = ps[pividx] + ps[pividx] = j + end + pivot = a[ps[k]*n+k] + for i in (k+1)...n do + psin = ps[i]*n + a[psin+k] = mult = a[psin+k].div(pivot,prec) + if mult!=zero then + pskn = ps[k]*n + for j in (k+1)...n do + a[psin+j] -= mult.mult(a[pskn+j],prec) + end + end + end + end + raise "Singular matrix" if a[ps[n1]*n+n1] == zero + ps + end + + # Solves a*x = b for x, using LU decomposition. + # + # a is a matrix, b is a constant vector, x is the solution vector. + # + # ps is the pivot, a vector which indicates the permutation of rows performed + # during LU decomposition. + def lusolve(a,b,ps,zero=0.0) + prec = BigDecimal.limit(nil) + n = ps.size + x = [] + for i in 0...n do + dot = zero + psin = ps[i]*n + for j in 0...i do + dot = a[psin+j].mult(x[j],prec) + dot + end + x <<= b[ps[i]] - dot + end + (n-1).downto(0) do |i| + dot = zero + psin = ps[i]*n + for j in (i+1)...n do + dot = a[psin+j].mult(x[j],prec) + dot + end + x[i] = (x[i]-dot).div(a[psin+i],prec) + end + x + end +end diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/math.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/math.rb new file mode 100644 index 0000000..1b3357c --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/math.rb @@ -0,0 +1,927 @@ +# frozen_string_literal: true +require 'bigdecimal' + +# +#-- +# Contents: +# sqrt(x, prec) +# cbrt(x, prec) +# hypot(x, y, prec) +# sin (x, prec) +# cos (x, prec) +# tan (x, prec) +# asin(x, prec) +# acos(x, prec) +# atan(x, prec) +# atan2(y, x, prec) +# sinh (x, prec) +# cosh (x, prec) +# tanh (x, prec) +# asinh(x, prec) +# acosh(x, prec) +# atanh(x, prec) +# log2 (x, prec) +# log10(x, prec) +# log1p(x, prec) +# expm1(x, prec) +# erf (x, prec) +# erfc(x, prec) +# gamma(x, prec) +# lgamma(x, prec) +# frexp(x) +# ldexp(x, exponent) +# PI (prec) +# E (prec) == exp(1.0,prec) +# +# where: +# x, y ... BigDecimal number to be computed. +# prec ... Number of digits to be obtained. +#++ +# +# Provides mathematical functions. +# +# Example: +# +# require "bigdecimal/math" +# +# include BigMath +# +# a = BigDecimal((PI(49)/2).to_s) +# puts sin(a,100) # => 0.9999999999...9999999986e0 +# +module BigMath + module_function + + # call-seq: + # sqrt(decimal, numeric) -> BigDecimal + # + # Computes the square root of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # BigMath.sqrt(BigDecimal('2'), 32).to_s + # #=> "0.14142135623730950488016887242097e1" + # + def sqrt(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :sqrt) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sqrt) + x.sqrt(prec) + end + + + # Returns [sign, reduced_x] where reduced_x is in -pi/2..pi/2 + # and satisfies sin(x) = sign * sin(reduced_x) + # If add_half_pi is true, adds pi/2 to x before reduction. + # Precision of pi is adjusted to ensure reduced_x has the required precision. + private_class_method def _sin_periodic_reduction(x, prec, add_half_pi: false) # :nodoc: + return [1, x] if -Math::PI/2 <= x && x <= Math::PI/2 && !add_half_pi + + mod_prec = prec + BigDecimal::Internal::EXTRA_PREC + pi_extra_prec = [x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC + while true + pi = PI(mod_prec + pi_extra_prec) + half_pi = pi / 2 + div, mod = (add_half_pi ? x + pi : x + half_pi).divmod(pi) + mod -= half_pi + if mod.zero? || mod_prec + mod.exponent <= 0 + # mod is too small to estimate required pi precision + mod_prec = mod_prec * 3 / 2 + BigDecimal::Internal::EXTRA_PREC + elsif mod_prec + mod.exponent < prec + # Estimate required precision of pi + mod_prec = prec - mod.exponent + BigDecimal::Internal::EXTRA_PREC + else + return [div % 2 == 0 ? 1 : -1, mod.mult(1, prec)] + end + end + end + + private_class_method def _sin_binary_splitting(x, prec) # :nodoc: + return x if x.zero? + x2 = x.mult(x, prec) + # Find k that satisfies x2**k / (2k+1)! < 10**(-prec) + log10 = Math.log(10) + logx = BigDecimal::Internal.float_log(x.abs) + step = (1..).bsearch { |k| Math.lgamma(2 * k + 1)[0] - 2 * k * logx > prec * log10 } + # Construct denominator sequence for binary splitting + # sin(x) = x*(1-x2/(2*3)*(1-x2/(4*5)*(1-x2/(6*7)*(1-x2/(8*9)*(1-...))))) + ds = (1..step).map {|i| -(2 * i) * (2 * i + 1) } + x.mult(1 + BigDecimal::Internal.taylor_sum_binary_splitting(x2, ds, prec), prec) + end + + private_class_method def _sin_around_zero(x, prec) # :nodoc: + # Divide x into several parts + # sin(x.xxxxxxxx...) = sin(x.xx + 0.00xx + 0.0000xxxx + ...) + # Calculate sin of each part and restore sin(0.xxxxxxxx...) using addition theorem. + sin = BigDecimal(0) + cos = BigDecimal(1) + n = 2 + while x != 0 do + partial_x = x.truncate(n) + x -= partial_x + s = _sin_binary_splitting(partial_x, prec) + c = (1 - s * s).sqrt(prec) + sin, cos = (sin * c).add(cos * s, prec), (cos * c).sub(sin * s, prec) + n *= 2 + end + sin.clamp(BigDecimal(-1), BigDecimal(1)) + end + + # call-seq: + # cbrt(decimal, numeric) -> BigDecimal + # + # Computes the cube root of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # BigMath.cbrt(BigDecimal('2'), 32).to_s + # #=> "0.12599210498948731647672106072782e1" + # + def cbrt(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :cbrt) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cbrt) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite? + return BigDecimal(0) if x.zero? + + x = -x if neg = x < 0 + ex = x.exponent / 3 + x = x._decimal_shift(-3 * ex) + y = BigDecimal(Math.cbrt(x.to_f), 0) + BigDecimal::Internal.newton_loop(prec + BigDecimal::Internal::EXTRA_PREC) do |p| + y = (2 * y + x.div(y, p).div(y, p)).div(3, p) + end + y._decimal_shift(ex).mult(neg ? -1 : 1, prec) + end + + # call-seq: + # hypot(x, y, numeric) -> BigDecimal + # + # Returns sqrt(x**2 + y**2) to the specified number of digits of + # precision, +numeric+. + # + # BigMath.hypot(BigDecimal('1'), BigDecimal('2'), 32).to_s + # #=> "0.22360679774997896964091736687313e1" + # + def hypot(x, y, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :hypot) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :hypot) + y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :hypot) + return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan? + return BigDecimal::Internal.infinity_computation_result if x.infinite? || y.infinite? + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + sqrt(x.mult(x, prec2) + y.mult(y, prec2), prec) + end + + # call-seq: + # sin(decimal, numeric) -> BigDecimal + # + # Computes the sine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is Infinity or NaN, returns NaN. + # + # BigMath.sin(BigMath.PI(5)/4, 32).to_s + # #=> "0.70710807985947359435812921837984e0" + # + def sin(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :sin) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sin) + return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan? + n = prec + BigDecimal::Internal::EXTRA_PREC + sign, x = _sin_periodic_reduction(x, n) + _sin_around_zero(x, n).mult(sign, prec) + end + + # call-seq: + # cos(decimal, numeric) -> BigDecimal + # + # Computes the cosine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is Infinity or NaN, returns NaN. + # + # BigMath.cos(BigMath.PI(16), 32).to_s + # #=> "-0.99999999999999999999999999999997e0" + # + def cos(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :cos) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cos) + return BigDecimal::Internal.nan_computation_result if x.infinite? || x.nan? + n = prec + BigDecimal::Internal::EXTRA_PREC + sign, x = _sin_periodic_reduction(x, n, add_half_pi: true) + _sin_around_zero(x, n).mult(sign, prec) + end + + # call-seq: + # tan(decimal, numeric) -> BigDecimal + # + # Computes the tangent of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is Infinity or NaN, returns NaN. + # + # BigMath.tan(BigDecimal("0.0"), 4).to_s + # #=> "0.0" + # + # BigMath.tan(BigMath.PI(24) / 4, 32).to_s + # #=> "0.99999999999999999999999830836025e0" + # + def tan(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :tan) + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + sin(x, prec2).div(cos(x, prec2), prec) + end + + # call-seq: + # asin(decimal, numeric) -> BigDecimal + # + # Computes the arcsine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.asin(BigDecimal('0.5'), 32).to_s + # #=> "0.52359877559829887307710723054658e0" + # + def asin(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :asin) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asin) + raise Math::DomainError, "Out of domain argument for asin" if x < -1 || x > 1 + return BigDecimal::Internal.nan_computation_result if x.nan? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + cos = (1 - x**2).sqrt(prec2) + if cos.zero? + PI(prec2).div(x > 0 ? 2 : -2, prec) + else + atan(x.div(cos, prec2), prec) + end + end + + # call-seq: + # acos(decimal, numeric) -> BigDecimal + # + # Computes the arccosine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.acos(BigDecimal('0.5'), 32).to_s + # #=> "0.10471975511965977461542144610932e1" + # + def acos(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :acos) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acos) + raise Math::DomainError, "Out of domain argument for acos" if x < -1 || x > 1 + return BigDecimal::Internal.nan_computation_result if x.nan? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + return (PI(prec2) / 2).sub(asin(x, prec2), prec) if x < 0 + return PI(prec2).div(2, prec) if x.zero? + + sin = (1 - x**2).sqrt(prec2) + atan(sin.div(x, prec2), prec) + end + + # call-seq: + # atan(decimal, numeric) -> BigDecimal + # + # Computes the arctangent of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.atan(BigDecimal('-1'), 32).to_s + # #=> "-0.78539816339744830961566084581988e0" + # + def atan(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan) + return BigDecimal::Internal.nan_computation_result if x.nan? + n = prec + BigDecimal::Internal::EXTRA_PREC + return PI(n).div(2 * x.infinite?, prec) if x.infinite? + + x = -x if neg = x < 0 + x = BigDecimal(1).div(x, n) if inv = x < -1 || x > 1 + + # Solve tan(y) - x = 0 with Newton's method + # Repeat: y -= (tan(y) - x) * cos(y)**2 + y = BigDecimal(Math.atan(x.to_f), 0) + BigDecimal::Internal.newton_loop(n) do |p| + s = sin(y, p) + c = (1 - s * s).sqrt(p) + y = y.sub(c * (s.sub(c * x.mult(1, p), p)), p) + end + y = PI(n) / 2 - y if inv + y.mult(neg ? -1 : 1, prec) + end + + # call-seq: + # atan2(decimal, decimal, numeric) -> BigDecimal + # + # Computes the arctangent of y and x to the specified number of digits of + # precision, +numeric+. + # + # BigMath.atan2(BigDecimal('-1'), BigDecimal('1'), 32).to_s + # #=> "-0.78539816339744830961566084581988e0" + # + def atan2(y, x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :atan2) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atan2) + y = BigDecimal::Internal.coerce_to_bigdecimal(y, prec, :atan2) + return BigDecimal::Internal.nan_computation_result if x.nan? || y.nan? + + if x.infinite? || y.infinite? + one = BigDecimal(1) + zero = BigDecimal(0) + x = x.infinite? ? (x > 0 ? one : -one) : zero + y = y.infinite? ? (y > 0 ? one : -one) : y.sign * zero + end + + return x.sign >= 0 ? BigDecimal(0) : y.sign * PI(prec) if y.zero? + + y = -y if neg = y < 0 + xlarge = y.abs < x.abs + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + if x > 0 + v = xlarge ? atan(y.div(x, prec2), prec) : PI(prec2) / 2 - atan(x.div(y, prec2), prec2) + else + v = xlarge ? PI(prec2) - atan(-y.div(x, prec2), prec2) : PI(prec2) / 2 + atan(x.div(-y, prec2), prec2) + end + v.mult(neg ? -1 : 1, prec) + end + + # call-seq: + # sinh(decimal, numeric) -> BigDecimal + # + # Computes the hyperbolic sine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.sinh(BigDecimal('1'), 32).to_s + # #=> "0.11752011936438014568823818505956e1" + # + def sinh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :sinh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :sinh) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + prec2 -= x.exponent if x.exponent < 0 + e = exp(x, prec2) + (e - BigDecimal(1).div(e, prec2)).div(2, prec) + end + + # call-seq: + # cosh(decimal, numeric) -> BigDecimal + # + # Computes the hyperbolic cosine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.cosh(BigDecimal('1'), 32).to_s + # #=> "0.15430806348152437784779056207571e1" + # + def cosh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :cosh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :cosh) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result if x.infinite? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + e = exp(x, prec2) + (e + BigDecimal(1).div(e, prec2)).div(2, prec) + end + + # call-seq: + # tanh(decimal, numeric) -> BigDecimal + # + # Computes the hyperbolic tangent of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.tanh(BigDecimal('1'), 32).to_s + # #=> "0.76159415595576488811945828260479e0" + # + def tanh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :tanh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :tanh) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal(x.infinite?) if x.infinite? + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + [-x.exponent, 0].max + e = exp(x, prec2) + einv = BigDecimal(1).div(e, prec2) + (e - einv).div(e + einv, prec) + end + + # call-seq: + # asinh(decimal, numeric) -> BigDecimal + # + # Computes the inverse hyperbolic sine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.asinh(BigDecimal('1'), 32).to_s + # #=> "0.88137358701954302523260932497979e0" + # + def asinh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :asinh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :asinh) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result * x.infinite? if x.infinite? + return -asinh(-x, prec) if x < 0 + + sqrt_prec = prec + [-x.exponent, 0].max + BigDecimal::Internal::EXTRA_PREC + log(x + sqrt(x**2 + 1, sqrt_prec), prec) + end + + # call-seq: + # acosh(decimal, numeric) -> BigDecimal + # + # Computes the inverse hyperbolic cosine of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.acosh(BigDecimal('2'), 32).to_s + # #=> "0.1316957896924816708625046347308e1" + # + def acosh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :acosh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :acosh) + raise Math::DomainError, "Out of domain argument for acosh" if x < 1 + return BigDecimal::Internal.infinity_computation_result if x.infinite? + return BigDecimal::Internal.nan_computation_result if x.nan? + + log(x + sqrt(x**2 - 1, prec + BigDecimal::Internal::EXTRA_PREC), prec) + end + + # call-seq: + # atanh(decimal, numeric) -> BigDecimal + # + # Computes the inverse hyperbolic tangent of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.atanh(BigDecimal('0.5'), 32).to_s + # #=> "0.54930614433405484569762261846126e0" + # + def atanh(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :atanh) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :atanh) + raise Math::DomainError, "Out of domain argument for atanh" if x < -1 || x > 1 + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result if x == 1 + return -BigDecimal::Internal.infinity_computation_result if x == -1 + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + (log(x + 1, prec2) - log(1 - x, prec2)).div(2, prec) + end + + # call-seq: + # BigMath.log2(decimal, numeric) -> BigDecimal + # + # Computes the base 2 logarithm of +decimal+ to the specified number of + # digits of precision, +numeric+. + # + # If +decimal+ is zero or negative, raises Math::DomainError. + # + # If +decimal+ is positive infinity, returns Infinity. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.log2(BigDecimal('3'), 32).to_s + # #=> "0.15849625007211561814537389439478e1" + # + def log2(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :log2) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log2) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1 + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2 + v = log(x, prec2).div(log(BigDecimal(2), prec2), prec2) + # Perform half-up rounding to calculate log2(2**n)==n correctly in every rounding mode + v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP) + v.mult(1, prec) + end + + # call-seq: + # BigMath.log10(decimal, numeric) -> BigDecimal + # + # Computes the base 10 logarithm of +decimal+ to the specified number of + # digits of precision, +numeric+. + # + # If +decimal+ is zero or negative, raises Math::DomainError. + # + # If +decimal+ is positive infinity, returns Infinity. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.log10(BigDecimal('3'), 32).to_s + # #=> "0.47712125471966243729502790325512e0" + # + def log10(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :log10) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log10) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal::Internal.infinity_computation_result if x.infinite? == 1 + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC * 3 / 2 + v = log(x, prec2).div(log(BigDecimal(10), prec2), prec2) + # Perform half-up rounding to calculate log10(10**n)==n correctly in every rounding mode + v = v.round(prec + BigDecimal::Internal::EXTRA_PREC - (v.exponent < 0 ? v.exponent : 0), BigDecimal::ROUND_HALF_UP) + v.mult(1, prec) + end + + # call-seq: + # BigMath.log1p(decimal, numeric) -> BigDecimal + # + # Computes log(1 + decimal) to the specified number of digits of precision, +numeric+. + # + # BigMath.log1p(BigDecimal('0.1'), 32).to_s + # #=> "0.95310179804324860043952123280765e-1" + # + def log1p(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :log1p) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :log1p) + raise Math::DomainError, 'Out of domain argument for log1p' if x < -1 + + return log(x + 1, prec) + end + + # call-seq: + # BigMath.expm1(decimal, numeric) -> BigDecimal + # + # Computes exp(decimal) - 1 to the specified number of digits of precision, +numeric+. + # + # BigMath.expm1(BigDecimal('0.1'), 32).to_s + # #=> "0.10517091807564762481170782649025e0" + # + def expm1(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :expm1) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :expm1) + return BigDecimal(-1) if x.infinite? == -1 + + exp_prec = prec + if x < -1 + # log10(exp(x)) = x * log10(e) + lg_e = 0.4342944819032518 + exp_prec = prec + (lg_e * x).ceil + BigDecimal::Internal::EXTRA_PREC + elsif x < 1 + exp_prec = prec - x.exponent + BigDecimal::Internal::EXTRA_PREC + else + exp_prec = prec + end + + return BigDecimal(-1) if exp_prec <= 0 + + exp(x, exp_prec).sub(1, prec) + end + + # call-seq: + # erf(decimal, numeric) -> BigDecimal + # + # Computes the error function of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.erf(BigDecimal('1'), 32).to_s + # #=> "0.84270079294971486934122063508261e0" + # + def erf(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :erf) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erf) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal(x.infinite?) if x.infinite? + return BigDecimal(0) if x == 0 + return -erf(-x, prec) if x < 0 + return BigDecimal(1) if x > 5000000000 # erf(5000000000) > 1 - 1e-10000000000000000000 + + if x > 8 + xf = x.to_f + log10_erfc = -xf ** 2 / Math.log(10) - Math.log10(xf * Math::PI ** 0.5) + erfc_prec = [prec + log10_erfc.ceil, 1].max + erfc = _erfc_asymptotic(x, erfc_prec) + return BigDecimal(1).sub(erfc, prec) if erfc + end + + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + x_smallprec = x.mult(1, Integer.sqrt(prec2) / 2) + # Taylor series of x with small precision is fast + erf1 = _erf_taylor(x_smallprec, BigDecimal(0), BigDecimal(0), prec2) + # Taylor series converges quickly for small x + _erf_taylor(x - x_smallprec, x_smallprec, erf1, prec2).mult(1, prec) + end + + # call-seq: + # erfc(decimal, numeric) -> BigDecimal + # + # Computes the complementary error function of +decimal+ to the specified number of digits of + # precision, +numeric+. + # + # If +decimal+ is NaN, returns NaN. + # + # BigMath.erfc(BigDecimal('10'), 32).to_s + # #=> "0.20884875837625447570007862949578e-44" + # + def erfc(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :erfc) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :erfc) + return BigDecimal::Internal.nan_computation_result if x.nan? + return BigDecimal(1 - x.infinite?) if x.infinite? + return BigDecimal(1).sub(erf(x, prec + BigDecimal::Internal::EXTRA_PREC), prec) if x < 0.5 + return BigDecimal::Internal.underflow_computation_result if x > 5000000000 # erfc(5000000000) < 1e-10000000000000000000 (underflow) + + if x >= 8 + y = _erfc_asymptotic(x, prec) + return y.mult(1, prec) if y + end + + # erfc(x) = 1 - erf(x) < exp(-x**2)/x/sqrt(pi) + # Precision of erf(x) needs about log10(exp(-x**2)/x/sqrt(pi)) extra digits + log10 = 2.302585092994046 + xf = x.to_f + high_prec = prec + BigDecimal::Internal::EXTRA_PREC + ((xf**2 + Math.log(xf) + Math.log(Math::PI)/2) / log10).ceil + BigDecimal(1).sub(erf(x, high_prec), prec) + end + + # Calculates erf(x + a) + private_class_method def _erf_taylor(x, a, erf_a, prec) # :nodoc: + return erf_a if x.zero? + # Let f(x+a) = erf(x+a)*exp((x+a)**2)*sqrt(pi)/2 + # = c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ... + # f'(x+a) = 1+2*(x+a)*f(x+a) + # f'(x+a) = c1 + 2*c2*x + 3*c3*x**2 + 4*c4*x**3 + 5*c5*x**4 + ... + # = 1+2*(x+a)*(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + ...) + # therefore, + # c0 = f(a) + # c1 = 2 * a * c0 + 1 + # c2 = (2 * c0 + 2 * a * c1) / 2 + # c3 = (2 * c1 + 2 * a * c2) / 3 + # c4 = (2 * c2 + 2 * a * c3) / 4 + # + # All coefficients are positive when a >= 0 + + scale = BigDecimal(2).div(sqrt(PI(prec), prec), prec) + c_prev = erf_a.div(scale.mult(exp(-a*a, prec), prec), prec) + c_next = (2 * a * c_prev).add(1, prec).mult(x, prec) + sum = c_prev.add(c_next, prec) + + 2.step do |k| + cn = (c_prev.mult(x, prec) + a * c_next).mult(2, prec).mult(x, prec).div(k, prec) + sum = sum.add(cn, prec) + c_prev, c_next = c_next, cn + break if [c_prev, c_next].all? { |c| c.zero? || (c.exponent < sum.exponent - prec) } + end + value = sum.mult(scale.mult(exp(-(x + a).mult(x + a, prec), prec), prec), prec) + value > 1 ? BigDecimal(1) : value + end + + private_class_method def _erfc_asymptotic(x, prec) # :nodoc: + # Let f(x) = erfc(x)*sqrt(pi)*exp(x**2)/2 + # f(x) satisfies the following differential equation: + # 2*x*f(x) = f'(x) + 1 + # From the above equation, we can derive the following asymptotic expansion: + # f(x) = (0..kmax).sum { (-1)**k * (2*k)! / 4**k / k! / x**(2*k)) } / x + + # This asymptotic expansion does not converge. + # But if there is a k that satisfies (2*k)! / 4**k / k! / x**(2*k) < 10**(-prec), + # It is enough to calculate erfc within the given precision. + # Using Stirling's approximation, we can simplify this condition to: + # sqrt(2)/2 + k*log(k) - k - 2*k*log(x) < -prec*log(10) + # and the left side is minimized when k = x**2. + prec += BigDecimal::Internal::EXTRA_PREC + xf = x.to_f + kmax = (1..(xf ** 2).floor).bsearch do |k| + Math.log(2) / 2 + k * Math.log(k) - k - 2 * k * Math.log(xf) < -prec * Math.log(10) + end + return unless kmax + + sum = BigDecimal(1) + # To calculate `exp(x2, prec)`, x2 needs extra log10(x**2) digits of precision + x2 = x.mult(x, prec + (2 * Math.log10(xf)).ceil) + d = BigDecimal(1) + (1..kmax).each do |k| + d = d.div(x2, prec).mult(1 - 2 * k, prec).div(2, prec) + sum = sum.add(d, prec) + end + sum.div(exp(x2, prec).mult(PI(prec).sqrt(prec), prec), prec).div(x, prec) + end + + # call-seq: + # BigMath.gamma(decimal, numeric) -> BigDecimal + # + # Computes the gamma function of +decimal+ to the specified number of + # digits of precision, +numeric+. + # + # BigMath.gamma(BigDecimal('0.5'), 32).to_s + # #=> "0.17724538509055160272981674833411e1" + # + def gamma(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :gamma) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :gamma) + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + if x < 0.5 + raise Math::DomainError, 'Numerical argument is out of domain - gamma' if x.frac.zero? + + # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) + pi = PI(prec2) + sin = _sinpix(x, pi, prec2) + return pi.div(gamma(1 - x, prec2).mult(sin, prec2), prec) + elsif x.frac.zero? && x < 1000 * prec + return _gamma_positive_integer(x, prec2).mult(1, prec) + end + + a, sum = _gamma_spouge_sum_part(x, prec2) + (x + (a - 1)).power(x - 0.5, prec2).mult(BigMath.exp(1 - x, prec2), prec2).mult(sum, prec) + end + + # call-seq: + # BigMath.lgamma(decimal, numeric) -> [BigDecimal, Integer] + # + # Computes the natural logarithm of the absolute value of the gamma function + # of +decimal+ to the specified number of digits of precision, +numeric+ and its sign. + # + # BigMath.lgamma(BigDecimal('0.5'), 32) + # #=> [0.57236494292470008707171367567653e0, 1] + # + def lgamma(x, prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :lgamma) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, prec, :lgamma) + prec2 = prec + BigDecimal::Internal::EXTRA_PREC + if x < 0.5 + return [BigDecimal::INFINITY, 1] if x.frac.zero? + + loop do + # Euler's reflection formula: gamma(z) * gamma(1-z) = pi/sin(pi*z) + pi = PI(prec2) + sin = _sinpix(x, pi, prec2) + log_gamma = BigMath.log(pi, prec2).sub(lgamma(1 - x, prec2).first + BigMath.log(sin.abs, prec2), prec) + return [log_gamma, sin > 0 ? 1 : -1] if prec2 + log_gamma.exponent > prec + BigDecimal::Internal::EXTRA_PREC + + # Retry with higher precision if loss of significance is too large + prec2 = prec2 * 3 / 2 + end + elsif x.frac.zero? && x < 1000 * prec + log_gamma = BigMath.log(_gamma_positive_integer(x, prec2), prec) + [log_gamma, 1] + else + # if x is close to 1 or 2, increase precision to reduce loss of significance + diff1_exponent = (x - 1).exponent + diff2_exponent = (x - 2).exponent + extremely_near_one = diff1_exponent < -prec2 + extremely_near_two = diff2_exponent < -prec2 + + if extremely_near_one || extremely_near_two + # If x is extreamely close to base = 1 or 2, linear interpolation is accurate enough. + # Taylor expansion at x = base is: (x - base) * digamma(base) + (x - base) ** 2 * trigamma(base) / 2 + ... + # And we can ignore (x - base) ** 2 and higher order terms. + base = extremely_near_one ? 1 : 2 + d = BigDecimal(1)._decimal_shift(1 - prec2) + log_gamma_d, sign = lgamma(base + d, prec2) + return [log_gamma_d.mult(x - base, prec2).div(d, prec), sign] + end + + prec2 += [-diff1_exponent, -diff2_exponent, 0].max + a, sum = _gamma_spouge_sum_part(x, prec2) + log_gamma = BigMath.log(sum, prec2).add((x - 0.5).mult(BigMath.log(x.add(a - 1, prec2), prec2), prec2) + 1 - x, prec) + [log_gamma, 1] + end + end + + # Returns sum part: sqrt(2*pi) and c[k]/(x+k) terms of Spouge's approximation + private_class_method def _gamma_spouge_sum_part(x, prec) # :nodoc: + x -= 1 + # Spouge's approximation + # x! = (x + a)**(x + 0.5) * exp(-x - a) * (sqrt(2 * pi) + (1..a - 1).sum{|k| c[k] / (x + k) } + epsilon) + # where c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(a - k) / (k - 1)! + # and epsilon is bounded by a**(-0.5) * (2 * pi) ** (-a - 0.5) + + # Estimate required a for given precision + a = (prec / Math.log10(2 * Math::PI)).ceil + + # Calculate exponent of c[k] in low precision to estimate required precision + low_prec = 16 + log10f = Math.log(10) + x_low_prec = x.mult(1, low_prec) + loggamma_k = 0 + ck_exponents = (1..a-1).map do |k| + loggamma_k += Math.log10(k - 1) if k > 1 + -loggamma_k - k / log10f + (k - 0.5) * Math.log10(a - k) - BigDecimal::Internal.float_log(x_low_prec.add(k, low_prec)) / log10f + end + + # Estimate exponent of sum by Stirling's approximation + approx_sum_exponent = x < 1 ? -Math.log10(a) / 2 : Math.log10(2 * Math::PI) / 2 + x_low_prec.add(0.5, low_prec) * Math.log10(x_low_prec / x_low_prec.add(a, low_prec)) + + # Determine required precision of c[k] + prec2 = [ck_exponents.max.ceil - approx_sum_exponent.floor, 0].max + prec + + einv = BigMath.exp(-1, prec2) + sum = (PI(prec) * 2).sqrt(prec).mult(BigMath.exp(-a, prec), prec) + y = BigDecimal(1) + (1..a - 1).each do |k| + # c[k] = (-1)**k * (a - k)**(k - 0.5) * exp(-k) / (k-1)! / (x + k) + y = y.div(1 - k, prec2) if k > 1 + y = y.mult(einv, prec2) + z = y.mult(BigDecimal((a - k) ** k), prec2).div(BigDecimal(a - k).sqrt(prec2).mult(x.add(k, prec2), prec2), prec2) + # sum += c[k] / (x + k) + sum = sum.add(z, prec2) + end + [a, sum] + end + + private_class_method def _gamma_positive_integer(x, prec) # :nodoc: + return x if x == 1 + numbers = (1..x - 1).map {|i| BigDecimal(i) } + while numbers.size > 1 + numbers = numbers.each_slice(2).map {|a, b| b ? a.mult(b, prec) : a } + end + numbers.first + end + + # Returns sin(pi * x), for gamma reflection formula calculation + private_class_method def _sinpix(x, pi, prec) # :nodoc: + x = x % 2 + sign = x > 1 ? -1 : 1 + x %= 1 + x = 1 - x if x > 0.5 # to avoid sin(pi*x) loss of precision for x close to 1 + sign * sin(x.mult(pi, prec), prec) + end + + # call-seq: + # frexp(x) -> [BigDecimal, Integer] + # + # Decomposes +x+ into a normalized fraction and an integral power of ten. + # + # BigMath.frexp(BigDecimal(123.456)) + # #=> [0.123456e0, 3] + # + def frexp(x) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :frexp) + return [x, 0] unless x.finite? + + exponent = x.exponent + [x._decimal_shift(-exponent), exponent] + end + + # call-seq: + # ldexp(fraction, exponent) -> BigDecimal + # + # Inverse of +frexp+. + # Returns the value of fraction * 10**exponent. + # + # BigMath.ldexp(BigDecimal("0.123456e0"), 3) + # #=> 0.123456e3 + # + def ldexp(x, exponent) + x = BigDecimal::Internal.coerce_to_bigdecimal(x, 0, :ldexp) + x.finite? ? x._decimal_shift(exponent) : x + end + + # call-seq: + # PI(numeric) -> BigDecimal + # + # Computes the value of pi to the specified number of digits of precision, + # +numeric+. + # + # BigMath.PI(32).to_s + # #=> "0.31415926535897932384626433832795e1" + # + def PI(prec) + # Gauss–Legendre algorithm + prec = BigDecimal::Internal.coerce_validate_prec(prec, :PI) + n = prec + BigDecimal::Internal::EXTRA_PREC + a = BigDecimal(1) + b = BigDecimal(0.5, 0).sqrt(n) + s = BigDecimal(0.25, 0) + t = 1 + while a != b && (a - b).exponent > 1 - n + c = (a - b).div(2, n) + a, b = (a + b).div(2, n), (a * b).sqrt(n) + s = s.sub(c * c * t, n) + t *= 2 + end + (a * b).div(s, prec) + end + + # call-seq: + # E(numeric) -> BigDecimal + # + # Computes e (the base of natural logarithms) to the specified number of + # digits of precision, +numeric+. + # + # BigMath.E(32).to_s + # #=> "0.27182818284590452353602874713527e1" + # + def E(prec) + prec = BigDecimal::Internal.coerce_validate_prec(prec, :E) + exp(1, prec) + end +end diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/newton.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/newton.rb new file mode 100644 index 0000000..bb70d47 --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/newton.rb @@ -0,0 +1,82 @@ +# frozen_string_literal: false +require "bigdecimal/ludcmp" +require "bigdecimal/jacobian" + +warn "'bigdecimal/newton' is deprecated and will be removed in a future release." + +# +# newton.rb +# +# Solves the nonlinear algebraic equation system f = 0 by Newton's method. +# This program is not dependent on BigDecimal. +# +# To call: +# n = nlsolve(f,x) +# where n is the number of iterations required, +# x is the initial value vector +# f is an Object which is used to compute the values of the equations to be solved. +# It must provide the following methods: +# +# f.values(x):: returns the values of all functions at x +# +# f.zero:: returns 0.0 +# f.one:: returns 1.0 +# f.two:: returns 2.0 +# f.ten:: returns 10.0 +# +# f.eps:: returns the convergence criterion (epsilon value) used to determine whether two values are considered equal. If |a-b| < epsilon, the two values are considered equal. +# +# On exit, x is the solution vector. +# +module Newton + include LUSolve + include Jacobian + module_function + + def norm(fv,zero=0.0) # :nodoc: + s = zero + n = fv.size + for i in 0...n do + s += fv[i]*fv[i] + end + s + end + + # See also Newton + def nlsolve(f,x) + nRetry = 0 + n = x.size + + f0 = f.values(x) + zero = f.zero + one = f.one + two = f.two + p5 = one/two + d = norm(f0,zero) + minfact = f.ten*f.ten*f.ten + minfact = one/minfact + e = f.eps + while d >= e do + nRetry += 1 + # Not yet converged. => Compute Jacobian matrix + dfdx = jacobian(f,f0,x) + # Solve dfdx*dx = -f0 to estimate dx + dx = lusolve(dfdx,f0,ludecomp(dfdx,n,zero,one),zero) + fact = two + xs = x.dup + begin + fact *= p5 + if fact < minfact then + raise "Failed to reduce function values." + end + for i in 0...n do + x[i] = xs[i] - dx[i]*fact + end + f0 = f.values(x) + dn = norm(f0,zero) + end while(dn>=d) + d = dn + end + nRetry + end +end diff --git a/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/util.rb b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/util.rb new file mode 100644 index 0000000..7f2232d --- /dev/null +++ b/vendor/bundle/ruby/3.4.0/gems/bigdecimal-4.1.2/lib/bigdecimal/util.rb @@ -0,0 +1,186 @@ +# frozen_string_literal: true +# +#-- +# bigdecimal/util extends various native classes to provide the #to_d method, +# and provides BigDecimal#to_d and BigDecimal#to_digits. +#++ + +require 'bigdecimal' + +class Integer < Numeric + # call-seq: + # int.to_d -> bigdecimal + # + # Returns the value of +int+ as a BigDecimal. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # 42.to_d # => 0.42e2 + # + # See also Kernel.BigDecimal. + # + def to_d + BigDecimal(self) + end +end + + +class Float < Numeric + # call-seq: + # float.to_d -> bigdecimal + # float.to_d(precision) -> bigdecimal + # + # Returns the value of +float+ as a BigDecimal. + # The +precision+ parameter is used to determine the number of + # significant digits for the result. When +precision+ is set to +0+, + # the number of digits to represent the float being converted is determined + # automatically. + # The default +precision+ is +0+. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # 0.5.to_d # => 0.5e0 + # 1.234.to_d # => 0.1234e1 + # 1.234.to_d(2) # => 0.12e1 + # + # See also Kernel.BigDecimal. + # + def to_d(precision=0) + BigDecimal(self, precision) + end +end + + +class String + # call-seq: + # str.to_d -> bigdecimal + # + # Returns the result of interpreting leading characters in +str+ + # as a BigDecimal. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # "0.5".to_d # => 0.5e0 + # "123.45e1".to_d # => 0.12345e4 + # "45.67 degrees".to_d # => 0.4567e2 + # + # See also Kernel.BigDecimal. + # + def to_d + BigDecimal.interpret_loosely(self) + end +end + + +class BigDecimal < Numeric + # call-seq: + # a.to_digits -> string + # + # Converts a BigDecimal to a String of the form "nnnnnn.mmm". + # This method is deprecated; use BigDecimal#to_s("F") instead. + # + # require 'bigdecimal/util' + # + # d = BigDecimal("3.14") + # d.to_digits # => "3.14" + # + def to_digits + if self.nan? || self.infinite? || self.zero? + self.to_s + else + i = self.to_i.to_s + _,f,_,z = self.frac.split + i + "." + ("0"*(-z)) + f + end + end + + # call-seq: + # a.to_d -> bigdecimal + # + # Returns self. + # + # require 'bigdecimal/util' + # + # d = BigDecimal("3.14") + # d.to_d # => 0.314e1 + # + def to_d + self + end +end + + +class Rational < Numeric + # call-seq: + # rat.to_d(precision) -> bigdecimal + # + # Returns the value as a BigDecimal. + # + # The +precision+ parameter is used to determine the number of + # significant digits for the result. When +precision+ is set to +0+, + # the number of digits to represent the float being converted is determined + # automatically. + # The default +precision+ is +0+. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # Rational(22, 7).to_d(3) # => 0.314e1 + # + # See also Kernel.BigDecimal. + # + def to_d(precision=0) + BigDecimal(self, precision) + end +end + + +class Complex < Numeric + # call-seq: + # cmp.to_d -> bigdecimal + # cmp.to_d(precision) -> bigdecimal + # + # Returns the value as a BigDecimal. + # If the imaginary part is not +0+, an error is raised + # + # The +precision+ parameter is used to determine the number of + # significant digits for the result. When +precision+ is set to +0+, + # the number of digits to represent the float being converted is determined + # automatically. + # The default +precision+ is +0+. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # Complex(0.1234567, 0).to_d(4) # => 0.1235e0 + # Complex(Rational(22, 7), 0).to_d(3) # => 0.314e1 + # Complex(1, 1).to_d # raises ArgumentError + # + # See also Kernel.BigDecimal. + # + def to_d(precision=0) + BigDecimal(self) unless self.imag.zero? # to raise error + + BigDecimal(self.real, precision) + end +end + + +class NilClass + # call-seq: + # nil.to_d -> bigdecimal + # + # Returns nil represented as a BigDecimal. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # nil.to_d # => 0.0 + # + def to_d + BigDecimal(0) + end +end |