/** 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[n, 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

    For better code samples, see the new (2023) paper:
    http://www.hvks.com/Numerical/papers.html
    and especially the paper
    http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Exp()%20calculation%20for%20arbitrary%20precision.pdf

    which contains implementations of the Taylor series, binary splitting,
    and sinh algorithms to quickly calculate exp(x).  Implementing it in terms
    of sinh is apparently fastest and would give you a sinh algorithm...
*/

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 = abs[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
         if reps < 1
            reps = 1
         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)]
}