Sample code without deprecated modules (#480)

* sample/pi Remove shebang and `include BigMath`

* Rewrite sample/linear.rb using gem matrix instead of bigdecimal/ludcmp

* Implement sample/nlsolve.rb without deprecated bigdecimal/newton

Author: tompng

sample/linear.rb

@@ -1,4 +1,3 @@
-#!/usr/local/bin/ruby
 # frozen_string_literal: false
 
 #
@@ -13,62 +12,101 @@
 
 # :stopdoc:
 require "bigdecimal"
-require "bigdecimal/ludcmp"
 
-#
-# NOTE:
-#   Change following BigDecimal.limit() if needed.
-BigDecimal.limit(100)
-#
+# Requires gem matrix
+require "matrix"
+
+class PrecisionSpecifiedValue
+  # NOTE:
+  #   Change following PREC if needed.
+
+  attr_reader :value
+  def initialize(value, prec)
+    @value = BigDecimal(value)
+    @prec = prec
+  end
+
+  def unwrap(value)
+    PrecisionSpecifiedValue === value ? value.value : value
+  end
+
+  def coerce(other)
+    [self.class.new(unwrap(other), @prec), self]
+  end
+
+  def abs
+    self.class.new(@value.abs, @prec)
+  end
+
+  def >(other)
+    @value > unwrap(other)
+  end
+
+  def <(other)
+    @value < unwrap(other)
+  end
+
+  def -(other)
+    self.class.new(@value.sub(unwrap(other), @prec), @prec)
+  end
+
+  def +(other)
+    self.class.new(@value.add(unwrap(other), @prec), @prec)
+  end
+
+  def *(other)
+    self.class.new(@value.mult(unwrap(other), @prec), @prec)
+  end
+
+  def quo(other)
+    self.class.new(@value.div(unwrap(other), @prec), @prec)
+  end
+end
+
+return if __FILE__ != $0
 
-include LUSolve
 def rd_order(na)
-   printf("Number of equations ?") if(na <= 0)
-   n = ARGF.gets().to_i
+  printf("Number of equations ?") if(na <= 0)
+  ARGF.gets().to_i
 end
 
-na   = ARGV.size
-zero = BigDecimal("0.0")
-one  = BigDecimal("1.0")
+na = ARGV.size
 
 while (n=rd_order(na))>0
   a = []
-  as= []
   b = []
   if na <= 0
      # Read data from console.
      printf("\nEnter coefficient matrix element A[i,j]\n")
      for i in 0...n do
-       for j in 0...n do
+       a << n.times.map do |j|
          printf("A[%d,%d]? ",i,j); s = ARGF.gets
-         a  << BigDecimal(s)
-         as << BigDecimal(s)
+         BigDecimal(s)
        end
        printf("Contatant vector element b[%d] ? ",i)
        b << BigDecimal(ARGF.gets)
      end
   else
-     # Read data from specified file.
-     printf("Coefficient matrix and constant vector.\n")
-     for i in 0...n do
-       s = ARGF.gets
-       printf("%d) %s",i,s)
-       s = s.split
-       for j in 0...n do
-         a  << BigDecimal(s[j])
-         as << BigDecimal(s[j])
-       end
-       b << BigDecimal(s[n])
-     end
+    # Read data from specified file.
+    printf("Coefficient matrix and constant vector.\n")
+    for i in 0...n do
+      s = ARGF.gets
+      printf("%d) %s",i,s)
+      s = s.split
+      a << n.times.map {|j| BigDecimal(s[j]) }
+      b << BigDecimal(s[n])
+    end
   end
-  x = lusolve(a,b,ludecomp(a,n,zero,one),zero)
+
+  prec = 100
+  matrix = Matrix[*a.map {|row| row.map {|v| PrecisionSpecifiedValue.new(v, prec) } }]
+  vector = b.map {|v| PrecisionSpecifiedValue.new(v, prec) }
+  x = matrix.lup.solve(vector).map(&:value)
+
   printf("Answer(x[i] & (A*x-b)[i]) follows\n")
   for i in 0...n do
      printf("x[%d]=%s ",i,x[i].to_s)
-     s = zero
-     for j in 0...n do
-       s = s + as[i*n+j]*x[j]
-     end
-     printf(" & %s\n",(s-b[i]).to_s)
+     diff = a[i].zip(x).sum {|aij, xj| aij*xj }.sub(b[i], 10)
+     printf(" & %s\n", diff.to_s)
   end
 end

sample/nlsolve.rb

@@ -1,4 +1,3 @@
-#!/usr/local/bin/ruby
 # frozen_string_literal: false
 
 #
@@ -7,34 +6,54 @@
 #
 
 require "bigdecimal"
-require "bigdecimal/newton"
-include Newton
-
-class Function # :nodoc: all
-  def initialize()
-    @zero = BigDecimal("0.0")
-    @one  = BigDecimal("1.0")
-    @two  = BigDecimal("2.0")
-    @ten  = BigDecimal("10.0")
-    @eps  = BigDecimal("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
+require_relative "linear"
+
+# Requires gem matrix
+require "matrix"
+
+# :stopdoc:
+
+def func((x, y)) # defines functions solved
+  [
+    x**2 + y**2 - 2,
+    (x - 1)**2 + (y + 1)**2 - 3
+  ]
+end
+
+def jacobian(x, f, delta, prec)
+  dim = x.size
+  dim.times.map do |i|
+    xplus = Array.new(dim) {|j| x[i] + (j == i ? delta : 0) }
+    xminus = Array.new(dim) {|j| x[i] - (j == i ? delta : 0) }
+    yplus = f.call(xplus)
+    yminus = f.call(xminus)
+    yplus.zip(yminus).map {|p, m| (p - m).div(2 * delta, prec) }
+  end.transpose
+end
+
+def nlsolve(initial_x, prec:, max_iteration: 100, &f)
+  initial_x = initial_x.map {|v| BigDecimal(v) }
+  x = initial_x
+  delta = BigDecimal(0.01)
+  calc_prec = prec + 10
+  max_iteration.times do |iteration|
+    # Newton step
+    jacobian = jacobian(x, f, delta, calc_prec)
+    matrix = Matrix[*jacobian.map {|row| row.map {|v| PrecisionSpecifiedValue.new(v, calc_prec) } }]
+    y = f.call(x)
+    vector = y.map {|v| PrecisionSpecifiedValue.new(v, calc_prec) }
+    dx = matrix.lup.solve(vector).map(&:value)
+    x_prev = x
+    x = x.zip(dx).map {|xi, di| xi.sub(di, prec) }
+    movement = x_prev.zip(x).map {|xn, xi| (xn - xi).abs }.max
+    delta = [movement, delta].min.mult(1, 10)
+    break if movement.zero? || movement.exponent < -prec
   end
+  x
 end
 
-f = BigDecimal.limit(100)
-f = Function.new
-x = [f.zero,f.zero]      # Initial values
-n = nlsolve(f,x)
-p x
+initial_value = [1, 1]
+ans = nlsolve(initial_value, prec: 100) {|x| func(x) }
+diff = func(ans).map {|v| v.mult(1, 10) }
+p(ans:)
+p(diff:)

sample/pi.rb

@@ -1,4 +1,3 @@
-#!/usr/local/bin/ruby
 # frozen_string_literal: false
 
 #
@@ -11,11 +10,9 @@
 require "bigdecimal"
 require "bigdecimal/math.rb"
 
-include BigMath
-
 if ARGV.size == 1
     print "PI("+ARGV[0]+"):\n"
-    p PI(ARGV[0].to_i)
+    p BigMath.PI(ARGV[0].to_i)
 else
     print "TRY: ruby pi.rb 1000 \n"
 end