• 【模板】多项式指数函数(多项式 exp)


    IV.【模板】多项式指数函数(多项式 exp)

    本题有两种解法,一种比较好理解,一种比较通用(并且速度快)。

    首先法一便是分治FFT解法。

    我们有

    \[B=e^A \]

    于是两边求导,得到

    \[B'=A'e^A \]

    因为又有\(B=e^A\),代入得

    \[B'=A'B \]

    我们再积分回去,得到

    \[\int B'=\int A'B \]

    也即

    \[B=\int A'B \]

    我们设\(A'B=T\),于是

    \[B=\int T \]

    所以

    \[B_i=\dfrac{T_{i-1}}{i} \]

    我们回到公式\(A'B=T\),将其展开得到

    \[T_i=\sum\limits_{j=0}^{i}B_jA_{i-j}' \]

    于是

    \[B_i=\dfrac{\sum\limits_{j=0}^{i-1}B_jA_{i-j-1}'}{i} \]

    为同分治NTT的形式统一,我们设一个\(C_i=A_{i-1}'=iA_i\)

    于是

    \[B_i=\dfrac{\sum\limits_{j=0}^{i-1}B_jC_{i-j}}{i} \]

    很明显有\(C_0=0\),因此我们仍然可以写成

    \[B_i=\dfrac{\sum\limits_{j=0}^iB_jC_{i-j}}{i} \]

    注意到一定有\(B_0=1\),这可以从式子中直接得出。

    然后就是分治NTT的模板题了。只需要记得当分治区间大小为\(1\)时,将其所代表的数除以\(i\)即可。

    时间复杂度\(O(n\log^2n)\)

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    const int N=1<<20;
    const int mod=998244353;
    const int G=3;
    int n,f[N],g[N],rev[N],lim,LG,invlim,all;
    int ksm(int x,int y){
    	int z=1;
    	for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
    	return z;
    }
    void NTT(int *a,int tp){
    	for(int i=0;i<lim;i++)if(rev[i]>i)swap(a[i],a[rev[i]]);
    	for(int md=1;md<lim;md<<=1){
    		int rt=ksm(G,(mod-1)/(md<<1));
    		if(tp==-1)rt=ksm(rt,mod-2);
    		for(int stp=(md<<1),pos=0;pos<lim;pos+=stp){
    			int w=1;
    			for(int i=0;i<md;i++,w=1ll*w*rt%mod){
    				int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
    				a[pos+i]=(x+y)%mod;
    				a[pos+md+i]=(x+mod-y)%mod;
    			}
    		}
    	}
    	if(tp==-1)for(int i=0;i<lim;i++)a[i]=1ll*a[i]*invlim%mod;
    }
    int a[N],b[N];
    void func(int l,int k){
    	for(int i=0;i<lim;i++)a[i]=b[i]=0;
    	for(int i=0;i<(1<<k);i++)a[i]=g[l+i];
    	for(int i=0;i<(2<<k);i++)b[i]=f[i];
    	lim=(4<<k),LG=k+2,invlim=ksm(lim,mod-2);
    	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
    	NTT(a,1),NTT(b,1);
    	for(int i=0;i<lim;i++)a[i]=1ll*a[i]*b[i]%mod;
    	NTT(a,-1);
    }
    void CDQ(int l,int k){
    	if(!k){if(l)g[l]=1ll*g[l]*ksm(l,mod-2)%mod;return;}
    	CDQ(l,k-1);
    	func(l,k-1);
    	for(int i=0;i<(1<<(k-1));i++)(g[l+(1<<(k-1))+i]+=a[(1<<(k-1))+i])%=mod;
    	CDQ(l+(1<<(k-1)),k-1);
    }
    int main(){
    	scanf("%d",&n);
    	for(int i=0;i<n;i++)scanf("%d",&f[i]),f[i]=1ll*f[i]*i%mod;
    	g[0]=1;
    	while((1<<all)<n)all++;
    	CDQ(0,all);
    	for(int i=0;i<n;i++)printf("%d ",g[i]);
    	return 0;
    }
    

    现在我们回到正解。

    首先是牛顿迭代法。这玩意是用来求方程\(G\Big(F(x)\Big)\equiv0\pmod {x^n}\)的根的,其中\(G(x)\)已知,要求出\(F(x)\)

    如同求逆一般,我们考虑递推。

    我们已经求出一个\(f(x)\)使得

    \[G\Big(f(x)\Big)\equiv0\pmod{x^n} \]

    现在我们要求出一个\(F(x)\)使得

    \[G\Big(F(x)\Big)\equiv0\pmod{x^{2n}} \]

    我们将\(G\Big(F(x)\Big)\)\(f(x)\)处泰勒展开,就有

    \[G\Big(F(x)\Big)\equiv\sum\limits_{i=0}^{\infty}\dfrac{G^{(i)}\big(f(x)\big)\Big(F(x)-f(x)\Big)^i}{i!}\pmod{x^{2n}} \]

    注意到\(f\equiv0\pmod{x^n}\Leftrightarrow f\text{前}n\text{项的系数皆为}0\),而\(F(x)-f(x)\equiv0\pmod{x^n}\)又是显然成立的,所以当它乘了\(2\)次以上的幂以后,前\(2n\)位都会为\(0\),也即

    \[\Big(F(x)-f(x)\Big)^i(i\geq 2)\equiv0\pmod{x^{2n}} \]

    所以当你把它带回泰勒展开式后,就会发现保留前\(2\)项即可。

    于是

    \[G\Big(F(x)\Big)\equiv\sum\limits_{i=0}^1\dfrac{G^{(i)}\big(f(x)\big)\Big(F(x)-f(x)\Big)^i}{i!}\pmod{x^{2n}} \]

    也即

    \[G\Big(F(x)\Big)\equiv G\big(f(x)\big)+G'\big(f(x)\big)\Big(F(x)-f(x)\Big)\pmod{x^{2n}} \]

    注意到依据定义,是有\(G\Big(F(x)\Big)\equiv0\pmod{x^{2n}}\)

    所以

    \[G\big(f(x)\big)+G'\big(f(x)\big)\Big(F(x)-f(x)\Big)\equiv 0\pmod{x^{2n}} \]

    两边同除\(G'\big(f(x)\big)\),得到

    \[\dfrac{G\big(f(x)\big)}{G'\big(f(x)\big)}+F(x)-f(x)\equiv 0\pmod{x^{2n}} \]

    所以移项,我们得到

    \[F(x)\equiv f(x)-\dfrac{G\big(f(x)\big)}{G'\big(f(x)\big)}\pmod{x^{2n}} \]

    这就是多项式牛顿迭代公式,它可以由\(f\)推至\(F\)

    我们将它框起来特别展示一下:

    \[\boxed{F(x)\equiv f(x)-\dfrac{G\big(f(x)\big)}{G'\big(f(x)\big)}\pmod{x^{2n}}} \]

    我们现在回到本题。

    我们有

    \[B(x)\equiv e^{A(x)}\pmod{x^n} \]

    两边求\(\ln\)之后,就得到了

    \[\ln B(x)\equiv A(x)\pmod{x^n} \]

    移项,有

    \[\ln B(x)-A(x)\equiv0\pmod{x^n} \]

    我们现在设\(G\Big(B(x)\Big)\equiv\ln B(x)-A(x)\pmod{x^n}\)

    按照上面的牛顿迭代公式,假设我们已经知道\(b(x)\),则有

    \[B(x)\equiv b(x)-\dfrac{G\big(b(x)\big)}{G'\big(b(x)\big)}\pmod{x^{2n}} \]

    代入\(G\)的定义,我们有(注意到\(G'\Big(b(x)\Big)=\dfrac{1}{b(x)}\)

    \[B(x)\equiv b(x)-\dfrac{\ln b(x)-A(x)}{\dfrac{1}{b_i}}\pmod{x^{2n}} \]

    变换得

    \[B(x)\equiv b(x)\Big(A(x)-\ln b(x)+1\Big)\pmod{x^{2n}} \]

    这就是多项式\(\exp\)的递推式。

    通过模板化\(\ln f\),我们可以实现它。

    代码:

    #include<bits/stdc++.h>
    using namespace std;
    const int N=1<<20;
    const int mod=998244353;
    const int G=3;
    int n,rev[N],f[N],g[N],all;
    int ksm(int x,int y){
    	int rt=1;
    	for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)rt=(1ll*rt*x)%mod;
    	return rt;
    }
    void NTT(int *a,int tp,int LG){
    	int lim=(1<<LG),invlim=ksm(lim,mod-2);
    	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
    	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    	for(int md=1;md<lim;md<<=1){
    		int rt=ksm(G,(mod-1)/(md<<1));
    		if(tp==-1)rt=ksm(rt,mod-2);
    		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
    			int w=1;
    			for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
    				int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
    				a[pos+i]=(x+y)%mod;
    				a[pos+md+i]=(x-y+mod)%mod;
    			}
    		}
    	}
    	if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
    }
    int A[N],B[N],C[N],D[N];
    void mul(int *a,int *b,int *c,int LG){
    	int lim=(1<<LG);
    	for(int i=0;i<lim;i++)A[i]=B[i]=0;
    	for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
    	NTT(A,1,LG),NTT(B,1,LG);
    	for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
    	NTT(A,-1,LG);
    	for(int i=0;i<lim;i++)c[i]=A[i];
    }
    void inv(int *a,int *b,int LG){
    	b[0]=ksm(a[0],mod-2);
    	for(int k=1;k<=LG+1;k++){
    		mul(b,a,C,k);
    		for(int i=0;i<(1<<k);i++)C[i]=(mod-C[i])%mod;
    		(C[0]+=2)%=mod;
    		mul(C,b,b,k);
    	}
    }
    void diff(int *a,int *b,int lim){
    	for(int i=0;i<lim;i++)b[i]=1ll*a[i+1]*(i+1)%mod;
    	b[lim-1]=0;
    }
    void inte(int *a,int *b,int lim){
    	for(int i=lim-1;i;i--)b[i]=1ll*a[i-1]*ksm(i,mod-2)%mod;
    	b[0]=0;
    }
    void ln(int *a,int *b,int LG){
    	inv(a,b,LG);
    	diff(a,C,1<<LG);
    	mul(b,C,b,LG+1);
    	inte(b,b,1<<LG);
    }
    void exp(int *a,int *b,int LG){
    	b[0]=1;
    	for(int k=1;k<=LG+1;k++){
    		ln(b,D,k-1);
    		for(int i=0;i<(1<<(k-1));i++)D[i]=(a[i]-D[i]+mod)%mod;
    		D[0]=(D[0]+1)%mod;
    //		for(int i=0;i<(1<<(k-1));i++)printf("%d ",D[i]);puts("");
    		mul(b,D,b,k);
    	}
    }
    int main(){
    	scanf("%d",&n);
    	for(int i=0;i<n;i++)scanf("%d",&f[i]);
    	while((1<<all)<n)all++;
    	exp(f,g,all);
    	for(int i=0;i<n;i++)printf("%d ",g[i]);puts("");
    	return 0;
    }
    

  • 相关阅读:
    解决方案下显示的网站名称被追加编号的问题解决方法
    sql server 2008 报表服务器(sql server 2008 report service)中如何固定行或者列
    如何通过参数来切换图表和数据
    饼图中如何将一定比例的小切片收集在一起
    饼图图例中显示百分比值
    如何根据条件隐藏列
    饼图上如何显示百分比值
    负载均衡
    在UTF8页中接收和使用以GB2312方式进行URL编码的中文数据
    ASP 格式化显示时间为几个月,几天前,几小时前,几分钟前,或几秒前
  • 原文地址:https://www.cnblogs.com/Troverld/p/14607818.html
Copyright © 2020-2023  润新知