summaryrefslogtreecommitdiff
path: root/ext/bigdecimal/lib
diff options
context:
space:
mode:
Diffstat (limited to 'ext/bigdecimal/lib')
-rw-r--r--ext/bigdecimal/lib/bigdecimal-rational.rb30
-rw-r--r--ext/bigdecimal/lib/delcr7
-rw-r--r--ext/bigdecimal/lib/jacobian.rb63
-rw-r--r--ext/bigdecimal/lib/linear.rb46
-rw-r--r--ext/bigdecimal/lib/ludcmp.rb75
-rw-r--r--ext/bigdecimal/lib/newton.rb75
-rw-r--r--ext/bigdecimal/lib/nlsolve.rb38
-rw-r--r--ext/bigdecimal/lib/pai.rb49
8 files changed, 383 insertions, 0 deletions
diff --git a/ext/bigdecimal/lib/bigdecimal-rational.rb b/ext/bigdecimal/lib/bigdecimal-rational.rb
new file mode 100644
index 0000000000..de999d19cf
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal-rational.rb
@@ -0,0 +1,30 @@
+#
+# BigDecimal <-> Rational
+#
+class BigDecimal
+ # Convert BigDecimal to Rational
+ def to_r
+ sign,digits,base,power = self.to_parts
+ numerator = sign*digits.to_i
+ denomi_power = power - digits.size # base is always 10
+ if denomi_power < 0
+ denominator = base ** (-denomi_power)
+ else
+ denominator = base ** denomi_power
+ end
+ Rational.new(numerator,denominator)
+ end
+end
+
+class Rational
+ # Convert Rational to BigDecimal
+ # to_d returns an array [quotient,residue]
+ def to_d(nFig=0)
+ num = self.numerator.to_s
+ if nFig<=0
+ nFig = BigDecimal.double_fig*2+1
+ end
+ BigDecimal.new(num).div(self.denominator,nFig)
+ end
+end
+
diff --git a/ext/bigdecimal/lib/delcr b/ext/bigdecimal/lib/delcr
new file mode 100644
index 0000000000..df4483ce34
--- /dev/null
+++ b/ext/bigdecimal/lib/delcr
@@ -0,0 +1,7 @@
+tr -d '\r' < bigdecimal-rational.rb > unix/bigdecimal-rational.rb
+tr -d '\r' <jacobian.rb > unix/jacobian.rb
+tr -d '\r' <linear.rb > unix/linear.rb
+tr -d '\r' <ludcmp.rb > unix/ludcmp.rb
+tr -d '\r' <newton.rb > unix/newton.rb
+tr -d '\r' <nlsolve.rb > unix/nlsolve.rb
+tr -d '\r' <pai.rb > unix/pai.rb
diff --git a/ext/bigdecimal/lib/jacobian.rb b/ext/bigdecimal/lib/jacobian.rb
new file mode 100644
index 0000000000..34a60ae67a
--- /dev/null
+++ b/ext/bigdecimal/lib/jacobian.rb
@@ -0,0 +1,63 @@
+#
+# jacobian.rb
+#
+# Computes Jacobian matrix of f at x
+#
+module Jacobian
+ 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
+
+ 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
+ s = f.zero
+ deriv = []
+ if(nRetry>100) then
+ raize "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
+
+ 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/ext/bigdecimal/lib/linear.rb b/ext/bigdecimal/lib/linear.rb
new file mode 100644
index 0000000000..f93404fb6f
--- /dev/null
+++ b/ext/bigdecimal/lib/linear.rb
@@ -0,0 +1,46 @@
+#!/usr/local/bin/ruby
+
+#
+# linear.rb
+#
+# Solves linear equation system(A*x = b) by LU decomposition method.
+# where A is a coefficient matrix,x is an answer vector,b is a constant vector.
+#
+require "bigdecimal"
+require "ludcmp"
+
+include LUSolve
+
+def rd_order
+ printf("Number of equations ?")
+ n = gets().chomp.to_i
+end
+
+zero = BigDecimal::new("0.0")
+one = BigDecimal::new("1.0")
+
+while (n=rd_order())>0
+ a = []
+ as= []
+ b = []
+ printf("\nEnter coefficient matrix element A[i,j]\n");
+ for i in 0...n do
+ for j in 0...n do
+ printf("A[%d,%d]? ",i,j); s = gets
+ a <<=BigDecimal::new(s);
+ as<<=BigDecimal::new(s);
+ end
+ printf("Contatant vector element b[%d] ? ",i);b<<=BigDecimal::new(gets);
+ end
+ printf "ANS="
+ x = lusolve(a,b,ludecomp(a,n,zero,one),zero)
+ p x
+ printf "A*x-b\n"
+ for i in 0...n do
+ s = zero
+ for j in 0...n do
+ s = s + as[i*n+j]*x[j]
+ end
+ p s-b[i]
+ end
+end
diff --git a/ext/bigdecimal/lib/ludcmp.rb b/ext/bigdecimal/lib/ludcmp.rb
new file mode 100644
index 0000000000..c36f0dea5b
--- /dev/null
+++ b/ext/bigdecimal/lib/ludcmp.rb
@@ -0,0 +1,75 @@
+#
+# ludcmp.rb
+#
+module LUSolve
+ def ludecomp(a,n,zero=0.0,one=1.0)
+ 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/nrmrow
+ 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]/pivot
+ if mult!=zero then
+ pskn = ps[k]*n
+ for j in (k+1)...n do
+ a[psin+j] -= mult*a[pskn+j]
+ end
+ end
+ end
+ end
+ raise "Singular matrix" if a[ps[n1]*n+n1] == zero
+ ps
+ end
+
+ def lusolve(a,b,ps,zero=0.0)
+ 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]*x[j] + 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]*x[j] + dot
+ end
+ x[i] = (x[i]-dot)/a[psin+i]
+ end
+ x
+ end
+end
diff --git a/ext/bigdecimal/lib/newton.rb b/ext/bigdecimal/lib/newton.rb
new file mode 100644
index 0000000000..67a92474ac
--- /dev/null
+++ b/ext/bigdecimal/lib/newton.rb
@@ -0,0 +1,75 @@
+#
+# newton.rb
+#
+# Solves 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 solution vector.
+# f is the object to be solved which must have following methods.
+#
+# f ... Object to compute Jacobian matrix of the equation systems.
+# [Methods required for f]
+# f.values(x) returns values of all functions at x.
+# f.zero returns 0.0
+# f.one returns 1.0
+# f.two returns 1.0
+# f.ten returns 10.0
+# f.eps convergence criterion
+# x ... initial values
+#
+require "ludcmp"
+require "jacobian"
+
+module Newton
+ include LUSolve
+ include Jacobian
+
+ def norm(fv,zero=0.0)
+ s = zero
+ n = fv.size
+ for i in 0...n do
+ s += fv[i]*fv[i]
+ end
+ s
+ end
+
+ 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
+ raize "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/ext/bigdecimal/lib/nlsolve.rb b/ext/bigdecimal/lib/nlsolve.rb
new file mode 100644
index 0000000000..08f17f9ecd
--- /dev/null
+++ b/ext/bigdecimal/lib/nlsolve.rb
@@ -0,0 +1,38 @@
+#!/usr/local/bin/ruby
+
+#
+# nlsolve.rb
+# An example for solving nonlinear algebraic equation system.
+#
+
+require "bigdecimal"
+require "newton"
+include Newton
+
+class Function
+ def initialize()
+ @zero = BigDecimal::new("0.0")
+ @one = BigDecimal::new("1.0")
+ @two = BigDecimal::new("2.0")
+ @ten = BigDecimal::new("10.0")
+ @eps = BigDecimal::new("1.0e-16")
+ end
+ def zero;@zero;end
+ def one ;@one ;end
+ def two ;@two ;end
+ def ten ;@ten ;end
+ def eps ;@eps ;end
+ def values(x) # <= defines functions solved
+ f = []
+ f1 = x[0]*x[0] + x[1]*x[1] - @two # f1 = x**2 + y**2 - 2 => 0
+ f2 = x[0] - x[1] # f2 = x - y => 0
+ f <<= f1
+ f <<= f2
+ f
+ end
+end
+ f = BigDecimal::limit(100)
+ f = Function.new
+ x = [f.zero,f.zero] # Initial values
+ n = nlsolve(f,x)
+ p x
diff --git a/ext/bigdecimal/lib/pai.rb b/ext/bigdecimal/lib/pai.rb
new file mode 100644
index 0000000000..4e2dd7104a
--- /dev/null
+++ b/ext/bigdecimal/lib/pai.rb
@@ -0,0 +1,49 @@
+#!/usr/local/bin/ruby
+
+#
+# pai.rb
+#
+
+require "bigdecimal"
+#
+# Calculates 3.1415.... using J. Machin's formula.
+#
+def pai(sig) # sig: Number of significant figures
+ exp = -sig
+ pi = BigDecimal::new("0")
+ two = BigDecimal::new("2")
+ m25 = BigDecimal::new("-0.04")
+ m57121 = BigDecimal::new("-57121")
+
+ u = BigDecimal::new("1")
+ k = BigDecimal::new("1")
+ w = BigDecimal::new("1")
+ t = BigDecimal::new("-80")
+ while (u.exponent >= exp)
+ t = t*m25
+ u,r = t.div(k,sig)
+ pi = pi + u
+ k = k+two
+ end
+
+ u = BigDecimal::new("1")
+ k = BigDecimal::new("1")
+ w = BigDecimal::new("1")
+ t = BigDecimal::new("956")
+ while (u.exponent >= exp )
+ t,r = t.div(m57121,sig)
+ u,r = t.div(k,sig)
+ pi = pi + u
+ k = k+two
+ end
+ pi
+end
+
+if $0 == __FILE__
+ if ARGV.size == 1
+ print "PAI("+ARGV[0]+"):\n"
+ p pai(ARGV[0].to_i)
+ else
+ print "TRY: ruby pai.rb 1000 \n"
+ end
+end
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-27 03:50:36 +0000