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*) = *c*_{1} *x*^{n − 1} + ⋅⋅⋅ + *c*_{n}
found by the Vandermonde method to interpolate *n* points
may be evaluated at a point *x* in a reasonably straightforward manner:

Set *s* = *c*_{1},

For *i* = 2, 3, ..., *n* do

*s* = *s*⋅*x* + *c*_{i}.

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*) = *c*_{1} + *c*_{2}(*x* − *x*_{1}) + ⋅⋅⋅ +
*c*_{n}(*x* − *x*_{1})⋅⋅⋅(*x* − *x*_{n − 1})
This requires *n* − 1 multiplications but 2*n* − 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 (*x*_{1}, *y*_{1}), ...,(*x*_{n}, *y*_{n}), and
the coefficients **c** = (*c*_{1}, ..., *c*_{n})^{T}
which define the Newton polynomial:

p(*x*) = *c*_{1} + *c*_{2}(*x* − *x*_{1}) + ⋅⋅⋅ +
*c*_{n}(*x* − *x*_{1})⋅⋅⋅(*x* − *x*_{n − 1})
we may efficiently evaluate this polynomial at a given point *x* as follows:

*s* = *c*_{n},

for *i* = *n* − 1, *n* − 2, ..., 3, 2, 1 do

*s* = *s*⋅(*x* − *x*_{i}) + *c*_{i}.

# 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(*n*^{2}) 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 x^{3} - 1.82 x^{2} + 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.