ガウス・ルジャンドルのアルゴリズムで円周率を求める

gauss.rb

require "gmp"

def pi(n)
  prec_bits = Math.log2(10).*(n).ceil

  an = GMP.F(1, prec_bits)
  bn = GMP.F(2, prec_bits).rec_sqrt
  tn = GMP.F(0.25, prec_bits)
  (0...Math.log2(n) - 1).each do |n|
    an1 = (an + bn) * GMP.F(0.5, 1)
    an, bn, tn = an1, an.*(bn).sqrt, tn - (an - an1) ** 2 * GMP.F(2, 1) ** n
  end

  (an + bn) ** 2 / (tn * 4)
end

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