[Tutorial] Polynomial Sample Shift

Today, we will be solving the Polynomial Sample Shift problem with in $$$O(n \log n)$$$ time using one convolution.

Problem: We are told that $$$f(x)$$$ is a polynomial with degree $$$d$$$ and we are given the samples $$$f(0), f(1), \dots, f(d)$$$, as well as a number $$$m$$$, and are asked to compute $$$f(m), f(m + 1), \dots, f(m + d)$$$.

Since the polynomial is degree $$$d$$$ and we are given $$$d + 1$$$ samples, $$$f(x)$$$ is uniquely defined. Using polynomial interpolation and then multi-point evaluation we can compute this in $$$O(d \log^2 d)$$$ using FFTs, but for this specific problem we can do better.

We start off with the Lagrange Interpolation Formula:

$$$\displaystyle f(x) = \sum_{i=0}^d f(i) \prod_{j = 0, j \neq i}^d \frac{x - j}{i - j}$$$

The goal is to plug in $$$x = m + k$$$ and get a formula in terms of $$$f(0), f(1), \dots, f(d)$$$.

$$$ \begin{align*} f(m + k) &= \sum_{i=0}^d f(i) \prod_{j = 0, j \neq i}^d \frac{m+k-j}{i-j} \\ &= \sum_{i=0}^d \frac{f(i)}{\prod_{j=0,j\neq i}^d (i-j)} \prod_{j=0,j\neq i}^d (m + k - j) \end{align*} $$$

We can precompute the product parts of this formula. First, the product in the denominator, letting

$$$\displaystyle p[i] = \prod_{j=0,j\neq i}^d (i - j)$$$

We can precompute this for $$$i = 0, 1, \dots, d$$$ by directly computing $$$p[0] = \prod_{j=0,j\neq i}^d (0-j)$$$, then using the formula $$$p[i] = \frac i{i-d-1}\ {p[i-1]}$$$ for the rest.

Next, we can remove the $$$j\neq i$$$ from the product on the right allowing that term to be in the product, but dividing by it, i.e.

$$$\displaystyle \prod_{j=0,j\neq i}^d (m + k - j) = \frac1{m + k - i}\ \prod_{j=0}^d (m + k - j) $$$

We also want to precompute this product. Let

$$$ q[k] = \prod_{j=0}^d (m + k - j) $$$

We can similarly precompute this for $$$k = 0, 1, \dots, d$$$ by directly computing $$$q[0] = \prod_{j=0}^d (m - j)$$$ then using the formula $$$q[k] = \frac{m+k}{m+k-d-1}\ q[k -1]$$$ for the rest.

So we have

$$$ \begin{align*} f(m + k) &= \sum_{i=0}^d \frac{f(i)}{p[i]}\ \frac1{m+k-i}\ q[k] \\ &= q[k] \sum_{i=0}^d \frac{f(i)}{p[i]}\ \frac1{m+k-i} \end{align*} $$$

Now each term sum is the product of two parts: one in terms of $$$i$$$, and one in term of $$$-i$$$, so we can express this as a convolution.

Let

$$$ \displaystyle a(x) = \sum_{i=0}^d \frac{f(i)}{p[i]} x^i \ \ \text{ and } \ \ b(x) = \sum_{i=0}^{2d} \frac1{m+i-d} x^i $$$

be polynomials with degree $$$d$$$ and $$$2d$$$ respectively.

Let $$$[x^i]p[x]$$$ denote the $$$i$$$-th coefficient of polynomial $$$p[x]$$$. The $$$c$$$-th term of the product is

$$$ \begin{align*} &\sum_{i=0}^{\min(c, d)} [x^i]a(x)\ [x^{c-i}]b(x)\ x^c\\ =&\sum_{i=0}^{\min(c, d)} \frac{f(i)}{p[i]} \frac1{m+(c-i)-d}\ x^c \end{align*} $$$

In particular, the $$$(k + d)$$$-th coefficient of the product is

$$$ \begin{align*} &\sum_{i=0}^{\min(k + d, d)} \frac{f(i)}{p[i]} \frac1{m+((k + d) -i)-d}\\ &=\sum_{i=0}^d \frac{f(i)}{p[i]} \frac1{m+ k - i} \end{align*} $$$

which matches exactly. Therefore, with $$$O(d)$$$ precomputation of $$$p[i]$$$ and $$$q[k]$$$ plus one FFT for the polynomial product, we have computed

$$$\displaystyle f(m + k) = q[k] \sum_{i=0}^d \frac{f(i)}{p[i]}\ \frac1{m+k-i} $$$

for all $$$k = 0, 1, \dots, d$$$, which is what we wanted.

Summary

The algorithm is as follows:

Input: An array $$$[f(0), f(1), \dots, f(d)]$$$, and a number $$$m$$$.

Output: An array $$$[f(m), f(m + 1), \dots, f(m + d)]$$$.

Compute $$$p[0] = \prod_{j = 0, j \neq i}^d (-j)$$$ and $$$p[i]$$$ for all $$$i = 1, \dots, d$$$ with $$$p[i] = \frac i{i - d - 1}p[i - 1]$$$
Compute $$$q[0] = \prod_{j=0}^d (m - j)$$$ and $$$q[k]$$$ for all $$$k = 1, \dots, d$$$ with $$$q[k] = \frac{m+k}{m+k-d-1}q[k-1]$$$
Compute $$$a[i] = \frac{f(i)}{p[i]}$$$ for $$$i = 0, 1, \dots, d$$$ and $$$b[i] = \frac1{m+i-d}$$$ for $$$i = 0, 1, \dots, 2d$$$.
Compute $$$c[i]$$$, the convolution of $$$a[i]$$$ and $$$b[i]$$$.
Pointwise multiply the coefficients of $$$q[k]$$$ and the $$$k$$$-th to $$$k + d$$$-th term of $$$c[i]$$$, and output.

Applications

Evenly Spaced Samples

It's worth noting that the samples we're given don't have to be sequential, just evenly spaced apart. Suppose we are told that $$$g(x)$$$ is a polynomial with degree $$$d$$$ and we are given the samples $$$g(a), g(a + b), g(a + 2b), \dots, g(a + db)$$$, and a number $$$m$$$ and are asked to compute $$$g(a + m), g(a + b + m), g(a + 2b + m), \dots, g(a + db + m)$$$.

We define $$$f(x) = g(a + bx)$$$ so that the samples given to us can be interpreted as $$$f(0), f(1), f(2), \dots, f(d)$$$. Performing a shift by $$$\frac mb$$$, we get $$$f(\frac mb), f(1 + \frac mb), f(2 + \frac{m}b), \dots, f(d + \frac mb)$$$ which are $$$g(a + m), g(a + b + m), g(a + 2b + m), \dots, g(a + db + m)$$$, as desired.

Fast Factorial Modulo p

Calculating factorials can be done with multi-point evaluation for a polynomial whose roots and sampling points are evenly spaced. Using this technique, we can calculate $$$n! \bmod p$$$ in $$$O(\sqrt n \log n)$$$ time. You can read about this technique in my blog post, or from here or here.

Problems

Shift of Sampling Points of Polynomial

Additional Notes

Technically the polynomial we are given is promised to have degree at most $$$d$$$, as we do not care whether e.g. the $$$x^d$$$, $$$x^{d-1}, \dots$$$ terms are 0.
If we attempt to shift by a value such that the output samples overlaps the given samples, which happens if $$$m$$$ is one of $$$0, 1, \dots, d$$$, then there is a division-by-zero in the above formula. In such a case, we should choose $$$m = d + 1$$$ and discard the extra values.
Looking closer at $$$p[i] = \prod_{j = 0, j \neq i}^d (i - j)$$$, this product has $$$d$$$ terms, the $$$i$$$ positive ones span from $$$1$$$ to $$$i$$$ and the $$$d - i$$$ negative ones span from $$$-1$$$ to $$$-(d - i)$$$, so we can alternatively write $$$p[i] = i! (d - i)! (-1)^{d-i}$$$. Either way, it is $$$O(d)$$$ time to compute all $$$d + 1$$$ values.

Sources

The math and derivations are based on Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator by Alin Bostan, Pierrick Gaudry, and Eric Schost.