【数学】FFT蒟蒻的研究历程

本来很久以前就学了，但是我一直没有写博客，直到最近要用的时候才发现忘了。于是就有了这篇博客.
我们知道，如果直接暴力将两个多项式比如
$A(x)=a_0+a_{1}x+a_2{x}^2+...+a_n{x}^n$
$B(x)=b_0+b_{1}x+b_2{x}^2+...+b_n{x}^n$
相乘将会花费$O(n^2)$的时间，有没有更快的方法呢？
我们注意到，耗时的原因是每一项都要花费$O(n)$的时间来乘，合起来乘的难度较大，说不定换一种表示方法就可以不一样了。
我们要确定一个有$n$次项的方程，不仅可以用上面的方法，我们也可以用$n$个点(向量)来表示，而且这样一来方程就可以唯一确定了。
很明显如果是用$n$个点来表示$A(x)$的话,$A(x),B(x)$可以表示为：
$(x_1,A(x_1)),(x_2,A(x_2)),...,(x_n,A(x_n))$
$(x_1,B(x_1)),(x_2,B(x_2)),...,(x_n,B(x_n))$
由于我们要求的$C(x)=A(x)\ast B(x)$因此$C(x)$可以表示为：
$(x_1,A(x_1)\ast B(x_1)),(x_2,A(x_2)\ast B(x_2)),...,(x_n,A(x_n)\ast B(x_n))$
在这种表示方法下我们只需要$O(n)$时间就可以将两个方程相乘,但是我们已有的是$A(x)$和$B(x)$的系数表示法的式子。
但如果暴力转成点值表示法，需要$O(n^2)$的时间，时间复杂度并没有丝毫的减少，而且常数反而变大了。我们需要找到更快的能在系数表示法和点值表示法之间转换的方法。
很明显，我们在带入n个值求点是在做操作重复的运算，而且由于这n个值的取值我们现在还没有确定，这里应该可以通过巧妙的取值来达到减少运算的目的，要是我们选的数的次方之间可以相互轮回，那么应该可以简化运算。一维的实数相乘的话只能在数轴上面走，二维的虚数就比较方便了，因为$(cos(a)+isin(a){)}^n=(cos(na)+isin(na))$我们如果用这样的数，那么就它的n次方就是相当于在平面上转圈，如下图：

因此如上图$(w_8^1{)}^n=w_8^n$.
我们可以发现它有这些性质：
$w_{2n}^{2m}=w_n^m$
$w_n^{m+\frac{n}{2}}=-w_n^m$
因此各次方的关系就有了，我们可以把最后的式子排在一起看看
比如现在只有4个:
$a_0+a_1\ast w_4^0+a_2\ast (w_4^0{)}^2+a_3\ast (w_4^0{)}^3$
$a_0+a_1\ast w_4^1+a_2\ast (w_4^1{)}^2+a_3\ast (w_4^1{)}^3$
$a_0+a_1\ast w_4^2+a_2\ast (w_4^2{)}^2+a_3\ast (w_4^2{)}^3$
$a_0+a_1\ast w_4^3+a_2\ast (w_4^3{)}^2+a_3\ast (w_4^3{)}^3$
即：
$a_0+a_1\ast w_4^0+a_2\ast w_4^0+a_3\ast w_4^0$
$a_0+a_1\ast w_4^1+a_2\ast w_4^2+a_3\ast w_4^3$
$a_0+a_1\ast w_4^2+a_2\ast w_4^4+a_3\ast w_4^6$
$a_0+a_1\ast w_4^3+a_2\ast w_4^6+a_3\ast w_4^9$
由于$w_n^{m+n}=w_n^m$,我们可以把次数模n加强联系即：
$a_0+a_1\ast w_4^0+a_2\ast w_4^0+a_3\ast w_4^0$
$a_0+a_1\ast w_4^1+a_2\ast w_4^2+a_3\ast w_4^3$
$a_0+a_1\ast w_4^2+a_2\ast w_4^0+a_3\ast w_4^2$
$a_0+a_1\ast w_4^3+a_2\ast w_4^2+a_3\ast w_4^1$
由于$w_4^2=-w_4^0,w_4^3=-w_4^1$我们可以进一步把次数降低加强联系即：
$a_0+a_1\ast w_4^0+a_2\ast w_4^0+a_3\ast w_4^0$
$a_0+a_1\ast w_4^1-a_2\ast w_4^0-a_3\ast w_4^1$
$a_0-a_1\ast w_4^0+a_2\ast w_4^0-a_3\ast w_4^0$
$a_0-a_1\ast w_4^1-a_2\ast w_4^0+a_3\ast w_4^1$
可以发现其中的(0,2)(1,3)式的次数是完全一样的，系数的符号有一半一样，另一半互为相反数，其实这很好证明第$m$式和第$m+\frac{n}{2}$式的关系。
$f_m=a_0+a_1\ast w_n^m+a_2\ast w_n^{2m}+...+a_{n-1}\ast w_n^{(n-1)\ast m}$
$f_{m+\frac{n}{2}}=a_0+a_1\ast w_n^{m+\frac{n}{2}}+a_2\ast w_n^{2m}+...+a_{n-1}\ast w_n^{(n-1)\ast (m+\frac{n}{2})}$
即
$f_{m+\frac{n}{2}}=a_0-a_1\ast w_n^m+a_2\ast w_n^{2m}+...-a_{n-1}\ast w_n^{(n-1)\ast m}$(假设n是偶数)
那么我们可以把这之间相同的部分拿出来
令
$A_0=a_0+a_2\ast w_n^{2m}+a_4\ast w_n^{4m}+...+a_{n-2}\ast w_n^{(n-2)\ast m}$
$A_1=a_1+a_3\ast w_n^{2m}+a_5\ast w_n^{4m}+...+a_{n-1}\ast w_n^{(n-2)\ast m}$
则
$A_0=a_0+a_2\ast w_{\frac{n}{2}}^m+a_4\ast w_{\frac{n}{2}}^{2m}+...+a_{n-2}\ast w_{\frac{n}{2}}^{(n-2)\ast m/2}$
$A_1=a_1+a_3\ast w_{\frac{n}{2}}^m+a_5\ast w_{\frac{n}{2}}^{2m}+...+a_{n-1}\ast w_{\frac{n}{2}}^{(n-2)\ast m/2}$
那么很明显
$f_m=A_0+A_1\ast w_n^m$
$f_{m+\frac{n}{2}}=A_0+A_1\ast w_n^{m+\frac{n}{2}}=A_0-A_1\ast w_n^m$
因此我们只要把$f_m$和$f_{m+\frac{n}{2}}$求出来就可以一下求出来$A_0$和$A_1$ 而且$f$的数据范围只有$A$的一半。也就是说只要$lo{g}_2(n)$层运算就可以求出解了，说以总的时间复杂度只有$nlo{g}_2(n)$.

$a_0+a_1{w}_4^0+a_2{w}_4^0+a_3{w}_4^0$	$a_0+a_1{w}_4^1+a_2{w}_4^2+a_3{w}_4^3$	$a_0+a_1{w}_4^2+a_2{w}_4^4+a_3{w}_4^6$	$a_0+a_1{w}_4^3+a_2{w}_4^6+a_3{w}_4^9$
$a_0+a_2{w}_2^0$	$a_0+a_2{w}_2^1$	$a_1+a_3{w}_2^0$	$a_1+a_3{w}_2^1$
$a_0$	$a_2$	$a_1$	$a_3$

为了方便处理我们将对应位置的f对着对应位置的A，表示成数字就是:
0 1 2 3 4 5 6 7
0 2 4 6 1 3 5 7
0 4 2 6 1 5 3 7
不难发现之后得到的位置的二进制就是按反过来之后排序，因为我们第一次把&1=0的移动到了前面，而吧&1=1的移到了后面,所以之后第二次把&2=0的移动到了前面，而吧&2=1的移到了后面，之后依次这样就可以了。这就是蝴蝶变换。
因此我们在实际操作的时候就反着来，从下至上依次求出每一层，这里不递归的原因是因为递归太慢了。
说以我们现在可以用$O(nlo{g}_2(n))$的时间复杂度正着把系数表示法变成点值表示法，那么怎么倒着装呢？
很明显为了求出上面的值，我们采取的方法是倒着来，既然
$f_m=A_0+A_1\ast w_n^m$
$f_{m+\frac{n}{2}}=A_0-A_1\ast w_n^m$
那么
$2\ast A_0=f_m+f_{m+\frac{n}{2}}$
$2\ast A_1=(f_m-f_{m+\frac{n}{2}})\ast w_n^{-m}$
我们先暴力把最终得到的式子求出来再统一除$n(n=2^{某次方})$,得到:

$4\ast a_0$	$4\ast a_1$	$4\ast a_2$	$4\ast a_3$
$2\ast (a_0+a_2{w}_4^0)$	$2\ast (a_1+a_3{w}_4^0)$	$2\ast (a_0+a_2{w}_4^2)$	$2\ast (a_1+a_3{w}_4^2)$
$a_0+a_1{w}_4^0+a_2{w}_4^0+a_3{w}_4^0$	$a_0+a_1{w}_4^2+a_2{w}_4^4+a_3{w}_4^6$	$a_0+a_1{w}_4^1+a_2{w}_4^2+a_3{w}_4^3$	$a_0+a_1{w}_4^3+a_2{w}_4^6+a_3{w}_4^9$

比如这里n就是4，FFT只能解决n是二的整次幂的情况，如果开始的n不是2的整次幂那么加一堆系数为0的项就可以了，反正加的项又不到一倍，时间也多不了多少。
这样式子就可以反过来了
不过貌似再写逆变换有点麻烦，大佬们搞了个逆矩阵至今没搞懂，留坑了。
这里是代码(洛谷3803)

#include<cstdio>
#include<complex>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int MAXN=4000005;
typedef complex<double> cpx;
const long double PI=acos(0.0)*2.0;
int N,M;
void FFT(cpx a[],int n,double op)
{
	for(int i=1,j=0;i<n-1;i++)
	{
		for(int s=n;j^=(s>>=1),(~j)&s;);
		if(i<j)
			swap(a[i],a[j]);
	}
	for(int len=2,mid=1;len<=n;len<<=1,mid<<=1)
	{
		cpx wm(cos(2.0*PI/len),op*sin(2.0*PI/len));
		for(int now=0;now<n;now+=len)
		{
			cpx w(1,0);
			for(int k=now;k<now+mid;k++,w=w*wm)
			{
				cpx l=a[k],r=w*a[k+mid];
				a[k]=l+r,a[k+mid]=l-r;
			}
		}
	}
	if(op==-1)
	{
		for(int i=0;i<n;i++)
			a[i]/=n;
	}
}
cpx A[MAXN],B[MAXN];
void mul(int a[],int n,int b[],int m,int c[])
{
	//memset(A,0,sizeof(A));
	//memset(B,0,sizeof(B));
	for(int i=0;i<n;i++)
		A[i]=cpx(a[i],0);
	for(int i=0;i<m;i++)
		B[i]=cpx(b[i],0);
	int l=n+m-1,len=1;
	while(len<l)
		len<<=1;
	FFT(A,len,1);
	FFT(B,len,1);
	for(int i=0;i<len;i++)
		A[i]=A[i]*B[i];
	FFT(A,len,-1);
	for(int i=0;i<len;i++)
		c[i]=int(A[i].real()+0.5);
}
int F[MAXN],G[MAXN],an[MAXN];
int main()
{
	scanf("%d %d",&N,&M);
	N++,M++;
	for(int i=0;i<N;i++)
		scanf("%d",&F[i]);
	for(int i=0;i<M;i++)
		scanf("%d",&G[i]);
	mul(F,N,G,M,an);
	for(int i=0;i<N+M-1;i++)
		printf("%d%c",an[i],i==N+M-2?'\n':' ');
}