题解P6613

原题链接

题意:
求满足\(\frac{d F(x)}{d x}=A(x)\times e^{F(x)-1}+B(x)\)\(F(x)\)的前\(n\)项系数,且\(F(0)=1\)

强解微分方程!

\(y=F(x),A1(x)=\int^{}{}{A(x)d x},B1(x)=\int^{}{}{B(x)d x}\)
先求\(\frac{d y}{d x}=A(x)\times e^{y-1}\)的解:
\(\int^{}{}{e^{1-y}dy}\)=\(\int^{}{}{A(x)d x}-C\)(\(C为常数\))
\(-e^{1-y}=A1(x)-C\)
\(y=1+\ln{\frac{1}{C-A1(x)}}\)

再由常数变易法设\(C=u(x)\)
\(\frac{d y}{d x}=(u(x)-A1(x))(\frac{-1}{(u(x)-A1(x))^2})(\frac{du(x)}{dx}-A(x))\)
\(=\frac{A(x)-\frac{du(x)}{dx}}{u(x)-A1(x)}\)
\(A(x)e^{F(x)-1}+B(x)=\frac{A(x)}{u(x)-A1(x)}+B(x)\)
所以\(\frac{du(x)}{dx}=A1(x)B(x)-u(x)B(x)\)

再一次,先求\(\frac{du(x)}{dx}=-u(x)B(x)\)的解:
\(\int^{}{}{\frac{1}{u}du}=-\int^{}{}{B(x)d x}+D\)(\(D\)为常数)
\(\ln(u(x))=-B1(x)+D\)
\(u(x)=e^{-B1(x)}\times E\)(常数\(E=e^D\)

再设\(E=v(x)\)
\(\frac{du(x)}{dx}=-v(x)B(x)e^{-B1(x)}+\frac{dv(x)}{dx}e^{-B1(x)}\)
\(A1(x)B(x)-u(x)B(x)=-v(x)B(x)e^{-B1(x)}+A1(x)B(x)\)
于是得到\(\frac{dv(x)}{dx}=e^{B1(x)}A1(x)B(x)\)
因此\(v(x)=\int^{}{}{e^{B1(x)}A1(x)B(x)dx}+G\)(\(G\)为常数)

由于\(C=u(x)=e^{-B1(x)}v(x)\)
所以\(F(x)=y=1+\ln{\frac{1}{e^{-B1(x)}v(x)-A1(x)}}=1+B1(x)-\ln(v(x)-A1(x)e^{B1(x)})\)
\(=1+B1(x)-\ln(\int^{}{}{e^{B1(x)}A1(x)B(x)dx}-A1(x)e^{B1(x)}+G)\)
\(F(0)=1\),取\(G=1\)即可

代码:

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+10;
const int mod=998244353;
int n,g;
int a[N],b[N];
int a1[N],b1[N],eb[N];
int c1[N],c2[N],c[N],d[N],dd[N];
int tmp[N],tmp1[N];
int tt[N],t1[N];
int tln[N],tinv[N];
int rev[N];
int inv[N];
void getinv(int n){
	inv[1]=1;
	for(int i=2;i<n;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
}
int getrev(int n){
	int n1=1,len=0;
	while(n1<(n<<1))n1<<=1,len++;
	for(int i=0;i<n1;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	return n1;
}
int qpow(int x,int k){
	int cnt=1;
	while(k){
		if(k&1)cnt=(cnt*x)%mod;
		x=(x*x)%mod;
		k>>=1;
	}
	return cnt;
}
void integral(int n,int *x,int *y){
	for(int i=1;i<n;i++)y[i]=x[i-1]*inv[i]%mod;
	y[0]=0;
} 
void ntt(int n,int *a,int t){
	for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	if(t==1)g=3;
	else g=332748118;
	for(int mid=1;mid<n;mid<<=1){
		int wn=qpow(g,(mod-1)/(mid<<1));
		for(int j=0;j<n;j+=(mid<<1)){
			int w=1;
			for(int k=0;k<mid;k++){
				int x=a[j+k],y=w*a[j+k+mid]%mod;
				a[j+k]=(x+y)%mod;
				a[j+k+mid]=(x-y+mod)%mod;
				w=w*wn%mod;
			}
		}
	}
	if(t==1)return;
	int inv=qpow(n,mod-2);
	for(int i=0;i<n;i++)a[i]=a[i]*inv%mod;
}
void polyinv(int n,int *f,int *g){
	if(n==1){
		g[0]=qpow(f[0],mod-2);
		return;
	}
	polyinv((n+1)>>1,f,g);
	int n1=getrev(n);
	for(int i=0;i<n;i++)tmp[i]=f[i];
	for(int i=n;i<n1;i++)tmp[i]=0;
	ntt(n1,tmp,1);
	ntt(n1,g,1);
	for(int i=0;i<n1;i++)g[i]=(2ll-tmp[i]*g[i]%mod+mod)%mod*g[i]%mod;
	ntt(n1,g,-1);
	for(int i=n;i<n1;i++)g[i]=0;
}
void polyln(int n,int *ff,int *gg){
	int n1=getrev(n);
	for(int i=0;i<n;i++)tt[i]=(i+1)*ff[i+1]%mod;
	for(int i=n;i<n1;i++)tt[i]=0,gg[i]=0;
	polyinv(n,ff,gg);
	ntt(n1,gg,1);
	ntt(n1,tt,1);
	for(int i=0;i<n1;i++)tt[i]=tt[i]*gg[i]%mod;
	ntt(n1,tt,-1);
	integral(n,tt,gg);
	for(int i=n;i<n1;i++)gg[i]=0;
}
void polyexp(int n,int *f1,int *g1){
	if(n==1){
		g1[0]=1;
		return;
	}
	polyexp((n+1)>>1,f1,g1);
	int n1=getrev(n);
	for(int i=0;i<n;i++)tmp1[i]=f1[i],tln[i]=0;
	for(int i=n;i<n1;i++)tmp1[i]=0,tln[i]=0,g1[i]=0;
	polyln(n,g1,tln);
	tmp1[0]=1;
	for(int i=0;i<n;i++)tmp1[i]=(tmp1[i]-tln[i]+mod)%mod;
	ntt(n1,g1,1);
	ntt(n1,tmp1,1);
	for(int i=0;i<n1;i++)g1[i]=g1[i]*tmp1[i]%mod;
	ntt(n1,g1,-1);
	for(int i=n;i<n1;i++)g1[i]=0;
}
signed main(){
	scanf("%d",&n),n++;
	for(int i=0;i<n;i++)scanf("%d",&a[i]);
	for(int i=0;i<n;i++)scanf("%d",&b[i]);
	int n1=getrev(n);
	getinv(n);
	integral(n,a,a1),integral(n,b,b1);
	polyexp(n,b1,eb);
	ntt(n1,eb,1),ntt(n1,a1,1);
	for(int i=0;i<n1;i++)c2[i]=a1[i]*eb[i]%mod;
	ntt(n1,c2,-1);
	for(int i=n;i<n1;i++)c2[i]=0;
	ntt(n1,c2,1);
	ntt(n1,b,1);
	for(int i=0;i<n1;i++)c1[i]=c2[i]*b[i]%mod;
	ntt(n1,c1,-1);
	integral(n,c1,d);
	ntt(n1,c2,-1);
	for(int i=0;i<n;i++)d[i]=(d[i]-c2[i]+mod)%mod;
	d[0]=1;
	polyln(n,d,dd);
	for(int i=0;i<n;i++)c[i]=(b1[i]-dd[i]+mod)%mod;
	c[0]=(c[0]+1)%mod;
	for(int i=0;i<n;i++)printf("%d ",c[i]);
}
posted @ 2021-01-16 00:23  Y_B_X  阅读(91)  评论(0编辑  收藏  举报