From fcbf63e62c627deae76c1b8cb8c0876c536ed811 Mon Sep 17 00:00:00 2001 From: Jari Vetoniemi Date: Mon, 16 Mar 2020 18:49:26 +0900 Subject: Fresh start --- jni/ruby/ext/bigdecimal/lib/bigdecimal/jacobian.rb | 87 ++++++++ jni/ruby/ext/bigdecimal/lib/bigdecimal/ludcmp.rb | 88 ++++++++ jni/ruby/ext/bigdecimal/lib/bigdecimal/math.rb | 231 +++++++++++++++++++++ jni/ruby/ext/bigdecimal/lib/bigdecimal/newton.rb | 79 +++++++ jni/ruby/ext/bigdecimal/lib/bigdecimal/util.rb | 127 +++++++++++ 5 files changed, 612 insertions(+) create mode 100644 jni/ruby/ext/bigdecimal/lib/bigdecimal/jacobian.rb create mode 100644 jni/ruby/ext/bigdecimal/lib/bigdecimal/ludcmp.rb create mode 100644 jni/ruby/ext/bigdecimal/lib/bigdecimal/math.rb create mode 100644 jni/ruby/ext/bigdecimal/lib/bigdecimal/newton.rb create mode 100644 jni/ruby/ext/bigdecimal/lib/bigdecimal/util.rb (limited to 'jni/ruby/ext/bigdecimal/lib') diff --git a/jni/ruby/ext/bigdecimal/lib/bigdecimal/jacobian.rb b/jni/ruby/ext/bigdecimal/lib/bigdecimal/jacobian.rb new file mode 100644 index 0000000..d56caab --- /dev/null +++ b/jni/ruby/ext/bigdecimal/lib/bigdecimal/jacobian.rb @@ -0,0 +1,87 @@ +# +# 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/jni/ruby/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/jni/ruby/ext/bigdecimal/lib/bigdecimal/ludcmp.rb new file mode 100644 index 0000000..6cbe29b --- /dev/null +++ b/jni/ruby/ext/bigdecimal/lib/bigdecimal/ludcmp.rb @@ -0,0 +1,88 @@ +require 'bigdecimal' + +# +# 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/jni/ruby/ext/bigdecimal/lib/bigdecimal/math.rb b/jni/ruby/ext/bigdecimal/lib/bigdecimal/math.rb new file mode 100644 index 0000000..4a4fcc2 --- /dev/null +++ b/jni/ruby/ext/bigdecimal/lib/bigdecimal/math.rb @@ -0,0 +1,231 @@ +require 'bigdecimal' + +# +#-- +# Contents: +# sqrt(x, prec) +# sin (x, prec) +# cos (x, prec) +# atan(x, prec) Note: |x|<1, x=0.9999 may not converge. +# PI (prec) +# E (prec) == exp(1.0,prec) +# +# where: +# x ... BigDecimal number to be computed. +# |x| must be small enough to get convergence. +# prec ... Number of digits to be obtained. +#++ +# +# Provides mathematical functions. +# +# Example: +# +# require "bigdecimal/math" +# +# include BigMath +# +# a = BigDecimal((PI(100)/2).to_s) +# puts sin(a,100) # => 0.10000000000000000000......E1 +# +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.new('2'), 16).to_s + # #=> "0.1414213562373095048801688724E1" + # + def sqrt(x, prec) + x.sqrt(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, 5).to_s + # #=> "0.70710678118654752440082036563292800375E0" + # + def sin(x, prec) + raise ArgumentError, "Zero or negative precision for sin" if prec <= 0 + return BigDecimal("NaN") if x.infinite? || x.nan? + n = prec + BigDecimal.double_fig + one = BigDecimal("1") + two = BigDecimal("2") + x = -x if neg = x < 0 + if x > (twopi = two * BigMath.PI(prec)) + if x > 30 + x %= twopi + else + x -= twopi while x > twopi + end + end + x1 = x + x2 = x.mult(x,n) + sign = 1 + y = x + d = y + i = one + z = one + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + sign = -sign + x1 = x2.mult(x1,n) + i += two + z *= (i-one) * i + d = sign * x1.div(z,m) + y += d + end + neg ? -y : y + 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(4), 16).to_s + # #=> "-0.999999999999999999999999999999856613163740061349E0" + # + def cos(x, prec) + raise ArgumentError, "Zero or negative precision for cos" if prec <= 0 + return BigDecimal("NaN") if x.infinite? || x.nan? + n = prec + BigDecimal.double_fig + one = BigDecimal("1") + two = BigDecimal("2") + x = -x if x < 0 + if x > (twopi = two * BigMath.PI(prec)) + if x > 30 + x %= twopi + else + x -= twopi while x > twopi + end + end + x1 = one + x2 = x.mult(x,n) + sign = 1 + y = one + d = y + i = BigDecimal("0") + z = one + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + sign = -sign + x1 = x2.mult(x1,n) + i += two + z *= (i-one) * i + d = sign * x1.div(z,m) + y += d + end + y + 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.new('-1'), 16).to_s + # #=> "-0.785398163397448309615660845819878471907514682065E0" + # + def atan(x, prec) + raise ArgumentError, "Zero or negative precision for atan" if prec <= 0 + return BigDecimal("NaN") if x.nan? + pi = PI(prec) + x = -x if neg = x < 0 + return pi.div(neg ? -2 : 2, prec) if x.infinite? + return pi / (neg ? -4 : 4) if x.round(prec) == 1 + x = BigDecimal("1").div(x, prec) if inv = x > 1 + x = (-1 + sqrt(1 + x**2, prec))/x if dbl = x > 0.5 + n = prec + BigDecimal.double_fig + y = x + d = y + t = x + r = BigDecimal("3") + x2 = x.mult(x,n) + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + t = -t.mult(x2,n) + d = t.div(r,m) + y += d + r += 2 + end + y *= 2 if dbl + y = pi / 2 - y if inv + y = -y if neg + y + end + + # call-seq: + # PI(numeric) -> BigDecimal + # + # Computes the value of pi to the specified number of digits of precision, + # +numeric+. + # + # BigMath.PI(10).to_s + # #=> "0.3141592653589793238462643388813853786957412E1" + # + def PI(prec) + raise ArgumentError, "Zero or negative precision for PI" if prec <= 0 + n = prec + BigDecimal.double_fig + zero = BigDecimal("0") + one = BigDecimal("1") + two = BigDecimal("2") + + m25 = BigDecimal("-0.04") + m57121 = BigDecimal("-57121") + + pi = zero + + d = one + k = one + t = BigDecimal("-80") + while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + t = t*m25 + d = t.div(k,m) + k = k+two + pi = pi + d + end + + d = one + k = one + t = BigDecimal("956") + while d.nonzero? && ((m = n - (pi.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + t = t.div(m57121,n) + d = t.div(k,m) + pi = pi + d + k = k+two + end + pi + end + + # call-seq: + # E(numeric) -> BigDecimal + # + # Computes e (the base of natural logarithms) to the specified number of + # digits of precision, +numeric+. + # + # BigMath.E(10).to_s + # #=> "0.271828182845904523536028752390026306410273E1" + # + def E(prec) + raise ArgumentError, "Zero or negative precision for E" if prec <= 0 + BigMath.exp(1, prec) + end +end diff --git a/jni/ruby/ext/bigdecimal/lib/bigdecimal/newton.rb b/jni/ruby/ext/bigdecimal/lib/bigdecimal/newton.rb new file mode 100644 index 0000000..db1a5ad --- /dev/null +++ b/jni/ruby/ext/bigdecimal/lib/bigdecimal/newton.rb @@ -0,0 +1,79 @@ +require "bigdecimal/ludcmp" +require "bigdecimal/jacobian" + +# +# 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/jni/ruby/ext/bigdecimal/lib/bigdecimal/util.rb b/jni/ruby/ext/bigdecimal/lib/bigdecimal/util.rb new file mode 100644 index 0000000..82c82c8 --- /dev/null +++ b/jni/ruby/ext/bigdecimal/lib/bigdecimal/util.rb @@ -0,0 +1,127 @@ +# BigDecimal extends the native Integer class to provide the #to_d method. +# +# When you require the BigDecimal library in your application, this methodwill +# be available on Integer objects. +class Integer < Numeric + # call-seq: + # int.to_d -> bigdecimal + # + # Convert +int+ to a BigDecimal and return it. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # 42.to_d + # # => # + # + def to_d + BigDecimal(self) + end +end + +# BigDecimal extends the native Float class to provide the #to_d method. +# +# When you require BigDecimal in your application, this method will be +# available on Float objects. +class Float < Numeric + # call-seq: + # flt.to_d -> bigdecimal + # + # Convert +flt+ to a BigDecimal and return it. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # 0.5.to_d + # # => # + # + def to_d(precision=nil) + BigDecimal(self, precision || Float::DIG) + end +end + +# BigDecimal extends the native String class to provide the #to_d method. +# +# When you require BigDecimal in your application, this method will be +# available on String objects. +class String + # call-seq: + # string.to_d -> bigdecimal + # + # Convert +string+ to a BigDecimal and return it. + # + # require 'bigdecimal' + # require 'bigdecimal/util' + # + # "0.5".to_d + # # => # + # + def to_d + BigDecimal(self) + end +end + +# BigDecimal extends the native Numeric class to provide the #to_digits and +# #to_d methods. +# +# When you require BigDecimal in your application, this method will be +# available on BigDecimal objects. +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' + # require 'bigdecimal/util' + # + # d = BigDecimal.new("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. + def to_d + self + end +end + +# BigDecimal extends the native Rational class to provide the #to_d method. +# +# When you require BigDecimal in your application, this method will be +# available on Rational objects. +class Rational < Numeric + # call-seq: + # r.to_d(precision) -> bigdecimal + # + # Converts a Rational to a BigDecimal. + # + # The required +precision+ parameter is used to determine the amount of + # significant digits for the result. See BigDecimal#div for more information, + # as it is used along with the #denominator and the +precision+ for + # parameters. + # + # r = (22/7.0).to_r + # # => (7077085128725065/2251799813685248) + # r.to_d(3) + # # => # + def to_d(precision) + if precision <= 0 + raise ArgumentError, "negative precision" + end + num = self.numerator + BigDecimal(num).div(self.denominator, precision) + end +end -- cgit v1.2.3