多项式简陋入门

多项式全家桶

然而并没有多点求值,快速插值,转下降/上升幂,复合,复合逆

疯狂多项式,v我50

namespace efX_poly
{
	const int maxlen=(1<<23)+1,maxSqrt=1e5+1;
	inline int add(int x,int y,int m) { x+=y; return x>=m?x-m:x; }
	inline int sub(int x,int y,int m) { x-=y; return x<0?x+m:x; }
	inline int mul(int x,int y,int m) { ll z=1ll*x*y; return z>=m?z%m:z; }
	inline void Add(int &x,int y,int m) { x=add(x,y,m); }
	inline void Sub(int &x,int y,int m) { x=sub(x,y,m); }
	inline void Mul(int &x,int y,int m) { x=mul(x,y,m); }
	inline int qpow(int a,int x,int m)
	{
		int r=1;
		while(x)
		{
			if(x&1) Mul(r,a,m);
			Mul(a,a,m),x>>=1;
		} return r;
	} inline int _inv(int a,int m) { return qpow(a,m-2,m); }
	namespace set_mod
	{
		static vector<int> d;
		inline void init(int p)
		{
			int tmp=p-1;
			if(!d.empty()) d.clear();
			for(int i=2;1ll*i*i<=tmp;++i)
				if(tmp%i==0)
			{
				d.emplace_back(i);
				while(tmp%i==0) tmp/=i;
			}
			if(tmp>1) d.emplace_back(tmp);
		}
		inline bool check(int x,int p)
		{
			for(int v:d)
				if(qpow(x,(p-1)/v,p)==1)
					return false;
			return true;
		}
		inline int getg(int p)
		{
			init(p);
			irep(i,2,p)
			if(check(i,p))
				return i;
			return -1;
		}
	}
	int r[maxlen];
	inline void make_rev(int lim,int l)
	{
		static int his=-1;
		if(his==lim) return;
		for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
		his=lim;
	}
	inline pii getLimit(int n)
	{
		int lim=1,l=0;
		while(lim<n) lim<<=1,++l;
		return make_pair(lim,l);
	}
	
	int inv[maxlen];
	int tmp[maxlen],tmpb[maxlen];
	inline void view_arr(int *a,int lim) { cerr<<"## "; for(int i=0;i<lim;++i) {cerr<<a[i]<<' ';} cerr_out(); }
	class poly
	{
	public:
		const int p,g,gi;
	private:
		int G[2][maxSqrt+1],GI[2][maxSqrt+1];
		inline void init()
		{
			memset(inv,0,sizeof inv);//can't use it any more
			G[0][0]=GI[0][0]=1;
			for(int i=1;i<=maxSqrt;++i)
			{
				G[0][i]=mul(G[0][i-1],g,p);
				GI[0][i]=mul(GI[0][i-1],gi,p);
			}
			G[1][0]=GI[1][0]=1;
			for(int i=1;i<=maxSqrt;++i)
			{
				G[1][i]=mul(G[1][i-1],G[0][maxSqrt],p);
				GI[1][i]=mul(GI[1][i-1],GI[0][maxSqrt],p);
			}
		}
		inline int very_quick_pow(bool o,int x)
		{
			if(o) return mul(G[1][x/maxSqrt],G[0][x%maxSqrt],p);
			return mul(GI[1][x/maxSqrt],GI[0][x%maxSqrt],p);
		}
	public:
		inline void NTT(int *a,int lim,int op)
		{
			for(int i=0;i<lim;++i) if(i<r[i]) swap(a[i],a[r[i]]);
			for(int mid=1,wn;mid<lim;mid<<=1)
			{
				wn=very_quick_pow(op==1,(p-1)/(mid<<1));
				for(int j=0,w;j<lim;j+=(mid<<1))
				{
					w=1;
					for(int k=0,x,y;k<mid;++k,w=mul(w,wn,p))
						x=a[j+k],
						y=mul(w,a[j+k+mid],p),
						a[j+k]=add(x,y,p),
						a[j+k+mid]=sub(x,y,p);
				}
			}
			if(op==-1)
			{
				int iv=_inv(lim,p);
				for(int i=0;i<lim;++i) Mul(a[i],iv,p);
			}
		}
		inline void clear(int *a,int st,int ed) { for(int i=st;i<ed;++i) a[i]=0; }
		inline void arrSub(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Sub(a[i],b[i],p); }
		inline void arrSub(int *a,int c,int lim) { for(int i=0;i<lim;++i) Sub(a[i],c,p); }
		inline void arrAdd(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Add(a[i],b[i],p); }
		inline void arrAdd(int *a,int c,int lim) { for(int i=0;i<lim;++i) Add(a[i],c,p); }
		inline void Reverse(int *a,int *b,int n) { for(int i=0;i<n;++i) b[i]=a[n-i-1]; }
		inline void arrMul(int *a,int *b,int lim) { for(int i=0;i<lim;++i) Mul(a[i],b[i],p); }
		inline void arrMul(int *a,int c,int lim) { for(int i=0;i<lim;++i) Mul(a[i],c,p); }
		inline void Inv(int *a,int *b,int n)
		{
			static int tmp[maxlen];//开成全局就gg了
			if(n==1) { b[0]=_inv(a[0],p); return; }
			Inv(a,b,(n+1)>>1);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memcpy(tmp,a,sizeof(int)*n);
			for(int i=n;i<lim;++i) tmp[i]=0;
			NTT(tmp,lim,1),NTT(b,lim,1);
			for(int i=0;i<lim;++i) b[i]=mul(sub(2,mul(b[i],tmp[i],p),p),b[i],p);
			NTT(b,lim,-1);//tmp就不用管了,以后用的时候会刷新的
			for(int i=n;i<lim;++i) b[i]=0;
		}
		inline void Sqrt(int *a,int *b,int n)
		{
			if(n==1) { b[0]=1; return; }
			if(!inv[2]) inv[2]=_inv(2,p);
			Sqrt(a,b,(n+1)>>1);
			memset(tmp,0,sizeof(int)*(n<<1));
			Inv(b,tmp,n);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memcpy(tmpb,a,sizeof(int)*n);
			clear(tmpb,n,lim);
			NTT(tmpb,lim,1);
			NTT(tmp,lim,1);
			NTT(b,lim,1);
			for(int i=0;i<lim;++i)
				Add(b[i],mul(tmpb[i],tmp[i],p),p),
				Mul(b[i],inv[2],p);
			NTT(b,lim,-1);
			clear(b,n,lim);
		}//be very very careful!
		inline void Int(int *a,int *b,int n)
		{
			for(int i=n-1;i;--i)
			{
				if(!inv[i]) inv[i]=_inv(i,p);
				b[i]=mul(a[i-1],inv[i],p);
			} b[0]=0;
		}
		inline void Der(int *a,int *b,int n)
		{
			for(int i=1;i<n;++i) b[i-1]=mul(a[i],i,p); 
			b[n-1]=0;//attention to the bound
		}
		inline void Ln(int *a,int *b,int n)
		{
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			memset(tmp,0,sizeof(int)*lim);
			memset(tmpb,0,sizeof(int)*lim);
			memcpy(tmp,a,sizeof(int)*n);
			Der(a,tmp,n);
			Inv(a,tmpb,n);
			NTT(tmp,lim,1),NTT(tmpb,lim,1);
			for(int i=0;i<lim;++i) tmp[i]=mul(tmp[i],tmpb[i],p);
			NTT(tmp,lim,-1);
			Int(tmp,tmp,n);
			memcpy(b,tmp,sizeof(int)*n);
		}
		inline void Exp(int *a,int *b,int n)
		{
			if(n==1) { b[0]=1; return; }
			Exp(a,b,(n+1)>>1);
			Ln(b,tmp,n);
			int l=0,lim=1;
			while(lim<(n<<1)) lim<<=1,++l;
			make_rev(lim,l);
			for(int i=0;i<n;++i) tmp[i]=sub(a[i],tmp[i],p);
			clear(tmp,n,lim);
			clear(b,n,lim);
			tmp[0]++;
			NTT(tmp,lim,1);
			NTT(b,lim,1);
			arrMul(b,tmp,lim);
			NTT(b,lim,-1);
			clear(b,n,lim);
		}
		inline void Pow(int *a,int k,int *b,int n)
		{
			static int tmp[maxlen];
			memset(tmp,0,sizeof(int)*(n<<1));
			Ln(a,tmp,n);
			arrMul(tmp,k,n);
			Exp(tmp,b,n);
		}//it's not totally safe
		inline void Mod(int *a,int *b,int *P,int *Q,int n,int m)
		{
			--n,--m;//强行钦定使用左闭右开区间
			static int ar[maxlen],br[maxlen];
			memset(ar,0,sizeof(int)*(n));
			memset(br,0,sizeof(int)*(n));
			Reverse(a,ar,n+1);
			Reverse(b,br,m+1);
			irep(i,n-m+1,m) br[i]=0;
			Inv(br,P,n-m+1);
			pii lm=getLimit(2*n-m+1);//ar*p=br^-1
			make_rev(lm.first,lm.second);
			NTT(ar,lm.first,1),
			NTT(P,lm.first,1);
			arrMul(ar,P,lm.first);
			NTT(ar,lm.first,-1);
			Reverse(ar,P,n-m+1);//所以p不需要ntt回来
			lm=getLimit(n+1);
			make_rev(lm.first,lm.second);
			clear(P,n-m+1,lm.first);
			memcpy(ar,b,sizeof(int)*lm.first);
			NTT(P,lm.first,1),
			NTT(ar,lm.first,1);
			arrMul(ar,P,lm.first);
			NTT(P,lm.first,-1);//p is needed
			NTT(ar,lm.first,-1);
			for(int i=0;i<lm.first;++i) Q[i]=sub(a[i],ar[i],p);
			clear(Q,m,lm.first);
		}//[a]=n,[b]=m,a%b=p...q
		poly():p(998244353),g(3),gi(332748118) { init(); }
		poly(int mod):p(mod),g(set_mod::getg(mod)),gi(_inv(g,mod)) { init(); }
	};
}

\(\text{poly}\) 提供了一个变幻模数的接口,不过使得常数变大了一倍,如果没有必要写任意模数的话应该拆开

函数的规范是 传入参数,传出参数,数组长度

比如求 \(g(x)=f^-1(x)\pmod {x^n}\) 就是 Inv(f,g,n),不难注意到这样写我们的数组是左闭右开的,事实上,我们要求所有函数这样(所以取模里面硬减了一下来

更详细的约定是

  • NTT(a,b,c)a为需要变幻的多项式,b表示次数,c\(-1\) 时进行逆变换,其他时候进行正变幻

  • Inv(a,b,c)a为需要求逆的多项式,b表示逆元(在 \(\mathbb{Z}_p\) 的意义下),c为次数,之后默认最后一个参数表示数列长度(也意味着是在 \(\bmod {x^n}\) 的意义下)

  • clear(a,b,c),清空a\([l,r)\) 的值

  • arrAdd(a,b,c),如果b为数列,则将对应位相加,如果b为常数,则a的每一位加上b

  • arrSub(a,b,c),同上,加法变为减法

  • arrMul(a,b,c),同上,加法变为乘法

  • Sqrt(a,b,c),求解 \(b\equiv\sqrt a\pmod{x^n}\)

  • Der(a,b,c),求解 \(b=a^\prime\)

  • Int(a,b,c),求解 \(b=\int a\),常数项置 \(0\)

  • Ln(a,b,c),求解 \(b=\ln a\),应当有 \(a_0=1\)

  • Exp(a,b,c),求解 \(b=e^x\),应当有 \(a_0=0\)

  • Pow(a,k,b,c),求解 \(b=a^k\),应当有 \(a_0=1\)

  • Mod(a,b,p,q,n,m)a的长度为 \(n\)b的长度为 \(m\),求解 \(c=\lfloor\dfrac{a}{b}\rfloor,d=a\bmod b\),就是小学学的大除法

原理都不难,主要就是多项式牛顿迭代

多项式牛迭是为了解决下面这个问题

\[G(F(x))\equiv 0\pmod{x^n} \]

我们期望用 \(G(x)\) 表示或者求出 \(F(x)\),当然,它们都是多项式

考虑一个倍增的做法

我们假设已经求出了 \(F_0(x)\) 满足

\[G(F_0(x))\equiv 0\pmod{x^{\frac{n}{2}}} \]

显然有

\[F_0(x)\equiv F(x)\pmod{x^{\frac n2}} \]

\(G(F(x))\) 看做关于 \(F(x)\) 的函数,在 \(F_0(x)\) 处展开

可得

\[\begin{align*} G(F(x))&=G(F_0(x))+\frac{G^\prime(F_0(x))}{1!}(F(x)-F_0(x))+\cdots\\ &=\sum_{i=0}^\infty\frac{G^{(i)}(F(x))}{i!}(F(x)-F_0(x))^i \end{align*} \]

我们只关心 \(F(x)\)\(\bmod x^n\) 的意义下的值,而显然有 \(F(x)-F_0(x)\equiv0\pmod{x^\frac{n}{2}}\)

也就是说,\((F(x)-F_0(x))^2\equiv0\pmod {x^n}\)

于是只需要保留泰勒展开的前两项,即有

\[F(x)=F_0(x)-\frac{G(F_0(x)}{G^\prime(F_0(x))} \]

最后 \(n=1\) 的边界一般都可以转化为数论问题,然后就可以倍增了

复杂度一般都是 \(O(n\log n)\),用主定理看一眼,都是 \(T(n)=O(n\log n)+T(\dfrac{n}{2})\),当然,还有一些这里没有的操作,有些就不是这个复杂度,如果你写半在线卷积或者分治 \(\text{NTT}\) 的话也不是

符号化方法

实例

或许先看一些例子对于理解理论是有好处的

直观的理解

我们期望用生成函数解决组合问题,于是需要关联生成函数和组合意义,假定需要求解的组合对象构成的集合为 \(\mathcal{A}\),我们定义一个函数把组合对象映射到非负整数,这样就完成了连接生成函数和组合概念的第一步

我们规定 \(\alpha\in\mathcal{A}\)\(|\alpha|\) 为其大小函数,表示从组合对象到自然数的映射

比如我们考虑求解 0/1 背包问题拿总重量为 \(V_0\) 的方案数。假设对于单个物品,它的重量为 \(V\),那么我们可以直接把重量当做这个物品的大小函数,如果没有选择,那么大小自然是 \(0\)。于是我们可以写出单个物品的生成函数 \(1+x^V\),考虑每个组合对象(这里就是每个物品),它们的组合类应该是 \(\mathcal{A}_i=1+x^{V_i}\),最后对于整个组合对象,也就是我们要求解的问题,应当有

\[\mathcal{A}=\mathcal{A}_1\times\mathcal{A}_2\times\cdots\times\mathcal{A}_n \]

其中,\(\times\) 表示笛卡尔积,可以感性的理解为把两个组合对象直接拼接起来,而对应的代数含义就是生成函数的卷积

引入中性对象 \(\xi=1\) 和元对象 \(\circ=x\),可以类比自然数中的 \(0\)\(1\)

形式化地,我们用 \(\mathcal{A}\cong \mathcal{B}\) 表示 \(\mathcal{A}\)\(\mathcal{B}\) 是同构的组合对象,显然我们有

\[\mathcal{A}\cong \mathcal{A}\times\xi\cong\xi\times \mathcal{A} \]

SEQ

终于开始有用的东西了...

我们考虑 \(\operatorname{SEQ}(S)\) 生成了包括 \(\varnothing\) 在内的所有用 \(S\) 中的元素构成的序列

比如 \(\text{SEQ}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,a\},\{a,b\},\{b,a\},\{b,b\}\cdots\}\)

很显然,我们可以发现 \(\text{SEQ}\) 只是把若干个集合中选出任意元素有序地拼接在一起,那么考虑把这个自然语言表述的过程符号化

\[\text{SEQ}(S)=1+S(x)+S(x)^2+S(x)^3+\cdots=\frac{1}{1-S(x)} \]

回想一下上面的过程,我们推导出 \(\text{SEQ}(\mathcal{A})\) 的封闭形式是 \(\dfrac{1}{1-\mathcal{A}(x)}\) 的过程只依赖于代数技巧,这启发我们这种函数可以应用于更多的场景,而不是局限于某种特定的组合概念(这也是符号化方法的优点,它简洁地指明了不同组合意义的共性)

求有序有根树的生成函数

树一定有根,儿子之间的顺序有意义

\[\mathcal{T}=\circ\times(\varnothing+\mathcal{T}+\mathcal{T}\times\mathcal{T}+\cdots)=\circ\times\text{SEQ}(\mathcal{T})\iff T(x)=\frac{x}{1-T(x)} \]

之后的处理方式有很多,个人比较推荐像链接里那样用拉格朗日反演

MSET

比起 \(\text{SEQ}\),无序的情形似乎更加常见,于是就有了 \(\text{MSET}(S)\),其包含所有可能的无序的组合

比如 \(\text{MSET}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,a\},\{a,b\},\{b,b\}\cdots\}\)

注意到 \(\text{SEQ}\) 中出现了 \(\{a,b\},\{b,a\}\),而在 \(\text{MSET}\) 中,它们被认为是一样的

为了求解 \(\text{MSET}\),我们考虑递推,递推就是用我们已经知道的去求解不知道的,我们注意到 \(\text{MSET}(\{a\})=\text{SEQ}(\{a\})=\{\varnothing,\{a\},\{a.a\},\cdots\}\)

于是可以有

\[\begin{align*} \text{MSET}(\{a_1,a_2,a_3\cdots a_n\})&=\text{MSET}(\{a_1,a_2,a_3\cdots a_{n-1}\})\times\text{SEQ}(\{a_n\})\\ &=\prod_{i=1}^n \text{SEQ}(\{a_i\})\\ &=\prod_{i=1}^n\frac{1}{1-x^{|a_i|}} \end{align*} \]

累乘的形式当然不能算作封闭形式,考虑高中数学常见的套路先取 \(\ln\) 再取 \(\exp\),设 \(c_i\) 表示 \(|a_k|=i\)\(a\) 的数量,那么

\[\begin{align*} \text{MSET}(\{a_1,a_2,a_3\cdots a_n\})&=e^{-\sum\limits_{i=1}^n\ln(1-x^{|a_i|})}\\ &=e^{-\sum\limits_{i\ge 1}c_i\ln(1-x^i)} \end{align*} \]

利用麦克劳林公式 \(\ln(1+z)=\sum\limits_{i\ge 1}\dfrac{(-1)^{i-1}z^i}{i}\)

\[\begin{align*} \text{MSET}(\{a_1,a_2,a_3\cdots a_n\})&=e^{-\sum\limits_{i\ge 1}c_i\ln(1-x^i)}\\ &=e^{-\sum\limits_{i\ge 1}c_i\sum\limits_{j\ge 1}\frac{-x^{ij}}{j}}\\ \end{align*} \]

此时,我们发现 \(S(x)=-\sum\limits_{i\ge 1}c_i\ln(1-x^i)\)

整理得

\[\text{MSET}(\mathcal{A}(x))=e^{\sum\limits_{i\ge 1}\frac{A(x^i)}{i}}=e^{\frac{A(x)}{1}+\frac{A(x^2)}{2}+\frac{A(x^3)}{3}+\cdots} \]

较为简单的例子

\(\text{MSET}\) 的一个应用实例是求解数的拆分,如果考虑用 DP 解决,比如 \(f_{i,j}\) 表示用到了 \(i\),和为 \(j\) 的方案数

那么不难列出 \(f_{i,j}=\sum f_{i-1,j-ki}\)

发现这和 \(\text{MSET}\) 的递推式没有区别,于是要求的就是 \(\text{MSET}(\mathbb{N})\)

更为直接的方式是直接写出目标的生成函数 \(A(x)=\prod\limits_{i\ge 1}\dfrac{1}{1-x^i}\),不过这样不够机械,有种组合意义的感觉

更为困难的例子

求大小为 \(n\) 的无编号无根树的数量 \(a_n\)

先求无编号有根树 \(\mathcal{T}\)

无编号有根树的子树也是无编号有根树,因为没有编号,所以儿子之间的顺序不重要,这完全等价于儿子间存在严格的偏序关系

于是

\[\mathcal{T}=\circ\times(\varnothing+\{\alpha\mid\alpha\in\mathcal{T}\}+\{(\alpha,\beta)\mid\alpha<\beta,\alpha,\beta\in\mathcal{T}\}+\cdots)=\circ\times\text{MSET}(\mathcal{T}) \]

论文中可知,无编号无根树的生成函数是 \(T(z)-\frac 12T(z)^2+\frac12T(z^2)\)

PSET

定义 \(\text{PSET}(\mathcal{S})\)\(\mathcal{S}\) 的所有子集构成的集合

比如 \(\text{PSET}(\{a,b\})=\{\varnothing,\{a\},\{b\},\{a,b\}\}\)

考虑 \(\text{PSET}(\{a\})=1+x^{|a|}\)

由定义知

\[\begin{align*} \text{PSET}(\mathcal{S})&=\prod_{\iota\in\mathcal{S}}\text{PSET}(\{\iota\})\\ &=\prod_{i\ge 1}(1+x^i)^{c_i}\\ &=e^{\sum_{i\ge 1}c_i\ln(1+x^i)}\\ &=e^{\sum_{i\ge 1}c_i\sum_{j\ge 1}\frac{(-1)^{j-1}x^{ij}}{j}}\\ &=e^{\frac{S(x)}{1}-\frac{S(x)^2}{2}+\frac{S(x)^3}{3}-\cdots}\\ &=e^{\sum_{i\ge1}\frac{(-1)^{i-1}S(x)^i}{i}} \end{align*} \]

CYC

我们定义 \(CYC(\mathcal{C})\) 表示用 \(\mathcal{C}\) 中元素可以表示的不同的环排列的集合,其中 \(\mathcal{C}\) 一般是可重集合

\[A(x)=\sum_{k\ge1}\frac{\varphi(k)}{k}\log\frac{1}{1-B(x^k)} \]

或者

\[\log(A(x))=\sum_{k\ge1}\frac{\varphi(k)}{k}\ln\frac{1}{1-A(x^i)} \]

证明有点难,我看懂了再补,咕咕咕咕咕咕

posted @ 2022-10-13 16:47  嘉年华_efX  阅读(24)  评论(0编辑  收藏  举报