# 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 } ```