Monday, November 10, 2014

Evaluating polynomials at regular-space grids

Remember the successive differences approach we learned in school to analyze series of numbers? Consider the set of 4th powers: 0, 1, 16, 81, 256, 625, 1296, 2401, 4096, 6561, ...
Finding the differences between successive terms we get another sequence:
1, 15, 65, 175, 369, 671, 1105, 1695, 2465, ...
Repeating this procedure, we eventually get:
0,    1,  16,  81,  256,  625,  1296,  2401,  4096,    6561, ...
    1,  15,  65,  175, 369,   671,  1105,   1695,   2465, ...
       14, 50, 110,  194,  302,   434,   590,     770, ...
          36,  60,    84,   108,   132,  156,     180, ...
              24,   24,    24,     24,     24,     24, ...

Each row is obtained by adding the numbers in the row below to obtain the next term. We notice that the fifth row is constant, i.e. we can obtain the fourth row by simply adding a constant to a previous term. This shows that we can compute each row simply by addition.  This approach is valid for any polynomial and is the essence of Nuttall's algorithm (Ref. [1]) to evaluate polynomials on a equally spaced grid. Consider the 3rd order polynomial $p(x) = x^3$ and we want to evaluate it at $x = x_0+n\cdot dx$ for $n=0,1,2, ...$. For simplicity (and without loss of generality, as we can always translate and scale the polynomial), assume that $x_0 = 0$ and $dx = 1$.  Using successive differences, we obtain:
$p(n) = p(n-1) + d_2(n)$ where $d_2(n) = 3n^2-3n+1$.
$d_2(n) = d_2(n-1) + d_1(n)$ where $d_1(n) = 6n-6$.
$d_1(n) = d_1(n-1) + d_0(n)$ where $d_0(n) = 6$.

The main idea is that $d_2(n) = p(n)-p(n-1)$ is a polynomial of a lower degree. We can then compute $d_i(n)$ iteratively from $d_i(n-1)$ by adding $d_{i+1}(n)$. Thus $p(n)$ can be computed using the $3$ additions above with the initial conditions $d_0(0) = 6, d_1(0) = -6, d_2(0) = 1, p(0) = 0$.

Given the initial conditions at $n=0$, for a $m$-degree polynomial, this algorithm will require $mk$ additions to compute $p(n)$ for $n = 1, ..., k$.

For the On-Line Encyclopedia of Integer Sequences (OEIS), there are several sequences where such polynomial evaluations are needed.  I have written several Python programs to compute these sequences using Nuttall's algorithm.  See for example: A161710A161715A235810A138091A075904, A243298, A007487.

Consider the interpolating polynomial for the divisors of 24 (OEIS A161710):
$\frac{-6n^7 + 154n^6 - 1533n^5 + 7525n^4 - 18879n^3 + 22561n^2 - 7302n+2520}{2520}$

Since the polynomial evaluates to integers for all integers $n$, given the initial conditions, Nuttall's algorithm evaluates this polynomial at integers $n=1,2,...$ using only integer arithmetic, without the need for multiplication or division:

A161710_list, m = [1], [-12, 80, -223, 333, -281, 127, -23, 1]
for _ in range(1, 10**3):
    for i in range(7):
        m[i+1] += m[i]
    A161710_list.append(m[-1]) 

In Ref. [2], it was shown that the initial conditions take some computation and Nuttall's algorithm is more efficient than the direct approach in terms of the number of multiplication operations needed if the number of gridpoints is larger than $o(m^2)$ where $m$ is the degree of the polynomial. However, for the case above the savings in the number of multiplications is more substantial for 2 reasons.  First, the initial conditions can be precomputed and hardcoded as in the code above.
Second, in the above example only integer arithmetic is used and the numbers have increasing number of digits as $n\rightarrow\infty$. The initial conditions have smaller number of bits than the numbers that are multiplied in the direct approach to obtain $p(n)$ when $n$ is large and the number of multiplication operations needed on a fixed wordlength computer increases for larger numbers.

In Ref. [1], Nuttall also considered functions of the form $f(x) = e^{p(x)}$ where $p$ is a polynomial.  In this case $f(n)/f(n-1) = e^{p(n)-p(n-1)}$, so the successive differences in the polynomial becomes successive ratios in $f$ and the additions become multiplications.

Extensions of this algorithm to polynomials in multiple variables can be found in Ref. [3].

References:
[1] Albert H. Nuttall, "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1486-1487, October 1987.
[2] S. C. Datta Roy and Shailey Minocha, "A Note on "Efficient Evaluation of Polynomials and Exponentials of Polynomails for Equispaced Arguments",", IEEE Transactions on Signal Processing, vol. 39, no. 11, pp. 2554-2556, November 1991.
[3] N. K. Bose, "Efficient Evaluation of Multivariate Polynomials Over Multidimensional Grids," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, no. 10, pp. 1630-1631, November 1987.