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 20152 days, 17 hours, 18 minutes ago.