limit 完成 快速傅裏葉變換 而已 lin 循環 http 思想 !=

1、多項式的兩種表示法

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

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

2、復數

當\(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\)為弧長。

則根據歐拉公式，可將一個復數表示為一個二元組\((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^k\ne 1\]

3、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\)|
這樣就可以把奇偶合並轉化為左右歸並，叠代實現了。

5、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。

5、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)\)。

6、其他擴展

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)\)

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

這樣就可以了。

時間復雜度：\(O(n\log n)\)

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)\)

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)\)

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}x^{k+1}=\sum_{k=1}^n\frac{a_{k-1}}{k-1}x^k\]

時間復雜度：\(O(n)\)

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)\)

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)\)

8.多項式求冪

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

在冪數很大的時候，為了確定出系數，時間復雜度為\(O(k n \log (nk))\)，我們必須進行優化。

我們發現：\[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)\)的時間內求出。

7、代碼實現

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_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_i^k y_{j+k+i}\)
      • 2.2.1.2：\(y_{j+k+i}=y_{j+k}-\omega_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

註意：

1、由於FFT要求\(n\)為\(2\)的冪且結果的次數界較大，所以要把兩個因式的系數補到\(l\)位，\(l\)滿足\(l=2^t\)且\(l/2\)大於等於因式的次數界。

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);
}

8、拉格朗日反演

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

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