快速傅里叶变换(FFT)——从入门到入土

FFT (快速傅里叶变换)——从入门到入土

请不要转载，谢谢。版权于江苏省前黄高级中学 刘成宇 手中

前置知识

多项式系数/点值表达法

对于一个多项式 $A(x)=\sum\limits_{i=0}^{n-1}{a_i*x^i}$

明显是一个 $n-1$ 次函数/多项式，给这种对于 对应关系 $A$ 的种类一个名称——多项式的系数表达法。

还有什么表示方式呢？下面给出一条引理

用 $n+1$ 个点可以确定一条 $n$ 次函数图像。

感性证明：

$\begin{cases}y_1=a_0+a_1*x_1+...+a_{n-1}*x_1^{n-1} \\ y_2=a_0+a_1*x_2+...+a_{n-1}*x_2^{n-1} \\ y_3=a_0+a_1*x_3+...+a_{n-1}*x_3^{n-1} \\ ...\end{cases}$

这时候注意到，如果我现在 $x_1 \not= x_2 \not= x_3 ... \not= x_n+1$ 那么我是不是用加减消元法最终可以确定 $a_0,a_1,a_2,...,a_{n-1}$

那么这时候就简单了，可以用点刻画函数。即为点值表示法。用形如 $A(x)=\{(x_1,y_1),(x_2,y_2),...,(x_{n},y_{n})\}$ 表示

复数计算

如果两个复数相乘 $(a+bi)*(c+di)=ac-bd+(bc+ad)i$

如果两个复数相加或相减比较简单，这里就不给出了。

复平面

本文令 $i = \sqrt{-1}$

快速傅里叶变换是基于复数的，复平面是其中一个基本的知识。

注意，原点是不在虚轴上的，因为 $0i=0$ 而 $0$ 是实数

就类似于上面这样一个平面。对于平面上任意一个点 $(a,b)$ 其都表示一个复数 $a+bi$ 。

如果我现在要对复平面上的一个点 $(a,b)$ 平方。

因为在单位圆上 所以有 $|(a,b)|=\sqrt{a^2+b^2}=1$ 即 $a^2+b^2=1$

首先计算它的值 $(a+bi)^2=a^2-b^2+2abi$

新的点 $(a^2-b^2,2ab)$ 计算一下模长

$|(a^2-b^2,2ab)|=\sqrt{(a^2-b^2)^2+4a^2b^2}=\sqrt{a^4+b^4-2a^2b^2+4a^2b^2}=\sqrt{a^4+b^4+2a^2b^2}=\sqrt{(a^2+b^2)^2}=1$

所以，现在的这个点仍然在单位圆上。

接下来考虑一下它与原点连线后的旋转角 $\theta_2$ 与原旋转角 $\theta_1$ 的关系。

由于原来的点坐标是 $(a,b)$ 所以有 $\tan{\theta_1}=\frac{b}{a}$ 不妨令 $\frac{b}{a}=\lambda$

那么综合现在的坐标就有 $\tan{\theta_2}=\frac{2ab}{a^2-b^2}$ 类似于 whk 中处理的套路，齐次式将上边的除下去。

$\tan{\theta_2}=\frac{2}{\frac{a}{b}-\frac{b}{a}}=\frac{2\lambda}{1-\lambda^2}$ 不难看出二倍角公式，那么有

$\theta_2=2\theta_1$ 就相当于多转了一个 $\theta_1$ 这种特殊性就是 FFT 为什么所撷取的。

真正的FFT

FFT 正变换

快速傅里叶变换，处理多项式乘法时非常有用。既然是多项式乘法，那不可能是用一般的系数表达式一个一个乘，那么得想一个更好的办法，由于:

对于点值表达法计算 $A(i)\times B(i)$ 答案 $C(i)=\{(x_0,A_{x0}B_{x0}),(x_1,A_{x1}B_{x1}),(x_2,A_{x2}B_{x2}),...,(x_n,A_nB_n)\}$

明显发现对于点值表达法计算乘法只需要 $\mathbf{O}(n)$ 如此优秀的时间就 ok 了，那么直接使用点值表达法就 ok 啦。这时候遇到的最大的一个问题就是，对于直观的系数表达法与抽象但好用的点值表达法之间，我如何去转化，这一点是值得我们去思考的。这也是 FFT 区别于或者说优于 DFT 的一个特点。

现在就是要计算 $n$ 个 $A(x)$ 且这些 $x$ 互不相同，总时间复杂度还不能有 $\mathbf{O}(n^2)$ 怎么办。

先把 $n$ 转化成一个大于 $n$ 二的整数次幂的形式，不够的高位上系数补 $0$ ，然后进行操作。

复数就派上用场了，选取的就是 $\omega^k_n\quad \forall k\in[0,n)$ 这个时候因为这里的 $n$ 肯定是大于原来的 $n$ 的，所以即使取 $n$ 个点也是可以计算出原来的值的。

现在就是要计算 $A(\omega^k_n)$ 了。

首先，稍微的将 $A(x)=\sum\limits_{i=0}^{n-1}{a_i*x^i}$ 展开看看 $A(x)=a_0+a_1*x+a_2*x^2+...+a_{n-1}x^{n-1}$ 这时候做这样一步构造

令 $A_1(x)=a_0+a_2*x+a_4*x^2+...+a_{n-2}*x^{\frac{n}{2}-1}$
$A_2(x)=a_1+a_3*x+a_5*x^2+...+a_{n-1}*x^{\frac{n}{2}-1}$

显然有 $A(x)=A_1(x^2)+x*A_2(x^2)$
将 $A(\omega^k_n)$ 写出来看看

$A(\omega^k_n)=A_1(\omega^{2k}_n)+\omega^k_n*A_2(\omega^{2k}_n)=A_1(\omega^k_{\frac{n}{2}})+\omega^k_n*A_2(\omega^k_{\frac{n}{2}})$

好的，但是问题并没有解决，即使化成这种形式，还是找不到任何优化的点。这时候分类讨论一波

1. 如果当 $0\leq k<\frac{n}{2}$ 时，我就用$A(\omega^k_n)=A_1(\omega^k_{\frac{n}{2}})+\omega^k_n*A_2(\omega^k_{\frac{n}{2}})$算就 ok 了。
2. 如果当 $n>k\ge\frac{n}{2}$ 时，不妨令 $t=k+\frac{n}{2}$ 那么这个 $t$ 又回到了上一个情况了，就有 $A(k)=A_1(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_n*A_2(\omega^{2k+n}_n)$

首先 $\omega^{2k+n}_n=\omega^{2k}_n*\omega^n_n=\omega^{2k}_n$
并且 $\omega^{k+\frac{n}{2}}_n=-\omega^{k}_n$
那么就有 $A(k)=A_1(\omega^{2k}_n)-\omega^{k}_n*A_2(\omega^{2k}_n)=A_1(\omega^k_{\frac{n}{2}})-\omega^k_n*A_2(\omega^k_{\frac{n}{2}})$

连起来

$A(\omega^k_n)=A_1(\omega^k_{\frac{n}{2}})+\omega^k_n*A_2(\omega^k_{\frac{n}{2}})\qquad k\in[0,\frac{n}{2})$
$A(\omega^{k+\frac{n}{2}}_n)=A_1(\omega^k_{\frac{n}{2}})-\omega^k_n*A_2(\omega^k_{\frac{n}{2}})\qquad k\in[\frac{n}{2},n)$

发现了特殊的地方，两个的差距只有一个 $\omega^k_n*A_2(\omega^k_{\frac{n}{2}})$ ,那么就其实我们可以通过计算 $A_1(\omega^k_{\frac{n}{2}})$ 以及 $A_2(\omega^k_{\frac{n}{2}})*\omega^k_n$ 就能求出 $A(\omega^k_n)$ 以及 $A(\omega^{k+\frac{n}{2}}_n)$

那么求 $A(1),A(2),...,A(n)$ 就被切成了一半，只需要求一半就 ok 了。

这个时候观察一下所要求的 $A_1(x),A_2(x)$ 也是一个可以使用这种方法的式子，那么使用这种同样的方法就可以了。

有点像线段树分治一样。

FFT 负变换

上面的正变换只是一半，是把这个东西转化成点值表达法的方式，由于 正常 的题目最后肯定还是要你转化成系数表达式的，那么还需还需要进行将点值表达式转化为系数表达式。

怎么操作呢，确实很难，但是总有大佬能想到。是一种构造法。

令 $f(x)= \sum\limits_{i=0}^{n}y_i*x^i$ 其中 $y_i$ 就是上面的点值表达式中的所有纵坐标

先给出结论把

$a_k=\frac{f(\omega_n^{-k}\ )}{n}$

证明一波：

$f(\omega_n^{-k})=\sum\limits_{i=0}^{n}y_i*\omega_n^{-ik}$

由于 $y_i=\sum\limits_{j=0}^{n}\omega_n^{ij}*a_j$

代入

$f(\omega_n^{-k})=\sum\limits_{i=0}^{n}[\sum\limits_{j=0}^{n}\omega_n^{ij}*a_j]*\omega_n^{-ik}$

像莫比乌斯反演一样套路调换顺序

$f(\omega_n^{-k})=\sum\limits_{j=0}^{n}\sum\limits_{i=0}^{n}[\omega_n^{ij}*a_j]*\omega_n^{-ik}$

整理整理

$f(\omega_n^{-k})=\sum\limits_{j=0}^{n}a_j\sum\limits_{i=0}^{n}(\omega_n^{j-k})^i$

这时候思考一下后面的东西

令 $S(x)=1+x+x^2+...+x^{n-1}$
$xS(x)=x+x^2+x^3+...+x^n$
两式相减 $(x-1)S(x)=x^n-1$
$S(x)=\frac{x^n\ -1}{x-1}$

将 $x=\omega_n^u$ 代入

$\\S(\omega_n^u)=\frac{(\omega_n^n)^u-1}{\omega_n^u-1}=\frac{1-1}{\omega_n^u-1}\\$

显然，上式为 $0$ ,但是，除非 $\omega_n^u-1 = 0$ 即分母为 $0$ 时，思考一下什么意思就是 $n=u$ 时呗。根据定义，此时 $S(\omega_n^u)=n$

好的，处理了这些，将 $S(x)$ 代回去

$f(\omega_n^{-k})=\sum\limits_{j=0}^{n}a_j\sum\limits_{i=0}^{n}(\omega_n^{j-k})^i$

发现当且仅当 $j=k$ 时原式才有值。 $f(\omega_n^{-k})=a_k*n$

证到这里，两边再除以 $n$ ，答案就出来了。

当然，为了求值，再回眸看一眼这个 $f(\omega_n^{-k})$ 有无何特殊的地方，这不就是刚才算的正变换把 $\omega_n^{k}$ 换成了 $\omega_n^{-k}$ 么？

所以，正变换和逆变换是两个相反的运算，从某种方面理解，就是多转了半圈。

FFT 递归版 die 码

#include <bits/stdc++.h>
#define debug puts("I ak IOI several times");
#define pb push_back
using namespace std;
template <typename T>inline void read(T& t){
    t=0; register char ch=getchar(); register int fflag=1;
    while(!('0'<=ch&&ch<='9')){if(ch=='-') fflag=-1;ch=getchar();}
    while(('0'<=ch&&ch<='9')){t=((t<<1)+(t<<3))+ch-'0'; ch=getchar();} t*=fflag;
}
template <typename T,typename... Args> inline void read(T& t, Args&... args){
    read(t);read(args...);
}
template <typename T>inline void write(T x){
    if(x<0) putchar('-'),x=~(x-1); int s[40],top=0;
    while(x) s[++top]=x%10,x/=10; if(!top) s[++top]=0;
    while(top) putchar(s[top--]+'0');
}
const double PI=acos(-1.0);
const int MAXN=2500005;
struct Complex{
	double x,y;
	Complex(double _x=0,double _y=0){x=_x,y=_y;}
}a[MAXN],b[MAXN],c[MAXN];
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
void FFT(Complex *y,int len,int type){
	if(len<=1) return;
	int mid=len>>1;
	Complex a1[mid],a2[mid];
	for(int i=0;i<mid;++i) a1[i]=y[i<<1],a2[i]=y[i<<1|1];
	FFT(a1,mid,type);FFT(a2,mid,type);
	Complex Wn=Complex(cos(2*PI/len),type*sin(2*PI/len)),w(1,0);
	for(int i=0;i<mid;++i){
		Complex wt=w*a2[i];
		y[i]=a1[i]+wt;
		y[i+mid]=a1[i]-wt;
		w=w*Wn;
	}
}
int n,m,limit;
int main(){
	read(n,m);
	for(int i=0;i<=n;++i) cin>>a[i].x;
	for(int i=0;i<=m;++i) cin>>b[i].x;
	limit=1;while(limit<=n+m) limit<<=1;
	FFT(a,limit,1); FFT(b,limit,1);
	for(int i=0;i<=limit;++i) c[i]=a[i]*b[i];
	FFT(c,limit,-1);
	for(int i=0;i<=n+m;++i) cout<<(int)(c[i].x/limit+0.5)<<' ';
    return 0;
}
//Welcome back,Chtholly.

递推优化 FFT

蝴蝶操作

观察到在原来的操作中有这样一步。

Complex wt=w*a2[i];
y[i]=a1[i]+wt;
y[i+mid]=a1[i]-wt;
w=w*Wn;

现在称其为蝴蝶操作
类似于下图

发现规律

对于原来的递归版 FFT ，不难发现，很多一大部分是将数组拆分的过程，那这个数组拆分的有没有规律呢？

还真有，这里继续引用一下谷里的题解图片。

这时候把最后一行全部转成二进制：

$000,100,010,110,001,101,011,111$

reverse一下

$000,001,010,011,100,101,110,111$

转成十进制发现就是 $0,1,2,3,4,5,6,7$ 那么这个规律好用了呀。

只需要提前构造一下数组就 ok 了，将数组变成 $0,4,2,6,1,5,3,7$ 的模样似乎就不用拆分了？ 是的。

这里想一下如何 reserve 。

令 $i$ 号反转成 $rev_i$ 。

如果我正着做，显然 $rev_{\left\lfloor\frac{i}{2}\right\rfloor}$ 时已经处理好了的。

想想 $rev_{\left\lfloor\frac{i}{2}\right\rfloor}$ 和 $rev_i$ 有没有什么关系。

$\left\lfloor\frac{i}{2}\right\rfloor$ 不就是 $i$ 右移一位嘛，那么它如果全部反转的话应该在第一位。好了，那么只需要看原来的 $i\&1$ 是不是 $1$ 就行了，是 $1$ 把开头一位换成 $1$ .

for(int i=0;i<limit;++i) rev[i]=(rev[i>>1]>>1)|((i&1)*(1<<L-1));

die 码一行，简洁易懂常数小

迭代 FFT

知道了这些，迭代 FFT 只不过是一步之遥。只需要想想 die 码的细节就好了，这里给出一张宝贵的图。

我感觉这张图已经非常好理解力。

下面上 die 码：

#include <bits/stdc++.h>
#define debug puts("I ak IOI several times");
#define pb push_back
using namespace std;
template <typename T>inline void read(T& t){
    t=0; register char ch=getchar(); register int fflag=1;
    while(!('0'<=ch&&ch<='9')){if(ch=='-') fflag=-1;ch=getchar();}
    while(('0'<=ch&&ch<='9')){t=((t<<1)+(t<<3))+ch-'0'; ch=getchar();} t*=fflag;
}
template <typename T,typename... Args> inline void read(T& t, Args&... args){
    read(t);read(args...);
}
template <typename T>inline void write(T x){
    if(x<0) putchar('-'),x=~(x-1); int s[40],top=0;
    while(x) s[++top]=x%10,x/=10; if(!top) s[++top]=0;
    while(top) putchar(s[top--]+'0');
}
const double PI=acos(-1.0);
const int MAXN=2500005;
struct Complex{
	double x,y;
	Complex(double _x=0,double _y=0){x=_x,y=_y;}
}a[MAXN],b[MAXN],c[MAXN];
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
int rev[MAXN],limit,L,n,m;
void FFT(Complex *y,int type){
	for(int i=0;i<limit;++i) if(rev[i]>i) swap(y[i],y[rev[i]]);
	for(int len=2;len<=limit;len<<=1){
		Complex Wn(cos(2*PI/len),type*sin(2*PI/len));
		int mid=len>>1;
		for(int L=0,R=len-1;R<=limit;L+=len,R+=len){
			Complex w(1,0);
			for(int p=L;p<L+mid;++p){
				Complex X=y[p]+w*y[p+mid],Y=y[p]-w*y[p+mid];
				y[p]=X;y[p+mid]=Y;
				w=w*Wn;
			}
		}
	}
}
int main(){
	read(n,m);
	for(int i=0;i<=n;++i) cin>>a[i].x;
	for(int i=0;i<=m;++i) cin>>b[i].x;
	limit=1;while(limit<=n+m) limit<<=1,++L;
	for(int i=0;i<limit;++i) rev[i]=(rev[i>>1]>>1)|((i&1)*(1<<L-1));
	FFT(a,1);FFT(b,1);
	for(int i=0;i<=limit;++i) c[i]=a[i]*b[i];
	FFT(c,-1);
	for(int i=0;i<=n+m;++i) cout<<(int)(c[i].x/limit+0.5)<<' ';
    return 0;
}
//Welcome back,Chtholly.