bbp.cr

require "math"
macro dcl(*args)
  {% aa = [] of (Path|Call) %}
  {% for x in args %}
    {% if x.is_a?(TypeDeclaration) %}
      {% aa.push(x.var) %}
      {% t = x.type %}
    {% else %}
      {% aa.push(x) %}
    {% end %}
  {% end %}
  {% for v in aa %}
    {{v.id}} : {{t.id}}
  {% end %}
##  {% debug %}
end

def mul_mod(a : Int64, b : Int64, m : Int64)
  (a * b) % m
end

def inv_mod(x : Int64, y : Int64)
   dcl q, u, v, a, c, t : Int64

   u = x
   v = y
   c = 1
   a = 0
   loop do
     q = (v / u).to_i64
     t = c
     c = a - q * c
     a = t

     t = u
     u = v - q * u
     v = t
     break if u == 0
   end
   a %= y
   a += y if a < 0
   a
end

def pow_mod(a : Int64, b : Int64, m : Int64)
  dcl r, aa : Int64

  r = 1
  aa = a
  loop do
    if b & 1 > 0
      r = mul_mod(r, aa, m)
    end
    b >>= 1

    break if b == 0
    aa = mul_mod(aa, aa, m)
  end
  r
end

def is_prime(n : Int64)
    return false if n & 1 == 0

    i = 3_i64

    while i <= n
      res, rem = n.divmod(i)
      return !(rem == 0) if i > res || rem == 0
      i += 2
    end
    true
end

def next_prime(n : Int64)
  cn = n + 1 + n & 1
  until is_prime(cn)
    cn += 1
  end
  cn
end

def main
  dcl av, a, vmax, nbig, n, num, den, k, kq, kq2, t, v, s, i : Int64
  sum : Float64
  n = ARGV[0]?.try &.to_i64 || raise "This program computes the n'th decimal digit of PI\nUsage: PROGRAM n, where \"n\" is the digit you want"
  nbig = ((n + 20) * Math.log(10) / Math.log(2)).to_i64

  sum = 0
  a = 3
  while a <= (2 * nbig)
    vmax = (Math.log(2 * nbig) / Math.log(a)).to_i64;
    av = a ** vmax

    s = 0
    num = 1
    den = 1
    v = 0
    kq = 1
    kq2 = 1
    (1_i64..nbig).each do |k|
        t = k
        if kq >= a
          loop do
            t = (t / a).to_i64
            v -= 1
            break unless (t % a) == 0
          end
          kq = 0
        end
        kq += 1
        num = mul_mod(num, t, av)

        t = (2_i64 * k - 1_i64);
        if kq2 >= a
            if kq2 == a
              loop do
                t = (t / a).to_i64
                v += 1
                break unless (t % a) == 0
              end
            end
            kq2 -= a
        end
        den = mul_mod(den, t, av)
        kq2 += 2

        if v > 0
            t = inv_mod(den, av)
            t = mul_mod(t, num, av)
            t = mul_mod(t, k, av)
            (v..vmax-1).each do
              t = mul_mod(t, a, av)
            end
            s += t
            s -= av if s >= av
        end
    end

    t = pow_mod(10, n - 1, av)
    s = mul_mod(s, t, av)
    sum = (sum + s.to_f64 / av.to_f64) % 1.0_f64

    a = next_prime(a)
  end # while a <= (2 * nbig)

  printf("Decimal digits of pi at position %d: %09d\n", n, (sum * 1e+9).to_i64 )
end

main