pi.rb

require "bigdecimal"
require "bigdecimal/math"

def measure(name, &block)
  s = Process.clock_gettime(Process::CLOCK_MONOTONIC)
  r = block.call
  t = Process.clock_gettime(Process::CLOCK_MONOTONIC)

  time = t - s
  STDERR.puts "#{name}: #{time} s"

  r
end

DIGITS_PER_COLUMN = 10
DIGITS_PER_LINE = 100
DIGITS_PER_BLOCK = 1000

A, B, C, D, E = 13591409, 545140134, 640320, 426880, 10005
C3_24 = C * C * C / 24
DIGITS_PER_TERM = 14.1816474627254776555 # log(53360^3) / log(10)

def compute_pqt(n1, n2)
  if n1 + 1 == n2
    p = 2 * n2 - 1
    p *= 6 * n2 - 1
    p *= 6 * n2 - 5
    q = C3_24 * n2 * n2 * n2
    t = (A + B * n2) * p
    t = -t if n2.odd?
  else
    m = (n1 + n2) / 2
    p1, q1, t1 = compute_pqt(n1, m)
    p2, q2, t2 = compute_pqt(m, n2)
    p = p1 * p2
    q = q1 * q2
    t = t1 * q2 + p1 * t2
  end
  return p, q, t
end

target_prec = ARGV[0].to_i

prec = target_prec + 15

n = (prec / DIGITS_PER_TERM).to_i
p, q, t = measure("compute_pqt") { compute_pqt(0, n) }
pi = measure("pi") do
  pi = BigDecimal(D, prec) * BigMath.sqrt(BigDecimal(E, prec), prec) * BigDecimal(q, prec)
  pi /= BigDecimal(A, prec) * BigDecimal(q, prec) + BigDecimal(t, prec)
  pi
end

s = measure("to_s") { pi.to_s("F") }

puts "π = #{s[0...2]}"
index = 2
(target_prec.to_f / DIGITS_PER_BLOCK).ceil.times do |i|
  (DIGITS_PER_BLOCK / DIGITS_PER_LINE).times do |j|
    break if index >= target_prec + 2
    (DIGITS_PER_LINE / DIGITS_PER_COLUMN).times do |k|
      break if index >= target_prec + 2
      print ' '
      next_index = [index + DIGITS_PER_COLUMN, s.size].min
      print s[index...next_index]
      index = next_index
    end
    puts
  end
  puts unless index >= target_prec + 2
end