// Program to calculate pi using binary splitting. // See http://numbers.computation.free.fr/Constants/Algorithms/splitting.html use karatsuba.frink class Pi { class var digitsPerIteration = 14.1816474627254776555 class var largestDigits = undef class var cachePi = undef /** This is the main public method to get the value of pi to a certain number of digits, calculating it if need be. If pi has already been calculated to a sufficient number of digits, this returns it from the cache. */ class getPi[digits = getPrecision[]] := { origPrec = getPrecision[] try { setPrecision[digits] if (largestDigits != undef) and (largestDigits >= digits) return 1. * cachePi else return 1. * calcPi[digits] } finally setPrecision[origPrec] } /** This is the main public method to get the value of 2 * pi to a certain number of digits, calculating it if need be. If pi has already been calculated to a sufficient number of digits, this returns it from the cache. */ class get2Pi[digits = getPrecision[]] := { origPrec = getPrecision[] try { setPrecision[digits] if (largestDigits != undef) and (largestDigits >= digits) return 2 * cachePi else return 2 * calcPi[digits] } finally setPrecision[origPrec] } /** This is an internal method that calculates the digits of pi if necessary. */ class calcPi[digits] := { oldPrec = getPrecision[] // Find number of terms to calculate k = max[floor[digits/digitsPerIteration], 1] try { setPrecision[digits+5] sum = 0 p = p[0,k] q = q[0,k] sqC = karatsubaSqrt[640320, digits+8, true] piprime = p * 53360 / (q + 13591409 * p) piRational = piprime * sqC // Convert to floating-point setPrecision[digits] pi = 1. * piRational largestDigits = digits cachePi = pi return pi } finally setPrecision[oldPrec] } /** Internal method for binary splitting. */ class 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] } /** Internal method for binary splitting. */ class p[a,b] := { if (b-a) == 1 return 10939058860032000 b^3 m = (a+b) div 2 return p[a,m] p[m,b] } /** Internal method for binary splitting. */ class 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] } }