Download or view piChudnovsky.frink in plain text format
/** 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
}
Download or view 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 20217 days, 23 hours, 17 minutes ago.