diff options
Diffstat (limited to 'ruby_1_8_6/ext/bigdecimal/lib/bigdecimal')
| -rw-r--r-- | ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/jacobian.rb | 85 | ||||
| -rw-r--r-- | ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb | 84 | ||||
| -rw-r--r-- | ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/math.rb | 235 | ||||
| -rw-r--r-- | ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/newton.rb | 77 | ||||
| -rw-r--r-- | ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/util.rb | 65 |
5 files changed, 546 insertions, 0 deletions
diff --git a/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/jacobian.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/jacobian.rb new file mode 100644 index 0000000000..d80eeab901 --- /dev/null +++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/jacobian.rb @@ -0,0 +1,85 @@ +# +# 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 1.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 + #-- + 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 + 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 + + # 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/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb new file mode 100644 index 0000000000..8f4888725e --- /dev/null +++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/ludcmp.rb @@ -0,0 +1,84 @@ +# +# Solves a*x = b for x, using LU decomposition. +# +module LUSolve + # 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/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/math.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/math.rb new file mode 100644 index 0000000000..f3248a3c5c --- /dev/null +++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/math.rb @@ -0,0 +1,235 @@ +# +#-- +# Contents: +# sqrt(x, prec) +# sin (x, prec) +# cos (x, prec) +# atan(x, prec) Note: |x|<1, x=0.9999 may not converge. +# exp (x, prec) +# log (x, prec) +# 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" +# require "bigdecimal/math" +# +# include BigMath +# +# a = BigDecimal((PI(100)/2).to_s) +# puts sin(a,100) # -> 0.10000000000000000000......E1 +# +module BigMath + + # Computes the square root of x to the specified number of digits of + # precision. + # + # BigDecimal.new('2').sqrt(16).to_s + # -> "0.14142135623730950488016887242096975E1" + # + def sqrt(x,prec) + x.sqrt(prec) + end + + # Computes the sine of x to the specified number of digits of precision. + # + # If x is infinite or NaN, returns NaN. + 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") + 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 + y + end + + # Computes the cosine of x to the specified number of digits of precision. + # + # If x is infinite or NaN, returns NaN. + 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") + 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 + + # Computes the arctangent of x to the specified number of digits of precision. + # + # If x is infinite or NaN, returns NaN. + # Raises an argument error if x > 1. + def atan(x, prec) + raise ArgumentError, "Zero or negative precision for atan" if prec <= 0 + return BigDecimal("NaN") if x.infinite? || x.nan? + raise ArgumentError, "x.abs must be less than 1.0" if x.abs>=1 + 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 + end + + # Computes the value of e (the base of natural logarithms) raised to the + # power of x, to the specified number of digits of precision. + # + # If x is infinite or NaN, returns NaN. + # + # BigMath::exp(BigDecimal.new('1'), 10).to_s + # -> "0.271828182845904523536028752390026306410273E1" + def exp(x, prec) + raise ArgumentError, "Zero or negative precision for exp" if prec <= 0 + return BigDecimal("NaN") if x.infinite? || x.nan? + n = prec + BigDecimal.double_fig + one = BigDecimal("1") + x1 = one + y = one + d = y + z = one + i = 0 + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + x1 = x1.mult(x,n) + i += 1 + z *= i + d = x1.div(z,m) + y += d + end + y + end + + # Computes the natural logarithm of x to the specified number of digits + # of precision. + # + # Returns x if x is infinite or NaN. + # + def log(x, prec) + raise ArgumentError, "Zero or negative argument for log" if x <= 0 || prec <= 0 + return x if x.infinite? || x.nan? + one = BigDecimal("1") + two = BigDecimal("2") + n = prec + BigDecimal.double_fig + x = (x - one).div(x + one,n) + x2 = x.mult(x,n) + y = x + d = y + i = one + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + x = x2.mult(x,n) + i += two + d = x.div(i,m) + y += d + end + y*two + end + + # Computes the value of pi to the specified number of digits of precision. + def PI(prec) + raise ArgumentError, "Zero or negative argument 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 + w = 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 + w = 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 + + # Computes e (the base of natural logarithms) to the specified number of + # digits of precision. + def E(prec) + raise ArgumentError, "Zero or negative precision for E" if prec <= 0 + n = prec + BigDecimal.double_fig + one = BigDecimal("1") + y = one + d = y + z = one + i = 0 + while d.nonzero? && ((m = n - (y.exponent - d.exponent).abs) > 0) + m = BigDecimal.double_fig if m < BigDecimal.double_fig + i += 1 + z *= i + d = one.div(z,m) + y += d + end + y + end +end diff --git a/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/newton.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/newton.rb new file mode 100644 index 0000000000..59ac0f7f04 --- /dev/null +++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/newton.rb @@ -0,0 +1,77 @@ +# +# 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 1.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. +# +require "bigdecimal/ludcmp" +require "bigdecimal/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 + 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/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/util.rb b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/util.rb new file mode 100644 index 0000000000..09e926acd5 --- /dev/null +++ b/ruby_1_8_6/ext/bigdecimal/lib/bigdecimal/util.rb @@ -0,0 +1,65 @@ +# +# BigDecimal utility library. +# +# To use these functions, require 'bigdecimal/util' +# +# The following methods are provided to convert other types to BigDecimals: +# +# String#to_d -> BigDecimal +# Float#to_d -> BigDecimal +# Rational#to_d -> BigDecimal +# +# The following method is provided to convert BigDecimals to other types: +# +# BigDecimal#to_r -> Rational +# +# ---------------------------------------------------------------------- +# +class Float < Numeric + def to_d + BigDecimal(self.to_s) + end +end + +class String + def to_d + BigDecimal(self) + end +end + +class BigDecimal < Numeric + # Converts a BigDecimal to a String of the form "nnnnnn.mmm". + # This method is deprecated; use BigDecimal#to_s("F") instead. + def to_digits + if self.nan? || self.infinite? || self.zero? + self.to_s + else + i = self.to_i.to_s + s,f,y,z = self.frac.split + i + "." + ("0"*(-z)) + f + end + end + + # Converts a BigDecimal to a Rational. + def to_r + sign,digits,base,power = self.split + numerator = sign*digits.to_i + denomi_power = power - digits.size # base is always 10 + if denomi_power < 0 + Rational(numerator,base ** (-denomi_power)) + else + Rational(numerator * (base ** denomi_power),1) + end + end +end + +class Rational < Numeric + # Converts a Rational to a BigDecimal + 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 |