Skip to the content of the web site.

5.4 Horner's Rule

Introduction Theory HOWTO Error Analysis Examples Questions Applications in Engineering Matlab Maple BLAS

Introduction

Horner's rule provides a fast and efficient method of evaluating polynomials, including those in expanded form and those in the form of Newton polynomials. We will focus on Horner's rule as applied to Newton polynomials.

Background

Useful background for this topic includes:

References

Theory

The polynomial p(x) = c1 xn − 1 + ⋅⋅⋅ + cn found by the Vandermonde method to interpolate n points may be evaluated at a point x in a reasonably straightforward manner:

Set s = c1,
For i = 2, 3, ..., n do
    s = s⋅x + ci.

This requires n − 1 multiplications and additions. Unfortunately, the numerical stability of such an expanded polynomial is difficult to determine, for example, subtractive cancellation could occur at any time and it is difficult to determine apriori at which points this problem will occur.

The Newton polynomial polynomial does not appear to be as easy to evaluate:

p(x) = c1 + c2(x − x1) + ⋅⋅⋅ + cn(x − x1)⋅⋅⋅(x − xn − 1)

This requires n − 1 multiplications but 2n − 2 additions. The benefit, however, is more significant. The Newton polynomial may be found using O(n) memory and the numerical stability is much better defined.

Both of the two methods above are termed Horner's rule, but we will focus on the application of Horner's rule using Newton polynomials.

HOWTO

Given the n points (x1, y1), ...,(xn, yn), and the coefficients c = (c1, ..., cn)T which define the Newton polynomial:

p(x) = c1 + c2(x − x1) + ⋅⋅⋅ + cn(x − x1)⋅⋅⋅(x − xn − 1)

we may efficiently evaluate this polynomial at a given point x as follows:

s = cn,
for i = n − 1, n − 2, ..., 3, 2, 1 do
    s = s⋅(x − xi) + ci.

Error Analysis

Suppose we are evaluating a Newton polynomial at a known point x. In this case, it is desirable to order the points being interpolated relative to their distance from x.

For example, suppose we are interpolating the points (2, 8.5), (3, 8.6), (4, 8.9) and we want to evaluate the Newton interpolating polynomial at x = 3.1 . In this case, we would order the points:

  3   8.6
              0.3
  4   8.9             0.1
              0.2
  2   8.5

and hence the interpolating Newton polynomial is 8.6 + 0.3(x − 3) + 0.1(x − 3)(x − 4).

Reference: J. Stoer and R. Bulirsch, p.47.

Examples

1. Use Horner's rule to evaluate the Newton polynomial defined by the points x = (0.5, 5.9, 1.3, 4.7, 3.5)T with corresponding coefficients c = (0.39, 0.47, 0.63, −0.53, 1.23)T at the points x = 3.7 and x = 4.2.

Starting with x = 3.7, we calculate:

1.23
1.23⋅(3.7 − 4.7) − 0.53 = −1.76
−1.76⋅(3.7 − 1.3) + 0.63 = −3.594
−3.594⋅(3.7 − 5.9) + 0.47 = 8.3768
8.3768⋅(3.7 − 0.5) + 0.39 = 27.196

Next, continuing with x = 4.2, we calculate:

1.23
1.23⋅(4.2 − 4.7) − 0.53 = −1.145
−1.145⋅(4.2 − 1.3) + 0.63 = −2.6905
−2.6905⋅(4.2 − 5.9) + 0.47 = 5.0439
5.0439⋅(4.2 − 0.5) + 0.39 = 19.052

In both cases, we used 8 additions and 4 multiplications.

2. Use Horner's rule to evaluate the Newton polynomial defined by the points x = (1, 2, 3, 4)T with corresponding coefficients c = (4, −3, 2, −1)T at the points x = 2.5 and x = 3.5.

−1
−1⋅(2.5 − 3) + 2 = 2.5
2.5⋅(2.5 − 2) − 3 = −1.75
−1.75*(2.5 − 1) + 4 = 1.375

−1
−1⋅(3.5 − 3) + 2 = 1.5
1.5⋅(3.5 − 2) − 3 = −0.75
−0.75*(3.5 − 1) + 4 = 2.125

3. Use Horner's rule to evaluate the Newton polynomial define by the points x = (1, 2, 4, 5)T with corresponding coefficients c = (11, 12, 13, 14)T at the point x = 3.

14
14⋅(3 − 4) + 13 = −1
−1⋅(3 − 2) + 12 = 11
11⋅(3 − 1) + 11 = 33

4. Compare the number of operations required if you were to evaluate the long form of the Newton polynomial given in Example 3 and compare to the number of operations used with Horner's rule.

The calculation 11 + 12 (x − 1) + 13 (x − 1)(x − 2) + 14 (x − 1)(x − 2)(x − 3) requires 9 additions and 6 multiplications. Horner's rule, from Example 3, requires 6 additions and 3 multiplications. More telling is the observation that the number of operations not using Horner's rule is O(n2) whereas the number of operations using Horner's rule is O(n), so simply trying to evaluate Newton polynomials is worse than the most efficient means of evaluating an expanded polynomial.

Questions

1. Use Horner's rule to evaluate the polynomial defined by the points x = (0.5, 0.7, 1.5, 1.8)T and the corresponding Newton polynomial coefficients c = (0.3, 0.8, -0.2, 0.6)T at the points x = 1.3, 1.5, and 2.5.

Answer: 0.78640, 0.94000, 3.3400

2. The expanded Newton polynomial in Question 1 is 0.6 x3 - 1.82 x2 + 2.33 x - 0.485. Write an algorithm in Matlab which evaluates this polynomial. Test your algorithm on the points given in Question 1.

3. Use Horner's rule to evaluate the polynomial defined by the points x = (8, 2, 6, 4)T and the corresponding Newton polynomial coefficients c = (-2, 2, 1, -1)T at the points x = 3, 5, 7.

Answer: -32, -26, -4

4. Use Horner's rule to evaluate the polynomial defined by the points x = (1.5, 2, 2.5)T and the corresponding Newton polynomial coefficients c = (-2, 0, 2)T at the points x = 1.75, 2.25.

Answer: -2.125, -1.625

Applications to Engineering

Horner's rule is the most efficient method of evaluating a dense polynomial at a particular value, both in terms of number of operations and even in terms of the number of registers. Thus, in any application where such evaluations are required, it is fast and efficient, and usually overlooked.

Matlab

The implementation of Horner's rule is quite straight-forward in Matlab:

>> xpts =      [3.5 4.7 5.1  6.9 8.3]'
>> coefflist = [5.3 3.2 2.1 -1.8 1.7]'
>> x = 5.5235423
>> s = xpts(end)
>> for i = (length(xpts) - 1):-1:1
    s = s*(x - xpts(i)) + coefflist(i)
end

Maple

This is most quickly done in Maple by doing:

xpts := [3.5, 4.7, 5.1,  6.9, 8.3];
coefflist := [5.3, 3.2, 2.1, -1.8, 1.7];
s := xpts[-1];
x := 5.5235423;
for i from -2 to -nops( xpts ) by -1 do
    s := s*(x - xpts[i]) + coefflist[i];
end do:
s;

BLAS

The basic linear algebra subroutine (BLAS) for for evaluating a Newton polynomial at a point x using Horner's rule is

  int gsl_poly_dd_eval( double dd[], double xs[], int N, double x );

This example finds the coefficients of the Newton polynomial interpolating the stated points

#include <stdio.h>
#include <gsl/gsl_poly.h>

/*************************************************
 * Evaluate a Newton polynomial at a number of
 * points using Horner's rule
 *************************************************/

#define N 5

int main() {
        /* The points
             (0.1, 1.2), (0.5, 2.7), (0.7, 3.4), (1.2, 4.7), (1.5, 6.0) */

        double xs[N] = {0.1, 0.5, 0.7, 1.2, 1.5};
        double ys[N] = {1.2, 2.7, 3.8, 4.7, 6.0};

        /* The coefficients of the newton polynomial. */
        double dd[N];
        gsl_poly_dd_init( dd, xs, ys, N );

        double x;

        for ( x = -2.0; x <= 2.0; x += 0.5 ) {
                printf( "(%+2.1lf, %12.8lf)\n",
                        x, gsl_poly_dd_eval( dd, xs, N, x ) );
        }

        return 0;
}

To compile and execute this code, you must have the GNU Scientific Library installed on your system, either in cygwin or in Linux. The example below looks at how you compile this code in Cygwin:

$ gcc -Wall -I/usr/include horner.c -lgsl -lgslcblas -lm
$ ./a.exe

The output is

(-2.0, 629.79090909)
(-1.5, 273.02857143)
(-1.0,  92.92857143)
(-0.5,  19.97792208)
(+0.0,   1.70909091)
(+0.5,   2.70000000)
(+1.0,   4.57402597)
(+1.5,   6.00000000)
(+2.0,  22.69220779)

As you may note, when the Newton polynomial is evaluated at the points x = 0.5 and x = 1.5, the results are y = 2.7 and y = 6.0, respectively, as expected.

Reference: GNU Scientific Library, 6.2 Divided Difference Representation of Polynomials.