/**  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]
   }
}