eaiovnaovbqoebvqoeavibavo lib/bigdecimal/newton.rb000064400000003435147634530120011243 0ustar00require "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) 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 lib/bigdecimal/jacobian.rb000064400000004073147634530120011476 0ustar00# # 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 s = f.zero 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 lib/bigdecimal/util.rb000064400000004376147634530120010713 0ustar00class Integer < Numeric # call-seq: # int.to_d -> bigdecimal # # Convert +int+ to a BigDecimal and return it. # # require 'bigdecimal' # require 'bigdecimal/util' # # 42.to_d # # => # # def to_d BigDecimal(self) end end 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 # # => # # def to_d(precision=nil) BigDecimal(self, precision || Float::DIG+1) end end 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 # # => # # def to_d BigDecimal(self) end end 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 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) # # => # def to_d(precision) if precision <= 0 raise ArgumentError, "negative precision" end num = self.numerator BigDecimal(num).div(self.denominator, precision) end end lib/bigdecimal/ludcmp.rb000064400000004142147634530120011211 0ustar00require '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 lib/bigdecimal/math.rb000064400000011742147634530120010662 0ustar00require '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" # require "bigdecimal/math" # # include BigMath # # a = BigDecimal((PI(100)/2).to_s) # puts sin(a,100) # -> 0.10000000000000000000......E1 # module BigMath module_function # 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") 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 # 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") 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 # Computes the arctangent of x to the specified number of digits of precision. # # If x is NaN, returns NaN. 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 # 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