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 was born 19963 days, 19 hours, 19 minutes ago.