FFT是个啥？
简单来说就是一个计算多项式乘法的东西呀..
以下内容基本都是在大黑书《算法导论》上的..
总述
对于项数为$n$的多项式$A(x)$和项数为$m$的多项式$B(x)$，可以如此表达：
$$A(x)=A_0+A_1x+A_2x^2+A_3x^3+...+A_{n-1}x^{n-1}$$
$$B(x)=B_0+B_1x+B_2x^2+B_3x^3+...+B_{m-1}x^{m-1}$$
把这两个多项式相乘可得到项数为$(n+m-1)$的多项式$C(x)$：
其中对于其任意系数$C_i$有$C_i=\sum\limits_{j=0}^{i}A_j\cdot B_{i-j}$
然后的话具体步骤如下：
系数表达->点值表达->相乘->点值表达->系数表达
把系数表达式转为点值表达式以及其逆运算都是$O(nlogn)$的
相乘不用说，$O(n)$
总复杂度为$O(nlogn)$
点值表达式
简单来说就是把多项式$A(x)$代入若干个$x$的值得到若干个点
比如说$A(x)$的点值表达式为${(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),...,(x_{n-1},A(x_{n-1}))}$
我们把从点值表达式转化为系数表达式的操作称为插值
插值多项式的唯一性定理
对于有$n$个点组成的集合仅存在一个项数为$n$的多项式与之对应
FFT即为取若干个特殊的值进行点值表达
n次单位复数根
从这里开始，所有提到的$n$将满足$(n=2^k,k\in N)$
n次单位复数根：即$\omega$使得$\omega^n=1$
复数与$e$的关系：$$e^{iu}=\cos u+i\sin u$$
（上述公式由于复数运算性质与指数运算性质大致相同，可以用泰勒公式证明，这里挖一个坑）
由于$e^{2\pi i}=\cos 2\pi+i\sin 2\pi=1$
所以定义主n次单位根为$\omega_n=e^{2\pi i/n}$
那么从上面可得 $\omega_n^i\ (i=0,1,2,...,n-1)$在平面直角坐标系上均匀分布且模长为$1$
也就有$\omega_n^i\cdot\omega_n^j=\omega_n^{(i+j)\ mod\ n}$
一些定理
消去引理： $$\omega_{dn}^{dk}=\omega_n^k$$
证明$$\omega_{dn}^{dk}=e^{2\pi i\cdot dk/dn}=e^{2\pi i\cdot k/n}=\omega_n^k$$
折半引理：如果$n>0$且$n$为偶数，那么$n$个n次单位复数根的平方的集合就是$n/2$个n/2次单位复数根的集合
证明
$$(\omega_n^k)^2=\omega_{n/2}^k$$
$$(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\cdot\omega_n^n=\omega_{n/2}^k$$
由此我们可以得到$$(\omega_n^k)^2=(\omega_n^{k+n/2})^2$$
即$$\omega_n^k=-\omega_n^{k+n/2}$$
这将会很有用
求和引理：对于$n\in Z,n\geq 1$，且$n\nmid k$，有$$\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0$$
由等比数列求和公式可得$$\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=\dfrac{(\omega_n^k)^n-1}{\omega_n^k-1}=0$$
DFT
我们计算多项式$A(x)=\sum\limits_{i=0}^{n-1}a_ix^i$在$\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}$处的结果
也就是其点值表达式$(y_0,y_1,y_2,...,y_{n-1})$
称其为$(a_0,a_1,a_2,...,a_{n-1})$的离散傅里叶变换
$$y=DFT_n(a)$$
FFT
利用复数单位根的特殊性质即可快速求出（即在$O(nlogn)$时间内）求出$DFT_n(a)$
这是一种分治策略，先把偶数系数和奇数系数的分开
即
$$A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}$$
$$A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}$$
所以可以得到
$$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$
假设通过递归分治前两者已经求好
对于$k\in [0,n/2-1]$
$$y_k=y^{[0]}_k+\omega_n^ky^{[1]}_k$$
$$y_{k+n/2}=y^{[0]}_k-\omega_n^ky^{[1]}_k$$
证明：
$y_k=y^{[0]}_k+\omega_n^ky^{[1]}_k$
$=A^{[0]}(\omega_{n/2}^k)+\omega_n^kA^{[1]}(\omega_{n/2}^k)$
$=A^{[0]}((\omega_n^k)^2)+\omega_n^kA^{[1]}((\omega_n^k)^2)$
$=A(\omega_n^k)$
同样的
$y_{k+n/2}=y^{[0]}_k-\omega_n^ky^{[1]}_k$
$=A^{[0]}(\omega_{n/2}^k)+\omega_n^{k+n/2}A^{[1]}(\omega_{n/2}^k)$
$=A^{[0]}((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}A^{[1]}((\omega_n^{k+n/2})^2)$
$=A(\omega_n^{k+n/2})$
证明结束，我们此时把$\omega_n^k$称作旋转因子
插值
由前面的式子，可以把$DFT$写成矩阵乘积的形式$y=V_na$
其中$V_n$是由$\omega_n$的若干次幂填充的矩阵
$$\left[\begin{array}{aaa}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1}\end{array}\right]=\left[\begin{array}{ccc}1 & 1 & 1 & 1 \ldots & 1 \\1 & \omega_n & \omega_n^2 & \omega_n^3 & \ldots & \omega_n^{n-1} \\1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \ldots & \omega_n^{2(n-1)} \\1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \ldots & \omega_n^{3(n-1)} \\\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \ldots & \omega_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{ppp}a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ y_{n-1}\end{array}\right]$$
对于$j,k\in [0,n-1]$，$V_n$的$(j,k)$元素为$\omega_n^{kj}$
所以$$a=y\cdot V_n^{-1}$$
定理：对于$j,k\in [0,n-1]$，$V_n^{-1}$的$(j,k)$元素为$\omega_n^{-kj}/n$
证明：我们只需证明$V_n^{-1}V_n=I_n$，即$n\times n$的单位矩阵
对于该矩阵$(j,j')$的元素
$$\left[V_n^{-1}V_n\right]_jj'=\sum\limits_{k=0}^{n-1}(\omega_n^{-kj}/n)(\omega_n^{kj'})=\sum\limits_{k=0}^{n-1}\omega_n^{k(j'-j)}/n$$
那么如果$j=j'$此和为$1$，否则根据求和引理，此和为$0$
证明完成
由此可得$DFT_n^{-1}(y)$
$$a_j=\dfrac{1}{n}\sum\limits_{k=0}^{n-1}y_k\cdot \omega_n^{-kj}$$
高效实现FFT
也就是使用迭代代替递归
比如$n=8$，以如下排列顺序即可实现
$${a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7}$$
也称作位逆序置换
Code
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cmath>
using namespace std;
const int Maxn = 400010;
struct cp {
double r, i;
cp (double r = 0.0, double i = 0.0) : r(r), i(i) {}
}A[Maxn], B[Maxn]; int An, Bn, N;
cp operator +(cp a, cp b) { return cp(a.r+b.r, a.i+b.i); }
cp operator -(cp a, cp b) { return cp(a.r-b.r, a.i-b.i); }
cp operator *(cp a, cp b) { return cp(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r); }
const double PI = acos(-1);
void fft(cp *a, int op) {
int i, j, k;
j = 0;
for(i = 0; i < N; i++){
if(i < j) swap(a[i], a[j]);
k = N >> 1;
while(j & k){ j -= k; k >>= 1; }
j += k;
}
for(i = 2; i <= N; i <<= 1){
cp wn = cp(cos(2.0*PI/i), op*sin(2.0*PI/i));
for(j = 0; j < N; j += i){
cp w = cp(1, 0);
for(k = j; k < j+i/2; k++){
cp x = a[k], y = w*a[k+i/2];
a[k] = x+y; a[k+i/2] = x-y;
w = w*wn;
}
}
}
if(op == -1) for(i = 0; i < N; i++) A[i].r /= N;
}
int main() {
int i, j, k;
scanf("%d%d", &An, &Bn);
N = 1;
while(N-1 <= An+Bn) N <<= 1;
for(i = 0; i <= An; i++) scanf("%lf", &A[i].r);
for(i = 0; i <= Bn; i++) scanf("%lf", &B[i].r);
fft(A, 1); fft(B, 1);
for(i = 0; i < N; i++) A[i] = A[i]*B[i];
fft(A, -1);
for(i = 0; i <= An+Bn; i++) printf("%d ", (int)(A[i].r+0.5));
printf("\n");
return 0;
}
完结撒花..
拖了差不多半年的坑终于补上了..
