Arbitrary-precision elementary mathematical functions (Python)

From LiteratePrograms

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.

Contents 1 General notes 2 Exponential function 3 Natural logarithm 4 Trigonometric functions 5 Hyperbolic functions 6 Constants 7 Wrapping it up

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,

<<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

Trigonometric functions

Hyperbolic functions

Constants

Wrapping it up

<<hpfunc.py>>=
exp
log
if __name__ == "__main__":
    test_exp()
    test_log()