FFT快速傅里叶变换

DFT

把系数表示把转化为点值表示法。

有多项式$f_n(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$，其中$n$是一个$2$的幂，现在要快速求出$\forall i,0\le i<n$的$f({\omega }_n^i)$。

首先把$f_n(x)$的系数按照奇偶性分成两组，分别为$h_{\frac{n}{2}}(x)={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}{x}^i$和$g_{\frac{n}{2}}(x)={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{x}^i$

那么把$f_n({\omega }_n^i)$拆一下，就有如下：
$$\begin{array}{rl}{f}_n({\omega }_n^i)& =\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j}({\omega }_n^i{)}^{2j}+\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j+1}({\omega }_n^i{)}^{2j+1}\\ & =\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j}({\omega }_n^{2i}{)}^j+\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j+1}({\omega }_n^{2i}{)}^j{\omega }_n^i\\ & =\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j}({\omega }_{\frac{n}{2}}^i{)}^j+\sum _{j=0}^{\frac{n}{2}-1}{a}_{2j+1}({\omega }_{\frac{n}{2}}^i{)}^j{\omega }_n^i\\ & =h_{\frac{n}{2}}({\omega }_{\frac{n}{2}}^i)+{\omega }_n^i{g}_{\frac{n}{2}}({\omega }_{\frac{n}{2}}^i)\end{array}$$
发现$g,h$​和$f$​是类似的，那么可以分治处理。

然后发现位置的变换是蝴蝶变换，就可以非递归计算了。

蝴蝶变换

可以发现在分治过程中，位置$i$的值会唯一的变换到$i^{\prime }$的位置，那么可否快速计算，然后从低到高合并呢？这样就可以降低常数。

当然是可以的，$i^{\prime }$的二进制表达式和$i$的二进制表达式是前后对调的，证明就不写了。

IDFT

把点值化为系数表示法。

假设现在要做多项式乘法$C(x)=A(x)B(x)$，那么可以把他们都拓展到$(n-1)$次的多项式，其中$n$是$2$的幂。接着对$A(x),B(x)$做$DFT$，将其转化为点值表示法，那么对于带入的$n$个${\omega }_n^i$，都有$C({\omega }_n^i)=A({\omega }_n^i)B({\omega }_n^i)$，就是把对应位置的点值乘到一块就好了，然后对其做$IDFT$即可。

假设$C(x)={\sum }_{i=0}^{n-1}{a}_i{x}^i$，$y_i=C({\omega }_n^i)$，那么有如下的矩阵乘法：
$$\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]$$
只要两边都乘上这个方阵的逆矩阵就好了，然后我们发现它的逆矩阵就是对应每一项取倒数并再除以$n$，如下：
$$\frac{1}{n}\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n^{-1}& {\omega }_n^{-2}& {\omega }_n^{-3}& \cdots & {\omega }_n^{-(n-1)}\\ 1& {\omega }_n^{-2}& {\omega }_n^{-4}& {\omega }_n^{-6}& \cdots & {\omega }_n^{-2(n-1)}\\ 1& {\omega }_n^{-3}& {\omega }_n^{-6}& {\omega }_n^{-9}& \cdots & {\omega }_n^{-3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{-(n-1)}& {\omega }_n^{-2(n-1)}& {\omega }_n^{-3(n-1)}& \cdots & {\omega }_n^{-(n-1{)}^2}\end{array}\right]$$
那么现在来证明：

对于乘完之后的方阵，第$i$行第$i$列的值可任意这样计算：${\sum }_{k=0}^{n-1}{\omega }_n^i{\omega }_n^{-i}={\sum }_{k=0}^{n-1}{\omega }_n^0=n$

对于第$i$行$j$列，其中$i\ne j$，可以这样计算：${\sum }_{k=0}^{n-1}({\omega }_n^{i-j}{)}^k=\frac{({\omega }_n^{i-j}{)}^n-({\omega }_n^{i-j}{)}^0}{{\omega }_n^{i-j}-1}=0$

再乘上矩阵前面的$\frac{1}{n}$，那么逆矩阵就得证了，于是又是一个$DFT$。这次的系数是计算出来的点值，带入的是${\omega }_n^{-i}$，最后把$\frac{1}{n}$乘进去就好了。

#include<bits/stdc++.h>
using namespace std;
const int N=1<<22;
typedef complex<double> comp;
comp f[N],g[N],h[N];
int r[N];
int nf,ng,nh=1;
void FFT(comp *f,int rev){//rev=1->DFT  rev=-1->IDFT
	for (int i=0;i<nh;++i) if (r[i]<i) swap(f[i],f[r[i]]);
	for (int l=2;l<=nh;l<<=1){
		for (int i=0;i<nh;i+=l){
			comp step=comp(cos(2.0*M_PI/l),sin(2.0*M_PI/l)*rev),mul=comp(1,0);
			for (int j=0;j<(l>>1);++j){
				comp tmp=f[i+j+(l>>1)];
				f[i+j+(l>>1)]=f[i+j]-mul*tmp;
				f[i+j]=f[i+j]+mul*tmp;
				mul=mul*step;
			}
		}
	}
}
int main(){
	scanf("%d%d",&nf,&ng);
	while (nh-1<nf+ng)
		nh<<=1;
	for (int i=0;i<=nf;++i){double x;scanf("%lf",&x);f[i]=comp(x,0);}
	for (int i=0;i<=ng;++i){double x;scanf("%lf",&x);g[i]=comp(x,0);}
	r[0]=0; 
	for (int i=1;i<nh;++i) r[i]=(r[i>>1]>>1)|(i&1?(nh>>1):0);
	FFT(f,1);
	FFT(g,1);
	for (int i=0;i<nh;++i) h[i]=f[i]*g[i];
	FFT(h,-1);
	for (int i=0;i<=nf+ng;++i) printf("%d ",(int)round(h[i].real()/nh));
}