多项式乘法

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

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

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

证明：设这个多项式$f=a_n{x}^n+a_{n-1}{x}^{n-1}+...+a_{1}x+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_{1}x+a_0$

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

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

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

$w_n^k$具有如下的性质：

1.$\forall i\ne j,w_n^i\ne w_n^j$

2.$w_n^k=\cos (\frac{2k\pi }{n})+i\sin (\frac{2k\pi }{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_2{x}^2+a_0,h=a_{n-1}{x}^{n-1}+a_{n-3}{x}^{n-3}+...+a_{1}x$

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

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

也就是说我们要求出将$w_n^0,w_n^1,...,w_n^{n-1}$代入$f$的值，只需要递归求解${\phi }_1,{\phi }_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}$

证明：
$$\begin{gathered} \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 \\ =n{a}_k \end{gathered}$$
，最后一步成立是因为
$$\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=\{\begin{array}{l}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\ne k\end{array}$$

于是我们设$\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$，则每一层的划分如下

$$\begin{gathered} a_0{a}_1{a}_2{a}_3{a}_4{a}_5{a}_6{a}_7 \\ {a}_0{a}_2{a}_4{a}_6|a_1{a}_3{a}_5{a}_7 \\ {a}_0{a}_4|a_2{a}_6|a_1{a}_5|a_3{a}_7 \\ {a}_0|a_4|a_2|a_6|a_1|a_5|a_3|a_7 \end{gathered}$$

我们现在需要快速获得最后一行的顺序。观察，第一行的二进制位与最后一行对应的二进制位刚好为翻转关系。比如第一行的$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;
}