快速傅里叶变换（FFT）略解

前言

如果我们能用一种时间上比 $O(n^2)$ 更优秀的方法来计算大整数（函数）的乘法，那就好了。快速傅里叶变换（FFT） 可以帮我们在 $O(n\log n)$ 的时间内解决问题。

函数乘积

计算两个大整数之积时，我们发现
$$\begin{gathered} (2x+3)(4x+5)=8x^2+22x+15...(\ast ) \\ 23\times 45=1035 \end{gathered}$$
而如果我们把 $(\ast )$ 式右边的每一位的系数看做一个数每位上的数码，正好得到了 $1035$。事实上，对于所有的多项式乘法，以上规律同样成立。
证明： （提示）考虑竖式乘法的过程，和多项式乘法的过程，它们的本质都是一样的。

这样，我们就把问题转换为：计算两个已知函数之积的函数的解析式。

复平面、单位圆

考虑 $\sqrt{-9}$ 的值。
$$\begin{array}{rl}\sqrt{-9}& =\sqrt{-1}\times \sqrt{9}=3\sqrt{-1}.\end{array}$$
类似地，$\forall N\in {\mathbb{Z}}_{-}$ 我们都可以用类似的方法得到 $$\sqrt{N}=\sqrt{-N}\times \sqrt{-1}$$
引入虚数单位 $\text{i}$，使 ${\text{i}}^2=-1.$ 这样我们就重新认识了数的范围，从实数扩充到复数。

一复数 $a+b\text{i}$ 中的 $a,b\in \mathbb{R}$，$a$ 是它的实数部分，$b\text{i}$ 是虚数部分。若 $b=0$，则它是实数。复数服从实数的大部分运算法则。
若两个复数，它们的实数部分相等，虚数部分之和为 $0$，我们称它们互为 共轭复数。

我们知道，数轴上的每个点与每个实数一一对应。类似地，我们可以使用 复平面 上的点表示复数。复平面与平面直角坐标系类似，它的 $x$ 轴单位长度为 $1$，$y$ 轴单位长度为 $\text{i}$。复平面上的点之横纵坐标表示的数之和即为该点表示的数。比如，$(1,2)$ 表示 $1+2\text{i}$。

以圆点为圆心，$1$ 为半径，在复平面上作圆，如图所示。这个圆叫 单位圆。

在圆上任取一点 $A$，过此点作 $x$ 轴的垂线段，垂足为 $B$。设 $\angle AOB=\theta$。易知$$\begin{gathered} OB=OA\times \cos \theta =\cos \theta \\ AB=OB\times \sin \theta =\sin \theta \end{gathered}$$$\therefore A(\cos \theta ,\sin \theta ).$ 这个 $\theta$ 叫做 $A$ 的 辐角。

单位根及其性质

1. $\forall k,n\in {\mathbb{N}}^{\ast }$ 有 ${\omega }_n^k\times {\omega }_n^1={\omega }_n^{k+1}$；
   证明：
   
   如图所示，过点 $A(1,0)$ 作圆内接 $n$ 边形。设 $\angle AO{A}^{\prime }=\theta .$ 则
   $$\begin{gathered} A^{\prime }(\cos \theta ,\sin \theta ) \\ {A}^{\prime \prime }(\cos 2\theta ,\sin 2\theta ) \\ \begin{array}{rl}& (\cos \theta +\sin \theta \times \text{i})\times (\cos 2\theta +\sin 2\theta \times \text{i})\\ & =\cos \theta \cos 2\theta +(\cos \theta \sin 2\theta +\sin \theta \cos 2\theta )\times \text{i}-\sin \theta \sin 2\theta \\ & =(\sin \theta \sin 2\theta +\cos \theta \cos 2\theta )+(\sin \theta \cos 2\theta +\cos \theta \cos 2\theta )\times \text{i}\end{array} \end{gathered}$$
   由三角函数恒等变换式
   $$\begin{gathered} \sin (\alpha +\beta )=\sin \alpha \sin \beta -\cos \alpha \cos \beta , \\ \cos (\alpha +\beta )=\sin \alpha \cos \beta +\sin \beta \cos \alpha \end{gathered}$$
   知
   $$\begin{array}{rl}原式& =\sin (\theta +2\theta )+\cos (\theta +2\theta )\times \text{i}\\ & =\sin 3\theta +\cos 3\theta \times \text{i}\end{array}$$
   换句话说，$A^{\prime }$ 和 $A^{\prime \prime }$ 表示的数相乘后得到了对应的 $A^{\prime \prime \prime }$。
   设这个多边形的顶点为 { ω i 0 } ,   i ∈ [ 0 , n − 1 ] \{\omega_i^0\},\ i\in[0,n-1] {ωi0​}, i∈[0,n−1]。那么我们一定有$${\omega }_n^k\times {\omega }_n^1={\omega }_n^{k+1},k\in [0,n-1]$$
   特别地，我们有 ${\omega }_n^n={\omega }_n^0$。它们叫做 $n$ 次 单位根。

2. ${\omega }_n^0={\omega }_n^n=1$；

3. ${\omega }_n^0,{\omega }_n^1,...,{\omega }_n^{n-1}$ 互不相同；

4. $\forall k,n\in {\mathbb{N}}^{\ast }$ 有 ${\omega }_n^k={\omega }_{2n}^{2k}$；

证明一： 感性理解。一个 $n$ 边形把单位圆分成了 $n$ 个部分，取这 $n$ 个圆弧的中点，顺次连接这 $2n$ 个点，得到一个新多边形。则 ${\omega }_n^k$ 即表示 $n$ 变形的第 $k$ 个单位根，也表示 $2n$ 变形的第 $2k$ 个点。

证明二：
$\begin{array}{rl}{\omega }_{2n}^{2k}& =\cos 2k\times \frac{2\pi }{2n}+\text{i}\times \sin 2k\times \frac{2\pi }{2n}\\ & =\cos k\times \frac{2\pi }{n}+\text{i}\times \sin k\times \frac{2\pi }{n}\\ & ={\omega }_n^k.\end{array}$

5. $\forall k,n\in {\mathbb{N}}^{\ast }$ 有 ${\omega }_n^{k+\frac{2}{n}}=-{\omega }_n^k$；

证明一： 感性理解。乘以 ${\omega }_n^2$ 的意思其实就是在单位圆逆时针转半圈。单位圆上的一个点，逆时针转了半圈后到达的点，与原来的点关于原点对称。

证明二： $$\begin{array}{rl}{\omega }_n^{k+\frac{2}{n}}& ={\omega }_n^k\times {\omega }_n^{\frac{2}{n}}\\ & ={\omega }_n^k\times (\cos \pi +\text{i}\times \sin \pi )\\ & ={\omega }_n^k\times (-1+0)\\ & =-{\omega }_n^k.\end{array}$$

6. $\forall k,n\in {\mathbb{N}}^{\ast }$ 有 $$\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=\{\begin{array}{l}0,k̸=0\\ n,k=0\end{array}$$
   证明： $(\text{I})$ 若 $k̸=0$，设 $$S=\sum _{i=0}^{n-1}({\omega }_n^k{)}^i...(1)$$ 则 $${\omega }_n^k\times S=\sum _{i=1}^n({\omega }_n^k{)}^i...(2)$$
   $(2)-(1)$ 得$$\begin{array}{rl}({\omega }_n^k-1)\times S& =({\omega }_n^k{)}^n-1\\ S& =\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}\\ & =\frac{1-1}{{\omega }_n^k-1}\\ & =0.\end{array}$$

$(\text{II})$ 当 $k=0$ 时，
$$\begin{array}{rl}原式& =\sum _{i=0}^{n-1}1\\ & =n.\end{array}$$

快速傅里叶变换

显然地，一个 $n$ 次多项式可以被 $n+1$ 个点唯一确定。

那么，我们可以在单位圆上取 $n+1$ 个单位根，代入已知的两个函数，得到 $n+1$ 对点，再把每对点相乘，得到结果函数上的 $n+1$ 个点，再求出结果函数。

设已知函数（合并后）为 $$\begin{array}{rl}f(x)& =\sum _{i=0}^{n-1}{a}_i{x}^i\\ & =(a_0+a_{2}x+...+a_{n-2}{x}^{n-2})\\ & +(a_{1}x+a_3{x}^3+...+a_{n-1}{x}^{n-1})\end{array}$$
设 $$\begin{array}{rl}g(x)& =a_0+a_{2}x+...+a_{n-2}{x}^{\frac{n}{2}-1}\\ h(x)& =a_1+a_{3}x+...+a_{n-1}{x}^{\frac{n}{2}-1}\end{array}$$
则 $$f(x^2)=g(x)+x\times h(x^2)$$
令 $x={\omega }_n^k$，由 单位根的性质1 得 $$\begin{array}{rl}f(x^2)& =g({\omega }_n^{2k})+{\omega }_n^k\times h({\omega }_n^{2k})\\ & =g({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\times h({\omega }_{\frac{n}{2}}^k)...(1)\end{array}$$
令 $x=-{\omega }_n^k={\omega }_n^{k+\frac{n}{2}}$ 得 $$\begin{array}{rl}f(x^2)& =g({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}\times h({\omega }_n^{2k+n})\\ & =g({\omega }_n^{2k})+{\omega }_n^{k+\frac{n}{2}}\times h({\omega }_n^{2k})\\ & =g({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\times h({\omega }_{\frac{n}{2}}^k)...(2)\end{array}$$
不难发现，$(1)$ 式与 $(2)$ 式 互轭！换句话说，当你求出 $(1)$ 式的值时，只需将虚数部分取相反数（时间复杂度为 $O(1)$）即得到了 $(2)$ 式。

这样，我们就将问题的规模减为一半。同理，剩下的一半也可以使用类似方法，分为更小的两半……没错，这就是 分治，这样就将时间复杂度降为了 $O(n\log n)$。

快速傅里叶逆变换

设 { y i } \{y_i\} {yi​} 是 { a i } \{a_i\} {ai​} 的傅里叶变换，即在 { ω n i } \{\omega_n^i\} {ωni​} 处的值；
设 { c i } \{c_i\} {ci​} 是 { y i } \{y_i\} {yi​} 在 { ω n − i } \{\omega_n^{-i}\} {ωn−i​} 的值。
则有 $$\begin{array}{rl}{c}_k& =\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^k{)}^j)\cdot ({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}{a}_j\times \sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i\end{array}$$
根据 单位根的性质6 ，有且仅有一个 $j\in [0,n-1]$ 使得 $j=k$。
$\therefore \forall j$ 有且仅有一个$$({\omega }_n^{j-k}{)}^i=n$$有且仅有 $n-1$ 个 $$({\omega }_n^{j-k}{)}^i=0$$
$$\begin{array}{rl}\therefore \sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i& =n,\\ c_k& =\sum _{i=0}^{n-1}{a}_j\times \sum _{j=0}^{n-1}({\omega }_n^{j-k}{)}^i\\ & =a_k\times n,\\ a_k& =\frac{c_k}{n}.\end{array}$$
换句话说，经过快速傅里叶逆变换后的数组 { c i } \{c_i\} {ci​}，除以 $n$ 后就得到了结果 { a i } \{a_i\} {ai​}。（这里的 $n$ 是指 $c$ 数组的长度）

小结

使用快速傅里叶变换（FFT）计算两个函数之积的步骤如下：

1. 分别对两个函数进行快速傅里叶变换；
2. 将两组结果相乘；
3. 对结果进行快速傅里叶逆变换，并将结果除以 $n$。
   
   

迭代快速傅里叶变换

经过上文的学习，容易写出以下代码（感谢 linjiayang2016 大佬的代码）

void fast_fast_tle(int limit,complex *a,int type){
    if(limit==1)
        return ;
    complex a1[limit>>1],a2[limit>>1];	//*
    for(int i=0;i<=limit;i+=2)
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    fast_fast_tle(limit>>1,a1,type);
    fast_fast_tle(limit>>1,a2,type);
    complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
    for(int i=0;i<(limit>>1);i++,w=w*Wn)
        a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}

发现 $\ast$ 行定义了两个不小的数组，而 FFT 函数是被递归调用的，所以会造成爆栈。

原序列	1	2	3	4	5	6	7	8
快速傅里叶变换后序列	1	5	3	7	2	6	4	8
原下标的二进制	000	100	010	110	001	101	011	111
新下标的二进制	000	001	010	011	100	101	110	111

不难发现，快速傅里叶变换后的序列中，每个数的新下标的二进制，在数值上等于原下标二进制的反转。

设原下标为 $x$ 的数的新下标为 $r[x]$。

以十进制数 $184$ 的二进制表示为例。设这个二进制数长度为 $l$。
$$\begin{array}{ccccccccccccccc}1& & 0& & 1& & 1& & 1& & 0& & 0& |& 0\end{array}$$
我们把这个数分为两部分，第一部分是前 $l-1$ 位，第二部分是最后一位。不难得到，它的翻转就是 在第二部分后面接上第一部分的翻转形成的数。
即

	r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));

至此，快速傅里叶变换的相关内容已全部结束。

$\text{luogu P3803}$

给你两个正整数 $n,m\le 10^6$，和 $n+1$ 次多项式 $f(x)$、$m+1$ 次多项式 $g(x)$。求 $f(x)$ 和 $g(x)$ 的卷积。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>

using namespace std;
#define reg register
const int MAXN=400010;
const double pi=acos(-1.0);

struct compex{
	double x,y;
	friend compex operator +(const compex a,const compex b){
		compex c;
		c.x=a.x+b.x;c.y=a.y+b.y;
		return c;
	}
	friend compex operator -(const compex a,const compex b){
		compex c;
		c.x=a.x-b.x;c.y=a.y-b.y;
		return c;
	}
	friend compex operator *(const compex a,const compex b){
		compex c;
		c.x=a.x*b.x-a.y*b.y;
		c.y=a.x*b.y+a.y*b.x;
		return c;
	}
}a[MAXN],b[MAXN];
int r[MAXN];
int n,m,l=0;
int limit=1;

void FFT(compex*t,int type){
	for(reg int i=0;i<limit;++i)
		if(i<r[i]) swap(t[i],t[r[i]]);
	for(reg int mid=1;mid<limit;mid<<=1){
		compex wn=(compex){cos(pi/mid),type*sin(pi/mid)};
		for(reg int j=0,R=(mid<<1);j<limit;j+=R){
			compex w={1.0,0};
			for(reg int k=0;k<mid;++k,w=w*wn){
				compex x=t[j+k],y=w*t[j+k+mid];
				t[j+k]=x+y;
				t[j+k+mid]=x-y;
			}
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);
	for(reg int i=0;i<=n;++i)
		scanf("%lf",&a[i].x);
	for(reg int i=0;i<=m;++i)
		scanf("%lf",&b[i].x);
	while(limit<=n+m) limit<<=1,++l;
	for(reg int i=1;i<limit;++i)
		r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
	FFT(a,1);FFT(b,1);
	for(reg int i=0;i<limit;++i)
		a[i]=a[i]*b[i];
	FFT(a,-1);
	for(reg int i=0;i<=n+m;++i)
		printf("%d ",(int)(a[i].x/limit+0.5));
}