[拉格朗日反演][FFT][NTT][多項式大全]詳解

一、多項式的兩種表示法

1.係數表示法

咱們最經常使用的多項式表示法就是係數表示法，一個次數界爲\(n\)的多項式\(S(x)\)能夠用一個向量\(s=(s_0,s_1,s_2,\cdots,s_n-1)\)係數表示以下：\[S(x)=\sum_{k=0}^{n-1}s_kx^k\]

係數表示法很適合作加法，能夠在\(O(n)\)的時間複雜度內完成，表達式爲：\[S(x)=A(x)+B(x)=\sum_{k=0}^{n-1}(a_k+b_k)x^k\]

當中\[s_k=a_k+b_k\]

可是，係數表示法不適合作乘法，時間複雜度爲\(O(n^2)\)，表達式爲：\[S(x)=A(x)B(x)=\sum_{k=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j b_{k-j}\right)x^k\]

當中\[s_k=\sum_{j=0}^k a_j b_{k-j}\]

這就是卷積的通常形式，記\(s=a\otimes b\)，咱們要想辦法加速這個過程。

2.點值表示法

顧名思義，點值就是多項式在一個點處的值。多項式\(A(x)\)的點值表達是一個集合：\[\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_{n-1},y_{n-1})\}\]

使得對於\(k=0,1,2,\cdots,n-1\)有\(x_k\)兩兩不相同且\(y_k=A(x_k)\)。

\(n\)個點能夠肯定惟一一個\(n\)次多項式。

點值表達有不少優良的性質，加法和乘法均可以在\(O(n)\)的時間複雜度內完成。

現有\(A(x)\)的點值表達\[\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_{n-1},y_{n-1})\}\]和\(B(x)\)的點值表達\[\{(x_0,y'_0),(x_1,y'_1),(x_2,y'_2),\cdots,(x_{n-1},y'_{n-1})\}\]

則\(C(x)=A(x)+B(x)\)的點值表達爲：\[\{(x_0,y_0+y'_0),(x_1,y_1+y'_1),(x_2,y_2+y'_2),\cdots,(x_{n-1},y_{n-1}+y'_{n-1})\}\]

\(C(x)=A(x)B(x)\)的點值表達爲：\[\{(x_0,y_0 y'_0),(x_1,y_1 y'_1),(x_2,y_2 y'_2),\cdots,(x_{n-1},y_{n-1} y'_{n-1})\}\]

可見，點值表示能夠幫助咱們更快地進行卷積，但是如何在係數表示法和點值表示法之間相互轉化呢？

二、複數

當\(x\)爲實數時，沒法很好地對轉換方法進行優化。爲了優化計算\(x^n\)所浪費的時間，咱們須要\(x\)有循環的性質。但點值表示法須要\(n\)個兩兩不一樣的值，而在實數域中只有\(1\)和\(-1\)，所以，咱們須要複數的幫助。

1.複數、複平面的定義

咱們把形如\(a+bi\)的數稱爲複數\(z\)，其中\(a\)爲實部(Real)，記爲\(\Re z\)；\(b\)爲虛部(Imaginary)，記爲\(\Im z\)。

每一點都對應惟一複數的平面叫複平面，至關於一個把\(\Re z\)做爲橫座標，把\(\Im z\)做爲縱座標的笛卡爾座標系。如圖：

模長：複平面上原點到複數\(z\)的距離，記爲\(|z|\)。根據勾股定理，\(|z|=|a+bi|=\sqrt{a^2+b^2}\)

輻角：複平面上\(x\)軸與複數\(z\)所對應向量之間的夾角，在\((-\frac{\pi}{2},\frac{\pi}{2})\)之間的記爲輻角主值\(\arg z\)。

2.歐拉公式

大名鼎鼎的歐拉公式：\[e^{i t}=\cos t+i\sin t\]

根據三角函數在單位圓上的幾何意義，公式是容易理解的。

幾何意義：

當中\(\theta\)爲角度，\(t\)爲弧長。

證實：

對\(\exp(x)\)在\(x=0\)泰勒展開：\[\exp(x)=\sum_{k=0}^\infty\frac{x^k}{k!}\]

對\(\sin(x)\)在\(x=0\)泰勒展開：\[\sin(x)=\sum_{k=0}^\infty\frac{x^{4 k+1}}{(4 k+1)!}-\sum_{k=0}^\infty\frac{x^{4 k+3}}{(4 k+3)!}\]

對\(\cos(x)\)在\(x=0\)泰勒展開：\[\cos(x)=\sum_{k=0}^\infty\frac{x^{4 k}}{(4 k)!}-\sum_{k=0}^\infty\frac{x^{4 k+2}}{(4 k+2)!}\]

把\(\exp(x)\)中的\(x\)替換爲\(i x\)：\[\exp(i x)=\sum_{k=0}^\infty\frac{x^{4 k}}{(4 k)!}+i\sum_{k=0}^\infty\frac{x^{4 k+1}}{(4 k+1)!}-\sum_{k=0}^\infty\frac{x^{4 k+2}}{(4 k+2)!}-i\sum_{k=0}^\infty\frac{x^{4 k+3}}{(4 k+3)!}\]

則顯然有：\[e^{i x}=\cos x+i\sin x\]

證畢。

則根據歐拉公式，可將一個複數表示爲一個二元組\((a,\theta)\)，即模長和輻角（至關於複平面上極座標系的表示方法）。值爲:\(a(\cos\theta+i\sin\theta)\)

特殊狀況：歐拉恆等式\(e^{i\pi}+1=0\)

3.複數的運算

（1）複數加法

運算規則：實部、虛部分別相加\[(a+bi)+(c+di)=a+c+bi+di=(a+c)+(b+d)i\]

幾何意義：如圖

結果至關於兩個向量所構成的平行四邊形的對角線。若是把一個複數所對應的向量視爲一個移動的變換，那麼向量加法就是連續運用這兩個變換至關於的新變換。

（2）複數乘法

運算規則：展開\[(a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+bc)i\]

幾何意義：如圖

如圖，\(\arg a+\arg b=\arg a\times b\),\(|a|\times|b|=|a\times b|\)

總結就是：模長相乘，輻角相加。

所以，若是模長爲\(1\)，那麼它的\(n\)次方必定還在單位圓上。

證實：

根據歐拉公式，已知\(x=(a_1,\theta_1)=a_1(\cos\theta_1+i\sin\theta_1),y=(a_2,\theta_2)=a_2(\cos\theta_2+i\sin\theta_2)\)

則\[\begin{align*} x\times y&=a_1 a_2(\cos\theta_1+i\sin\theta_1)(\cos\theta_2+i\sin\theta_2)\\ &=a_1 a_2\left[(\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2)+i(\cos\theta_1\sin\theta_2+\sin\theta_1\cos\theta_2)\right]\\ &=a_1 a_2\left[\left(\frac{\cos(\theta_1+\theta_2)+\cos(\theta_1-\theta_2)}{2}+\frac{\cos(\theta_1+\theta_2)-\cos(\theta_1-\theta_2)}{2}\right)\right.\\ &\left.+i\left(\frac{\sin(\theta_1+\theta_2)-\sin(\theta_1-\theta_2)}{2}+\frac{\sin(\theta_1+\theta_2)+\sin(\theta_1-\theta_2)}{2}\right)\right]\tag{積化和差公式}\\ &=a_1 a_2\left[\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2)\right] \end{align*}\]

\(\therefore |x\times y|=|x|\times |y|,\arg(x\times y)=\arg x+\arg y\)

證畢。

4.單位複數根

（1）基本性質

單位複數根是方程\(\omega^n=1\)的解，第\(k\)個解記爲\(\omega_n^k\)(這裏的\(k\)事實上是乘方的含義)

\(n=16\)的解在複平面上的位置以下：

能夠看到，\(n\)個解把單位圓分紅了\(n\)等弧，交點即爲根。並且，\(\omega_n^k\)其實是\(\omega_n\)的\(n\)次方，模長仍爲\(1\),輻角翻倍！

爲何呢？

\(\because |x^n|=|x|^n,\arg x^n=n\arg x\)

\(\therefore |\omega|^n=|\omega^n|,\arg\omega^n=n\arg\omega\)

\(\therefore |\omega|^n=1(|\omega|\in\mathbb{R}^+),\arg\omega=\frac{360^\circ}{n}\)

\(\therefore |\omega|=1,\arg\omega=\frac{360^\circ}{n}\)

這就很明顯了。

因此，\(\omega_n^k\)事實上表示的是\(\omega_n\)的\(k\)次冪。爲何選擇單位複數根呢？由於它有循環的優良性質，即\(\omega_n^n=1\)。因爲其餘的均可以由\(\omega_n^1\)獲得，所以稱爲主\(n\)次單位根，又記爲\(\omega_n\)。

根據單位複數根的平分圓的意義和歐拉公式，\(\omega_n^k=e^\frac{2\pi i k}{n}=\cos\frac{2\pi k}{n}+i\sin\frac{2\pi k}{n}\)。

（2）計算引理

顯然，因爲單位複數根循環(\(\omega_n^{zn}=e^{2\pi iz}=\left[\left(e^{\pi i}\right)^2\right]^z=1^z=1\))，有變換恆等式：\[\omega_n^k=\omega_n^{k+wn}(w\in\mathbb{Z})\]

每一份再分紅\(k\)份，編號也變成\(k\)倍，位置天然不變(\(\omega_{d n}^{d k}=e^\frac{2\pi i d k}{dn}=e^\frac{2\pi i k}{n}=\omega_n^k\))，因此有消去引理：\[\omega_{d n}^{d k}=\omega_n^k\]

因爲過了\(n/2\)就會繞過半圈(\(\omega_n^{n/2}=e^\frac{\pi i n}{n}=e^{\pi i}=-1\))，因此有折半引理：\[\omega_n^k=-\omega_n^{k\pm n/2}\]

對單位複數根求和，根據幾何級數(等差數列求和公式)，可得(\(\sum\limits_{k=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^n)^j-1}{\omega_n^1-1}=0\))，即有求和引理（要注意公式的使用條件）：\[\sum_{k=0}^{n-1}(\omega_n^k)^j=0,\omega_n\ne 1\]

三、DFT&FFT

1.DFT

DFT就是求多項式\(A(x)\)在點\((\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1})\)處取值的過程。即：\[y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}a_j\omega_n^{kj}\]

結果\(y=(y_0,y_1,y_2,\cdots,y_{n-1})\)就是\(a\)的離散傅里葉變換(DFT),記爲\(y=DFT_n(a)\)

2.FFT

（1）遞歸

DFT的\(O(n^2)\)算法太慢了，FFT使用分治策略優化速度到\(O(n\log n)\)。

考慮奇偶分治。

如今，咱們假設\(n=2^t\)，設原係數\(a=(a_0,a_1,a_2,\cdots,a_{n-1})\)，分治爲偶數部分\(a_1=(a_0,a_2,a_4,\cdots,a_{n-2})\)，奇數部分\(a_2=(a_1,a_3,a_5,\cdots,a_{n-1})\)，已經遞歸求出\(y_1=DFT_{n/2}(a_1)\)，\(y_2=DFT_{n/2}(a_2)\)，如今咱們要合併\(y_1,y_2\)，獲得\(y=DFT_n(a)\)(蝴蝶操做)。

對於\(n=1\)的邊界狀況，結果是顯然的：由於\(k=0\)，故\(\omega_1^0=1\)，即結果等於原係數。

對於\(n>1\)，如今咱們枚舉\(k\in[1,n]\)要合併出\(y_k\)：
\[\begin{align*} y_k&=A(\omega_n^k)\\ &=\sum_{j=0}^{n-1}a_j\omega_n^{kj}\\ &=\sum_{j=0}^{n/2-1}a_{2 j}\omega_n^{2 k j}+\sum_{j=0}^{n/2-1}a_{2 j+1}\omega_n^{2 k j+k}\\ &=\sum_{j=0}^{n/2-1}a_{2 j}\omega_n^{2 k j}+\omega_n^k\sum_{j=0}^{n/2-1}a_{2 j+1}\omega_n^{2 k j}\\ &=\sum_{j=0}^{n/2-1}a_{1_j}\omega_{n/2}^{k j}+\omega_n^k\sum_{j=0}^{n/2-1}a_{2_j}\omega_{n/2}^{k j}\tag{消去引理} \end{align*}\]

• 對於\(k<n/2\)：
  \[\begin{align*} y_k&=y_{1_k}+\omega_n^k y_{2_k} \end{align*}\]

• 對於\(k\geq n/2\)：
  \[\begin{align*} y_k&=\sum_{j=0}^{n/2-1}a_{1_j}\omega_{n/2}^{(k-n/2)j}+\omega_n^k\sum_{j=0}^{n/2-1}a_{2_j}\omega_{n/2}^{(k-n/2)j}\tag{變換恆等式}\\ &=y_{1_{k-n/2}}+\omega_n^k y_{2_{k-n/2}}\\ &=y_{1_{k-n/2}}-\omega_n^{k-n/2}y_{2_{k-n/2}}\tag{折半引理} \end{align*}\]

咱們用\(k+n/2\)替代\(k\)，就獲得\(y_{k+n/2}=y_{1_k}-\omega_n^k y_{2_k}\)

結合在一塊兒就獲得\(\begin{cases}y_k&=y_{1_k}+\omega_n^k y_{2_k}\\y_{k+n/2}&=y_{1_k}-\omega_n^k y_{2_k}\end{cases}\)
這樣咱們就能夠把兩個\(n/2\)長的序列合併爲一個\(n\)長的序列了。

如下圖的遞歸序，就能夠在\(O(n\log n)\)的時間複雜度內完成求解了。

（2）迭代

遞歸方法消耗內存、時間過大，沒法承受。咱們每次把下標分爲奇數部分和偶數部分，是否有辦法直接求出最後的遞歸運算順序，以免遞歸呢？

這樣想：
第一次奇偶劃分，咱們按照二進制的倒數第一位排序；
第二次奇偶劃分，咱們按照二進制的倒數第二位排序；
第\(n\)次奇偶劃分，咱們按照二進制的倒數第\(n\)位排序；
依此類推。

所以，結果順序就是原序列按照二進制位翻轉的大小排序的結果。只要依次交換\(a_k,a_{rev(k)}\)，求出序列，就能夠用迭代方法相鄰歸併實現快速傅里葉變換。

或者，咱們也能夠用更加代數的方法來發現這個結論。
已知如今位置爲\(k=(b_1 b_2 b_3 \cdots b_n)_2\)，按照奇偶重排：

• 若\(b_n=0\)，則位置變爲\(\frac{k}{2}=(0 b_1 b_2 \cdots b_{n-1})_2=(b_n b_1 b_2 \cdots b_{n-1})_2\)，即爲把最後一位提到第一位。
• 若\(b_n=1\)，則位置變爲\(2^{n-1}-1+\frac{k+1}{2}=\frac{2^n+k-1}{2}=\frac{(1 b_1 b_2 \cdots b_{n-1} 0)_2}{2}=(b_n b_1 b_2 \cdots b_{n-1})_2\)，一樣是把最後一位提到第一位。

則反覆\(n\)次以後，就至關於二進制反轉了。

如\(n=8\)時，求出二進制：

0	1	2	3	4	5	6	7
\(000_2\)	\(001_2\)	\(010_2\)	\(011_2\)	\(100_2\)	\(101_2\)	\(110_2\)	\(111_2\)

翻轉：

0	1	2	3	4	5	6	7
\(000_2\)	\(100_2\)	\(010_2\)	\(110_2\)	\(001_2\)	\(101_2\)	\(011_2\)	\(111_2\)

按翻轉後的值排序：

0	4	2	6	1	5	3	7
\(000_2\)	\(001_2\)	\(010_2\)	\(011_2\)	\(100_2\)	\(101_2\)	\(110_2\)	\(111_2\)

這樣就能夠把奇偶合並轉化爲左右歸併，迭代實現了。

五、IDFT&IFFT

何把點值表達變回係數表達呢？若是把求值寫成矩陣形式，就是：
\[\begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}\]

若是咱們要把\(y\)變成\(a\)，就須要求出第一個矩陣的逆，即：
\[\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix}^{-1} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}\]

這個範德蒙德矩陣極爲特殊，它的逆矩陣是：
\[\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}=\frac{1}{n} \begin{bmatrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0 \\ \omega_n^0 & \omega_n^{-1} & \omega_n^{-2} & \omega_n^{-3} & \cdots & \omega_n^{-(n-1)} \\ \omega_n^0 & \omega_n^{-2} & \omega_n^{-4} & \omega_n^{-6} & \cdots & \omega_n^{-2(n-1)} \\ \omega_n^0 & \omega_n^{-3} & \omega_n^{-6} & \omega_n^{-9} & \cdots & \omega_n^{-3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \omega_n^{-3(n-1)} & \cdots & \omega_n^{-(n-1)^2} \end{bmatrix} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}\]

只是把每項取倒數，並將結果除以\(n\)便可。

證實：

記原矩陣爲\(A_{n n}\),咱們給出的矩陣爲\(B_{n n}\)。

\(\therefore A_{i j}=\omega_n^{i j},B_{i j}=\frac{1}{n}\omega_n^{-i j}(0\le i,j\le n-1)\)

\[\begin{align*} AB_{i j}&=\sum_{k=0}^{n-1} A_{i k}B_{k j}\\ &=\frac{1}{n}\sum_{k=0}^{n-1} \omega_n^{i k}\omega_n^{-k j}\\ &=\frac{1}{n}\sum_{k=0}^{n-1} \omega_n^{i k-k j}\\ &=\frac{1}{n}\sum_{k=0}^{n-1} (\omega_n^{i-j})^k \end{align*}\]

• 當\(i=j\)時：

\[\begin{align*} AB_{i j}&=\frac{1}{n}\sum_{k=0}^{n-1}1^k\\ &=\frac{1}{n}\times n\\ &=1 \end{align*}\]

• 當\(i\ne j\)時：

\[\begin{align*} AB_{i j}&=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{i-j})^k\\ &=0\tag{求和引理} \end{align*}\]

綜上所述，\(AB=I_n\)，即\(B=A^{-1}\)

以上，離散傅里葉逆變換(IDFT)的表達式爲：\[a_k=\frac{1}{n}\sum_{j=0}^{n-1}y_k\omega_n^{-k j}\]
記\(a=IDFT_n(y)\)。

同理，能夠用相同的方法把IDFT加速到\(O(n\log n)\),稱爲IFFT。

五、NTT

1.定義

有時候咱們會想要模素數\(p\)意義下的多項式乘法。此時，由次數界爲\(n\)的多項式\(A(x),B(x)\)的係數表達\(a,b\)求\(S(x)=A(x)B(x)\)的係數表達\(s\)的公式爲：\[s_k=\sum_{j=0}^k a_j b_{k-j}\mod p\]

FFT無能爲力，咱們須要一種新的DFT，以數論的辦法進行，這就是快速數論變換(NTT)。

2.原根

（1）定義

咱們須要一種有數論循環性質的新數，原根剛好知足咱們的要求。

設\(m\)爲正整數，\(a\)爲整數，若\(a\)模\(m\)的階等於\(\varphi(m)\)，則稱\(a\)爲模\(m\)的一個原根。

假設\(g\)是素數\(p\)的原根，有\(1<g<p\)，且對於\(k=0,1,2,\cdots,p-1\)，有\(g^k\mod p\)的結果兩兩不一樣，且\(g^{p-1}\equiv 1 \pmod p\)。

能夠發現，原根一樣有循環的性質。所以，咱們類比\(\omega_n^k\)的定義，把原來的\(\omega_n^k=e^\frac{2\pi ik}{n}\)替換爲\(g^\frac{k(p-1)}{n}\)。

（2）性質

咱們來證實一些相似單位複數根的性質。

變換恆等式：
由於：\[g^{p-1}\equiv 1\]

因此：\[g^\frac{(k+n)(p-1)}{n}\equiv g^\frac{k(p-1)}{n}g^{p-1}\equiv g^\frac{k(p-1)}{n}\]

---

消去引理：
顯然有：\[g^\frac{d k(p-1)}{d n}\equiv g^\frac{k(p-1)}{n}\]

---

折半引理：
由於：\[g^\frac{(k+\frac{n}{2})(p-1)}{n}\equiv g^\frac{k(p-1)}{n}g^\frac{p-1}{2}\]

因此若要使：\[g^\frac{k(p-1)}{n}+g^\frac{(k+\frac{n}{2})(p-1)}{n}\equiv 0\pmod p\]

成立，即：\[g^\frac{k(p-1)}{2}\left(1+g^\frac{p-1}{2}\right)\equiv 0\]

只需證：\[g^\frac{p-1}{2}\equiv p-1\]

因爲：\[g^{k}\equiv 0,1,2,\cdots,p-1\]

那麼咱們能夠設：\[g^\frac{p-1}{2}\equiv x,x=0,1,2,\cdots,p-1\]

由於：\[\left(g^\frac{p-1}{2}\right)^2\equiv g^{p-1}\equiv 1\]

因此：\[x^2-1\equiv 0\]

即：\[(x+1)(x-1)\equiv 0\]

又由於\(p\)爲素數，因此有：\[x+1\equiv 0\;或\;x-1\equiv 0\]

因此：\[x=1,p-1\]

又由於：\[g^{p-1}\equiv 1,g^k\mod p兩兩不一樣\]

因此：\[x=p-1\]

即：\[g^\frac{p-1}{2}\equiv p-1\]

得證：\[g^\frac{k(p-1)}{n}+g^\frac{(k+\frac{n}{2})(p-1)}{n}\equiv 0\]

---

求和引理：
\[\sum\limits_{k=0}^{n-1}g^\frac{k j(p-1)}{n}\equiv\sum\limits_{k=0}^{n-1}\left(g^\frac{j(p-1)}{n}\right)^k\equiv\frac{g^\frac{n j(p-1)}{n}-1}{g^\frac{j(p-1)}{n}-1}\equiv\frac{\left(g^{p-1}\right)^j-1}{g^\frac{j(p-1)}{n}-1}\equiv 0\]

這樣，咱們就證實了原根因爲複數單位根相同的性質。

3.公式

咱們用原根替換複數單位根，獲得：\[y_k=A(g^\frac{k(p-1)}{n})=\sum_{j=0}^{n-1}a_k g^\frac{k j(p-1)}{n}\]

即\(y=NTT_n(a)\)。逆變換：\[a_k=\frac{1}{n}\sum_{j=0}^{n-1}y_k g^\frac{-k j(p-1)}{n}\]

即\(a=INTT_n(y)\)。

六、其餘擴展

1.任意模數FFT

有的時候咱們須要卷積後模上一個數，這個數不是NTT模數，甚至可能不是個質數。那咱們該怎麼作呢？

這裏使用拆係數FFT，本質是以時間換精度。

如今給定次數界爲\(m\)的兩個多項式\(A(x),B(x)\)，要求\(A(x)B(x) \mod P\)。

首先，令\(M=\lfloor\sqrt{P}\rfloor\)，再對於每一個\(a_i\)或\(b_i\)，把其變爲\(k_i M+r_i(r_i<M)\)的形式。這樣，\(k_i\)和\(r_i\)就都小於等於\(M\)了。

那麼卷積就能夠寫成：\[c_i=\sum_{k=0}^i a_k b_{i-k}=\sum_{k=0}^i(k_{a_i}M+r_{a_i})(k_{b_{i-k}}M+r_{b_{i-k}})=M^2\sum_{k=0}^i k_{a_i}k_{b_{i-k}}+M\sum_{k=0}^i (k_{a_i}r_{b_{i-k}}+r_{a_i}k_{b_{i-k}})+\sum_{k=0}^i r_{a_i}r_{b_{i-k}}\]

那麼咱們看到，\(c_i\)能夠由\(k\)和\(r\)合併獲得。那麼咱們對\(k_a,k_b,r_a,r_b\)分別作FFT，再對應算出\(x_1=k_a k_b,x_2=k_a r_b+r_a k_b,x_3=r_a r_b\)，對它們再分別作IFFT，就能夠由\(c=M^2 x_1+M x_2+x_3\)獲得答案了。

這麼作的好處究竟在哪裏呢？事實上，計算後的最大值由\(m P\)變爲了\(m \lfloor\sqrt{P}\rfloor\)，避免了溢出。

時間複雜度：\(O(n\log n)\)，常數爲\(7\)。

2.多項式求逆

如今咱們有一個次數界爲\(n\)的多項式\(A(x)\)，要求\(B(x)\)知足\(A(x)B(x)\equiv 1\pmod{x^n}\)。

咱們考慮倍增實現。

• 當\(n=1\)時，直接求逆元求得\(B(x)\equiv A(x)^{p-2}\)。
• 當\(n>1\)時，已有\(A(x)G(x)\equiv 1\pmod{x^\frac{n}{2}}\)：

由於：\[A(x)B(x)\equiv 1\pmod{x^n}\]

又由於除了\(0\)次項以外到\(n-1\)次都爲\(0\)，所以到\(\frac{n}{2}-1\)次項也爲零：\[A(x)B(x)\equiv 1\pmod{x^\frac{n}{2}}\]

又\[A(x)G(x)\equiv 1\pmod{x^\frac{n}{2}}\]

兩式相減：\[A(x)[B(x)-G(x)]\equiv 0\pmod{x^\frac{n}{2}}\]

由於：\[A(x)\ne 0\]

因此：\[B(x)-G(x)\equiv 0\pmod{x^\frac{n}{2}}\]

既然\(0\)至\(\frac{n}{2}-1\)次項都爲零，那麼平方以後\(0\)至\(n-1\)次項也爲零：\[B(x)^2-2 B(x)G(x)+G(x)^2\equiv 0\pmod{x^n}\]

又\[A(x)B(x)\equiv 1\pmod{x^n}\]

兩邊同時乘上\(A(x)\)，得：\[B(x)-2 G(x)+A(x)G(x)^2\equiv 0\pmod{x^n}\]

移項，得：\[B(x)\equiv 2 G(x)-A(x)G(x)^2\pmod{x^n}\]

這樣就能夠了。

時間複雜度證實：

顯然有遞歸式：
\[\begin{cases}T(0)=1\\T(n)=T(\frac{n}{2})+n\log_2(n)\end{cases}\]

展開可得：\[\begin{align*} T(n)&=\sum_{i=0}^{\log_2(n)}2^i \log_2(2^i)\\ &=\sum_{i=0}^{\log_2(n)}2^i i \end{align*}\]

即咱們要求和式：\[S_n=\sum_{k=0}^n 2^k k\]

的值。

對和式用擾動法，得：\[\begin{align*} S_n+(n+1)2^{n+1}&=\sum_{k=0}^n (k+1)2^{k+1}\\ S_n+(n+1)2^{n+1}&=2\sum_{k=0}^n k 2^k+2\sum_{k=0}^n 2^k\\ S_n+(n+1)2^{n+1}&=2S_n+2^{n+2}-2\\ -S_n&=2^{n+2}-2-(n+1)2^{n+1}\\ S_n&=2((n+1)2^n-2^{n+1}+1)\\ S_n&=2(2^n n-2^n+1) \end{align*}\]

代入，得：\[\begin{align*} T(n)&=S_{\log_2(n)}\\ &=2(2^{\log_2(n)}\log_2(n)-2^{\log_2(n)}+1)\\ &=2(n\log_2(n)-n+1)\\ \end{align*}\]

則時間複雜度爲：\[O(T(n))=O(n \log n)\]

由於多項式乘法的常數爲\(3\)，所以多項式求逆的常數爲\(2\times3=6\)。

3.多項式開根

前置：多項式求逆。

如今咱們有一個次數界爲\(n\)的多項式\(A(x)\)，要求\(B(x)\)知足\(B(x)^2\equiv A(x)\pmod{x^n}\)。

仍是倍增。

• 當\(n=1\)時，\(B(x)\)等於\(A(x)\)在模意義下的開根。
• 當\(n>1\)時，已有\(G(x)^2\equiv A(x)\pmod{x^\frac{n}{2}}\)：

由於：\[G(x)^2\equiv A(x)\pmod{x^\frac{n}{2}}\]

移項，得：\[G(x)^2-A(x)\equiv 0\pmod{x^\frac{n}{2}}\]

兩邊平方，同理可得：\[G(x)^4-2 G(x)^2 A(x)+A(x)^2\equiv 0\pmod{x^n}\]

因此：\[[G(x)^2+A(x)]^2-4 G(x)^2 A(x)\equiv 0\pmod{x^n}\]

即：\[[G(x)^2+A(x)]^2\equiv 4 G(x)^2 A(x)\pmod{x^n}\]

除過去：\[\frac{[G(x)^2+A(x)]^2}{4 G(x)^2}\equiv A(x)\pmod{x^n}\]

獲得：\[A(x)\equiv\left(\frac{G(x)^2+A(x)}{2 G(x)}\right)^2\equiv B(x)^2\pmod{x^n}\]

即：\[B(x)\equiv\frac{G(x)^2+A(x)}{2 G(x)}\equiv\frac{A(x)}{2 G(x)}+\frac{G(x)}{2}\pmod{x^n}\]

這就能夠了。

時間複雜度：\(O(n\log n)\)，常數爲\(2\times(6+3)=18\)。

4.多項式求導

根據導數的可加性和冪函數求導法則\(\frac{\mathbb{d}(cx^k)}{\mathbb{d}x}=c k x^{k-1}\)，有多項式的導數爲：\[\frac{\mathbb{d}(\sum\limits_{k=0}^{n-1}a_k x^k)}{\mathbb{d} x}=\sum_{k=0}^{n-1}\frac{\mathbb{d}(a_k x^k)}{\mathbb{d} x}=\sum_{k=1}^{n-1}k a_k x^{k-1}=\sum_{k=0}^{n-2}(k+1)a_{k+1}x^k\]

時間複雜度：\(O(n)\)，常數爲\(1\)。

5.多項式積分

根據積分的可加性和冪函數的不定積分\(\displaystyle\int c x^k\mathbb{d}x=\frac{c}{k}x^{k+1}\)，有多項式的不定積分爲：\[\int\sum_{k=0}^{n-1}a_k x^k \mathbb{d}x=\sum_{k=0}^{n-1}\int a_k x^k\mathbb{d}x=\sum_{k=0}^{n-1}\frac{a_k}{k+1}x^{k+1}=\sum_{k=1}^n\frac{a_{k-1}}{k}x^k\]

時間複雜度：\(O(n)\)，常數爲\(1\)。

6.多項式求ln

前置：多項式求導&積分&求逆

如今已知多項式\(A(x)\)，要求\(B(x)=\ln A(x)\)。咱們兩邊微分，獲得：\[B'(x)=\frac{\mathbb{d}(\ln A(x))}{\mathbb{d} A(x)}\frac{\mathbb{d} A(x)}{\mathbb{d} x}\]
\[B'(x)=\frac{A'(x)}{A(x)}\]

再兩邊積分，就獲得：\[B(x)=\int\frac{A'(x)}{A(x)}\mathbb{d}x\]

所以，咱們直接多項式求導+求逆+積分解決。

時間複雜度：\(O(n\log n)\)，常數爲\(6+3=9\)。

7.多項式求exp

前置：多項式求ln

如今，咱們已知多項式\(A=A(x)\)（這樣寫是由於在這裏把\(A\)視爲是與\(x\)無關的），求\(F(x)=\exp(A)=e^{A}\)。只要咱們設\(G(x)=\ln x-A\)，就獲得：\[G(F(x))=\ln F(x)-A=0\]

咱們考慮用牛頓迭代法來倍增解這個方程。

對於牛頓迭代法的初始解，即結果的常數項，咱們並不知道具體值。可是若是不對的話，也只是缺乏了一個常數罷了，那咱們不妨設\(F(x)=1\)。

倍增：如今設咱們已經求出了\(F(x)\)的前\(n\)項\(F_0(x)\)，即：\[F_0(x)\equiv F(x)\pmod{x^n}\]

根據牛頓迭代法，咱們求出下一個近似解爲：\[F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-\frac{\ln F_0(x)-A}{F_0(x)^{-1}}\equiv F_0(x)(1-\ln F_0(x)+A(x))\]

如此，就能夠倍增實現了。

時間複雜度：\(O(n\log n)\)，常數：\(2\times(9+3)=24\)。

8.多項式求冪

已知多項式\(A(x)\)和指數\(k\)，求\(A(x)^k\)。

在冪數很大的時候，若是使用快速冪的思想，時間複雜度爲\(O(n \log n\log k)\)，常數爲\(6\)。當\(k\)很大時速度很慢，咱們必須進行優化。

咱們發現：\[A(x)^k=\left(e^{\ln A(x)}\right)^k=e^{k \ln A(x)}=\exp(k \ln A(x))\]

因而咱們能夠用多項式求ln+多項式求exp求出。

時間複雜度：\(O(n\log n)\)，常數：\(9+24=33\)。

9.時間複雜度總結

內容	時間複雜度	常數
多項式卷積	\(O(n\log n)\)	\(3\)
多項式求逆	\(O(n\log n)\)	\(6\)
多項式開根	\(O(n\log n)\)	\(18\)
多項式求導	\(O(n)\)	\(1\)
多項式積分	\(O(n)\)	\(1\)
多項式求ln	\(O(n\log n)\)	\(9\)
多項式求exp	\(O(n\log n)\)	\(24\)
多項式求冪	\(O(n\log n)\)	\(33\)

七、代碼實現

1.二進制反轉

能夠用一種相似dp的思想計算。

邊界：\(0\)的二進制翻轉爲\(0\)

遞歸式：對於\(a\)，已經算出了\(rev_{\lfloor\frac{a}{2}\rfloor}\)，\(a\)就是除去最後一位的二進制翻轉(\(rev_{\lfloor\frac{a}{2}\rfloor}\))向後移動一位再補上第一位。即：

rev[a]=(rev[a>>1]>>1)|((a&1)<<(l-1))

\(l\)爲要翻轉的位數。

2.複數類

套公式便可。

struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
};

3.FFT

• 1：根據二進制翻轉交換\(a_r\)和\(a_{rev(r)}\)
• 2：枚舉歸併步長\(i\in[1,n)\)，\(i\)爲二的冪；
  • 2.1： 根據歐拉公式求出\(\omega_{2 i}^1=\cos\frac{\pi}{i}+i\sin\frac{\pi}{i}\)
  • 2.2：枚舉歸併位置\(j\)，歸併\([j,j+i)\)和\([j+i,j+2 i)\)，步長爲\(2 i\)
    • 2.2.1：枚舉\(x\)的冪數\(k\in[0,i)\)進行蝴蝶操做計算\(y\)，根據單位根計算\(\omega_i^k\)
      • 2.2.1.1：\(y_{j+k}=y_{j+k}+\omega_{2 i}^k y_{j+k+i}\)
      • 2.2.1.2：\(y_{j+k+i}=y_{j+k}-\omega_{2 i}^k y_{j+k+i}\)

注意：因爲在C++中值會被更改，所以須要引入臨時變量。

void FFT(complex c[]){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}

4.IFFT

因爲公式中只差一個負號而已，所以引入一個參數\(type\)，在歐拉公式的地方乘上去，再除以\(n\)就能夠了。

void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=m;i++)a[i].re/=n;
}

5.多項式乘法

模板：洛谷3083

注意：

一、因爲FFT要求\(n\)爲\(2\)的冪且結果的次數界較大，因此要把兩個因式的係數補到\(l\)位，\(l\)知足\(l=2^t\)且\(l/2\)大於等於因式的次數界。

二、\(FFT\)雖然在數學上精準，但在C++中偏差巨大，所以虛部不會爲\(0\)，忽略便可。實部也不爲正數，能夠加上\(0.1\)再向下取整。

代碼：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const double pi=acos(-1);
struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}a[2097153],b[2097153];
int n,m,l,r[2097153];
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=m;i++)a[i].re/=n;
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&a[i].re);
    for(int i=0;i<=m;i++)scanf("%lf",&b[i].re);
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a);
    FFT(b);
    for(int i=0;i<=n;i++)a[i]=a[i]*b[i];
    IFFT(a);
    for(int i=0;i<=m;i++)printf("%d ",int(a[i].re+0.5));
}

6.(I)FFT+(I)NTT

給出一個\(n\)次多項式和一個\(m\)次多項式的係數表達(\(1\le n,m\le 4000000\))，求乘積。\(type=0\)時，直接計算；\(type=1\)時，結果取模\(998244353\)(原根爲\(3\))。

注：爲了方便閱讀，代碼效率不高。若要提速，能夠把單位根打表。並且，因爲\(g^\frac{p-1}{n}\)中\(\frac{p-1}{n}\)必須爲整數，故僅在第一個比\(n+m\)大的\(2\)的整數次冪是\(p-1\)的約數時纔可行，此處\(998244353-1=998244352=119\times2^{23}\)，\(2^{23}=8388608>n+m\)。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const double pi=acos(-1);
const int p=998244353;
typedef long long ll;
struct complex{
    double re,im;
    complex(double re=0,double im=0):re(re),im(im){}
    complex operator+(complex b)const{return complex(re+b.re,im+b.im);}
    complex operator-(complex b)const{return complex(re-b.re,im-b.im);}
    complex operator*(complex b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}af[131073],bf[131073];
struct modp{
    int a;
    modp(int a=0):a(a){}
    modp operator+(modp b)const{return (a+b.a)%p;}
    modp operator-(modp b)const{return (a-b.a+p)%p;}
    modp operator*(modp b)const{return ll(a)*b.a%p;}
    modp operator/(modp b)const{return (b^(p-2))*a;}
    modp operator^(int b)const{
        modp ans=1,bs=a;
        while(b){
            if(b&1)ans=ans*bs;
            bs=bs*bs;
            b>>=1;
        }
        return ans;
    }
}an[131073],bn[131073];
const modp g=3;
int n,m,l,r[131073],type;
modp gn[18];
void init(){
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<=17;i++)gn[i]=g^((p-1)/(1<<i));
}
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=n;i++)c[i].re/=n;
}
void NTT(modp c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1,id=1;i<n;i<<=1,id++){
        modp ggn=gn[id]^(tp==1?1:p-2);
        for(int j=0;j<n;j+=(i<<1)){
            modp gg=1;
            for(int k=0;k<i;k++,gg=gg*ggn){
                modp x=c[j+k],y=gg*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void INTT(modp c[]){
    NTT(c,-1);
    for(int i=0;i<=n;i++)c[i]=c[i]/n;
}
int main(){
    scanf("%d%d%d",&n,&m,&type);
    if(type==0){
        for(int i=0;i<=n;i++)scanf("%lf",&af[i].re);
        for(int i=0;i<=m;i++)scanf("%lf",&bf[i].re);
    }else{
        for(int i=0;i<=n;i++)scanf("%d",&an[i].a);
        for(int i=0;i<=m;i++)scanf("%d",&bn[i].a);
    }
    init();
    if(type==0){
        FFT(af);
        FFT(bf);
        for(int i=0;i<=n;i++)af[i]=af[i]*bf[i];
        IFFT(af);
        for(int i=0;i<=m;i++)printf("%d ",int(af[i].re+0.5));
    }else{
        NTT(an);
        NTT(bn);
        for(int i=0;i<=n;i++)an[i]=an[i]*bn[i];
        INTT(an);
        for(int i=0;i<=m;i++)printf("%d ",an[i].a);
    }
    printf("\n");
}

常數只有上面那個的三分之一的NTT（正式考試請務必採用這種寫法）：

PS:有一道題上面那個3700ms,這個1000ms

inline ll pow(ll a,int b){
    ll ans=1;
    while(b){
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
inline ll add(ll a,ll b){return a+b>p?a+b-p:a+b;}
inline ll cut(ll a,ll b){return a-b<0?a-b+p:a-b;}
void init(){
    for(n=1;n<=s;n<<=1);
    nn=n;
    gn[0][0]=gn[1][0]=1;
    gn[0][1]=pow(g,(p-1)/(n<<1));
    gn[1][1]=pow(gn[0][1],p-2);
    for(int i=2;i<(n<<1);i++){gn[0][i]=gn[0][i-1]*gn[0][1]%p;gn[1][i]=gn[1][i-1]*gn[1][1]%p;}
    inv[1]=1;
    for(int i=2;i<=(n<<1);i++)inv[i]=inv[p%i]*(p-p/i)%p;
}
void NTT(ll c[],int n,int tp=1){
    for(int i=0;i<n;i++){
        r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        if(i<r[i])swap(c[i],c[r[i]]);
    }
    for(int i=1;i<n;i<<=1){
        for(int j=0;j<n;j+=(i<<1)){
            for(int k=0;k<i;k++){
                ll x=c[j+k],y=gn[tp!=1][nn/i*k]*c[j+k+i]%p;
                c[j+k]=add(x,y);
                c[j+k+i]=cut(x,y);
            }
        }
    }
}
void INTT(ll c[],int n){
    NTT(c,n,-1);
    for(int i=0;i<n;i++)c[i]=c[i]*inv[n]%p;
}

7.任意模數FFT

給定\(n,m,P\)，再給定次數界爲\(n\)的第一個多項式和次數界爲\(m\)的第二個多項式，求兩個多項式的卷積模\(P\)。

注意：拆係數FFT精度損失極大，須要使用long double保證正確性。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long ll;
int x,n,m,l,r[262145],pm,p;
const long double pi=acos(-1);
struct complex{
    long double re,im;
    complex(long double re=0,long double im=0):re(re),im(im){}
    complex operator+(const complex &b)const{return complex(re+b.re,im+b.im);}
    complex operator-(const complex &b)const{return complex(re-b.re,im-b.im);}
    complex operator*(const complex &b)const{return complex(re*b.re-im*b.im,re*b.im+im*b.re);}
}k1[262145],r1[262145],k2[262145],r2[262145],c1[262145],c2[262145],c3[262145];
void FFT(complex c[],int tp=1){
    for(int i=0;i<n;i++)if(i<r[i])swap(c[i],c[r[i]]);
    for(int i=1;i<n;i<<=1){
        complex omn(cos(pi/i),tp*sin(pi/i));
        for(int j=0;j<n;j+=(i<<1)){
            complex om(1,0);
            for(int k=0;k<i;k++,om=om*omn){
                complex x=c[j+k],y=om*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void IFFT(complex c[]){
    FFT(c,-1);
    for(int i=0;i<=n;i++)c[i].re/=n;
}
void init(){
    m+=n;
    for(n=1;n<=m;n<<=1)l++;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
int main(){
    scanf("%d%d%d",&n,&m,&p);
    pm=sqrt(p);
    for(int i=0;i<=n;i++){
        scanf("%d",&x);
        k1[i]=x/pm;
        r1[i]=x%pm;
    }
    for(int i=0;i<=m;i++){
        scanf("%d",&x);
        k2[i]=x/pm;
        r2[i]=x%pm;
    }
    init();
    FFT(k1);
    FFT(r1);
    FFT(k2);
    FFT(r2);
    for(int i=0;i<=n;i++){
        c1[i]=k1[i]*k2[i];
        c2[i]=k1[i]*r2[i]+r1[i]*k2[i];
        c3[i]=r1[i]*r2[i];
    }
    IFFT(c1);
    IFFT(c2);
    IFFT(c3);
    for(int i=0;i<=m;i++){
        ll s1=ll(c1[i].re+0.5)%p*pm%p*pm%p,s2=ll(c2[i].re+0.5)%p*pm%p,s3=ll(c3[i].re+0.5)%p;
        printf("%lld ",((s1+s2)%p+s3)%p);
    }
}

7.多項式的運算

依賴關係：

直接按照公式打就行了。

咱們先修改NTT：

void NTT(modp c[],int n,int tp=1){
    for(int i=0;i<n;i++){
        r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
        if(i<r[i])swap(c[i],c[r[i]]);
    }
    for(int i=1,id=1;i<n;i<<=1,id++){
        modp ggn=gn[id]^(tp==1?1:p-2);
        for(int j=0;j<n;j+=(i<<1)){
            modp gg=1;
            for(int k=0;k<i;k++,gg=gg*ggn){
                modp x=c[j+k],y=gg*c[j+k+i];
                c[j+k]=x+y;
                c[j+k+i]=x-y;
            }
        }
    }
}
void INTT(modp c[],int n){
    NTT(c,n,-1);
    for(int i=0;i<n;i++)c[i]=c[i]/n;
}

求逆：

void inverse(modp c[],int n=n){
    static modp t[262145],tma[262145];
    t[0]=c[0]^(p-2);
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<(k<<1);i++)tma[i]=(i<k?c[i]:0);
        for(int i=(k>>1);i<(k<<1);i++)t[i]=0;
        NTT(tma,k<<1);
        NTT(t,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=t[i]*2-t[i]*t[i]*tma[i];
        INTT(t,k<<1);
    }
    memcpy(c,t,sizeof(int)*n);
}

開根(\(c[0]=1\))：

void sqrt(modp c[],int n=n){
    static modp t[262145],tma[262145],tmb[262145];
    t[0]=1;
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<k;i++)tma[i]=t[i]*2;
        inverse(tma,k);
        for(int i=0;i<(k<<1);i++)tmb[i]=(i<k?c[i]:0);
        NTT(tma,k<<1);
        NTT(tmb,k<<1);
        for(int i=0;i<(k<<1);i++){
            modp tmp=tma[i];
            tma[i]=t[i];
            t[i]=tmp*tmb[i];
        }
        INTT(t,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=(i<k?t[i]+tma[i]/2:0);
    }
    memcpy(c,t,sizeof(int)*n);
}

求導：

void derivative(modp c[],int n=n){for(int i=0;i<n;i++)c[i]=c[i+1]*(i+1);}

積分：

void integrate(modp c[],int n=n){for(int i=n-1;i>=1;i--)c[i]=c[i-1]*inv[i];c[0]=0;}

求ln：

void ln(modp c[],int n=n){
    static modp t[262145];
    for(int i=0;i<(n<<1);i++)t[i]=(i<n?c[i]:0);
    derivative(t,n);
    inverse(c,n);
    NTT(t,n<<1);
    NTT(c,n<<1);
    for(int i=0;i<(n<<1);i++)c[i]=c[i]*t[i];
    INTT(c,n<<1);
    for(int i=n;i<(n<<1);i++)c[i]=0;
    integrate(c,n);
}

求exp：

void exp(modp c[]){
    static modp t[262145],ta[262145];
    t[0]=1;
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<(k<<1);i++)ta[i]=t[i];
        ln(ta,k);
        for(int i=0;i<k;i++)ta[i]=c[i]-ta[i];
        ta[0].a++;
        NTT(t,k<<1);
        NTT(ta,k<<1);
        for(int i=0;i<(k<<1);i++)t[i]=t[i]*ta[i];
        INTT(t,k<<1);
        for(int i=k;i<(k<<1);i++)t[i]=0;
    }
    memcpy(c,t,sizeof(modp)*n);
}

求冪：
咱們在多項式求exp中假定了\(c[0]=1\)，那麼若是常數項不是\(1\)的話咱們就把常數項變爲\(1\)在運算後再用快速冪變回來便可。

void pow(modp c[],int k){
    ln(c);
    for(int i=0;i<n;i++)c[i]=c[i]*k;
    exp(c);
}

八、拉格朗日反演

1.形式冪級數

對於任意一個域\(F\)咱們定義在其上的形式冪級數爲：\[f(x)=\sum_{k=0}^\infty a_k x^k,a_k\in F\]
記全部的形式冪級數爲\(F[[x]]\)。

2.反演公式

拉格朗日反演是求關於函數方程的冪級數展開係數很是重要的工具，能夠用於組合計數函數的係數提取。

公式內容：

這裏\([x^n]f(x)\)指取\(f(x)\)中\(x^n\)的係數。

若\(f(x),g(x)\in F[[x]]\)且\(f(g(x))=x\)，則稱\(f(x)\)爲\(g(x)\)的複合逆。知足：\[[x^n]g(x)=\frac{1}{n}[x^{-1}]\frac{1}{f(x)^n}\tag{1}\]
特別的，若是\(f(x)=\displaystyle\frac{x}{\phi(x)}\)，那麼有：\[[x^n]g(x)=\frac{1}{n}[x^{n-1}]\phi(x)^n\tag{2}\]

公式的推導：
由式\(f(x)=\displaystyle\frac{x}{\phi(x)}\)，得\(\phi(x)=\displaystyle\frac{x}{f(x)}\)，代入\((2)\)式可得：\[[x^n]g(x)=\frac{1}{n}[x^{n-1}]\left(\frac{x}{f(x)}\right)^n\]

公式的推廣：
\[[x^n]h(g(x))=\frac{1}{n}[x^{n-1}]h'(x)\left(\frac{x}{f(x)}\right)^n\]