【ACM数论】从“快速傅里叶变换 FFT”到“快速数论变换 NTT”

1、FFT前言

FFT全称（Fast Fourier Transformation）快速傅里叶变换

NTT全称（Fast Number Theoretical Transformation）快速数论变换

1.1、引入问题

如何计算多项式$f(x)$与$g(x)$的乘积？

例如：计算$f(x)=1+x^2+x^3$与$g(x)=1+x^5+x^6+x^7$的乘积

通常而言，我们思考暴力的算法：
$$\begin{array}{rl}f(x)\times g(x)& =1\times (1+x^5+x^6+x^7)+x^2\times (1+x^5+x^6+x^7)+x^3\times (1+x^5+x^6+x^7)\\ & =1+x^5+x^6+x^7+x^2+x^7+x^8+x^9+x^3+x^8+x^9+x^{10}\\ & =1+x^2+x^3+x^5+x^6+2x^7+2x^8+2x^9+x^{10}\end{array}$$
显而易见，这个算法的时间复杂度是$O(n^2)$

那么有快速的计算方法吗？

当然是有的，于1965年由J.W.库利和T.W.图基提出FFT算法，可以在$O(nlogn)$的时间复杂度下计算出答案。

1.2、定义一些记号：

# $f(x)$表示一个多项式

**#**对于一个$n$次多项式$f(x)$的第$m$项，记为$f[m]$;换言之，$f(x)={\sum }_{i=0}^{n}f[i]x^i$

1.3、多项式乘法理解

例如：现有两个$n$次多项式$f(x)$与$g(x)$,需要计算这两个多项式的乘积.

不妨设$h(x)=f(x)\times g(x)$

那么通过1.1中的例子，我们可以发现，其实$h[k]={\sum }_{i=0}^{k}f[i]g[k-i]$

换言之，我们可以理解为：要想乘出$x^k$,需要$f(x)$的$x^i$项和$g(x)$的$x^{k-i}$项相乘，然后将他们的相乘的贡献累加便是$h[k]$

1.4、卷积

将形如:$h[k]={\sum }_{i\oplus j}f[i]g[j]$的式子称为卷积，其中⊕在数学中称为假设符号（假设什么运算）

例如，在多项式乘法中，为加法卷积：$h[k]={\sum }_{i+j}f[i]g[j]$

因此，对于卷积可以利用FFT构造多项式进行计算！

2、FFT思想

2.1、多项式的点值表达式

我们可以从

得出推论：在平面直角坐标系中，$n+1$个点就能唯一确定一个$n$次多项式。

因此，对于一个$n$次多项式$f(x)$而言，我们可以用$n+1$个点来描述一个多项式；换言之，两个多项式相乘，等价于两个多项式的对应点值相乘。(点值的乘法对应着整个多项式的乘法,也就是浓缩了整个多项式)

将能唯一表达一个$n$次多项式$f(x)$的$n+1$个点,称为$f(x)$的点值表达式。

$$
例如：n次多项式f(x)的点值表达式为：\

\left{(x_1,f(x_1)),(x_2,f(x_2)),(x_3,f(x_3))…(x_{n+1},f(x_{n+1}))\right}
$$
由于$n$次多项式$f(x)$与$n$次多项式$g(x)$,他们的乘积一定是一个$2n+1$次的多项式，所以我们只需要对$f(x)$和$g(x)$找到$2n+1$个对应的点，然后将他们对应点的纵坐标相乘，便可以确定多项式$f(x)\times g(x)$的点值表达式；

例如： 2 n + 1 次多项式 f ( x ) × g ( x ) 的点值表达式为： { ( x 1 , f ( x 1 ) g ( x 1 ) ) , ( x 2 , f ( x 2 ) g ( x 2 ) ) , ( x 3 , f ( x 3 ) g ( x 3 ) ) . . . . . . ( x 2 n + 1 , f ( x 2 n + 1 ) g ( x 2 n + 1 ) ) } 例如：2n+1次多项式f(x)\times g(x)的点值表达式为： \\ \left\{(x_1,f(x_1)g(x_1)),(x_2,f(x_2)g(x_2)),(x_3,f(x_3)g(x_3))......(x_{2n+1},f(x_{2n+1})g(x_{2n+1}))\right\} 例如：2n+1次多项式f(x)×g(x)的点值表达式为：{(x1​,f(x1​)g(x1​)),(x2​,f(x2​)g(x2​)),(x3​,f(x3​)g(x3​))......(x2n+1​,f(x2n+1​)g(x2n+1​))}
最后我们将多项式$f(x)\times g(x)$的点值表达式转化为多项式的一般表达式（系数表达式）便得到了我们想要的多项式乘积了。

而这也是FFT算法的流程

2.2、DFT&IDFT

DFT：把多项式的系数表达式转化为点值表达式的变换过程

IDFT：把多项式的点值表达式转化为系数表达式的变换过程

下面这张图便对比了FFT算法与暴力算法的区别

3、FFT的实现

3.1、DFT

​ 由于我们想要多项式的点值表达式，所以我们需要利用一些方法获得多项式的上的点；但是细想一下会发现，
$$\begin{gathered} 如果对于一个n次多项式f(x)，我们任选n+1个点的横坐标x_1,x_2,x_3....x_n,x_{n+1}, \\ 再分别计算其对应的纵坐标f(x_1),f(x_2),f(x_3)...f(x_n),f(x_{n+1})的时间复杂度是O(n^2) \\ (由于计算每个横坐标的时间复杂度是O(n),但是有n个点)； \end{gathered}$$
​ 这样一来，我们通过随便找$n+1$个点的方式，会让我们的时间复杂度和暴力没什么区别！

​ 那么有快速找到$n+1$个点的方法么？答案当然是有的！

​ 我们可以这样考虑问题：

3.1.1、Good idea 1：函数的奇偶性

​ 首先，我们不考虑一般多项式，只考虑一个简单的多项式$f(x)=x^2$

那么，如果我们要从$f(x)$中选择8个点，我们如何选择？

这里我们需要注意$f(x)$的一个优美的性质————$f(x)$是偶函数，所以我们计算出一个点$(x,f(x))$后，便可以立马知道点$(-x,f(x))$也一定是$y(x)$上的点。

那么，如果$f(x)=x^3$,应该如何处理喃？

显而易见，我们会发现$f(x)=x^3$是一个奇函数，它与偶函数的区别在于，奇函数上的点是关于原点对称的，换言之，如果我们计算出点$(x,f(x))$,那就可以知道点$(-x,-f(x))$也是多项式上的点！

既然这样，对于多项式$f(x)=1+2x^2+3x^3+4x^4+5x^5$

那么我可以将$f(x)$中按幂数分离，因此可得$f(x)=(1+2x^2+4x^4)+(3x^3+5x^5)$
$$\begin{gathered} 不妨设f_o(x)=1+2x+4x^2,f_e(x)=3x+5x^2 \\ 那么有：\begin{array}{rl}f(x)& =(1+2x^2+4x^4)+x(3x^2+5x^4)\\ & =f_o(x^2)+x{f}_e(x^2)\end{array} \end{gathered}$$
这样，对于$f_o(x^2),f_e(x^2)$都是关于$x^2$的多项式，我们便可以利用函数的奇偶性快速找到6个点，然后利用点值表达式的计算方式计算出该函数的点值表达式。

然后我们可以推广到一般多项式$f(x)=a_0+a_{1}x+a_2{x}^2+.....a_n{x}^n$
$$\begin{gathered} 可以设f(x)按照幂数分离，将偶幂数部分设为f_e(x^2),将奇幂数部分提出一个x后的部分设为f_o(x^2) \\ 所以：f(x)=f_e(x^2)+x{f}_o(x^2) \\ 然而对于，f_e(x)与f_o(x)而言，最高次都变成了\frac{n}{2}-1 \\ 同理可以推出重要结论\{\begin{array}{l}f(x_i)=f_e(x_i^2)+x_i{f}_o(x_i^2)\\ f(-x_i)=f_e(x_i^2)-x_i{f}_o(x_i^2)\end{array} \\ 这样一来，问题得到了简化： \\ 分别求f_e(x^2),f_o(x^2)在x_1^2,x_2^2....x_{\frac{n}{2}}^2的值,因为这样我们便可以得到\pm x_1,\pm x_2...\pm x_{\frac{n}{2}}这n个点 \end{gathered}$$
而对于$f_e(x)与f_o(x)$,我们可以重复$f(x)$的奇偶分离的方式；这样我们可以把对于$n$次多项式转化为求$\frac{n}{2}$次多项式（类似分治的思想），所以这个思路下，求多项式点值表达式的时间复杂度是$O(nlogn)$

因此，只要我们最后找到的点为合法的点，那么算法就是合法的！

但是这里出现了一个问题：对于$f(x)$中的$[\pm x_1,\pm x_2...\pm x\frac{n}{2}]$而言，每一个$\pm x_i$都是$\pm$配对的；但是对于$f_o(x^2)与f_e(x^2)$中的$[x_1^2,x_2^2..x_n^2]$而言,没有一个$x_i^2$是配对的;所以在实数域下，是找不到满足条件的点的！

不过解决问题的的方法总比问题多，所以我们考虑把数域扩展到复数域！

3.1.2、Good idea 2：复数域下的单位根

3.1.2.1解决疑惑：为何复数域下，满足3.1.1中的分治找点

一个简单例子：求多项式$f(x)=1+2x+3x^3$的点值表达式
$$\begin{gathered} 在这里我们不妨设这四个点的横坐标为:x_1,x_2,x_3,x_4 \\ 那么根据分治的规则:那么必须保证x_1=-x_2,x_3=-x_4 \\ \therefore x_1^2=x_2\dots \dots 2,x_3^2=x_4^2 \\ \because x_1^2=-x_3^2,那么可令x_1=1 \\ \therefore x_3^2=-1 \\ \because 在复数域下有i^2=-1 \\ \therefore x_3=i,x_4=-i \\ \therefore 对于横坐标x_1=1,x_2=-1,x_3=i,x_4=-i是满足条件的 \end{gathered}$$

3.1.2.2引入单位根

单位根：在复平面上，以单位圆点为起点，单位圆的 $n$ 等分点为终点，在单位圆上可以得到 $n$ 个复数，设幅角为正且最小的复数为$n$次单位根，记作${\omega }_n$

**#**推导${\omega }_n$
$$\begin{gathered} 由于在复平面当中，故俯角最小为\frac{2\pi }{n},所以可知在复平面中，{\omega }_n=cos\frac{2\pi }{n}+isin\frac{2\pi }{n} \\ 又\because 欧拉公式为e^{i\theta }=cos\theta +isin\theta \\ \therefore {\omega }_n=e^{i\frac{2\pi }{n}} \end{gathered}$$
**#**推导单位根的一些性质:

(1)$w_n^a\times w_n^b=w_n^{i+j}$
$$\begin{array}{rl}{w}_n^a\times w_n^b& =e^{ai\frac{2\pi }{n}}\times e^{bi\frac{2\pi }{n}}\\ & =e^{(a+b)i\frac{2\pi }{n}}\\ & =(e^{i\frac{2\pi }{n}}{)}^{a+b}\\ & =w_n^{a+b}\end{array}$$
(2)$w_{rn}^{rk}=w_n^k$
$$w_{rn}^{rk}=(e^{i\frac{2\pi }{rn}}{)}^{rk}=(e^{i\frac{2\pi }{rn}r}{)}^k=(e^{i\frac{2\pi }{n}}{)}^k=w_n^k$$
(3)$w_n^{2k}=w_{\frac{n}{2}}^k$
$$w_n^{2k}=e^{2ki\frac{2\pi }{n}}=e^{ki\frac{2\pi }{\frac{n}{2}}}=w_{\frac{n}{2}}^k$$
(4)如果$n$是偶数, $w_n^{k+\frac{n}{2}}=-w_n^k$
$$\begin{array}{rl}{w}_n^{k+\frac{n}{2}}& =(e^{i\frac{2\pi }{n}}{)}^{k+\frac{n}{2}}\\ & =e^{ki\frac{2\pi }{n}+i\pi }\\ & =e^{i\pi }{e}^{ki\frac{2\pi }{n}}\\ & =-1\times e^{ki\frac{2\pi }{n}}\\ & =-w_n^k\end{array}$$

3.1.2.3代入单位根

现已知$f(x)=f_e(x^2)+x{f}_o(x^2)$

(1)代入$x=w_n^k,k<\frac{n}{2}$
$$\begin{array}{rl}f(w_n^k)& =f_e(w_n^{2k})+w_n^k{f}_o(w_n^{2k})\\ & =f_e(w_{\frac{n}{2}}^k)+w_n^k{f}_o(w_{\frac{n}{2}}^k)\end{array}$$
(2)代入$x=w_n^{k+\frac{n}{2}},k<\frac{n}{2}$
$$\begin{array}{rl}f(w_n^{k+\frac{n}{2}})& =f_e(w_n^{2(k+\frac{n}{2})})+w_n^{k+\frac{n}{2}}{f}_o(w_n^{2(k+\frac{n}{2})})\\ & =f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})+w_n^{k+\frac{n}{2}}{f}_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ & =f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})-w_n^k{f}_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ & =f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})-w_n^k{f}_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ & =f_e(w_n^{2k+n})-w_n^k{f}_o(w_n^{2k+n})\\ & =f_e(w_n^{2k})-w_n^k{f}_o(w_n^{2k})\\ & =f_e(w_{\frac{n}{2}}^k)-w_n^k{f}_o(w_{\frac{n}{2}}^k)\end{array}$$

$$\begin{gathered} 那么得到（1）（2）两个结论，我们可以： \\ 如果知道多项式f_e(x),f_o(x)在w_{\frac{n}{2}}^0,w_{\frac{n}{2}}^1,w_{\frac{n}{2}}^2....w_{\frac{n}{2}}^{\frac{n}{2}-1}的点值表达式， \\ 利用（1）可以O(n)求出f(x)在w_{\frac{n}{2}}^0,w_{\frac{n}{2}}^1,w_{\frac{n}{2}}^2....w_{\frac{n}{2}}^{\frac{n}{2}-1}的点值表达式; \\ 利用（2）可以O(n)求出f(x)在w_{\frac{n}{2}}^{\frac{n}{2}},w_{\frac{n}{2}}^{\frac{n}{2}+1},w_{\frac{n}{2}}^{\frac{n}{2}+2}....w_{\frac{n}{2}}^{n-1}的点值表达式; \end{gathered}$$

既然这样，那么我们再3.1.1中$O(nlogn)$求多项式的思路是可行的！！！！

3.1.3：优化分治（蝴蝶变换）

我们用下面一张图来表示3.1.1中的分治过程：

显然,对于上述过程是可以使用递归的写法,但是递归的写法常数过大，

但是我们需要观察每个数最终的位置,最终发现对于幂数$i$而言，其实$i$最终在序列中的位置其实是$i$在二进制下翻转（蝴蝶变换）得到的数
$$\begin{gathered} 设i通过蝴蝶变换后的位置为reverse[i],i在二进制下的位数为len(i) \\ 那么有reverse[i]=\lfloor \frac{reverse\lfloor \frac{x}{2}\rfloor }{2}\rfloor +(i\times (mod2))\times (\frac{len(i)}{2}) \end{gathered}$$
因此我们可以把每个数都放在其最终位置上，不断往上迭代即可！

**#**附一个蝴蝶变换的小结论：

蝴蝶变换过后的序列，任取三个数都不构成等差数列。

3.2、IDFT

现在我们已经知道$k$次多项式$f(x)$的点值表达式：$\{(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),(w_n^2,f(w_n^2)),....,(w_n^n,f(w_n^n))\}(k\le n,n=2^m-1)$

那么可一把点值表达式改写成与之对应的矩阵的形式：
$$\left[\begin{array}{cccccc}& 1& (w_n^0{)}^1& (w_n^0{)}^2& \dots & (w_n^0{)}^n\\ & 1& (w_n^1{)}^1& (w_n^1{)}^2& \dots & (w_n^1{)}^n\\ & 1& (w_n^2{)}^1& (w_n^2{)}^2& \dots & (w_n^2{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^n{)}^1& (w_n^n{)}^2& \dots & (w_n^n{)}^n\end{array}\right]\times \left[\begin{array}{cc}& a_0\\ & a_1\\ & a_2\\ & ⋮\\ & a_n\end{array}\right]=\left[\begin{array}{cc}& f(w_n^0)\\ & f(w_n^1)\\ & f(w_n^2)\\ & ⋮\\ & f(w_n^n)\end{array}\right]$$
那么在这个矩阵乘法中：
$$\left[\begin{array}{cc}& a_0\\ & a_1\\ & a_2\\ & ⋮\\ & a_n\end{array}\right]即为多项式的系数表达式$$

$$那么我们只需要两边同时乘以\left[\begin{array}{cccccc}& 1& (w_n^0{)}^1& (w_n^0{)}^2& \dots & (w_n^0{)}^n\\ & 1& (w_n^1{)}^1& (w_n^1{)}^2& \dots & (w_n^1{)}^n\\ & 1& (w_n^2{)}^1& (w_n^2{)}^2& \dots & (w_n^2{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^n{)}^1& (w_n^n{)}^2& \dots & (w_n^n{)}^n\end{array}\right]的逆矩阵，便可以得到多项式的系数表达式$$

$$因为矩阵\left[\begin{array}{cc}& f(w_n^0)\\ & f(w_n^1)\\ & f(w_n^2)\\ & ⋮\\ & f(w_n^n)\end{array}\right]是一个关于f(w_n^i)已知的矩阵$$

$$\begin{gathered} 所以最后IDFT这一块的问题等价于求矩阵\left[\begin{array}{cccccc}& 1& (w_n^0{)}^1& (w_n^0{)}^2& \dots & (w_n^0{)}^n\\ & 1& (w_n^1{)}^1& (w_n^1{)}^2& \dots & (w_n^1{)}^n\\ & 1& (w_n^2{)}^1& (w_n^2{)}^2& \dots & (w_n^2{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^n{)}^1& (w_n^n{)}^2& \dots & (w_n^n{)}^n\end{array}\right]的逆矩阵 \\ 在这里我们不妨考虑，该矩阵的每一行： \\ 对于每j行有:[1(w_n^j{)}^1(w_n^j{)}^2\dots (w_n^j{)}^n]肯定有逆矩阵的第j列与之相乘 \\ \left[\begin{array}{cccccccc}& 1& (w_n^j{)}^1& (w_n^j{)}^2& \dots & (w_n^j{)}^p& \dots & (w_n^j{)}^n\end{array}\right]\times \left[\begin{array}{cc}& x_0\\ & x_1\\ & x_2\\ & ⋮\\ & x_q\\ & ⋮\\ & x_n\end{array}\right] \\ \therefore 对于每一项有w_n^{jp}\times x_q \\ 这个我们考虑单位根性质：w_n^j\times w_n^{n-j}=w_n^n=1 \\ \therefore 可以构造x_q=w_n^{(n-j)p} \\ \therefore w_n^{jp}\times x_q=w_n^{np}=1 \\ \therefore \left[\begin{array}{cc}& x_0\\ & x_1\\ & x_2\\ & ⋮\\ & x_q\\ & ⋮\\ & x_n\end{array}\right]=\left[\begin{array}{cc}& w_n^{(n-j)0}\\ & w_n^{(n-j)1}\\ & w_n^{(n-j)2}\\ & ⋮\\ & w_n^{(n-j)p}\\ & ⋮\\ & w_n^{(n-j)n}\end{array}\right] \end{gathered}$$

$$\begin{gathered} \therefore \left[\begin{array}{cccccccc}& 1& (w_n^j{)}^1& (w_n^j{)}^2& \dots & (w_n^j{)}^p& \dots & (w_n^j{)}^n\end{array}\right]\times \left[\begin{array}{cc}& w_n^{(n-j)0}\\ & w_n^{(n-j)1}\\ & w_n^{(n-j)2}\\ & ⋮\\ & w_n^{(n-j)p}\\ & ⋮\\ & w_n^{(n-j)n}\end{array}\right]=n \\ 又\because \forall d\in Z,d\ne 0且n-j+d>0有\left[\begin{array}{cccccccc}& 1& (w_n^j{)}^1& (w_n^j{)}^2& \dots & (w_n^j{)}^p& \dots & (w_n^j{)}^n\end{array}\right]\times \left[\begin{array}{cc}& w_n^{(n-j+d)0}\\ & w_n^{(n-j+d)1}\\ & w_n^{(n-j+d)2}\\ & ⋮\\ & w_n^{(n-j+d)p}\\ & ⋮\\ & w_n^{(n-j+d)n}\end{array}\right]=0 \end{gathered}$$

$$\therefore \left[\begin{array}{cccccc}& 1& (w_n^0{)}^1& (w_n^0{)}^2& \dots & (w_n^0{)}^n\\ & 1& (w_n^1{)}^1& (w_n^1{)}^2& \dots & (w_n^1{)}^n\\ & 1& (w_n^2{)}^1& (w_n^2{)}^2& \dots & (w_n^2{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^n{)}^1& (w_n^n{)}^2& \dots & (w_n^n{)}^n\end{array}\right]\times \left[\begin{array}{cccccc}& 1& (w_n^{n-0}{)}^1& (w_n^{n-0}{)}^2& \dots & (w_n^{n-0}{)}^n\\ & 1& (w_n^{n-1}{)}^1& (w_n^{n-1}{)}^2& \dots & (w_n^{n-1}{)}^n\\ & 1& (w_n^{n-2}{)}^1& (w_n^{n-2}{)}^2& \dots & (w_n^{n-2}{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^{n-n}{)}^1& (w_n^{n-n}{)}^2& \dots & (w_n^{n-n}{)}^n\end{array}\right]=\left[\begin{array}{cccccc}& n& 0& 0& \dots & 0\\ & 0& n& 0& \dots & 0\\ & 0& 0& n& \dots & 0\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 0& 0& 0& \dots & n\end{array}\right]$$

$$\therefore 矩阵的逆矩阵为：\frac{1}{n}\left[\begin{array}{cccccc}& 1& (w_n^{n-0}{)}^1& (w_n^{n-0}{)}^2& \dots & (w_n^{n-0}{)}^n\\ & 1& (w_n^{n-1}{)}^1& (w_n^{n-1}{)}^2& \dots & (w_n^{n-1}{)}^n\\ & 1& (w_n^{n-2}{)}^1& (w_n^{n-2}{)}^2& \dots & (w_n^{n-2}{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^{n-n}{)}^1& (w_n^{n-n}{)}^2& \dots & (w_n^{n-n}{)}^n\end{array}\right]$$

$$\therefore 系数\left[\begin{array}{cc}& a_0\\ & a_1\\ & a_2\\ & ⋮\\ & a_n\end{array}\right]=\frac{1}{n}\left[\begin{array}{cccccc}& 1& (w_n^{n-0}{)}^1& (w_n^{n-0}{)}^2& \dots & (w_n^{n-0}{)}^n\\ & 1& (w_n^{n-1}{)}^1& (w_n^{n-1}{)}^2& \dots & (w_n^{n-1}{)}^n\\ & 1& (w_n^{n-2}{)}^1& (w_n^{n-2}{)}^2& \dots & (w_n^{n-2}{)}^n\\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ & 1& (w_n^{n-n}{)}^1& (w_n^{n-n}{)}^2& \dots & (w_n^{n-n}{)}^n\end{array}\right]\times \left[\begin{array}{cc}& f(w_n^0)\\ & f(w_n^1)\\ & f(w_n^2)\\ & ⋮\\ & f(w_n^n)\end{array}\right]$$

最后我们注意到对于$w_n^{n-j}$，其实它为$w_n^j$的共轭复数；因此在IDFT中，只需将DFT中的乘以单位根变为乘以其共轭复数即可，但是分治结束后，需要除以$n$

其实上述需要求逆的矩阵是一个范德蒙德矩阵，也可以利用范德蒙德矩阵求逆的过程

好了！耗时三天终于敲完FFT的公式推导了！！！下面的代码实现就比较容易了！

4、FFT代码实现

直接上模版：

double PI = acos(-1.0);
//复数类
struct Complex
{
    double r, i;
    Complex() {}
    Complex(double r, double i) : r(r), i(i) {}
    inline void real(const double &x) { r = x; }
    inline double real() { return r; }
    inline Complex operator+(const Complex &rhs) const
    {
        return Complex(r + rhs.r, i + rhs.i);
    }
    inline Complex operator-(const Complex &rhs) const
    {
        return Complex(r - rhs.r, i - rhs.i);
    }
    inline Complex operator*(const Complex &rhs) const
    {
        return Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
    }
    inline Complex operator*(const double &rhs) const
    {
        return Complex(r * rhs, i * rhs);
    }
    inline void operator/=(const double &x)
    {
        r /= x, i /= x;
    }
    inline void operator*=(const Complex &rhs)
    {
        *this = Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
    }
    inline void operator+=(const Complex &rhs)
    {
        r += rhs.r, i += rhs.i;
    }
    inline Complex conj()
    {
        return Complex(r, -i);
    }
};
// i的二进制翻转数组
int bitrev[maxn];
// FFT预处理
void fft_prepare(int n)
{
    int bits = 0;
    //由于FFT需要项数为2的整数次方倍，所以多项式（形式）次数为第一个大于n的二的正整数次方，bits表示多少位
    while ((1 << bits) < n)
    {
        bits++;
    }
    for (int i = 0; i < n; i++)
    {
        //递推式：二进制翻转
        // 0 1 2 3 4 5 6 7
        // 0|4|2|6|1|5|3|7
        bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
    }
}

//快速傅里叶变换系数表达式转换为点值表达式
void fft(Complex *a, int n, int type = 1)
{
    //求出要迭代的序列
    for (int i = 0; i < n; i++)
    {
        if (i < bitrev[i])
        {
            swap(a[i], a[bitrev[i]]);
        }
    }
    //待合并区间的中点
    for (int mid = 1; mid < n; mid <<= 1)
    {
        //单位根
        Complex Wn{cos(PI / mid), type * sin(PI / mid)};
        // R是区间的右端点，j表示已经到哪个位置
        for (int R = mid << 1, j = 0; j < n; j += R)
        {
            Complex w{1, 0}; //幂
            //枚举左半部分
            for (int k = 0; k < mid; k++, w *= Wn)
            {
                //蝴蝶效应
                Complex x = a[j + k], y = w * a[j + mid + k];
                a[j + k] = x + y;
                a[j + k + mid] = x - y;
            }
        }
    }
    // type为-1则进行逆变换，即IDFT
    if (type == -1)
    {
        for (int i = 0; i < n; i++)
        {
            a[i].r /= n;
        }
    }
}

代码有注释，就不细说了！

5、FFT的瓶颈

虽然FFT的在时间复杂度上非常优秀，但是计算机在计算时仍然会遇到一些问题！

例如：由于复数类中使用的是double类型，而且在计算单位根存在除法，除此之外计算机计算下$\pi$始终是一个逼近值，因此必然也必然会存在误差！所以在特殊情况下计算机计算会出现一定的误差！

但是有没有解决办法喃？

可以说在一定情况下是可以避免误差的！这个时候就要引出原根

6、To NTT

原根：
$$\begin{gathered} 设m\in N^{\ast },a\in Z,若gcd(a,m)=1,且{\delta }_m(a)=\phi (m),则称a为模m的原根 \\ 注：{\delta }_m(a)为a模m的阶，\phi 为m的欧拉函数 \end{gathered}$$
通过定义我们可以知道，原根是一种定义在模意义下的定义，因此对于FFT遇到的的精度问题，我们也只能保证在模意义下的答案正确

我们通过推导可以发现原根具有和单位根几乎一样的性质，因此我们可以通过原根推出基本相同的结论

这样一来我们就不必考虑复数域了！！

也可以理解为NTT其实就是把FFT中的单位根换成了原根

6.1 NTT代码实现

const int g = 3, gi = 332748118, mod = 998244353;
// 998244353的一个原根为3且998244353-1=2^23*119，3在模998244353意义下的逆元为332748118
int qp(int a, int n)
{
    int ans = 1;
    while (n)
    {
        if (n & 1)
        {
            ans = ans * a % mod;
        }
        a = a * a % mod;
        n >>= 1;
    }
    return ans;
}
// i的二进制翻转数组
int bitrev[maxn];
// NTT预处理
void ntt_prepare(int n)
{
    int bits = 0;
    //由于NTT需要项数为2的整数次方倍，所以多项式（形式）次数为第一个大于n的二的正整数次方，bits表示多少位
    while ((1 << bits) < n)
    {
        bits++;
    }
    for (int i = 0; i < n; i++)
    {
        //递推式：二进制翻转
        // 0 1 2 3 4 5 6 7
        // 0|4|2|6|1|5|3|7
        bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
    }
}
// NTT，type=1时系数表示法转点值表示法，否则点值表示法转系数表示法
void NTT(int *a, int n, int type = 1)
{
    for (int i = 0; i < n; i++)
    {
        if (i < bitrev[i])
        {
            swap(a[i], a[bitrev[i]]);
        }
    }
    //在这之前NTT与FFT无异
    for (int mid = 1; mid < n; mid <<= 1)
    {
        //单位原根g_n
        int gn = qp(type == 1 ? g : gi, (mod - 1) / (mid << 1));
        //枚举具体区间，j也就是区间右端点
        for (int R = mid << 1, j = 0; j < n; j += R)
        {
            int g0 = 1;
            //合并，记得要取模
            for (int k = 0; k < mid; k++, g0 = gn * g0 % mod)
            {
                int x = a[j + k], y = g0 * a[mid + j + k] % mod;
                a[j + k] = (x + y) % mod;
                a[j + k + mid] = (x - y + mod) % mod;
            }
        }
    }
}

可以发现，NTT和FFT代码大体逻辑是相同的！！！只是从单位根换成了原根

6.2 tips

在抽象代数中，原根就是循环群的生成元。这个概念只在模 $m$缩剩余系关于乘法形成的群中有“原根”这个名字，在一般的循环群中都称作“生成元”。

并非每个模$m$缩剩余系关于乘法形成的群都是循环群，存在原根就表明它同构于循环群，如果不存在原根就表明不同构。