【luogu P3803】【模板】多项式乘法（FFT）

【模板】多项式乘法（FFT）

题目链接：luogu P3803

题目大意

给你两个多项式，要你求这两个多项式乘起来得到的多项式。（卷积）

思路

系数表示法

就是我们一般来表示一个多项式的方法：
$A(x)=a_1x^k+a_2x^{k-1}+...+a_k$
或者可以这样表示：
$A(x)=\sum\limits_{i=1}^{k}a_i\times x_i$

那你很容易看到，用来做这道题用系数表示法来做是 $O(n^2)$ 的。

点值表示法

假设我们已经知道了这个多项式，那我们把 $n$ 个不同的 $x$ 带入，会得到 $n$ 个取值 $y$。
由于 $n$ 个点可以确定一个 $n$ 次多项式，那我们就可以根据这 $n$ 个点确定这个多项式。

那你如果选了 $n$ 个点 $(x_1,y_1),...,(x_n,y_n)$，那就有：
$y_i=\sum\limits_{j=0}^{n-1}a_j\times x_i^j$

但用它计算还是 $O(n^2)$，你选点 $O(n)$，对于每个点计算也是 $O(n)$。

考虑优化

至于系数表示法，它的系数固定，你一改其它都要改，似乎很难弄。

但是第二种点值法似乎没有特别大的限制让它难以优化。
但是怎么优化呢？这就要看点前置知识才可以继续了。

前置知识

向量

普通的量只有大小，但是向量就是有方向又有大小的量。
在几何里面它其实就是一个箭头。

圆的弧度制

等于半径长的圆弧所对的圆心角叫做1弧度的角，叫做弧度，用 rad 表示。
用它做单位来量角的制度就是弧度制。

$1^{\circ}=\frac{\pi}{180}rad$
$180^{\circ}=\pi rad$

复数

设 $a,b$ 为实数，$i^2=-1$，那形如 $a+bi$ 的数就叫做复数。
$i$ 就是其中的复数单位。

复数其实可以在一个二维平面表示出来，$x$ 轴表示实数 $a$，$y$ 轴（当然在原点的时候不是虚数）表示虚数 $bi$。
然后从 $(0,0)$ 到 $(a,b)$ 的向量就表示了复数 $a+bi$。

模长：就是那个点到原点的距离，勾股就可以得到。
幅角：从 x 正半轴以逆时针为正方向到一直向量的转角。

有关复数的运算法则

那由于复数在平面上是向量，那我们其实可以用向量的加减法则来弄。（就是用那个平行四边形定则）

然后乘法就稍微复杂一点。
几何的意义就是，复数相乘，模长也相乘，然后幅角就会相加。
其实我们这里要的是代数意义，我们化一下式子：
$(a+bi)*(c+di)$
$=ac+a*di+b*ci+bi*di$
$=ac+i(ad+bc)-bd$（$i^2=-1$）
$=(ac-bd)+(ad+bc)i$

单位根

在复平面上，以原点为圆心，$1$ 为半径作圆，所得的圆叫单位圆。

以圆点为起点，圆的 $n$ 等分点为终点，做 $n$ 个向量，设幅角为正且最小的向量对应的复数为 $\omega_n$，称为 $n$次单位根。

那剩下的那 $n-1$ 个就是 $\omega_n^2,\omega_n^3,...,\omega_n^n$。

然后要注意的就是 $\omega_n^0$ 和 $\omega_n^n$ 都是 $1$，对于的是 $x$ 正半轴的向量。

那怎么计算其它的呢？
我们可以用欧拉公式来求：$\omega_n^k=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac{2\pi}{n}$
（大概就是把模长那条线往 $x$ 轴做垂线得到直角三角形，然后用三角函数得到两条直角边的长度，从而得到坐标，也就是复数）

当然显而易见，单位根幅角就是 $\frac{1}{n}$

还有就是如果 $z^n=1$，那 $z$ 就是 $n$ 次单位根。

单位根的一些性质

$\omega_{2n}^{2k}=\omega_n^k$
证明：
$\omega_{2n}^{2k}=\cos 2k\dfrac{2\pi}{2n}+i\sin 2k\dfrac{2\pi}{2n}$
$=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac{2\pi}{n}=\omega_n^k$

$\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$
证明：
$\omega_{n}^{\frac{n}{2}}=\cos \dfrac{n}{2}*\dfrac{2\pi}{n}+i\sin \dfrac{n}{2}*\dfrac{2\pi}{n}$
$=\cos \pi+i\sin\pi=-1$
（$\cos \pi=\cos 180^\circ=-1,\sin\pi=\sin 180^\circ=0$）
那就有 $\omega_{n}^{k+\frac{n}{2}}=\omega_{n}^{k}*\omega_{n}^{\frac{n}{2}}=\omega_{n}^{k}*(-1)=-\omega_{n}^{k}$

OK，现在前置知识结束，开始正文。

快速傅里叶变换（FFT）

我们设 $A(x)$ 的系数为 $(a_0,a_1,...,a_{n-1})$（$n$ 次多项式）
那就有 $A(x)=a_0+a_1x+a_2{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}$

我们可以把它奇偶分开，$A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1+a_3x^3+...+a_{n-1}x^{n-1})$
那我们可以设 $A_1(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1},A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}$

那我们容易看出 $A(x)=A_1(x^2)+xA_2(x^2)$。

把 $\omega_n^k$ 带入，$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$

把 $\omega_n^{k+\frac{n}{2}}$ 带入，$A(\omega_n^k)=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})$，然后把这个式子化一下：
$=A_1(\omega_n^{2k}*\omega_n^n)-\omega_n^{k}A_2(\omega_n^{2k}*\omega_n^n)$
$=A_1(\omega_n^{2k})-\omega_n^{k}A_2(\omega_n^{2k})$

你会发现这两个东西不同的只有一个常数项。
那我们可以枚举得到第一个答案，然后通过第一个答案得到第二个。
那你会发现你只用枚举前面的一半，那问题规模就减半。

那分出的两个问题可以继续分半解决，那我们可以用递归来搞。
那这个其实就是类似分治的感觉，是 $O(nlogn)$ 的。

快速傅里叶逆变换（IFFT）

前面我们在弄的时候都是用点值表示法。

但一般我们（题目）给的多半都是系数表示法。
从系数表示法得到点值表示法我们已经学会了，但是怎么转回来呢？

这个时候就要用到 IFFT 了。
$(y_0,y_1,...,y_{n-1})$ 为 $(a_0,a_1,...,a_{n-1})$ 的傅里叶变换（就是点值表示）
那另外有一个 $(c_0,c_1,...,c_{n-1})$ 满足：
$c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i$
（其实就是把 $(y_0,y_1,...,y_{n-1})$ 当做多项式，求它在 $\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}$ 的点值表示）

然后来推推公式：
$c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i$
$=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i$
$=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i$
$=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)$
$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(\omega_n^{j-k})^i$
$=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i$

我们设 $S(x)=\sum\limits_{i=0}^{n-1}x^i$。
带入 $\omega_n^k$：$S(\omega_n^k)=\sum\limits_{i=0}^{n-1}(\omega_n^k)^i$

然后分类讨论：

当 $k=0$ 时，$\omega_n^k=\omega_n^0=1$。
那式子就变成了 $S(\omega_n^k)=\sum\limits_{i=0}^{n-1}1^i=n$

当 $k\neq0$ 时，我们考虑怎么求。
观察到没项都乘了一个值，我们考虑用这么一种方法：
乘 $\omega_n^k$：$\omega_n^kS(\omega_n^k)=\sum\limits_{i=1}^{n}(\omega_n^k)^i$
然后两个相减：$\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-(w_n^k)^0$
$(\omega_n^k-1)S(\omega_n^k)=(\omega_n^k)^n-1$
$(\omega_n^k-1)S(\omega_n^k)=(\omega_n^n)^k-1$
$(\omega_n^k-1)S(\omega_n^k)=1^k-1$
$S(\omega_n^k)=\dfrac{1-1}{\omega_n^k-1}=0$

OK，分析完了 $S(x)$，我们回到之前卡克的地方：
$c_k=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i$
$=\sum\limits_{j=0}^{n-1}a_jS(\omega_n^{j-k})$
那只会有一个 $j=k$，使得 $j-k=0$，贡献是 $n$，其它的时候，贡献都是 $0$。
那就是 $c_k=a_kn$
那你就会发现，$a_k=\dfrac{c_k}{n}$

那我们就可以通过求出 $(c_0,c_1,...,c_{n-1})$，然后再根据这个一弄就可以得到 $a_k$ 了。
（而且你会发现求 $(c_0,c_1,...,c_{n-1})$ 和求 FFT 的很像，只是由 $\omega_n^i$ 变成了 $\omega_n^{-i}$，那我们只要再搞一个 $op$，记录它单位的改变即可以了）

算法总结

它的主要思想就是把系数表示法转成点值表示法，然后用点值表示法的分支方法来快速得到答案，然后再将答案转回系数表示法。

实现

妹想到把，还有。

别的地方都没有问题，我们来讲讲递归的做法。

递归

（我一定不会告诉你这个是过不了的）

递归其实就是根据我们前面的奇偶分开做法，把一个序列弄成两个，然后递归处理之后合并。

递归到只剩一个常数项的时候，就可以直接返回。

然后我们先来一个小小的优化。
它有一个很牛逼的名字——蝴蝶效应。
（其实就是记忆化）
因为向量的乘法耗比较多时间，我们可以先把要乘的乘出来放一个变量里，然后要用的时候直接就是这个变量。
（这样如果这个乘法值要用多次的话就可以节省时间，原本要乘很多次，现在只用乘一次）

然而就算你加了这个优化，也是过不了的。
因为你递归，加上常数比较大，在 $10^6$ 就爆了。

那我们就要找一个更优的方法来弄，那就是迭代实现。

迭代实现

我们考虑观察一下原序列和翻转之后的序列。

我们观察一下它是怎么变的：
$0,1,2,3,4,5,6,7$
$0,2,4,6,1,3,5,7$
$0,4,2,6,1,5,3,7$

似乎没有什么想法，转成二进制：
$000,100,010,110,001,101,011,111$
再把一开始的也转成二进制：
$000,001,010,011,100,101,110,111$
你会发现它就是把它们的二进制翻转了。

那你可以用一个 $O(n)$ 的方法得到对应关系——DP。
$f_i=(f_{i/2}/2)|((i\&1)^{n-1})$
大概就是把之前翻转好的往左边移一个，流出位置给原序列中右边新出现的一位放。

接着就是怎么通过这个翻转对应得到我们迭代的序列。
现有一个显然的东西，就是 $f_{f_i}=i$

那也就是说，它们之间是两两对应的。
那就需要把每对都只翻转一次，就可以了。
那简单，我们只要找一个在一对里面只会发生一次的事，就比如 $f_i>i$。（当然你小于也可以，只要让一对只有一个会发生就可以了）

然后我们来讲讲迭代要怎么搞。
其实就想到与从下段直接网上合并。（因为你已经确定了最下的序列）
就枚举合并的大小，然后枚举区间，然后就合并。

小声哔哔

终于写完了，好累啊。
awa

代码

#include<cstdio>
#include<cmath>
#include<iostream>

using namespace std;

struct complex {
	double x, y;
	complex (double xx = 0, double yy = 0) {
		x = xx;
		y = yy;
	}
}a[5000005], b[5000005];
int n, m, limit, l_size;
int an[5000005];
double Pi = acos(-1.0);

complex operator +(complex x, complex y) {//定义向量的加减乘
	return complex(x.x + y.x, x.y + y.y);
}

complex operator -(complex x, complex y) {
	return complex(x.x - y.x, x.y - y.y);
}

complex operator *(complex x, complex y) {
	return complex(x.x * y.x - x.y * y.y, x.x * y.y + x.y * y.x);
}

//op=1则系数变点值，为-1则点值边系数
//至于为啥看看求的两个公式就知道了
void FFT(complex *now, int op) {
	for (int i = 0; i < limit; i++)
		if (i < an[i]) swap(now[i], now[an[i]]);//求出要迭代的序列
	
	for (int mid = 1; mid < limit; mid <<= 1) {//枚举合并序列的大小
		complex Wn(cos(Pi / mid), op * sin(Pi / mid));//单位根
		for (int R = mid << 1, j = 0; j < limit; j += R) {//R是区间右端点，j表示左端点
			complex w(1, 0);//幂
			for (int k = 0; k < mid; k++, w = w * Wn) {//枚举区间的每个位置
				complex x = now[j + k], y = w * now[j + mid + k];//先求出来减少时间
				now[j + k] = x + y;//得到一个求出另外一个的
				now[j + mid + k] = x - y;
			}
		}
	}
}

int main() {
	scanf("%d %d", &n, &m);
	for (int i = 0; i <= n; i++) scanf("%lf", &a[i].x);
	for (int i = 0; i <= m; i++) scanf("%lf", &b[i].x);
	
	limit = 1;
	while (limit <= n + m) {
		limit <<= 1;
		l_size++;
	}
	
	for (int i = 0; i < limit; i++)
		an[i] = (an[i >> 1] >> 1) | ((i & 1) << (l_size - 1));
	
	FFT(a, 1);
	FFT(b, 1);
	
	for (int i = 0; i <= limit; i++)
		a[i] = a[i] * b[i];
	
	//将点值表示法转换为系数表示法再输出
	FFT(a, -1);
	
	for (int i = 0; i <= n + m; i++)
		printf("%d ", (int)(a[i].x / limit + 0.5));
	
	return 0;
}

__EOF__

本文作者：あおいSakura
本文链接：https://www.cnblogs.com/Sakura-TJH/p/luogu_P3803_FFT.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！