多项式全家桶

不知道什么时候写封装，有空就写。这里是我使用的封装好的 modint。

FWT：link。

FWT 同时带给我们一个深切的见解：一切多项式算法都是把多项式卷积转化成可以点积，再逆回去。

FFT

处理实数意义下的多项式 \((+,\times)\) 卷积。也就是多项式乘法。

点值表示法和系数表示法不必赘述了。

对于两个以系数表示法表示的多项式，显然我们可以 \(\mathcal O(n)\) 点起来。如果我们存在一个较快的办法把系数表示法的多项式转化为点值表示法，并且支持转回去就好了。

使用朴素的 DFT 和 IDFT（代入和拉插）可以做到 \(\mathcal O(n^2)-\mathcal O(n)\)，和直接暴力卷没什么区别。

下面介绍 FFT。

考虑先优化代入的操作。考虑比较优秀的点值可以带来一些性质，帮助我们快速算出多项式的值。

根据代数基本定理，任何 \(x^n=1\) 的方程都有 \(n\) 个根。我们不妨用这 \(n\) 个根作为点值表示一个 \(n\) 次多项式。根据傅里叶对代数基本定理的证明，这 \(n\) 个根在复平面上平分了单位圆。

类似平面直角坐标系上的单位圆，复平面上的单位圆上与原点连线相对于 \(x\) 轴正半轴旋转了 \(\theta\) 的点表示为 \(\cos\theta+i\sin\theta\)。

那么显而易见的，我们从 \((1,0)\) 开始沿着单位圆逆时针旋转，找到的第一个根应该是 \(\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\)。我们把这个根称为 \(n\) 处的单位根，记为 \(\omega_n\)。

根据欧拉公式 \(\cos\theta+i\sin\theta=e^{i\theta}\)，有 \(\omega_n=e^{i\frac{2\pi}{n}}\)。

从而所有根都是 \(\omega_n\) 的幂次。也就是说，取 \(p\in[0,n)\)，\(\omega_n^p\) 可以表示 \(x^n=1\) 的所有根，并且显然这个东西是循环的，即 \(\forall i\ge n,\omega_n^i=\omega_n^{i\bmod n}\)。

并且利用指数的性质和 \(e^{i\pi}=-1\) 我们可以轻易证明：

\[\begin{cases} \omega_{dn}^{dk}=w_{n}^k & ,\forall n,k\in\N,d\in\N^+\\ (\omega_{n}^{k+\frac{n}{2}})^2=\omega_n^{2k+n}=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^k & ,\forall n,k\in\N,n\text{ is even} \end{cases} \]

利用这些性质进行 FFT。

假设多项式项数 \(n\) 是 \(2\) 的次幂，令 \(n\gets 2^n\)，下文中 \(n\) 均表示指数。显然不够 \(2\) 的次幂在 \((+,\times)\) 卷积下用 \(0\) 补齐就好。

我们维护点值 \(\omega_n^0,\omega_n^1\cdots \omega_n^{n-1}\) 处多项式的值。

考虑多项式 \(\sum\limits_{i=0}^{n-1}a_ix^i\)，进行下面的一通操作：

pE9Xbh6.md.png

总而言之就是对于 \(k\in[0,\frac{n}{2})\)，有 \(A_1(\omega_{\frac{n}{2}}^k)\) 和 \(A_2(\omega_{\frac{n}{2}}^k)\) 我们就能求出所有 \(A(\omega_n^k),A(\omega_{n}^{k+\frac{n}{2}})\)。考虑这显然是一个子问题，分治下去做就可以了。

IFFT

考虑我们的 FFT 可以写成一个矩阵乘法状物：

\[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{1}_n & \omega^{2}_n & \cdots & \omega^{n - 1}_n \\ 1 & \omega^{2}_n & \omega^{4}_n & \cdots & \omega^{2n - 2}_n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{n - 1}_n & \omega^{2n - 2}_n & \cdots & \omega^{(n-1)^2}_n \end{bmatrix} \times \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ A(\omega_n^2) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \]

左侧每一行都表示一个取值，第 \(i\) 行表示取 \(\omega_n^{i-1}\) 的点值多项式的每一项，右侧就是系数矩阵。乘出来当然就是多项式的点值。

考虑怎么倒回去。注意到：

\[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{-1}_n & \omega^{-2}_n & \cdots & \omega^{-n + 1}_n \\ 1 & \omega^{-2}_n & \omega^{-4}_n & \cdots & \omega^{-2n + 2}_n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-n + 1}_n & \omega^{-2n + 2}_n & \cdots & \omega^{-(n-1)^2}_n \end{bmatrix} \times \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ A(\omega_n^2) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} = n \times \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{bmatrix} \]

所以把 \(\omega_{n}^k\) 转成 \(\omega_{n}^{-k}\)，做一遍 FFT 最后除以 \(n\) 就好了。

写法上需要维护复数类，namespace::std 中自带的 complex<> 类非常好使。我们声明一个 complex <double> 类变量，正常使用即可：

.real() 获取实部，.imag() 获取虚部。
可以赋值为整数类。
使用 ={x,y} 可以赋值为 \(x+iy\)。
其余四则运算均可照常。

于是我们可以写出美丽的分治版 FFT，实现上参考了 Hch 的：

const double pi=acos(-1.0);
#define complex complex<double>
void FFT(vector <complex> &A,double coef){
	int n=A.size();
	if(n==1) return;
	vector <complex> A1,A2;
	fr1(i,0,n-1) (i&1?A2:A1).pb(A[i]);
	FFT(A1,coef);
	FFT(A2,coef);
	double theta=coef*2.0*pi/(double)n;
	complex wn={cos(theta),sin(theta)};
	complex wp=1;
	fr1(i,0,n/2-1){
		A[i]=A1[i]+wp*A2[i];
		A[i+n/2]=A1[i]-wp*A2[i];
		wp*=wn;
	}
}
vector <int> multi(vector <int> a,vector <int> b){
	vector <complex> A(a.begin(),a.end());
	vector <complex> B(b.begin(),b.end());
	int t=A.size()+B.size()-1;
	int st=t;
	while(lowbit(t)!=t) t++;
	A.resize(t),B.resize(t);
	FFT(A,1);
	FFT(B,1);
	fr1(i,0,t-1) A[i]=A[i]*B[i];
	FFT(A,-1);
	vector <int> c(st);
	fr1(i,0,st-1) c[i]=(A[i].real()/(double)t+0.5);
	return c;
}

注意在多项式中所有系数表示多项式的 vector 中系数都是从低次到高次的。此外需要注意：\(a,b\) 项多项式有 \(a+b-1\) 项，这里 \(t\) 表示乘出来的多项式的项数，也就是 FFT 的项数，这个东西需要补全到 \(2^k\)。而 \(a+b-2\) 是多项式的次数。

然而这玩意常数太大了！作为多项式全家桶我们应该卡常，下面介绍非递归写法。

蝴蝶变换

考虑类似 FWT，我们如果可以把分治树从下层往上层处理就好了。我们先直接跑到最后一层，然后向上合并。

FWT 非常简单，就是直接两半边，但是这里需要按照奇偶分治，怎么办呢？

考虑这个例子：

      0 1 2 3 4 5 6 7
1st:  0 2 4 6|1 3 5 7
2nd:  0 4|2 6|1 5|3 7
3rd:  0|4|2|6|1|5|3|7

注意到任何一个数在第 \(i\) 次分治后所在的位置就是它二进制表示循环右移移位到达的位置（如果最高是 \(2^k-1\)，那么只考虑 \(k\) 位），特别地，不计算最后一次。比如 \(3=(011)_2\)，那么第一次分治到达了 \((101)_2=5\)，第二次分治就到达了 \((110)_2=6\)，第三次分治不管。

于是我们可以轻易先搞到最后一层。只需求出最终位置 \(pos\) 后，把 \(i<pos_i\) 的交换即可（显然要么是两个数交换，要么是不动）。然后我也不知道为什么……只需要在 FFT 的时候做一遍这个，再在 IFFT 的时候做一遍这个，就能把奇偶分治看成折半分治了。所以我们可以直接类似 FWT 的写法进行操作：

const double pi=acos(-1.0);
#define complex complex<double>
vector <int> pos;
void FFT(vector <complex> &A,double coef){
	int n=A.size();
	fr1(i,0,n-1) if(i<pos[i]) swap(A[i],A[pos[i]]);
	for(int p=2,q=1;p<=n;p<<=1,q<<=1){
		double theta=coef*2*pi/(double)p;
		complex wp={cos(theta),sin(theta)};
		for(int i=0;i<n;i+=p){
			complex wn=1;
			fr1(j,i,i+q-1){
				complex x=A[j],y=A[j+q]*wn;
				A[j]=x+y;
				A[j+q]=x-y;
				wn*=wp;
			}
		}
	}
}
vector <int> multi(vector <int> a,vector <int> b){
	vector <complex> A(a.begin(),a.end());
	vector <complex> B(b.begin(),b.end());
	int t=A.size()+B.size()-1;
	int st=t;
	while(lowbit(t)!=t) t++;
	A.resize(t),B.resize(t);
	pos.resize(t);
	int lg=__lg(t);
	fr1(i,1,t-1) pos[i]=(pos[i>>1]>>1)|((i&1)<<lg-1);
	FFT(A,1);
	FFT(B,1);
	fr1(i,0,t-1) A[i]=A[i]*B[i];
	FFT(A,-1);
	vector <int> c(st);
	fr1(i,0,st-1) c[i]=(A[i].real()/(double)t+0.5);
	return c;
}

可以比递归版快一倍。

NTT

在模意义下进行 FFT。

结论：考虑模数 \(p\) 的原根 \(g\)，那么 \(\omega_n\equiv g^{\frac{p-1}{n}}\pmod p\)。所以直接写一个 FFT，一切照常，所有幂次 \(-1\) 和除法处使用逆元（INTT 需要使用），乘法加法取模就完事了。

此处要求 \(p-1\) 必须是 \(n\) 的整数倍，这也意味着 \(p-1\) 中需要有大量的 \(2\) 因子。所以人们常说 \(998244353\) 是 NTT 模数，因为 \(998244352=7\times 17\times 2^{23}\)，这意味着它可以处理项数小于等于 \(2^{23}\) 的多项式，这已经是了大多数多项式题不敢达到的数据范围了。

寻找原根是简单的：

原根存在定理：\(2,4,p^\alpha,2p^\alpha\) 存在原根，其中 \(p\) 为奇素数，\(\alpha\in\N^+\)。
原根判定定理：对于 \(\varphi(p)\) 的所有因子 \(i\)，均有 \(g^{\frac{\varphi(p)}{i}}\not\equiv 1\pmod p\)，证明 \(g\) 是原根。这里用了原根满足 \(g^i,i\in[0,\varphi(p))\) 及任意一个等长的连续幂集等于整个模 \(p\) 剩余系的性质。

论文拼盘证明了原根的级别为 \(\mathcal O(p^{0.25})\)，所以暴力枚举原根就行了。

任意模数 NTT 暂时没有研究。一般来说可以被三模 NTT 平替掉。找三个乘积大于最大系数的 NTT 模数，做三遍 NTT 之后 CRT。~~现代算法竞赛谁卡三倍常数啊。~~

为了方便封装，我们直接全部一步到位（找原根也封装在一起了）。为了方便没有考虑模数非质数的情况，直接取 \(\varphi(p)=p-1\) 了。

多项式全家桶

多项式求逆

求解 \(f(x)g(x)\equiv 1\pmod {x^n}\)（\(n\) 项多项式），系数对 \(p=998244353\) 取模。

如果没有模数可以分治 FFT。具体来说考虑 \(f(x)\) 的系数 \(\{a_n\}\)，\(g(x)\) 的系数 \(\{b_n\}\)，那么 \(a_0b_0=1\)，并且 \(\sum\limits_{i=0}^{k} a_ib_{k-i}=0\)，类似约数容斥地把 \(b_k\) 项拆出来，那么有：

\[\sum\limits_{i=1}^{k}a_ib_{k-i}=-a_0b_k\iff b_k=\sum\limits_{i=1}^k(-\frac{a_i}{a_0})b_{k-i} \]

观察到右侧是一个 cdq 分治状物，我们只需知道 \(b_{0\sim k-1}\) 就能推出 \(b_k\)。所以先分治计算 \(b_{l\sim mid}\)，它对右侧的贡献是一个卷积，用 FFT 算一下就完了。

其实有点蠢。

常规地，我们使得 \(n\) 补齐为 \(2\) 的幂次。假设我们现在知道 \(f(x)g(x)\equiv1\pmod {x^{\frac{n}{2}}}\)，那么：

\[[f(x)g(x)-1]^2\equiv 0\pmod{x^n} \]

这是常见的折半再升高模数的技巧，有助于转化为子问题。

考虑展开，那么：

\[f^2(x)g^2(x)-2f(x)g(x)+1\equiv 0\pmod {x^n} \]

不妨设 \(x^n\) 下 \(f(x)\) 的逆是 \(g'(x)\)，那么显然应该左右两边各乘一个 \(g'(x)\) 来引入它：

\[\begin{aligned} &f(x)^2g(x)-2g(x)+g'(x)\equiv 0\pmod {x^n}\\ \iff&g'(x)\equiv 2g(x)-f(x)g^2(x)\pmod{x^n}\\ \iff&g'(x)\equiv g(x)(2-f(x)g(x))\pmod{x^n} \end{aligned} \]

那么我们直接倍增就好了。算出 \(g(x)\) 后进行两次乘法就可以获取 \(g'(x)\)。

struct modint{
	int x;
	static int Mod;
	static void setmod(int _){Mod=_;}
	static int getmod(){return Mod;}
	int qpow(int b,int p){
		if(!p) return 1;
		int d=qpow(b,p>>1);
		if(p&1) return 1ll*d*d%Mod*b%Mod;
		else return 1ll*d*d%Mod;
	}
	modint(int o=0){x=o;}
	modint &operator=(int o){return x=o,*this;}
	modint &operator+=(modint o){return x=(x+o.x>=Mod?x+o.x-Mod:x+o.x),*this;}
	modint &operator-=(modint o){return x=(x>=o.x?x-o.x:x-o.x+Mod),*this;}
	modint &operator*=(modint o){return x=1ll*x*o.x%Mod,*this;}
	modint &operator^=(modint o){return x=qpow(x,o.x),*this;}
	modint &operator/=(modint o){return *this*=(o^=(Mod-2));}
	modint &operator++() & {return (*this)+=1;}
	modint &operator--() & {return (*this)-=1;}
	modint &operator++(int) & {
		modint temp=*this;
		(*this)+=1;
		return temp;
	}
	modint &operator--(int) & {
		modint temp=*this;
		(*this)-=1;
		return temp;
	}
	friend modint operator+(modint a,modint b){return a+=b;}
	friend modint operator-(modint a,modint b){return a-=b;}
	friend modint operator*(modint a,modint b){return a*=b;}
	friend modint operator/(modint a,modint b){return a/=b;}
	friend modint operator^(modint a,modint b){return a^=b;}
	friend bool operator==(modint a,modint b){return a.x==b.x;}
	friend bool operator!=(modint a,modint b){return a.x!=b.x;}
	bool operator<(const modint &b) const{return x<b.x;}
	bool operator>(const modint &b) const{return x>b.x;}
	bool operator<=(const modint &b) const{return x<b.x||x==b.x;}
	bool operator>=(const modint &b) const{return x>b.x||x==b.x;}
	bool operator!(){return !x;}
	modint operator-(){return x?Mod-x:0;};
	friend std::istream &operator>>(std::istream &is,modint &k) {
		int val;
		is>>val;
		k=val;
		return is;
	}
	friend std::ostream &operator<<(std::ostream &os,modint k){return os<<k.x;}
};
int modint::Mod=998244353;
modint g,ig;
modint Find_proot(int p){
	vector <int> d;
	int lim=1;
	while(1){
		lim++;
		if(1ll*lim*lim>p-1) break;
		if((p-1)%lim==0) d.pb(lim),d.pb((p-1)/lim);
	}
	modint g=0;
	while(1){
		label:
		g++;
		for(auto i:d) if((g^((p-1)/i))==1) goto label;
		return g;
	}
}
#define VM vector <modint>
vector <int> pos;
void NTT(VM &A,int coef){
	int n=A.size();
	int mod=modint::getmod();
	fr1(i,0,n-1) if(i<pos[i]) swap(A[i],A[pos[i]]);
	for(int p=2,q=1;p<=n;p<<=1,q<<=1){
		modint wp=(coef?g:ig)^((mod-1)/p);
		for(int i=0;i<n;i+=p){
			modint wn=1;
			fr1(j,i,i+q-1){
				modint x=A[j],y=wn*A[j+q];
				A[j]=x+y;
				A[j+q]=x-y;
				wn*=wp;
			}
		}
	}
	if(!coef){
		modint x=1/(modint)n;
		fr1(i,0,n-1) A[i]*=x;
	}
}
VM multi(VM A,VM B){
	int a=A.size(),b=B.size();
	int t=a+b-1;
	int st=t;
	while(t!=lowbit(t)) t++;
	int lg=__lg(t);
	pos.resize(t);
	fr1(i,1,t-1) pos[i]=(pos[i>>1]>>1)|((i&1)<<lg-1);
	A.resize(t),B.resize(t);
	NTT(A,1),NTT(B,1);
	fr1(i,0,t-1) A[i]*=B[i];
	NTT(A,0);
	A.resize(st);
	return A;
}
VM inv(VM &F){
	VM G,nowF;
	int t=F.size();
	while(t!=lowbit(t)) t++;
	F.resize(t),G.resize(1);
	nowF.reserve(t);
	G[0]=1/F[0];
	nowF.pb(F[0]),nowF.pb(F[1]);
	for(int p=2;p<=t;p<<=1){
		VM temp=multi(nowF,G);
		temp.resize(p);
		fr1(i,0,p-1) temp[i]=-temp[i];
		temp[0]+=2;
		G=multi(G,temp);
		G.resize(p);
		if(p!=t) fr1(i,p,(p<<1)-1) nowF.pb(F[i]);
	}
	return G;
}
void init(int M){
	modint::setmod(M);
	g=Find_proot(M);
	ig=1/g;
}
/*
Remember to init()! Remember to init()! Remember to init()! Remember to init()! Remember to init()!
*/
struct Poly{
	VM p;
	Poly(){}
	Poly(int n){p.resize(n);}
	Poly(VM &a){p=a;}
	Poly(const std::initializer_list <modint> &a){p=a;}
	void push_back(modint x){p.pb(x);}
	int size(){return p.size();}
	Poly &operator*=(Poly o){return p=multi(p,o.p),*this;}
	Poly &operator/=(Poly o){return p=multi(p,inv(o.p)),*this;}
	friend Poly operator*(Poly a,Poly b){return a*=b;}
	friend Poly operator/(Poly a,Poly b){return a/=b;}
};

posted @ 2025-01-08 10:44 Shunpower 阅读(11) 评论(0) 收藏举报

刷新页面返回顶部

Shunpower 的小顺摩尔庄园 o(￣▽￣)ブ

May the force be with you and me.

多项式全家桶

FFT

FFT

IFFT

蝴蝶变换

NTT

多项式全家桶

多项式求逆

公告

Shunpower 的小顺摩尔庄园 o(*￣▽￣*)ブ

May the force be with you and me.

多项式全家桶

FFT

FFT

IFFT

蝴蝶变换

NTT

多项式全家桶

多项式求逆

公告

Shunpower 的小顺摩尔庄园 o(￣▽￣)ブ