多项式全家桶

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

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} ceil}))

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

注意到这个式子可以平方一下(G(x)^2-2A(x)G(x)+A(x)^2equiv0(mod x^n))
两边乘上(F(x),G(x)-2A(x)+F(x)A(x)^2equiv0(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} ceil})).

考虑将(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)的意义下,若(kgeq 2), 那么((G(x)-G_1(x))^kequiv 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)^2equiv 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;
原文地址:https://www.cnblogs.com/encodetalker/p/12663969.html