浅谈 FFT 和 NTT
Lazy_Labs
·
·
科技·工程
离散傅里叶变换 (DFT)
可以干一些频谱分析的事，有时间再补充（坑）
省略一些推导，反正可知我们要干的事情：
\text{DFT}(a)_i=\sum_{j=0}^{n-1} a_i\omega_n^{ij}
显然可以用矩阵形式表示：
\vec{\text{DFT}(a)}=F\vec{a}
其中：
F=
\begin{bmatrix}
\omega_n^{00} &\omega_n^{01} &\cdots&\omega_n^{0(n-1)}\\
\omega_n^{10} &\omega_n^{11} &\cdots&\omega_n^{1(n-1)}\\
\vdots &\vdots &\ddots&\vdots\\
\omega_n^{(n-1)0}&\omega_n^{(n-1)1}&\cdots&\omega_n^{(n-1)(n-1)}
\end{bmatrix}
DIT-FFT
实际上，上述过程暴力做是 $O(n^2)$ 的。$\text{FFT}$ 所做的事情只不过是将这个复杂度优化到和 $\text{FWT}$ 一样的 $O(n\log n)$ 罢了。
## 介绍：分治！
这里介绍的是一般大众所写的 $\text{FFT}$。
为了方便，一般 $n$ 取为 $2^k$。
我们令：
$$f_{k,i,\{a\}}=\sum_{j=0}^{n-1}\omega_n^{ij}a_j$$
此处 $n=2^k$。
本式子表示的含义即为将 $\omega_n^i$ 带入 $\{a\}$ 为系数的多项式所得到的结果。
容易得知 $\text{DFT}(a)_i=f_{k,i,\{a\}}$。
$$
\begin{align*}
f_{k,i,\{a\}}&=\sum_{j=0}^{n'-1}\omega_n^{i(2j)}a_{2j}+\omega_n^{i(2j+1)}a_{2j+1}\\
&=\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j}+\omega_n^i\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j+1}\\
&=f_{k-1,i,\{b\}}+\omega_n^if_{k-1,i,\{c\}}
\end{align*}
$$
且有：
$$
\begin{align*}
f_{k,i+n',\{a\}}&=\sum_{j=0}^{n'-1}\omega_n^{(i+n')(2j)}a_{2j}+\omega_n^{(i+n')(2j+1)}a_{2j+1}\\
&=\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j}+\omega_n^{i+n'}\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j+1}\\
&=\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j}-\omega_n^i\sum_{j=0}^{n'-1}\omega_{n'}^{ij}a_{2j+1}\\
&=f_{k-1,i,\{b\}}-\omega_n^if_{k-1,i,\{c\}}
\end{align*}
$$
此处 $n'=2^{k-1}$，同时 $b_i=a_{2i},c_i=a_{2i+1}$。
## 线代推导？启动！
我们成功地发现了 $f_{k,i,\{a\}}$ 和 $f_{k,i+n',\{a\}}$ 均可以用 $f_{k-1,i,\{b\}}$ 和 $f_{k-1,i,\{c\}}$ 的线性组合表示。
即：
$$
\begin{bmatrix}
f_{k,i,\{a\}}\\
f_{k,i+n',\{a\}}
\end{bmatrix}
=
\begin{bmatrix}
1&\omega_n^i\\
1&-\omega_n^i
\end{bmatrix}
\begin{bmatrix}
f_{k-1,i,\{b\}}\\
f_{k-1,i,\{c\}}
\end{bmatrix}
$$
### 整合起来：
假设我们将本矩阵补全（以 $n=8,k=3$ 为例）。
$$
\begin{bmatrix}
\text{DFT}(a)_0\\
\text{DFT}(a)_1\\
\text{DFT}(a)_2\\
\text{DFT}(a)_3\\
\text{DFT}(a)_4\\
\text{DFT}(a)_5\\
\text{DFT}(a)_6\\
\text{DFT}(a)_7\\
\end{bmatrix}
=
\begin{bmatrix}
f_{3,0,\{a\}}\\
f_{3,1,\{a\}}\\
f_{3,2,\{a\}}\\
f_{3,3,\{a\}}\\
f_{3,4,\{a\}}\\
f_{3,5,\{a\}}\\
f_{3,6,\{a\}}\\
f_{3,7,\{a\}}\\
\end{bmatrix}
=
\begin{bmatrix}
1&0&0&0&\omega_8^0&0&0&0\\
0&1&0&0&0&\omega_8^1&0&0\\
0&0&1&0&0&0&\omega_8^2&0\\
0&0&0&1&0&0&0&\omega_8^3\\
1&0&0&0&\omega_8^4&0&0&0\\
0&1&0&0&0&\omega_8^5&0&0\\
0&0&1&0&0&0&\omega_8^6&0\\
0&0&0&1&0&0&0&\omega_8^7\\
\end{bmatrix}
\begin{bmatrix}
f_{2,0,\{b\}}\\
f_{2,1,\{b\}}\\
f_{2,2,\{b\}}\\
f_{2,3,\{b\}}\\
f_{2,0,\{c\}}\\
f_{2,1,\{c\}}\\
f_{2,2,\{c\}}\\
f_{2,3,\{c\}}\\
\end{bmatrix}
$$
注：
1. 在这里选出任意相隔为 $4$ 的行和列的交（如 第二、五行 和 第二、五列 的交）即为上述矩阵，否则没有贡献。整合起来便是如此。
2. $\omega_n^5=-\omega_n^1$，其余同理，这里这样写是因为好看。
故我们发现**合起来写的操作和我们对于每个小矩阵进行操作是等价的**，我们下面均对于那个小矩阵进行讨论。
同时我们发现左侧的我们要求的式子被一个矩阵转化为了右侧要求的式子，容易发现上下各一半显然互不关联，且是左侧的一个子问题，可以递归计算（但是大家都知道显然可以不用递归计算，接下来马上讲到）。
### 递归下去：
继续对右侧进行子矩阵的乘法：
$$
\begin{bmatrix}
f_{2,0,\{b\}}\\
f_{2,1,\{b\}}\\
f_{2,2,\{b\}}\\
f_{2,3,\{b\}}\\
f_{2,0,\{c\}}\\
f_{2,1,\{c\}}\\
f_{2,2,\{c\}}\\
f_{2,3,\{c\}}\\
\end{bmatrix}
=
\begin{bmatrix}
F'_2&O\\
O&F'_2\\
\end{bmatrix}
\begin{bmatrix}
f_{1,0,\{bb\}}\\
f_{1,1,\{bb\}}\\
f_{1,2,\{bc\}}\\
f_{1,3,\{bc\}}\\
f_{1,0,\{cb\}}\\
f_{1,1,\{cb\}}\\
f_{1,2,\{cc\}}\\
f_{1,3,\{cc\}}\\
\end{bmatrix}
$$
其中：
$$
F'_2=
\begin{bmatrix}
1&0&\omega_4^0&0\\
0&1&0&\omega_4^1\\
1&0&\omega_4^2&0\\
0&1&0&\omega_4^3\\
\end{bmatrix}
$$
$O$ 为全 $0$ 矩阵，$\{?b\},\{?c\}$ 则是形如先前一般的变换，具体地：
$$
\{bb\}=\{b\}_{2i+0}=\{a_0,a_2,a_4,a_6\}_{2i+0}=\{a_0,a_4\}\\
\{bc\}=\{b\}_{2i+1}=\{a_0,a_2,a_4,a_6\}_{2i+1}=\{a_2,a_6\}\\
\{cb\}=\{c\}_{2i+0}=\{a_1,a_3,a_5,a_7\}_{2i+0}=\{a_1,a_5\}\\
\{cc\}=\{c\}_{2i+1}=\{a_1,a_3,a_5,a_7\}_{2i+1}=\{a_3,a_7\}\\
$$
同理，可以继续写下：
$$
\begin{bmatrix}
f_{1,0,\{bb\}}\\
f_{1,1,\{bb\}}\\
f_{1,0,\{bc\}}\\
f_{1,1,\{bc\}}\\
f_{1,0,\{cb\}}\\
f_{1,1,\{cb\}}\\
f_{1,0,\{cc\}}\\
f_{1,1,\{cc\}}\\
\end{bmatrix}
=
\begin{bmatrix}
F'_1&O&O&O\\
O&F'_1&O&O\\
O&O&F'_1&O\\
O&O&O&F'_1\\
\end{bmatrix}
\begin{bmatrix}
f_{0,0,\{bbb\}}\\
f_{0,0,\{bbc\}}\\
f_{0,0,\{bcb\}}\\
f_{0,0,\{bcc\}}\\
f_{0,0,\{cbb\}}\\
f_{0,0,\{cbc\}}\\
f_{0,0,\{ccb\}}\\
f_{0,0,\{ccc\}}\\
\end{bmatrix}
\\
F'_1=
\begin{bmatrix}
1&\omega_2^0\\
1&\omega_2^1\\
\end{bmatrix}
$$
其中右侧的向量根据类似上述的推导，显然等于 $\begin{bmatrix}a_0&a_4&a_2&a_6&a_1&a_5&a_3&a_7\end{bmatrix}^\top$。
其中这个 $\top$ 表示转置，$A^\top_{i,j}=A_{j,i}$。在这里，就是把行向量变为了列向量。
### 最终 & 蝴蝶变换
我们显然可以最终得到一系列等式，最终构成了下面的形式：
$$
\begin{bmatrix}
\text{DFT}(a)_0\\
\text{DFT}(a)_1\\
\text{DFT}(a)_2\\
\text{DFT}(a)_3\\
\text{DFT}(a)_4\\
\text{DFT}(a)_5\\
\text{DFT}(a)_6\\
\text{DFT}(a)_7\\
\end{bmatrix}
=
F_3F_2F_1
\begin{bmatrix}
a_0\\
a_4\\
a_2\\
a_6\\
a_1\\
a_5\\
a_3\\
a_7
\end{bmatrix}
$$
其中，$F_3$ 就是上面的那个大矩阵，$F_2,F_1$ 则是继续递归下去的产物。我们应该会发现，最终右侧的向量的每个数都应当为 $f_{0,0,\{?\}}$ 的形式，特别地，本处的 $\{?\}$ 显然是一个长为 $1$ 的“数组”，值为 $a_{0\sim8}$。通过 $a\to b,c$ 可以得到最终的式子形如右侧的形式，我们称这个向量为 $\vec{r}$。
更加普遍地，我们可以证明此时 $r_i$ 为 $a_{rev(i)}$，这里的 $rev(i)$ 表示 $i$ 在二进制下反转后的结果（如 $rev(3)=rev(011_2)=110_2=6$，则 $r_3=a_{rev(3)}=a_6$）。
我们称上述变换为“蝴蝶变换”。
这个证明并不是重点，因此我们不再加以过多赘述（如证明等可以去 `oiwiki` 等或者随便找一篇 $\text{FFT}$ 的题解应当都有写），我们应当要注意到的是：
$$\vec{r}=R\vec{a}$$
在 $n=8$ 时，$R$ 为：
$$
\begin{bmatrix}
1&0&0&0&0&0&0&0\\
0&0&0&0&1&0&0&0\\
0&0&1&0&0&0&0&0\\
0&0&0&0&0&0&1&0\\
0&1&0&0&0&0&0&0\\
0&0&0&0&0&1&0&0\\
0&0&0&1&0&0&0&0\\
0&0&0&0&0&0&0&1\\
\end{bmatrix}
$$
于是我们就可以将我们的一次正向 $\text{FFT}$ 写成矩阵的形式形如：
$$
\vec{\text{DFT}(a)}=(F_3F_2F_1R)\vec{a}
$$
的矩阵形式了，更为普遍地：
$$
\vec{\text{DFT}(a)}=(F_kF_{k-1}\cdots F_1R)\vec{a}
$$
中间括号所括起来的式子就恰好等于之前提到的 $F$（什么你忘了？快去最前面看看吧！此处为何相等呢？因为 $\vec{\text{DFT}(a)}=F\vec{a}$ 啊），我们的程序也是像这样从右至左依次执行这样的矩阵操作完成一次 $\text{FFT}$ 的。
这样，我们就用线性代数的方法生动形象地理解了大众所写的 $\text{FFT}$ 的全过程，这中间甚至只用到了最简单的矩阵乘法知识，相信大家都能看懂！
## 代码
这里我们贺一下我以前写的构思代码，说明一下上述原理：
```cpp
void FFT(comp A[])
{
//相当于矩阵 R
for(int i=0;i<N;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
//每次最外层循环分别为 F_1,F_2,...,F_K
for(int ns=1;ns<N;ns<<=1)
{
//当前循环乘的矩阵为 F_k,ns 即 n'=2^(k-1)
//单位根 \omega_n,你看不出来是因为单位根是 2Pi/n = Pi/ns
comp w(cos(Pi/ns),sin(Pi/ns));
//可以看出，矩阵 F_i 就需相隔 2^(i-1) 行/列。
for(int n=ns<<1,j=0;j<N;j+=n)
{
//n=2^k,j 为当前乘 F'_i 的位移
comp W(1,0);
for(int i=0;i<ns;i++,W=W*w)
{
//所乘小矩阵 F'_i 对应的变换，W=\omega_n^i
comp x=A[j+i],y=W*A[j+ns+i];
A[j+i]=x+y;A[j+ns+i]=x-y;
}
}
}
}
```
注释很详细了，同时也充分展现了为何我们将 $F_i$ 写成 $F'_i$ 构成的分块矩阵的样式。
# DIT-FFT 的逆变换
我们先不加证明的说出：
$$\vec{\text{IDFT}(a)}=F^{-1}\vec{a}$$
## 主流写法
一种合理的想法是将上述的 $\text{FFT}$ 过程每个矩阵都取逆，即：
$$
\vec{\text{DFT}(a)}=(F_kF_{k-1}\cdots F_1R)\vec{a}
$$
变为
$$
\vec{\text{IDFT}(a)}=(F_1^{-1}\cdots F_{k-1}^{-1}F_k^{-1}R^{-1})\vec{a}
$$
但是如果你真的尝试去取逆的话（还记得吗，**等价于**，所以只要对每个小矩阵取逆，将结果并入大矩阵即可。当然，在代码中，我们并不需要多一步“并入大矩阵”操作），你会发现结果异常得难看，多带了一个 $\frac{1}{2}$。
于是，神秘的巧合出现了：
我们考虑原矩阵的逆，通过细心观察，容易发现其逆矩阵等于原矩阵的共轭矩阵！（即每个元素均在复数意义下取共轭）即：
$$
F^{-1}=
\begin{bmatrix}
\omega_n^{00} &\omega_n^{01} &\cdots&\omega_n^{0(n-1)}\\
\omega_n^{10} &\omega_n^{11} &\cdots&\omega_n^{1(n-1)}\\
\vdots &\vdots &\ddots&\vdots\\
\omega_n^{(n-1)0}&\omega_n^{(n-1)1}&\cdots&\omega_n^{(n-1)(n-1)}
\end{bmatrix}^{-1}
=
\cfrac{1}{n}
\begin{bmatrix}
\omega_n^{-00} &\omega_n^{-01} &\cdots&\omega_n^{-0(n-1)}\\
\omega_n^{-10} &\omega_n^{-11} &\cdots&\omega_n^{-1(n-1)}\\
\vdots &\vdots &\ddots&\vdots\\
\omega_n^{-(n-1)0}&\omega_n^{-(n-1)1}&\cdots&\omega_n^{-(n-1)(n-1)}
\end{bmatrix}
$$
读者自证不难。
这意味着什么？我们发现：
$$\text{IDFT}(a)_i=\sum_{j=0}^{n-1} a_i\omega_n^{-ij}$$
然后如此如此，如同上方一同推导一般，容易发现仅有转移矩阵（即先前所说的“小矩阵”）不同，其余完全一致！
转移矩阵的变化如下：
$$
\begin{bmatrix}
1&\omega_n^i\\
1&-\omega_n^i
\end{bmatrix}
\rightarrow
\begin{bmatrix}
1&\omega_n^{-i}\\
1&-\omega_n^{-i}
\end{bmatrix}
$$
所以我们仅仅需要给上述代码加一个简单的改动，就能同时做到 $\text{FFT}$ 和 $\text{IFFT}$。当然，这也正是当前主流的，只加了蝴蝶变换优化的快速傅里叶变换写法：
```cpp
void FFT(int N,comp A[],int id)//id 为 1 表示正变换，-1 为逆变换
{
//相当于矩阵 R
for(int i=0;i<N;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
//每次最外层循环分别为 F_1,F_2,...,F_K
for(int ns=1;ns<N;ns<<=1)
{
//当前循环乘的矩阵为 F_k,ns 即 n'=2^(k-1)
//单位根 \omega_n,你看不出来是因为单位根是 2Pi/n = Pi/ns
comp w(cos(Pi/ns),id*sin(Pi/ns));//在这里，多加了正负的判断
//可以看出，矩阵 F_i 就需相隔 2^(i-1) 行/列。
for(int n=ns<<1,j=0;j<N;j+=n)
{
//n=2^k,j 为当前乘 F'_i 的位移
comp W(1,0);
for(int i=0;i<ns;i++,W=W*w)
{
//所乘小矩阵（转移矩阵）对应的变换，W=\omega_n^i
comp x=A[j+i],y=W*A[j+ns+i];
A[j+i]=x+y;A[j+ns+i]=x-y;
}
}
}
}
```
当然，不能忘了还有一个 $\dfrac{1}{n}$ 的常数，在函数外乘上即可。
实际上呢，还有一种连代码都不需要改的方法完成“逆变换”：
## 重要优化
读者不难发现：
$$
F^2=
\cfrac{1}{n}
\begin{bmatrix}
1 & 0 & 0 & \cdots & 0 & 0 \\
0 & 0 & 0 & \cdots & 0 & 1 \\
0 & 0 & 0 & \cdots & 1 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 1 & \cdots & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0
\end{bmatrix}
$$
而右侧矩阵（我们令其为 $P$）的平方为单位矩阵 $I$（对角线全 $1$），所以不难得到：
$$
F^{-1}=
\cfrac{1}{n}
P
F
$$
这个等式的寓意是：先乘一遍 $F$，再乘一遍 $P$，最后全部除以 $n$ 就等于做了一次逆变换。
转换为代码语言就是：先直接用原来的函数 $\text{FFT}$ 执行一遍后，将 $1\sim n-1$ 倒序，最后全部除以 $n$ 即可做完逆变换。
这样子，我们甚至只需要在结尾 `reverse` 一下就可以完成先前的操作！
### 新代码出炉~
```cpp
#include<bits/stdc++.h>
using namespace std;
inline int read()
{
int x(0),f(1);char c=getchar();
while(c>'9'||c<'0')f=c=='-'?-1:1,c=getchar();
while(c>='0'&&c<='9')x=x*10+c-48,c=getchar();
return f*x;
}
const double Pi=acos(-1.0);
const int N=2000007;
struct comp
{
double r,i;
comp(double R=0,double I=0){r=R,i=I;}
}a[N<<1],b[N<<1];
int n,m,rev[N<<1];
comp operator+(comp a,comp b){return comp(a.r+b.r,a.i+b.i);}
comp operator-(comp a,comp b){return comp(a.r-b.r,a.i-b.i);}
comp operator*(comp a,comp b){return comp(a.r*b.r-a.i*b.i,a.r*b.i+b.r*a.i);}
void FFT(int N,comp A[])
{
//相当于矩阵 R
for(int i=0;i<N;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
//每次最外层循环分别为 F_1,F_2,...,F_K
for(int ns=1;ns<N;ns<<=1)
{
//当前循环乘的矩阵为 F_k,ns 即 n'=2^(k-1)
//单位根 \omega_n,你看不出来是因为单位根是 2Pi/n = Pi/ns
comp w(cos(Pi/ns),sin(Pi/ns));
//可以看出，矩阵 F_i 就需相隔 2^(i-1) 行/列。
for(int n=ns<<1,j=0;j<N;j+=n)
{
//n=2^k,j 为当前乘 F'_i 的位移
comp W(1,0);
for(int i=0;i<ns;i++,W=W*w)
{
//所乘小矩阵（转移矩阵）对应的变换，W=\omega_n^i
comp x=A[j+i],y=W*A[j+ns+i];
A[j+i]=x+y;A[j+ns+i]=x-y;
}
}
}
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++)a[i]=read();for(int i=0;i<=m;i++)b[i]=read();
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)?n>>1:0);
FFT(n,a);FFT(n,b);
for(int i=0;i<n;i++)a[i]=a[i]*b[i];
FFT(n,a);reverse(a+1,a+n);//即乘上矩阵 P
for(int i=0;i<=m;i++)printf("%d ",int(a[i].r/n+0.5));//别忘了除以 n
return 0;
}
```
## 补充说明
上述问题中我们不加证明得给出结论：离散傅里叶逆变换相当于乘上其逆矩阵，下面我们来说明这件事的合理性：
我们考虑把先前的结论推倒重来，我们需要得到一个快速计算多项式乘法的方法。更加一般地，我们考虑循环卷积，即：
$$c_k=\sum\nolimits_{i+j \bmod n=k}a_ib_j$$
仿照 $\text{FWT}$，我们需要一个矩阵 $F$ 使得：
$$F\vec{a} \cdot F\vec{b}=F\vec{c}$$
和 $\text{FWT}$ 类似地推导，可以得到：$F_{k,i}F_{k,j}=F_{k,i+j\bmod n}$，解得：$F_{k,i}=F_{k,1}^i,F_{k,1}^n=1$。为了满秩（使得矩阵可逆），我们令 $F_{k,1}=\omega_n^k$，故得到傅里叶矩阵：
$$
F=
\begin{bmatrix}
\omega_n^{00} &\omega_n^{01} &\cdots&\omega_n^{0(n-1)}\\
\omega_n^{10} &\omega_n^{11} &\cdots&\omega_n^{1(n-1)}\\
\vdots &\vdots &\ddots&\vdots\\
\omega_n^{(n-1)0}&\omega_n^{(n-1)1}&\cdots&\omega_n^{(n-1)(n-1)}
\end{bmatrix}
$$
正和先前提到的定义一致！
所以，我们所需要的“逆变换”，则和 $\text{FWT}$ 一样，是 $F$ 的逆矩阵。
没学过 $\text{FWT}$ 也没事，只要知道 $F\vec{a} \cdot F\vec{b}=F\vec{c}$ 即可。
也就是说，我们选中傅里叶矩阵的原因是它本身就是做“循环卷积”这件事的。
所以说如果 $n$ 过于小，调用 $\text{FFT}$ 会返回循环卷积的原因也在此了。
# DIF-FFT
这是本篇文章的重头戏，是一种新型（存疑）的 $\text{FFT}$。
你是否还在因为蝴蝶变换太难写，蝴蝶变换太慢导致每次都得卡常而痛心疾首？
那就快来看看这里，这里有最新最热的 DIF-FFT，可以帮助您省去无用的蝴蝶变换。
## 优化蝴蝶变换？
蝴蝶变换太无用了！我们要想办法把它优化掉。
在优化它之前，我们要先来看看蝴蝶变换矩阵 $R$ 的一些性质：
$$R_{i,j}=[i=rev(j)]$$
由 $rev(rev(i))=i$ 可得，$R$ 是一个对称矩阵，即 $R=R^\top$，即 $R_{i,j}=R_{j,i}$。
当然，同样有 $F=F^\top$。
同时，由 $R$ 矩阵对应的操作为：交换 $i$ 和 $rev(i)$ 可知，作用两遍 $R$ 矩阵相当于没有作用。那么写成矩阵乘法则是 $R^2=I$（单位矩阵）。
我们再来看看基于 DIT-FFT 的式子：
$$\vec{DFT(a)}=(F_kF_{k-1}\cdots F_1R)\vec{a}$$
我们将 $F_kF_{k-1}\cdots F_1$ 称为 $F'$，则有 $F=F'R$。
### 插播一则转置原理：
$(AB)^\top=B^\top A^\top$，多个的置换也同理，比如：
$(ABC)^\top=C^\top B^\top A^\top$。
读者自证不难。
有了这一工具，消去蝴蝶变换的 $R$ 矩阵则是触手可及了！
### 正式推导！
$$F=F'R\rightarrow F=F^\top=R^\top (F')^\top$$
我们先考虑 $\text{IDFT}(\text{DFT}(a))=a$ 的过程：
$$
\vec{\text{DFT(a)}}=F\vec{a}=R^\top (F')^\top\vec{a}\\
\vec{\text{IDFT(a)}}=F^{-1}\vec{a}=\dfrac{1}{n}PF\vec{a}=\dfrac{1}{n}PF'R\vec{a}\\
\downarrow\\
\text{IDFT}(\text{DFT}(a))=\frac{1}{n}PF'RR^\top(F')^\top\vec{a}
$$
发现蝴蝶变换可以抵消！于是我们尝试让：
$$
\vec{\text{nDFT(a)}}=(F')^\top\vec{a}\\
\vec{\text{nIDFT(a)}}=\dfrac{1}{n}PF'\vec{a}
$$
### 证明可行
通过观察发现，$(F')^\top$ 即为正常进行傅里叶变换之后再进行一次蝴蝶变换（即 $RF$）。此时点乘后依然是对应位置相乘。所以其实蝴蝶变换并没有什么软用，于是就可以直接丢掉了！
### 具体实现
我们成功将蝴蝶变换丢掉了，但是具体实现呢？
对于 $\text{nIDFT}$，我们发现它仅仅是丢掉了蝴蝶变换后的 DIT-FFT（即去掉蝴蝶变换后的原函数）后 `reverse` 再除以 $n$。和先前的代码并没有太大的差别（除了少了讨人厌的蝴蝶变换）。
然而对于 $\text{nDFT}$ 就不一样了，$(F')^\top$ 仍然是我们所不知道的一个东西，肯定得新写点东西来实现它。但是无妨，毕竟少了三次甚至两次蝴蝶变换，多写一点也无所谓。
### nDFT
我们继续分解 $F'$，由很早前我们所推得的结论可知：
$$F'=F_kF_{k-1}\cdots F_1$$
那就容易了，我们再次运用转置原理，得到：
$$(F')^\top=F_1^\top\cdots F_{k-1}^\top F_k^\top$$
显然可以发现，乘矩阵的方向反了，同时矩阵变为了其转置。
别忘了我们所说的 **等价于**，所以我们只需要考虑转移矩阵（小矩阵）的置换，就能推知组合为大矩阵了（当然，代码中并不需要我们实现“组合”的过程）。
再看一眼我们的转移矩阵：
$$
\begin{bmatrix}
1&\omega_n^i\\
1&-\omega_n^i
\end{bmatrix}^\top
=
\begin{bmatrix}
1&1\\
\omega_n^i&-\omega_n^i
\end{bmatrix}
$$
即此时的转移发生变化。
于是，我们就完成了 $\text{nDFT}$ 的具体实现。
总结一下：
1. 在转移时，需要 $F_k$ 先乘，从大往小乘。
2. 在转移时，转移矩阵改变，需要更改。
至此，DIF-FFT 的理论实现就讲完了。
接下来是代码时间：
### new-新代码出炉~
```cpp
#include<bits/stdc++.h>
using namespace std;
inline int read()
{
int x(0),f(1);char c=getchar();
while(c>'9'||c<'0')f=c=='-'?-1:1,c=getchar();
while(c>='0'&&c<='9')x=x*10+c-48,c=getchar();
return f*x;
}
const double Pi=acos(-1.0);
const int N=2000007;
struct comp
{
double r,i;
comp(double R=0,double I=0){r=R,i=I;}
}a[N<<1],b[N<<1];
int n,m,rev[N<<1];
comp operator+(comp a,comp b){return comp(a.r+b.r,a.i+b.i);}
comp operator-(comp a,comp b){return comp(a.r-b.r,a.i-b.i);}
comp operator*(comp a,comp b){return comp(a.r*b.r-a.i*b.i,a.r*b.i+b.r*a.i);}
void FFT(int N,comp A[])
{
//同样没有蝴蝶变换
//每次最外层循环分别为 F_k^\top,...,F_2^\top,F_1^\top
for(int ns=N>>1;ns;ns>>=1)
{
//当前循环乘的矩阵为 F_k^\top,ns 即 n'=2^(k-1)
//单位根 \omega_n,你看不出来是因为单位根是 2Pi/n = Pi/ns
comp w(cos(Pi/ns),sin(Pi/ns));
//可以看出，矩阵 F_i^\top 就需相隔 2^(i-1) 行/列。
for(int n=ns<<1,j=0;j<N;j+=n)
{
//n=2^k,j 为当前乘 F'_i^\top 的位移
comp W(1,0);
for(int i=0;i<ns;i++,W=W*w)
{
//所乘小矩阵（转移矩阵）对应的变换，W=\omega_n^i
comp x=A[j+i]+A[j+ns+i],y=A[j+i]-A[j+ns+i];
A[j+i]=x;A[j+ns+i]=W*y;
}
}
}
}
void IFFT(int N,comp A[])
{
//删去了蝴蝶变换
// for(int i=0;i<N;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
//每次最外层循环分别为 F_1,F_2,...,F_K
for(int ns=1;ns<N;ns<<=1)
{
//当前循环乘的矩阵为 F_k,ns 即 n'=2^(k-1)
//单位根 \omega_n,你看不出来是因为单位根是 2Pi/n = Pi/ns
comp w(cos(Pi/ns),sin(Pi/ns));
//可以看出，矩阵 F_i 就需相隔 2^(i-1) 行/列。
for(int n=ns<<1,j=0;j<N;j+=n)
{
//n=2^k,j 为当前乘 F'_i 的位移
comp W(1,0);
for(int i=0;i<ns;i++,W=W*w)
{
//所乘小矩阵（F'_i）对应的变换，W=\omega_n^i
comp x=A[j+i],y=W*A[j+ns+i];
A[j+i]=x+y;A[j+ns+i]=x-y;
}
}
}
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++)a[i]=read();for(int i=0;i<=m;i++)b[i]=read();
for(m+=n,n=1;n<=m;n<<=1);
FFT(n,a);FFT(n,b);
for(int i=0;i<n;i++)a[i]=a[i]*b[i];
IFFT(n,a);reverse(a+1,a+n);//即乘上矩阵 P
for(int i=0;i<=m;i++)printf("%d ",int(a[i].r/n+0.5));//别忘了除以 n
return 0;
}
```
实战上至少比加了蝴蝶优化的快非常多，但是还是很慢！毕竟有复数，同时还怕精度丢失。
这样，就不得不介绍：NTT 了。
# 基于 DIF-FFT 的朴素 NTT
我们发现模意义下的原根也有很类似单位根的性质，于是：
## 原根
我们将复数的 单位根 换为，模一个大质数 $p$ 意义下的“原根”。但是同样需要满足一些性质：
$$g_n^n=1,g_n^{n'}=-1,g_n^{2i}=g_{n'}^i$$
因为我们已经选取了 $n=2^k$，故对于 $p=r\times 2^q+1$ 的质数来说，需要找到一个数 $g$ 使得其阶为 $2^k$ 即可（即 $g^{2^k}=1$）。
容易发现，原阶为 $r\times 2^q$，故我们取 $\omega_n=g^{r2^{q-k}}=g^\frac{p-1}{n}$ 即可（此处需满足 $n\le 2^q$ 即 $k \le q$）。
容易发现这样选取可以轻松满足先前的要求。
## 经典模数
$998244353=119\times 2^{23}+1,g=3$，可以做到 $n\le 8388608$。一般选取这个即可。
$167772161=5\times 2^{25}+1,g=3$，可以做到 $n\le 33554432$。基本能应对大部分题目。
其它模数在参考文章中有提到。
## 具体实现
只需将复数的加减乘除换为模意义下的加减乘除，单位根换为原根的幂次即可。
## 常数优化
1. 最后除以的 $\dfrac{1}{n}$ 不需要用快速幂计算。我们发现直接取 $p-\dfrac{p-1}{n}$ 即可，这必然能被整除。
2. 可以将我们所用到的所有“单位根” $\omega_n^i$ 预处理到数组 `w[ns+i]` 的位置上（因为 $i\in [0,n')$），这样我们访问单位根便是取址连续的，有优秀的常数。
3. 处于记忆，在模数 $p=998244353$ 时，我们不取 $3^{119}\bmod p$，而取 $31$。容易验证，其也满足 $31^{2^{23}}\equiv 1\pmod p$。那么 $\omega_n=31^{2^{23}/n}$。
## 代码
```cpp
#include<bits/stdc++.h>
using namespace std;
inline int read()
{
int x(0),f(1);char c=getchar();
while(c>'9'||c<'0')f=c=='-'?-1:1,c=getchar();
while(c>='0'&&c<='9')x=x*10+c-48,c=getchar();
return f*x;
}
const double Pi=acos(-1.0);
const int N=2000007,mod=998244353;
int w[N<<1];
void InitOmega(int N)
{
static int Bs[23];long long rt=31;
for(int i=22;i>=0;i--)Bs[i]=rt,rt=rt*rt%mod;
for(int n=1,tot=0;n<N;n<<=1)
{
long long nw=1;const int B=Bs[tot++];
for(int i=0;i<n;i++)w[n+i]=nw,nw=nw*B%mod;
}
}
int n,m,a[N<<1],b[N<<1];
void NTT(int N,int A[])
{
for(int ns=N>>1;ns;ns>>=1)for(int j=0;j<N;j+=(ns<<1))for(int i=0;i<ns;i++)
{
//w[ns+i] 即为 \omega_n^i
const int x=A[j+i]+A[j+ns+i],y=A[j+i]-A[j+ns+i]+mod;
A[j+i]=(x>=mod?x-mod:x);A[j+ns+i]=1ll*w[ns+i]*y%mod;
}
}
void INTT(int N,int A[])
{
for(int ns=1;ns<N;ns<<=1)for(int j=0;j<N;j+=(ns<<1))for(int i=0;i<ns;i++)
{
//w[ns+i] 即为 \omega_n^i
const int x=A[j+i],y=1ll*w[ns+i]*A[j+ns+i]%mod;
A[j+i]=(x+y>=mod?x+y-mod:x+y);A[j+ns+i]=(x-y<0?x-y+mod:x-y);
}
reverse(a+1,a+N);const long long inv=mod-(mod-1)/N;
for(int i=0;i<N;i++)a[i]=a[i]*inv%mod;
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++)a[i]=read();
for(int i=0;i<=m;i++)b[i]=read();
for(m+=n,n=1;n<=m;n<<=1);InitOmega(n);
NTT(n,a);NTT(n,b);for(int i=0;i<n;i++)a[i]=1ll*a[i]*b[i]%mod;
INTT(n,a);for(int i=0;i<=m;i++)printf("%d ",a[i]);
return 0;
}
```
还有一个优化，之后再说，咕咕咕。
# 参考文章
1. [快速傅里叶变换](https://www.cnblogs.com/JerryTcl/p/17155756.html) by [JerryTcl](https://www.luogu.com.cn/user/370352)
1. [FFT/FWT 相关理论推导复习](https://www.cnblogs.com/HaHeHyt/p/18190014) by HaHeHyt
1. [快速傅里叶变换](https://oi-wiki.org/math/poly/fft/) by [OIwiki](https://oi-wiki.org/)
1. [素数原根表](https://www.cnblogs.com/LzyRapx/p/9109496.html) by LzyRapx
1. [常用素数原根表（NTT快速数论变换）](https://blog.csdn.net/weixin_45697774/article/details/115732914) by 繁凡さん