pi2.frink

View or download 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]
   }
}


View or download 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 18492 days, 22 hours, 26 minutes ago.