/**  Program to calculate e to arbitrary precision.  This uses binary
     splitting for speed.

     You usually use this by calling E.getE[digits] which will return e
     to the number of digits you specify.

     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.

     This differs from the basic algorithm in eBinarySplitting.frink
     because when you ask for more digits of e than you previously had,
     this *resumes* the calculation rather than starting from scratch.

     See:
     http://numbers.computation.free.fr/Constants/Algorithms/splitting.html
     http://www.ginac.de/CLN/binsplit.pdf
*/
class E
{
   class var largestDigits = -1
   class var cacheE = undef

   /** These variables help us resume the state of the binary-splitting
       algorithm. */
   class var largestP = 1
   class var largestQ = 1
   class var largestK = 1
   class var logFactorial = 0.

   /** This is the main public method to get the value of e to a certain
       number of digits, calculating it if need be.  If e has already been
       calculated to a sufficient number of digits, this returns it from the
       cache.
   */
   class getE[digits = getPrecision[]] :=
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[digits]
         if (largestDigits >= digits)
            return 1. * cacheE
         else
            return 1. * calcE[digits]
      }
      finally
         setPrecision[origPrec]
   }


   /** This is an internal method that calculates the digits of e if
       necessary. */
   class calcE[digits] :=
   {
      origPrec = getPrecision[]
      try
      {
         setPrecision[15]
         // Find number of terms to calculate.
         // ln[x!] = ln[1] + ln[2] + ... + ln[x]
         logMax = digits * ln[10]
         k = largestK
         
         while (logFactorial < logMax)
         {
            logFactorial = logFactorial + ln[k]
            k = k + 1
         }

         setPrecision[digits+3]
         if largestK < k    // We may have already calculated enough!
         {
            if (largestK == 1)  // This is the first iteration
            {
               p = P[0,k]
               q = Q[0,k]
            } else
            {
               // Continuing iterations
               pmb = P[largestK, k]

               qmb = Q[largestK, k]

               p = largestP qmb + pmb
               q = largestQ qmb
            }

            // Store the biggest values back in the cache
            largestP = p
            largestQ = q
            largestK = k
            
            cacheE = 1 + (1. * p)/q
         }
         
         setPrecision[digits]

         return 1. * cacheE
      }
      finally
         setPrecision[origPrec]
   }

   /** Private method for binary splitting. */
   class P[a,b] :=
   {
      if (b-a) == 1
         return 1
      
      m = (a+b) div 2
      r = P[a,m] Q[m,b] + P[m,b]
      return r
   }

   /** Private method for binary splitting. */
   class Q[a,b] :=
   {
      if (b-a) == 1
         return b
      
      m = (a+b) div 2
      return Q[a,m] Q[m,b]
   }
}