Taylor series and Horner’s rule

How does a computer work out cos 1? It has no circle to measure. What it can do is add and multiply, and cos x can be written, as closely as we like, using only those:

That is the start of its Taylor series; the sum goes on for ever, and we keep the first six terms. The question for this page is how to add them up with as little work as possible. The instrument does it three ways and counts the multiplications.

⋯9static const double coef[6] = { 1.0, -1.0 / 2, 1.0 / 24, -1.0 / 720, 1.0 / 40320, -1.0 / 3628800 };10 11double horner_rec(int k, int n, double u)12{13    if (k == n - 1)14        return coef[k];15    double rest = horner_rec(k + 1, n, u);16    mults++;17    return coef[k] + u * rest;18}19 20double horner_loop(int n, double u)21{22    double r = coef[n - 1];23    for (int k = n - 2; k >= 0; k--) {24        r = r * u + coef[k];25        mults++;26    }27    return r;28}⋯48int main(void)49{50    double x = 1.0;51    double u = x * x;52    double r = horner_rec(0, TERMS, u);53    printf("recursive: %.10f in %d multiplications\n", r, mults);54    mults = 0;55    double l = horner_loop(TERMS, u);56    printf("loop:      %.10f in %d multiplications\n", l, mults);57    mults = 0;58    double s = cos_naive(x, TERMS);59    printf("naive:     %.10f in %d multiplications\n", s, mults);60    return 0;61}

A polynomial, nested

Every term has an even power of x. So write u for x², and the six terms become a polynomial in u:

So a₀ = 1, a₁ = −1/2, a₂ = 1/24, a₃ = −1/720, a₄ = 1/40320, a₅ = −1/3628800. These are the six numbers in the program’s array coef.

Now take u out of every term after the first, then again, and again:

This is Horner’s rule. To evaluate it, start in the innermost bracket and work outwards. Each bracket is its own coefficient plus u times the bracket inside it: one multiplication and one addition. Watch the column in the clip fill from the bottom: the innermost bracket is just a₅, and every row above uses the value of the row below. At x = 1, so u = 1, the top row ends at 0.5403023038.

The code and its cost

By recursion

Each bracket is “my coefficient plus u times the rest”, and the rest is the same kind of thing, one bracket further in. That is a recursive function:

double horner_rec(int k, int n, double u)
{
    if (k == n - 1)
        return coef[k];
    double rest = horner_rec(k + 1, n, u);
    mults++;
    return coef[k] + u * rest;
}

horner_rec(k, n, u) evaluates the brackets from $a_k$ inwards. The base case is the innermost bracket, $a_{n-1}$ alone. main calls horner_rec(0, 6, u). As in Recursion (4.1), the calls go down first, and the work happens on the way back, which is exactly the inside-out order Horner’s rule needs.

By a loop

The work on the way back is the same at every step, so a loop can do it from the inside out directly:

double horner_loop(int n, double u)
{
    double r = coef[n - 1];
    for (int k = n - 2; k >= 0; k--) {
        r = r * u + coef[k];
        mults++;
    }
    return r;
}

r starts at the innermost coefficient and passes through the same values as the column: −0.0000002756, 0.0000245260, and so on up to 0.5403023038. It uses one frame instead of six.

Counting the multiplications

With n terms, both Horner versions make n − 1 multiplications: one for each bracket but the innermost. Six terms: 5. (Squaring x to get u is one more, done once in main; the program does not count it.)

The naive way, cos_naive in the program, builds each term from scratch: for term k it multiplies x by itself 2k times for the power and multiplies 1 × 2 × … × 2k for the factorial, then divides. That is 4k multiplications for term k, and over n terms 4(0 + 1 + … + (n − 1)) = 2n(n − 1). Six terms: 60. Drag the slider and read the counts: Horner’s grow by one per term, the naive count grows with the square of n.

How good is it?

The six terms give 0.5403023038; the true cos 1 is 0.5403023059. The difference is about 2.1 × 10⁻⁹. The first term left out is 1/12! ≈ 2.09 × 10⁻⁹, nearly the same. That is no accident: once the terms of a series like this one, with alternating signs, are getting smaller, stopping early leaves an error no larger than the first term left out.

More terms do not make it good everywhere at once. Look at the curve with 6 terms: near 0 it lies on top of cos x, but towards ±2π it peels away, because x¹⁰ grows huge there. Drag down to 3 or 4 terms and watch the good region shrink towards 0. A maths library first moves x into a small range near 0, using cos’s repeats, and only then uses a few terms.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic cosine.c -o cosine
./cosine
gcc -std=c11 -Wall -Wextra -pedantic -DTERMS=3 cosine.c -o cosine3
./cosine3

The first prints the three lines with 0.5403023038, in 5, 5 and 60 multiplications.

Worked example

Six terms at x = 1, by hand. With u = 1, each step is just “add the next coefficient”. Start inside:

Before you read on: what does the program print with two terms? (0.5: the polynomial is 1 − u/2, and u = 1.)

Three terms at x = 0.5. Now u = 0.25, so 1 − 0.25/2 + 0.25²/24 = 1 − 0.125 + 0.0026041667 = 0.8776041667. The true cos 0.5 is 0.8775825619, an error of 2.2 × 10⁻⁵. The first term left out is 0.5⁶/720 ≈ 2.17 × 10⁻⁵: again the error is just under it.

A polynomial of degree 20. Horner’s rule takes 20 multiplications. Computing each power of u separately, u¹ needs none, u² one, and so on up to u²⁰, nineteen: 0 + 1 + … + 19 = 190, and then 20 more to multiply each power by its coefficient, 210 in all.

Where it’s used

Maths libraries. The cos, sin and exp in a C library compute a short polynomial for x near 0, evaluated by Horner’s rule, after moving x into that small range. The polynomial’s coefficients are chosen a little differently from the Taylor ones, to spread the error evenly, but the evaluation is the same nesting.

Reading numbers and hashing text. Turning the digits “4”, “0”, “3” into the number 403 is Horner’s rule with u = 10: ((4 × 10) + 0) × 10 + 3. String hashing (18.1) and CRC checksums use the same nesting over the characters of a text.

One DFT value. The Goertzel algorithm computes a single frequency of a signal by the same nesting, with the signal’s samples as the coefficients; see also Sinusoids.

Reference card

Platform: C11, gcc on x86-64 Linux, double 8 bytes; the program prints 10 decimal places.