diff options
Diffstat (limited to 'ruby_1_8_5/ext/bigdecimal/lib/bigdecimal')
| -rw-r--r-- | ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/jacobian.rb | 85 | ||||
| -rw-r--r-- | ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/ludcmp.rb | 84 | ||||
| -rw-r--r-- | ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/math.rb | 235 | ||||
| -rw-r--r-- | ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/newton.rb | 77 | ||||
| -rw-r--r-- | ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/util.rb | 65 |
5 files changed, 0 insertions, 546 deletions
diff --git a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/jacobian.rb b/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/jacobian.rb deleted file mode 100644 index d80eeab901..0000000000 --- a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/jacobian.rb +++ /dev/null @@ -1,85 +0,0 @@ -# -# 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_5/ext/bigdecimal/lib/bigdecimal/ludcmp.rb b/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/ludcmp.rb deleted file mode 100644 index 8f4888725e..0000000000 --- a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/ludcmp.rb +++ /dev/null @@ -1,84 +0,0 @@ -# -# 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_5/ext/bigdecimal/lib/bigdecimal/math.rb b/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/math.rb deleted file mode 100644 index f3248a3c5c..0000000000 --- a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/math.rb +++ /dev/null @@ -1,235 +0,0 @@ -# -#-- -# 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_5/ext/bigdecimal/lib/bigdecimal/newton.rb b/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/newton.rb deleted file mode 100644 index 59ac0f7f04..0000000000 --- a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/newton.rb +++ /dev/null @@ -1,77 +0,0 @@ -# -# 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_5/ext/bigdecimal/lib/bigdecimal/util.rb b/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/util.rb deleted file mode 100644 index 09e926acd5..0000000000 --- a/ruby_1_8_5/ext/bigdecimal/lib/bigdecimal/util.rb +++ /dev/null @@ -1,65 +0,0 @@ -# -# 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 |