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

cal_pi_arctan_2.rb

#! /usr/local/bin/ruby
#*********************************************
# 円周率計算 by Arctan 系の公式
# ( 各項の Arctan を個別に計算後に加減算する方法 )
#
#    1: Machin
#    2: Klingenstierna
#    3: Euler
#    4: Euler(2)
#    5: Gauss
#    6: Stormer
#    7: Stormer(2)
#    8: Takano
#*********************************************
#
class CalcPiArctan
  # 公式
  FORMULA = [
    "Machin", "Klingenstierna", "Euler",    "Euler2",
    "Gauss",  "Stormer",        "Stormer2", "Takano"
  ]  # 公式名
  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)
    [ 48, 1, 18,  32, 1,  57, -20, 1, 239               ],  #   5: Gauss
    [ 24, 1,  8,   8, 1,  57,   4, 1, 239               ],  #   6: Stormer
    [176, 1, 57,  28, 1, 239, -48, 1, 682, 96, 1,  12943],  #   7: Stormer(2)
    [ 48, 1, 49, 128, 1,  57, -20, 1, 239, 48, 1, 110443]   #   8: Takano
  ]
  # 保存ファイル名
  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 = Array.new
    0.upto(@ks - 1) do |i|              # 計算項数
      n  = @l / (Math::log10(@kk[i * 3 + 2]) - Math::log10(@kk[i * 3 + 1])) + 1
      @n << (n / 2).truncate + 1
    end
  end

  # 計算・結果出力
  def calc
    # ==== 計算開始時刻
    t0 = Time.now
    # ==== 配列宣言・初期化
    a  = Array.new(@ks, nil)    # 各項の(2k-1)部分を除いた部分
    a.each_index {|i| a[i] = Array.new(@l1 + 2, 0)}
    q  = Array.new(@ks, nil)    # 各項の(2k-1)部分を含む部分
    q.each_index {|i| q[i] = Array.new(@l1 + 2, 0)}
    sk = Array.new(@ks, nil)    # 各項の総和
    sk.each_index {|i| sk[i] = Array.new(@l1 + 2, 0)}
    s  = 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
    # ==== 計算
    0.upto(@ks - 1) do |i|
      1.upto(@n[i]) do |k|
        # 各項の分数の部分
        a[i] = long_div(a[i], @kk[i * 3 + 2])
        a[i] = long_div(a[i], @kk[i * 3 + 2])
        unless @kk[i * 3 + 1] == 1  # 分子が 1 でない時
          a[i] = long_mul(a[i], @kk[i * 3 + 1])
          a[i] = long_mul(a[i], @kk[i * 3 + 1])
        end
        # 各項の 1 / (2 * k - 1) の部分
        q[i] = long_div(a[i], 2 * k - 1)
        # 各項の総和に加減算
        sk[i] = k % 2 == 0 ? long_sub(sk[i], q[i]) : long_add(sk[i], q[i])
      end
    end
    # 各項の総和の加減算
    0.upto(@ks - 1) do |i|
      s = (@kk[i * 3] >= 0 ? long_add(s, sk[i]) : long_sub(s, sk[i]))
    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)\n"
    print "5:Gauss, 6:Stormer, 7:Stormer(2), 8:Takano : "
    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