多项式全家桶

众所周知,生成函数是一个十分强大的东西,许多与多项式相关的算法也就应运而生了,在这里选取几种较为简单的算法做一个介绍.

p.s. 这篇文章在去年NOI前已经完成了一半,现在笔者将其补充完整后发出,同时也为了纪念那一段美好的时光。

多项式乘法(FFT/NTT)

https://www.cnblogs.com/encodetalker/p/10212036.html

https://www.cnblogs.com/encodetalker/p/10285657.html

多项式求逆

已知\(F(x)\),求\(G(x)\)使得\(F(x)G(x)\equiv 1(mod\ x^n)\)

倍增,假设已经求出\(A(x)F(x)\equiv1(mod\ x^{\lceil \frac{n}{2}\rceil})\)

\(G(x)-A(x)\equiv0(mod\ x^{\lceil \frac{n}{2}\rceil})\)

注意到这个式子可以平方一下\(G(x)^2-2A(x)G(x)+A(x)^2\equiv0(mod\ x^n)\)
两边乘上\(F(x),G(x)-2A(x)+F(x)A(x)^2\equiv0(mod\ x^n)\)

最后得到\(G(x)\equiv A(x)(2-F(x)A(x))(mod\ x^n)\)

多项式取模

咕咕咕

多项式牛顿迭代

假设已知\(n\)次多项式\(F(x)\),求多项式\(G(x)\)使得\(F(G(x))\equiv 0 (\mod x^n)\).

同样考虑倍增,假设\(F(G_1(x))\equiv 0(\mod x^{\lceil \frac{n}{2}\rceil})\).

考虑将\(F(G(x))\)\(G_1(x)\)处Taylor展开,有:

\[F(G(x))=F(G_1(x))+\frac{F'(G(x))}{1!}(G(x)-G_1(x))+\frac{F^{(2)}(G(x)}{2!}(G(x)-G_1(x))^2+\cdots \]

注意到在\(\mod n\)的意义下,若\(k\geq 2\), 那么\((G(x)-G_1(x))^k\equiv 0(\mod x^n)\).

于是有

\[F(G(x))\equiv F(G_1(x))+F'(G(x))(G(x)-G_1(x))(\mod x^n) \]

移项之后得到:

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

多项式开方

已知\(n\)次多项式\(F(x)\),求多项式\(G(x)\)满足\(G(x)^2\equiv F(x)(\mod x^n)\).

考虑牛顿迭代,构造\(A(G(x))=G(x)^2-F(x).\)显然最后需要\(A(G(x))\equiv 0(\mod x^n)\).

带入牛顿迭代的结论中可以得到:

\[G(x)\equiv G_1(x)-\frac{G_1(x)^2-F(x)}{2G_1(x)}=\frac{G_1(x)^2+F(x)}{2G_1(x)}(\mod x^n) \]

多项式求导&&积分

\(F(x)=\sum_{i=0}^n a_ix^i\)

求导:\(F'(x)=\sum_{i=1}^na_iix^{i-1}\)

积分:\(\int F(x)=\sum_{i=0}^n(\frac{a_i}{i+1}x^{i+1})+C\)(一个常数)

多项式Ln

复合函数求导:记\(H(x)=F(G(x))\),则\(H'(x)=F'(G(x))G'(x)\)

已知多项式\(F(x)\),求\(G(x)=ln(F(x))\)

对两边求导:\(G'(x)=\frac{F'(x)}{F(x)}\)

之后再积分回去:\(G(x)=\int \frac{F'(x)}{F(x)}\),多项式求逆之后再积分即可。

多项式Exp

已知\(n\)次多项式\(F(x)\),求\(G(x)\)满足\(G(x)\equiv e^{F(x)}(\mod x^n)\).

两边同时取\(\ln\)\(\ln G(x)\equiv F(x)(\mod x^n)\).

\(A(x)=\ln G(x)-F(x)\), 继续套用牛顿迭代得到:

\[G(x)\equiv G_1(x)(1-\ln G_1(x)+F(x))(\mod x^n) \]

总代码如下:

namespace Polynomial{
	
	int A[N<<2],r[N<<2],B[N<<2],C[N<<2],D[N<<2];
	
	int calcr(int len)
	{
		int lim=1,cnt=0;
		while (lim<len) {lim<<=1;cnt++;}
		rep(i,0,lim-1)
			r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
		return lim;
	}
	
	void ntt(int lim,int *a,int typ)
	{
		rep(i,0,lim-1)
			if (i<r[i]) swap(a[i],a[r[i]]);
		for (int mid=1;mid<lim;mid<<=1)
		{
			int len=(mid<<1),gn=qpow(3,(maxd-1)/len);
			if (typ==-1) gn=getinv(gn);
			for (int sta=0;sta<lim;sta+=len)
			{
				int g=1;
				for (int j=0;j<mid;j++,g=mul(g,gn))
				{
					int x=a[sta+j],y=mul(a[sta+j+mid],g);
					a[sta+j]=add(x,y);a[sta+j+mid]=dec(x,y);
				}
			}
		}
		if (typ==-1)
		{
			int inv=getinv(lim);
			rep(i,0,lim-1) a[i]=mul(a[i],inv);
		}
	}
	
	void Inv(int *a,int *b,int n)
	{
		if (n==1)
		{
			b[0]=getinv(a[0]);
			return;
		}
		Inv(a,b,(n+1)>>1);
		int lim=calcr(n<<1);
		rep(i,0,n-1) A[i]=a[i];
		ntt(lim,A,1);ntt(lim,b,1);
		rep(i,0,lim-1) b[i]=mul(b[i],dec(2,mul(A[i],b[i])));
		ntt(lim,b,-1);
		rep(i,0,lim-1) A[i]=0;
		rep(i,n,lim-1) b[i]=0;
	}
	
	void Direv(int *a,int *b,int n)
	{
		rep(i,1,n-1) b[i-1]=mul(a[i],i);
		b[n-1]=0;
	}
	
	void Inter(int *a,int *b,int n)
	{
		rep(i,1,n-1) b[i]=mul(a[i-1],getinv(i));
		b[0]=0;
	}
	
	void Ln(int *a,int *b,int n)
	{
		rep(i,0,(n<<1)-1) B[i]=C[i]=0;
		Direv(a,B,n);Inv(a,C,n);
		int lim=calcr(n<<1);
		ntt(lim,B,1);ntt(lim,C,1);
		rep(i,0,lim-1) B[i]=mul(B[i],C[i]);
		ntt(lim,B,-1);Inter(B,b,n);
		rep(i,n,lim-1) b[i]=0;
	}
	
	void Exp(int *a,int *b,int n)
	{
		if (n==1) {b[0]=1;return;}
		Exp(a,b,(n+1)>>1);
		Ln(b,D,n);
		rep(i,0,n-1)
		{
			if (i) D[i]=add(maxd-D[i],a[i]);
			else D[i]=add(maxd-D[i],a[i]+1);
		}
		int lim=calcr(n<<1);
		ntt(lim,D,1);ntt(lim,b,1);
		rep(i,0,lim-1) b[i]=mul(b[i],D[i]);
		ntt(lim,b,-1);
		rep(i,n,lim-1) b[i]=0;
	}
	
	void Sqrt(int *a,int *b,int n)
	{
		if (n==1) {b[0]=1;return;}
		Sqrt(a,b,(n+1)>>1);
		rep(i,0,(n<<1)-1) B[i]=C[i]=0;
		Inv(b,B,n);
		int lim=calcr(n<<1);
		rep(i,0,n-1) C[i]=a[i];
		ntt(lim,B,1);ntt(lim,C,1);
		rep(i,0,lim-1) B[i]=mul(C[i],B[i]);
		ntt(lim,B,-1);
		rep(i,0,n-1) b[i]=mul(add(b[i],B[i]),inv2);
	}
}
using namespace Polynomial;
posted @ 2020-04-08 23:52  EncodeTalker  阅读(159)  评论(0编辑  收藏  举报