Polynomial Shift | Continuous Point-Value Shift

Polynomial Shift

Polynomial shift is the simplest polynomial composition transform: given the coefficients of $f(x)=\sum _ {i=0}^nf_ix^i$ and a constant $c$, compute the coefficients of $f(x+c)$, namely $f(x)\mapsto f(x+c)$.

Divide-and-Conquer Method

Let

$$f(x)=f_0(x)+x^{\left\lfloor n/2\right\rfloor}f_1(x)$$

Then

$$f(x+c)=f_0(x+c)+(x+c)^{\left\lfloor n/2\right\rfloor}f_1(x+c)$$

The coefficient of $(x+c)^{\left\lfloor n/2\right\rfloor}$ is a binomial coefficient, so

$$T(n)=2T(n/2)+O(n\log n)=O(n\log^2 n)$$

where $O(n\log n)$ is the time complexity of polynomial multiplication.

Taylor Formula Method

Applying the Taylor formula to $f(x)$ at $c$ gives

$$f(x)=f(c)+\frac{f'(c)}{1!}(x-c)+\frac{f''(c)}{2!}(x-c)^2+\cdots +\frac{f^{(n)}(c)}{n!}(x-c)^n$$

Therefore

$$f(x+c)=f(c)+\frac{f'(c)}{1!}x+\frac{f''(c)}{2!}x^2+\cdots +\frac{f^{(n)}(c)}{n!}x^n$$

For $t\geq 0$, we have

$$\begin{aligned} t!\lbrack x^t\rbrack f(x+c)&=f^{(t)}(c)\\ &=\sum _ {i=t}^nf_ii!\frac{c^{i-t}}{(i-t)!}\\ &=\sum _ {i=0}^{n-t}f _ {i+t}(i+t)!\frac{c^i}{i!} \end{aligned}$$

Let

$$\begin{aligned} A_0(x)&=\sum _ {i=0}^nf _ {n-i}(n-i)!x^i\\ B_0(x)&=\sum _ {i=0}^n\frac{c^i}{i!}x^i \end{aligned}$$

Then

$$\begin{aligned} \lbrack x^{n-t}\rbrack (A_0(x)B_0(x))&=\sum _ {i=0}^{n-t} (\lbrack x^{n-t-i}\rbrack A_0(x))(\lbrack x^i\rbrack B_0(x))\\ &=\sum _ {i=0}^{n-t}f _ {i+t}(i+t)!\frac{c^i}{i!}\\ &=t!\lbrack x^t\rbrack f(x+c) \end{aligned}$$

Binomial Theorem Method

Using the binomial theorem $\displaystyle (a+b)^n=\sum _ {i=0}^n\binom{n}{i}a^ib^{n-i}$, we obtain

$$\begin{aligned} f(x+c)&=\sum _ {i=0}^nf_i(x+c)^i\\ &=\sum _ {i=0}^nf_i\left(\sum _ {j=0}^i\binom{i}{j}x^jc^{i-j}\right)\\ &=\sum _ {i=0}^nf_ii!\left(\sum _ {j=0}^i\frac{x^j}{j!}\frac{c^{i-j}}{(i-j)!}\right)\\ &=\sum _ {i=0}^n\frac{x^i}{i!}\left(\sum _ {j=i}^{n}f_jj!\frac{c^{j-i}}{(j-i)!}\right) \end{aligned}$$

This gives the same result as the method above.

Continuous Point-Value Shift

Example: LOJ 166 — Lagrange interpolation 2

Given a polynomial $f$ of degree at most $n$, together with its consecutive point values $f(0),f(1),\dots ,f(n)$, compute $f(c),f(c+1),\dots ,f(c+n)$ modulo $998244353$, where $1\leq n\leq 10^5,n < m\leq 10^8$.

Lagrange Interpolation Formula Method

Consider Lagrange interpolation formula.

$$\begin{aligned} f(x)&=\sum _ {0\leq i\leq n}f(i)\prod _ {0\leq j\leq n\,\land \,j\neq i}\frac{x-j}{i-j}\\ &=\sum _ {0\leq i\leq n}f(i)\frac{x!}{(x-n-1)!(x-i)}\frac{(-1)^{n-i}}{i!(n-i)!}\\ &=\frac{x!}{(x-n-1)!}\sum _ {0\leq i\leq n}\frac{f(i)}{(x-i)}\frac{(-1)^{n-i}}{i!(n-i)!} \end{aligned}$$

Although the formula above has a convolution form, the denominator $x-i\neq 0$ is not always invertible. We therefore consider only the case $c > n$ below; other cases (for example, when denominators of $B_0(x)$ coefficients vanish modulo a prime) can be handled by case analysis. Let

$$\begin{aligned} A_0(x)&=\sum _ {0\leq i\leq n}\frac{f(i)(-1)^{n-i}}{i!(n-i)!}x^i\\ B_0(x)&=\sum _ {i\geq 0}\frac{1}{c-n+i}x^i \end{aligned}$$

Then, for $t\geq 0$,

$$\begin{aligned} \lbrack x^{n+t}\rbrack (A_0(x)B_0(x))&=\sum _ {i=0}^{n+t}(\lbrack x^i\rbrack A_0(x))(\lbrack x^{n+t-i}\rbrack B_0(x))\\ &=\sum _ {i=0}^{n}\frac{f(i)(-1)^{n-i}}{i!(n-i)!}\frac{1}{c+t-i}\\ &=\frac{(c+t-n-1)!}{(c+t)!}f(c+t) \end{aligned}$$

In the implementation, truncate the required part of $B_0(x)$ to obtain more point values; cyclic convolution can also be used.

With a slight modification of the problem, suppose the point values of $d$ at $f(d),f(d+k),\dots ,f(d+nk)$ are given. We can compute $f(c+d),f(c+d+k),\dots ,f(c+d+nk)$, viewing them as the point values $g(x)=f(d+kx)$ of the shifted polynomial $g(0),g(1),\dots ,g(n)$ at $g(c/k),g(c/k+1),\dots ,g(c/k+n)$.

The Lagrange interpolation formula also gives a linear-time method for computing one point value by maintaining several prefix and suffix products.

Applications

Unsigned Stirling Numbers of the First Kind in One Row

Example: P5408 — unsigned Stirling numbers of the first kind

Compute $\displaystyle {n\brack 0},{n\brack 1},\dots ,{n\brack n}$ modulo the prime $167772161$, where $1\leq n< 262144$.

Consider

$$x^{\overline{n}}=\sum _ {i=0}^n{n\brack i}x^i,\quad n\geq 0$$

where $x^{\overline{n}}=x\cdot (x+1)\cdots (x+n-1)$ is the rising factorial, and let $f_n(x)=x^{\overline{n}}$. Then

$$f_{2n}(x)=x^{\overline{n}}\cdot (x+n)^{\overline{n}}=f_n(x)f_n(x+n)$$

Polynomial shift computes $O(n\log n)$ in $f_n(x+n)$, reducing the problem to half its original size: compute the coefficients of $f_n(x)$. Thus

$$T(n)=T(n/2)+O(n\log n)=O(n\log n)$$

Factorial Modulo a Prime

Example: P5282 [Template] — Fast factorial algorithm

Compute $n!\bmod p$, where $p$ is prime and $1\leq n< p\leq 2^{31}-1$.

Let $v=\lfloor\sqrt{n}\rfloor$ and $g(x)=\prod _ {i=1}^v(x+i)$. Then

$$n!\equiv \left(\prod _ {i=0}^{v-1}g(iv)\right)\cdot \prod _ {i=v^2+1}^n i\pmod{p}$$

The product $\prod _ {i=v^2+1}^n i$ can be computed in $O(\sqrt{n})$ time. We want to compute the first product in the formula efficiently.

Polynomial Multipoint Evaluation

The coefficients of $g(x)$ can be computed by the polynomial-shift algorithm in $O(n\log n)$, but evaluating it at the multiple points $g(0),g(v),g(2v),\dots ,g(v^2-v)$ takes $O(\sqrt{n}\log ^2n)$.

Continuous Point-Value Shift

Let $g_d(x)=\prod _ {i=1}^d(x+i)$. We can uniquely determine this degree-$d$ polynomial from its $d+1$ point values $g_d(0),g_d(v),\dots ,g_d(dv)$. Also,

$$g _ {2d}(x)=g_d(x)g_d(x+d)$$

Thus only $2d+1$ point values are needed to uniquely determine $g _ {2d}(x)$. Applying continuous point-value shift to compute $g_d((d+1)v),g_d((d+2)v),\dots ,g_d(2dv)$ (that is, shift the point values $h(x)=g_d(vx)$, $h(0),h(1),\dots ,h(d)$ to $h(d+1),h(d+2),\dots ,h(2d)$) and $g_d(d),g_d(v+d),\dots ,g_d(2dv+d)$ (that is, shift the point values $h(x)=g_d(vx)$, $h(0),h(1),\dots ,h(d)$ to $h(d/v),h(d/v+1),h(d/v+2),\dots ,h(d/v+2d)$), then multiplying the corresponding point values gives $g _ {2d}(0),g _ {2d}(v),\dots ,g _ {2d}(2dv)$.

To compute $g_d(0),g_d(v),\dots ,g_d(dv)$ from $g _ {d+1}(0),g _ {d+1}(v),\dots ,g _ {d+1}(dv),g _ {d+1}((d+1)v)$, consider

$$g _ {d+1}(x)=(x+d+1)\cdot g_d(x)$$

The one additional point value can be obtained by a linear-time algorithm. Therefore, starting with $g_1(0)=1,g_1(v)=v+1$, we maintain these point values by repeated doubling with continuous point-value shift, obtaining

$$T(n)=T(n/2)+O(n\log n)=O(n\log n)$$

We need only about $\sqrt{n}$ point values, so the time complexity is $O(\sqrt{n}\log n)$.

Prefix Sums of Binomial Coefficients Modulo a Prime

Example: LOJ 6386 — prefix sums of binomial coefficients

Compute $\displaystyle \sum _ {i=0}^m\binom{n}{i}\bmod 998244353$, where $0\leq m\leq n\leq 9\times 10^8$.

Describe the recurrence $n!=n\cdot (n-1)!$ using matrices. We have

$$\begin{bmatrix} n! \end{bmatrix} =\left( \prod _ {i=0}^{n-1} \begin{bmatrix}i+1\end{bmatrix} \right) \begin{bmatrix} 1 \end{bmatrix}$$

Similarly, the recurrence for the prefix sum of binomial coefficients can be written as

$$\begin{bmatrix} \binom{n}{m+1}\\ \sum _ {i=0}^m\binom{n}{i} \end{bmatrix}= \begin{bmatrix} (n-m)/(m+1)&0\\ 1&1 \end{bmatrix} \begin{bmatrix} \binom{n}{m}\\ \sum _ {i=0}^{m-1}\binom{n}{i} \end{bmatrix}$$

Taking the order of matrix multiplication into account,

$$\begin{aligned} \begin{bmatrix} \binom{n}{m+1}\\ \sum _ {i=0}^m\binom{n}{i} \end{bmatrix} &=\left( \prod _ {i=0}^{m} \begin{bmatrix} (n-i)/(i+1)&0\\1&1 \end{bmatrix} \right) \begin{bmatrix} 1\\0 \end{bmatrix}\\ &= \frac{1}{(m+1)!} \left( \prod _ {i=0}^{m} \begin{bmatrix} n-i&0\\i+1&i+1 \end{bmatrix} \right) \begin{bmatrix} 1\\0 \end{bmatrix} \end{aligned}$$

Let $v=\lfloor\sqrt{m}\rfloor$ and consider maintaining the matrix

$$\begin{aligned} M _ d(x)&= \prod _ {i=1}^d \begin{bmatrix} -x+n+1-i&0\\ x+i&x+i \end{bmatrix}\\ &= \begin{bmatrix} f_d(x)&0\\ g_d(x)&h_d(x) \end{bmatrix} \end{aligned}$$

The point values $M _ d(0),M _ d(v),\dots ,M_d(dv)$, namely $f_d(0),f_d(v),\dots ,f_d(dv)$, $h_d(0),\dots ,h_d(dv)$, and $g_d(0),\dots ,g_d(dv)$, together with

$$\begin{aligned} M _ {2d}(x)&= \prod _ {i=1}^{2d} \begin{bmatrix} -x+n+1-i&0\\ x+i&x+i \end{bmatrix}\\ &= \left( \prod _ {i=1}^d \begin{bmatrix} -x-d+n+1-i&0\\ x+d+i&x+d+i \end{bmatrix} \right) \left( \prod _ {i=1}^d \begin{bmatrix} -x+n+1-i&0\\ x+i&x+i \end{bmatrix} \right) \\ &= \begin{bmatrix} f_d(x+d)&0\\ g_d(x+d)&h_d(x+d) \end{bmatrix} \begin{bmatrix} f_d(x)&0\\ g_d(x)&h_d(x) \end{bmatrix} \\ &= \begin{bmatrix} f_d(x+d)f_d(x)&0\\ g_d(x+d)f_d(x)+h_d(x+d)g_d(x)&h_d(x+d)h_d(x) \end{bmatrix} \end{aligned}$$

The lower-right entry of the matrix is exactly the quantity maintained in the factorial algorithm, so

$$\begin{aligned} \prod _ {i=0}^{m} \begin{bmatrix} n-i&0\\i+1&i+1 \end{bmatrix}= \left( \prod _ {i=(k+1)v}^m \begin{bmatrix} n-i&0\\ i+1&i+1 \end{bmatrix} \right) \begin{bmatrix} f_v(kv)&0\\ g_v(kv)&h_v(kv) \end{bmatrix} \cdots \begin{bmatrix} f_v(0)&0\\ g_v(0)&h_v(0) \end{bmatrix} \end{aligned}$$

The computation can be completed in $O(\sqrt m\log m)$ time.

Harmonic Numbers Modulo a Prime

Example: P5702 — summing a harmonic series

Compute $\sum _ {i=1}^ni^{-1}\bmod p$, where $p$ is prime and $1\leq n< p< 2^{30}$.

Let $H_n=\sum _ {k=1}^nk^{-1}$. One step of the recurrence is

$$\begin{bmatrix} (n+1)!\\(n+1)!H _ {n+1} \end{bmatrix}= \begin{bmatrix} n+1&0\\1&n+1 \end{bmatrix} \begin{bmatrix} n!\\n!H_n \end{bmatrix}$$

Thus

$$\begin{bmatrix} {n+1\brack 1}\\{n+1\brack 2} \end{bmatrix}= \begin{bmatrix} n!\\n!H_n \end{bmatrix}= \left( \prod _ {i=0}^{n-1} \begin{bmatrix} i+1&0\\1&i+1 \end{bmatrix} \right) \begin{bmatrix} 1\\0 \end{bmatrix}$$

Here $\displaystyle {n+1\brack 1}$ and $\displaystyle {n+1\brack 2}$ are unsigned Stirling numbers of the first kind. The matrix method for maintaining point values is the same as above.

Polynomial Recurrences

For the more general case, analogous to the fast-factorial examples above, what kind of algorithm can we obtain?

Example: P6115 [Template] — Polynomial Recurrences

Given a sequence $a$ satisfying $\forall n\ge m,\sum_{k=0}^ma_{n-k}P_k(n)=0$, where $P_k$ is a polynomial of degree at most $d$. Given the coefficients of all $P_k$ and the initial values $a_0,a_1,\dots,a_{m-1}$, find $a_n$. Compute modulo $998244353$; $n\le6\times10^8$, $1\le m,d\le7$, $7s$, with a time limit of {11}.

To describe the matrix constructions in the examples above more systematically, we introduce the concept of lambda matrix ($\lambda$).

For polynomial recurrences, note that the point values maintained in the fast-factorial algorithm are not actually $n!$, but $\prod_{i=0}^{T-1}(aT+i)$: the ratio between a pair of point values.

Because the order of a polynomial recurrence, $m$, is no longer just $1$, we cannot directly maintain the ratio between two numbers. Instead, we maintain a linear transformation between a pair of $m$-dimensional vectors, namely an $m\times m$ matrix; each matrix entry is the point value of a polynomial.

For a general polynomial-recurrence problem asking for a distant term, we can construct

$$-{\frac{1}{P_0(n)}}\begin{bmatrix}P_1(n)&P_2(n)&P_3(n)&\cdots&P_{m-1}(n)&P_m(n)\\-P_0(n)\\&-P_0(n)\\&&-P_0(n)\\&&&\ddots\\&&&&-P_0(n)\\\end{bmatrix} \begin{bmatrix}a_{n-1}\\a_{n-2}\\a_{n-3}\\\vdots\\a_{n-m+1}\\a_{n-m}\end{bmatrix} =\begin{bmatrix}a_n\\a_{n-1}\\a_{n-2}\\\vdots\\a_{n-m+2}\\a_{n-m+1}\end{bmatrix}$$

Let

$$B(\lambda)=\begin{bmatrix} P_1(\lambda)&P_2(\lambda)&P_3(\lambda)&\cdots&P_{m-1}(\lambda)&P_m(\lambda)\\ -P_0(\lambda)\\ &-P_0(\lambda)\\ &&-P_0(\lambda)\\ &&&\ddots\\ &&&&-P_0(\lambda)\\ \end{bmatrix}$$

Ignoring the leading factor $-\frac1{P_0(n)}$ for now, we need to maintain a quantity of the form $\prod_{i=0}^{T-1}B(aT+m+i)$, where multiplication is performed from right to left.

It is easy to see that $B_T(\lambda)=\prod_{i=0}^{T-1}B(\lambda+i)$ is a $\lambda$ matrix whose entries have degree at most $dT$; it can be maintained using only $dT+1$ values.

We therefore maintain the point values of the $B_T(m)$ matrices $B_T(m+T)$, $B_T(m+2T)$, $\dots$, $B_T(m+(dT-1)T)$, $B_T(m+dT^2)$, and $\lambda$, and use polynomial point-value shifting and doubling in the same way as in the fast-factorial algorithm.

Specifically, to lift $t=\log_2T$ to $1$, we proceed as follows:

1. In $O(m^2dT\log(dT))$ time, obtain $B_T(p+dT^2)$, $B_T(p+(dT+1)T)$, $B_T(p+(dT+2)T)$, $\cdots$, $B_T(p+(2dT-1)T)$, and $B_T(p+(2dT)dT)$.
2. In $O(m^2dT\log(dT))$ time, obtain $B_T(p+2dT^2)$, $B_T(p+(2dT+1)T)$, $B_T(p+(2dT+2)T)$, $\cdots$, $B_T(p+(3dT-1)T)$, and $B_T(p+(3dT)dT)$.
3. In $O(m^2dT\log(dT))$ time, obtain $B_T(p+3dT^2)$, $B_T(p+(3dT+1)T)$, $B_T(p+(3dT+2)T)$, $\cdots$, $B_T(p+(4dT-1)T)$, and $B_T(p+(4dT)dT)$.
4. Compute $B_{2T}(v)=B_{T}(v+T)B_{T}(v)$.

Each round spends $O(m^2dT\log(dT))$ complexity on shifting. Each round performs only $\Theta(dT)$ matrix multiplication, whose complexity can be regarded as $O(m^3dT)$.

Finally, it is enough to reach $T\ge\sqrt{n/d}$.

The preceding $-\frac1{P_0(n)}$ factor can be handled similarly.

Thus the preprocessing complexity is $\Theta(\sqrt{nd}(m^3+m^2\log(nd)))$.

For a query, we need only $\Theta(n/T)$ vector-matrix multiplications and $O(T)$ brute-force transitions.

This part is clearly not the complexity bottleneck.

Therefore, the total complexity of the algorithm is $\Theta(\sqrt{nd}(m^3+m^2\log(nd)))$.

When coding, cyclic convolution can be used to reduce the constant factor of the NTT.

In applications, we often extract a distant coefficient of a known D-finite GF, so $m,d$ are all constants; this gives distant-coefficient evaluation in $\Theta(\sqrt n\log n)$.

References

• Alin Bostan, Pierrick Gaudry, and Eric Schost. Linear recurrences with polynomial coefficients and application to integer factorization and Cartier–Manin operator.
• Min_25's blog
• ZZQ's blog - factorial modulo a large prime