Polynomials Polynomials are one of the oldest classes of mathematical functions with an important history in mathematics, science, engineering, and other quantitative fields [1]. Part of the polynomial's appeal comes from the fact that it may be numerically evaluated using a finite number of multiplications and additions. In the processes of factoring, evaluating, and deflating polynomials, Horner's methods are central and are the focus of this note.Any Nth degree polynomial can be written in coefficient form as: fN a z a0 a1 z a2 z2 ... aN zN can be written in a nested form as: fN a z a0 z a1 z a2 ... z a N - 1 z a N ... or nested in a different order as: fN a z z ... z z z aN a N - 1 a N - 2 ... a0 and in a factored form as: f z a N d 1 D z z d Q d where Qd is the multiplicity of the dth zero, D is the number of distinct zeros, and d 1 D Qd N . The polynomial can be written in a first order remainder form as: fN a z qN-1 b z z z0 R where qN-1 b z b0 b1 z b2 z2 ... b N - 1 zN1 is the quotient polynomial obtained when fN a z is divided by zz0 and the remainder is a constant easily seen to be the value R fN az0 If z0 is a zero of fN z , then R0 and q N - 1 z contains the same zeros as fN z except for the one at z0 . A more general formulation of (5) is: fN az q N - M bz dM cz R z There are four polynomial problems that are often posed: Factor: Given the coefficients ak , find the roots or zeros zr . Matlab command: roots() Unfactor: Given the zeros zr , find the coefficients ak . Matlab command: poly() Evaluate: Given the coefficients or zeros and z0 , find the value of the polynomial at z z0 , f z0 (or the value of the derivative of the polynomial at z z0 , f z0 . Matlab command: polyval() Deflate: Given a polynomial and one of its roots, find the polynomial containing the other roots. Horner's methods are important for evaluation and deflation, therefore, for factoring. For many high degree polynomial factoring schemes, it is important to use stable evaluation and deflation and to deflate in an order that maximizes the conditioning of the quotient. Unfactoring is simply the multiplying of the factors to obtain the polynomial but, because of round-off error, the order of multiplication is important. As an alternative, the FFT/DFT can be used to unfactor and evaluate. Of the four, only factoring is nonlinear and it is the most difficult. In many strategies, factoring uses the other three. For implementation of any of the four, the number of arithmetic operations required will determine the speed of execution and stability and conditioning will affect the error characteristics.

Evaluate the Polynomial and Its Derivatives From the structure of the nested forms in (2) and (3), one can formulate a recursive relation which is also a linear first order difference equation: x k1 z xk a N-1-k , x0 an for k 0,1,..., N-1 with the value of the polynomial at z given by fN a z xn [3,4,5]. Here the xk are the values of the successively evaluated bracketed epressions in (2). Matlab code to evaluate fN a z0 with coefficients a is given by:


   m = length(a);     % poly degree plus one: N+1
   f = a(1);          % initial condition
   for k = 2:m        % iterative Horner's algorithm
     f=z*f + a(k);    % recursive evaluation of f(z)
   end

Program 1. Forward Evaluation of Polynomial Remember that Matlab stores polynomial coefficients in the reverse order of the notation used by many writers and the reverse of that shown in (1) and starts indexing with 1 rather than 0. This formulation is very efficient if implemented with the "multiplicity-accumulate" single-cycle instruction of several DSP chips and some microprocessors. The nested forms of (2) and (3), the remainder form of (5), and the recursive from of (9) are all versions of Horner's method [2,5,6] or synthetic division. There are also similarities to the Goertzel algoirthm to evaluate the DFT, the Padé or Prony method and inversion of z-transforms by division. Indeed, all of these methods convert an evaluation into a recursion and then into a difference equation which has considerable literature [7]. While in some cases it is possible to reduce the number of multiplications used by Horner in polynomial evaluation, it is generally a small gain and accompanied by an increase in the number of additions and the complexity of the algorithm. If preprocessing of the polynomial coefficients is not used, Horner gives a minimum arithmetic evaluation [3,5]. If we differentiate (5), we get: fN az qN-1 bz qN-1 bz z z0 which, if evaluated at zz0 , gives: fN az0 qN-1 bz0 which can be evaluated by a minor variation of (9). Thus Horner's method can evaluate an Nth degree polynomial with N multiplications and additions or evaluate it and its derivative with 2N multiplications and additions or avaluate it and its derivative with 2N multiplications and additions. Note qN-1 az is not the derivative, but is the value of the derivative at z0 . The simple Matlab code for this is:



m = length(a);          % N+1
f = a(1); fp = 0;       % initial conditions
for i = 2:m             % iterative Horner's algorithms
    fp = z * fp + f;    % recursive evaluation of f'(z)
    f = z * f + a(i);   % recursive evaluation of f(z)
end

Program 2. Forward Evaluation of Polynomial and Its Derivative This can be made much faster by using the filter() function in Matlab as:


N = length(a) - 1;
b = filter(1,[1 -z], a);       % Horner by filter: f(z)
f = b(N+1);
c = filter(1,[1 -z], b(1:N));  % Horner by filter: f'(z)
fp = c(N);

Program 3. Forward Evaluation of Polynomial and Its Derivative by Filters From (11) and (6), we see that the value of the derivative of the Nth degree polynomial fN az is the value of the N1th degree polynomial qN-1 bz with coefficients bk as solutions to the difference equation b k1 z bk a N-1-k which has the same form as (9) but saving the intermediate values of bk . This means that the solution to the difference equation (12) with the N input values of ak gives N1 output values of bk followed by the remainder R1 which is the value of fN az . A similar argument shows that solving (12) with an input of bk will give N2 output values of ck followed by R2 which is the value of the derivative of fN az . Using the ck as an input will give N3 values of dk and R3 which is the second derivative. This can be continued N times to evaluate the function and all of its N1 derivatives (see Fig 1). These values of Rk allow the original polynomial to be written as: fN a z R1 R2 zz0 R3 zz0 2 ... R N1 zz0 N which is a power series expansion [5] around z z0 with the coefficients Rk being defined in terms of derivatives evaluated at z z0 . The following program evaluates the function and its first and second derivatives in a very compact single loop solution of the basic first order difference equation without saving the coefficients bk or ck .


m = length(a);              % poly degree + 1
f = 0; fp = 0; fpp = 0;     % initial values
for i = 1:m-2               % Horner's algorithm
    f = z*f + a(i);         % recursive evaluation of f(z)
    fp = z*fp + f;          % recursive evaluation of f'(z)
    fpp = z*fpp + fp;       % recursive evaluation of f''(z)
end
f = z*f + a(m-1);
fp = z*fp + f;
f = z*f + a(m);             % Finish evaluations
fpp = 2*fpp;

Program 4. Forward Evaluation of Polynomial and Two Derivatives Higher order derivatives can be evaluated by a simple extension of this idea. Describing the algorithm as a flow-graph illustrates the recursion as a feedback path with the feedback parameter being z and shows the structure of evaluating derivatives.

Flow Graph for Horner's Algorithm

Polynomial Deflation If z0 is a zero of the polynomial fN a z , then R0 and (5) becomes fN az qN-1 bz z z0 with qN-1 being the N1 degree polynomial (6) having the same zeros at fN except for the one at z z0 . In other words, Horner's method can also be used to deflate a polynomial with a known zero. A Matlab program that will deflate fN az is given by:


m = length(a);               % order = one: N+1
b = zeros(1, m-1);
b(1) = a(1);                 % initial conditions
for k = 1:m-2                % iterative Horner's algorithm 
  b(k+1) = z*b(k) + a(k+1);  % recursive deflation of f(z)
end

Program 5. Forward Deflation of Polynomial or using the more efficient Matlab "filter" command:


N = length(a) - 1;
b = filter(1,[1 -z, a);   % Deflation by filter
b = b(1:N);

Program 6. Forward Deflation of Polynomial using Filter Because multiplication of two polynomials is the same operation as the convolution of their coefficients, one can get the same results by multiplying the discrete Fourier transform (DFT)'s of the polynomial coefficients and taking the inverse DFT of the product. This gives an alternative to Horner's algorithm for deflating polynomials and is the technique used in the Lindsey-Fox polynomial factoring scheme [8, 2]. The use of the FFT for deflation or unfactoring of polynomials has very different error characteristics from Horner's method and is often superior. A Matlab program that deflates a polynomial with coefficients ak by a divisor with coefficients dk to give a quotient with coefficients qk using the FFT is:


m = length(a); md = length(d);
q = ifft(ifft(a)./fft([d,zeros(1,(m-md))]));

Program 7. Deflation of Polynomial usng the FFT A more general formulaiton of the deflation problem is of the form: fN az qN-M bz dM cz R z If the factors of dM z are all roots of fN z , then q NM contains the others when R z 0 . This allows the easy deflation by quadratic factors or other factors that are already known.

Stability Horner's method is implemented by the recursive equation (9) or (12) which, in this form, is seen to be a linear first order constant coefficient difference equation. There is considerable literature on the stability of linear differential and difference equations [7], and an equation of this form is known to be stable [9] for z 1 and unstable for z 1 . To reduce error accumulation when z 1 we nest (1) and (2) in still a different way: fN a z zN z-1 ... z-1 z-1 z-1 a0 a1 ... a N1 aN to give an alternate recursive relationship to (9) which is the following difference equation: x k1 z-1 xk a k1 , x0 a0 for k 0,1,..., N-1 with the value of the polynomial at z given by fN a z zN xN . This equation is stable [7,9] for z 1 . Again, remembering the Matlab index convention, the program for this form of Horner's algorithm is:


m = length(a);        % N+1
f = a(m);             % initial condition
for k = 1:m-1         % iterative Horner's method 
  f = f/z + a(m-k);   % recursive evaluation of f(z)
end
f = (z^(m-1))*f;      % remove the factor of z^{-m+1}

Program 8. Backward Evaluation of Polynomial which is stable for |z|>1 The same can be done for deflation.This algorithm is equivalent to reversing the order ("flipping") the polynomial coefficients and applying (9) which will give zeros that are the reciprocal of those of the original polynomial [6,2]. Using the standard Matlab commands, this can be done by:


m = length(a);
f = (z^(m-1))*polyval(a(m:-1:1),1/z);

Program 9. Backward Evaluation of Polynomial by using Reversed Order Coefficients It is easy to see that stability is very important when evaluating or deflating high degee polynomials by considering the effects if the degree is as high as one million [2]. Using this difference equation formulation of Horner's methods allows a more general deflation using a quadratic or higher degree divisor polynomial based on (8) and (15) to be easily analysed for stability.

Conclusions We have seen that polynomial evaluation and deflation can be done by Horner's approach which is the conversion of the problem into a linear difference equation. Indeed, evaluation and first order deflation are accomplished by the same algorithm and the stability is easily controlled. Horner's methods are used in several basic Newton type factoring programs such as "broots". The alternative use of the DFT (FFT) for deflation should always be considered. Having forwards and backwards algorithms allows one to always use a stable evaluation or deflation scheme [10,11,5,12,6,2,9]. If the magnitude of z where the polynomial is being evaluated (or where the polynomial is being deflated) is less than one, use (9), if it is greater than one, use (17).

References V. Y. Pan. Solving a polynomial equation: some history and recent progress. SIAM Review, 39(2):187–220, June 1997. Gary A. Sitton, C. Sidney Burrus, James W. Fox, and Sven Treitel. Factoring very high degree polynomials. IEEE Signal Processing Magazine, 20(6):27–42, November 2003. Donald E. Knuth. The Art of Computer Programming, Vol. 2, Seminumerical Algorithms. Addison-Wesley, Reading, MA, third edition, 1997. E. J. Barbeau. Polynomials. Springer-Verlag, New York, 1989. N. J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM, 1996. Chapter 5 on polynomials. W. H. Press, B. P.Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in Fortran: The Art of Scientiﬁc Computing. Cambridge University Press, Cambridge, second edition, 1992. Saber N. Elaydi.An Introduction to Diﬀerence Equations. Springer-Verlag, New York, second edition, 1999. First edition in 1996. J. P. Lindsey and James W. Fox. A method of factoring long z-transform polynomials. Computational Methods in Geosciences, SIAM, 78–90, 1992. Reprinted in Seismic Source Signature Estimation and Measurement, edited by Osman Osman and Enders Robinson, Society of Exploration Geophysicists, Geophysics Reprint Series no. 18, 712–724, 1996. James W. Fox. Improving the accuracy of polynomial evaluation, deﬂation, and root ﬁnding. draft manuscript, January 14 2002. Daniela Calvetti and Lothar Reichel. On the evaluation of polynomial coeﬃcients. Numerical Algorithms, 33:153–161, 2003. J. H. Wilkinson.Rounding Errors in Algebraic Processes. Prentice-Hall, Englewood Cliﬀs, NJ, 1963. Now published by Dover, 1994. Forman S. Acton. Numerical Methods that Work. The Mathematical Association of America, Washington, DC, 1990.