How to use chebyshev polynomials to calculate exponents (antilogs)?

Question

I am trying to write an algorithm to accurately calculate exponents (antilogs) for a variable precision floating point library I am working on. The base is not relevant since I can convert between them.

I was able to manually calculate log10() using repetitive application of x^10. This is a digit by digit calculation and requires 4 multiplies per digit. I can reverse the algorithm to calculate exp10(), but this requires repeated application of a 10th root. Calculating the 10th root is significantly more CPU costly than 10th power.

I searched the web and a lot of people suggested using a Taylor Series to calculate exp_e(). I did that and found that it requires about 2 iterations per digit for accurate results. Only two multiplies and one divide per iteration. This is still a bit steep in terms of CPU cycles especially when some FP numbers can be 100 digits long.

Now, I also found the algorithm that was used to calculate EXP in the old Sinclair ZX81. The author claimed that it was Chebyshev polynomials. I mention this because when I tested it, the algorithm was calculating accurately to one digit per iteration - much better than the Taylor Series.

I would use the algorithm as-is if it weren't for the fact that the floating point library has to be accurate to an arbitrary number of digits. The ZX81 EXP code is only accurate to 8 digits. There is no explanation as to how to extend the number of iterations to get more accuracy.

So does anyone know how to calculate EXP() using Chebyshev Polynomials? Can they be expanded like the Taylor Series for more accuracy? Anything better than either?

[Please no long math proofs. That's over my head. I just want the algorithms.]

UPDATE The test results for the Taylor Series are as follows - LOG(25):

EXP10(1.397940) Round   1: = Taylor=19.16290731874155394000     
EXP10(1.397940) Round   2: = Taylor=23.36085084533393060000     
EXP10(1.397940) Round   3: = Taylor=24.64302976078316452000      
EXP10(1.397940) Round   4: = Taylor=24.93674192499081185000     
EXP10(1.397940) Round   5: = Taylor=24.99056707177124531000     
EXP10(1.397940) Round   6: = Taylor=24.99878698562735818000     
EXP10(1.397940) Round   7: = Taylor=24.99986296146780963000     
EXP10(1.397940) Round   8: = Taylor=24.99998619980410040000     
EXP10(1.397940) Round   9: = Taylor=24.99999874670913982000     
EXP10(1.397940) Round  10: = Taylor=24.99999989637041996000     
EXP10(1.397940) Round  11: = Taylor=24.99999999213623594000     
EXP10(1.397940) Round  12: = Taylor=24.99999999944868008000     
EXP10(1.397940) Round  13: = Taylor=24.99999999996408967000     
EXP10(1.397940) Round  14: = Taylor=24.99999999999782289000     
EXP10(1.397940) Round  15: = Taylor=24.99999999999988352000     
EXP10(1.397940) Round  16: = Taylor=25.00000000000000153000     
EXP10(1.397940) Round  17: = Taylor=25.00000000000000789000     
EXP10(1.397940) Round  18: = Taylor=25.00000000000000821000     
EXP10(1.397940) Round  19: = Taylor=25.00000000000000822000     

Taylor Series 16 digits of accuracy with long double and 19 iterations. Note that the optimal iterations in this example is 16 since those that follow are actually farther from the mark.

Sinclair EXP(1.397940)= 25.00000001907205060000 

Sinclair 9 digits of accuracy with double and 8 iterations

Here are the actual functions used: Taylor Series (iterates until no change):

// Taylor series to figure exp10^x

void TaylorEx(long double x)
{
    int i, j, intpart;
    long double a, frac;
    long double factorial;
    long double power, inp, old, out;

    // separate the int part from the frac part
    intpart = (int) x;
    frac = x - intpart;

    // The Taylor series operates on base E.
    // To convert base e to base 10, multiply input by ln(10)
    inp = frac * 2.3025850929940456840179914546844;

    factorial = 1;
    power = inp;
    a = 1;
    for (i = 1; i < 50; i++)
    {
        factorial *= i;

        old = a;
        a += (power / factorial);
        if (a == old) break;

        // for display, add base 10 exponent to A
        out = a;
        if (intpart > 0)
            for (j = 0; j < intpart; j++) out *= 10;
        else if (intpart < 0)
            for (j = 0; j < intpart; j++) out /= 10;

        printf("EXP10(%Lf) Round %3d: = Taylor=%.20Lf\tFPU=%.20Lf\n", x, i, out, powl(10,x));

        power *= inp;
    }
}

Sinclair (fixed at 8 iterations):

double Sinclair_Exp(double C)
{
    int N;
    double T, D, Z, BERG, M0, M1, M2, I, U;
    union {unsigned ui[2]; double f; } u;
    double A[8] =
    {
        0.000000001,        // A1   1 / 1000000000.0
        0.000000053,        // A2   1 / 18867924.528301886792452830188679
        0.000001851,        // A3   1 / 540248.51431658562938951917882226
        0.000053453,        // A4   1 / 18708.023871438459955474903185976
        0.001235714,        // A5   1 / 809.24874202283052550994809478569
        0.021446556,        // A6   1 / 46.627533110677537223225957584985
        0.248762434,        // A7   1 / 4.0198995640957589279738274308733
        1.456999875         // A8   1 / 0.68634185709864937359723520909705
    };

    // DEMONSTRATION FOR EXP X

    //D = C * 1.4426950408889634073599246810019;    // Log2(e) - uncomment for input base E
    D = C * 3.3219280948873623478703194294894;    // log2(10) - uncomment for inout base 10
    N = (int) D;
    Z = D - N;
    Z = 2 * Z - 1;

    // USE "SERIES CALCULATOR"
    // SERIES CALCULATOR
    // FIRST VALUE IN Z

    M0 = 2 * Z;
    M2 = 0;
    T = 0;
    for (I = 0; I < 8; I++)
    {
        M1 = M2;
        U = T * M0 - M2 + A[I];
        M2 = T;
        T = U;
    }
    T = T - M1;
    // LAST VALUE IN T

    // get original exponent of T
    u.f = T;
    u.ui[1] >>= 20;         // shift out mantissa
    u.ui[1] &= 0x7FF;       // mask off sign

    // Add correction
    N += u.ui[1];
    if (N > 2048) printf("Exponent Overflow!\n");

    if (N < 0.0) T = 0.0;
    else
    {
        // modify exponent
        u.f = T;
        u.ui[1] &= 0x800FFFFF;  // clear old exponent
        N <<= 20;               // shift new exponent into place
        u.ui[1] |= N;           // replace exponent
        T = u.f;
    }

    printf("Sinclair EXP(%lf)= %.20lf\tFPU EXP(%lf)=%.20lf\n", C, T, C, pow(10,C));


    return(T);
}

So am I stuck with the Taylor Series or is there a way to extend the Sinclair Chebyshev algorithm to n arbitrary precision?

Answer 1

Step 1 is range reduction. Exp(log(2)) = 2, and exp(k log 2 + y) = 2^k * exp(y). 2^k is easy to calculate, and you use that to reduce y to the range +/- log (2) / 2. So you get about -0.35 < y < 0.35, and even the Taylor formula works much better in this range.

Chebyshev polynomials are the optimal polynomials over a fixed range. For example draw the exp function in the reduced range and the Taylor polynomial of degree 1, which is 1+x. Just by moving the polynomial up or down a bit you reduce the maximum error. The best polynomial of degree n will interpolate the function at n+1 points, and drawing the error, it gives a maximum error +/- E with alternating sign in n+1 points.

PS. There seems to be some misconceptions that you can just “extend the range of a Taylor polynomial”. For most functions, you can’t. Try calculating exp(-100) and you’ll get an absolute nonsense result. Or try calculating sin x using a Taylor polynomial for x = 0, 0.1, 0.2 etc and tell me at what point the results have no relation to reality anymore.

Answer 2

Chebyshev approximations have fallen out of favor in recent decades. For fixed-precision libraries one typically uses the closely related minimax approximations, while for arbitrary-precision libraries, exp is indeed commonly computed via a Taylor series expansion. The original argument to the exponential is reduced to make the input to the Taylor series very small, which means fewer terms of the Taylor series have to be considered.

In a first step, j = floor (x * log2(e) and r = x - j * log(2). This allows us to compute exp(x) as 2^j * exp(r), with r in [0, log(2)]. Note that in the computation of r one needs to represent the constant log(2) with twice the precision of x to account for subtractive cancellation. The multiplication with 2^j is error free.

For further reduction, one can compute t = r / 2^k, thus exp (r) = (exp(t)) ^2k with t in [0, log(2)/2^k]. This requires repeated squarings, which can be faster than general multiplication by some constant factor. One needs about k additional bits of precision to counter the computational error.

It is also possible to accelerate this argument reduction via lookup tables of exp(i/2ⁿ), with i in [0, 2ⁿ-1], which is even done for normal double-precision computation by some standard math libraries.

At sufficiently high precision, we can compute u = sinh(t) via Taylor series expansion, then exp(t) = u + sqrt (1 + u²). This halves the required number of terms in the Taylor series at the expense of an arbitrary-precision square root. A square root is typically computed using division-free iterations for reciprocal square root with quadratic or cubic convergence.

As for the evaluation of the Taylor series itself, Smith showed how to re-arrange the computation to achieve a lower asymptotic complexity of O(t^1/3 M(t)) versus the complexity of O(t^1/2 M(t)) of a naive method, where M(t) is the time required to multiply two t-digit numbers:

David M. Smith, "Efficient Multiple-Precision Evaluation of Elementary Functions", Mathematics of Computation, Vol. 52, No. 185, January 1989, pp. 131-134

For optimal performance, you will probably need to develop heuristics for how to balance work in the argument reduction with work in evaluating the Taylor series, and when to switch to the Taylor series expansion for sinh.