As discussed in ruby-dev ML:

  lib directory moved.
  util.rb created instead of bigdecimal-rational.rb


git-svn-id: svn+ssh://ci.ruby-lang.org/ruby/trunk@4093 b2dd03c8-39d4-4d8f-98ff-823fe69b080e

5 files changed, 325 insertions, 0 deletions

diff --git a/ext/bigdecimal/lib/bigdecimal/jacobian.rb b/ext/bigdecimal/lib/bigdecimal/jacobian.rb
new file mode 100644
index 0000000000..34a60ae67a
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/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/bigdecimal/ludcmp.rb b/ext/bigdecimal/lib/bigdecimal/ludcmp.rb
new file mode 100644
index 0000000000..c36f0dea5b
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/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/bigdecimal/newton.rb b/ext/bigdecimal/lib/bigdecimal/newton.rb
new file mode 100644
index 0000000000..67a92474ac
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/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/bigdecimal/nlsolve.rb b/ext/bigdecimal/lib/bigdecimal/nlsolve.rb
new file mode 100644
index 0000000000..08f17f9ecd
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/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/bigdecimal/util.rb b/ext/bigdecimal/lib/bigdecimal/util.rb
new file mode 100644
index 0000000000..d9b7ae3131
--- /dev/null
+++ b/ext/bigdecimal/lib/bigdecimal/util.rb
@@ -0,0 +1,74 @@
+#
+# BigDecimal utility library.
+# ----------------------------------------------------------------------
+# Contents:
+#
+#   String#
+#     to_d      ... to BigDecimal
+#
+#   Float#
+#     to_d      ... to BigDecimal
+#
+#   BigDecimal#
+#     to_digits ... to xxxxx.yyyy form digit string(not 0.zzzE?? form).
+#     to_r      ... to Rational
+#
+#   Rational#
+#     to_d      ... to BigDecimal
+#
+# ----------------------------------------------------------------------
+#
+class Float < Numeric
+  def to_d
+    BigFloat::new(selt.to_s)
+  end
+end
+
+class String
+  def to_d
+    BigDecimal::new(self)
+  end
+end
+
+class BigDecimal < Numeric
+  # to "nnnnnn.mmm" form digit string
+  def to_digits
+     if self.nan? || self.infinite?
+        self.to_s
+     else
+       s,i,y,z = self.fix.split
+       s,f,y,z = self.frac.split
+       if s > 0
+         s = ""
+       else
+         s = "-"
+       end
+       s + i + "." + f
+     end
+  end
+
+  # 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