Horner's method

Horner's method evaluates a polynomial $p (x) = a_{n} x^{n} + \dots + a_{1} x + a_{0}$ by nesting the multiplications, $p (x) = (\dots ((a_{n} x + a_{n - 1}) x + a_{n - 2}) \dots) x + a_{0}$ , which takes $n$ multiplications and $n$ additions. Ostrowski and Pan proved that no method uses fewer for a general polynomial^[4]^[5]. The intermediate values of the computation are the coefficients of the quotient of $p (x)$ by $x - c$ , so the same loop performs synthetic division; repeating it gives the derivatives of $p$ at $c$ and the shifted polynomial $p (x + c)$ . The method is named after William Horner, who published a root-finding method built on it in 1819^[1], but it had been described by Ruffini in 1804^[2] and in thirteenth-century China^[10].

§ 01

The method

Rewriting the polynomial with nested brackets turns evaluation into a recurrence:

b_{n} = a_{n}, b_{k} = b_{k + 1} c + a_{k} (k = n - 1, \dots, 0), p (c) = b_{0}

Computing each power $x^{k}$ separately takes $k - 1$ multiplications for the power and one more for the coefficient, $n (n + 1) /2$ in all; keeping a running power takes $2 n$ ; Horner's recurrence takes $n$ . It also needs only one register for the running value, reads the coefficients in order, and is the usual way polynomials are evaluated in numerical libraries and in the kernels of elementary functions such as exp and sin, which are approximated by polynomials on small intervals.

Optimality

Ostrowski proved in 1954 that $n$ additions are necessary to evaluate a general polynomial of degree $n$ and that Horner's method is optimal in multiplications for degree at most 4^[4]. Pan proved in 1966 that $n$ multiplications are necessary for every degree, when the coefficients are arbitrary^[5]. If the same polynomial is evaluated many times, its coefficients can be preprocessed into a form that needs only about $n /2$ multiplications per evaluation^[6].

§ 02

Synthetic division

The numbers $b_{k}$ have a meaning of their own. Dividing $p (x)$ by $x - c$ gives

p (x) = (x - c) (b_{n} x^{n - 1} + b_{n - 1} x^{n - 2} + \dots + b_{1}) + b_{0}

which can be checked by comparing coefficients: the coefficient of $x^{k}$ on the right is $b_{k} - c b_{k + 1}$ , which equals $a_{k}$ by the recurrence. So the remainder of the division is $p (c)$ , the polynomial remainder theorem, and the quotient comes free. Dividing the quotient by $x - c$ again gives $p^{'} (c)$ as its remainder, and $k$ divisions give the Taylor coefficients $p^{(k)} (c) / k!$ . Carrying out all $n$ divisions computes $p (x + c)$ , the Taylor shift, in $O (n^{2})$ operations; it is the basic step of root isolation by Descartes' rule of signs.

Horner's method with a count of multiplications, synthetic division, value and derivative together, the Taylor shift, and Newton's method.

(* Coefficients are listed from the highest degree down:
   [| a_n; ...; a_1; a_0 |]. *)
let mults = ref 0
let ( *: ) a b = incr mults; a *. b
 
(* Each term computed from scratch: x^k takes k - 1 multiplications. *)
let naive a x =
  let n = Array.length a - 1 in
  let s = ref 0. in
  Array.iteri (fun i c -> let k = n - i in let pw = ref 1. in for _ = 1 to k do pw := !pw *: x done; s := !s +. (c *: !pw)) a;
  !s
 
(* Horner: a_n x^n + ... + a_0 = (...((a_n x + a_(n-1)) x + a_(n-2)) ...) x + a_0 *)
let horner a x =
  let acc = ref a.(0) in
  for i = 1 to Array.length a - 1 do acc := (!acc *: x) +. a.(i) done;
  !acc
 
(* The same recurrence yields the quotient on division by (x - c), with
   the remainder p(c) as the last value: synthetic division. *)
let synthetic a c =
  let b = Array.make (Array.length a) 0. in
  Array.iteri (fun i ai -> b.(i) <- (if i = 0 then ai else (b.(i - 1) *. c) +. ai)) a;
  (Array.sub b 0 (Array.length a - 1), b.(Array.length a - 1))
 
(* Value and derivative together, for Newton's method. *)
let horner2 a x = Array.fold_left (fun (p, d) c -> ((p *. x) +. c, (d *. x) +. p)) (0., 0.) a
 
(* p(x + c) by n rounds of synthetic division: the Taylor shift. *)
let taylor_shift a c =
  let a = Array.copy a and n = Array.length a - 1 in
  for k = 0 to n - 1 do for i = 1 to n - k do a.(i) <- a.(i) +. (a.(i - 1) *. c) done done;
  a
 
let show a = "[" ^ String.concat "; " (Array.to_list (Array.map (Printf.sprintf "%g") a)) ^ "]"

Running it.

degree 10 at x = 1.5: naive 490.985352 with 66 multiplications, Horner 490.985352 with 10
x^3 - 6x^2 + 11x - 6 divided by x - 4: quotient [1; -2; 3], remainder 6 = p(4)
divided by x - 2: quotient [1; -4; 3], remainder 0
Taylor shift p(x + 2): [1; 0; -1; 0]   (roots move to -1, 0, 1)
Newton step 1: x = 2.100000000000000
Newton step 2: x = 2.094568121104185
Newton step 3: x = 2.094551481698199
Newton step 4: x = 2.094551481542327
Newton step 5: x = 2.094551481542327

The last lines apply Newton's method to $x^{3} - 2 x - 5 = 0$ , the example Newton used to illustrate it^[11], with horner2 computing $p (x)$ and $p^{'} (x)$ in one pass. Horner's own paper used the division step to shift the polynomial so that the next digit of a root could be read off, a digit-by-digit method that was taught in schools well into the twentieth century^[1]^[3].

§ 03

Rounding errors

In floating-point arithmetic, Horner's method is backward stable: the computed value is the exact value of a polynomial whose coefficients differ from the $a_{k}$ by relative amounts of order $2 n u$ , where $u$ is the unit roundoff. The forward error is therefore bounded by^[7]:

p (x) - p (x); \leq; γ_{2 n} k = 0 \sum n ∣ a_{k} ∣ ∣ x ∣^{k}, γ_{2 n} = \frac{2 n u}{1 - 2 n u}

The relative error is small when $∣ p (x) ∣$ is comparable to $\sum ∣ a_{k} ∣∣ x ∣^{k}$ , and can be enormous when the terms cancel, as they do near a multiple root. Graillat, Langlois and Louvet's compensated Horner scheme computes the rounding error of every multiplication and addition exactly, with error-free transformations, and evaluates those errors with a second Horner recurrence; the result is as accurate as if Horner's method had been run in twice the working precision^[8]:

(x - 1)^7, expanded, evaluated near 1 by Horner's method and by the compensated scheme.

(* Evaluating (x - 1)^7, expanded, near its multiple root. The true
   values are tiny and the rounding errors are not. A compensated Horner
   scheme recovers the lost digits by computing each rounding error
   exactly and evaluating them alongside. *)
let p = [| 1.; -7.; 21.; -35.; 35.; -21.; 7.; -1. |]
 
let horner a x = Array.fold_left (fun acc c -> (acc *. x) +. c) 0. a
 
(* Error-free transformations: s + e = a + b and p + e = a * b exactly. *)
let two_sum a b = let s = a +. b in let z = s -. a in (s, (a -. (s -. z)) +. (b -. z))
let two_prod a b = let p = a *. b in (p, Float.fma a b (-. p))
 
let comp_horner a x =
  let s = ref a.(0) and c = ref 0. in
  for i = 1 to Array.length a - 1 do
    let pr, pe = two_prod !s x in
    let sm, se = two_sum pr a.(i) in
    s := sm;
    c := (!c *. x) +. (pe +. se)
  done;
  !s +. !c

Running it.

x             (x - 1)^7         Horner    compensated
0.9900    -1.000000e-14  -5.884182e-15  -1.000000e-14
0.9950    -7.812500e-17   1.776357e-15  -7.812500e-17
0.9990    -1.000000e-21   3.108624e-15  -1.000000e-21
1.0010     1.000000e-21  -1.998401e-15   1.000000e-21
1.0050     7.812500e-17   6.661338e-16   7.812500e-17
1.0100     1.000000e-14   7.993606e-15   1.000000e-14

Near $x = 1$ the true values are as small as $1 0^{- 21}$ , and the plain evaluation returns rounding noise of order $1 0^{- 15}$ , with the wrong sign in several cases. The compensated scheme agrees with $(x - 1)^{7}$ to all digits shown. The first column is itself computed in floating point from $x - 1$ , which is exact for these inputs up to the representation of $x$ .

§ 04

Parallel evaluation

Horner's recurrence is sequential: each step needs the previous one, so its latency is $n$ multiply-adds even on hardware that could do several at once. Estrin's scheme splits the polynomial into pairs and combines them with powers $x^{2}, x^{4}, x^{8}, \dots$ :

p (x) = (a_{0} + a_{1} x) + (a_{2} + a_{3} x) x^{2} + ((a_{4} + a_{5} x) + (a_{6} + a_{7} x) x^{2}) x^{4}

which has depth $O (lo g n)$ at the cost of a few extra multiplications^[9]. Vectorized implementations of elementary functions often use Estrin's scheme or a mixture of the two.

§ 05

In computer algebra

Over exact coefficient rings the same recurrence evaluates polynomials at integers or rationals, divides by linear factors, and substitutes one polynomial into another, with the $b_{k}$ computed by polynomial arithmetic. Multipoint evaluation at many points at once is faster with subproduct trees and FFT-based multiplication, and pseudo-division generalizes the division step to divisors of higher degree over rings that are not fields.

§ 06

History

Qin Jiushao described the method, as part of a procedure for extracting roots of polynomial equations, in his Mathematical Treatise in Nine Sections of 1247^[10]. Ruffini described synthetic division and a root-approximation method based on it in 1804^[2], and Horner published his method in the Philosophical Transactions in 1819^[1]; Cajori pointed out Ruffini's priority in 1911^[3]. Ostrowski's 1954 paper began the study of the complexity of polynomial evaluation^[4], and Pan settled the number of multiplications in 1966^[5].