/** Program to calculate pi to arbitrary precision. This algorithm uses binary splitting. This version has been updated to be fastest with Frink: The Next Generation. Due to limitations in Java's BigDecimal class, this is limited to a little over 101 million digits. You usually use this by calling Pi.getPi[digits] which will return pi to the number of digits you specify. You can also call Pi.get2Pi[digits] to return 2 pi. See http://numbers.computation.free.fr/Constants/Algorithms/splitting.html */ use sqrtWayne.frink class Pi { class var digitsPerIteration = 14.1816474627254776555 class var largestDigits = -1 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 >= 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 = sqrt[640320, digits+8] piprime = p * 53360. / (q + 13591409 * p) piFull = piprime * sqC // Truncate to the desired number of digits setPrecision[digits] pi = 1. * piFull 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] } }