• 多项式全家桶


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

    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;
    
  • 相关阅读:
    JAVA_WEB--jsp概述
    npr_news英语新闻听力——每日更新
    词根词缀高效背单词技巧--词霸天下完整版
    python刷LeetCode:1071. 字符串的最大公因子
    python刷LeetCode:1013. 将数组分成和相等的三个部分
    python刷LeetCode:543. 二叉树的直径
    python刷LeetCode:121. 买卖股票的最佳时机
    python刷LeetCode:38. 外观数列
    python刷LeetCode:35. 搜索插入位置
    python刷LeetCode:28. 实现 strStr()
  • 原文地址:https://www.cnblogs.com/encodetalker/p/12663969.html
Copyright © 2020-2023  润新知