zoukankan      html  css  js  c++  java
  • 多项式进阶操作

    看到成爷爷写了一篇就跟风写一篇吧,慢慢更

    一些前置技能

    泰勒展开

    ddy:这里用大家小学二年级学的泰勒展开。。。
    

    对于一个不是很好表示的函数(f(x)),我们可以用一个多项式函数来拟合它,用神奇的泰勒展开即可

    [f(x)=f(x_0)+frac{f'(x_0)}{1!}(x-x_0)+frac{f''(x_0)}{2!}(x-x_0)^2+...+frac{f^n(x_0)}{n!}(x-x_0)^n+... ]

    这里的(n)越大,误差就越小;由于常规的多项式操作都是在({ m mod} x^n)意义下,所以我们只需要展开前(n)项即可

    多项式牛顿迭代

    定义一个多项式(F(x)),求一个多项式(G(x))满足(F(G(x))equiv 0({ m mod} x^n))

    考虑现在有(F(A(x))equiv 0({ m mod} x^n)),推出(F(B(x))equiv 0({ m mod} x^{2n}))

    不难发现有(B-Aequiv 0({ m mod} x^n)),则((B-A)^2equiv 0({ m mod} x^n))

    我们在(A)处对(F(B))做泰勒展开

    [F(B)equiv F(A)+F'(A)(B-A)+F''(A)(B- A)^2+...({ m mod} x^{2n})]

    由于((B-A)^2equiv 0({ m mod} x^n)),所以后面的那些东西相当于没有,只剩下了

    [F(B)equiv F(A)+F'(A)(B-A) ({ m mod} x^{2n}) ]

    整理一下即

    [Bequiv A+frac{F(B)-F(A)}{F'(A)} ({ m mod} x^{2n}) ]

    注意一下(F(B)equiv 0({ m mod} x^{2n})),于是

    [Bequiv A-frac{F(A)}{F'(A)}({ m mod} x^{2n}) ]

    1.多项式求逆

    就是给定多项式(F(x)),求一个多项式(G(x))

    使得

    [F(x) imes G(x)equiv 1(mod x^n) ]

    (mod x^n)就是只考虑前(n)

    这是一个基于倍增的算法,就是推一下如何从(mod x^{frac{n}{2}})推到(mod x^n)

    设多项式(B'(x))满足

    [F(x)B'(x)equiv 1(mod x^{frac{n}{2}}) ]

    求多项式(B(x))满足

    [F(x)B(x)equiv 1(mod x^{n}) ]

    后面那个忽略前(n)项,所以前(frac{n}{2})项肯定会满足

    [F(x)B(x)equiv 1(mod x^{frac{n}{2}}) ]

    和第一个柿子一减

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

    显然(F)不可能是(0),于是只能是((B-B')=0)

    于是

    [B-B'equiv 0 (mod x^{frac{n}{2}}) ]

    两边平方并且乘上(F)

    [F(B-B')^2equiv 0(mod x^n) ]

    [F imes B^2-F imes 2BB'+F imes B'^2equiv 0(mod x^n) ]

    别忘了(F(x)B(x)equiv 1(mod x^{frac{n}{2}}))就有

    [B-2B'+FB'^2equiv 0(mod x^n) ]

    [Bequiv 2B'-FB'^2(mod x^n) ]

    发现这个东西可以倍增来求,就是需要一个(NTT)

    时间复杂度(T(n)=T(n/2)+O(nlogn)=O(nlogn))

    代码

    #include<algorithm>
    #include<iostream>
    #include<cstring>
    #include<cstdio>
    #define maxn 300005
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const LL mod=998244353;
    const LL g[2]={332748118,3};
    int n,rev[maxn],len;
    LL a[maxn],b[maxn],C[maxn];
    inline LL quick(LL a,LL b) {LL S=1;while(b) {if(b&1ll) S=S*a%mod;b>>=1ll;a=a*a%mod;}return S;}
    inline void NTT(LL *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1;LL og1=quick(g[o],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			LL og=1,t;
    			for(re int x=l;x<l+ln;x++) {
    				t=(og*f[ln+x])%mod;
    				f[ln+x]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og=(og*og1)%mod;
    			}
    		}
    	}
    	if(o) return;
    	LL inv=quick(len,mod-2);
    	for(re int i=0;i<len;i++) f[i]=(f[i]*inv)%mod;
    }
    inline void Inv(int n,LL *A,LL *B) {
    	if(n==1) {B[0]=quick(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);len=1;
    	while(len<=n+n-2) len<<=1;
    	for(re int i=0;i<n;i++) C[i]=A[i];
    	for(re int i=n;i<len;i++) C[i]=0;
    	for(re int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
    	NTT(C,1),NTT(B,1);
    	for(re int i=0;i<len;i++) B[i]=(2ll*B[i]%mod-C[i]*B[i]%mod*B[i]+mod)%mod,B[i]=(B[i]+mod)%mod;
    	NTT(B,0);for(re int i=n;i<len;i++) B[i]=0;
    	
    } 
    int main()
    {
    	n=read();
    	for(re int i=0;i<n;i++) a[i]=read();
    	Inv(n,a,b);
    	for(re int i=0;i<n;i++) printf("%lld ",b[i]);
    	return 0;
    }
    

    2.分治(NTT)

    拿luogu上的板子吧

    给定(G={g_1,g_2,g_3...})

    [f_i=sum_{j=1}^if_{j}g_{i-j} ]

    原先没学生成函数,自然做不动了,现在用生成函数来求就很休闲了

    (f)的生成函数为(F(x))

    就有

    [F(x)=sum_{i=0}(sum_{j=1}^if_jg_{i-j}+[i=0])x^i ]

    [=1+sum_{j=1}^if_jx^j imes g_{i-j}x^{i-j} ]

    发现里面是(F(x) imes G(x))

    于是就有

    [F(x)=1+F(x)G(x) ]

    解得

    [F(x)=frac{1}{1-G(x)} ]

    这是生成函数的关系,从多项式来看,就是对(1-G(x))求一个逆

    于是还是套多项式求逆的板子

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    const int maxn=262144+10;
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const LL G[2]={3,332748118};
    const LL mod=998244353;
    LL g[maxn],b[maxn],C[maxn];
    int rev[maxn];
    int n,len;
    inline LL quick(LL a,LL b) {LL S=1;while(b){if(b&1) S=S*a%mod;b>>=1;a=a*a%mod;}return S;}
    inline void NTT(LL *f,int v) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1;LL og1=quick(G[v],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			LL t,og=1;
    			for(re int x=l;x<l+ln;x++) {
    				t=(og*f[ln+x])%mod;
    				f[x+ln]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og*=og1;og%=mod;
    			}
    		}
    	}
    	if(!v) return;
    	LL inv=quick(len,mod-2);
    	for(re int i=0;i<len;i++) f[i]=(f[i]*inv)%mod;
    }
    inline void Inv(int n,LL *A,LL *B) {
        if(n==1) {B[0]=quick(A[0],mod-2);return;}
        Inv((n+1)>>1,A,B);len=1;
        while(len<=n+n-2) len<<=1;
        for(re int i=0;i<n;i++) C[i]=A[i];
        for(re int i=n;i<len;i++) C[i]=0;
        for(re int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?len>>1:0);
        NTT(C,0),NTT(B,0);
        for(re int i=0;i<len;i++) B[i]=(2ll*B[i]%mod-C[i]*B[i]%mod*B[i]+mod)%mod,B[i]=(B[i]+mod)%mod;
        NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    } 
    int main() {
    	n=read();g[0]=1;
    	for(re int i=1;i<n;i++) g[i]=mod-read();
    	Inv(n,g,b);
    	for(re int i=0;i<n;i++) printf("%d ",(int)b[i]);
    	return 0;
    }
    

    3.多项式开根

    就是给你一个多项式(F(x)),让你找到一个(B(x)),满足(B^2(x)equiv F(x) [n])

    和求逆一样,也是考虑倍增假设我们现在求得了(B'^2(x)equiv F [frac{n}{2}]),考虑推出(B^2(x)equiv F [n])

    (n)项相同肯定也满足前(frac{n}{2})项相同

    于是

    [B^2(x)equiv F(x) [frac{n}{2}] ]

    (B'^2(x)equiv F(x) [frac{n}{2}])做一个差

    [B^2(x)-B'^2(x)equiv 0 [frac{n}{2}] ]

    [B^4(x)-2B^2(x)B'^2(x)+B'^4(x)equiv 0 [n] ]

    左右两边加上(4B^2(x)B'^2(x))

    [B^4(x)+2B^2(x)B'^2(x)+B'^4(x)equiv 4B^2(x)B'^2(x) [n] ]

    [(B^2(x)+B'^2(x))^2equiv (2B(x)B'(x))^2 [n] ]

    [B^2(x)+B'^2(x)equiv 2B(x)B'(x) [n] ]

    [B(x)equivfrac{F(x)+B'^2(x)}{2B'(x)} [n] ]

    发现开根还需要套上一个求逆,但是复杂度依旧是(O(nlogn))

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    const int maxn=262144+1005;
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const LL mod=998244353;
    const LL G[2]={3,332748118};
    int rev[maxn],len,n;
    LL A[maxn],B[maxn],C[maxn],D[maxn],K[maxn],H[maxn];
    inline LL ksm(LL a,LL b) {LL S=1;while(b) {if(b&1) S=S*a%mod;b>>=1;a=a*a%mod;}return S;}
    inline void NTT(LL *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1;LL og1=ksm(G[o],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			LL t,og=1;
    			for(re int x=l;x<l+ln;x++) {
    				t=(f[ln+x]*og)%mod;
    				f[x+ln]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og=(og*og1)%mod;
    			}
    		}
    	}
    	if(!o) return;
    	LL inv=ksm(len,mod-2);
    	for(re int i=0;i<len;i++) f[i]=(f[i]*inv)%mod;
    }
    inline void Inv(int n,LL *A,LL *B) {
    	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);len=1;
    	while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	for(re int i=0;i<n;i++) C[i]=A[i];
    	for(re int i=n;i<len;i++) C[i]=0;
    	NTT(C,0),NTT(B,0);
    	for(re int i=0;i<len;i++) B[i]=(2ll*B[i]-C[i]*B[i]%mod*B[i]%mod+mod)%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    inline void mul(LL *A,LL *B,int n) {
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	NTT(A,0),NTT(B,0);
    	for(re int i=0;i<len;i++) A[i]=(A[i]*B[i])%mod;
    	NTT(A,1);for(re int i=n;i<len;i++) A[i]=0;
    }
    inline void Sqrt(int n,LL *A,LL *B) {
    	if(n==1) {B[0]=1;return;}
    	Sqrt((n+1)>>1,A,B);
    	memset(K,0,sizeof(K));
    	memset(D,0,sizeof(D));
    	memset(H,0,sizeof(H));
    	for(re int i=0;i<n;i++) K[i]=(B[i]*2ll)%mod;Inv(n,K,D);
    	for(re int i=0;i<n;i++) H[i]=B[i];mul(B,H,n);
    	for(re int i=0;i<n;i++) B[i]=(B[i]+A[i])%mod;mul(B,D,n);
    }
    int main() {
    	n=read();
    	for(re int i=0;i<n;i++) A[i]=read();
    	Sqrt(n,A,B);
    	for(re int i=0;i<n;i++) printf("%lld ",B[i]);
    	return 0;
    }
    

    4.多项式求导和积分

    不会,于是就背过吧

    如果有

    [F(x)=sum_{i=0}^nf_ix^i ]

    那么我们对(F)求导就有

    [F'(x)=sum_{i=1}^nif_ix^{i-1} ]

    显然,求导和积分是逆操作,于是我们对(F)求积分就有

    [int Fdx=sum_{i=1}^{n+1}frac{f_{i-1}}{i}x^i ]

    5.多项式对数函数

    板子

    我们求一个多少项式(G(x)),满足(G(x)=ln F(x))

    我们先两边求导

    [G(x)'=ln'F(x) ]

    对于(x)来说满足(ln'x=frac{1}{x}),于是

    [G'(x)=frac{F'(x)}{F(x)} ]

    我们这样求出来是(G'(x)),之后求个积分就好了

    于是

    [ln F(x)=int frac{F'(x)}{F(x)} ]

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    const int maxn=262144+5;
    const int mod=998244353;
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const int G[2]={3,332748118};
    int n,len,rev[maxn],F[maxn],K[maxn],T[maxn],g[maxn],H[maxn],inv[maxn];
    inline int ksm(int a,int b) {
    	int S=1;
    	while(b) {if(b&1) S=(1ll*S*a)%mod;b>>=1;a=(1ll*a*a)%mod;}
    	return S;
    }
    inline void NTT(int *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1,og1=ksm(G[o],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			int t,og=1;
    			for(re int x=l;x<l+ln;++x) {
    				t=(1ll*f[x+ln]*og)%mod;
    				f[x+ln]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og=(1ll*og*og1)%mod;
    			}
    		}
    	}
    	if(!o) return;
    	int Inv=ksm(len,mod-2);
    	for(re int i=0;i<len;i++) f[i]=(1ll*f[i]*Inv)%mod;
    }
    void Inv(int n,int *A,int *B) {
    	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	for(re int i=0;i<n;i++) K[i]=B[i],T[i]=A[i];
    	for(re int i=n;i<len;i++) K[i]=T[i]=0;
    	NTT(K,0),NTT(T,0);
    	for(re int i=0;i<len;i++) 
    		B[i]=(2ll*K[i]-1ll*T[i]*K[i]%mod*K[i]%mod+mod)%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    int main() {
    	n=read();
    	for(re int i=0;i<n;i++) F[i]=read();
    	for(re int i=1;i<n;i++) g[i-1]=(1ll*i*F[i])%mod;
    	Inv(n,F,H);len=1;while(len<n+n) len<<=1;
    	NTT(H,0),NTT(g,0);
    	for(re int i=0;i<len;i++) H[i]=(1ll*H[i]*g[i])%mod;
    	NTT(H,1);putchar('0');putchar(' ');
    	inv[1]=1;
    	for(re int i=2;i<n;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    	for(re int i=1;i<n;i++) printf("%d ",1ll*H[i-1]*inv[i]%mod);
    	return 0;
    }
    

    6.多项式指数函数

    板子

    首先前置的东西是牛顿迭代

    是这样的,如果我们有一个多项式(F(x)),需要求一个多项式(G(x))满足

    [F(G(x))equiv 0 [n] ]

    假设我们已经求出了(F(G_0(x))equiv 0 [frac{n}{2}])

    那么

    [G=G_0-frac{F(G_0)}{F'(G_0)} ]

    现在我们要求的东西是(G(x)=e^{F(x)})

    两边取一个(ln)

    [ln G=F ]

    [ln G-F=0 ]

    (A(G)=ln G-F),把(F)看成一个常数项,显然有(A'(G)=frac{1}{G})

    假设我们已经求出了

    [A(G_0)equiv 0 [frac{n}{2}] ]

    需要求

    [A(G)equiv 0 [n] ]

    直接套上牛顿迭代

    [G=G_0-frac{A(G_0)}{A'(G_0)}=G_0-G_0(ln G_0-F)=G_0(1-ln G_0+F) ]

    于是我们倍增就好啦

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const int maxn=262144+5;
    const int mod=998244353;
    const int G[2]={3,332748118};
    int n,inv[maxn],a[maxn],b[maxn],K[maxn];
    int H[maxn],T[maxn],g[maxn],C[maxn],rev[maxn],len;
    inline int ksm(int a,int b) {
    	int S=1;
    	while(b) {if(b&1) S=(1ll*S*a)%mod;b>>=1;a=(1ll*a*a)%mod;}
    	return S;
    }
    inline void NTT(int *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1,og1=ksm(G[o],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			int t,og=1;
    			for(re int x=l;x<l+ln;++x) {
    				t=(1ll*f[ln+x]*og)%mod;
    				f[x+ln]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og=(1ll*og*og1)%mod;
    			}
    		}
    	}
    	if(!o) return;
    	int Inv=inv[len];
    	for(re int i=0;i<len;i++) f[i]=(1ll*f[i]*Inv)%mod;
    }
    void Inv(int n,int *A,int *B) {
    	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	for(re int i=0;i<n;i++) K[i]=A[i];
    	for(re int i=n;i<len;i++) K[i]=0;
    	NTT(K,0),NTT(B,0);
    	for(re int i=0;i<len;i++)
    		B[i]=(2ll*B[i]-1ll*K[i]*B[i]%mod*B[i]%mod+mod)%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    void Ln(int n,int *A,int *B) {
    	for(re int i=0;i<n;i++) g[i]=0,B[i]=0;
    	for(re int i=1;i<n;i++) g[i-1]=(1ll*i*A[i])%mod;
    	memset(C,0,sizeof(C));Inv(n,A,C);
    	len=1;while(len<n+n) len<<=1;
    	NTT(C,0),NTT(g,0); 
    	for(re int i=0;i<len;i++) g[i]=(1ll*g[i]*C[i])%mod;
    	NTT(g,1);
    	for(re int i=1;i<n;i++) B[i]=(1ll*inv[i]*g[i-1])%mod;
    }
    void Exp(int n,int *A,int *B) {
    	if(n==1) {B[0]=1;return;}
    	Exp((n+1)>>1,A,B);
    	Ln(n,B,T);len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	for(re int i=0;i<n;i++) T[i]=(A[i]-T[i]+mod)%mod;
    	for(re int i=n;i<len;i++) T[i]=0;T[0]++;
    	NTT(B,0),NTT(T,0);
    	for(re int i=0;i<len;i++) B[i]=(1ll*B[i]*T[i])%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    int main() {
    	n=read();
    	for(re int i=0;i<n;i++) a[i]=read();
    	inv[1]=1;
    	for(re int i=2;i<=262144;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod; 
    	Exp(n,a,b);
    	for(re int i=0;i<n;i++) printf("%d ",b[i]);
    	return 0;
    }
    

    7.多项式快速幂

    板子

    显然直接裸上快速幂是两个(log)

    我们考虑要求的东西是(A^k(x))

    直接取一个(ln)就能把(k)从指数上拿下来啦,之后我们得到的是(ln A^k),再(exp)回去就好了

    于是就是

    [A^k(x)=exp(kln A(x)) ]

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    const int maxn=262144+5;
    const int mod=998244353;
    const int G[2]={3,332748118};
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    int n,len,rev[maxn],inv[maxn],T[maxn],C[maxn];
    int a[maxn],K[maxn],g[maxn],H[maxn],O[maxn],m;
    inline int getPow() {
    	int x=0;char c=getchar();
    	while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(10ll*x+c-48)%mod,c=getchar();
    	return x;
    }
    inline int ksm(int a,int b) {
    	int S=1;
    	while(b) {if(b&1) S=(1ll*S*a)%mod;b>>=1;a=(1ll*a*a)%mod;}
    	return S;
    }
    inline void NTT(int *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1,og1=ksm(G[o],(mod-1)/i);
    		for(re int l=0;l<len;l+=i) {
    			int t,og=1;
    			for(re int x=l;x<l+ln;++x) {
    				t=(1ll*f[ln+x]*og)%mod;
    				f[ln+x]=(f[x]-t+mod)%mod;
    				f[x]=(f[x]+t)%mod;
    				og=(1ll*og*og1)%mod;
    			}
    		}
    	}
    	if(!o) return;
    	int Inv=inv[len];
    	for(re int i=0;i<len;i++) f[i]=(1ll*f[i]*Inv)%mod;
    }
    void Inv(int n,int *A,int *B) {
    	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);
    	for(re int i=0;i<n;i++) T[i]=A[i];
    	for(re int i=n;i<len;i++) T[i]=0;
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	NTT(T,0),NTT(B,0);
    	for(re int i=0;i<len;i++)
    		B[i]=(2ll*B[i]-1ll*T[i]*B[i]%mod*B[i]%mod+mod)%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    void Ln(int n,int *A,int *B) {
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) g[i]=B[i]=0;
    	for(re int i=1;i<n;i++) g[i-1]=(1ll*i*A[i])%mod;
    	memset(C,0,sizeof(C));Inv(n,A,C);
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	NTT(C,0),NTT(g,0);
    	for(re int i=0;i<len;i++) C[i]=(1ll*g[i]*C[i])%mod;
    	NTT(C,1);for(re int i=1;i<n;i++) B[i]=(1ll*inv[i]*C[i-1])%mod;
    }
    void Exp(int n,int *A,int *B) {
    	if(n==1) {B[0]=1;return;}
    	Exp((n+1)>>1,A,B);Ln(n,B,K);
    	len=1;while(len<n+n) len<<=1;
    	for(re int i=0;i<n;i++) K[i]=(A[i]-K[i]+mod)%mod;
    	for(re int i=n;i<len;i++) K[i]=0;K[0]++;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	NTT(B,0),NTT(K,0);
    	for(re int i=0;i<len;i++) B[i]=(1ll*B[i]*K[i])%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    int main() {
    	inv[1]=1;
    	for(re int i=2;i<=262144;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    	n=read(),m=getPow();
    	for(re int i=0;i<n;i++) a[i]=read();
    	Ln(n,a,H);
    	for(re int i=0;i<n;i++) H[i]=(1ll*H[i]*m)%mod;
    	Exp(n,H,O);
    	for(re int i=0;i<n;i++) printf("%d ",O[i]);
    	return 0;
    }
    

    8.多项式带余除法

    板子

    就是给定两个多项式(F(x),G(x)),求满足如下条件的两个多项式

    [F(x)=G(x) imes Q(x)+R(x) ]

    其中(F(x))(n)次多项式,(G(x))(m)次多项式,要求(Q(x))(n-m)次多项式,(R(x))最高次项小于(m)

    有一个非常牛逼的东西,定义(F_r(x))为多项式(F(x))的转置

    那么就有

    [F_r(x)=x^nF(frac{1}{x})(mod x^n) ]

    这非常显然

    于是我们来化柿子了

    [F(x)=G(x) imes Q(x)+R(x) ]

    [F(frac{1}{x})=G(frac{1}{x}) imes Q(frac{1}{x})+R(frac{1}{x}) ]

    [x^nF(frac{1}{x})=x^mG(frac{1}{x}) imes x^{n-m}Q(frac{1}{x})+x^{n-m+1} imes x^{m-1}R(frac{1}{x}) ]

    [F_r(x)=G_r(x) imes Q_r(x)+x^{n-m+1} imes R_r(x) ]

    [F_r(x)equiv G_r(x) imes Q_r(x)(mod x^{n-m+1}) ]

    [Q_r(x)equiv F_r(x) imes G_r^{-1}(x)(mod x^{n-m+1}) ]

    于是多项式求逆即可

    求出(G_r(x))显然(R(x))就很好求了

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    #include<algorithm>
    #define re register
    #define LL long long
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    inline int read() {
    	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
    	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
    }
    const int mod=998244353;
    const int G[2]={3,(mod+1)/3};
    const int maxn=262144+5;
    int n,m,len;
    int C[maxn],a[maxn],b[maxn],res[maxn],rev[maxn],ar[maxn],br[maxn],T[maxn],K[maxn];
    inline int ksm(int a,int b) {
    	int S=1;for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) S=1ll*S*a%mod;return S;
    }
    inline void NTT(int *f,int o) {
    	for(re int i=0;i<len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
    	for(re int i=2;i<=len;i<<=1) {
    		int ln=i>>1,og1=ksm(G[o],(mod-1)/i);
    		for(re int t,og=1,l=0;l<len;l+=i,og=1)
    			for(re int x=l;x<l+ln;++x) {
    				t=1ll*og*f[x+ln]%mod,og=1ll*og*og1%mod;
    				f[x+ln]=(f[x]-t+mod)%mod,f[x]=(f[x]+t)%mod;
    			}
    	}
    	if(!o) return;
    	int Inv=ksm(len,mod-2);
    	for(re int i=0;i<len;i++) f[i]=1ll*f[i]*Inv%mod;
    }
    inline void Rev(int n,int *A,int *B) {
    	for(re int i=0;i<n;i++) A[i]=B[n-i-1];
    }
    inline void Inv(int n,int *A,int *B) {
    	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
    	Inv((n+1)>>1,A,B);
    	len=1;while(len<n+n-1) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	for(re int i=0;i<n;i++) C[i]=A[i];
    	for(re int i=n;i<len;i++) C[i]=0;
    	NTT(C,0),NTT(B,0);
    	for(re int i=0;i<len;i++) B[i]=(2ll*B[i]-1ll*C[i]*B[i]%mod*B[i]%mod+mod)%mod;
    	NTT(B,1);for(re int i=n;i<len;i++) B[i]=0;
    }
    int main() {
    	n=read();m=read();
    	for(re int i=0;i<=n;i++) a[i]=read();
    	for(re int i=0;i<=m;i++) b[i]=read();
    	Rev(m+1,br,b),Rev(n+1,ar,a);
    	Inv(n-m+1,br,T);
    	len=1;while(len<n+1+n-m+1) len<<=1;
    	for(re int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    	NTT(T,0),NTT(ar,0);
    	for(re int i=0;i<len;i++) ar[i]=1ll*ar[i]*T[i]%mod;
    	NTT(ar,1);for(re int i=n-m;i>=0;i--) printf("%d ",ar[i]);
    	for(re int i=n-m+1;i<len;i++) ar[i]=0;
    	Rev(n-m+1,K,ar);puts("");
    	NTT(K,0),NTT(b,0);
    	for(re int i=0;i<len;i++) b[i]=1ll*b[i]*K[i]%mod;
    	NTT(b,1);
    	for(re int i=0;i<m;i++) printf("%d ",(a[i]-b[i]+mod)%mod);
    	return 0;
    }
    
  • 相关阅读:
    WHU-1551-Pairs(莫队算法+分块实现)
    JS日历控件 灵活设置: 精确的时分秒.
    java集群优化——多线程下的单例模式
    课程设计——银行系统
    互联网金融,巨头天下还是创业者天堂?
    Android 使用图片异步载入框架Universal Image Loader的问题
    程序C++ to C#交互
    浅谈asp.net通过本机cookie仿百度(google)实现搜索input框自己主动弹出搜索提示
    全栈JavaScript之路(十六)HTML5 HTMLDocument 类型的变化
    推荐安卓开发神器(里面有各种UI特效和实例)
  • 原文地址:https://www.cnblogs.com/asuldb/p/10543247.html
Copyright © 2011-2022 走看看