简单来说就是一个计算多项式乘法的东西呀..
以下内容基本都是在大黑书《算法导论》上的..
总述
对于项数为$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=sumlimits_{j=0}^{i}A_jcdot 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,kin N)$
n次单位复数根：即$omega$使得$omega^n=1$
复数与$e$的关系：$$e^{iu}=cos u+isin u$$
（上述公式由于复数运算性质与指数运算性质大致相同，可以用泰勒公式证明，这里挖一个坑）
由于$e^{2pi i}=cos 2pi+isin 2pi=1$
所以定义主n次单位根为$omega_n=e^{2pi i/n}$
那么从上面可得 $omega_n^i (i=0,1,2,...,n-1)$在平面直角坐标系上均匀分布且模长为$1$
也就有$omega_n^icdotomega_n^j=omega_n^{(i+j) mod n}$
一些定理
消去引理： $$omega_{dn}^{dk}=omega_n^k$$
证明$$omega_{dn}^{dk}=e^{2pi icdot dk/dn}=e^{2pi icdot 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}cdotomega_n^n=omega_{n/2}^k$$
由此我们可以得到$$(omega_n^k)^2=(omega_n^{k+n/2})^2$$
即$$omega_n^k=-omega_n^{k+n/2}$$
这将会很有用
求和引理：对于$nin Z,ngeq 1$，且$n mid k$，有$$sumlimits_{j=0}^{n-1}(omega_n^k)^j=0$$
由等比数列求和公式可得$$sumlimits_{j=0}^{n-1}(omega_n^k)^j=dfrac{(omega_n^k)^n-1}{omega_n^k-1}=0$$
DFT
我们计算多项式$A(x)=sumlimits_{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)$$
假设通过递归分治前两者已经求好
对于$kin [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[egin{array}{aaa}y_0 \ y_1 \ y_2 \ y_3 \ vdots \ y_{n-1}end{array} ight]=left[egin{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} ight]left[egin{array}{ppp}a_0 \ a_1 \ a_2 \ a_3 \ vdots \ y_{n-1}end{array} ight]$$
对于$j,kin [0,n-1]$，$V_n$的$(j,k)$元素为$omega_n^{kj}$
所以$$a=ycdot V_n^{-1}$$
定理：对于$j,kin [0,n-1]$，$V_n^{-1}$的$(j,k)$元素为$omega_n^{-kj}/n$
证明：我们只需证明$V_n^{-1}V_n=I_n$，即$n imes n$的单位矩阵
对于该矩阵$(j,j')$的元素
$$left[V_n^{-1}V_n ight]_jj'=sumlimits_{k=0}^{n-1}(omega_n^{-kj}/n)(omega_n^{kj'})=sumlimits_{k=0}^{n-1}omega_n^{k(j'-j)}/n$$
那么如果$j=j'$此和为$1$，否则根据求和引理，此和为$0$
证明完成
由此可得$DFT_n^{-1}(y)$
$$a_j=dfrac{1}{n}sumlimits_{k=0}^{n-1}y_kcdot 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("
");
return 0;
}
完结撒花..
拖了差不多半年的坑终于补上了..