跳转至

多項式平移|連續點值平移

多項式平移

多項式平移是簡單情況的多項式複合變換，給出 \(f(x)=\sum _ {i=0}^nf_ix^i\) 的係數和一個常數 \(c\)，求 \(f(x+c)\) 的係數，即 \(f(x)\mapsto f(x+c)\)。

分治法

令

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

那麼

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

\((x+c)^{\left\lfloor n/2\right\rfloor}\) 的係數為二項式係數，那麼

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

其中 \(O(n\log n)\) 為多項式乘法的時間。

Taylor 公式法

對 \(f(x)\) 在 \(c\) 處應用 Taylor 公式，有

\[ 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 \]

那麼

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

觀察到對於 \(t\geq 0\) 有

\[ \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} \]

令

\[ \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} \]

那麼

\[ \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} \]

二項式定理法

考慮二項式定理 \(\displaystyle (a+b)^n=\sum _ {i=0}^n\binom{n}{i}a^ib^{n-i}\) 那麼

\[ \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} \]

得到的結果與上述方法相同。

連續點值平移

例題 LOJ 166 拉格朗日插值 2

給出度數小於等於 \(n\) 的多項式 \(f\) 的連續點值 \(f(0),f(1),\dots ,f(n)\)，在模 \(998244353\) 意義下計算 \(f(c),f(c+1),\dots ,f(c+n)\)，其中 \(1\leq n\leq 10^5,n < m\leq 10^8\)。

Lagrange 插值公式法

考慮 Lagrange 插值公式

\[ \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} \]

上式雖然是卷積形式但不能保證分母上 \(x-i\neq 0\)，所以下面僅考慮 \(c > n\) 的情況，其他情況（如係數在模素數意義下時須避免 \(B_0(x)\) 係數的分母出現零）可以分類討論解決，令

\[ \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} \]

那麼對於 \(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} \]

實現中取 \(B_0(x)\) 需要的部分截斷可求出更多點值，且可利用循環卷積。

對問題稍加修改，假設對於某個 \(d\) 給出的點值為 \(f(d),f(d+k),\dots ,f(d+nk)\)，我們可以計算 \(f(c+d),f(c+d+k),\dots ,f(c+d+nk)\)，視作平移 \(g(x)=f(d+kx)\) 的點值 \(g(0),g(1),\dots ,g(n)\) 為 \(g(c/k),g(c/k+1),\dots ,g(c/k+n)\)。

Lagrange 插值公式也給出了通過維護一些前後綴積的線性計算單個點值的方法。

應用

同一行第一類無符號 Stirling 數

例題 P5408 第一類斯特林數·行

在模素數 \(167772161\) 意義下求 \(\displaystyle {n\brack 0},{n\brack 1},\dots ,{n\brack n}\)，其中 \(1\leq n< 262144\)。

考慮

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

其中 \(x^{\overline{n}}=x\cdot (x+1)\cdots (x+n-1)\) 為上升階乘冪，令 \(f_n(x)=x^{\overline{n}}\) 那麼

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

通過多項式平移可在 \(O(n\log n)\) 求出 \(f_n(x+n)\)，問題被縮小為原先的一半即求出 \(f_n(x)\) 的係數，那麼

\[ T(n)=T(n/2)+O(n\log n)=O(n\log n) \]

模素數意義下階乘

例題 P5282【模板】快速階乘算法

求 \(n!\bmod p\)，其中 \(p\) 為素數且 \(1\leq n< p\leq 2^{31}-1\)。

令 \(v=\lfloor\sqrt{n}\rfloor\) 和 \(g(x)=\prod _ {i=1}^v(x+i)\) 那麼

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

其中 \(\prod _ {i=v^2+1}^n i\) 可在 \(O(\sqrt{n})\) 時間計算，我們希望可以快速計算上式的前半部分。

多項式多點求值

\(g(x)\) 係數的計算可用上述多項式平移算法在 \(O(n\log n)\) 時間得到，但多點求值計算 \(g(0),g(v),g(2v),\dots ,g(v^2-v)\) 需要 \(O(\sqrt{n}\log ^2n)\) 時間。

連續點值平移

令 \(g_d(x)=\prod _ {i=1}^d(x+i)\)，我們可以用 \(d+1\) 個點值 \(g_d(0),g_d(v),\dots ,g_d(dv)\) 唯一確定這個度數為 \(d\) 的多項式，又

\[ g _ {2d}(x)=g_d(x)g_d(x+d) \]

所以只需 \(2d+1\) 個點值可以唯一確定 \(g _ {2d}(x)\)，那麼使用連續點值平移計算 \(g_d((d+1)v),g_d((d+2)v),\dots ,g_d(2dv)\)（即平移 \(h(x)=g_d(vx)\) 的點值 \(h(0),h(1),\dots ,h(d)\) 為 \(h(d+1),h(d+2),\dots ,h(2d)\)）和 \(g_d(d),g_d(v+d),\dots ,g_d(2dv+d)\)（即平移 \(h(x)=g_d(vx)\) 的點值 \(h(0),h(1),\dots ,h(d)\) 為 \(h(d/v),h(d/v+1),h(d/v+2),\dots ,h(d/v+2d)\)）後將這兩者的對應點值相乘即得 \(g _ {2d}(0),g _ {2d}(v),\dots ,g _ {2d}(2dv)\)。

由 \(g_d(0),g_d(v),\dots ,g_d(dv)\) 計算 \(g _ {d+1}(0),g _ {d+1}(v),\dots ,g _ {d+1}(dv),g _ {d+1}((d+1)v)\) 考慮

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

額外增加的一個點值使用線性時間的算法即可。那麼在開始時維護 \(g_1(0)=1,g_1(v)=v+1\) 後使用連續點值平移來倍增地維護這些點值，有

\[ T(n)=T(n/2)+O(n\log n)=O(n\log n) \]

而我們只需要約 \(\sqrt{n}\) 個點值，所以時間複雜度為 \(O(\sqrt{n}\log n)\)。

模素數意義下二項式係數前綴和

例題 LOJ 6386 組合數前綴和

求 \(\displaystyle \sum _ {i=0}^m\binom{n}{i}\bmod 998244353\)，其中 \(0\leq m\leq n\leq 9\times 10^8\)。

考慮使用矩陣描述 \(n!=n\cdot (n-1)!\) 這一步遞推，我們有

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

類似的可以將二項式係數前綴和的遞推描述為

\[ \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} \]

注意矩陣乘法的順序，那麼

\[ \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} \]

令 \(v=\lfloor\sqrt{m}\rfloor\)，考慮維護矩陣

\[ \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} \]

的點值 \(M _ d(0),M _ d(v),\dots ,M_d(dv)\) 即 \(f_d(0),f_d(v),\dots ,f_d(dv)\)、\(h_d(0),\dots ,h_d(dv)\) 和 \(g_d(0),\dots ,g_d(dv)\)，又

\[ \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} \]

且矩陣右下角元素恰為我們在階乘算法中所維護的，那麼

\[ \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} \]

可在 \(O(\sqrt m\log m)\) 時間完成計算。

模素數意義下調和數

例題 P5702 調和級數求和

求 \(\sum _ {i=1}^ni^{-1}\bmod p\)，其中 \(p\) 為素數且 \(1\leq n< p< 2^{30}\)。

記 \(H_n=\sum _ {k=1}^nk^{-1}\)，一步遞推為

\[ \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} \]

那麼

\[ \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} \]

在這裏 \(\displaystyle {n+1\brack 1}\) 和 \(\displaystyle {n+1\brack 2}\) 為第一類無符號 Stirling 數。維護點值矩陣的方法同上。

整式遞推

對於更一般的情況，類似於上述快速階乘算法的案例，我們期望得到一個怎麼樣的算法？

例題 P6115【模板】整式遞推

現有數列 \(a\) 滿足 \(\forall n\ge m,\sum_{k=0}^ma_{n-k}P_k(n)=0\)，其中 \(P_k\) 為不超過 \(d\) 次的多項式。
給定所有 \(P_k\) 的係數，和 \(a_0,a_1,\dots,a_{m-1}\)，求 \(a_n\)。 對 \(998244353\) 取模。\(n\le6\times10^8\)，\(1\le m,d\le7\)，時限 \(7s\)。

為了更系統地描述上述幾道例題中構造矩陣的過程，我們引入 \(\lambda\) 矩陣 的概念。

為了實現整式遞推，我們應當注意到快速階乘算法過程中，我們維護的點值其實並不是 \(n!\)，而是 \(\prod_{i=0}^{T-1}(aT+i)\)，即 一對點值之間的倍數關係。

由於整式遞推階數 \(m\) 不止是 \(1\) 了，我們就 不能直接維護一對數之間的倍數關係了；而是維護出 一對 \(m\) 維向量之間的線性變換，即 \(m\times m\) 的一個矩陣，矩陣的每一項對應於某個多項式的一個點值。

容易發現，對於一般的整式遞推遠處係數求值問題，我們可以構造

\[ -{\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} \]

設

\[ 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} \]

我們先撇開前面的 \(-\frac1{P_0(n)}\) 因子不論，我們現在要維護 \(\prod_{i=0}^{T-1}B(aT+m+i)\) 這種形式的量，其中乘法自右往左。

容易發現 \(B_T(\lambda)=\prod_{i=0}^{T-1}B(\lambda+i)\) 是一個各項次數不高於 \(dT\) 的 \(\lambda\) 矩陣，只用 \(dT+1\) 個值即可維護。

於是我們維護出 \(B_T(m)\)，\(B_T(m+T)\)，\(B_T(m+2T)\)，\(\dots\)，\(B_T(m+(dT-1)T)\)，\(B_T(m+dT^2)\) 這幾個 \(\lambda\) 矩陣的 點值，然後用類似於快速階乘算法的方式暴力進行多項式點值平移和倍增就好了。

具體地，為了讓 \(t=\log_2T\) 抬高 \(1\)，我們這麼幹：

  1. 在 \(O(m^2dT\log(dT))\) 時間內獲取 \(B_T(p+dT^2)\)，\(B_T(p+(dT+1)T)\)，\(B_T(p+(dT+2)T)\)，\(\cdots\)，\(B_T(p+(2dT-1)T)\)，\(B_T(p+(2dT)dT)\)。
  2. 在 \(O(m^2dT\log(dT))\) 時間內獲取 \(B_T(p+2dT^2)\)，\(B_T(p+(2dT+1)T)\)，\(B_T(p+(2dT+2)T)\)，\(\cdots\)，\(B_T(p+(3dT-1)T)\)，\(B_T(p+(3dT)dT)\)。
  3. 在 \(O(m^2dT\log(dT))\) 時間內獲取 \(B_T(p+3dT^2)\)，\(B_T(p+(3dT+1)T)\)，\(B_T(p+(3dT+2)T)\)，\(\cdots\)，\(B_T(p+(4dT-1)T)\)，\(B_T(p+(4dT)dT)\)。
  4. 計算 \(B_{2T}(v)=B_{T}(v+T)B_{T}(v)\)。

我們每輪花費 \(O(m^2dT\log(dT))\) 的複雜度進行平移；同時，我們每輪只用做 \(\Theta(dT)\) 次矩陣乘法，複雜度可以認為是 \(O(m^3dT)\)。

最後，我們只用做到 \(T\ge\sqrt{n/d}\) 即可。

之前的 \(-\frac1{P_0(n)}\) 因子可以用類似的方法解決。

這樣，我們預處理的複雜度即為 \(\Theta(\sqrt{nd}(m^3+m^2\log(nd)))\)。

考慮查詢，我們只用 \(\Theta(n/T)\) 次向量與矩陣的乘法，以及 \(O(T)\) 次暴力轉移。

容易發現這部分計算並不是複雜度瓶頸。

因此，該算法的總複雜度為 \(\Theta(\sqrt{nd}(m^3+m^2\log(nd)))\)。

編碼時，我們可以使用循環卷積的技巧來減小 NTT 的常數。

在實際應用時，我們往往是對一個已知的微分有限的 GF 提取其遠處係數，從而 \(m,d\) 均為常數，也即做到了 \(\Theta(\sqrt n\log n)\) 的遠處係數求值。

參考文獻

  • Alin Bostan, Pierrick Gaudry, and Eric Schost. Linear recurrences with polynomial coefficients and application to integer factorization and Cartier–Manin operator.
  • Min_25 的博客
  • ZZQ 的博客 - 階乘模大質數
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用