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[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)]
}
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 20139 days, 6 hours, 54 minutes ago.