[笔记]快速傅里叶变换（FFT）

模板题：P3803 【模板】多项式乘法（FFT）

快速傅里叶变换（Fast Fourier Transform,FFT）在算法竞赛中主要用于求卷积，或者说多项式乘法。如果我们枚举两数的各系数相乘，时间复杂度是\(O(n^2)\)，而FFT可以将这一过程优化到\(O(n\log n)\)。

流程

整个FFT算法分\(3\)个过程：

• 将\(2\)个多项式的系数表示法转换为点值表示法（FFT）。
• 将\(2\)个多项式的点值相乘，得到结果的点值表示。
• 将结果的点值表示转换为系数表示法（IFFT）。

其中系数表示法就是形如\(\sum\limits_{i=0}^{n-1}a_i x^i\)的表示。
将其转换为点值表示，就是选了\(n\)个具体的点\(x_0,x_1,\dots,x_{n-1}\)，对应在该点的值\(y_1,y_2,\dots,y_{n-1}\)。
可以证明，任意互不相同的\(n\)个点值可以确定一个唯一的系数表示。

其中，我们发现步骤\(2\)仅仅需要将每个\(x\)对应的的\(2\)个\(y\)相乘，就可以得到结果的点值表示了，这一步骤是\(O(n)\)的。

所以我们重点分析步骤\(1\)和\(3\)。

注：以下步骤假设\(n\)是\(2\)的正整数次幂，具体代码实现中，为了保证正常运行，也会将\(n\)补成\(2\)的正整数次幂。

系数转点值

现在我们要选取\(n\)个点，将\(A,B\)两个多项式转为点值表示，这里拿\(A\)说明。

如果我们直接代数进去，时间复杂度仍然是\(O(n^2)\)。

接下来就要上FFT的核心了。我们可以将\(A\)拆成\(1\)个偶函数和\(1\)个奇函数的和，也就是把\(A\)按奇偶次项划分成\(2\)部分。

\[\begin{aligned} A(x)&=a_0+a_1x^1+a_2x^2+a_3x^3+\dots+a_{n-1}x^{n-1}\\ &=(a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2})+(a_1+a_3x_3+a_5x^5+\dots+a_{n-1}x^{n-1}) \end{aligned}\]

设

\[A_1(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{\frac{n-2}{2}} \]

\[A_2(x)=a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{\frac{n-2}{2}} \]

则

\[A(x)=A_1(x^2)+x\times A_2(x^2) \]

\[A(-x)=A_1(x^2)-x\times A_2(x^2) \]

为什么呢？因为\(A_1\)是偶函数，而\(x\times A_2\)是奇函数。

我们发现如果我们令\(x_0,x_1,x_2,\dots,x_{n-1}\)是正负配对的，其中\([x_0,\dots,x_{\frac{n}{2}-1}]\)和\([x_\frac{n}{2},\dots,x_{n-1}]\)对应相反，而它们的平方也对应相等……

我们就仅需算出\(A_1(x_0^2),A_1(x_1^2),\dots,A_1(x_{\frac{n}{2}-1}^2)\)和\(A_2(x_0^2),A_2(x_1^2),\dots,A_2(x_{\frac{n}{2}-1}^2)\)，然后用上面的式子就可以\(O(n)\)计算出每个\(A(x_0),A(x_1),\dots,A(x_{n-1})\)了。

这相当于我们将该问题分成了\(2\)个规模减半的子问题。对于每个子问题，再重复上面的操作。

如果这个方法可行，我们就能再\(O(n\log n)\)的时间复杂度内求出点值表达。

但我们发现对于初始状态，我们可以令\(x_0,x_1,x_2,\dots,x_{n-1}\)正负配对。但是在实数范围内它们的平方\([x_0^2,x_1^2,\dots,x_{\frac{n}{2}-1}^{2}]\)都是非负数，怎么能保证正负配对呢？

没错，就是复数。

如上图，上一层的平方是下一层，而下一层可用于推出上一层，和我们刚才的过程相同。

我们发现，最终要求的\(x_0,x_1,\dots,x_{n-1}\)，就是\(x^8=k\)在复数域上的解。

这里不妨取\(k=1\)，这样我们解出的\(8\)个根，就是“\(n\)次单位复数根”。

如上图，根据欧拉公式，有\(e^{i\theta}=\cos \theta+i\sin \theta\)。那么我们可以得出对于\(k\in[0,n)\)，\(n\)次单位根\(\omega_n^{k}=\omega^{2\pi i*\frac{k}{n}}\)（因为\(2\pi\)是\(360^{\circ}\)嘛，\(k\)表示的就是\(n\)份中的第几份）

单位根\(\omega\)有如下性质（可以结合上面的图来理解）：
\(\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}}\)
\(\omega_{an}^{ak}=\omega_{n}^{k}\)

发现了么？第\(1\)个性质和\(x\)的对应关系\(x_{i}=x_{i+\frac{n}{2}}\)正好相同！

而第\(2\)个性质说明平方之后，这\(\frac{n}{2}\)个单位根仍然平均分布在单位圆上。

也就是说，我们可以用\(\omega\)代替\(x\)：

至此，我们可以写出代码了。不过写之前我们发现有一个小优化：

如果我们每层额外用一个\(O(n)\)来进行奇偶分类，虽然不影响总时间复杂度\(O(n\log n)\)，但是常数较大。我们试着观察奇偶分类后的形态（图片来源自为风月马前卒的博客）。

奇偶分类后的序列，下标的二进制表示正是原序列下标二进制表示的翻转。所以我们根据这个规则，先把\(A\)的系数数组按这样的规则重排一下，然后按正常步骤进行就可以了。

点值转系数

其实我们的多项式也可以表示为矩阵乘的形式：

\[\def\mat#1{\textbf{#1}} \mat{Y}=\mat{X}\times \mat{A}\]

即

\[\begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}=\begin{bmatrix} 1&x_0&x_0^1&\cdots&x_0^{n-1}\\ 1&x_1&x_1^1&\cdots&x_1^{n-1}\\ 1&x_2&x_2^1&\cdots&x_2^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_{n-1}&x_{n-1}^1&\cdots&x_{n-1}^{n-1} \end{bmatrix}\times \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]

带入\([x_0,x_1,\dots,x_{n-1}]=[1,\omega_n^1,\dots,\omega_n^{n-1}]\)：

\[\begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}=\begin{bmatrix} 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)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)(n-1)} \end{bmatrix}\times \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]

现在我们要倒推\(\textbf{A}\)，就两边同时左乘\(\textbf{X}^{-1}\)：

\[\begin{bmatrix} 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)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\cdots&\omega_n^{(n-1)(n-1)} \end{bmatrix}^{-1}\times \begin{bmatrix} y_0\\y_1\\y_2\\ \vdots\\y_{n-1} \end{bmatrix}= \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots\\a_{n-1} \end{bmatrix}\]

中间的矩阵\(\textbf{X}\)是一个范德蒙德矩阵的转置矩阵，其行列式是\(\prod\limits_{1\le i < j \le n }(x_j-x_i)\)（因为是线性代数的范畴所以就不给证明了，逆阵也是）。

由于\(x_i\)互不相等，所以行列式不为\(0\)，存在逆阵：

\[\textbf{X}^{-1}=\frac{1}{n}\begin{bmatrix} 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)}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\omega_n^{-(n-1)}&\omega_n^{-2(n-1)}&\cdots&\omega_n^{-(n-1)(n-1)} \end{bmatrix}^{-1} \]

可以发现，与原矩阵的区别，就是每个\(\omega\)都取了倒数，最后再乘一个\(\frac{1}{n}\)。

所以代码不用重新写了，只需要额外加一个参数type。FFT提供1，IFFT提供-1。IFFT中，从\(\omega_n^0\)逆时针走\(1\)个单位，要变成顺时针。所以单位的复数部分需要\(\times\)type。具体见代码。

Code

点击查看代码

#include<bits/stdc++.h>
#define N 4000010
#define Pi 3.1415926535897932384626
using namespace std;
struct complx{
	double x,y;
	complx(double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
int l,r[N],limit=1;
inline complx operator+(complx a,complx b){return complx(a.x+b.x,a.y+b.y);}
inline complx operator-(complx a,complx b){return complx(a.x-b.x,a.y-b.y);}
inline complx operator*(complx a,complx b){return complx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void fft(complx *a,int type){//采用非递归实现 
	for(int i=0;i<limit;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	//先按照变换规则提前处理出交换后的数组 
	for(int mid=1;mid<limit;mid<<=1){//mid表示分割点（相对位置） 
		complx Wn(cos(Pi/mid),type*sin(Pi/mid));
		//计算单位，每遍历一个单位根，都需要乘一下这个单位 
		for(int R=mid<<1,j=0;j<limit;j+=R){//j表示每组的开始位置 
			complx w(1,0);//从(1,0)开始计算单位根，每乘一次Wn就到了下一个单位根 
			for(int k=0;k<mid;k++,w=w*Wn){//k枚举的是组内元素 
				complx x=a[j+k],y=w*a[j+mid+k];
				//a[j~j+mid-1]是偶数部分，a[j+mid~j+2*mid-1]是奇数部分 
				a[j+k]=x+y;
				a[j+mid+k]=x-y;//原理见FFT过程的说明 
			}
		}
	}
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	int n,m;
	cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i].x;
	for(int i=0;i<=m;i++) cin>>b[i].x;
	while(limit<=n+m) limit<<=1,l++;//把长度变成2的正整数次幂 
	for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//处理出每个数交换后的位置  
	//原理：模拟二进制加法，注意这里的二进制加法是反转意义下的，可以自己试试来理解 
	fft(a,1);//1.FFT 
	fft(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i];//2.计算结果点值 
	fft(a,-1);//2.IFFT
	for(int i=0;i<=n+m;i++) cout<<(int)(a[i].x/limit+0.5)<<" ";//答案需要除以n，+0.5是四舍五入 
	return 0;
}