Horner's method

In mathematics, Horner's method (also known as Horner scheme in the UK or Horner's rule in the U.S.) is either of two things: (i) an algorithm for calculating polynomials, which consists in transforming the monomial form into a computationally efficient form; or (ii) a method for approximating the roots of a polynomial. The latter is also known as Ruffini–Horner's method.

These methods are named after the British mathematician William George Horner. It is now established that both methods were known before him.

Description of the algorithm

Given the polynomial

$p(x) = \sum_{i=0}^n a_i x^i = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n,$

where $a_0, \ldots, a_n$ are real numbers, we wish to evaluate the polynomial at a specific value of x, say x₀.

To accomplish this, we define a new sequence of constants as follows:

$\begin{align} b_n & := a_n \\ b_{n-1} & := a_{n-1} + b_n x_0 \\ & {}\ \ \vdots \\ b_0 & := a_0 + b_1 x_0. \end{align}$

Then b₀ is the value of p(x₀).

To see why this works, note that the polynomial can be written in the form

$p(x) = a_0 + x(a_1 + x(a_2 + \cdots + x(a_{n-1} + a_n x)+\cdots)). \,$

Thus, by iteratively substituting the $b_i$ into the expression,

$\begin{align} p(x_0) & = a_0 + x_0(a_1 + x_0(a_2 + \cdots + x_0(a_{n-1} + b_n x_0)+ \cdots)) \\ & = a_0 + x_0(a_1 + x_0(a_2 + \cdots + x_0(b_{n-1})+ \cdots)) \\ & {} \ \ \vdots \\ & = a_0 + x_0(b_1) \\ & = b_0. \end{align}$

Examples

Evaluate

$f(x)=2x^3-6x^2+2x-1\,$ for $x=3\;$.

We use synthetic division as follows:

 x₀│   x³    x²    x¹    x⁰
 3 │   2    −6     2    −1
   │         6     0     6    
   └────────────────────────
       2     0     2     5

The entries in the third row are the sum of those in the first two. Each entry in the second row is the product of the x-value (3 in this example) with the third-row entry immediately to the left. The entries in the first row are the coefficients of the polynomial to be evaluated. Then the remainder of $f(x)$ on division by $x-3$ is 5.

But by the remainder theorem, we know that the remainder is $f(3)$. Thus $f(3) = 5$

In this example, if $a_3 = 2, a_2 = -6, a_1 = 2, a_0 = -1$ we can see that $b_3 = 2, b_2 = 0, b_1 = 2, b_0 = 5$, the entries in the third row. So, synthetic division is based on Horner's method.

As a consequence of the polynomial remainder theorem, the entries in the third row are the coefficients of the second-degree polynomial, the quotient of $f(x)$ on division by $x-3$. The remainder is 5. This makes Horner's method useful for polynomial long division.

Divide $x^3-6x^2+11x-6\,$ by $x-2\,$:

 2 │   1    -6    11    -6
   │         2    -8     6    
   └────────────────────────
       1    -4     3     0

The quotient is $x^2-4x+3\,$.

Let $f_1(x)=4x^4-6x^3+3x-5\,$ and $f_2(x)=2x-1\,$. Divide $f_1(x)\,$ by $f_2\,(x)$ using Horner's method.

  2 │  4    -6    0    3   │   -5
────┼──────────────────────┼───────
  1 │        2   -2   -1   │    1
    │                      │  
    └──────────────────────┼───────
       2    -2    -1   1   │   -4    

The third row is the sum of the first two rows, divided by 2. Each entry in the second row is the product of 1 with the third-row entry to the left. The answer is

$\frac{f_1(x)}{f_2(x)}=2x^3-2x^2-x+1-\frac{4}{2x-1}.$

Floating point multiplication and division

Main article: multiplication algorithm#Shift and add

Horner's method is a fast, code-efficient method for multiplication and division of binary numbers on a microcontroller with no hardware multiplier. One of the binary numbers to be multiplied is represented as a trivial polynomial, where, (using the above notation): a_i = 1, and x = 2. Then, x (or x to some power) is repeatedly factored out. In this binary numeral system (base 2), x = 2, so powers of 2 are repeatedly factored out.

Example

For example, to find the product of two numbers, (0.15625) and m:

$\begin{align} ( 0.15625) m & = (0.00101_b) m = ( 2^{-3} + 2^{-5}) m = (2^{-3})m + (2^{-5})m \\[4pt] & = 2^{-3} (m + (2^{-2})m) = 2^{-3} (m + 2^{-2} (m)). \end{align}$

Method

To find the product of two binary numbers, "d" and "m".

• 1. A register holding the intermediate result is initialized to (d).
• 2. Begin in (m) with the least significant (rightmost) non-zero bit,
  • 2b. Count (to the left) the number of bit positions to the next most significant non-zero bit. If there are no more-significant bits, then take the value of the current bit position.
  • 2c. Using that value, perform a right-shift operation by that number of bits on the register holding the intermediate result
• 3. If all the non-zero bits were counted, then the intermediate result register now holds the final result. Otherwise, add (d) to the intermediate result, and continue in step #2 with the next most significant bit in (m).

Derivation

In general, for a binary number with bit values: ($d_3 d_2 d_1 d_0$) the product is:

$(d_3 2^3 + d_2 2^2 + d_1 2^1 + d_0 2^0)m = d_3 2^3 m + d_2 2^2 m + d_1 2^1 m + d_0 2^0 m$

At this stage in the algorithm, it is required that terms with zero-valued coefficients are dropped, so that only binary coefficients equal to one are counted, thus the problem of multiplication or division by zero is not an issue, despite this implication in the factored equation:

$= d_0(m + 2 \frac{d_1}{d_0} (m + 2 \frac{d_2}{d_1} (m + 2 \frac{d_3}{d_2} (m))))$

The denominators all equal one (or the term is absent), so this reduces to:

$= d_0(m + 2 {d_1} (m + 2 {d_2} (m + 2 {d_3} (m))))$

or equivalently (as consistent with the "method" described above):

$= d_3(m + 2^{-1} {d_2} (m + 2^{-1}{d_1} (m + {d_0} (m))))$

In binary (base 2) math, multiplication by a power of 2 is merely a register shift operation. Thus, multiplying by 2 is calculated in base-2 by an arithmetic shift. The factor (2⁻¹) is a right arithmetic shift, a (0) results in no operation (since 2⁰ = 1, is the multiplicative identity element), and a (2¹) results in a left arithmetic shift. The multiplication product can now be quickly calculated using only arithmetic shift operations, addition and subtraction.

The method is particularly fast on processors supporting a single-instruction shift-and-addition-accumulate. Compared to a C floating-point library, Horner's method sacrifices some accuracy, however it is nominally 13 times faster (16 times faster when the "canonical signed digit" (CSD) form is used), and uses only 20% of the code space.

Polynomial root finding

Using Horner's method in combination with Newton's method, it is possible to approximate the real roots of a polynomial. The algorithm works as follows. Given a polynomial $p_n(x)$ of degree $n$ with zeros $z_n < z_{n-1} < \cdots < z_1$, make some initial guess $x_0$ such that $x_0 > z_1$. Now iterate the following two steps:

1. Using Newton's method, find the largest zero $z_1$ of $p_n(x)$ using the guess $x_0$.

2. Using Horner's method, divide out $(x-z_1)$ to obtain $p_{n-1}$. Return to step 1 but use the polynomial $p_{n-1}$ and the initial guess $z_1$.

These two steps are repeated until all real zeros are found for the polynomial. If the approximated zeros are not precise enough, the obtained values can be used as initial guesses for Newton's method but using the full polynomial rather than the reduced polynomials.

Example

Polynomial root finding using Horner's method

Consider the polynomial,

$p_6(x) = (x-3)(x+3)(x+5)(x+8)(x-2)(x-7)$

which can be expanded to

$p_6(x) = x^6 + 4x^5 - 72x^4 -214x^3 + 1127x^2 + 1602x -5040.$

From the above we know that the largest root of this polynomial is 7 so we are able to make an initial guess of 8. Using Newton's method the first zero of 7 is found as shown in black in the figure to the right. Next $p(x)$ is divided by $(x-7)$ to obtain,

$p_5(x) = x^5 + 11x^4 + 5x^3 - 179x^2 - 126x + 720 \,$

which is drawn in red in the figure to the right. Newton's method is used to find the largest zero of this polynomial with an initial guess of 7. The largest zero of this polynomial which corresponds to the second largest zero of the original polynomial is found at 3 and is circled in red. The degree 5 polynomial is now divided by $(x-3)$ to obtain,

$p_4(x) = x^4 + 14x^3 + 47x^2 - 38x - 240 \,$

which is shown in yellow. The zero for this polynomial is found at 2 again using Newton's method and is circled in yellow. Horner's method is now used to obtain,

$p_3(x) = x^3 + 16x^2 + 79x + 120 \,$

which is shown in green and found to have a zero at −3. This polynomial is further reduced to

$p_2(x) = x^2 + 13x + 40 \,$

which is shown in blue and yields a zero of −5. The final root of the original polynomial may be found by either using the final zero as an initial guess for Newton's method, or by reducing $p_2(x)$ and solving the linear equation. As can be seen, the expected roots of −8, −5, −3, 2, 3, and 7 were found.

Octave implementation

The following Octave code was used in the example above to implement Horner's method.

function [y b] = horner(a,x)
  % Input a is the polynomial coefficient vector, x the value to be evaluated at.
  % The output y is the evaluated polynomial and b the divided coefficient vector.
  b(1) = a(1);
  for i = 2:length(a)
    b(i) = a(i)+x*b(i-1);
  end
  y = b(length(a));
  b = b(1:length(b)-1);
end

Python implementation

The following Python code implements Horner's method.

def horner(x, *args):
    """A function that implements the Horner Scheme for evaluating a
    polynomial of coefficients *args in x."""
    result = 0
    for coefficient in args:
        result = result * x + coefficient
    return result

Application

Horner's method can be used to convert between different positional numeral systems – in which case x is the base of the number system, and the a_i coefficients are the digits of the base-x representation of a given number – and can also be used if x is a matrix, in which case the gain in computational efficiency is even greater. In fact, when x is a matrix, further acceleration is possible which exploits the structure of matrix multiplication, and only $\sqrt{n}$ instead of n multiplies are needed (at the expense of requiring more storage) using the 1973 method of Paterson and Stockmeyer.

Efficiency

Evaluation using the monomial form of a degree-n polynomial requires at most n additions and (n² + n)/2 multiplications, if powers are calculated by repeated multiplication and each monomial is evaluated individually. (This can be reduced to n additions and 2n − 1 multiplications by evaluating the powers of x iteratively.) If numerical data are represented in terms of digits (or bits), then the naive algorithm also entails storing approximately 2n times the number of bits of x (the evaluated polynomial has approximate magnitude xⁿ, and one must also store xⁿ itself). By contrast, Horner's method requires only n additions and n multiplications, and its storage requirements are only n times the number of bits of x. Alternatively, Horner's method can be computed with n fused multiply–adds. Horner's method can also be extended to evaluate the first k derivatives of the polynomial with kn additions and multiplications.

Horner's method is optimal, in the sense that any algorithm to evaluate an arbitrary polynomial must use at least as many operations. Alexander Ostrowski proved in 1954 that the number of additions required is minimal. Victor Pan proved in 1966 that the number of multiplications is minimal. However, when x is a matrix, Horner's method is not optimal.

This assumes that the polynomial is evaluated in monomial form and no preconditioning of the representation is allowed, which makes sense if the polynomial is evaluated only once. However, if preconditioning is allowed and the polynomial is to be evaluated many times, then faster algorithms are possible. They involve a transformation of the representation of the polynomial. In general, a degree-n polynomial can be evaluated using only ${\scriptstyle{\left\lfloor n/2 \right\rfloor + 2}}$ multiplications and n additions (see Knuth: The Art of Computer Programming, Vol.2).

History

Qin Jiushao algorithm for solving quadratic polymial equation$-x^4+763200x^2-40642560000=0$
result:x=840

According to Augustus De Morgan, Horner made his method known in a paper titled "A new method of solving numerical equations of all orders, by continuous approximation", which was read on July 1, 1819 by Davies Gilbert before the Royal Society of London. This paper was subsequently published in the same year in the "Philosophical Transactions".

De Morgan's advocacy of Horner's priority in discovery led to "Horner's method" being so called in textbooks. However, this method was known long before. In reverse chronological order, Horner's method was already known to:

• Paolo Ruffini in 1809 (see Ruffini's rule)
• Isaac Newton in 1669
• the Chinese mathematician Zhu Shijie in the 14th century
• the Chinese mathematician Qin Jiushao in his Mathematical Treatise in Nine Sections in the 13th century
• the Persian mathematician Sharaf al-Dīn al-Tūsī in the 12th century
• the Chinese mathematician Jia Xian in the 11th century (Song Dynasty)
• The Nine Chapters on the Mathematical Art, a Chinese work of the Han Dynasty (202 BC – 220 AD) edited by Liu Hui (fl. 3rd century).

See also

• Clenshaw algorithm to evaluate polynomials in Chebyshev form
• De Boor's algorithm to evaluate splines in B-spline form
• De Casteljau's algorithm to evaluate polynomials in Bézier form
• Estrin's scheme to facilitate parallelization on modern computer architectures
• Lill's method to approximate roots graphically
• Ruffini's rule to divide a polynomial by a binomial of the form x − r

References

Bibliography

• Horner, William George (July 1819). "A new method of solving numerical equations of all orders, by continuous approximation". Philosophical Transactions (Royal Society of London): pp. 308–335. 
• Spiegel, Murray R. (1956). Schaum's Outline of Theory and Problems of College Algebra. McGraw-Hill Book Company. 
• Knuth, Donald (1997). The Art of Computer Programming. Vol. 2: Seminumerical Algorithms (3rd ed.). Addison-Wesley. pp. 486–488 in section 4.6.4. ISBN 0-201-89684-2. 
• Kripasagar, Venkat (March 2008). "Efficient Micro Mathematics – Multiplication and Division Techniques for MCUs". Circuit Cellar magazine (212): p. 60. 

External links

• Weisstein, Eric W., "Horner's method" from MathWorld.
• Module for Horner's Method by John H. Mathews