Polynomial evaluation explained

In mathematics and computer science, polynomial evaluation refers to computation of the value of a polynomial when its indeterminates are substituted for some values. In other words, evaluating the polynomial

P(x _{1,} x _{2)} = 2 x _{1x} _{2} + x ^{3} _{1} + 4 at x _{1=2,} x _{2=3} consists of computing P( 2 , 3 ) = 2 ⋅ 2 ⋅ 3 + 2 ^{3+4=24.} See also a ^{n + a} _{n -1} x ^{n-1} + … + a _{0,} the most naive method would use n multiplications to compute a _{n} x ^{n} , use n -1 multiplications to compute a _{n-1} x ^{n-1} and so on for a total of \tfrac{n(n +1 )}{ 2 } multiplications and n additions.Using better methods, such as Horner's rule, this can be reduced to n multiplications and n additions. If some preprocessing is allowed, even more savings are possible.

Background

This problem arises frequently in practice. In computational geometry, polynomials are used to compute function approximations using Taylor polynomials. In cryptography and hash tables, polynomials are used to compute k-independent hashing.

In the former case, polynomials are evaluated using floating-point arithmetic, which is not exact. Thus different schemes for the evaluation will, in general, give slightly different answers. In the latter case, the polynomials are usually evaluated in a finite field, in which case the answers are always exact.

General methods

Horner's rule

See also: Horner's method.

Horner's method evaluates a polynomial using repeated bracketing:This method reduces the number of multiplications and additions to just

n

Horner's method is so common that a computer instruction "multiply–accumulate operation" has been added to many computer processors, which allow doing the addition and multiplication operations in one combined step.

Multivariate

If the polynomial is multivariate, Horner's rule can be applied recursively over some ordering of the variables.E.g.

P(x, y) = 4 + x + 2 x y + 2 x ^{2} y + x ^{2} y ^{2} can be written as \begin{align} P(x, y) & = 4 + x ( 1 + y( 2 ) + x (y( 2 + y))) or \\ P(x, y) & = 4 + x + y(x( 2 + x( 2 )) + y(x ^{2)). \end{align}}

An efficient version of this approach was described by Carnicer and Gasca.^[1]

Estrin's scheme

See also: Estrin's scheme.

While it's not possible to do less computation than Horner's rule (without preprocessing), on modern computers the order of evaluation can matter a lot for the computational efficiency.A method known as Estrin's scheme computes a (single variate) polynomial in a tree like pattern:

$$\backslash beginP(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.\backslash end$$

Combined by Exponentiation by squaring, this allows parallelizing the computation.

Evaluation with preprocessing

Arbitrary polynomials can be evaluated with feweroperations than Horner's rule requires if we first "preprocess"the coefficients

a _{n,} ... , a _{0} .

An example was first given by Motzkin^[2] who noted that

P(x) = x ^{4} + a _{3} x ^{3} + a _{2} x ^{2} + a _{1} x + a _{0} can be written asy = (x + \beta _{0)x+\beta} _{1,} P(x) = (y + x + \beta _{2)y+\beta} _{3,} where the values \beta _{0,} ... , \beta _{3} are computed in advanced, based on a _{0,} ... , a _{3} .Motzkin's method uses just 3 multiplications compared to Horner's 4.

The values for each

\beta _{i} can be easily computed by expanding P(x) and equating the coefficients:\begin{align} \beta _{0&=\tfrac12(a} _{3-1), &z&=a} _{2-\beta} _{0(\beta} _{0+1), &\beta} _{1&=a} _{1-\beta} _{0} z,\\ \beta _{2&=z-2\beta} _{1,} &\beta _{3&=a} _{0-\beta} _{1(\beta} _{1+\beta} _{2).\end{align}}

Example

\exp(x) ≈ 1+ x + x ^{2/2+x} ^{3/6+x} ^{4/24} ,we can upscale by a factor 24, apply the above steps, and scale back down.That gives us the three multiplication computationy = (x +1 . 5 )x +11 . 625 , P(x) = (y + x -15 )y/ 24+2 . 63477 . Improving over the equivalent Horner form (that is P(x) = 1 + x ( 1 + x ( 1 / 2 + x( 1 / 6 + x/ 24 ))) ) by 1 multiplication.

Some general methods include the Knuth–Eve algorithm and the Rabin–Winograd algorithm.

Multipoint evaluation

Evaluate of a

n -degree polynomial P(x) in multiple points x _{1,} ... , x _{m} can be done with mn multiplications by using Horner's method m times. Using above preprocessing approach, this can be reduced that by a factor of two, that is, to mn/ 2 multiplications. However, it is possible to do better.

It is possible to reduce the time requirement to just

O ( (n + m) log ^{2(n} + m) ) .^[3] The idea is to define two polynomials that are zero in respectively the first and second half of the points: m _{0(x)=(x-x} _{1) … (x-x} _{n/2} ) and m _{1(x)=(x-x} _{n/2+1} ) … (x - x _{n} ) .We then compute R _{0} = P \bmod m _{0} and R _{1} = P \bmod m _{1} using the Polynomial remainder theorem, which can be done in O(n log n) time using a fast Fourier transform.This means P(x) = Q(x)m _{0(x)} + R _{0(x)} and P(x) = Q(x)m _{1(x)} + R _{1(x)} by construction, where R _{0} and R _{1} are polynomials of degree at most n/ 2 .Because of how m _{0} and m _{1} were defined, we have\begin{align} R _{0(x} _{i)} & = P(x _{i)} for i \le n/ 2 and \\ R _{1(x} _{i)} & = P(x _{i)} for i > n/ 2 . \end{align} Thus to compute P on all n of the x _{i} , it suffices to compute the smaller polynomials R _{0} and R _{1} on each half of the points.This gives us a divide-and-conquer algorithm with T(n) = 2 T(n/ 2 ) + n log n , which implies T(n) = O(n( log n) ^{2)} by the master theorem.

In the case where the points in which we wish to evaluate the polynomials have some structure, simpler methods exist.For example, Knuth^[4] section 4.6.4gives a method for tabulating polynomial values of the type

P(x _{0} + h), P(x _{0} + 2 h), ... .

Dynamic evaluation

In the case where

x _{1,} ... , x _{m} are not known in advance,Kedlaya and Umans^[5] gave a data structure for evaluating polynomials over a finite field of size F _{q} in time ( log n) ^{O(1)} ( log _{2} q) ^{1+o(1)} per evaluation after some initial preprocessing.This was shown by Larsen^[6] to be essentially optimal.

The idea is to transform

P(x) of degree n into a multivariate polynomial f(x _{1,} x _{2,} ... , x _{m)} , such that P(x) = f(x, x ^{d,} x ^{d ^{2}} , ... , x ^{d ^{m}} ) and the individual degrees of f is at most d .Since this is over \bmod q , the largest value f can take (over Z ) is M = d ^{m} (q -1 ) ^{dm} .Using the Chinese remainder theorem, it suffices to evaluate f modulo different primes p _{1,} ... , p _{\ell} with a product at least M .Each prime can be taken to be roughly log M = O(dm log q) , and the number of primes needed, \ell , is roughly the same.Doing this process recursively, we can get the primes as small as loglog q .That means we can compute and store f on all the possible values in T = ( loglog q) ^{m} time and space.If we take d = log q , we get m = \tfrac{ log n}{ loglog q} , so the time/space requirement is just n ^{(loglog q)/(logloglog q)} .

Kedlaya and Umans further show how to combine this preprocessing with fast (FFT) multipoint evaluation.This allows optimal algorithms for many important algebraic problems, such as polynomial modular composition.

Specific polynomials

While general polynomials require

\Omega(n) operations to evaluate, some polynomials can be computed much faster.For example, the polynomial P(x) = x ^{2+2x+1} can be computed using just one multiplication and one addition since P(x) = (x +1 ) ^{2}

Evaluation of powers

See main article: Exponentiation by squaring and Addition-chain exponentiation.

A particularly interesting type of polynomial is powers like

x ^{n} .Such polynomials can always be computed in O( log n) operations.Suppose, for example, that we need to compute x ^{16} ; we could simply start with x and multiply by x to get x ^{2} .We can then multiply that by itself to get x ^{4} and so on to get x ^{8} and x ^{16} in just four multiplications.Other powers like x ^{5} can similarly be computed efficiently by first computing x ^{4} by 2 multiplications and then multiplying by x .

The most efficient way to compute a given power

x ^{n} is provided by addition-chain exponentiation. However, this requires designing a specific algorithm for each exponent, and the computation needed for designing these algorithms are difficult (NP-complete^[7]), so exponentiation by squaring is generally preferred for effective computations.

Polynomial families

Often polynomials show up in a different form than the well known

a _{n} x ^{n} + ... + a _{1} x + a _{0} .For polynomials in Chebyshev form we can use Clenshaw algorithm.For polynomials in Bézier form we can use De Casteljau's algorithm,and for B-splines there is De Boor's algorithm.

Hard polynomials

The fact that some polynomials can be computed significantly faster than "general polynomials" suggests the question: Can we give an example of a simple polynomial that cannot be computed in time much smaller than its degree?Volker Strassen has shown^[8] that the polynomial

P(x) = \sum ^{n} _{k =0} 2 ^{2 ^{kn ^{3}}} x ^{k} cannot be evaluated by with less than \tfrac 12 n - 2 multiplications and n - 4 additions.At least this bound holds if only operations of those types are allowed, giving rise to a so-called "polynomial chain of length <n ^{2/log} n ".

The polynomial given by Strassen has very large coefficients, but by probabilistic methods, one can show there must exist even polynomials with coefficients just 0's and 1's such that the evaluation requires at least

\Omega(n/ log n) multiplications.

For other simple polynomials, the complexity is unknown.The polynomial

(x +1 )(x +2 ) … (x + n) is conjectured to not be computable in time ( log n) ^{c} for any c .This is supported by the fact, that if it can be computed fast then integer factorization can be computed in polynomial time, breaking the RSA cryptosystem.^[9]

Matrix polynomials

Sometimes the computational cost of scalar multiplications (like

ax ) is less than the computational cost of "non scalar" multiplications (like x ^{2} ).The typical example of this is matrices.If M is an m x m matrix, a scalar multiplication aM takes about m ^{2} arithmetic operations, while computing M ^{2} takes about m ^{3} (or m ^{2.3} using fast matrix multiplication).

Matrix polynomials are important for example for computing the Matrix Exponential.

Paterson and Stockmeyer ^[10] showed how to compute a degree

n polynomial using only O(\sqrt n) non scalar multiplications and O(n) scalar multiplications.Thus a matrix polynomial of degree can be evaluated in O(m ^{2.3} \sqrt{n} + m ^{2n)} time. If m = n this is O(m ^{3)} , as fast as one matrix multiplication with the standard algorithm.

This method works as follows: For a polynomial

P(M) = a _{n-1} M ^{n-1} + ... + a _{1} M + a _{0} I, let be the least integer not smaller than \sqrt{n}. The powers M, M ^{2,} ... , M ^{k} are computed with k matrix multiplications, and M ^{2k} , M ^{3k} , ... , M ^{k ^{2-k}} are then computed by repeated multiplication by M ^{k.} Now,\begin{align}P(M) = &(a _{0} I + a _{1} M + ... + a _{k-1} M ^{k-1} ) \\ + &(a _{k} I + a _{k+1} M + ... + a _{2k-1} M ^{k-1} )M ^{k \\ + & ... \\ + &(a} _{n - k} I + a _{n-k+1} M + ... + a _{n-1} M ^{k-1} )M ^{k ^{2-k}} , \end{align} ,where a _{i=0} for .This requires just k more non-scalar multiplications.

We can write this succinctly using the Kronecker product:

P(M) = \begin{bmatrix}I\\M\\\vdots\\M ^{k-1} \end{bmatrix} ^{T \left(\begin{bmatrix} a} _{0} & a _{1} & a _{2} & ... \\ a _{k} & a _{k+1} & \ddots \\ a _{2k} & \ddots \\ \vdots\end{bmatrix} ⊗ I\right) \begin{bmatrix}I\\M ^{k\\M} ^{2k} \\\vdots\end{bmatrix} .

The direct application of this method uses

2 \sqrt{n} non-scalar multiplications, but combining it with Evaluation with preprocessing, Paterson and Stockmeyer show you can reduce this to \sqrt{ 2 n} .

Methods based on matrix polynomial multiplications and additions have been proposed allowing to save nonscalar matrix multiplications with respect to the Paterson-Stockmeyer method.^[11]

See also

• Estrin's scheme to facilitate parallelization on modern computer architectures
• Arithmetic circuit complexity theory studies the computational complexity of evaluating different polynomials.

Notes and References

1. Carnicer . J. . Gasca . M. . 1990 . Evaluation of Multivariate Polynomials and Their Derivatives . . 54 . 189 . 231–243 . 10.2307/2008692 . 2008692 . free.
2. Motzkin. T. S.. 1955. Evaluation of polynomials and evaluation of rational functions. Bulletin of the American Mathematical Society. 61. 163. 10.
3. Book: Von Zur Gathen . Joachim . Modern computer algebra . Jürgen . Gerhard . . 2013 . 9781139856065 . Chapter 10 . Joachim von zur Gathen.
4. Book: Knuth, Donald . . . 2005 . 9780201853926 . 2: Seminumerical Algorithms . Donald Knuth.
5. Kedlaya. Kiran S.. Kiran Kedlaya. Umans. Christopher. Chris Umans. 2011. Fast Polynomial Factorization and Modular Composition. SIAM Journal on Computing. 40. 6. 1767–1802. 10.1137/08073408x. 1721.1/71792. 412751 . free.
6. Book: Larsen, K. G.. 2012 IEEE 53rd Annual Symposium on Foundations of Computer Science . Higher Cell Probe Lower Bounds for Evaluating Polynomials . Kasper Green Larsen. 2012. IEEE. 53. 293–301. 10.1109/FOCS.2012.21. 978-0-7695-4874-6 . 7906483 .
7. Downey . Peter . Leong . Benton . Sethi . Ravi . Computing Sequences with Addition Chains . SIAM Journal on Computing . 1981 . 10 . 3 . 27 January 2024.
8. Strassen. Volker. Volker Strassen. 1974. Polynomials with Rational Coefficients Which are Hard to Compute. SIAM Journal on Computing. en. 3. 2. 128–149. 10.1137/0203010.
9. Chen, Xi, Neeraj Kayal, and Avi Wigderson. Partial derivatives in arithmetic complexity and beyond. Now Publishers Inc, 2011.
10. Paterson. Michael S.. Mike Paterson. Stockmeyer. Larry J.. Larry J. Stockmeyer. 1973. On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials. SIAM Journal on Computing. en. 2. 1. 60–66. 10.1137/0202007.
11. Fasi . Massimiliano . Optimality of the Paterson–Stockmeyer method for evaluating matrix polynomials and rational matrix functions . Linear Algebra and its Applications . 1 August 2019 . 574 . 185 . 10.1016/j.laa.2019.04.001 . 0024-3795. free .