uoj50【UR#3】链式反应

uoj50【UR#3】链式反应

  • 题解:

    • 令$a(x)$为破坏死光的$EFG$,$f(x)$为方案的$EGF$:$f(x) = x + \int \  \frac{1}{2} f^2(x) a(x) \  dt$;
    • 注意到$f(0)=0$,所以考虑如何解:$f'(x) =  \frac{1}{2} a(x) f(x)^2 + 1$
    • 设$g(f) = 1 + \frac{1}{2}af^2$,即求解$f' = g(f)$;
    • 主要思想是牛顿迭代,假设已经求得$f  \equiv f_{0}   \pmod {x^{n}} $:
    • 泰勒展开得:$f' \equiv g(f_{0}) + g'(f_{0})(f - f_{0}) \pmod {x^{2n}} $
    • $f' - g'(f_{0})f \equiv g(f_{0}) - g'(f_{0})f_{0} $
    • 进一步令:$r \equiv e ^ {\int \  - g'(f_{0}) \ dt}$,(注意有$r' = -g'(f_{0}) r$)
    • 两边同时乘以$r$ 得到:
    • $f'r - fr' \equiv  (g(f_{0}) - g'(f_{0})f_{0}) r$
    • $(fr)' \equiv (g(f_{0}) - g'(f_{0})f_{0}) r$
    • $ f \equiv \frac{  \int \  (1 - \frac{1}{2}af_{0}^{2})r \ dt }{r} \pmod {x^{2n}}$
    • 就可以求了;
 #include<bits/stdc++.h>
#define mod 998244353
#define Run(i,l,r) for(int i=l;i<=r;++i)
#define Don(i,l,r) for(int i=l;i>=r;--i)
const int N=<<;
using namespace std;
int n,a[N],fac[N],ans[N];
char s[N],ps[],*pp=ps;
int pw(int x,int y){
if(y<)y+=mod-;
int re=;
for(;y;y>>=,x=1ll*x*x%mod)if(y&)re=1ll*re*x%mod;
return re;
}
void push(char x){
if(pp==ps+)fwrite(ps,,,stdout),pp=ps;
*pp++=x;
}
void write(int x){
static int top,sta[];
if(!x){push(''),push('\n');return;}
while(x)sta[++top]=x%,x/=;
while(top)push(sta[top--]^'');
push('\n');
}
void flush(){fwrite(ps,,pp-ps,stdout);}
namespace poly{
int g=,iv[N],rev[N],L;
void init(int l){for(int i=;i<=l;++i)iv[i]=pw(i,mod-);}
void cls(int*A,int l,int r){for(int i=l;i<r;++i)A[i]=;}
void cpy(int*A,int*B,int l){for(int i=;i<l;++i)A[i]=B[i];}
void der(int*A,int l){Run(i,,l-)A[i]=1ll*(i+)*A[i+]%mod;A[l-]=;}
void dif(int*A,int l){Don(i,l-,)A[i]=1ll*iv[i]*A[i-]%mod;A[]=;}
void ntt(int*A,int l,int f){
for(L=;(<<L)<l;++L);
for(int i=;i<l;++i){
rev[i]=(rev[i>>]>>)|((i&)<<(L-));
if(i<rev[i])swap(A[i],A[rev[i]]);
}
for(int i=;i<l;i<<=){
int wn=pw(g,f*(mod-)/(i<<));
for(int j=;j<l;j+=i<<){
int w=;
for(int k=;k<i;++k,w=1ll*w*wn%mod){
int x=A[j+k],y=1ll*w*A[j+k+i]%mod;
A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod;
}
}
}
if(!~f){for(int i=;i<l;++i)A[i]=1ll*A[i]*iv[l]%mod;}
}
void inv(int*A,int*B,int l){
static int t[N];
if(l==){B[]=pw(A[],mod-);return;}
int len=l<<;
inv(A,B,l>>);cls(B,l,len);
cpy(t,A,l);cls(t,l,len);
ntt(t,len,);ntt(B,len,);
for(int i=;i<len;++i)B[i]=1ll*B[i]*(-1ll*t[i]*B[i]%mod+mod)%mod;
ntt(B,len,-);cls(B,l,len);
}
void ln(int*A,int*B,int l){
static int t[N];
int len=l<<;
inv(A,B,l);
cpy(t,A,l);cls(t,l,len);
der(t,l);
ntt(B,len,);ntt(t,len,);
for(int i=;i<len;++i)B[i]=1ll*B[i]*t[i]%mod;
ntt(B,len,-);cls(B,l,len);
dif(B,l);
}
void exp(int*A,int*B,int l){
static int t[N];
if(l==){B[]=;return;}
int len=l<<;
exp(A,B,l>>);cls(B,l,len);
ln(B,t,l);
for(int i=;i<l;++i)t[i]=(A[i]-t[i]+mod)%mod;
t[]++;
ntt(B,len,);ntt(t,len,);
for(int i=;i<len;++i)B[i]=1ll*B[i]*t[i]%mod;
ntt(B,len,-);cls(B,l,len);
}
void solve(int*A,int l){
static int t[N],r[N];
if(l==){A[]=;return;}
int len=l<<;
solve(A,l>>);cls(A,l,len);
cpy(t,a,l);cls(t,l,len);
ntt(A,len,);ntt(t,len,);
for(int i=;i<len;++i){
int tmp=A[i];
A[i]=(mod-1ll*iv[]*t[i]%mod*A[i]%mod*A[i]%mod)%mod;
t[i]=(mod-1ll*t[i]*tmp%mod)%mod;
}
ntt(A,len,-);cls(A,l,len);A[]++;
ntt(t,len,-);cls(t,l,len);dif(t,l);
exp(t,r,l);inv(r,t,l);
ntt(A,len,);ntt(r,len,);
for(int i=;i<len;++i)A[i]=1ll*A[i]*r[i]%mod;
ntt(A,len,-);cls(A,l,len);
dif(A,l);
ntt(A,len,);ntt(t,len,);
for(int i=;i<len;++i)A[i]=1ll*A[i]*t[i]%mod;
ntt(A,len,-);cls(A,l,len);
}
}
int main(){
// freopen("uoj50.in","r",stdin);
// freopen("uoj50.out","w",stdout);
scanf("%d%s",&n,s+);
for(int i=fac[]=;i<=n;++i){
fac[i]=1ll*fac[i-]*i%mod;
if(s[i]-'')a[i-]=pw(fac[i-],mod-);
}
int len=;for(;len<=n;len<<=);
poly::init(len<<);
poly::solve(ans,len);
for(int i=;i<=n;++i)write(1ll*ans[i]*fac[i]%mod);
flush();
return ;
}
上一篇:使用JS对字符串进行MD5加密


下一篇:poj 2888 Magic Bracelet(Polya+矩阵快速幂)