多项式乘法

FFT主要用于快速求多项式的乘积。多项式的乘积就叫做卷积

对\(F\)和\(G\)来说，显然暴力算法的复杂度是\(O(nm)\)，而FFT的时间复杂度为\(O(nlogn)\)

多项式的性质：用任意\(n+1\)个横坐标不同的点，可以唯一确定一个\(n\)次多项式。这个性质叫做多项式的点表示法

证明：设这个多项式\(f=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0\)，那么我们代入\(n+1\)个点可以获得\(n+1\)个方程，这是\(n+1\)元一次方程，写出系数行列式就会发现是范德蒙德行列式，由于横坐标不同，所以行列式的值不为\(0\)，于是方程有唯一解

设\(H=FG\)，于是\(H\)为\(n+m\)次多项式，我们选取\(n+m+1\)个点代入\(F,G\)之后即可确定\(H\)。这就是点表示法的好处，时间复杂度从\(O(nm)\)降低到了\(O(n+m)\)

于是我们现在要解决的问题是如何快速将一个多项式的系数表示法（也就是我们通常写出来的多项式）与点表示法互相转换

下文的\(f(x)=a_{n-1}x^{n-1}+...+a_1x+a_0\)

系数表示法转点表示法：我们要在取\(n+1\)个特殊点上下功夫，我们取复数域上的单位根

假设\(n=2^c\)（如果不是的话我们补齐），画出复数域上的单位圆：

设\(w_n^k\)表示将单位圆等分成\(n\)份，将\(x\)轴逆时针旋转\(\frac{2πk}{n}\)所得到的复数，其中\(k\)的取值为\([0,n-1]\)。比如上图，\(w_n^0=1\)，\(w_n^3=-i\)

\(w_n^k\)具有如下的性质：

1.\(\forall i\neq j,w_n^i\neq w_n^j\)

2.\(w_n^k=\cos(\frac{2kπ}{n})+i\sin(\frac{2kπ}{n})\)

3.\(w_{2n}^{2k}=w_n^k\)

4.\(w_n^{k+\frac{n}{2}}=-w_n^k\)

现在，我们将\(f\)奇数项和偶数项分开，有\(f=g+h\)，其中\(g(x)=a_{n-2}x^{n-2}+a_{n-4}x^{n-4}+...+a_2x^2+a_0,h=a _ {n-1} x ^ {n-1} + a _ {n-3}x^{n-3}+...+a_1x\)

将\(g\)中每个\(x\)用\(x^2\)代替可设\(φ_1(x)=a_{n-2}x^{\frac{n}{2}-1}+a_{n-4}x^{\frac{n}{2}-2}+...+a_2x+a_0\)，将\(h\)提取一个\(x\)并将剩下的式子的\(x\)用\(x^2\)代替可得\(φ_2(x)=a_{n-1}x^{\frac{n}{2}-1}+a_{n-3}x^{\frac{n}{2}-2}+...+a_1\)，于是\(f(x)=φ_1(x^2)+xφ_2(x^2)\)

现在，设\(k∈[0,\frac{n}{2}-1]\)，则\(f(w_n^k)=φ_1(w_n^{2k})+w_n^kφ_2(w_n^{2k})=φ_1(w_{\frac{n}{2}}^{k})+w_n^kφ_2(w_{\frac{n}{2}}^{k}),f(w_n^{k+\frac{n}{2}})=φ_1(w_n^{2k})-w_n^kφ_2(w_n^{2k})=φ_1(w_{\frac{n}{2}}^{k})-w_n^kφ_2(w_{\frac{n}{2}}^{k})\)

也就是说我们要求出将\(w_n^0,w_n^1,...,w_n^{n-1}\)代入\(f\)的值，只需要递归求解\(φ_1,φ_2\)即可；一共会划分\(O(\log n)\)层，每层的总复杂度为\(O(n)\)，所以时间复杂度为\(O(n\log n)\)

点表示法转换为系数表示法：设我们现在已经知道了\((w_n^0,f(w_n^0)),...,(w_n^{n-1},f(w_n^{n-1}))\)，我们要求\(f(x)\)，则有\(a_k=\frac{\sum_{i=0}^{n-1}f(w_n^{i})(w_n^{-k})^i}{n}\)

证明：

\[\sum_{i=0}^{n-1}f(w_n^{i})(w_n^{-k})^i\\=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^{i})^j)(w_n^{-k})^i\\=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(w_n^{j})^i)(w_n^{-k})^i\\=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(w_n^{j-k})^i\\=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(w_n^{j-k})^i\\=na_k \]

，最后一步成立是因为

\[\sum_{i=0}^{n-1}(w_n^{j-k})^i=\begin{cases} n,j=k \\ \frac{1-(w_n^{j-k})^n}{1-w_n^{j-k}}=\frac{1-(w_n^n)^{j-k}}{1-w_n^{j-k}}=\frac{1-1}{1-w_n^{j-k}}=0,j\neq k \end{cases} \]

于是我们设\(\rho(x)=\sum_{i=0}^{n-1}f(w_n^{i})x^{i}\)，则我们求出\(\rho(w_n^{0}),\rho(w_n^{-1}),...,\rho(w_n^{-(n-1)})\)即可，这就是上文讲的傅里叶正变换（系数表示法转化为点表示法）

最后讲一下实现。实践证明，如果用递归来实现的话，常数是非常大的，所以我们一般利用迭代来实现

假设现在\(n=8\)，则每一层的划分如下

\[a_0\space a_1\space a_2\space a_3\space a_4\space a_5\space a_6\space a_7\\a_0\space a_2\space a_4\space a_6\space |\space a_1\space a_3\space a_5\space a_7\\a_0\space a_4|a_2\space a_6 |a_1\space a_5|a_3\space a_7\\a_0|a_4|a_2|a_6|a_1|a_5|a_3|a_7 \]

我们现在需要快速获得最后一行的顺序。观察，第一行的二进制位与最后一行对应的二进制位刚好为翻转关系。比如第一行的\(a_4\)对应最后一行的\(a_1\)，两者的二进制位分别是\(100\)和\(001\)

证明：考虑一个数如何走到最后一层。从第一层到第二层，如果其末尾是\(1\)，那么就会被放到右边，否则就会被放到左边；从第二层到第三层，如果其右移一位后末尾是\(1\)那么就会被放到右边，否则就会被放到左边；依次类推，每次走完一层之后都会右移一位。于是一个数的每个二进制位决定了其每一步是向左走还是向右走。而每往右走一步，当前下标就会加上当前组数的个数的一半，向左走则不变。比如\(a_5\)，先向右走，而当前组（指第一层）有\(8\)个数，于是下标加上\(4\)，再向左走，下标不变，再向右走，而当前组（第三层第三组的{\(a_1,a_5\)}）有两个数，于是下标加上\(1\)，最终的下标就是\(0+4+1=5\)。由上面的过程可知，每次下标加的数就是将当前二进制位放到翻转的位置。于是结论成立

代码见下

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e5+10;
const double PI=acos(-1);
int n,m,len;
complex<double> y[3][N<<2];
int rev[N<<2];
void change(int op)//交换 
{
	for(int i=0;i<len;i++)
	if(i<rev[i]) swap(y[op][i],y[op][rev[i]]);//防止交换两次 
}
void fft(int op,int on)//像BFS一样从下往上一层一层地合并 
{
	change(op);
	for(int h=2;h<=len;h<<=1)//h为当前这一层的区间长度 
	{
		complex<double> wn(cos(2*PI/h),sin(on*2*PI/h));//将单位圆分成h份 
		for(int j=0;j<len;j+=h)//合并起始点 
		{
			complex<double> w(1,0);
			for(int k=j;k<j+h/2;k++)
			{
				complex<double> u=y[op][k],t=w*y[op][k+h/2];
				y[op][k]=u+t,y[op][k+h/2]=u-t;
				//蝴蝶变换，见OI-wiki 
				w*=wn;
			}
		}
	}
}
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n+1;i++)//千万注意这里要加一 
	{
		int a;
		scanf("%d",&a);
		complex<double> u(a,0);//complex是C++自带的一个类 
		y[0][i-1]=u;
	}
	for(int i=1;i<=m+1;i++)
	{
		int a;
		scanf("%d",&a);
		complex<double> u(a,0);
		y[1][i-1]=u;
	}
	len=n+m+1;
	for(int i=0;;i++)//补齐 
	if(len==(1<<i)) break;
	else if(len<(1<<i))
	{
		len=1<<i;
		break;
	}
	for(int i=0;i<len;i++)//这一步操作见OI-wiki的位逆序置换O(n)实现 
	{
		rev[i]=rev[i>>1]>>1;
		if(i&1) rev[i]|=len>>1;
	}
	fft(0,1),fft(1,1);
	for(int i=0;i<len;i++)
	y[2][i]=y[0][i]*y[1][i];
	fft(2,-1);
	for(int i=0;i<=n+m;i++)
	printf("%d ",(int)(real(y[2][i])/len+0.5));//这里都要这么写，不然有精度误差 
    return 0;
}