题解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]);
}