Arbitrary-precision elementary mathematical functions (Python)

From LiteratePrograms

Jump to: navigation, search

This program is under development.
Please help to debug it. When debugging
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,

from decimal import *
def exp(x):
    getcontext().prec += 3
    if x < 0:
        s = 1 / exp(abs(x))
        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:
            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 main
log test

Trigonometric functions

Hyperbolic functions


Wrapping it up

if __name__ == "__main__":