数学黑科技1——FFT

用途

• 让我们看一下下面这一个问题：
• 对于$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$,$B(x)=b_0+b_{1}x+b_2{x}^2+...+b_{m-1}{x}^{m-1}$
• 求$C(x)=A(x)\ast B(x)$

暴力

• 我们显然可以枚举A的每一次项，B的每一次项。
• 然后一一乘起来，就是答案了！
• $$c_i=\sum _{j=0}^{j=i}{a}_j\ast b_{i-j}$$
• 但是这样子做的时间复杂度$O(n\ast m)$的，显然很慢！
• 所以，我们要想办法优化

点值与插值

• 在介绍更优的算法之前，就让我们先了解另一个多项式的表示方式
• 在通常情况下，我们是用“系数表示法”，也就是刚才的形如：
• $A(x)=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$的形式
• 但其实，还有另一种表示法，叫做“点值表示法”。
• 对于一个最高次数为$n-1$的多项式，我们可以选$n$个不同的数$(x_0,x_1,...,x_{n-1})$代入次多项式，得出$n$个结果$(y_0,y_1,...,y_{n-1})$
• 可以证明，只要给出这$2\ast n$个数，$x$与$y$，就可以得出有且仅有一个多项式
• 证明，考虑变成矩阵乘法：
• 只要证明“x矩阵”是可逆的，就可以用$a=y\ast x^{-1}$来得出解。
• 我们把x矩阵变成范德蒙德行列式，可以通过证明得出右图的等式（数学归纳法）
• 我们知道这是不为0的，因为这些数两两不相同。
• 所以x矩阵是可逆的，a也就是唯一的
• 而我们把点值变回原多项式的操作叫做插值

瓶颈

• 既然我们只要知道了点值与插值，我们就可以把$A(x)$与$B(x)$用相同的$x$带入，然后得出两组不同的$y$，再把这两组$y$一一乘起来，最后用插值求出$C(x)$就行了。
• 但很可惜的是，朴素的点值和插值不管怎么样，似乎都是$O(n^2)$的，这样子还不如直接暴力呢！！！（说那么多有什么用？）
• 但不要就此放弃，就在这时，傅里叶站了出来，他找到了突破口。

突破口

• 他提出了一个关键的思想：如果把点值时取的$x$特殊化，会不会有什么新发现呢？
• 顺着这个思想，他找到了n次单位复根。

n次单位复根

• 这个名词听起来很高深，实际也挺高深的。
• 说白了，就是如果$x^n=1$，它就是n次单位复根。
• 在实数数域内，可能符合条件的只有1和-1，而且-1得在偶数次方的情况下才成立。
• 但如果推广到复数，满足条件的就有n个数了。
• 首先来简单地讲一下复数是什么。
• 其实，实数已经包含了所有存在的数，但就在这时，人们发现又有一类数，它们虽然不存在，但在解决一些数学问题中有重要的作用，这类数的核心就是——$\sqrt{-1}$，也就是$i$。
• 在很多情况下，我们把复数数域看成是一个平面直角坐标系，x轴表示实数部分，y轴表示虚数部分。例如一个数$3+5\sqrt{-1}=3+5i$，它在这个坐标系上对应的点就是$（3,5）$。
• 复数作为一种数，自然也有乘法
• $(a+bi)\ast (c+di)=(ac-bd)+(ad+bc)i$，显然
• 但如果把它放在坐标系上，乘积也是有规律的。就是模长相乘，幅角相加。
• 模长指的是这个复数所在的点与原点的距离，而幅角就是这一条连线与x轴正半轴的夹角。
• 具体证明这里就不赘述了。
• 那么，对于一个复数，它的n次方的模长是原数的n次方，而幅角就是原数的n倍。
• 好的，下面是关键部分：
• 我们知道，如何一个非零数的0次方都为1，所以$(\sqrt[n]{1}{)}^0=1$（这里的次方根是在复数数域内的）
• 而$(\sqrt[n]{1}{)}^n=1$也是成立的。
• 也就是说在乘了$n$次$\sqrt[n]{1}$之后，模长没有变（1），幅角则是转了360度之后又回到了原位。
• 所以，$w_n^1$（我们一般把$(\sqrt[x]{1}{)}^y$表示为$w_x^y$）在平面坐标系上模长为1，幅角为$\frac{2\pi }{n}$。
• 那么，剩下的就是一些简单的三角函数了，这里就不赘述了，最终$w_n^1$的坐标就是$(cos(\frac{2\pi }{n}),sin(\frac{2\pi }{n}))$。

有关的定理

群的性质

• 因为$w_n^0=w_n^n=1$
• 所以$w_n^i\ast w_n^j=w_n^{i+j}=w_n^{(i+j)\%n}$
• 于是有$w_n^{-1}=w_n^{n-1}$，$w_n^k=w_n^{k+n}$，$w_n^k=w_n^{2k}$

消去引理

• $w_{dn}^{di}=w_n^i$

折半引理

• 对于任意正整数n，2*n次单位根的平方的集合与n次单位根集合相同。
• 也就是 { ( w 2 n ) 2 } = { w n } \{(w_{2n})^2\}=\{w_n\} {(w2n​)2}={wn​}

求和引理

• $$\sum _{j=0}^{n-1}(w_n^k{)}^j=0$$
• 证明，通过等比数列求和可得到原式等于$\frac{(w_n^k{)}^n-1}{w_n^k-1}$
• 其分子=$w_n^{nk}-1$=$1-1$=$0$
• 得证。

FFT

• 好了，铺垫了那么久，终于来到了最终的正题。
• 我们把这特殊的n次单位根带入，就会有特殊的结果。
• 这特殊的结果就是多项式A的离散傅里叶变换（DFT）
• 紧接着，我们考虑分治。
• 首先，我们把n补成形如2^x（x为正整数）的数，原本没有的项就把系数补为0。
• 然后，我们考虑这样的一个多项式$A(x)$：
• $$a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1}$$
• 我们把奇数次项和偶数次项分开：
• $$(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+...+a_{n-1}{x}^{n-1})$$
• 于是，我们把这个多项式分成了两个多项式：
• $$A_0(x^2)+x\ast A_1(x^2)$$
• 在这里，项数也就变为了$n/2$，一直到变为1为止，这时候，只剩下1个常数项，就返回1个代入$w_1^0$（因为有多少项就要有多少个点值），也就是当前的这个系数。
• 假设我们已经处理好了$A_0$与$A_1$两个$w_{n/2}^0$~$w_{n/2}^{n/2-1}$的点值，考虑如何合并。
• 对于代入$w_n^k$，有以下两种情况：
• 1. $k<\frac{n}{2}$，那么$A(w_n^k)=A_0(w_{n/2}^k)+w_n^k\ast A_1(w_{n/2}^k)$，直接相加即可（非常简单）
• 2.否则，$A(w_n^{k+n/2})=A_0(w_{n/2}^{k+n/2})+w_n^{k+n/2}\ast A_1(w_{n/2}^{k+n/2})=A_0(w_{n/2}^k)-w_n^k\ast A_1(w_{n/2}^k)$
• 然后，就直接合并就好了（是不是比想象中的要简单呢？）
• 然后就这样合并${\log }_{2}n$次，就可以得到点值了！！！！！
• 这样点积的时间复杂度就是$O(nlo{g}_{2}n)$

插值

• 点积已经被我们KO了，那么插积怎么搞呢？
• 实际上，DFT有一个非常巧妙的结论。
• 对于一个多项式A，我们求出它的DFT，然后把得出来的$(y_0,y_1,...,y_{n-1})$当做是一个新的多项式B的系数，然后用n次单位复根的的倒数作为代入值再对B求一遍DFT，得出一个$z$
• 可证$a_i=\frac{z_i}{n}$
• 证明，我们先把$z_k$用式子求出来
• $$z_k=\sum _{i=0}^{n-1}{y}_i\ast (w_n^{-k}{)}^i=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j\ast (w_n^i{)}^j)\ast (w_n^{-k}{)}^i$$
• $$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j\ast (w_n^{-k}{)}^i\ast (w_n^i{)}^j=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j\ast w_n^{(j-k)\ast i}$$
• $$=\sum _{j=0}^{n-1}{a}_j\ast \sum _{i=0}^{n-1}{w}_n^{(j-k)\ast i}$$
• 若$j=k$，${\sum }_{i=0}^{n-1}{w}_n^{(j-k)\ast i}$显然为$n$，不然，就是${\sum }_{i=0}^{n-1}(w_n^i{)}^{j-k}$，通过求和引理我们可以知道这个式子为0.
• 所以$a_k=\frac{z_k}{n}$，得证
• 这样子的话，就相当于把插值变成了点值，时间复杂度为$O(nlo{g}_{2}n)$

一些关于实现的东西

• 直接把所有分治时点值存下来是十分不现实的。
• 所以我们最好开滚动数组，这样子合并起来就很方便，但注意复数的运算要打好点，不要让常数太大（毕竟是复数嘛）。

代码（多项式乘法）（洛谷3803 模板题）

#include<cstdio>
#include<cstring>
#include<cmath>
#define pi M_PI
using namespace std;
struct xushu{
	double x,y;
	xushu operator +(const xushu &p){
	  return {x+p.x,y+p.y};
	}
	xushu operator -(const xushu &p){
	  return {x-p.x,y-p.y};
	}
	xushu operator *(const xushu &p){
	  return {x*p.x-y*p.y,x*p.y+y*p.x};
	}
	xushu operator /(const int &p){
	  return {x/p,y/p};
	}
};
xushu a[2100000][2],b[2100000];
int read(){
	int x=0;char ch=getchar();
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
	return x;
}
int mymax(int x,int y) {return x>y?x:y;}
int abs(int x) {return x>0?x:-x;}
inline xushu mi(int a,int b){
	if(a<0) a+=b;
	xushu t={cos(2*pi/b*a),sin(2*pi/b*a)};
	return t;
}
int g[2100000],t;
void find(int len,int bk){
	for(int i=0;i<len;i++){
		if(g[i]>i)
			{xushu t=a[i][0];a[i][0]=a[g[i]][0];a[g[i]][0]=t;}
	}
	int d=1;t=0;
	while(1){
		len/=2;d*=2;
		if(!len) break;t=1-t;
		for(int i=1;i<=len;i++){
			xushu one=mi(bk,d),temp={1,0};
			for(int j=0;j<d/2;j++){
				int t2=(i-1)*d+j,t3=t2+d/2;
				a[t2][t]=a[t2][1-t]+temp*a[t3][1-t];
				a[t3][t]=a[t2][1-t]-temp*a[t3][1-t];
				temp=temp*one;
			}
		}
	}
}
xushu c[2100000],d[2100000];
xushu ans[2100000];
int main()
{
	int n,m,xx;scanf("%d%d\n",&n,&m);n++;m++;
	for(int i=0;i<n;i++) a[i][0].x=read();
	for(int i=0;i<m;i++) b[i].x=read();
	int len=n+m-1,d1=1,d2=0;
	while(d1<len) d1*=2,d2++;len=d1;
	for(int i=0;i<len;i++){
		int e=i,h=0;
		for(int j=1;j<=d2;j++) h=h*2+e%2,e/=2;
		g[i]=h;
	}
	find(len,1);
	for(int i=0;i<len;i++) c[i]=a[i][t],a[i][0]=b[i];
	find(len,1);
	for(int i=0;i<len;i++) d[i]=a[i][t];
	for(int i=0;i<len;i++) c[i]=c[i]*d[i];
	for(int i=0;i<len;i++) a[i][0]=c[i];
	find(len,-1);
	for(int i=0;i<len;i++) ans[i]=a[i][t]/len;
	for(int i=0;i<n+m-1;i++) printf("%.0lf ",ans[i].x+0.001);
	printf("\n");
	return 0;
}