/** This library implements a general algorithm for performing binary splitting
    calculations of arbitrary-precision calculations.  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

    To make this work for a particular binary splitting algorithm, you will
    need to implement the interface ParamProvider below in one of your own
    classes.  See BinarySplittingPi.frink for an example of calculating pi
    using binary splitting or BinarySplittingExp.frink for an (incomplete)
    implementation of the exp[x] function.
*/
class BinarySplitting
{
   /** This is usually the method you want to call.  It performs the binary
       splitting and the final "fixup" calculation S = T / (B Q) to obtain the
       final result.  If you need to manage significant digits better, you
       may need to call the split[] method yourself and then do that fixup
       calculation on the result.
   */
   class calc[arg is ParamProvider, n1, n2, digits = getPrecision[]] :=
   {
      result = split[arg, n1, n2]

      origPrec = getPrecision[]
      try
      {
         setPrecision[digits+3]
         // This is a final floating-point result
         S = result.T / (result.B * result.Q)

         setPrecision[digits]
         return 1. * S   // Round to the specified number of digits
      }
      finally
         setPrecision[origPrec]
   }

   /** This performs recursive binary splitting on particular arguments.
       You usually don't want to call this directly, but call the calc[]
       function.

       returns:  a SplitResult object containing the results of the calculation.
   */
   class split[arg is ParamProvider, n1, n2] :=
   {
      diff = n2 - n1
      if diff == 0
         return "split:  Error.  Difference n2 - n1 should be greater than 0."
      
      result = new SplitResult
      
      if diff == 1
      {
         // The result at the point n1
         p = arg.p[n1]
         result.P = p
         result.Q = arg.q[n1]
         result.B = arg.b[n1]
         result.T = arg.a[n1] * p
         return result
      }

      /* The paper above *very helpfully* says "cases 2, 3, 4 left out for
         simplicity".  Augh.  But these are the hard-to-derive and tricky parts.
         I've dug the implementations out of their source code for the CLN
         libraries (written in C) but haven't implemented them here yet.
         The libraries are available at:
         https://ginac.de/CLN/

         From their source, you probably want the very-obviously-named file
         src/float/transcendental/cl_LF_ratseries_pqab.cc
      
         and in that, maybe the function eval_pqab_series_aux
         to see how more of the recursive cases are defined.
      */

      // Default recursive case
      // find the middle of the interval
      nm = (n1 + n2) div 2

      // sum left side
      L = split[arg, n1, nm]

      // sum right side
      R = split[arg, nm, n2]

      // Put together
      result.P = L.P * R.P
      result.Q = L.Q * R.Q
      result.B = L.B * R.B
      result.T = R.B * R.Q * L.T + L.B * L.P * R.T
      return result
   }
}

/** This is an interface that you will implement in a class that provides
    the values of a[n], b[n], p[n], and q[n].  Since these are usually a
    simple function of n (and maybe some other parameters,) it's implemented
    as an interface with basic functions rather than by forcing the user to
    create a (possibly infinite) array of all the possible values of each
    variable for each n.  These functions should all return integers.
*/
interface ParamProvider
{
   // You must implement these functions to return values for integer values
   // of n, n>=0
   a[n]
   b[n]
   p[n]
   q[n]
}

/** This class is a simple data structure that holds the current values of
    P, Q, B, and T.  It may someday be used to *resume* (or somday checkpoint or
    parallelize) calculations.
*/
class SplitResult
{
   var P
   var Q
   var B
   var T
}