前置知识
FFT（Fast Fourier Transformation），中文名快速傅里叶变换，OI中用来加速多项式乘法。实际上是离散傅里叶变换 (Discrete,Fourier Transformation，DFT) 的计算机快速计算方法的统称
朴素高精度乘法时间 \(O(n^2)\)，但 FFT 能 \(O(n\log n)\) 内解决
点值表达式
对于多项式 \(A(x)=a_0+a_1x+a_2x^2+a_3x^3\dots a_nx^n\)。
任意 \(n+1\) 个不同点 \((x_1,A(x_1)),(x_2,A(x_2)),\dots ,(x_{n+1},A(x_{n+1}))\) 可以唯一确定一组 \(\{a_i\}\) 。
不严谨验证：
我们将所有的点代入，试图求 \(\{a_i\}\) ，得到一个如下图的方程组：
它的系数矩阵行列式显然是一个转置的 \(n+1\) 阶 \(Vandermonde\) 行列式：
由于点之间两两不同，所以行列式值一定不为 \(0\)，则原矩阵满秩，故原方程组有唯一解。(对复数域成立)
那么这个东西方便在哪里呢？
对于两个多项式，我们各取 \(n+m+1\) 个点 \(x_1,x_2,\dots ,x_{n+m+1}\)。
得到 \((x_1,A(x_1))\ ,\ (x_2,A(x_2))\ ,\ (x_3,A(x_3))\ \dots\ (x_{n+m+1},A(x_{n+m+1}))\ ,\ (x_1,B(x_1))\ ,\ \dots (x_{n+m+1},B(x_{n+m+1}))\)。
我们将对应点相乘：\((x_1,A(x_1)B(x_1))\ ,\ (x_2,A(x_2)B(x_2))\ ,\ (x_3,A(x_3)B(x_3))\ ,\ \dots\ ,\ (x_{n+m+1}\ ,\ A(x_{n+m+1})B(x_{n+m+1}))\)。
发现这些点唯一确定了一个 \(C(x)=A(x)B(x)\) 。
使用系数表示法的时候计算乘积的复杂度是 \(\Theta(nm)\) ，而上面的方法是 \(\Theta (n+m)\) 。
但是我们一般使用的都是系数表示，还要进行点值与系数表示的互相转换。
朴素的系数转点值算法叫做离散傅里叶变换 (DFT)。逆向变换叫做离散傅里叶逆变换 (IDFT)。
复数
讲 DFT 之前先提亿下复数。
复数定义
称形如 \(a+bi,a,b\in R\) 的数，其中 \(i\) 是虚数单位，满足 \(i^2=-1\)。
复平面
复数平面即是 \(z=a+bi\) ,它对应的坐标为 \((a,b)\)。其中，\(a\) 表示的是复平面内的横坐标，\(b\) 表示的是复平面内的纵坐标，表示实数 \(a\) 的点都在 \(x\) 轴上，所以 \(x\) 轴又称为“实轴”；表示纯虚数 \(bi\) 的点都在 \(y\) 轴上，所以 \(y\) 轴又称为“虚轴”。\(y\) 轴上有且仅有一个实点即为原点 "\(0\)"。
一个复数可以唯一对应在复平面上的一个点或向量。
复数的模
对于复数 \(z=a+bi\) ，它的模 \(|\ z\ |=\sqrt{a^2+b^2}\)。
共轭复数
对于复数 \(z=a+bi\) ，我们记 \(\overline{z}=a-bi\) 为 \(z\) 的共轭复数。
在复平面上，表示两个共轭复数的点或向量关于 \(x\) 轴对称。
共轭复数的性质：
- \(|a+bi|=|a-bi|\)
- \((a+bi)\times (a-bi)=a^2+b^2\)。
复数运算
-
加法：
代数形式：\((a+bi)+(c+di)=(a+c)+(b+d)i\)
-
乘法
代数形式：\((a+bi)\times(c+di)=(ac-bd)+(bc+ad)i\)
几何意义："模长相乘，幅角相加"。
-
除法
代数形式：\(\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{(a+bi)(c-di)}{c^2+d^2}\) (即分母实数化)
单位根
在复平面上以原点为中心，以 \(1\) 为半径作一个圆，这个圆即是单位圆。在圆内接一个正 \(n\) 边形，其中一个顶点对应的复数为 \(1\) ，则这些顶点对应的复数就是 \(n\) 次单位根。
\(n\) 次单位根的严格定义为 \(z^n=1\) 的复数解。
从这里开始，下文单位根的次数 \(n\) 均为 \(2\) 的整数次幂。
不妨设 \(\omega_n\) 为这些复数解中幅角为正且最小的复数解。则任意一个 \(n\) 次单位根可以表示为 \(\omega_n^k\)。
特殊的，\(\omega_n^0=\omega_n^n=1\)。
那么如何得到其他单位根的值呢？
先给公式： \(\omega_n^k=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}\)
图形推导：这些单位根在圆上对应的点将圆 \(n\) 等分了，考虑最小的复数解 \(\omega_n\) ，画图：
其中 \(\theta=\frac{2\pi}n\) ，用三角函数表示 \(\omega_n=\cos \theta+i\sin \theta\)。
所以根据复数乘法的几何意义 \(\omega_n^k=\cos k\theta+i\sin k\theta\)
根据欧拉公式简单证明：
同时易得单位根幅角是圆周角的 \(\frac1n\)。
单位根还有一些性质：
-
\(\omega_n^i\ne \omega_n^j\)
-
\(\omega_{2n}^{2k}=\omega^k_n\)
-
\(\omega_{n}^k =- \omega_n^{k+\frac n2}\)
证明均只需要代入公式即可。
理论部分
快速傅里叶变换
接下来的思路参考自bilibili 快速傅里叶变换(FFT)——有史以来最巧妙的算法？
考虑将下面这个多项式转成点值表达式
首先想到暴力解决：我们随便带入 \(n\) 个互不相同的初值，直接计算对应的结果。这样的时间复杂度是 \(O(n^2)\) 的，显然不太优秀。
那么如果我们不随便找代入的点，而是去找一些特殊点呢？
比如在高中数学中，我们总是从对称点对 \(x\) 和 \(-x\) 考虑关于函数的一些性质：若一个函数 \(f(x)\) 是偶函数则 \(f(x)=f(-x)\) ；如果是奇函数则 \(f(-x)+f(x)=0\)。
那么考虑将一组相反数 \(x_1\) 与 \(x_2=-x_1\) 代入 \(A(x)\) 得到：
观察两个式子的差别，我们将多项式中的奇偶项分类：
分别设：
于是
也就是说，若是要求 \(A(x)\) 在 \(\pm x_1,\pm x_2,\ \dots\ \pm x_{n/2}\) 位置上的取值，我们需要求 \(A_1(x)\ ,\ A_2(x)\) 在 \(x_1^2,x_2^2,\ \dots\ x_{n/2}^2\) 处的取值。
一个 \(n-1\) 阶多项式 \(n\) 点取值问题被转化为了两个 \(\frac n2-1\) 阶多项式在 \(\frac n2\) 个点的取值问题，我们会联想到递归。
但是真的能递归吗？
细想一下我们会发现，\(A_1(x_1^2)\) 和 \(A_2(x_1^2)\) 的取值点都必须是正数，但是上面缩小数据范围的方式的关键在于取互为相反数的两个点，而如果我们只能取正点的话，求解单个点时，我们就只能暴力硬算，单次复杂度 \(O(n)\) 。总的时间复杂度仍然为 \(O(n^2)\)，只是少了点常数罢了。
继续往下细想：为什么不能递归？因为有平方的存在，我们不能取到负数点了。
那么什么样的数能够让我们在平方之后还能取到相反值呢？自然的，我们会想到复数。当然单纯 “是复数” 这一个条件是不够的，我们还要找出什么样的复数才能使得递归正常进行。
在递归的过程中，我们需要将两个相反点 \(x_1,-x_1\) 相互配对，再往下递归。递归到下一层 \(x_1^2\) 之后，我们需要另外一组 \(x_2,-x_2\) 递归到 \(x_2^2=-x_1^2\) ，这两个数配对在向下递归得到 \(x_1^4\) ……
当它最终递归到只剩一项时，会得到 \(x_1^n\)。我们不妨设 \(x_1^n=1\) ，那么这就是单位根的定义！也就是说，我们最开始取点的时候，取 \(n\) 阶单位根即可！
现在我们来看代入单位根之后递归过程有什么变化。
我们要求 \(A(\omega_n^k)\)，即求 \(A_1(\omega_{n/2}^k)+\omega_n^kA_2(\omega_{n/2}^k)\)
它对应的相反项：\(A(\omega_n^{k+n/2})=A_1(\omega_{n/2}^k)-\omega_n^kA_2(\omega_{n/2}^k)\)
也就是说我们每次将当前层数对应的点数 \(n\) 折成两半 \([0,\frac n2),[\frac n2,n]\)，其中 \(w_n^k=w_n^{k+n/2},k\in [0,\frac n2)\) 且为整数。
快速傅里叶变换的逆变换
设 \((y_0,y_1,y_2,\dots,y_{n-1})\) 是多项式 \((a_0,a_1,a_2,\dots, a_{n-1})\) 在 \((\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1})\) 的点值向量，有一个向量 \((c_0,c_1,\dots,c_{n-1})\) 满足
则有 \(a_k=\frac{c_k}n\)
小小证明一下：
设 \(S(x)=\sum_{i=0}^{n-1}x^i\)
代入 \(\omega_n^{k}\) ：
当 \(k\ne 0\) 时，
当 \(k=0\) 时，\(\omega_n^0=1,\ S(1)=n\)。
我们回到刚才的式子：
当 \(k\ne j\) 时，该项为 \(0\)；当 \(k=j\) 时该项为 \(na_k\)。
所以 \(c_k=na_k\)。
那么如何求 \(c_k\) 呢？
构造函数 \(B(x)=y_0+y_1x+y_2x^2+\cdots+y_{n-1}x^{n-1}\)。我们发现 \(c_k\) 就是这个多项式在 \(n\) 次单位根上的点值，仍然可以利用 FFT 的过程来做。
实现部分
首先我们可以想到一个按照流程的递归实现。但是很容易知道这样做的常数是巨大的。
所以我们一般用迭代来实现 FFT 。
首先我们观察递归的过程：
对于一个系数向量
首先我们会把他们奇偶分项，可以得到两个系数向量：
继续分直至每个向量内只有单个元素：
这构成了一棵树的关系，我们首先需要将最底层的顺序处理出来，然后计算最底层的值，并一层层归并上去。
那么如何得到最底层的顺序？
先给出答案：二进制翻转。
例如 \(a_3\) 的下标为 \(3=11_{(2)}\) ，我们将前缀 \(0\) 补齐到总位数与 \(n-1\) 相同，即 \(011\) ，然后翻转，得到 \(110\) ，那么 \(a_3\) 所在的下标就应该是 \(6\)。
怎么来的？我们考虑分组的过程。
对于第一次分组，奇偶分项的过程就是看 \(a_k\) 的下标 \(k\) 二进制最低位是否为 \(1\) 。当这一位为 \(0\) ，我们会把它放到前半段，此时它在数组中的下标二进制最高位就是 \(0\)；当这一位为 \(1\)，那么分组后它在数组中的位置下标最高位就位 \(1\)。
对于第二次分组，我们只看前半个数组，这里面的所有 \(a_k\)，\(k\) 的最低位都是 \(0\)。那么我们对它们的奇偶分组就是在看它的次低位是否为 \(1\)。
对于第三次分组，我们仍然是根据它的第三低位是否为 \(1\) 分组。
以此类推，到最后对于某个 \(a_k\) 所在位置的下标 \(i\)，\(k\) 的最高位对应 \(i\) 的最低位，\(k\) 的次高位对应 \(i\) 的次低位 …… \(k\) 的最低位对应 \(i\) 的最高位。
也就是说，\(i\) 就是 \(k\) 的二进制翻转。
全过程代码：
#include <bits/stdc++.h>
using namespace std;
const int N=5e5+10;
const double PI=acos(-1),eps=1e-7;
struct Comp//复数类
{
double x,y;//实部虚部
Comp operator + (const Comp &t) const
{
return {x+t.x,y+t.y};
}
Comp operator - (const Comp &t)const
{
return {x-t.x, y-t.y};
}
Comp operator * (const Comp &t)const
{
return {x*t.x-y*t.y,y*t.x+x*t.y};
}
} a[N],b[N];//两个多项式的系数
int n,m;
int rev[N],bit;//预处理顺序，最高位
int tot;
void FFT(Comp a[],int inv)
{
for(int i=0;i<tot;i++)
if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<tot;mid<<=1)
{
Comp w1=Comp({cos(PI/mid),inv*sin(PI/mid)});
for(int i=0;i<tot;i+=(mid<<1))
{
Comp wk=Comp({1,0});
for(int j=0;j<mid;j++,wk=wk*w1)
{
Comp x=a[i+j],y=wk*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=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);
while((1<<bit)<n+m+1) ++bit;
tot=1<<bit;
for(int i=0;i<tot;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
FFT(a,1), FFT(b,1);//FFT
for(int i=0;i<tot;i++)
a[i]=a[i]*b[i];//IFFT
FFT(a,-1);
for(int i=0;i<=n+m;i++)
printf("%d ",(int)(a[i].x/tot+0.5));
return 0;
}