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