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

chudnovsky.rb

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]
  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

class Integer
  def self.isqrt(n)
    return Integer.sqrt(n) if n.bit_length < 512
    shift = [n.bit_length / 4, 1].max
    x = Integer.isqrt(n >> (shift * 2))
    x = (x << (shift - 1)) + (n >> (shift + 1)) / x
    xx = x * x
    while xx > n
      xx -= 2 * x - 1
      x -= 1
    end
    x
  end
end

def int_pow_large(x, n)
  return x**n if n < 10000
  a = int_pow_large(x, n/2)
  n.even? ? a * a : a * a * x
end

def pi(n)
  ten_pow_n = int_pow_large(10, n)
  _, q, t = series(0, (n + 13) / 14)
  a = (4270934400 * q)
  b = (Integer.isqrt(10005 * int_pow_large(10, 2*n)) * (13591409 * q + t))
  m = ((b.bit_length - a.bit_length) / Math.log2(10)).ceil
  int_pow_large(10, n + m) * a / b
end

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