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