piChudnovskyNew.frink

Download or view piChudnovskyNew.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.

    This is a new version that's intended to be used with Frink: The Next
    Generation.  It uses mostly fast floating-point algorithms.

    See:  https://frinklang.org/experimental.html
*/


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

//use sqrtWayne.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["Time spent in binary splitting: " + format[end-start1, "s", 3]];

start = now[];
sqC = sqrt[640320, digits+2]
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]]

q[a,b] :=
{
   if (b-a)==1
      return (-1)^b * g[a,b] * (13591409 + 545140134 b)

   m = (a+b) div 2
   return q[a,m] p[m,b] + q[m,b] g[a,m]
}

p[a,b] :=
{
   if (b-a) == 1
      return 10939058860032000 b^3  // Constant is C^3/24

   m = (a+b) div 2
   return p[a,m] p[m,b]
}

g[a,b] :=
{
   if (b-a) == 1
      return (6b-5)(2b-1)(6b-1)

   m = (a+b) div 2
   return g[a,m] g[m,b]
}


Download or view piChudnovskyNew.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 19944 days, 20 hours, 15 minutes ago.