Arbitrary-precision elementary mathematical functions (Python)

This program is under development.
is complete, remove the {{develop}} tag.

The aim of this program is to implement elementary mathematical functions, in particular those all functions in Python's standard math module, for use with the Decimal class. The functions should ideally handle arbitrary precision and always return a correctly rounded result.

General notes

To avoid rounding errors, the functions must temporarily increase the working precision for intermediate calculations. This is done by incrementing getcontext().prec. After a calculation has finished, the precision is decremented and the value to be returned is normalized to this precision with +x.

Exponential function

The exponential function is most easily calculated from its Taylor series,

$e^x = 1 + \frac{x}{1!} + \frac{x^2}{2!} + \frac{x^3}{3!} + \dots$
<<exp>>=
from decimal import *

def exp(x):
getcontext().prec += 3
if x < 0:
s = 1 / exp(abs(x))
else:
xpower = 1
ns = 0
s = 1
n = 0
factorial = 1
while s != ns:
s = ns
term = Decimal(xpower) / factorial
ns = s + term
xpower *= x
n += 1
factorial *= n
getcontext().prec -= 3
return +s

def test_exp():
assert exp(4) == +Decimal("54.59815003314423907811026120286")
assert exp(0) == 1
assert exp(-8) == +Decimal("0.0003354626279025118388213891257809")
assert exp(Decimal("0.6931471805599453094172321214582")) == 2


Natural logarithm

With the exponential function available, the natural logarithm can be computed easily using Newton's method.

<<log main>>=
from math import log as _flog

def log(x):
if x < 0:
return Decimal("NaN")
if x == 0:
return Decimal("-inf")
getcontext().prec += 3
eps = Decimal("10")**(-getcontext().prec+2)
# A good initial estimate is needed
r = Decimal(repr(_flog(float(x))))
while 1:
r2 = r - 1 + x/exp(r)
if abs(r2-r) < eps:
break
else:
r = r2
getcontext().prec -= 3
return +r


To test the function, we can calculate the logarithms of a few values and feed them to exp to hopefully obtain the original values. Since two separate operations are performed, a tiny epsilon must be allowed to account for rounding error.

<<log test>>=
def test_log():
A = [4, 1, Decimal("0.67753892")]
for p in [5, 40, 28]:
getcontext().prec = p
eps = Decimal("10")**(-getcontext().prec + 2)
for a in A:
w = exp(log(a))
assert abs(w - a) < eps

<<log>>=
log main
log test


Wrapping it up

<<hpfunc.py>>=
exp
log

if __name__ == "__main__":
test_exp()
test_log()


hijacker
hijacker
hijacker
hijacker