BinarySplittingExp.frink

Download or view BinarySplittingExp.frink in plain text format


/** This class is an example of using the BinarySplitting.frink file to
    calculate exp[x] to arbitrary precision.  It follows the algorithm
    described in:

   "Fast multiprecision evaluation of series of rational numbers" by Bruno
    Hible and Thomas Papanikolaou:

    https://ginac.de/CLN/binsplit.pdf
  
    See the references in that paper for discussion of the parameters.  See
    section 2.2.1 for the algorithm for exp[x].

    Call BinarySplittingExp.exp[digits] to calculate the specified number of
    digits.

    TODO:  This is incomplete and needs a better algorithm to calculate the
    number of iterations because the equation given in the paper is not useful.

    Also see:
    https://stackoverflow.com/questions/57510825/binary-splitting-in-pari-gp
*/


use BinarySplitting.frink
    
class BinarySplittingExp implements ParamProvider
{
   // The numerator
   var u

   // The denominator
   var v

   new[num, denom] :=
   {
      u = num
      v = denom
   }
   
   /** Call this method with the number of digits of exp[x] you want to
   calculate.  */

   class exp[x, digits] :=
   {
      [u,v] = numeratorDenominator[toRational[x]]
      splitter = new BinarySplittingExp[u,v]

      origPrec = getPrecision[]
      try
      {
         setPrecision[18]
         d1 = log[digits]
         d2 = log[1/abs[x]]
         d = d1 + d2
         println["d1 is $d1"]
         println["d2 is $d2"]
         println["d  is $d"]
         reps = ceil[digits / d * 1.43]  // Found empirically for some numbers
         setPrecision[digits+3]
         println["reps is $reps"]
         return BinarySplitting.calc[splitter, 0, reps, digits]
      }
      finally
         setPrecision[origPrec]
   }

   /** Implementation of the ParamProvider interface. */
   a[n] := 1

   b[n] := 1

   p[n] :=
   {
      if n == 0
         return 1
      else
         return u
   }

   q[n] :=
   {
      if n <= 0
         return 1
      else
         return n v
   }
}


use ArbitraryPrecision.frink

test[x,digits] :=
{
   s = now[]
   c = arbitraryExp[x, digits]
   e = now[]
   println["calculated c  in " + format[e-s, "ms", 0]]

   s = now[]
   bs = BinarySplittingExp.exp[x, digits]
   e = now[]
   println["calculated bs in " + format[e-s, "ms", 0]]
   
   setPrecision[digits]
   println["diff is " + (c-bs)]
}


Download or view BinarySplittingExp.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 19760 days, 22 hours, 41 minutes ago.