2

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?

3
  • This seems to be more a Computer Science question than an Software Engineering one. You might get better answers over at Computer Science. – Bart van Ingen Schenau Jan 18 at 6:45
  • It’s numerical maths actually, considering that Chebyshev worked long before non-human computers were invented. – gnasher729 Jan 18 at 7:02
  • While heavy in maths, read this. Read the human parts and it'll help you understand what is going on and why. Highly recommend staring at the maths and teasing it apart to understand how it does it. – Kain0_0 Jan 18 at 22:23
2

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.

2
  • The Sinclair/Chebyshev exp() algorithm does split the input into int and fraction components. Only the fraction component is run through the series. Then the int portion is simply added to the exponent of the result. This works since the exp() is calculated in base 2. – user264480 Jan 18 at 13:21
  • When I said extend, I was referring to the number of fractional digits. With the Taylor series, this took about 18 iterations before the input and output stopped changing on a long double. I will try to post more details later. – user264480 Jan 18 at 13:32
1

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 2j * 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 2j is error free.

For further reduction, one can compute t = r / 2k, thus exp (r) = (exp(t)) 2k with t in [0, log(2)/2k]. 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/2n), with i in [0, 2n-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 + u2). 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(t1/3 M(t)) versus the complexity of O(t1/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.

10
  • I've seen "Chebychev polynomial" used in two different meanings: One, interpolation at nodes calculated according to some specific formula, but also exactly as the polynomial minimising the maximum error. – gnasher729 Feb 13 at 23:28
  • A truncated Chebyshev approximation basically constitutes a near-minimax approximation, but for the past 30 years or so math library implementers have preferred "true" minimax approximations generated numerically with the Remez algorithm (and possibly post processed with heuristics to adjust for quantization effects in finite-precision floating-point computation, e.g. Sollya's fpminimax command). – njuffa Feb 13 at 23:37
  • Minimax approximations follow the “Chebyshev theorem” (for any continuous function on an interval, a polynomial with n+2 maximum errors with opposite sign is optimal) and therefore are often called Chebyshev polynomials. – gnasher729 Feb 14 at 11:48
  • @gnasher729 If some people insist on using imprecise terminology, that is not something I can fix, but I don't intend to so the same. Relevant standard references do make the distinction I pointed out, e.g. Jean-Michel Muller, "Elementary Functions, Third Edition" (2016), p. 50: "As pointed out by Hart et. al. [225], if the function approximated is regular enough, then its Chebyshev approximation is very close to its minimax approximation (this is the case for the exponential function, see Figure 3.7) [...]". – njuffa Feb 14 at 12:08
  • I will research the Remez Algorithm. I have tried several approaches so far. I settled for using the sinclair chebyshev routine which only takes 11 multiplies. Then that result gets fed into the newton algorithm which usually resolves in a single iteration to the full precision. The bad thing is that each iteration has to call Log(x) which takes 4 multiplies per digit. The taylor series works well and only requires two multiplies per digit. Whatever method I use has to be extremely efficient computationally. – user264480 Feb 16 at 1:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.