【模板】多项式乘法（FFT）
题目链接：luogu P3803
题目大意
给你两个多项式，要你求这两个多项式乘起来得到的多项式。（卷积）
思路
系数表示法
就是我们一般来表示一个多项式的方法：
(A(x)=a_1x^k+a_2x^{k-1}+...+a_k)
或者可以这样表示：
(A(x)=sumlimits_{i=1}^{k}a_i imes x_i)
那你很容易看到，用来做这道题用系数表示法来做是 (O(n^2)) 的。
点值表示法
假设我们已经知道了这个多项式，那我们把 (n) 个不同的 (x) 带入，会得到 (n) 个取值 (y)。
由于 (n) 个点可以确定一个 (n) 次多项式，那我们就可以根据这 (n) 个点确定这个多项式。
那你如果选了 (n) 个点 ((x_1,y_1),...,(x_n,y_n))，那就有：
(y_i=sumlimits_{j=0}^{n-1}a_j imes 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 kdfrac{2pi}{n}+isin kdfrac{2pi}{n})
（大概就是把模长那条线往 (x) 轴做垂线得到直角三角形，然后用三角函数得到两条直角边的长度，从而得到坐标，也就是复数）
当然显而易见，单位根幅角就是 (frac{1}{n})
还有就是如果 (z^n=1)，那 (z) 就是 (n) 次单位根。
单位根的一些性质
(omega_{2n}^{2k}=omega_n^k)
证明：
(omega_{2n}^{2k}=cos 2kdfrac{2pi}{2n}+isin 2kdfrac{2pi}{2n})
(=cos kdfrac{2pi}{n}+isin kdfrac{2pi}{n}=omega_n^k)
(omega_{n}^{k+frac{n}{2}}=-omega_{n}^{k})
证明：
(omega_{n}^{frac{n}{2}}=cos dfrac{n}{2}*dfrac{2pi}{n}+isin dfrac{n}{2}*dfrac{2pi}{n})
(=cos pi+isinpi=-1)
（(cos pi=cos 180^circ=-1,sinpi=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=sumlimits_{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=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^i)^j)(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i)(omega_n^{-k})^i)
(=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i))
(=sumlimits_{i=0}^{n-1}sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i)
(=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{j-k})^i)
我们设 (S(x)=sumlimits_{i=0}^{n-1}x^i)。
带入 (omega_n^k)：(S(omega_n^k)=sumlimits_{i=0}^{n-1}(omega_n^k)^i)
然后分类讨论：
当 (k=0) 时，(omega_n^k=omega_n^0=1)。
那式子就变成了 (S(omega_n^k)=sumlimits_{i=0}^{n-1}1^i=n)
当 (k
eq0) 时，我们考虑怎么求。
观察到没项都乘了一个值，我们考虑用这么一种方法：
乘 (omega_n^k)：(omega_n^kS(omega_n^k)=sumlimits_{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=sumlimits_{j=0}^{n-1}a_jsumlimits_{i=0}^{n-1}(omega_n^{j-k})^i)
(=sumlimits_{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;
}