piChudnovsky.frink

View or download piChudnovsky.frink in plain text format


// Program to calculate pi using binary splitting.
// See http://numbers.computation.free.fr/Constants/Algorithms/splitting.html

//input["Hit Enter. "]    // This is so we can attach jvirtualvm profiler.

use karatsuba.frink

digits = million

if length[ARGS] >= 1
   digits = eval[ARGS@0]

digitsPerIteration = 14.1816474627254776555

// Find number of terms to calculate
k = floor[digits/digitsPerIteration]

setPrecision[digits+3]

sum = 0

println["k=$k"]

start1 = now[];

// pi = p[0,k] * (C/D) * sqrt[C] / (q[0,k] + A * p)
p = p[0,k]
q = q[0,k]
end = now[];

//println["p is $p: " + factor[p]]
//println["q is $q: " /* + factor[abs[q]]*/]

println["Time spent in binary splitting: " + (end-start1)];

start = now[];
sqC = karatsubaSqrt[640320, digits+2, true]
end = now[];
println["Time spent in square root: " + (end-start)];

start = now[];
piprime = p * 53360 / (q + 13591409 * p)
pi = piprime * sqC
end = now[];
println["Time spent in combining operations " + (end-start)]


println["Calculation complete."]

setPrecision[digits]
start = now[]
value = 1. * pi
end = now[]
println["Time spent in floating-point calculation: " + (end-start)]

start = now[]
println[value]
end = now[]
println["Time spent in printing: " + (end-start)]

end = now[]

println["Total time is " + (end-start1)]

//sqrtPi = sqrt[pi, digits]
//println["Square root of pi is:"]
//println[sqrtPi]
//println[sqrtPi^2 - pi]


q[a,b] :=
{
   if (b-a)==1
   {
      result = (-1)^b * g[a,b] * (13591409 + 545140134 b)
//      println["q[$a,$b]=" + factor[abs[result]]]
      return result
   }

   m = (a+b) div 2
//   println["m=$m"]
   r = q[a,m] p[m,b] + q[m,b] g[a,m]
//   println["q[$a,$b]=" + factor[abs[r]]]
//   println["r=$r"]
   return r
}

p[a,b] :=
{
   if (b-a) == 1
   {
      result = 10939058860032000 b^3  // Constant is C^3/24
//      println["p[$a,$b]=" + factor[result]]
      return result
   }
   m = (a+b) div 2
   result = p[a,m] p[m,b]
//   println["p[$a,$b]=" + factor[result]]
   return result
}

g[a,b] :=
{
   if (b-a) == 1
   {
      result = (6b-5)(2b-1)(6b-1)
//      println["g[$a,$b]=" + factor[result]]
      return result
   }
   m = (a+b) div 2
   result = g[a,m] g[m,b]
//   println["g[$a,$b]=" + factor[result]]
   return result
}


View or download piChudnovsky.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17683 days, 10 hours, 54 minutes ago.