zoukankan      html  css  js  c++  java
  • 高级(并不)多项式算法总结

    牛顿迭代

    说白了就是给你一个(F(x)),你需要求出一个(G(x)),使得(F(G(x)) equiv 0 mod x^n)

    假设我们已经求出了(H(x))满足(F(H(x)) equiv 0 mod x^n),我们需要推出(F(G(x)) equiv 0 mod x^{2n})。我们可以在(H(x))处使用正无穷次拉格朗日中值定理,可以得到泰勒展开式:

    [F(G(x)) equiv F(H(x))+frac{F'(H(x))(G(x)-H(x))}{1!}+frac{F''(H(x))(G(x)-H(x))^2}{2!}+... mod x^{2n} ]

    由于对于任意的((G(x)-H(x))^q,q geq 2)都满足((G(x)-H(x))^q equiv 0 mod x^{2n})

    所以有:

    [F(G(x)) equiv F(H(x))+F'(H(x))(G(x)-H(x)) mod x^{2n} ]

    因为需要(F(G(x)) equiv 0 mod x^{2n})

    所以:

    [F(H(x))+F'(H(x))(G(x)-H(x)) equiv 0 mod x^{2n} ]

    简单化一下就是:

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

    这个东西可以用来计算多项式开方,多项式(exp)

    多项式(ln)

    给你一个(A(x)),你需要求出一个(B(x) equiv ln A(x) mod x^n)

    两边同时求导,根据复合函数求导的链式法则,有:

    [B'(x) equiv ln'(A(x))A'(x) equiv frac{A'(x)}{A(x)} mod x^n ]

    注意这里(A(x))的常数项必须为(1),求出来的(B(x))常数项为(0)

    多项式(exp)

    给你一个(A(x)),你需要求出一个(B(x) equiv e^{A(x)} mod x^n)

    题中的式子可以转化为(ln B(x) equiv A(x) mod x^n)

    现在我们要求的东西相当于函数(F(B(x)) = A(x)-ln B(x))(x^n)意义下的一个零点,代入牛顿迭代的式子里,有:

    [B(x) equiv C(x)-frac{A(x)-ln C(x)}{-frac{1}{C(x)}} mod x^{2n} ]

    [B(x) equiv C(x)(1+A(x)-ln C(x)) mod x^{2n} ]

    注意这里(A(x))的常数项必须为(0),求出来的(B(x))常数项为(1)

    代码

    应该只放一个多项式(exp)的代码就够了吧。

    #include <bits/stdc++.h>
    
    #define rin(i,a,b) for(register int i=(a);i<=(b);++i)
    #define irin(i,a,b) for(register int i=(a);i>=(b);--i)
    #define trav(i,a) for(register int i=head[a];i;i=e[i].nxt)
    #define Size(a) (int)a.size()
    #define pb push_back
    typedef long long LL;
    
    using std::cerr;
    using std::endl;
    
    inline int read(){
    	int x=0,f=1;char ch=getchar();
    	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
    	return x*f;
    }
    
    const int MAXN=100005;
    const int MOD=998244353;
    const int G=3,INVG=332748118;
    
    int NTT,N,n,m,len,rev[MAXN<<2];
    int w[1048576],iw[1048576],inv[MAXN];
    int a[MAXN];
    int A[MAXN<<2],B[MAXN<<2],C[MAXN<<2];
    
    inline int qpow(int x,int y){
    	int ret=1,tt=x%MOD;
    	while(y){
    		if(y&1) ret=1ll*ret*tt%MOD;
    		tt=1ll*tt*tt%MOD;
    		y>>=1;
    	}
    	return ret;
    }
    
    void ntt(int *c,int dft){
    	rin(i,0,n-1) if(i<rev[i])
    		std::swap(c[i],c[rev[i]]);
    	for(register int mid=1;mid<n;mid<<=1){
    		int r=(mid<<1),u=NTT/r;
    		for(register int l=0;l<n;l+=r){
    			int v=0;
    			for(register int i=0;i<mid;++i,v+=u){
    				int x=c[l+i],y=1ll*c[l+mid+i]*(dft>0?w[v]:iw[v])%MOD;
    				c[l+i]=x+y<MOD?x+y:x+y-MOD;
    				c[l+mid+i]=x-y>=0?x-y:x-y+MOD;
    			}
    		}
    	}
    	if(dft<0){
    		int invn=qpow(n,MOD-2);
    		rin(i,0,n-1) c[i]=1ll*c[i]*invn%MOD;
    	}
    }
    
    void prepare(){
    	for(n=1,len=0;n<=m;n<<=1,++len);
    	rin(i,1,n-1) rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
    }
    
    // 给B求导到A
    void deriv(int idx){rin(i,0,idx-2)A[i]=1ll*B[i+1]*(i+1)%MOD;}
    // 给A积分到A
    void integ(int idx){irin(i,idx-1,1)A[i]=1ll*A[i-1]*inv[i]%MOD;}
    
    // 给B求逆到C利用A
    void inver(int idx){
    	if(idx==1){C[0]=qpow(B[0],MOD-2);return;}
    	inver((idx+1)/2);m=(idx-1)+((idx+1)/2-1)*2;prepare();
    	rin(i,0,idx-1)A[i]=B[i];
    	ntt(A,1);ntt(C,1);rin(i,0,n-1)C[i]=(2*C[i]-1ll*A[i]*C[i]%MOD*C[i]%MOD+MOD)%MOD;ntt(C,-1);
    	memset(A,0,n<<2);rin(i,idx,n-1)C[i]=0;
    }
    
    // 给B求ln到A利用C
    void polyln(int idx){
    	inver(idx);deriv(idx);
    	m=(idx-2)+(idx-1);prepare();
    	ntt(A,1);ntt(C,1);rin(i,0,n-1)A[i]=1ll*A[i]*C[i]%MOD;ntt(A,-1);
    	integ(idx);A[0]=0;rin(i,idx,n-1)A[i]=0;memset(C,0,n<<2);
    }
    
    // 给a求exp到B利用A
    void polyexp(int idx){
    	if(idx==1){B[0]=1;return;}
    	polyexp((idx+1)/2);polyln(idx);m=((idx+1)/2-1)+(idx-1);prepare();
    	rin(i,0,idx-1)A[i]=(a[i]-A[i]+MOD)%MOD;A[0]=(A[0]+1)%MOD;
    	ntt(A,1);ntt(B,1);rin(i,0,n-1)B[i]=1ll*A[i]*B[i]%MOD;ntt(B,-1);
    	memset(A,0,n<<2);rin(i,idx,n-1)B[i]=0;
    }
    
    void init(int N){
    	for(NTT=1;NTT<(N<<1);NTT<<=1);
    	int v=qpow(G,(MOD-1)/NTT),iv=qpow(INVG,(MOD-1)/NTT);w[0]=iw[0]=1;
    	rin(i,1,NTT-1) w[i]=1ll*w[i-1]*v%MOD,iw[i]=1ll*iw[i-1]*iv%MOD;
    	inv[1]=1;rin(i,2,N) inv[i]=(-1ll*(MOD/i)*inv[MOD%i]%MOD+MOD)%MOD;
    }
    
    int main(){
    	N=read();init(N);
    	rin(i,0,N-1) a[i]=read();
    	polyexp(N);
    	rin(i,0,N-1) printf("%d ",B[i]);putchar('
    ');
    	return 0;
    }
    
  • 相关阅读:
    hdoj 2041 超级楼梯
    hdoj 刚入门~把11页A了一些~~
    编写Powerdesigner脚本,快速生成数据库表
    网站
    duwamish
    面试题
    http://www.sqlclub.com
    Microsoft PetShop 3.0 设计与实现数据访问层
    Java代码查询站点
    从Blog上面去掉那该死的“狗狗订阅”的logo !
  • 原文地址:https://www.cnblogs.com/ErkkiErkko/p/10582531.html
Copyright © 2011-2022 走看看