Download or view ArbitraryPrecision.frink in plain text format
/** Functions for performing calculations to arbitrary precision.
A good reference is Ronald W. Potter,
"Arbitrary Precision Calculation of Selected Higher
Functions." References to "Potter" and equation numbers are from this book.
Also see Henrik Vestermark's excellent site which is a treasure trove:
http://www.hvks.com/
especially
http://www.hvks.com/Numerical/papers.html
*/
use pi2.frink
// See http://en.literateprograms.org/Arbitrary-precision_elementary_mathematical_functions_(Python)
/** This is a somewhat naive implementation of e^x which has a Maclaurin
series of
exp[x] = 1 + x + x^2/2! + x^3/3! + ...
It could probably be done faster with a binary splitting algorithm.
See BinarySplittingExp.frink for a sample.
*/
arbitraryExp[x, digits=getPrecision[], debug=false] :=
{
if debug
println["in arbitraryExp[$x, $digits]"]
origPrec = getPrecision[]
setPrecision[digits+3]
var s
try
{
if x < 0
{
s = 1.0 / arbitraryExp[abs[x],digits,debug]
} else
{
xpower = 1.
ns = 0.
s = 1
n = 0.
factorial = 1.
while s != ns
{
s = ns
term = xpower / factorial
if debug
println["term is $term"]
ns = s + term
xpower = xpower * x
n = n + 1.
factorial = factorial * n
}
if debug
println["s is $s"]
}
setPrecision[digits]
retval = 1.0 s
if debug
println["arbitraryExp[$x, $digits] returning $retval"]
return retval
}
finally
setPrecision[origPrec]
}
/** Natural log to arbitrary precision. This uses a cubic convergence
algorithm (that is, the number of correct digits in the result
approximately triple with each iteration) with adaptive precision using
equation 3.47 in Potter. It is significantly faster than the previous
algorithm that did not have adaptive precision and used a quadratic
Newton's method algorithm.
*/
arbitraryLn[x, digits=getPrecision[], debug=false] :=
{
if debug
println["in arbitraryLn2[$x]"]
if x <= 0
return "Error: Arbitrary logs of negative numbers not yet implemented."
origPrec = getPrecision[]
try
{
setPrecision[10]
eps = 10.^-(digits+1)
// A good initial estimate is needed
if (x < 10^308)
{
y = ln[x]
prec = 15
} else // TODO: Store ln[2] somewhere.
{
y = approxLog2[x] * ln[2]
prec = 10 // Is this reasonable? approxLog2 has a lot of latitude.
}
if debug
println["Epsilon is $eps"]
// Use Newton's method
do
{
setPrecision[prec]
y2 = y
le = arbitraryExp[-y, prec]
zn = 1 - x le
y = y - zn(1 + zn/2)
if debug
println["y is $y"]
prec = prec * 3
if (prec > digits + 3)
prec = digits + 3
} while (prec < digits) or (abs[y2-y] > eps)
setPrecision[digits]
retval = 1. y
if debug
println["arbitraryLn about to return $retval"]
return retval
}
finally
setPrecision[origPrec]
}
// Arbitrary-precision power x^p
// This uses the relationship that x^p = exp[p * ln[x]]
arbitraryPow[x, p, digits = getPrecision[], debug=false ] :=
{
// TODO: Make this work much faster for integer and rational powers.
prec = getPrecision[]
try
{
workingdigits = digits + 2
if digits <= 12
workingdigits = digits + 4
setPrecision[workingdigits]
ret = arbitraryExp[p * arbitraryLn[x, workingdigits, debug],
workingdigits, debug]
setPrecision[digits]
return 1. * ret
}
finally
setPrecision[prec]
}
// Arbitrary log to the base 10.
arbitraryLog[x, digits=getPrecision[], debug=false] :=
{
origPrec = getPrecision[]
try
{
setPrecision[digits+2]
// TODO: Store ln[10] somewhere.
retval = arbitraryLn[x, digits+2] / arbitraryLn[10, digits+2, debug]
setPrecision[digits]
return 1. retval
}
finally
setPrecision[origPrec]
}
// This method computes sine of a number to an arbitrary precision.
// This method is actually a dispatcher function which conditions the values
// and tries to dispatch to the appropriate method which will be most likely
// to converge rapidly.
arbitrarySin[x, digits=getPrecision[], debug=false] :=
{
origPrec = getPrecision[]
try
{
// If x is something like 10^50, we actually need to work with
// 50+digits at this point to get a meaningful result.
extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10
if debug
println["Extradigits is " + extradigits]
setPrecision[digits+extradigits+4]
if debug
println["Dividing by pi to " + (digits + extradigits + 4) + " digits"]
pi = Pi.getPi[digits+extradigits+4]
// Break up one period of a sinewave into octants, each with width pi/4
octant = floor[(x / (pi/4)) mod 8]
// Adjust x into [0, 2 pi]
x = x mod (2 pi)
if debug
println["Octant is $octant"]
if debug
println["Adjusted value is $x"]
if octant == 0
val = arbitrarySinTaylor[x, digits]
else
if octant == 1
val = arbitraryCosTaylor[-(x - pi/2), digits]
else
if octant == 2
val = arbitraryCosTaylor[x - pi/2, digits]
else
if octant == 3 or octant == 4
val = -arbitrarySinTaylor[x-pi, digits]
else
if octant == 5
val = -arbitraryCosTaylor[-(x - 3/2 pi), digits]
else
if octant == 6
val = -arbitraryCosTaylor[x - 3/2 pi, digits]
else
val = arbitrarySinTaylor[x - 2 pi, digits]
setPrecision[digits]
return 1. * val
}
finally
setPrecision[origPrec]
}
/* This method computes cosine of a number to an arbitrary precision.
This method actually just calls arbitrarySin[x + pi/2]
*/
arbitraryCos[x, digits=getPrecision[]] :=
{
origPrec = getPrecision[]
// If x is something like 10^50, we actually need to work with
// 50+digits at this point to get a meaningful result.
extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10
if debug
println["Extradigits is " + extradigits]
if debug
println["Dividing by pi to " + (digits + extradigits + 4) + " digits"]
pi = Pi.getPi[digits+extradigits+4]
try
{
setPrecision[digits+extradigits+4]
pi = Pi.getPi[digits+extradigits+4]
arg = x+pi/2
return arbitrarySin[arg, digits]
}
finally
setPrecision[origPrec]
}
// Arbitrary-precision sine
arbitrarySinTaylor[x, digits=getPrecision[], returnInterval = false] :=
{
origPrec = getPrecision[]
eps = 10.^-(digits+3)
terms = new array
try
{
setPrecision[digits+5]
x = x / radian // Factor out radians if we use them
pi = Pi.getPi[digits+5]
x = x mod (2 pi)
if x > pi
x = x - 2 pi
num = x
sum = x
term = 3
denom = 1
factor = -x*x
terms.push[sum]
do
{
prevSum = sum
num = num * factor
denom = denom * (term-1) * term
part = num/denom
sum = sum + part
term = term + 2
terms.push[part]
} while prevSum != sum
// println["terms for sin is $term"]
sum = sum[reverse[terms]]
setPrecision[digits]
return 1. * sum
}
finally
setPrecision[origPrec]
}
// Cosine for arbitrary digits. We could write this in terms of the sine
// function (cos[x] = sin[x + pi/2]) but it's faster and more accurate
// (especially around pi/2) to write it as the Taylor series expansion.
arbitraryCosTaylor[x, digits=getPrecision[]] :=
{
origPrec = getPrecision[]
eps = 10.^-(digits+4)
terms = new array
try
{
setPrecision[digits+4]
x = x / radian // Factor out radians if we use them
pi = Pi.getPi[digits+4]
x = x mod (2 pi)
// println["Effective x is $x"]
num = 1
sum = 1
term = 2
denom = 1
factor = -x*x
terms.push[sum]
do
{
prevSum = sum
num = num * factor
denom = denom * (term-1) * term
part = num/denom
sum = sum + part
// println["$term $part $sum"]
term = term + 2
terms.push[part]
} while prevSum != sum
// println["terms for cos is $term"]
sum = sum[reverse[terms]]
setPrecision[digits]
return 1. * sum
}
finally
setPrecision[origPrec]
}
// Tangent for arbitrary digits. This is written in terms of
// sin[x]/cos[x] but it seems to behave badly around pi/2 where cos[x] goes to
// zero.
//
// TODO: Make this a series expansion with the tangent numbers. This might
// be more efficient.
// See: http://mathworld.wolfram.com/TangentNumber.html
// also TangentNumbers.frink
// which calculate these numbers directly and efficiently.
//
// We could also try using Newton's method to invert arctan[x] which
// has a simple series expansion,
// arctan[x] = sum[(-1)^k x^(2k+1) / (2k + 1), {k, 0, infinity}]
// but this only converges for abs[x] <= 1, x != +/- i
//
// See:
// Fast Algorithms for High-Precision Computation of Elementary Functions,
// Richard P. Brent, 2006
// https://pdfs.semanticscholar.org/bf5a/ce09214f071251bfae3a09a91100e77d7ff6.pdf
arbitraryTan[x, digits=getPrecision[], debug=false] :=
{
// If x is something like 10^50, we actually need to work with
// 50+digits at this point to get a meaningful result.
extradigits = max[0, ceil[approxLog2[abs[x]]/ 3.219]] // Approx. log 10
if debug
println["Extradigits is " + extradigits]
origPrec = getPrecision[]
try
{
setPrecision[digits+extradigits+4]
retval = arbitrarySin[x, digits+4] / arbitraryCos[x, digits+4]
setPrecision[digits]
return 1. * retval
}
finally
setPrecision[origPrec]
}
// Polylogarithm. See:
// https://en.wikipedia.org/wiki/Polylogarithm
polylog[s, z, digits = getPrecision[]] :=
{
// if x <= 0
// return "Error: Arbitrary logs of negative numbers not yet implemented."
origPrec = getPrecision[]
try
{
setPrecision[digits+3]
eps = 10.^-(digits+1)
sum = 1. * z
term = 0
k = 2
do
{
term = z^k / k^s
sum = sum + term
k = k + 1
// println[sum]
} while abs[term] > eps
setPrecision[digits]
retval = 1. sum
return retval
}
finally
setPrecision[origPrec]
}
// Binary logarithm (that is, logarithm to base 2.)
binaryLog[x, digits = getPrecision[]] :=
{
origPrec = getPrecision[]
try
{
setPrecision[digits+3]
x = 1. x
y = 0
b = .5
while x < 1
{
x = 2 x
y = y - 1
}
while x >= 2
{
x = x / 2
y = y + 1
}
setPrecision[15]
epsilon = y * 10.^-(digits+3)
setPrecision[digits+3]
println["Epsilon is $epsilon"]
do
{
x = x^2
if x >= 2
{
x = x/2
y = y + b
}
b = b/2
//println["$y $x $b"]
} while b > epsilon
setPrecision[digits]
return 1. * y
}
finally
setPrecision[origPrec]
}
/** Performs a step of binarySplitting and returns [xp, p, q]
This is based on Henrik Vestermark's algorithm,
http://www.hvks.com/
especially
http://www.hvks.com/Numerical/papers.html
http://www.hvks.com/Numerical/arbitrary_precision.html
and specifically outlined in the paper:
http://www.hvks.com/Numerical/Downloads/HVE%20Fast%20Exp()%20calculation%20for%20arbitrary%20precision.pdf
*/
binarySplittingExp[x, xp, a, b] :=
{
// println["In binarySplittingExp a = $a b=$b"]
diff = b-a
if diff == 1
{
p = xp * x
return [p, p, b]
}
if diff == 2
{
xp = xp * x
p = (x+b) xp
xp = xp * x
q = b (b-1)
return [xp, p, q]
}
mid = (a+b) div 2
[xp, p, q] = binarySplittingExp[x, xp, a, mid] // Interval a...mid
[xp, pp, qq] = binarySplittingExp[x, xp, mid, b] // Interval a...mid
return [xp, p*qq + pp, q * qq]
}
arbitraryExpSplitting[x, digits=getPrecision[]] :=
{
origPrec = getPrecision[]
try
{
setPrecision[18]
prec = digits + 2 + ceil[log[digits]]
c1 = 1
c2 = 2
v = x
xp = 1
// Automatically calculate optimal reduction factor as a power of two
r = 8 * ceil[ln[2] * ln[prec]]
r = r + 1 // r += v.exponent() + 1
r = max[0, r]
// Adjust the precision
prec = prec + floor[log[prec]] * r
println["Final precision is $prec"]
setPrecision[prec]
// e^-x == 1/e^x
if realSignum[v] == -1
v = -v;
println["r is $r"]
println["v was $v"]
// v = v * 2^-r //v.adjustExponent(-r);
println["v is now $v"]
// Calculate needed Taylor terms
k = xstirlingApprox[prec]
if k < 2
k = 2 // Minimum 2 terms otherwise it can't split
[p1, q1] = numeratorDenominator[toRational[v]]
// q1 = q1 * 2^r
println["p1 = $p1, q1 = $q1"]
[xp, p, q] = binarySplittingExp[p1, q1, 0, k]
println["Out of binary splitting."]
println["xp=$xp"]
println["p=$p"]
println["q=$q"]
println["p/q is " + (p/q)]
// Adjust and calculate exp[x]
pp = q
p = p + pp
p = p / q
p = 1. p
println["Before adjust, p = $p"]
// Reverse argument reduction
// Brent enhancement avoids loss of significant digits when x is small.
/*
if r > 0
{
p = p - 1
while r > 0
{
println["r=$r"]
p = p *(p+2)
r = r - 1
}
p = p + 1
}
*/
println["Past reduction"]
if realSignum[x] == -1
p = 1/p
setPrecision[digits]
return 1. p
}
finally
{
setPrecision[origPrec]
}
}
// Stirling approximation for calculating e^x
xstirlingApprox[digits] :=
{
test = (digits + 1) * ln[10]
// x^n/k!<10^-p, where p is the precision of the number
// x^n~2^x’exponent
// Stirling approximation of k!~Sqrt(2*pi*k)(k/e)^k.
// Taking ln on both sides you get:
//
// -k*log(2^xexpo) + k*(log((k)-1)+0.5*log(2*pi*m)=test=>
// -k*xexpo*log(2) + k*(log((k)-1)+0.5*log(2*pi*m)=test
// Use the Newton method to find in less than 4-5 iteration
xold = 5
xnew = 0
NEWTONLOOP:
while true
{
f = -xold * ln[2] + xold * (ln[xold]-1) + 0.5 ln[2 pi xold]
f1 = 0.5 / xold + ln[xold]
xnew = xold - (f - test) / f1
if ceil[xnew] == ceil[xold]
break NEWTONLOOP
xold = xnew
}
println["xStirlingApprox returning " + ceil[xnew]]
return ceil[xnew]
}
/*
digits = 1
setPrecision[digits]
pi = Pi.getPi[digits]
num = pi/4
collapseIntervals[false]
inum = new interval[num, num, num]
println[arbitrarySin[num,digits]]
println[arbitrarySin[-num,digits]]
println[sin[num]]
println[sin[inum]]
println[]
println[arbitraryCos[num,digits]]
//println[arbitraryCosAround2Pi[num,digits]]
println[arbitraryCos[-num,digits]]
println[cos[num]]
println[cos[inum]]
println[]
println[arbitraryTan[num,digits]]
println[arbitraryTan[-num,digits]]
println[tan[num]]
println[tan[-num]]
println[tan[inum]]
setPrecision[2]
g = new graphics
ps = new polyline
pc = new polyline
for x=-20 pi to 20 pi step .1
{
ps.addPoint[x,-arbitrarySin[x]]
pc.addPoint[x,-arbitraryCos[x]]
}
g.add[ps]
g.color[0,0,1]
g.add[pc]
g.show[]
*/
Download or view ArbitraryPrecision.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 19762 days, 18 hours, 10 minutes ago.