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:

can be written in a nested form as:

or nested in a different order as:

and in a factored form as:

where QdQd is the multiplicity of the dthdth zero, DD is the number of distinct zeros, and ∑d=1DQd=N d 1 D Qd N .

The polynomial can be written in a first order remainder form as:

where

is the quotient polynomial obtained when fNaz fN a z is divided by z-z0zz0 and the remainder is a constant easily seen to be the value

If z0z0 is a zero of fNz fN z , then R=0R0 and q N - 1 z q N - 1 z contains the same zeros as fNz fN z except for the one at z0z0 . A more general formulation of (5) is:

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:

for k= 0,1,..., N-1 k 0,1,..., N-1 with the value of the polynomial at z given by fNaz=xn fN a z xn [3,4,5]. Here the xkxk are the values of the successively evaluated bracketed epressions in (2). Matlab code to evaluate fNaz0 fN a z0 with coefficients aa 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

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:

which, if evaluated at z=z0zz0 , gives:

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 qN-1 az is not the derivative, but is the value of the derivative at z0z0 . 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

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);

From (11) and (6), we see that the value of the derivative of the Nth degree polynomial fNaz fN az is the value of the N-1thN1th degree polynomial qN-1bz qN-1 bz with coefficients bkbk as solutions to the difference equation

which has the same form as (9) but saving the intermediate values of bkbk . This means that the solution to the difference equation (12) with the NN input values of akak gives N-1N1 output values of bkbk followed by the remainder R1R1 which is the value of fNaz fN az .

A similar argument shows that solving (12) with an input of bkbk will give N-2N2 output values of ckck followed by R2R2 which is the value of the derivative of fNaz fN az . Using the ckck as an input will give N-3N3 values of dkdk and R3R3 which is the second derivative. This can be continued NN times to evaluate the function and all of its N-1N1 derivatives (see Fig 1). These values of RkRk allow the original polynomial to be written as:

which is a power series expansion [5] around z=z0 z z0 with the coefficients RkRk being defined in terms of derivatives evaluated at z=z0 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 bkbk or ckck .

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;

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 zz and shows the structure of evaluating derivatives.

Figure 5: Flow Graph for Horner's Algorithm

Polynomial Deflation

If z0z0 is a zero of the polynomial fNaz fN a z , then R=0R0 and (5) becomes

with qN-1qN-1 being the N-1N1 degree polynomial (6) having the same zeros at fNfN except for the one at z=z0 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 fNaz 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

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);

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 akak by a divisor with coefficients dkdk to give a quotient with coefficients qkqk using the FFT is:

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

A more general formulaiton of the deflation problem is of the form:

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 z 1 and unstable for |z|>1 z 1 .

To reduce error accumulation when |z|>1 z 1 we nest (1) and (2) in still a different way:

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}

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);

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 zz 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