多项式乘法 FFT&NTT

hhh，最开心难绷的一集，多项式算法，博大精深啊（话说它怎么居然现在被加到NOI大纲了啊。。。

多项式乘法

对于多项式 \(f(x)=\sum_{i=0}^n a_ix^i\)

如果我们对两个多项式做乘法，会得到：\(f(x)g(x)=\sum_{i=0}^{2n}(\sum_{j=0}^i a_jb_{i-j})x^i\)

其中新多项式的系数可以被写成 \(c_n=\sum_{i=0}^na_ib_{n-i}\)

这样的形式被称为卷积，一般来说就是相乘的两个数的下标为定值。（注意到，这只是卷积的一种，卷积的严格定义是两个多项式经过一个特定的“卷积算子”得到的结果，而并不仅是乘法

快速傅立叶变换 FFT

(Fast Fourier Transform）

卷积的暴力计算，复杂度是 \(O(n^2)\) 的。而FFT，就是来快速求解卷积的。~~(它的得名就有一个是因为它非常fast~~

笔者在看了一堆玄妙的转来转去的东西之后深感无力，到头来竟发现，我们可以只本着信息奥赛的特点，从更为简单的角度理解FFT。

前置芝士

系数表示与点值表示

一个多项式形如 \(f(x)=\sum_{i=0}^n a_i x^i\) ，这就是它的系数表示

重要的是点值表示，我们把它看成一个函数 \(y=f(x)\) ，则我们可以用它经过的n+1个点表示这个函数，这个函数就是确定的了。\((x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\) 。

考虑两个多项式 \(f(x)\) 与 \(g(x)\) 相乘，我们惊喜的发现，点值表示上，他们的卷积结果是 \((x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),...,(x_n,f(x_n)g(x_n))\)。此时我们只需要 \(O(n)\) 的复杂度就可以直接求出两个多项式卷积的结果了。考虑如何在多项式的系数表示和点值表示间进行快速的转换。

复数

讲个笑话，你知道复数为什么叫“复数”吗？它的英文是 complex number，就是复杂的数hhh

对于 \(x^2=-1\) 的东西，显然它在实数范围内无解，但我们有时要解决这类问题。

因此定义 \(i^2=-1\) ，其中 \(i\) 被称为虚数单位。

复数=实数部分+虚数部分

即，任何一个复数都可以被表示为 \(a+bi\) 的形式（\(a,b\in R\)），全体复数集为 \(C\) .

其中 \(a\) 为实部，\(b\) 为虚部。

n次单位根

数1在复数范围内的n次方根，称为n次单位根。也就是\(x^n=1\)在复数域内的根。我们将这个方程的所有根对应到复平面上。以原点为圆心，做一个半径为一的圆，称为单位圆。则这些根会n等分单位圆，且始终有一个根为1。

其中这些根的表达式就是 \(\omega_n=cos(\frac{2k\pi}{n})+sin(\frac{2k\pi}{n})i\)

我们称这n个根为n次单位根，同时将从1开始逆时针找到的第一个根记为 \(\omega_n\)。它的性质就是，剩下的所有单位根都能用 \(\omega_n\) 的幂次表示，即n次单位根分别为 \(\omega_n^0,\omega_n^1,...,\omega_n^{n-1}\) 。
（\(\omega_n\) 就是表达式中 \(k=1\) 的情况）

其他n次单位根的重要性质：(可以简单画图理解)

\(\omega_n^k=\omega_n^{k+n}\)
\(\omega_n^k=\omega_{2n}^{2k}\)
\(\omega_{2n}^{n+k}=-\omega_{2n}^k\)

离散傅立叶变换DFT

实际上描述的就是将多项式从系数表示转化成点值表示的一种方法。

如对多项式 \(f(x)=\sum_{i=0}^{n-1} a_ix^i\) ，我们将 \(n\) 个不同的 \(n\) 次单位根带入 \(f(x)\) 中，就能得到它对应的点值表达，但这样暴力做的复杂度依然是 \(O(n^2)\) 的。FFT就提供了一种快速进行DFT的算法。

考虑一个\(2n-1\)次多项式\(f(x)\)，如果我们分别设：
\(g(x)=a_0+a_2x+a_4x^2+...\)
\(h(x)=a_1+a_3x+a_5x^2+...\)

那么我们有 \(f(x)=g(x^2)+xh(x^2)\)

即 \(f(\omega_{2n}^k)=g(\omega_{2n}^{2k})+\omega_{2n}^kh(\omega_{2n}^{2k})=g(\omega_n^k)+\omega_{2n}^kh(\omega_n^k)\)
\(f(\omega_{2n}^{n+k})=g(\omega_{2n}^{2n+2k})+\omega_{2n}^{n+k}h(\omega_{2n}^{2n+2k})=g(\omega_{2n}^{2k})-\omega_{2n}^kh(\omega_{2n}^{2k})=g(\omega_n^k)-\omega_{2n}^kh(\omega_n^k)\)

我们发现了一个分治的思路，这样递归下去复杂度是 \(O(nlogn)\) 的。

但分治常数过大，我们一般把多项式的项数补到一个2的整数次幂，然后使用倍增法实现。倍增实现有两个要点：~~（叫的名字还都怪高级的嘞~~

位逆序置换 对一个8项多项式，考虑分治过程，
初始下标为 \(\set{0,1,2,3,4,5,6,7}\)
第一次分治 \(\set{0,2,4,6},\set{1,3,5,7}\)
第二次分治 \(\set{0,4},\set{2,6},\set{1,5},\set{3,7}\)
第三次分治 \(\set{0},\set{4},\set{2},\set{6},\set{1},\set{5},\set{3},\set{7}\)

初始下标：\(000,001,010,011,100,101,110,111\)
结束下标：\(000,100,010,110,001,101,011,111\)

我们发现了一件神奇的事情，就是初始下标和结束下标都是互相对称的，这就是位逆序置换。我们可以提前 \(O(n)\) 预处理出每一个数字位逆序置换后的结果，然后交换每一对需要交换的数。

蝶形运算优化 我们发现在拆分系数之后，\(g(\omega_n^k)\) 存储在第 \(k\) 项，而 \(h(\omega_n^k)\) 存储在第 \(k+n\) 项。而合并之后 \(f(\omega_{2n}^k)\) 存储在第 \(k\) 项，\(f(\omega_{2n}^{n+k})\) 存储在第 \(k+n\) 项。我们完全可以对第 \(k\) 和 \(k+n\) 项进行覆盖，且把这两个放在一起算还能少一半常数（因为复数运算很慢）。（不知道为什么的可以翻着看上面推导的式子）

逆离散傅立叶变换IDFT

顾名思义，就是跟DFT相反的操作，把点值表示转化为系数表示。事实上我们可以把IDFT的过程转化为DFT的过程。

贴了 UKE_Automation 大佬写的博客的截图，~~因为不想打矩阵了~~

小结：FFT流程

多项式系数表示->点值表示->系数表示

注意到对于点值表示，只要它选择的横坐标相同（这里我们选择了n次方根的若干次方，是因为它的优秀性质），多项式的加减乘除就可以直接在纵坐标上计算。最后因为我们得到的是一个n+m次多项式，我们需要得到n+m+1个点的坐标，因此每个多项式选点都要选至少n+m+1个点才行。

我们要对两个多项式分别做一次DFT，然后再做一次IDFT，总共三次FFT操作。听说有神秘常数优化，可以只做两次FFT操作。

FFT模板

#include<bits/stdc++.h>
using namespace std;
const int maxn=4e6+5;//得开四倍，补全的2倍，两个多项式算下来是四倍 
const double pi=acos(-1);//反余弦函数，返回余弦值对应的弧度 
int n,m,r[maxn],mx=1,lg;
struct cplx{
	double x,y;
	cplx operator +(cplx b){ return {x+b.x,y+b.y};}
	cplx operator -(cplx b){ return {x-b.x,y-b.y};} 
	cplx operator *(cplx b){ return {x*b.x-y*b.y,x*b.y+y*b.x};}
	cplx operator /(int b){ return {x/b,y/b};}
};
struct polyn{
	int n;
	cplx a[maxn];
	void fft(int typ){
		for(int i=0;i<mx;i++){
			if(i<r[i]) swap(a[i],a[r[i]]);
		}
		for(int len=1;len<mx;len<<=1){
			cplx now={cos(pi/len),typ*sin(pi/len)};//2pi/2len 约了 
			for(int i=0;i<mx;i+=(len<<1)){
				cplx w={1,0};
				for(int j=0;j<len;j++,w=w*now){
					cplx x=a[i+j],y=w*a[i+j+len];
					a[i+j]=x+y,a[i+j+len]=x-y;
				}
			}
		}
		if(typ==-1){
			for(int i=0;i<mx;i++) a[i]=a[i]/mx;
		}
	}
}F,G;
int main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&F.a[i].x);
	for(int i=0;i<=m;i++) scanf("%lf",&G.a[i].x);
	while(mx<=n+m) mx<<=1,lg++;
	for(int i=0;i<mx;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(lg-1));//直接处理反转所对应的下标
	F.fft(1),G.fft(1);
	for(int i=0;i<=mx;i++) F.a[i]=F.a[i]*G.a[i];
	F.fft(-1);
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(F.a[i].x+0.5));//注意这里，直接输出会向下取整，为了避免精度问题，我们要四舍五入 
	return 0;
}

例题

「BZOJ2194」快速傅里叶之二

题目要求 \(c_k=\sum_{i=k}^{n-1} a_ib_{i-k}\) ，我们一般求的卷积都是两项和为定值，那么不妨考虑一个转化，把 \(a\) 倒序过来，然后 \(a_i\) 变成 \(a_{n-i}\) ，我们就能得到和为定值的形式

「ZJOI2014」力

发现那个 \(q_j\) 乘上了又除掉了，我们直接不考虑。考虑前面一坨减后面一坨可以直接拆成两部分算，还是通过上一题的套路把它转化为和一定的卷积形式。

注：推这种式子的时候注意一下，是和为多少加第几项，跟从哪里开始枚举关系不大。

「P10632」Normal

考虑 \(e=P\times w\) ，我们要求的相当于 \(\sum_{i=1}^n\sum_{j=1}^n P(\text{j在i为分治中心的子树内})\)，进一步转化为 \(\sum_{i=1}^n\sum_{j=1}^n \frac{1}{dis(i,j)}\)，也就是 \(i\) 到 \(j\) 路径上的点，要满足要求，\(i\) 必须是第一个被删掉的。

然后考虑用点分治求出每种长度的点有多少个。经典手法是子树容斥。先求出当前子树内所有点到分治中心的距离，然后考虑组合。设距离分治中心为i的点有 \(p_i\) 个，那么合并后长度为 \(n\) 的路径数有 \(cnt_n=\sum_{i=1}^{n-1}p_i\times p_{n-i}\) 。然后把在同一子树内的方案减掉，对每个子树再做一遍这东西。复杂度\(O(nlog^2n)\)。

这题会炸int。。。

「BZOJ2194」Triple

因为只选三个，所以可以分类讨论。这题类似生成函数，构造一个 \(a_ix^i\) 表示占用了 \(i\) 的空间的选择有 \(a_i\) 种方案数。

令选一个的多项式为 \(F\)。

选两个，\(F\times F- G\)，其中 \(G\) 表示选相同的两个的多项式。且最后要除以二。

选三个，\(F\times F\times F -3\times(G\times F-\times H)- H\)。其中 \(H\) 表示选相同的三个的多项式。最后除以六。

「BZOJ3160」万径人踪灭

答案等于位置对称的回文子序列数-回文子串数（这可以用manacher计算）

我们需要求位置对称的回文子序列数，就相当于寻找有多少个字符，在对称中心两侧对应位置且字符相同，然后根据对称中心在夹缝里或字符上分类讨论贡献方案数就行。

这题给我们的启示在于：FFT可以求解字符串问题。

同样是类似于生成函数（没见过的想不到多项式啊。。。），我们考虑构造两个多项式，其中 \(A\) 为所有原串中为 \(a\) 的位置系数为1，\(b\) 系数为0；\(B\) 相反。

然后 \(A\),\(B\) 自乘，就能求出答案了。

快速数论变换 NTT

有这样一类问题，要求对多项式卷积的结果取模，在数据范围很大的情况下，用FFT可能会不够，或者炸精度。这就需要NTT来解决。

我们考虑需要改动FFT——实际上只有单位根不能用了。考虑单位根的性质如下：

\(\omega_n^k=\omega_n^{k+n}\)
\(\omega_n^k=\omega_{2n}^{2k}\)
\(\omega_{2n}^{n+k}=-\omega_{2n}^k\)
\(\omega_n^k\) 对于 \(k\in [0,n-1]\) 互不相等，且 \(\omega_n^n=1\)

如果你学过原根，你就会难以发现，原根可以替代这个单位根。

最后一个性质就是原根的性质，一个质数 \(p\) 的原根 \(g\) 满足 \(g^k,k\in[1,p)\) 在模 \(p\) 意义下互不相等。为满足 \(\omega_n^n=1\) 的条件，我们要找出一个 \((g^r)^n\equiv 1 \pmod{p}\) ，根据欧拉定理，我们只需要让 \(rn=\varphi(p)\) 就行，\(r=\frac{\varphi(p)}{n}=\frac{p-1}{n}\)，因此 \(g^{\frac{p-1}{n}}\) 就是新的单位根。然后这几条性质都可以证明。

注意到我们新单位根的形式要求 \(n|p-1\) ，又 \(n\) 就是 \(2^k\) ，所以我们的 \(p\) 应该是 \(a2^k+1\) 的形式。常见 NTT 模数：

\(998244353=119\times 2^{23}+1,g=3\)
\(469762049=7\times 2^{26}+1,g=3\)
\(1004535809=479\times 2^{21}+1,g=3\)

换掉单位根之后我们就可以进行整数运算而不用到复数了，对应的，我们还需要把 IDFT 中 \(\omega_n^{-k}\) 的部分换成求逆元。

NTT 主要部分

const int mod=1004535809,g=3,invg=334845270;
void ntt(int typ){
	for(int i=0;i<mx;i++){
		if(i<r[i]) swap(a[r[i]],a[i]);
	}
	for(int len=1;len<mx;len<<=1){
		int now=qpow((typ==1)?g:invg,(mod-1)/(len<<1),mod);
		for(int i=0;i<mx;i+=(len<<1)){
			int w=1;
			for(int j=0;j<len;j++,w=w*now%mod){
				int x=a[i+j],y=w*a[i+j+len]%mod;
				a[i+j]=(x+y)%mod,a[i+j+len]=(x-y+mod)%mod;
			}
		}
	}
	if(typ==-1){
		int inv=qpow(mx,mod-2,mod);
		for(int i=0;i<mx;i++) a[i]=a[i]*inv%mod;
	}
}

ps.更有一类变态出题人，一定要让你用一个神秘大模数做ntt。那么就需要用任意模数ntt，有很多种写法，其中较为简单且常数小的一种是拆系数fft。

我们令每个数 \(x=a\times M+b\) ，这样两个数相乘就可以变为 \(x_1x_2=(a_1\times M+b_1)(a_2\times M+b_2)\) 了。也就是我们只需要分别维护 \(a_1,a_2,b_1,b_2\) 的数组，做四次 DFT ；然后分别对 \(a_1a_2,a_1b_2+b_1a_2,b_1b_2\) 做 IDFT 即可。一共是七次fft。我们的 \(M\) 一般取 \(2^15\) ，这样对于一个长为 \(1e5\) 左右的多项式，最大的过程值在 \(10^{14}\) 级别，是long double所能承受的。

还有一种优化是，我们发现“实数部 \(\times\) 实数部 \(\eq\) 实数部”“实数部 \(\times\) 虚数部 \(\eq\) 虚数部”，那就不放把 \(a_2,b_2\) 分别放在同一个函数的实数部和虚数部上，然后做三次 DFT 两次 IDFT 就行了。

例题

「SDOI2015」序列统计

权值不是我们擅长用生成函数处理的加法，而变成了乘法。不妨考虑一种经典转化，取对数。在模 \(m\) 意义下，就是指标了。

我们对于 \(s\) 中所有数取指标，那么原题就可以转化为在模 \(\varphi(m)\) 意义下n个数的指标之和等于 \(x\) 的指标的选择方案数。相当于原多项式的 \(n\) 次方中 \(x\) 次项对应的系数，对多项式做快速幂。然后对 \(m\) 取模就是直接把大于等于 \(\varphi(m)\) 的项放到前面并清空，这样我们对多项式做快速幂时就不会炸项数了。

「HAOI2018」染色

~~一看这种我没有一点头绪的题，一般就是数学了~~

满足要求的最多颜色数是 \(min(m,\frac{n}{s})\) ，我们想求恰好有 \(k\) 种颜色出现次数为 \(s\) 的方案数，是不好求的。~~（不要和我一样唐觉得直接拿生成函数能求）~~

考虑经典的二项式反演，我们先钦定 \(k\) 种颜色出现 \(s\) 次，记为 \(g(k)\) ；恰好 \(k\) 种颜色出现 \(s\) 次我们记为 \(f(k)\)。得到：

\[g(k)=\binom{m}{k}\times\frac{n!}{(s!)^k(n-s\times k)!}\times(m-k)^{n-s\times k} \]

就是选一些颜色 \(\times\) 可重排列数 \(\times\) 剩下任意。

考虑我们算重了多少次，枚举一下可以发现，\(g(k)=\sum_{i=k}^m\binom{i}{k}\times f(i)\)。也就是每次都计算了（钦定选择的东西）在（满足要求的东西）中的方案数次。

我们根据二项式反演可以得到：\(f(k)=\sum_{i=k}^m(-1)^{i-k}\binom{i}{k}g(i)\)

发现这个式子还是 \(O(n^2)\) 的，我们考虑能不能拆一下。

\(f(k)=\sum_{i=k}^m(-1)^{i-k}\times\frac{i!}{k!(i-k)!}\times g(i)=\frac{1}{k!}\sum_{i=k}^m\frac{(-1)^{i-k}}{(i-k)!}\times i! \times g(i)\)

我们发现后面的和式居然满足卷积的形式！！这给我们的启示是，组合数等东西也可能能用卷积来优化！！！

\(f(k)\) 求出来，答案也就出来了。

还有各种高级卷积，比如位运算卷积 FWT,FMT 等等，回头补坑？放大佬博客

分治 FFT/NTT

如果题目要求，两部分之间是多项式运算，且第一部分先得到，然后对第二部分做贡献，得到第二部分。这样我们就可以通过分治 FFT/NTT 把算法优化到 \(O(nlog^2n)\) .（暴力操作是 \(O(n^2)\) 的）

分治 FFT/NTT 的过程整体类似 cdq分治，只是注意边界条件和中间做卷积运算就行了。

这里给出一道例题和代码：CF1553I Stairs

先简化题意，每个标着a[i]的点都至少要有连着的a[i]个，把这样的一组划分成一个“段”。只需要排每个段的相对顺序并考虑升降序。

显然是要容斥的，因为如果升降序相同且前后给出的自然数相接，会连成一个大段。需要一个dp，因为长度为1的段升降序是一样的。我们只需要求出来G[i]表示至多i个不同连续段，就可以得到F[i]表示恰好i个不同连续段。

dp[i][j]表示到第i位，有至多j个连续段。这个dp经过前缀和优化仍是\(O(n^2)\)的，且j这维省不了。发现转移都是从 i-1,i-2,j-1 这样的位置转移过来，我们可以把它写成矩阵和生成函数的形式！！然后用分治NTT 就能到 \(O(nlog^2n)\) 。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=998244353,g=3,invg=332748118;
const int maxn=1e5+5;
int n,a[maxn],m,b[maxn],rv[maxn<<2],mx,lg;
ll dp[2005][2005],sum[2005][2005],fac[maxn],ans;
struct Polyn{
	vector<ll>a; int n;
	void reset(int len){ a.resize(len); n=len-1; }
	Polyn operator * (Polyn x){
		Polyn y; y.reset(max(n,x.n)+1);
		for(int i=0;i<=n;i++) y.a[i]=a[i]*x.a[i]%mod;
		return y;
	} 
	Polyn operator + (Polyn x){
		Polyn y; y.reset(max(n,x.n)+1);
		for(int i=0;i<=n;i++) y.a[i]=(a[i]+x.a[i])%mod;
		return y;
	}
};
ll qpow(ll x,ll y){
	ll res=1;
	while(y){
		if(y&1) res=res*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return res;
}
void ntt(Polyn &A,int typ){
	A.reset(mx);
	for(int i=0;i<mx;i++){
		if(i<rv[i]) swap(A.a[i],A.a[rv[i]]);
	}
	for(int len=1;len<mx;len<<=1){
		ll now=qpow((typ==1)?g:invg,(mod-1)/(len<<1));
		for(int i=0;i<mx;i+=(len<<1)){
			ll w=1;
			for(int j=0;j<len;j++,w=w*now%mod){
				ll x=A.a[i+j],y=w*A.a[i+j+len]%mod;
				A.a[i+j]=(x+y)%mod,A.a[i+j+len]=(x-y+mod)%mod;
			}
		}
	} 
	if(typ==-1){
		int iv=qpow(mx,mod-2);
		for(int i=0;i<mx;i++) A.a[i]=A.a[i]*iv%mod;
	}
}
struct mat{
	Polyn p[2][2];
	mat operator * (mat x){
		mat y;
		y.p[0][0]=p[0][0]*x.p[0][0]+p[0][1]*x.p[1][0];
		y.p[0][1]=p[0][0]*x.p[0][1]+p[0][1]*x.p[1][1];
		y.p[1][0]=p[1][0]*x.p[0][0]+p[1][1]*x.p[1][0];
		y.p[1][1]=p[1][0]*x.p[0][1]+p[1][1]*x.p[1][1];
		return y;
	}
};
mat cdq(int l,int r){
	if(l==r){
		mat ret;
		ret.p[0][0].reset(2),ret.p[0][1].reset(2),ret.p[1][0].reset(2),ret.p[1][1].reset(2);
		ret.p[0][1].a[0]=ret.p[1][1].a[0]=1,ret.p[1][0].a[1]=2,ret.p[0][0].a[1]=(b[l]==1)?1:2;
		return ret;
	}
	int mid=(l+r)>>1;
	mat pl=cdq(l,mid),pr=cdq(mid+1,r),ret;
	mx=1,lg=0;
	while(mx<=r-l+1) mx<<=1,lg++;
	for(int i=0;i<mx;i++) rv[i]=(rv[i>>1]>>1)|((i&1)<<(lg-1));
	for(int i=0;i<2;i++) for(int j=0;j<2;j++) ntt(pl.p[i][j],1),ntt(pr.p[i][j],1);
	ret=pl*pr;
	for(int i=0;i<2;i++) for(int j=0;j<2;j++) ntt(ret.p[i][j],-1);
	return ret;
}
int main(){
	scanf("%d",&n),fac[0]=1;
	for(int i=1;i<=n;i++){
		scanf("%d",&a[i]);
		fac[i]=fac[i-1]*i%mod;
	}
	for(int i=1;i<=n;i++){
		int j=i;
		while(j<n&&a[j+1]==a[i]) j++;
		if((j-i+1)%a[i]!=0){
			printf("0");
			return 0;
		}
		int k=(j-i+1)/a[i]; i=j;
		while(k--) b[++m]=a[i];
	}
	Polyn res=cdq(1,m).p[0][0]; ll fl=1;
	for(int i=m;i>0;i--){
		ans=(ans+fl*fac[i]*res.a[i])%mod,fl=-fl;
	}
	printf("%lld",(ans+mod)%mod);
	return 0;
}

倍增 FFT/NTT

实际上跟普通的倍增法差不多。就是考虑一个转移能写成从i到2i的形式。重点是关注它和分治的区别：倍增是要求dp转移有明显的顺序性或初始化条件，而分治则是关注每一个位置本身的限制条件。

例题：CF773F Test Data Generation

posted @ 2025-03-05 17:43 theWeeds'Defense 阅读(82) 评论(0) 收藏举报

刷新页面返回顶部

theWeeds'Defense

——We rise where tides forget to kiss the shore; between the cracks of sidewalks,cold and sore.