Ruby script to compute Pi with arctan's formula.(old ver.)

calc_pi_arctan_old.rb

#! /usr/local/bin/ruby
#*********************************************
# 円周率計算 by Arctan 系の公式
#
#    1: Machin
#    2: Klingenstierna
#    3: Euler
#    4: Euler(2)
#*********************************************
#
class CalcPiArctan
  # 公式
  FORMULA = ["Machin", "Klingenstierna", "Euler", "Euler2"]  # 公式名
  KK = [                                    # 公式内係数
    [16, 1,  5,  -4, 1, 239             ],  #   1: Machin
    [32, 1, 10,  -4, 1, 239, -16, 1, 515],  #   2: Klingenstierna
    [20, 1,  7,   8, 3,  79             ],  #   3: Euler
    [16, 1,  5,  -4, 1,  70,   4, 1,  99]   #   4: Euler(2)
  ]
  # 保存ファイル名
  FILE_PRE = "PI_"   # プリフィックス
  FILE_EXT = ".txt"  # 拡張子
  # 出力文字列 ( コンソール、テキストファイル共通 )
  STR_TITLE   = "** Pi Computation by Arctan method **"
  STR_FORMULA = "   Formula = [ %s ]"
  STR_DIGITS  = "   Digits  = %d"
  STR_TIME    = "   Time    = %f seconds"
  # 多桁計算
  MAX_DIGITS = 100000000  # 1つの int で8桁扱う

  def initialize(x, y)
    # 結果出力ファイル名生成
    @fname = "#{FILE_PRE}#{FORMULA[x - 1]}#{FILE_EXT}"
    puts
    puts "[ Output file : #{@fname} ]"
    puts
    @formula = FORMULA[x-1]             # 公式名取得
    @ks = KK[x-1].size / 3              # 計算対象公式の項数
    @kk = KK[x-1]                       # 計算対象公式の係数
    @l  = y                             # 計算桁数
    @l1 = (@l / 8) + 1                  # 配列サイズ
    n   = @l / Math::log10(@kk[2]) + 1
    @n  = (n / 2).truncate + 1          # 計算項数
  end

  # 計算・結果出力
  def calc
    # ==== 計算開始時刻
    t0 = Time.now
    # ==== 配列宣言・初期化
    s  = Array.new(@l1 + 2, 0)  # 総和
    a  = Array.new(@ks, nil)    # 公式の各項
    a.each_index {|i| a[i] = Array.new(@l1 + 2, 0)}
    q  = Array.new(@l1 + 2, 0)  # 各項の和
    # ==== 計算初期値
    0.upto(@ks - 1) do |i|
      a[i][0] = @kk[i * 3] * @kk[i * 3 + 2]
      a[i][0] *= -1 if @kk[i * 3] < 0
      # 分子が 1 で無ければ、さらに分子の値で除算しておく
      a[i] = long_div(a[i], @kk[i * 3 + 1]) unless @kk[i * 3 + 1] == 1
    end
    # ==== 計算
    1.upto(@n) do |k|
      # 各項の分数計算 ( １つ前の計算結果に追加で除算・乗算 )
      0.upto(@ks - 1) do |i|
        # 本来 x * x で除算したい（した方が高速だ）が、
        # x が 8桁以下でも x * x が8桁超になるケースがあるため２回に分けている。
        a[i] = long_div(a[i], @kk[i * 3 + 2])
        a[i] = long_div(a[i], @kk[i * 3 + 2])
        # 分子が 1 でない時
        #   本来 x * x を乗算したい（した方が高速だ）が、
        #   x が 8桁以下でも x * x が8桁超になるケースがあるため２回に分けている。
        unless @kk[i * 3 + 1] == 1
          a[i] = long_mul(a[i], @kk[i * 3 + 1])
          a[i] = long_mul(a[i], @kk[i * 3 + 1])
        end
      end
      # 各項の加減算
      1.upto(@ks - 1) do |i|
        if @kk[i * 3] >= 0  # 加算
          q = i == 1 ? long_add(a[i - 1], a[i]) : long_add(q, a[i])
        else                # 減算
          q = i == 1 ? long_sub(a[i - 1], a[i]) : long_sub(q, a[i])
        end
      end
      # 1 / (2 * k - 1) の部分
      q = long_div(q, 2 * k - 1)
      # 総和計算
      s = (k % 2) == 0 ? long_sub(s, q) : long_add(s, q)
    end
    # ==== 計算終了時刻
    t1 = Time.now
    # ==== 計算時間
    tt = t1 - t0
    # ==== 結果出力
    display(tt, s)
  rescue => e
    raise
  end

  # ロング + ロング
  def long_add(a, b)
    z = Array.new(@l1 + 2, 0)
    cr = 0
    (@l1 + 1).downto(0) do |i|
      z[i] = a[i] + b[i] + cr
      if z[i] < MAX_DIGITS
        cr = 0
      else
        z[i] -= MAX_DIGITS
        cr = 1
      end
    end
    return z
  rescue => e
    raise
  end

  # ロング - ロング
  def long_sub(a, b)
    z = Array.new(@l1 + 2, 0)
    br = 0
    (@l1 + 1).downto(0) do |i|
      z[i] = a[i] - b[i] - br
      if z[i] >= 0
        br = 0
      else
        z[i] += MAX_DIGITS
        br = 1
      end
    end
    return z
  rescue => e
    raise
  end

  # ロング * ショート
  def long_mul(a, b)
    z = Array.new(@l1 + 2, 0)
    cr = 0
    (@l1 + 1).downto(0) do |i|
      w = a[i]
      z[i] = (w * b + cr) % MAX_DIGITS
      cr = (w * b + cr) / MAX_DIGITS
    end
    return z
  rescue => e
    raise
  end

  # ロング / ショート
  def long_div(a, b)
    z = Array.new(@l1 + 2, 0)
    r = 0
    0.upto(@l1 + 1) do |i|
      w = a[i]
      z[i] = (w + r) / b
      r = ((w + r) % b) * MAX_DIGITS
    end
    return z
  rescue => e
    raise
  end

  # 結果出力
  def display(tt, s)
    puts STR_TITLE
    puts STR_FORMULA % @formula
    puts STR_DIGITS % @l
    puts STR_TIME % tt
    # ファイル出力
    out_file = File.open(@fname, "w")
    out_file.puts STR_TITLE
    out_file.puts STR_FORMULA % @formula
    out_file.puts STR_DIGITS % @l
    out_file.puts STR_TIME % tt
    out_file.puts "\n          %d." % s[0]
    1.upto(@l1 - 1) do |i|
      out_file.printf("%08d:", (i - 1) * 8 + 1) if i % 10 == 1
      out_file.printf(" %08d", s[i])
      out_file.puts if i % 10 == 0
    end
  rescue => e
    raise
  end
end

if __FILE__ == $0
  begin
    # ==== 使用公式番号入力
    print "1:Machin, 2:Klingenstierna, 3:Euler, 4:Euler(2) : "
    f = gets.to_i
    f = 1 if f < 1 || f > 8  # 範囲外なら 1:Machin と判断
    # ==== 計算桁数入力
    print "Please input number of Pi Decimal-Digits : "
    n = gets.to_i
    # 計算クラスインスタンス化
    obj = CalcPiArctan.new(f, n)
    # 円周率計算
    obj.calc
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}\n"
    e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
  end
end