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

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

题目链接：luogu P3803

题目大意

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

思路

系数表示法

就是我们一般来表示一个多项式的方法：
$A(x)=a_1{x}^k+a_2{x}^{k-1}+...+a_k$
或者可以这样表示：
$A(x)=\sum _{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 _{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)\ast (c+di)$
$=ac+a\ast di+b\ast ci+bi\ast 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\frac{2\pi }{n}+i\sin k\frac{2\pi }{n}$
（大概就是把模长那条线往 $x$ 轴做垂线得到直角三角形，然后用三角函数得到两条直角边的长度，从而得到坐标，也就是复数）

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

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

单位根的一些性质

${\omega }_{2n}^{2k}={\omega }_n^k$
证明：
${\omega }_{2n}^{2k}=\cos 2k\frac{2\pi }{2n}+i\sin 2k\frac{2\pi }{2n}$
$=\cos k\frac{2\pi }{n}+i\sin k\frac{2\pi }{n}={\omega }_n^k$

${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$
证明：
${\omega }_n^{\frac{n}{2}}=\cos \frac{n}{2}\ast \frac{2\pi }{n}+i\sin \frac{n}{2}\ast \frac{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\ast {\omega }_n^{\frac{n}{2}}={\omega }_n^k\ast (-1)=-{\omega }_n^k$

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

快速傅里叶变换（FFT）

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

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

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

把 ${\omega }_n^k$ 带入，$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_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}\ast {\omega }_n^n)-{\omega }_n^k{A}_2({\omega }_n^{2k}\ast {\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 _{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 _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i$
$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i$
$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i)({\omega }_n^{-k}{)}^i$
$=\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{a}_j({\omega }_n^j{)}^i({\omega }_n^{-k}{)}^i)$
$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^{j-k}{)}^i$
$=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i$

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

然后分类讨论：

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

当 $k\ne 0$ 时，我们考虑怎么求。
观察到没项都乘了一个值，我们考虑用这么一种方法：
乘 ${\omega }_n^k$：${\omega }_n^{k}S({\omega }_n^k)=\sum _{i=1}^n({\omega }_n^k{)}^i$
然后两个相减：${\omega }_n^{k}S({\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)=\frac{1-1}{{\omega }_n^k-1}=0$

OK，分析完了 $S(x)$，我们回到之前卡克的地方：
$c_k=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i$
$=\sum _{j=0}^{n-1}{a}_{j}S({\omega }_n^{j-k})$
那只会有一个 $j=k$，使得 $j-k=0$，贡献是 $n$，其它的时候，贡献都是 $0$。
那就是 $c_k=a_{k}n$
那你就会发现，$a_k=\frac{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;
}