summaryrefslogtreecommitdiff
path: root/ext/bigdecimal/lib
diff options
context:
space:
mode:
authorshigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-07-18 15:24:25 +0000
committershigek <shigek@b2dd03c8-39d4-4d8f-98ff-823fe69b080e>2003-07-18 15:24:25 +0000
commite242cf9defb5442ef223535abe93399352cf139e (patch)
treeb8eb591e57c326fd4e7cb49365f7e51783179237 /ext/bigdecimal/lib
parentb10272dc371a03844f685dc6db765c2bea11d29b (diff)
More pathes from Tadasi Saito.
As discussed in ruby-dev ML: E,PI, etc are disabled. BigDecimal op String disabled. to_f changed. lib directory moved. git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@4092 b2dd03c8-39d4-4d8f-98ff-823fe69b080e
Diffstat (limited to 'ext/bigdecimal/lib')
-rw-r--r--ext/bigdecimal/lib/bigdecimal-rational.rb30
-rw-r--r--ext/bigdecimal/lib/jacobian.rb63
-rw-r--r--ext/bigdecimal/lib/ludcmp.rb75
-rw-r--r--ext/bigdecimal/lib/newton.rb75
4 files changed, 0 insertions, 243 deletions
diff --git a/ext/bigdecimal/lib/bigdecimal-rational.rb b/ext/bigdecimal/lib/bigdecimal-rational.rb
deleted file mode 100644
index 041e8c7..0000000
--- a/ext/bigdecimal/lib/bigdecimal-rational.rb
+++ /dev/null
@@ -1,30 +0,0 @@
-#
-# BigDecimal <-> Rational
-#
-class BigDecimal < Numeric
- # Convert BigDecimal to 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
- denominator = base ** (-denomi_power)
- else
- denominator = base ** denomi_power
- end
- Rational(numerator,denominator)
- end
-end
-
-class Rational < Numeric
- # 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/jacobian.rb b/ext/bigdecimal/lib/jacobian.rb
deleted file mode 100644
index 34a60ae..0000000
--- a/ext/bigdecimal/lib/jacobian.rb
+++ /dev/null
@@ -1,63 +0,0 @@
-#
-# 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/ludcmp.rb b/ext/bigdecimal/lib/ludcmp.rb
deleted file mode 100644
index c36f0de..0000000
--- a/ext/bigdecimal/lib/ludcmp.rb
+++ /dev/null
@@ -1,75 +0,0 @@
-#
-# 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
deleted file mode 100644
index 67a9247..0000000
--- a/ext/bigdecimal/lib/newton.rb
+++ /dev/null
@@ -1,75 +0,0 @@
-#
-# 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