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}}$
- 就可以求了;
1 #include<bits/stdc++.h> 2 #define mod 998244353 3 #define Run(i,l,r) for(int i=l;i<=r;++i) 4 #define Don(i,l,r) for(int i=l;i>=r;--i) 5 const int N=1<<19; 6 using namespace std; 7 int n,a[N],fac[N],ans[N]; 8 char s[N],ps[1000000],*pp=ps; 9 int pw(int x,int y){ 10 if(y<0)y+=mod-1; 11 int re=1; 12 for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)re=1ll*re*x%mod; 13 return re; 14 } 15 void push(char x){ 16 if(pp==ps+1000000)fwrite(ps,1,1000000,stdout),pp=ps; 17 *pp++=x; 18 } 19 void write(int x){ 20 static int top,sta[20]; 21 if(!x){push('0'),push('\n');return;} 22 while(x)sta[++top]=x%10,x/=10; 23 while(top)push(sta[top--]^'0'); 24 push('\n'); 25 } 26 void flush(){fwrite(ps,1,pp-ps,stdout);} 27 namespace poly{ 28 int g=3,iv[N],rev[N],L; 29 void init(int l){for(int i=1;i<=l;++i)iv[i]=pw(i,mod-2);} 30 void cls(int*A,int l,int r){for(int i=l;i<r;++i)A[i]=0;} 31 void cpy(int*A,int*B,int l){for(int i=0;i<l;++i)A[i]=B[i];} 32 void der(int*A,int l){Run(i,0,l-2)A[i]=1ll*(i+1)*A[i+1]%mod;A[l-1]=0;} 33 void dif(int*A,int l){Don(i,l-1,1)A[i]=1ll*iv[i]*A[i-1]%mod;A[0]=0;} 34 void ntt(int*A,int l,int f){ 35 for(L=0;(1<<L)<l;++L); 36 for(int i=0;i<l;++i){ 37 rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 38 if(i<rev[i])swap(A[i],A[rev[i]]); 39 } 40 for(int i=1;i<l;i<<=1){ 41 int wn=pw(g,f*(mod-1)/(i<<1)); 42 for(int j=0;j<l;j+=i<<1){ 43 int w=1; 44 for(int k=0;k<i;++k,w=1ll*w*wn%mod){ 45 int x=A[j+k],y=1ll*w*A[j+k+i]%mod; 46 A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod; 47 } 48 } 49 } 50 if(!~f){for(int i=0;i<l;++i)A[i]=1ll*A[i]*iv[l]%mod;} 51 } 52 void inv(int*A,int*B,int l){ 53 static int t[N]; 54 if(l==1){B[0]=pw(A[0],mod-2);return;} 55 int len=l<<1; 56 inv(A,B,l>>1);cls(B,l,len); 57 cpy(t,A,l);cls(t,l,len); 58 ntt(t,len,1);ntt(B,len,1); 59 for(int i=0;i<len;++i)B[i]=1ll*B[i]*(2-1ll*t[i]*B[i]%mod+mod)%mod; 60 ntt(B,len,-1);cls(B,l,len); 61 } 62 void ln(int*A,int*B,int l){ 63 static int t[N]; 64 int len=l<<1; 65 inv(A,B,l); 66 cpy(t,A,l);cls(t,l,len); 67 der(t,l); 68 ntt(B,len,1);ntt(t,len,1); 69 for(int i=0;i<len;++i)B[i]=1ll*B[i]*t[i]%mod; 70 ntt(B,len,-1);cls(B,l,len); 71 dif(B,l); 72 } 73 void exp(int*A,int*B,int l){ 74 static int t[N]; 75 if(l==1){B[0]=1;return;} 76 int len=l<<1; 77 exp(A,B,l>>1);cls(B,l,len); 78 ln(B,t,l); 79 for(int i=0;i<l;++i)t[i]=(A[i]-t[i]+mod)%mod; 80 t[0]++; 81 ntt(B,len,1);ntt(t,len,1); 82 for(int i=0;i<len;++i)B[i]=1ll*B[i]*t[i]%mod; 83 ntt(B,len,-1);cls(B,l,len); 84 } 85 void solve(int*A,int l){ 86 static int t[N],r[N]; 87 if(l==1){A[0]=0;return;} 88 int len=l<<1; 89 solve(A,l>>1);cls(A,l,len); 90 cpy(t,a,l);cls(t,l,len); 91 ntt(A,len,1);ntt(t,len,1); 92 for(int i=0;i<len;++i){ 93 int tmp=A[i]; 94 A[i]=(mod-1ll*iv[2]*t[i]%mod*A[i]%mod*A[i]%mod)%mod; 95 t[i]=(mod-1ll*t[i]*tmp%mod)%mod; 96 } 97 ntt(A,len,-1);cls(A,l,len);A[0]++; 98 ntt(t,len,-1);cls(t,l,len);dif(t,l); 99 exp(t,r,l);inv(r,t,l); 100 ntt(A,len,1);ntt(r,len,1); 101 for(int i=0;i<len;++i)A[i]=1ll*A[i]*r[i]%mod; 102 ntt(A,len,-1);cls(A,l,len); 103 dif(A,l); 104 ntt(A,len,1);ntt(t,len,1); 105 for(int i=0;i<len;++i)A[i]=1ll*A[i]*t[i]%mod; 106 ntt(A,len,-1);cls(A,l,len); 107 } 108 } 109 int main(){ 110 // freopen("uoj50.in","r",stdin); 111 // freopen("uoj50.out","w",stdout); 112 scanf("%d%s",&n,s+1); 113 for(int i=fac[0]=1;i<=n;++i){ 114 fac[i]=1ll*fac[i-1]*i%mod; 115 if(s[i]-'0')a[i-1]=pw(fac[i-1],mod-2); 116 } 117 int len=1;for(;len<=n;len<<=1); 118 poly::init(len<<1); 119 poly::solve(ans,len); 120 for(int i=1;i<=n;++i)write(1ll*ans[i]*fac[i]%mod); 121 flush(); 122 return 0; 123 }