【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;
}
