piChudnovsky.frink

/** Program to calculate pi to a large number of digits using the Chudnovsky
algorithm with binary splitting.

see:
http://numbers.computation.free.fr/Constants/Algorithms/splitting.html

This program is a testbed for timing the performance of this algorithm and
prints extraneous output.  If you want a quieter, cached version that you
can use in your own library, see pi.frink instead. */

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

println["Calculating \$digits digits."]
//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: " + format[end-start1, "s", 3]];

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

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

println["Calculation complete."]

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

start = now[]
println[value]
end = now[]
println["Time spent in printing: " + format[end-start, "s", 3]]

end = now[]

println["Total time is " + format[end-start1, "s", 3]]

//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
}

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 18322 days, 0 hours, 33 minutes ago.