多项式半家桶模板

多项式半家桶

多项式乘法

考虑两个多项式 $ F(x)=\sum_{i=0}^{i<n} A_i x_i \ \ , G(x)=\sum_{i=0}^{i<m} B_i x_i $ ,求 $ F(x) \times G(x) $。

暴力做,复杂度为 $ O(NM) $ ,不够优美。我们有 $ O((N+M)\log(N+M)) $ 的做法。

一个 \(n\) 次多项式可以由 \(n+1\) 个互不相同的点值表示,即点值表示法,取 \(N+M\) 个点的点值分别得到两个多项式的点值表示法,那么两个多项式的乘积的点值表示法就是两者直接按位相乘,这一部分复杂度为 $ O(N+M) $ 。

那么这里的瓶颈在于快速求 $ F(x),G(x) $ 的点值和将点值转化为多项式。

快速傅里叶变换(FFT)

不妨认为 \(n\) 为2的整数次幂,不符的话就向高位补零。

前置知识:单位根

设 $ x^n=1 $ ,则在复数域下 \(x\)\(n\) 个解,第 \(k\) 个解记为第 \(k\)\(n\) 次单位根 $ w_n^k $ ,其值为 $ e^{i\frac{2k\pi}{n}} = \cos(\frac{2k\pi}{n}) + i \sin(\frac{2k\pi}{n}) $ 。可以在单位圆上直观的表示。

有如下性质:

  1. $ w^{k+\frac{n}{2}}_n = -w_n^k $
  2. $ ( w_{2n}^k )^2 = w_{2n}^{2k} = w_n^k $
  3. $ w_{\frac{n}{2}}^k = w_n^{2k} $

对于多项式 $ F(x)=\sum_{i=0}^{i<n} A_i x_i $ 按奇偶性分开,变成 $ F_0(x) = \sum_{i=0}^{i<\frac{n}{2}} A_{2i} x^{2i} , F_1(x) = x \times \sum_{i=0}^{i<\frac{n}{2}} A_{2i+1} x^{2i} $ ,进行分治求解,最后再合并。

假设我们已经得到 $ F_0,F_1 $ 的点值表示法。

\[F(x) = F_0(x^2) + x F_1(x^2) \]

对于 $ k < \frac{n}{2} $ ,有

\[ F(w_n^k) = F_0(w_n^{2k}) + w_n^k F_1(w_n^{2k}) \\ F(w_n^k) = F_0(w_{\frac{n}{2}}^k) + w_n^k F_1(w_{\frac{n}{2}}^k) \]

\[ F(w_n^{k+\frac{n}{2}}) = F_0(w_n^{2({k+\frac{n}{2}})}) + w_n^{k+\frac{n}{2}} F_1(w_n^{2({k+\frac{n}{2}})}) \\ F(w_n^{k+\frac{n}{2}}) = F_0(w_{\frac{n}{2}}^{k+\frac{n}{2}}) + w_n^{k+\frac{n}{2}} F_1(w_{\frac{n}{2}}^{k+\frac{n}{2}}) \\ F(w_n^{k+\frac{n}{2}}) = F_0(w_{\frac{n}{2}}^{k}) - w_n^k F_1(w_{\frac{n}{2}}^{k}) \]

这样就可以类似于线段树一样分治了。最后 $ F(x) $ 的第 \(i\) 项即为 $ x = w_n^i $ 对应的点值。

快速傅里叶逆变换

接下来考虑把点值转化成多项式的一般形式。

设 $ y_i = F(w_n^i) $ , 定义 $ c_k = \sum_{i=0}^{ i<n } y_i (w_n^{ -k })^i $

\[ c_k = \sum_{i=0}^{i<n} y_i (w_n^{-k})^i \\ c_k = \sum_{i=0}^{i<n} ( \sum_{j=0}^{j<n} a_j w_n^{ij} ) w_n^{-ik} \\ c_k = \sum_{j=0}^{j<n} ( a_j \times \sum_{i=0}^{i<n} w_n^{i(j-k)} ) \]

容易发现,当 $ j \neq k $ 时,有 $ \sum_{i=0}^{i<n} w_n^{i(j-k)} = 0 $ ,因此

\[c_k = n \times a_k \]

由上述定义可知,$ c_k $ 就是多项式 $ H(x) = \sum_{i=0}^{i<n} y_i x^i $ 在 $ w_n^{-k} $ 时的点值。

所以只要把得到的点值再做一遍FFT即可。

递归版实现

CODE
struct comp{
	double x,y;
	comp operator + (comp k) {return {x+k.x,y+k.y};}
	comp operator - (comp k) {return {x-k.x,y-k.y};}
	comp operator * (comp k) {return {x*k.x-y*k.y,x*k.y+y*k.x};}
}f[N],g[N],h[N];
void FFT(comp *x,int len,int tp){
	if(len==1)return ;
	for(int i=0;i<(len>>1);i++)h[i]=x[i<<1],h[i+(len>>1)]=x[i<<1|1];
	for(int i=0;i<len;i++)x[i]=h[i];
	FFT(x,len>>1,tp);
	FFT(x+(len>>1),len>>1,tp);
	comp wn={cos(pi*2/len),sin(pi*2/len)*tp},w={1,0};
	for(int i=0;i<len;i++)h[i]=x[i];
	for(int i=0;i<(len>>1);i++,w=w*wn){
		comp t=w*h[i+(len>>1)];
		x[i]=h[i]+t;
		x[i+(len>>1)]=h[i]-t;
	}
}
	while(len<=n+m)len<<=1;
	FFT(f,len,1);
	FFT(g,len,1);
	for(int i=0;i<len;i++)f[i]=f[i]*g[i];
	FFT(f,len,-1);
	for(int i=0;i<=n+m;i++)printf("%d ",int(f[i].x/len+0.5));

注意数组大小和精度误差。千万不要忘记除多项式长度

迭代版实现

从底层往上层迭代,复杂度不变,但是可以有效减小常数。

观察一下可以发现,第 \(i\) 项多项式系数在底层的位置是它的二进制翻转。(卡常以后再说)

CODE
void init(){
	int p=0,d=(len>>1);
	tax[p++]=0; tax[p++]=d;
	for(int i=2;i<len;i<<=1){
		d>>=1;
		for(int j=0;j<i;j++)tax[p++]=(tax[j]|d);
	}
}
void FFT(comp *x,int tp){
	for(int i=0;i<len;i++){
		if(tax[i]>i)swap(x[i],x[tax[i]]);
	}
	for(int i=2;i<=len;i<<=1){
		for(int j=i;j<=len;j+=i){
			comp wn={cos(pi*2/i),sin(pi*2/i)*tp},w={1,0};
			for(int k=0;k<(i>>1);k++,w=w*wn){
				comp t=w*x[k+(i>>1)+j-i];
				h[k+j-i]=x[k+j-i]+t;
				h[k+(i>>1)+j-i]=x[k+j-i]-t;
			}
			for(int k=j-i;k<j;k++)x[k]=h[k];
		}
	}
}

快速数论变换(NTT)

这个似乎更常用一点。

前置知识:原根

满足 $ a^{n} \equiv 1 \ \ (\mod m) $ 的最小正整数 \(n\) ,称为 \(a\)\(m\) 的阶,记作 $ \delta_m(a) $ 。

已知欧拉定理:$ gcd(a,m) = 1 , a^{\varphi(m)} \equiv 1 \ \ (\mod m) $ 。

正整数 \(a\)\(m\) 的原根,当且仅当 $ \delta_m(a) = \varphi(m) $

主要用到这一性质:$ a^1 , a^2 , \dots , a^{ \delta_m(a) } $ 在模 \(m\) 的意义下互不相同。(具体证明可以看 oi-wiki

建议掌握一下求原根的方法。

如果你已经学会FFT,那NTT就很简单,直接拿原根替换单位根即可。避免了精度误差而且常数更小。

模数一般取 \(998244353\) ,最小的原根是 \(3\) 。(稍后会说明为什么常用)

CODE
namespace Poly{
	int len,tax[N],L;
	ll h[N],dw[2][N];
	void init(){
		for(int i=0;i<len;i++)
			tax[i]=((tax[i>>1]>>1)|((i&1)<<(L-1)));
		dw[1][1]=dw[0][1]=1;
		for(int i=2;i<len;i<<=1){
			ll w1=qpow(G,(mod-1)/i/2),w0=qpow(Gi,(mod-1)/i/2);//?
			dw[1][i]=1;dw[0][i]=1;
			for(int j=1;j<i;j++){
				dw[1][i+j]=dw[1][i+j-1]*w1%mod;
				dw[0][i+j]=dw[0][i+j-1]*w0%mod;
			}
		}
	}
	void NTT(ll *x,int tp){
		static ull pa[N];
		for(int i=0;i<len;i++)pa[i]=x[tax[i]];
		for(int i=2;i<=len;i<<=1){
			if((i>>18)&1){
				for(int j=0;j<len;j++)pa[j]%=mod;
			}
			for(int j=i,bg=0;j<=len;j+=i,bg+=i){
				for(int k=bg;k<bg+(i>>1);k++){
					ull ch=pa[k+(i>>1)]*dw[tp][(i>>1)+k-bg]%mod;
					pa[k+(i>>1)]=pa[k]+mod-ch;
					pa[k]+=ch;
				}
			}
		}
		for(int i=0;i<len;i++)x[i]=pa[i]%mod;
	}
	void Mul(ll *aas,ll *fir,int lx,ll *sec,int ly){
		static ll pa[N],pb[N];
		memcpy(pa,fir,sizeof(ll)*lx);
		memcpy(pb,sec,sizeof(ll)*ly);
		len=1;L=0;
		while(len<lx+ly-1)len<<=1,L++;
		memset(pa+lx,0,sizeof(ll)*(len-lx));
		memset(pb+ly,0,sizeof(ll)*(len-ly));
		init();
		NTT(pa,1);NTT(pb,1);
		for(int i=0;i<len;i++)aas[i]=pa[i]*pb[i]%mod;
		NTT(aas,0);
		ll k=qpow(len,mod-2);
		for(int i=0;i<len;i++)aas[i]=aas[i]*k%mod;
	}
}

注意到代码中打问号的一行,不知道大家有没有疑惑,如果不能整除怎么办。然而事实是,如果真不能整除,那就得换模数。但一般不会这样,因为我们惊奇地发现:$ 998244352 = 2^{23} \times 7 \times 17 $ ,大多数数据都不会卡满所有2因子。

多项式乘法逆

给定 $ G(x) $ ,求 \(F(x)\) 使得 $ F(x) G(x) \equiv 1 \ \ (\mod x^n) $ 。保证 $ [x^0] G(x) \neq 0 $

倍增求解,不妨认为我们已经求出模 $ x^{\lceil \frac{n}{2} \rceil} $ 意义下的乘法逆。

设 $ H(x) G(x) \equiv 1 \ \ (\mod x^{\lceil \frac{n}{2} \rceil}) $

\[ F(x) G(x) \equiv H(x) G(x) \mod x^{\lceil \frac{n}{2} \rceil} \\ G(x) ( F(x) - H(x) ) \equiv 0 \mod x^{\lceil \frac{n}{2} \rceil} \\ ( F(x) - H(x) )^2 \equiv 0 \mod x^n \\ F^2(x) + H^2(x) - 2 F(x) H(x) \equiv 0 \mod x^n \\ F(x) + G(x) H^2(x) - 2 H(x) \equiv 0 \mod x^n \\ F(x) \equiv 2 H(x) - G(x) H^2(x) \mod x^n \]

可以用牛顿迭代推

多项式开根

给定 $ G(x) $ ,求 $ F(x) $ 使得 $ F^2(x) \equiv G(x) \mod x^n $

设 $ H^2(x) \equiv G(x) \mod x^{\lceil \frac{n}{2} \rceil} $

\[ F(x) - H(x) \equiv 0 \mod x^{\lceil \frac{n}{2} \rceil} \\ F^2(x) - 2 F(x) H(x) + H^2(x) \equiv 0 \mod x^n \\ G(x) - 2 F(x) H(x) + H^2(x) \equiv 0 \mod x^n \\ F(x) \equiv \frac{ G(x) + H^2(x) }{ 2 H(x) } \mod x^n \]

需要用到多项式求逆。

可以用牛顿迭代推

多项式对数函数(多项式ln)

求 $ F(x) $ 使得 $ F(x) \equiv \ln G(x) \mod x^n $

前置知识:基础微积分

只能讲一些基本的概念和用法,建议找几篇博客学习,效果可能更好。

微分,在我的理解中,就是对可导函数求导,表示的是函数的瞬时变化率。(如果你会导数就可以跳过了)

举个例子,设 $ f(x) = x^2 $ , \(dx\) 为一极小量, $ f(x) $ 的导数 $ f'(x) = \frac{ ( x+dx )^2 - x^2 }{dx} = \frac{ 2 x \ dx + dx^2 }{dx} = \lim_{ dx \to 0 } ( 2 x + dx ) = 2x $ 。同理,若 $ f(x) = x^3 $ ,则有 $ f'(x) = 3 x^2 $

给出一些常见函数的导数:

$ ( k x )' = k $

$ ( x^a )' = a \times x^{ a-1 } $

$ ( \ln x )' = \frac1{x} $

$ ( a^x )' = a^x \times \ln a $

导数符合以下运算法则:

  1. $ ( f(x) + g(x) )' = f'(x) + g'(x) $
  2. $ ( f(x) \ g(x) )' = f'(x) g(x) + f(x) g'(x) $
  3. $ f( g(x) )' = g'(x) \ f'( g(x) ) $

至于积分,就是微分的逆运算,在函数图像上理解为面积。

对于一个多项式 $ F(x) = \sum A_i x_i $ 来说, $ F'(x) = \sum i A_i x_{i-1} $ 。对其逆向运算即可得到原多项式。

已知 $ F(x) \equiv \ln G(x) \mod x^n $ ,求导得到

\[ F'(x) \equiv ( \ln G(x) )' \mod x^n \\ F'(x) \equiv \frac{ G'(x) }{ G(x) } \mod x^n \]

最后再对 $ F'(x) $ 积分一下就做完了。

可以用牛顿迭代推

多项式指数函数(多项式exp)

前置知识:多项式牛顿迭代

首先介绍泰勒展开:

\[f(x) = \sum_{i=0} \frac{ f^{(i)}(x_0) }{i!} ( x-x_0 )^i \]

其中 $ f^{(i)} $ 为f的 \(i\) 阶导数,上式称为 $ f(x) $ 在 $ x_0 $ 处的泰勒展开。

现在考虑如下问题:求一个多项式 $ f(x) $ ,使得函数 $ g(f(x)) \equiv 0 \mod x^n $

我们倍增地去解决这个问题,设已经得到模 $ x^{ \frac{n}{2} } $ 的解 $ f_0 $ ,对 $ g(f(x)) $ 在 $ f_0 $ 处泰勒展开,得到 $$ g(f(x)) = \sum_{i=0} \frac{ g^{(i)}(f_0) }{i!} ( f(x)-f_0 )^i \mod x^n $$

假设存在一个 $ h(x) $ 满足 $ f(x) = f_0 + x^{\frac{n}{2}} \cdot h(x) $ 。

代入得到:

\[g(f(x)) = \sum_{i=0} \frac{ g^{(i)}(f_0) }{i!} ( x^{\frac{n}{2}} \cdot h(x) )^i \mod x^n \]

注意到,当 $ i \ge 2 $ 时 ,$ ( x^{\frac{n}{2}} \cdot h(x) )^i \equiv 0 \mod x^n $ ,化简可知

\[g(f(x)) = g( f_0 ) + x^{\frac{n}{2}} \cdot g'(f_0) \cdot h(x) = 0 \]

\[h(x) = - \frac{ g(f_0) }{ x^{\frac{n}{2}} \cdot g'(f_0) } \]

代入 $ f(x) = f_0 + x^{\frac{n}{2}} \cdot h(x) $

\[f(x) = f_0 - \frac{ g(f_0) }{ g'(f_0) } \mod x^n \]

至于 $ h(x) $ 为什么总是存在,先感受一下。

求 $ F(x) \equiv e^{ G(x) } \mod x^n $

\[\ln F(x) - G(x) \equiv 0 \mod x^n \]

代入多项式牛顿迭代

\[F(x) = F_0(x) - \frac{ \ln F(x) - G(x) }{ ( \ln F(x) - G(x) )' } \mod x^n \]

\[F(x) = F_0(x) - ( \ln F(x) - G(x) ) \cdot F(x) \mod x^n \]

好玩吧

模板

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=8e5+100;
const ll mod=998244353;
const ll G=3,Gi=332748118;
int m,n;
inline int read(){
	int k=0;
	char ch=getchar();
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9'){
		k=k*10+ch-'0';
		ch=getchar();
	}
	return k;
}
inline void write(int x){
	if(!x)return putchar('0'),void();
	int tmp[20],len=0;
	while(x)tmp[++len]=x%10,x/=10;
	for(int i=len;i>=1;i--)putchar(char(tmp[i]+'0'));
}
inline 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;
}
ll f[N],t[N],ny[N];
namespace Poly{
	int len,tax[N],L;
	ll h[N],dw[2][N];
	void init(){
		for(int i=0;i<len;i++)
			tax[i]=((tax[i>>1]>>1)|((i&1)<<(L-1)));
		dw[1][1]=dw[0][1]=1;
		for(int i=2;i<len;i<<=1){
			ll w1=qpow(G,(mod-1)/i/2),w0=qpow(Gi,(mod-1)/i/2);
			dw[1][i]=1;dw[0][i]=1;
			for(int j=1;j<i;j++){
				dw[1][i+j]=dw[1][i+j-1]*w1%mod;
				dw[0][i+j]=dw[0][i+j-1]*w0%mod;
			}
		}
	}
	void NTT(ll *x,int tp){
		static ull pa[N];
		for(int i=0;i<len;i++)pa[i]=x[tax[i]];
		for(int i=2;i<=len;i<<=1){
			if((i>>18)&1){
				for(int j=0;j<len;j++)pa[j]%=mod;
			}
			for(int j=i,bg=0;j<=len;j+=i,bg+=i){
				for(int k=bg;k<bg+(i>>1);k++){
					ull ch=pa[k+(i>>1)]*dw[tp][(i>>1)+k-bg]%mod;
					pa[k+(i>>1)]=pa[k]+mod-ch;
					pa[k]+=ch;
				}
			}
		}
		for(int i=0;i<len;i++)x[i]=pa[i]%mod;
	}
	void Mul(ll *aas,ll *fir,int lx,ll *sec,int ly){
		static ll pa[N],pb[N];
		memcpy(pa,fir,sizeof(ll)*lx);
		memcpy(pb,sec,sizeof(ll)*ly);
		len=1;L=0;
		while(len<lx+ly-1)len<<=1,L++;
		memset(pa+lx,0,sizeof(ll)*(len-lx));
		memset(pb+ly,0,sizeof(ll)*(len-ly));
		init();
		NTT(pa,1);NTT(pb,1);
		for(int i=0;i<len;i++)aas[i]=pa[i]*pb[i]%mod;
		NTT(aas,0);
		ll k=qpow(len,mod-2);
		for(int i=0;i<len;i++)aas[i]=aas[i]*k%mod;
	}
	inline void Lim(ll *x,int lx,int ln){
		if(lx>ln)memset(x+ln,0,sizeof(ll)*(lx-ln));
	}
	inline void Add(ll *x,ll *y,int ln){
		for(int i=0;i<ln;i++)x[i]=(x[i]+y[i])%mod;
	}
	inline void Red(ll *x,ll *y,int ln){
		for(int i=0;i<ln;i++)x[i]=(x[i]+mod-y[i])%mod;
	}
	inline void mul(ll *x,ll k,int ln){
		for(int i=0;i<ln;i++)x[i]=x[i]*k%mod;
	}
	inline void Der(ll *x,ll *y,int ln){
		for(int i=1;i<ln;i++)x[i-1]=y[i]*i%mod;
	}
	inline void Int(ll *x,ll *y,int ln){
		for(int i=ln;i>0;i--)x[i]=y[i-1]*ny[i]%mod;
		x[0]=0;
	}
	void Inv(ll *aas,ll *x,int lx){
		assert(x[0]);
		aas[0]=qpow(x[0],mod-2);
		static ll pa[N],pb[N];
		memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
		for(int i=1;i<lx;i<<=1){
			memcpy(pa,x,sizeof(ll)*(i<<1));
			memcpy(pb,aas,sizeof(ll)*i);
			Mul(pa,pa,i<<1,pb,i<<1);Lim(pa,len,i<<1);
			Mul(pa,pa,i<<1,pb,i<<1);Lim(pa,len,i<<1);
			mul(pb,2,i);Red(pb,pa,i<<1);
			memcpy(aas,pb,sizeof(ll)*(i<<1));
		}
	}
	void Sqrt(ll *aas,ll *x,int lx){
		aas[0]=1;
		static ll pa[N],pb[N],pc[N];
		memset(pa,0,sizeof(pa));
		memset(pb,0,sizeof(pb));
		memset(pc,0,sizeof(pc));
		for(int i=1;i<lx*2;i<<=1){
			memcpy(pa,x,sizeof(ll)*(i<<1));
			memcpy(pb,aas,sizeof(ll)*i);
			Inv(pc,pb,i);
			Mul(pc,pa,i<<1,pc,i);Lim(pc,i<<2,i<<1);
			Add(pc,pb,i<<1);mul(pc,qpow(2,mod-2),i<<1);
			memcpy(aas,pc,sizeof(ll)*(i<<1));
		}
	}
	void Ln(ll *aas,ll *x,int lx){
		static ll pa[N],pb[N];
		memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
		Der(pa,x,lx);Inv(pb,x,lx);
		Mul(aas,pa,lx,pb,lx);
		Int(aas,aas,len);
	}
	void Exp(ll *aas,ll *x,int lx){
		aas[0]=1;
		static ll pa[N],pb[N];
		memset(pa,0,sizeof(pa));memset(pb,0,sizeof(pb));
		for(int i=1;i<lx*2;i<<=1){
			memcpy(pa,aas,sizeof(ll)*i);
			memcpy(pb,x,sizeof(ll)*(i<<1));
			Ln(pa,aas,i);
			Red(pb,pa,i<<1);pb[0]++;
			Mul(aas,aas,i,pb,i<<1);Lim(aas,len,i<<1);
		}
	}
}

关于常数

说真的,最开始的时候被多项式的常数震撼到了。写多项式开根的时候实现不精细直接T飞了,以上模板是卡常后的结果,跑的还是不快。

卡常可以看九日的博客 。实现好的多项式板子效率还是可以接受的。

posted @ 2025-06-20 21:51  Abnormal123  阅读(27)  评论(0)    收藏  举报