Download or view BinarySplitting.frink in plain text format
/** 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
}
Download or view BinarySplitting.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 19965 days, 22 hours, 35 minutes ago.