FFT学习笔记（快速傅里叶变换）

用途

快速傅里叶变换（Fast Fourier Transformation，简称FFT） 一般用来加速多项式乘法。求两个$n$次多项式相乘，朴素算法需要$O(n^2)$，但FFT只需要$O(n\log n)$就能解决。

多项式

系数表示法

一个$n$次$n+1$项多项式可表示为
$$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$$

这就是多项式的系数表示法。

点值表示法

将$n+1$个不同的$x$值代入$f(x)$，可以得到$n+1$个对应的$f(x)$值。这$n$对$(x_k,f(x_k))$可以唯一确定一个$n$次多项式。

为什么呢？可以将$n+1$个式子联立起来得到一个$n+1$元一次方程组，这个方程组有唯一解。

用$n+1$个点表示一个$n$次函数，称为点值表示法。

系数转为点值

对于一个$n$次多项式，在已知系数的情况下，将其转为点值表示法，称为DFT（离散傅里叶变换）。

使用秦九韶算法，即可$O(n)$转换一个$x$值，那么求$n$个点的时间复杂度为$O(n^2)$。

点值转为系数

对于一个由点值表示法表示的$n$次多项式，将其转化为系数表示法，称为IDFT（离散傅里叶逆变换）。

用高斯消元需要$O(n^3)$，用拉格朗日插值法可以达到$O(n^2)$。

多项式乘法

对于用系数表示的$n$次多项式$f(x)$和$g(x)$，求两式相乘后的多项式，时间复杂度为$O(n^2)$。

用点值表示的呢？

对于$f(x)$和$g(x)$，我们各取$2n+1$个点。

$$\begin{gathered} (x_0,f(x_0)),(x_1,f(x_1)),\dots ,(x_{2n},f(x_{2n})) \\ (x_0,g(x_0)),(x_1,g(x_1)),\dots ,(x_{2n},g(x_{2n})) \end{gathered}$$

上下相乘得到

$$(x_0,f(x_0)\times g(x_0)),(x_1,f(x_1)\times g(x_1)),\dots ,(x_{2n},f(x_{2n})\times g(x_{2n}))$$

显然这$2n+1$个点唯一地确定了一个$2n$次函数，这$2n+1$个点就是$f(x)$和$g(x)$的乘积的点值表示法的形式，所以用点值表示法只需要$O(n)$。

那如果将系数表示法转换为点值表示法相乘，再转回系数表示法，会不会更快呢？

用朴素算法是不会的，因为DFT 和IDFT 的时间复杂度都为$O(n^2)$，所以最终的时间复杂度为$O(n^2)$，和直接相乘一样。

有没有其他办法呢？那就要用到FFT 了。

快速傅里叶变换

复数

复数可以写成$a+bi$的形式。在复平面上为点$(a,b)$。

$x^n=1$有$n$个根，令${\omega }_n$表示$x^n=1$的单位复数根，则${\omega }_n=e^{2\pi i/n}$，则方程其他的根就是${\omega }_n$的若干次幂。

由欧拉公式可得
$$e^{i\theta }=\cos \theta +i\sin \theta$$

所以我们可以通过三角函数来求$e^{2\pi i/n}$的实部和虚部。也由此可得，$x^n=1$的根都均匀分布在复平面上以原点为圆心的单位圆的圆周上。

单位复数根有以下性质：

1. 消去引理

对于任意整数$n\ge 0,k\ge 0,d>0$，都有${\omega }_{dn}^{dk}={\omega }_n^k$

证明： ${\omega }_{dn}^{dk}=(e^{2\pi i/dn}{)}^{dk}=(e^{2\pi i/n}{)}^k={\omega }_n^k$

2. 折半引理

对于偶数$n\ge 0$和整数$k\ge 0$，$({\omega }_n^{k+n/2}{)}^2={\omega }_{n/2}^k$

证明： $({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}={\omega }_{n/2}^k$

3. 求和引理

对于任意整数$n\ge 1$和不能被$n$整除的整数$k$，有$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=0$

证明： $\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{1-({\omega }_n^k{)}^n}{1-{\omega }_n^k}=\frac{0}{1-{\omega }_n^k}=0$

FFT

我们考虑将${\omega }_n^0,{\omega }_n^1,{\omega }_n^2,\dots ,{\omega }_n^{n-1}$代入$n-1$次多项式求值。
（下文将$n$视作2的幂，次数大于等于$n$的项可看作系数为0）

有多项式
$$A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}$$

将$A(x)$根据下标奇偶性分开
$$\begin{gathered} A(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1}) \\ =(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+\cdots +a_{n-1}{x}^{n-2}) \end{gathered}$$

令

$$\begin{gathered} A_1(x)=a_0+a_{2}x+\cdots +a_{n-2}{x}^{n/2-1} \\ {A}_2(x)=a_1+a_{3}x+\cdots +a_{n-1}{x}^{n/2-1} \end{gathered}$$

可得
$$A(x)=A_1(x^2)+x{A}_2(x^2)$$

设$0\le k<n/2$，将$x={\omega }_n^k$代入$A(x)$，得
$$\begin{gathered} A({\omega }_n^k)=A_1(({\omega }_n^k{)}^2)+{\omega }_n^k{A}_2(({\omega }_n^k{)}^2) \\ {A}_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})=A_1({\omega }_{n/2}^k)+{\omega }_n^k{A}_2({\omega }_{n/2}^k) \end{gathered}$$

将$x={\omega }_n^{k+n/2}$代入$A(x)$，得
$$\begin{gathered} A({\omega }_n^{k+n/2})=A_1(({\omega }_n^{k+n/2}{)}^2)+{\omega }_n^{k+n/2}{A}_2(({\omega }_n^{k+n/2}{)}^2) \\ {A}_1({\omega }_n^{2k})+{\omega }_n^k\cdot {\omega }_n^{n/2}{A}_2({\omega }_n^{2k})=A_1({\omega }_{n/2}^k)-{\omega }_n^k{A}_2({\omega }_{n/2}^k) \end{gathered}$$

我们可以发现，$A({\omega }_n^k)$和$A({\omega }_n^{k+n/2})$用$A_1(x)$和$A_2(x)$表示之后只有符号不同。
$$\begin{gathered} A({\omega }_n^k)=A_1({\omega }_{n/2}^k)+{\omega }_n^k{A}_2({\omega }_{n/2}^k) \\ A({\omega }_n^{k+n/2})=A_1({\omega }_{n/2}^k)-{\omega }_n^k{A}_2({\omega }_{n/2}^k) \end{gathered}$$

于是，当我们求出$A_1({\omega }_{n/2}^k)$和$A_2({\omega }_{n/2}^k)$之后，就能求出$A({\omega }_n^k)$和$A({\omega }_n^{k+n/2})$的值。

那么问题就转化为求$A_1(x)$和$A_2(x)$在${\omega }_{n/2}^0,{\omega }_{n/2}^1,\dots ,{\omega }_{n/2}^{n/2-1}$处的值，可以用分治来解决问题。

以下是递归写法，用于理解，后面会介绍更优的迭代写法。

code

const double pi=acos(-1.0);
int a1[N],a2[N],y[N];
int fft(int *a,int l){
	if(l==1) return a;
	cp wn=(cp){cos(2*pi/l),sin(2*pi/l)};
	cp w=(cp){1,0};
	for(int i=0;i<l/2;i++){
		a1[i]=a[i*2];
		a2[i]=a[i*2+1];
	}
	int *y0=fft(a1,l/2);
	int *y1=fft(a2,l/2);
	for(int i=0;i<l/2;i++){
		y[i]=y0[i]+y1[i]*w;
		y[i+l/2]=y0[i]-y1[i]*w;
		w=w*wn;
	}
	return y;
}

其结构与$cdq$分治相似，相当于有$\log n$层，每层都是$n$个位置，所以时间复杂度为$O(n\log n)$。

IFFT

我们已经可以在$O(n\log n)$的时间复杂度下将系数表示法转化为点值表示法。接下来，我们要考虑如何将两个多项式的乘积的点值表示法再变回系数表示法。

对于前面的FFT，其过程可以用以下矩阵表示：
$$\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right]$$

左边的矩阵是范德蒙德矩阵，$V_{i,j}={\omega }_n^{ij}$，可以构造出它的逆矩阵$V^{-1}$，其中$(V^{-1}{)}_{i,j}=\frac{{\omega }_n^{-ij}}{n}$

$$\left[\begin{array}{ccccc}1/n& 1/n& 1/n& \cdots & 1/n\\ 1/n& {\omega }_n^{-1}/n& {\omega }_n^{-2}/n& \cdots & {\omega }_n^{-(n-1)}/n\\ 1/n& {\omega }_n^{-2}/n& {\omega }_n^{-4}/n& \cdots & {\omega }_n^{-2(n-1)}/n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1/n& {\omega }_n^{-(n-1)}/n& {\omega }_n^{-2(n-1)}/n& \cdots & {\omega }_n^{-(n-1)(n-1)}/n\end{array}\right]\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right]=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]$$

相当于在做FFT时在原来的基础上将${\omega }_n$变成${\omega }_n^{-1}$，然后对结果的每一项除以$n$即可，时间复杂度也为$O(n\log n)$。

蝴蝶变换

由上文可得知，可以将$A(x)$拆成$A_0(x)$和$A_1(x)$来求值。

但是用上述的递归代码跑的速度较慢，所以我们考虑用更高效的算法。

观察分治的过程，看看有什么规律。

我们发现，最底层的系数排列顺序，正好是其二进制的镜面翻转。

如第$1$个，$1$的二进制是$001$，翻转后是$100$，即十进制的$4$，所以从第$0$个位置开始的第一个位置是$a_4$。

为什么呢？因为在每次分治的时候，我们都将低位为$0$的放在左，为$1$的放在右。但如果将每个数的下标镜面翻转，我们就发现这是在将$a$序列按二进制位从高到低来分，最后自然就是从小到大。再镜面翻转回来，即可得到原来的系数。这个变换称为蝴蝶变换。

代码中的$j$是$i$翻转后的值，$j$在这个循环中其实就是在模拟$i$进位的过程。

code

void ch(cp *a,int l){
	for(int i=1,j=l/2,k;i<l-1;i++){
		if(i<j) swap(a[i],a[j]);
		k=l/2;
		while(j>=k){
			j-=k;k>>=1;
		}
		j+=k;
	}
}

以下是FFT的迭代写法。

void fft(cp *a,int l){
	for(int i=2;i<=l;i<<=1){
		wn=(cp){cos(2*pi/i),sin(2*pi/i)};
		for(int j=0;j<l;j+=i){
			w=(cp){1,0};
			for(int k=j;k<j+i/2;k++,w=w*wn){
				cp t=a[k],u=w*a[k+i/2];
				a[k]=t+u;
				a[k+i/2]=t-u;
			}
		}
	}
}

总结

解决多项式乘法的过程：

• 用FFT将系数表示法转化为点值表示法
• 用点值表示法来求多项式乘法
• 用IFFT将点值表示法再变回系数表示法

这样即可用$O(n\log n)$的时间复杂度求出多项式的乘积。

例题

多项式乘法（FFT）

模板题

令$l$表示最后多项式的次数，则$n$需满足$n\ge l$且$n$为$2$的若干次幂。

code

#include<bits/stdc++.h>
using namespace std;
const int N=5000005;
const double pi=acos(-1.0);
char s1[N],s2[N];
long long ans[N];
struct cp{
	double a,b;
	cp operator +(const cp ax)const{
		return (cp){a+ax.a,b+ax.b};
	}
	cp operator -(const cp ax)const{
		return (cp){a-ax.a,b-ax.b};
	}
	cp operator *(const cp ax)const{
		return (cp){a*ax.a-b*ax.b,b*ax.a+a*ax.b};
	}
}w,wn,a1[N],a2[N];
void ch(cp *a,int l){
	for(int i=1,j=l/2,k;i<l-1;i++){
		if(i<j) swap(a[i],a[j]);
		k=l/2;
		while(j>=k){
			j-=k;k>>=1;
		}
		j+=k;
	}
}
void fft(cp *a,int l,int fl){
	for(int i=2;i<=l;i<<=1){
		wn=(cp){cos(fl*2*pi/i),sin(fl*2*pi/i)};
		for(int j=0;j<l;j+=i){
			w=(cp){1,0};
			for(int k=j;k<j+i/2;k++,w=w*wn){
				cp t=a[k],u=w*a[k+i/2];
				a[k]=t+u;
				a[k+i/2]=t-u;
			}
		}
	}
}
int main()
{
	int n,l,l1,l2,x;
	scanf("%d%d",&l1,&l2);++l1;++l2;
	l=l1+l2;n=1;
	while(n<l) n<<=1;
	for(int i=0;i<l1;i++) scanf("%d",&x),a1[l1-i-1]=(cp){(double)(x),0};
	for(int i=l1;i<n;i++) a1[i]=(cp){0,0};
	for(int i=0;i<l2;i++) scanf("%d",&x),a2[l2-i-1]=(cp){(double)(x),0};
	for(int i=l2;i<n;i++) a2[i]=(cp){0,0};
	ch(a1,n);ch(a2,n);
	fft(a1,n,1);
	fft(a2,n,1);
	for(int i=0;i<n;i++) a2[i]=a1[i]*a2[i];
	ch(a2,n);
	fft(a2,n,-1);
	for(int i=0;i<l1+l2-1;i++){
		ans[i]=(long long)(a2[i].a/n+0.5);
	}
	while(ans[l]==0&&l>0) --l;
	for(int i=l;i>=0;i--) printf("%lld ",ans[i]);
	return 0;
}

A*B Problem 升级版（FFT 快速傅里叶变换）

将$a$和$b$看作两个多项式$f(x)$和$g(x)$，求多项式的乘积，再代入$x=10$即可。

code

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
char s1[5000005],s2[5000005];
long long ans[5000005];
struct cp{
	double a,b;
	cp operator +(const cp ax)const{
		return (cp){a+ax.a,b+ax.b};
	}
	cp operator -(const cp ax)const{
		return (cp){a-ax.a,b-ax.b};
	}
	cp operator *(const cp ax)const{
		return (cp){a*ax.a-b*ax.b,b*ax.a+a*ax.b};
	}
}w,wn,a1[5000005],a2[5000005];
void ch(cp *a,int l){
	for(int i=1,j=l/2,k;i<l-1;i++){
		if(i<j) swap(a[i],a[j]);
		k=l/2;
		while(j>=k){
			j-=k;k>>=1;
		}
		j+=k;
	}
}
void fft(cp *a,int l,int fl){
	for(int i=2;i<=l;i<<=1){
		wn=(cp){cos(fl*2*pi/i),sin(fl*2*pi/i)};
		for(int j=0;j<l;j+=i){
			w=(cp){1,0};
			for(int k=j;k<j+i/2;k++,w=w*wn){
				cp t=a[k],u=w*a[k+i/2];
				a[k]=t+u;
				a[k+i/2]=t-u;
			}
		}
	}
}
int main()
{
	int n,l,l1,l2;
	scanf("%s%s",s1,s2);
	l1=strlen(s1);
	l2=strlen(s2);
	l=l1+l2;n=1;
	while(n<l) n<<=1;
	for(int i=0;i<l1;i++) a1[l1-i-1]=(cp){(double)(s1[i]-'0'),0};
	for(int i=l1;i<n;i++) a1[i]=(cp){0,0};
	for(int i=0;i<l2;i++) a2[l2-i-1]=(cp){(double)(s2[i]-'0'),0};
	for(int i=l2;i<n;i++) a2[i]=(cp){0,0};
	ch(a1,n);ch(a2,n);
	fft(a1,n,1);
	fft(a2,n,1);
	for(int i=0;i<n;i++) a2[i]=a1[i]*a2[i];
	ch(a2,n);
	fft(a2,n,-1);
	for(int i=0;i<l1+l2-1;i++){
		ans[i]+=(long long)(a2[i].a/n+0.5);
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
	while(ans[l]==0&&l>0) --l;
	for(int i=l;i>=0;i--) printf("%lld",ans[i]);
	return 0;
}