题解 数学 FFT/NTT 省选/NOI− 分治 多项式 未完

多项式学习笔记

FFT/NTT

Posted by TH911 on July 17, 2025

多项式乘法

本部分主要介绍 傅里叶变换 及 快速傅里叶变换（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：

    1

    complex<T>x;

  • 引用 x 的实部：

    1

    x.real();

    将 x 的实部赋值为 y：

    1

    x.real(y);

  • 引用 x 的虚部：

    1

    x.imag();

    将 x 的虚部赋值为 y：

    1

    x.imag(y);

  • 运算复数 x、y：

    1
2
3

    x+y;x+=y;
x-y;x-=y;
x*y;x*=y;

实际上，STL 的 complex 类或许比你手写的 complex 类更快。这里提供一份手写 complex 类：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

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）

例题链接

模板题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

//#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$ 的情况，需要处理一遍进位。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

//#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 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

//#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$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

//#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;
}
CATALOG
    FEATURED TAGS
    树状数组 题解 提高+/省选− 省选/NOI− DP 普及+/提高 基础算法 数学 组合数学 普及/提高− 贪心 字符串 线段树 未完 图论 平衡树 数论 树形 DP 搜索 Tarjan Trie 二分 位运算 二分答案 状态压缩 FHQ Treap 前端 无聊时候的一点思考 最短路 树链剖分 NOI/NOI+/CTSC
    FRIENDS