Chudnovskyの公式で円周率を計算するプログラム

chudnovsky.rb

require "gmp"

def series(l, r)
  if r - l <= 1
    p = (2 * r - 1) * (6 * r - 5) * (6 * r - 1)
    q = r * r * r * 10939058860032000
    a = (r % 2 == 0 ? 1 : -1) * (13591409 + 545140134 * r)

    [p, q, a * p].map { GMP.Z(_1) }
  else
    m = (l + r) / 2
    p0, q0, t0 = series(l, m)
    p1, q1, t1 = series(m, r)

    [p0 * p1, q0 * q1, t0 * q1 + p0 * t1]
  end
end

def pi(n)
  _, q, t = series(0, (n + 13) / 14)
  GMP.F(10005, Math.log2(10).*(n).ceil).rec_sqrt * 4270934400 * q / (13591409 * q + t)
end

if $0 == __FILE__
  if ARGV.size == 1
    p pi(ARGV[0].to_i)
  else
    puts "TRY: ruby chudnovsky.rb 1000"
  end
end