Download or view pi2.frink in plain text format
/** Program to calculate pi using 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
This differs from the basic pi algorithm because when you ask for more
digits of pi than you previously had, this *resumes* the calculation
rather than starting from scratch. It also contains a resumable
sqrt[640320] to extend the results.
*/
class Pi
{
class var digitsPerIteration = 14.1816474627254776555
class var largestDigits = -1
class var cachePi = undef
/** These variables help us resume the state of the binary-splitting
algorithm. */
class var largestP = 1
class var largestQ = 1
class var largestK = 0
class var largestG = 1
/** The current best estimate of sqrt[640320] */
class var sqrt640320 = 800.199975006248
/** Current working precision of sqrt in digits */
class var sqrtPrecision = 14
/** 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-1 >= 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]
// println["Calculating [$largestK, $k]"]
if (largestK < k) // We may have already calculated enough!
{
if (largestK == 0) // This is the first iteration
{
p = p[0,k]
q = q[0,k]
g = g[0,k]
} else
{
// Continuing iterations
pmb = p[largestK, k]
// p[a,b] = p[a,m] p[m,b]
p = largestP * pmb
// We need to cache g[0,largestK]
g = largestG * g[largestK, k]
// We need a special rule here
// q[a,m] p[m,b] + q[m,b] g[a,m]
q = largestQ * pmb + q[largestK, k] * largestG
}
// Store the biggest values back in the cache
largestP = p
largestQ = q
largestG = g
largestK = k
}
// Do a resumable square root of 640320 (it may not need to do any
// more work)
sqC = sqrt640320[digits+8]
// At this point, largestP contains p and largestQ contains q
// for this number.
piprime = largestP * 53360. / (largestQ + 13591409 * largestP)
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]
}
/** Resumable sqrt of 640320. If we have already calculated sqrt[640320]
to the required number of digits, this will return quickly, otherwise,
it will resume the calculation efficiently and provide more digits. */
class sqrt640320[digits] :=
{
t = sqrt640320
// println["Resuming sqrt with precision $sqrtPrecision"]
origPrecision = getPrecision[]
try
{
while (sqrtPrecision <= 2 * digits)
{
setPrecision[1 + 2 * sqrtPrecision]
t = (640320/t + t) / 2.
sqrtPrecision = sqrtPrecision * 2
}
// println["Bailing with next precision at $sqrtPrecision"]
sqrt640320 = t
setPrecision[digits]
return 1. * sqrt640320
}
finally
setPrecision[origPrecision]
}
}
Download or view pi2.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, eliasen@mindspring.com