[多项式] FFT小计

引入

给出两个多项式 $A,B$ ，计算它们相乘的结果。

我们能轻易写出 code：

for(int i=0;i<=n;i++)
	for(int j=0;j<=n;j++)
		C[i+j]+=A[i]*B[j];

然后超时了。
FFT 是一种将多项式乘法优化成 $O(n\log n)$ 的神仙算法。
接下来会介绍傅里叶的四个神仙 trick。

分析

上面的式子没有任何优化空间。
什么意思呢？就是怎么乱搞都优化不了。于是我们看看天才算法是怎么搞的。

点值表示法

我们知道 $n$ 个点确定一条 $n-1$ 次多项式。
所以我们把多项式转到平面上。
我们随便取 $2n+1$ 个点，先在 $A$ 中计算出这些点，表示成 $(x,y_A)$，然后在 $B$ 中计算出这些点：$(x,y_B)$ ，然后我们就能很容易计算出 $y_C=y_A\times y_B$ 。然后我们用 $n$ 个点确定一条多项式。

天才思路！但是是不是有什么不对的地方？
计算一个点是 $O(n)$ 的，然后有 $n$ 个点要计算...时间复杂度 $O(n^2)$。

如果只是随便取点，就是低估了傅里叶的上限了。

递归计算

我们想想我们有没有知道一个点 $x_1$ 的值，就能马上另一点的值呢？
我们先来看二次函数 $f(x)=x^2$ ：
如果知道 $f(5)=25$，马上就能知道 $f(-5)=25$。
再来看看三次函数 $f(x)=x^3$，和二次函数一样，但是加个负号就行，说人话就是奇偶函数。

我们推广到多项式这一般形式：
$A(x)=16x^5+22x^4+7x^3+5x^2+13x+8$
我们把奇偶次数分开一下：

$$p=16x^5+7x^3+13x$$

$$q=22x^4+5x^2+8$$

假设我们算出 $x$ 点的取值为 $p+q$ ，那么 $-x$ 点的取值不就是 $q-p$ 吗！

也就是说，我们现在只需要取 $\frac{n}{2}$ 这么多点就行了。

然后能不能继续拓展？
这个看起来能递归的样子：
我们定义函数 $D(A)$ 表示给这个函数一个 $n$ 次的多项式，它会返回 $n+1$ 个点和它们的取值。
我们看看这个式子：

$$A(x)=x(16x^4+7x^2+13)+(22x^4+5x^2+8)$$

我们看看括号的式子，我们把 $x^2$ 看成 $x$ ，这不就是子问题吗！
所以递归求解就行了！
举个例子，假设我们计算 $x_1,x_2,x_3,x_4,-x_1,-x_2,-x_3,-x_4$，我们计算子问题 $x_1^2,x_2^2,x_3^2,x_4^2$ 的取值就行了。

等等，是不是有些问题？
新要计算的数全是平方项，都是正数，不存在相反数，所以递归不成立(悲)。所以难道就只能结束了吗？

取点

这里是 FFT 算法最天才的部分。
我们想想刚才的递归。
假设第一轮取的四个数是 $x_1,x_2,-x_1,-x_2$
那么第二轮的两个数是 $x_1^2,-x_1^2$
最后一轮取的数就是 $x^4$，好消息是，最后一轮这个 $x$ 我们可以随便取，所以我们令 $x=1$。
我们往下推，$x^4=x_1^2$，所以我们取 $x_1=1$，所以我们已经确定了两项：$x_1=1,-x_1=-1$。
我们看另一部分， $x_2^2=-x_1^2=-1$，所以 $x^2_2=-1$ 这不显然无解吗。

所以傅里叶最天才的部分就到了，实数部分没有，取虚数不就行了！
所以 $x_2=i$。

我们继续看上面的例子。
我们计算子问题 $x_1^2,x_2^2,x_3^2,x_4^2$ 的取值。
假设我们满足 $x_2^2=-x_1^2,x_4^2=-x_3^2$

这个时候我们计算 $x_1^4，x_3^4$ 就行。
最后，我们发现，其实我们就是解一个方程： $x^8=1$

整个流程是一棵二叉树。
我们想要对于每个节点，它的两个叶子节点成相反数，当前节点的值等于叶子节点值的平方。
看这个图，也可以说，FFT 一次性就是在处理每一层的信息，我们也能感受到 FFT 的时间复杂度：$O(n\log n)$ 。
问题来了：这些数怎么取呢？
这部分不是 FFT 重点，请看 这个。

知道单位根以后，这部分就很简单了。
我们看看这个图：

$A:\omega_2^0(=1),B:\omega_2^1$
$D:\omega_4^0(=1),C:\omega_4^2,F:\omega_4^1,E:\omega_4^3$
所以，对于任意一个点，若是 $\omega_n^k$，那么它的儿子分别是 $\omega_{2n}^k,\omega_{2n}^{k+n}$。
于是，我们在 $O(n\log n)$ 的优秀时间复杂度内完成了取点。

点值转多项式

这个怎么办呢？
拉插的 $O(n^2)$ 明显无法接受。
我们知道最终的点值表示：$G(k)=\sum\limits_{i=0}^{n-1}F(i)(\omega_n^k)^i$
emm这个正推很搞，直接给出结论： $F(k)\times n=\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^iG(i)$

我们反证。
我们把 $G(i)$ 给带入一下。
$\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^iG(i)$
$=\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^i\sum\limits_{j=0}^{n-1}F(j)(\omega_n^i)^j$
$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}\omega_n^{-ik}F(j)\omega_n^{ij}$
$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}F(j)\omega_n^{i(j-k)}$

我们分类讨论一下。

• $1.j=k$
  此时原式子等价于 $nF(k)$。

• $2.j\ne k$
  我们定义 $p=j-k$。
  原式等价于
  $=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}F(p+k)\omega_n^{ip}$
  我们把 $j$ 枚举提前，相当于枚举 $p$ 。
  $=\sum\limits_{j=0}^{n-1}\sum\limits_{i=0}^{n-1}F(p+k)\omega_n^{ip}$

我们看这个部分。
$\sum\limits_{i=0}^{n-1}F(p+k)\omega_n^{ip}$
$=F(p+k)\sum\limits_{i=0}^{n-1}\omega_n^{ip}$

玄学的地方就在这里：

$$\sum\limits_{i=0}^{n-1}\omega_n^{ip}=0(p\ne 0)$$

设 $S=\sum\limits_{i=0}^{n-1}(\omega_n^{p})^i$
则 $\omega_n^pS-S=\omega_n^n-\omega_n^0=0$
由于 $p\ne 0$，所以 $\omega_n^p\ne 1$，所以 $S=0$。所以这一部分没有贡献。

所以我们看到这个公式：
$F(k)\times n=\sum\limits_{i=0}^{n-1}(\omega_n^{-k})^iG(i)$
$\times n$ 那一项我们最后消掉即可。然后我们发现这个式子怎么长得和得到的式子这么像啊，所以我们直接对 $G$ 做转点值，但是系数是 $\omega_n^0,\omega_n^{-1},\omega_n^{-2}...\omega_n^{-(n-1)}$。这与 $\omega_n^0,\omega_n^{n-1},\omega_n^{n-2}...\omega_n^{1}$ 等价。所以可以进行 $FFT$ 加速。

还是详细写一下这里：
举个例子：
$\omega_4^{0},\omega_4^{-1},\omega_4^{-2},\omega_4^{-3}$
那么根据 DFT 能得到：
$\omega_8^{0},\omega_8^{-1},\omega_8^{-2},\omega_8^{-3},\omega_8^{4},\omega_8^{3},\omega_8^{2},\omega_8^1$
改造一下就是：
$\omega_8^{0},\omega_8^{-1},\omega_8^{-2},\omega_8^{-3},\omega_8^{-4},\omega_8^{-5},\omega_8^{-6},\omega_8^{-7}$

所以可以进行迭代。
至此，经过了四个神仙优化的锤炼后，理论部分终于完工。

代码部分

我们知道原理后，代码变得异常好写。

#include<bits/stdc++.h>
#define ld double
#define N 2000005
using namespace std;
const ld Pi=acos(-1);
int n,m;
struct CP{
	ld x,y;
}f[N*2],g[N*2],tmp[N*2];
CP operator +(CP a,CP b){return (CP){a.x+b.x,a.y+b.y};}
CP operator -(CP a,CP b){return (CP){a.x-b.x,a.y-b.y};}
CP operator *(CP a,CP b){return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(CP *f,int len,bool flag)
{
	if(len==1) return;
	CP *fl=f,*fr=f+len/2;
	for(int i=0;i<len;i++) tmp[i]=f[i];
	for(int i=0;i<len/2;i++)
		fl[i]=tmp[(i<<1)],fr[i]=tmp[(i<<1)|1];
	FFT(fl,len/2,flag),FFT(fr,len/2,flag);
	CP cur=(CP){1,0},w=(CP){cos(2*Pi/len),(flag?1:-1)*sin(2*Pi/len)};
	for(int i=0;i<len/2;i++)
	{
		tmp[i]=fl[i]+cur*fr[i];
		tmp[i+len/2]=fl[i]-cur*fr[i];
		cur=cur*w;
	}
	for(int i=0;i<len;i++) f[i]=tmp[i];
}
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	for(m=n+m+1,n=1;n<m;n<<=1);
	FFT(f,n,1);FFT(g,n,1);
	for(int i=0;i<n;i++) f[i]=f[i]*g[i];
	FFT(f,n,0);
	for(int i=0;i<m;i++)
		printf("%d ",(int)(f[i].x/n+0.5));
	return 0;
}

虽然没超时，但是很慢。
其中 $*f$ 是指向 $f$ 数组的指针。$f+x$ 可以取到 $f$ 往后 $x$ 位的指针。

接下来介绍迭代实现。

迭代实现

我们回到这个性质：

对于一个操作 $\omega_n^k$，它的左节点是 $\omega_{2n}^k$，右节点是 $\omega_{2n}^{n+k}$

画的太屎了

我们可以发现，在第 $x$ 层往右走，那么它会加上 $2^x$，往左则不变。
我们再仔细想想这个操作：

• 如果往左走，不变。
• 如果往右走，在最高位补 $1$。

我们分析最后一层第 $6$ 个数(不算 $0$)。
将 $6$ 进行二进制拆分，显然是 $110$。
那么我们能分析它的行走规律：第 $1$ 层右走，第 $2$ 层右走，第 $3$ 层左走。
所以它的值就是 $2^0+2^1=3$。
我们可以发现，每个数本质上就是当前的数翻转。
这个东西叫蝴蝶变换。

我们把最后一层的搞出来以后，往后的就直接搞就行了。
这样就把递归去了。

所以现在问题就是……
怎么快速求出一个数 $x$ 的二进制翻转呢？
$n\log n$ 固然可做，但是太慢了。

假设我们已经处理出 $0\to x-1$，怎么求 $x$ 的翻转？

举个例子：
我们求 $(00110)1$ 的翻转。
我们已经知道了 $0(00110)$ 的翻转。
所以我们先得到 $0(00110)$ 的翻转：$(01100)0$
我们知道，括号部分是我们想要的，我们右移一位就得到了。
最后判断第一位是不是 $1$，是的话在最高位补就 ok 了。

数列求完后，我们想想怎么改变原数组，我们知道开一个新数组是很浪费时间的。我们想要原地搞。

比如说，$x$ 的翻转是 $!x$ ，$!x$ 的翻转就是 $x$。所以我们交换一下这两个位置的数就行了。

我们的迭代 FFT 火热出炉！

#include<bits/stdc++.h>
#define ld double
#define N 2000005
using namespace std;
const ld Pi=acos(-1);
int n,m;
int tr[N*2];
struct CP{
	ld x,y;
}f[N*2],g[N*2];
CP operator +(CP a,CP b){return (CP){a.x+b.x,a.y+b.y};}
CP operator -(CP a,CP b){return (CP){a.x-b.x,a.y-b.y};}
CP operator *(CP a,CP b){return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(CP *f,int n,bool flag)
{
	for(int i=0;i<n;i++)
		if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(int k=2;k<=n;k<<=1)
	{
		int len=(k>>1);
		CP cur=(CP){1,0},w=(CP){cos(2*Pi/k),(flag?1:-1)*sin(2*Pi/k)};
		for(int i=0;i<n;i+=k)
		{
			cur=(CP){1,0};
			for(int j=i;j<i+len;j++)
			{
				CP tmp=cur*f[j+len];
				f[j+len]=f[j]-tmp;
				f[j]=f[j]+tmp;
				cur=cur*w;
			}
		}
	}
}
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
	for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
	for(m=n+m+1,n=1;n<m;n<<=1);
	for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
	FFT(f,n,1),FFT(g,n,1);
	for(int i=0;i<n;i++) f[i]=f[i]*g[i];
	FFT(f,n,0);
	for(int i=0;i<m;i++) printf("%d ",(int)(f[i].x/n+0.5));
	return 0;
}

所以差不多就结束了，这已经是最经典的 FFT 了。

附加：三次变两次 FFT

我们定义一个函数 $T(x)$，$T(x)=F(x)+G(x)i$。
然后给它平方一下。$T(x)=F(x)^2-G(x)^2+2F(x)G(x)i$
所以虚部就是我们想要的，这样两次就能搞了。

所以一般数满足的，多项式一般都有，这个时候多项式就产出来全家桶了。