多项式学习笔记

多项式乘法

本部分主要介绍 傅里叶变换 及 快速傅里叶变换（FFT）等在 OI范围内 的应用。

参考 Blog 参考 Blog 参考 Blog 参考文章

本文会简要介绍 FFT 的原理、思想，涉及诸如 复平面、平面向量 等数学内容，文中会简要介绍。

FFT 快速傅里叶变换

傅里叶变换的本质就是将一个多项式在 系数表示法 和 点值表示法 之间转换。

FFT 能优化乘法的本质是点值表示法下可以直接乘。那么将多项式 \(f(x),g(x)\) 转换为点值表示法后直接相乘，再转换回系数表示法即可求出答案。

FFT 可以在 \(\mathcal O\left(n\log n\right)\) 的时间复杂度内计算出两个 \(n\) 次多项式的乘积。同时，考虑大整数可视为关于 \(10\) 的多项式，因此 FFT 也可以用于优化高精度乘法。

前置知识

复数与复平面

我们对实数域进行扩域可以得到复数域，这是最大的已知域。

一个复数 \(z\) 的代数形式形如：

\[z=a+b\mathrm i \]

其中，\(\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\) 对应向量的夹角），满足：

\[\tan \theta=\frac{b}{a} \]

称 \(\theta\) 为 \(z\) 的幅角，记为：

\[\theta=\operatorname{Arg}z \]

此时，复数的乘除法变得简单：

乘法：模相乘，幅角相加。
除法：模相除，幅角相减。

欧拉公式

\[e^{\mathrm ix}=\cos x+\mathrm i\sin x \]

单位根

称单位圆为复平面内以 \((0,0)\) 为原点，\(1\) 为半径的圆。

有 \(n\) 次方程：

\[x^n=1 \]

此方程的根有 \(n\) 个，称这 \(n\) 个根均为 单位根，则 \(n\) 次单位根等分了单位圆（因为模均为 \(1\)，幅角相加）。

设幅角为 \(\dfrac{2\pi}n\) 的单位复数：

\[\omega_n=e^{\frac{2\pi\mathrm i}{n}}=\cos\dfrac{2\pi}n+\mathrm i\sin\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\)）：

\[\omega_n^n=1\\ \omega_n^k=\omega_{2n}^{2k}\\ \omega_{2n}^{k+n}=-\omega_{2n}^k \]

代码实现

在 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)\)：

\[f(x)=a_0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1} \]

然而，多项式还可以通过点值表示。即确定 \(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)\)：

\[A(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1} \]

不妨令 \(n\) 为 \(2\) 的正整数次幂。若 \(n\) 不是，则可以将高位补全至 \(2\) 的正整数次幂，系数均为 \(0\) 即可。

考虑将 \(A(x)\) 按照奇偶性分类：

\[\begin{aligned} A_0(x)&=a_0x^0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2}\\ A_1(x)&=a_1x^0+a_3x^2+a_5x^4+\cdots+a_{n-1}x^{n-1}\\ \end{aligned} \]

则有：

\[A(x)=A_0(x)+xA_1(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)\) 可得：

\[A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k}) \]

将 \(\omega_n^{\large\frac n2},\omega_n^{\large\frac n2+1},\cdots,\omega_n^{n-1}\) 代入 \(A(x)\) 可得：

\[\begin{aligned} A\left(\omega_n^{k+\large\frac n2}\right)&=A_0\left(\omega_n^{2k+n}\right)+\omega_n^{k+\large\frac n2}A_1\left(\omega_n^{2k+n}\right)\\ &=A_0\left(\omega_n^{2k}\right)-\omega_n^kA_1\left(\omega_n^{2k}\right) \end{aligned} \]

可以发现，差别仅仅是一个正负号。

因此可以通过将 \(\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 中的单位根进行运算。

单位根满足：

\[\omega_n^n=1\\ \omega_n^k=\omega_{2n}^{2k}\\ \omega_{2n}^{k+n}=-\omega_{2n}^k \]

而原根在模 \(p\) 意义下同样满足这些性质。

最常见的模数与其原根：

\[p=998244353=7\times17\times2^{23}+1,g=3 \]

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){
	if(n<0){
		return qpow(base,-1ll*n*(P-2)%(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}\) 即答案。

因为：

\[(a+b\mathrm i)^2=a^2-b^2+2ab\mathrm i \]

虚部即 \(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].real()/limit);
		cout<<pl<<' ';
	}
	
	cout.flush(); 
	
	/*fclose(stdin);
	fclose(stdout);*/
	return 0;
}

posted @ 2025-07-22 19:47 TH911 阅读(18) 评论(0) 收藏举报

刷新页面返回顶部

TH911 Blog

光阴不可测，岁月亦无歌。

多项式学习笔记

多项式乘法

FFT 快速傅里叶变换

前置知识

复数与复平面

几何意义

欧拉公式

单位根

代码实现

DFT/IDFT 离散傅里叶（逆）变换

多项式的点值表示法

IDFT 离散傅里叶逆变换

FFT 快速傅里叶变换

实现细节

FFT 例题

P3803 【模板】多项式乘法（FFT）

P1919 【模板】高精度乘法 | A*B Problem 升级版

NTT 快速数论变换

前置知识

原根与阶

NTT 快速数论变换

P3803 【模板】多项式乘法（FFT） AC 代码

FFT/NTT 优化技巧

FFT 三次变两次

公告