多项式乘法（FFT）
Fiendish
·
·
算法·理论
复数
引入
我们比较熟悉的是实数，比如 1,1.5,\sqrt 3,\frac{13}{14} 都是实数。它们都可以表示在数轴上。
但是你会发现，\sqrt{-1} 似乎和它们不太一样。你无法在数轴上找到它。
我们引入一个新数 \text{i}=\sqrt{-1}，\text{i} 也被称为 虚数单位。我们定义形如 a+b\text{i}(a,b\in\R) 的数为复数。全体复数的集合称为 复数集，记作 \mathbb C。
代数形式
复数通常用 z 表示，即 z=a+b\text{i}，这是复数的代数形式。其中 a 称为复数 z 的 实部，b 称为复数 z 的 虚部。无特殊说明，一般认为 a,b\in\R。
几何意义
首先，我们定义两个复数 z_1=a+b\text{i},z_2=c+d\text{i} 相等当且仅当 a=c,b=d。这个定义十分自然。
考虑一个平面直角坐标系。在这个平面直角坐标系中，每个复数 z=a+b\text{i} 唯一对应一个点 Z(a,b)。这个平面直角坐标系被称为 复平面。其 x 轴对应复数的实部，被称为 实轴；y 轴对应复数的虚部，被称为 虚轴。
进一步地，点 Z(a,b) 也可以表示向量 \overrightarrow{OZ}，所以 a+b\text{i} 也对应着向量 (a,b)。
于是我们可以将向量的知识迁移到复数上来。我们定义复数的模即为其对应的向量的模，也即 |a+b\text{i}|=\sqrt{a^2+b^2}。
由于向量无法比较大小，我们可以发现复数也是无法比较大小的。
复数的加减乘除
我们定义两个复数 z_1=a+b\text{i},z_2=c+d\text{i}。
加法和减法
减法同理。$z_1-z_2=(a-c)+(b-d)\text{i}$。可以对应到向量减法。即 $(a,b)-(c,d)=(a-c,b-d)$。
容易发现复数的加法满足结合律和交换律，而复数的减法是复数加法的逆运算。
#### 乘法和除法
乘法：
$$\begin{align*}
z_1z_2=& (a+b\text{i})(c+d\text{i}) & \\
=& ac+ad\text{i}+bc\text{i}+bd\text{i}^2 & \\
=& (ac-bd)+(ad+bc)\text{i} &
\end{align*}$$
除法：
$$\begin{align*}
\frac{z_1}{z_2}=& \frac{a+b\text{i}}{c+d\text{i}} & \\
=& \frac{(a+b\text{i})(c-d\text{i})}{(c+d\text{i})(c-d\text{i})} & \\
=& \frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}\text{i} &(z_2\neq0)
\end{align*}$$
可以发现，复数的乘法满足交换律、结合律和乘法分配律。复数的除法是复数乘法的逆运算。
### 辐角和辐角主值
复数 $z$ 也可以用极坐标 $(r,\theta)$ 来表示，其中 $r$ 为 $z$ 的模，$\theta$ 是实轴正向到 **非零** 复数 $z=a+b\text{i}$ 的夹角。
$\theta$ 满足：
$$\tan\theta=\frac{y}{x}$$
称为 $z$ 的辐角。记为：
$$\theta=\arg z$$
任意一个 **非零** 复数 $z$ 都有无穷多个辐角，因此 $\arg z$ 事实上是一个集合。令 $\text{Arg }z$ 表示其中一个特定值，满足 $-\pi<\text{Arg }z\le\pi$。容易发现对于一个 **非零** 复数 $z$，这样的值是唯一的。我们称 $\text{Arg }z$ 为 $z$ 的 **辐角主值** 或 **主辐角**。而 $z$ 的所有辐角都可以在 **辐角主值** 的基础上加若干个（可能是任意整数个） $2\pi$ 得来。也即
$$\arg z=\{\text{Arg }z+2k\pi|k\in\Z\}$$
这时我们来观察复数乘法和复数除法的几何意义。可以发现，复数的乘法的几何意义为：两复数相乘，模长相乘，辐角相加。除法类似：两复数相除，模长相除，辐角相减。
### 单位圆和单位圆周
我们称复平面上模小于 $1$ 的所有复数构成的几何图形为 **单位圆**，而模等于 $1$ 的复数构成的几何图形为 **单位圆周**。在不引起混淆的情况下，有时单位圆周也简称单位圆。
### 单位根
称 $x^n=1$ 在复数意义下的解是 $n$ 次复根。显然这样的解有 $n$ 个，称这 $n$ 个解都是 $n$ 次 **单位根** 或 $n$ 次 **单位复根**。这 $n$ 个复数将单位圆周平分。
设 $\omega_n$ 为 $n$ 次单位根中辐角主值为正且辐角主值最小的一个，那么显然 $\omega_n^0,\omega_n^1,\dots,\omega_n^{n-1}$ 取遍了这 $n$ 个单位根（每乘上一个 $\omega_n$ 相当于逆时针旋转了 $\frac{2\pi}{n}$ 且模长还是 $1$ 没变）。
那我们如何计算 $\omega_n^k$ 的值呢？我们可以借助欧拉公式。
根据欧拉公式，我们有
$$\omega_n^k=\cos\frac{2k\pi}{n}+\text{i}\sin\frac{2k\pi}{n}$$
根据这个公式，我们来推导一些性质。
$$\omega_{2n}^{2k}=\cos\frac{4k\pi}{2n}+\text{i}\sin\frac{4k\pi}{2n}=\cos\frac{2k\pi}{n}+\text{i}\sin\frac{2k\pi}{n}=\omega_n^k$$
$$\omega_n^{k+\frac{n}{2}}=\omega_n^k\omega_n^{\frac{n}{2}}=\omega_n^k(\cos\pi+\text{i}\sin\pi)$$
由 $\text{e}^{\text{i}\pi}=\cos \pi+\text{i}\sin \pi=-1$ 可得
$$\omega_n^{k+\frac{n}{2}}=\omega_n^k(\cos\pi+\text{i}\sin\pi)=-\omega_n^k$$
## 快速傅里叶变换
### 引入
进入正题。
给定一个 $n$ 次多项式 $A$ 和一个 $m$ 次多项式 $B$，求多项式 $C=A\times B$。
直接做时间复杂度是 $O(nm)$ 的，非常不优。
考虑一个比较奇葩的思路。对于一个 $n-1$ 次多项式 $f(x)=\sum\limits_{i=0}^{n-1}a_ix^i$，我们取 $n$ 个互不相同的值代入 $f$ 中可以获得 $n$ 个数，这对应着平面直角坐标系中的 $n$ 个点，可以证明这 $n$ 个点唯一对应一个 $n-1$ 次多项式。我们称这样表示一个多项式的方式为点值表示法。
那么我们取 $n+m+1$ 个互不相同的数，称为 $x_0,x_1,\dots,x_{n+m}$，将这 $n+m+1$ 个数代入 $A$ 和 $B$ 中，可以得到共 $2(n+m+1)$ 个点。显然 $C(x_k)=A(x_k)B(x_k)$。也就是说，如果我们能快速地求出 $A$ 的点值表示法和 $B$ 的点值表示法，就可以 $O(n+m)$ 地求出 $C$ 的点值表示法，然后如果我们能把 $C$ 的点值表示法再快速地转换回去，那么我们就能快速求出 $C$ 了！
事实证明，点值表示法和系数表示法之间的互相转换可以做到 $O(n\log n)$ 的复杂度。接下来我们详细讲述一下如何做到这一点。
### 快速傅里叶变换
首先，系数表示法转换成点值表示法。我们考虑一个 $n-1$ 次多项式 $f(x)=\sum\limits_{i=0}^{n-1}a_ix^i$。
我们把 $n$ 变为大于等于 $n$ 的最小的 $2$ 的整数次幂，多出的项系数赋成 $0$ 就好。
我们可以取 $\omega_n^0,\omega_n^1,\omega_n^2,\dots,\omega_n^{n-1}$，将这 $n$ 个数代入 $f(x)$ 中，获得 $f(x)$ 的点值表示法。
将 $f(x)$ 的系数按下标奇偶性分类，再写出两个多项式
$$f_1(x)=a_0+a_2x+\dots+a_{n-2}x^{\frac{n}{2}-1}$$
$$f_2(x)=a_1+a_3x+\dots+a_{n-1}x^{\frac{n}{2}-1}$$
然后容易发现
$$f(x)=f_1(x^2)+xf_2(x^2)$$
将 $\omega_n^k(k<\frac{n}{2})$ 代入
$$f(\omega_n^k)=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})=f_1(\omega_{\frac{n}{2}}^k)+\omega_n^kf_2(\omega_{\frac{n}{2}}^k)$$
将 $\omega_n^{k+\frac{n}{2}}(k<\frac{n}{2})$ 代入
$$f(\omega_n^{k+\frac{n}{2}})=f_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}f_2(\omega_n^{2k+n})=f_1(\omega_{\frac{n}{2}}^k\omega_n^n)-\omega_n^kf_2(\omega_{\frac{n}{2}}^k\omega_n^n)=f_1(\omega_{\frac{n}{2}}^k)-\omega_n^kf_2(\omega_{\frac{n}{2}}^k)$$
我们发现，两个式子只有一个常数项不同。
所以只要我们只要算出第一个式子的值，就可以 $O(1)$ 求出第二个式子的值。所以问题规模变为原来的 $\frac{1}{2}$。于是我们递归下去做即可，时间复杂度 $O(n\log n)$。
### 快速傅里叶逆变换
其次，点值表示法转换成系数表示法。
假设我们现在已经知道了 $f(x)=\sum\limits_{i=0}^{n-1}a_ix^i$ 的点值表示法 $y=(y_0,y_1,\dots,y_{n-1})$，其中 $y_k=f(\omega_n^k)$。
我们先求出 $g(x)=\sum\limits_{i=0}^{n-1}y_ix^i$ 的点值表示法 $c=(c_0,c_1,\dots,c_{n-1})$，其中 $c_k=g(\omega_n^{-k})$。
推一下式子。
$$\begin{align*}
c_k= & \sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i \\
=& \sum\limits_{i=0}^{n-1}(\omega_n^{-k})^i\sum\limits_{j=0}^{n-1}a_j(\omega_n^i)^j \\
=& \sum\limits_{i=0}^{n-1}(\omega_n^{-k})^i\sum\limits_{j=0}^{n-1}a_j(\omega_n^j)^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
\end{align*}$$
根据等差数列求和公式，我们知道
$$\sum\limits_{i=0}^{n-1}x^i=\frac{1-x^n}{1-x}(x\neq 1)$$
当 $x=1$ 时
$$\sum\limits_{i=0}^{n-1}x^i=n$$
那么
$$\sum\limits_{i=0}^{n-1}(\omega_n^k)^i=\frac{1-(\omega_n^k)^n}{1-\omega_n^k}=\frac{1-(\omega_n^n)^k}{1-\omega_n^k}=0(0<k<n)$$
特别地
$$\sum\limits_{i=0}^{n-1}(\omega_n^0)^i=n$$
也就是说，当且仅当 $j=k$ 时，
$$\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i=n$$
否则等于 $0$。
也就是说，$c_k=na_k$，$a_k=\frac{c_k}{n}$。
这样我们就可以以 $O(n\log n)$ 的复杂度求出两个多项式的乘积了。
### 实现
#### 递归实现
我们发现，不论是系数 $\to$ 点值，还是点值 $\to$ 系数，我们真正需要的都是系数 $\to$ 点值。于是我们只需要一个函数。可以递归实现。
```cpp
#include<bits/stdc++.h>
using namespace std;
const int N=1e7;
const double Pi=acos(-1.0);
struct com{//复数
double x,y;
com(double X=0,double Y=0){
x=X;
y=Y;
}
}a[N],b[N];
com operator +(com a,com b){
return com(a.x+b.x,a.y+b.y);
}
com operator -(com a,com b){
return com(a.x-b.x,a.y-b.y);
}
com operator *(com a,com b){
return com(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
void fft(int lim,com *a,int t){//t=1是系数变点值，t=-1是点值变系数
if(lim==1) return;
com a1[lim>>1],a2[lim>>1];
for(int i=0;i<lim;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];//按下标奇偶性分类
fft(lim>>1,a1,t);
fft(lim>>1,a2,t);
com Wn=com(cos(2.0*Pi/lim),sin(2.0*t*Pi/lim)),w=com(1,0);
for(int i=0;i<(lim>>1);i++,w=w*Wn) a[i]=a1[i]+w*a2[i],a[i+(lim>>1)]=a1[i]-w*a2[i];
}
int n,m;
int main(){
ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
cin>>n>>m;
for(int i=0;i<=n;i++) cin>>a[i].x;
for(int i=0;i<=m;i++) cin>>b[i].x;
int lim=1;
while(lim<=n+m) lim<<=1;
fft(lim,a,1);
fft(lim,b,1);
for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
fft(lim,a,-1);
for(int i=0;i<=n+m;i++) cout<<(int)(a[i].x/lim+0.5)<<' ';//最后系数要除以n
}
```
#### 倍增实现
我们还有其他的实现方法。观察多项式系数再递归过程中位置发生的变化。我们以一个 $7$ 次多项式为例：
1. $$\{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\}$$
1. $$\{a_0,a_2,a_4,a_6\}\{a_1,a_3,a_5,a_7\}$$
1. $$\{a_0,a_4\}\{a_2,a_6\}\{a_1,a_5\}\{a_3,a_7\}$$
1. $$\{a_0\}\{a_4\}\{a_2\}\{a_6\}\{a_1\}\{a_5\}\{a_3\}\{a_7\}$$
系数的位置变化有没有规律呢？我们以 $a_1$ 为例，最开始其下标为 $1$，最终为 $4$。而我们将 $1$ 二进制表示翻转之后就会得到 $4$（`001` $\to$ `100`）。这个结论是正确的。再选一个观察，$a_2$ 的位置没变，是因为 `010` 翻转之后还是 `010`，没有变化。
当然我们要注意，我们需要先把二进制的前导零补上该结论才正确。
我们称这个变换为位逆序置换，这个可以以 $O(n)$ 的时间复杂度完成。
#### 蝶形运算优化
假设现在我们已经算出了 $f_1(\omega_{\frac{n}{2}}^k)$ 和 $f_2(\omega_{\frac{n}{2}}^k)$ 后，我们要求出 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac{n}{2}})$。
在使用位逆序置换后：
- $f_1(\omega_{\frac{n}{2}}^k)$ 的值存储在数组下标为 $k$ 的位置，而 $f_2(\omega_{\frac{n}{2}}^k)$ 存储在数组下标为 $k+\frac{n}{2}$ 的位置；
- $f(\omega_n^k)$ 的值将被存储在数组下标为 $k$ 的位置，而 $f(\omega_n^{k+\frac{n}{2}})$ 将被存储在数组下标为 $k+\frac{n}{2}$ 的位置。
也就是说，我们不会用到或者影响数组其他位置的值。
所以我们可以直接在数组下标为 $k$ 和 $k+\frac{n}{2}$ 的位置修改，不需要额外开数组存，此方法被称为 **蝶形运算优化**。
```cpp
#include<bits/stdc++.h>
using namespace std;
struct com{//复数
double x,y;
com(double X=0,double Y=0){
x=X,y=Y;
}
};
com operator +(com &a,com &b){
return com(a.x+b.x,a.y+b.y);
}
com operator -(com &a,com &b){
return com(a.x-b.x,a.y-b.y);
}
com operator *(com &a,com &b){
return com(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
const int N=1e7+10;
const double PI=acos(-1);
int r[N];
void fft(int n,com *A,int t){
for(int i=0;i<n;i++) if(i<r[i]) swap(A[i],A[r[i]]);//位逆序置换
for(int mid=1;mid<n;mid<<=1){
com W(cos(PI/mid),t*sin(PI/mid));//单位根
for(int R=mid<<1,now=0;now<n;now+=R){
com w(1,0);
for(int k=0;k<mid;k++){
//mid是区间长度的一半，R是区间长度，now是当前区间的左端点，k是当前区间中的一个位置
com x=A[now+k],y=w*A[now+k+mid];//蝶形运算优化
A[now+k]=x+y;
A[now+k+mid]=x-y;
w=w*W;
}
}
}
}
int n,m,lim,l;
com a[N],b[N],c[N];
int main(){
cin>>n>>m;
for(int i=0;i<=n;i++) cin>>a[i].x;
for(int i=0;i<=m;i++) cin>>b[i].x;
lim=1;
while(lim<=n+m) lim<<=1,l++;
for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|(((i&1)<<(l-1)));//预处理二进制的翻转
fft(lim,a,1),fft(lim,b,1);
for(int i=0;i<lim;i++) c[i]=a[i]*b[i];
fft(lim,c,-1);
for(int i=0;i<=n+m;i++) cout<<(int)(c[i].x/lim+0.5)<<' ';
}
```
## 例题
[【模板】多项式乘法（FFT）](https://www.luogu.com.cn/problem/P3803)
模板题，没什么可说的。
[【模板】高精度乘法 | A*B Problem 升级版](https://www.luogu.com.cn/problem/P1919)
似乎是高精度乘法。
但是我们发现，高精度复杂度就炸了。
考虑优化高精度。可以发现，一个整数可以表示成一个多项式。比如说
$$A=\sum\limits_{i=0}^{\log_{10} A}\text{num}(i)\times10^i$$
其中 $\text{num}(i)$ 表示 $A$ 在十进制下从高往低数第 $i$ 位的数字（最低位是第 $0$ 位）。
于是，两个整数相乘就相当于两个多项式相乘。可以用 FFT 解决。最后要处理进位。
```cpp
#include<bits/stdc++.h>
using namespace std;
const int N=5000010;
const double Pi=acos(-1.0);
struct com{
double x,y;
com(double X=0,double Y=0){
x=X;
y=Y;
}
}a[N],b[N];
com operator +(com a,com b){
return com(a.x+b.x,a.y+b.y);
}
com operator -(com a,com b){
return com(a.x-b.x,a.y-b.y);
}
com operator *(com a,com b){
return com(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
void fft(int lim,com *a,int t){
if(lim==1) return;
com a1[lim>>1],a2[lim>>1];
for(int i=0;i<lim;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fft(lim>>1,a1,t);
fft(lim>>1,a2,t);
com Wn=com(cos(2.0*Pi/lim),sin(2.0*t*Pi/lim)),w=com(1,0);
for(int i=0;i<(lim>>1);i++,w=w*Wn) a[i]=a1[i]+w*a2[i],a[i+(lim>>1)]=a1[i]-w*a2[i];
}
int n,m;
string s1,s2;
int c[N];
int main(){
cin>>s1>>s2;
reverse(s1.begin(),s1.end());
reverse(s2.begin(),s2.end());
for(auto i:s1) a[n++].x=i-'0';
for(auto i:s2) b[m++].x=i-'0';
int lim=1;
while(lim<=n+m) lim<<=1;
fft(lim,a,1);
fft(lim,b,1);
for(int i=0;i<=lim;i++) a[i]=a[i]*b[i];
fft(lim,a,-1);
for(int i=0;i<lim;i++) c[i]=(int)(a[i].x/lim+0.5);
for(int i=0;i<lim;i++){
if(c[i]>=10){
c[i+1]+=c[i]/10;
c[i]%=10;
lim++;
}
}
while(lim>2&&!c[lim-1]) lim--;
for(int i=lim-1;i>=0;i--) cout<<c[i];
return 0;
}
```