From fcbf63e62c627deae76c1b8cb8c0876c536ed811 Mon Sep 17 00:00:00 2001
From: Jari Vetoniemi <jari.vetoniemi@indooratlas.com>
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
+ # # => #<BigDecimal:1008ef070,'0.42E2',9(36)>
+ #
+ 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
+ # # => #<BigDecimal:1dc69e0,'0.5E0',9(18)>
+ #
+ 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
+ # # => #<BigDecimal:1dc69e0,'0.5E0',9(18)>
+ #
+ 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)
+ # # => #<BigDecimal:1a44d08,'0.314E1',18(36)>
+ 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-70-g09d2