多项式学习笔记
多项式乘法
本部分主要介绍 傅里叶变换 及 快速傅里叶变换（FFT）等在 OI范围内 的应用。
本文会简要介绍 FFT 的原理、思想，涉及诸如 复平面、平面向量 等数学内容，文中会简要介绍。
FFT 快速傅里叶变换
傅里叶变换的本质就是将一个多项式在 系数表示法 和 点值表示法 之间转换。
FFT 能优化乘法的本质是点值表示法下可以直接乘。那么将多项式 \(f(x),g(x)\) 转换为点值表示法后直接相乘，再转换回系数表示法即可求出答案。
FFT 可以在 \(\mathcal O\left(n\log n\right)\) 的时间复杂度内计算出两个 \(n\) 次多项式的乘积。同时，考虑大整数可视为关于 \(10\) 的多项式，因此 FFT 也可以用于优化高精度乘法。
前置知识
复数与复平面
我们对实数域进行扩域可以得到复数域，这是最大的已知域。
一个复数 \(z\) 的代数形式形如：
其中，\(\mathrm i=\sqrt{-1}\) 为虚数单位，称 \(a\) 是 \(z\) 的实部，记为 \(\operatorname{Re}(z)\)，\(b\) 是 \(z\) 的虚部，记为 \(\operatorname{Im}(z)\)，且 \(a,b\in\R\)。
由此，可以有复平面（类似于平面直角坐标系）。例如点 \((3,2)\) 在复平面上就对应 \(3+2\mathrm i\)。
复数集为 \(\mathbb C\)。
复数的运算律与实数相同。
几何意义
对于复数 \(z=a+b\mathrm i\)，对应的复数点 \(Z(a,b)\)，对应平面向量 \(\overrightarrow{OZ}=(a,b)\) 可以使用极坐标 \((r,\theta)\) 表示。\(r=\sqrt{a^2+b^2}\) 为向量的模，\(\theta\) 为向量夹角（实轴正方向到 \(z\) 对应向量的夹角），满足：
称 \(\theta\) 为 \(z\) 的幅角，记为：
此时，复数的乘除法变得简单：
- 乘法：模相乘，幅角相加。
- 除法：模相除，幅角相减。
欧拉公式
单位根
称单位圆为复平面内以 \((0,0)\) 为原点，\(1\) 为半径的圆。
有 \(n\) 次方程：
此方程的根有 \(n\) 个，称这 \(n\) 个根均为 单位根，则 \(n\) 次单位根等分了单位圆（因为模均为 \(1\)，幅角相加）。
设幅角为 \(\dfrac{2\pi}n\) 的单位复数：
则 \(\omega_n\) 是一个 \(n\) 次单位根。因为 \((\omega_n)^n=e^{2\pi\mathrm i}=1\)。
则 \(\omega_n^k\) 也是一个 \(n\) 次单位根，因为 \((\omega_n^k)^n=1^k=1=\cos\dfrac{2\pi k}{n}+\mathrm i\sin\dfrac{2\pi k}n\)。
单位根具有性质（\(n,k\in \N\)）：
代码实现
在 C++ 的 STL 中，提供了标准库 <complex> 作为复数基本类型。
使用方法：
-
定义系数类型为
T的复数x：complex<T>x; -
引用
x的实部：x.real();将
x的实部赋值为y：x.real(y); -
引用
x的虚部：x.imag();将
x的虚部赋值为y：x.imag(y); -
运算复数
x、y：x+y;x+=y; x-y;x-=y; x*y;x*=y;
实际上，STL 的 complex 类或许比你手写的 complex 类更快。这里提供一份手写 complex 类：
struct complex{
double a,b;
complex(double A=0,double B=0){
a=A,b=B;
}
}f[Limit+1],g[Limit+1];
complex operator +(complex x,complex y){
return {x.a+y.a,x.b+y.b};
}
complex operator -(complex x,complex y){
return {x.a-y.a,x.b-y.b};
}
complex operator *(complex x,complex y){
return {x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a};
}
complex operator *=(complex &x,complex y){
return x=x*y;
}
手写用时 \(\text{3.48s}\)，使用 STL 用时 \(\text{2.52s}\)。
DFT/IDFT 离散傅里叶（逆）变换
多项式的点值表示法
多项式的常见表示方法是系数表示法。设 \(n-1\) 次多项式 \(f(x)\)：
然而，多项式还可以通过点值表示。即确定 \(n\) 个点 \((x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\)，满足 \(f(x_i)=y_i\)，从而确定 \(f(x)\)。
DFT 将 \(n\) 个复数点 \(\omega_n^0,\omega_n,\omega_2^2,\cdots,\omega_n^{n-1}\) 代入 \(f(x)\)，从而得到了 \(f(x)\) 的点值表示法。
IDFT 离散傅里叶逆变换
IDFT 将 \(n-1\) 次多项式 \(f(x)\) 的 DFT 结果当作另一个 \(n-1\) 次多项式 \(g(x)\) 的系数，取 \(\omega_n^{-0},\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-(n-1)}\) 再做一次 DFT 即可还原原来系数表示法的 \(f(x)\)。这个过程被称为 IDFT（离散傅里叶逆变换）。
具体而言，若还原得到的结果为 \(a_0,a_1,\cdots,a_{n-1}\)，则实际系数为 \(\dfrac{\operatorname{Re}(a_0)}n,\dfrac{\operatorname{Re}(a_{1})}n,\cdots,\dfrac{\operatorname{Re}(a_{n-1})}n\)。
也正因此，要取单位根 \(\omega_n^k\) 代入。
FFT 快速傅里叶变换
考虑到，朴素的 DFT/IDFT 的时间复杂度为 \(\mathcal O\left(n^2\right)\)，这需要一些优化，因此有了 FFT。FFT 的本质是基于分治思想的 DFT。（事实上，FFT 应当被称为 FDFT，但由于只有 DFT 能够快速地求，因此 FFT 默认指 FDFT。）
考虑 \(n-1\) 次多项式 \(A(x)\)：
不妨令 \(n\) 为 \(2\) 的正整数次幂。若 \(n\) 不是，则可以将高位补全至 \(2\) 的正整数次幂，系数均为 \(0\) 即可。
考虑将 \(A(x)\) 按照奇偶性分类：
则有：
那么，我们分治处理出 \(A_0(x),A_1(x)\) 的点值表示法，随后快速合并答案。这样，因为 \(A_0(x),A_1(x)\) 的规模都缩小了一半，就可以快速求解。
将 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{\large\frac n2-1}\) 代入 \(A(x)\) 可得：
将 \(\omega_n^{\large\frac n2},\omega_n^{\large\frac n2+1},\cdots,\omega_n^{n-1}\) 代入 \(A(x)\) 可得：
可以发现，差别仅仅是一个正负号。
因此可以通过将 \(\omega_{\large\frac n2}^0,\omega_{\large\frac n2}^1,\cdots,\omega_{\large\frac n2}^{\large\frac n2-1}\) 代入 \(A_1(x),A_2(x)\)，于是就能够在原来一半规模的时间内合并出答案。
这样，FFT 的时间复杂度便被优化到了 \(\mathcal O\left(n\log n\right)\)。
实现细节
-
对于整系数多项式相乘，取答案系数时，应当取实部四舍五入后的结果。
因为 FFT 使用的
double具有浮点误差。
FFT 例题
P3803 【模板】多项式乘法（FFT）
模板题。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<complex>
using namespace std;
constexpr const int N=1e6,M=N,Limit=1<<(int)(ceil(log2(N+M)))|1;
constexpr const double Pi=acos(-1);
#define complex complex<double>//便于书写
complex f[Limit+1],g[Limit+1];
void FFT(complex a[],int n,int op){
if(n==1){
return;
}
complex a0[n>>1|1],a1[n>>1|1];
for(int i=0;i<(n>>1);i++){
a0[i]=a[i<<1];
a1[i]=a[i<<1|1];
}
FFT(a0,n>>1,op);FFT(a1,n>>1,op);
complex w={1,0},omega={cos(2*Pi/n),sin(2*Pi/n)*op};
for(int i=0;i<(n>>1);i++,w*=omega){
a[i]=a0[i]+w*a1[i];
a[i+(n>>1)]=a0[i]-w*a1[i];
}
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,m;
cin>>n>>m;
for(int i=0;i<=n;i++){
double pl;
cin>>pl;
f[i].real(pl);
}
for(int i=0;i<=m;i++){
double pl;
cin>>pl;
g[i].real(pl);
}
int limit=1;
while(limit<=n+m){
limit<<=1;
}
FFT(f,limit,1);
FFT(g,limit,1);
for(int i=0;i<limit;i++){
f[i]*=g[i];
}
FFT(f,limit,-1);
for(int i=0;i<=n+m;i++){
int pl=round(f[i].real()/limit);
cout<<pl<<' ';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
P1919 【模板】高精度乘法 | A*B Problem 升级版
倒序存储 \(a,b\) 即可。注意最后乘出来的结果可能会出现某一位大于 \(9\) 的情况，需要处理一遍进位。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<complex>
using namespace std;
constexpr const int N=1000000,M=N,Limit=1<<(int)(ceil(log2(N+M)))|1;
constexpr const double Pi=acos(-1);
#define complex complex<double>
complex f[Limit+1],g[Limit+1];
int ans[Limit+1];
void FFT(complex a[],int limit,int op){
if(limit==1){
return;
}
complex a0[limit>>1|1],a1[limit>>1|1];
for(int i=0;i<(limit>>1);i++){
a0[i]=a[i<<1];
a1[i]=a[i<<1|1];
}
FFT(a0,limit>>1,op);FFT(a1,limit>>1,op);
complex w={1,0},omega={cos(2*Pi/limit),sin(2*Pi/limit)*op};
for(int i=0;i<(limit>>1);i++,w*=omega){
a[i]=a0[i]+w*a1[i];
a[i+(limit>>1)]=a0[i]-w*a1[i];
}
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
string fStr,gStr;
cin>>fStr>>gStr;
reverse(fStr.begin(),fStr.end());
reverse(gStr.begin(),gStr.end());
int n=fStr.size()-1,m=gStr.size()-1;
for(int i=0;i<=n;i++){
f[i].real(fStr[i]-'0');
}
for(int i=0;i<=m;i++){
g[i].real(gStr[i]-'0');
}
int limit=1;
while(limit<=n+m){
limit<<=1;
}
FFT(f,limit,1);
FFT(g,limit,1);
for(int i=0;i<limit;i++){
f[i]*=g[i];
}
FFT(f,limit,-1);
for(int i=0;i<limit;i++){
ans[i]+=round(f[i].real()/limit);
ans[i+1]+=ans[i]/10;
ans[i]%=10;
}
while(ans[limit]==0&&limit-1>0){
limit--;
}
for(int i=limit;i>=0;i--){
cout<<ans[i];
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
NTT 快速数论变换
事实上，「NTT」是指「数论变换」，而「FNTT」才是「快速数论变换」。但是，在算法竞赛中常提到的 NTT 一词，往往实际指的是快速数论变换，一般默认「数论变换」是指「快速数论变换」。
数论变换，即 DFT 在数论基础上的实现。快速数论变换，即 FFT 在数论基础上的实现。
NTT 可以在模意义下 \(\mathcal O\left(n\log n\right)\) 计算多项式乘法，并且不需要 FFT 的复数，也不需要浮点数运算，从而使其精度更高、常数更小。
前置知识
原根与阶
设 \(a\perp p\)。若最小的正整数 \(n\) 满足 \(a^n\equiv1\pmod p\)，称 \(n\) 为 \(a\) 模 \(p\) 的阶，记为 \(n=\delta_p(a)\)。
若 \(\delta_p(a)\equiv\varphi(a)\pmod p\)，则称 \(a\) 是模 \(p\) 的原根。
若 \(p\) 为质数，则有其原根 \(g_p=\Omega(\log p)\)，因此可以暴力找。
原根可以在 NTT 中代替 FFT 中的单位根进行运算。
单位根满足：
而原根在模 \(p\) 意义下同样满足这些性质。
最常见的模数与其原根：
NTT 快速数论变换
NTT 其实很简单，设 \(g\) 是模数 \(p\) 的原根，将 \(\omega_n\) 替换为 \(g^{\large\frac{p-1}{n}}\) 即可。
注意在 NTT 中也需要将 \(n\) 补全至 \(2\) 的正整数次幂，这可以保证 \((p-1)\bmod n=0\)。
P3803 【模板】多项式乘法（FFT） AC 代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<complex>
using namespace std;
constexpr const int N=1e7,M=N,P=998244353,G=3;
int f[N+1],g[M+1];
int qpow(int base,int n){
n=(n%(P-1)+P-1)%(P-1);
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
void NTT(int a[],int n,int op){
if(n==1){
return;
}
int a0[n>>1|1],a1[n>>1|1];
for(int i=0;i<(n>>1);i++){
a0[i]=a[i<<1];
a1[i]=a[i<<1|1];
}
NTT(a0,n>>1,op);NTT(a1,n>>1,op);
int w=1,omega=qpow(G,op*(P-1)/n);
for(int i=0;i<(n>>1);i++,w=1ll*w*omega%P){
a[i]=(a0[i]+1ll*w*a1[i]%P)%P;
a[i+(n>>1)]=(a0[i]-1ll*w*a1[i]%P)%P;
}
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,m;
cin>>n>>m;
for(int i=0;i<=n;i++){
cin>>f[i];
}
for(int i=0;i<=m;i++){
cin>>g[i];
}
int limit=1;
while(limit<=n+m){
limit<<=1;
}
NTT(f,limit,1);
NTT(g,limit,1);
for(int i=0;i<limit;i++){
f[i]=1ll*f[i]*g[i]%P;
}
NTT(f,limit,-1);
int inv=qpow(limit,P-2);
for(int i=0;i<=n+m;i++){
int pl=1ll*f[i]*inv%P;
if(pl<0){
pl+=P;
}
cout<<pl<<' ';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
FFT/NTT 优化技巧
尽管 FFT/NTT 的时间复杂度为 \(\mathcal O\left(n\log n\right)\)，但是因为涉及递归和较为庞大的内存操作（数组），因此常数较大。可以进行一定程度上的优化。
FFT 三次变两次
对于基于复数的 FFT，可以将三次 FFT 操作转化为两次，从而优化常数。
求 \(f(x),g(x)\) 的乘积，一般是分别进行 FFT，求得点值表示法后直接相乘，再转化回系数表示法。
然而，可以将 \(g(x)\) 放到 \(f(x)\) 的虚部上得到 \(h(x)\)，求出 \(\dfrac{\operatorname{Im}(h^2(x))}{2n}\) 即答案。
因为：
虚部即 \(2ab\mathrm i\)，除以 \(2\)，得到 \(ab\)。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
#include<complex>
using namespace std;
constexpr const int N=1e6,M=N,Limit=1<<(int)(ceil(log2(N+M)))|1;
constexpr const double Pi=acos(-1);
#define complex complex<double>//±ãÓÚÊéд
complex f[Limit+1];
void FFT(complex a[],int n,int op){
if(n==1){
return;
}
complex a0[n>>1|1],a1[n>>1|1];
for(int i=0;i<(n>>1);i++){
a0[i]=a[i<<1];
a1[i]=a[i<<1|1];
}
FFT(a0,n>>1,op);FFT(a1,n>>1,op);
complex w={1,0},omega={cos(2*Pi/n),sin(2*Pi/n)*op};
for(int i=0;i<(n>>1);i++,w*=omega){
a[i]=a0[i]+w*a1[i];
a[i+(n>>1)]=a0[i]-w*a1[i];
}
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,m;
cin>>n>>m;
for(int i=0;i<=n;i++){
double pl;
cin>>pl;
f[i].real(pl);
}
for(int i=0;i<=m;i++){
double pl;
cin>>pl;
f[i].imag(pl);
}
int limit=1;
while(limit<=n+m){
limit<<=1;
}
FFT(f,limit,1);
for(int i=0;i<limit;i++){
f[i]*=f[i];
}
FFT(f,limit,-1);
for(int i=0;i<=n+m;i++){
int pl=round(f[i].imag()/limit/2);
cout<<pl<<' ';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
迭代版 FFT
参考代码
void FFT(vector<complex<double>> &a,int n,int op){
for(int i=1,j=0;i<n;i++){
int bit=n>>1;
for(;j&bit;bit>>=1){
j^=bit;
}
j^=bit;
if(i<j){
swap(a[i],a[j]);
}
}
for(int len=2;len<=n;len<<=1){
complex<double> wlen(cos(2*Pi/len),sin(2*Pi/len)*op);
for(int i=0;i<n;i+=len){
complex<double> w(1,0);
for(int j=0;j<len/2;j++){
complex<double> u=a[i+j];
complex<double> v=a[i+j+len/2]*w;
a[i+j]=u+v;
a[i+j+len/2]=u-v;
w*=wlen;
}
}
}
}
浙公网安备 33010602011771号