From 77443517086c68cb51f26f68f5d25279f1be7daf Mon Sep 17 00:00:00 2001 From: shigek Date: Fri, 28 Mar 2003 05:00:21 +0000 Subject: Copied from rough/bigdecimal,documents & some sample programs added. git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@3625 b2dd03c8-39d4-4d8f-98ff-823fe69b080e --- ext/bigdecimal/lib/bigdecimal-rational.rb | 30 +++++++++++++ ext/bigdecimal/lib/delcr | 7 +++ ext/bigdecimal/lib/jacobian.rb | 63 ++++++++++++++++++++++++++ ext/bigdecimal/lib/linear.rb | 46 +++++++++++++++++++ ext/bigdecimal/lib/ludcmp.rb | 75 +++++++++++++++++++++++++++++++ ext/bigdecimal/lib/newton.rb | 75 +++++++++++++++++++++++++++++++ ext/bigdecimal/lib/nlsolve.rb | 38 ++++++++++++++++ ext/bigdecimal/lib/pai.rb | 49 ++++++++++++++++++++ 8 files changed, 383 insertions(+) create mode 100644 ext/bigdecimal/lib/bigdecimal-rational.rb create mode 100644 ext/bigdecimal/lib/delcr create mode 100644 ext/bigdecimal/lib/jacobian.rb create mode 100644 ext/bigdecimal/lib/linear.rb create mode 100644 ext/bigdecimal/lib/ludcmp.rb create mode 100644 ext/bigdecimal/lib/newton.rb create mode 100644 ext/bigdecimal/lib/nlsolve.rb create mode 100644 ext/bigdecimal/lib/pai.rb (limited to 'ext/bigdecimal/lib') 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' unix/jacobian.rb +tr -d '\r' unix/linear.rb +tr -d '\r' unix/ludcmp.rb +tr -d '\r' unix/newton.rb +tr -d '\r' unix/nlsolve.rb +tr -d '\r' 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 -- cgit v1.2.3