[数学记录]Uoj#50. 【UR #3】链式反应

· · 个人记录

题意 : 某种原子核 A ,受到中子撞击之后,有可能裂变,也有可能退出反应。

裂变会产生恰好 2 个中子，且发出光子照射使得 c 个原子核退出反应。

其中，c 要满足 c∈S。

每个原子核只能影响编号比他大的原子核，现在有一排 n 个原子核，第一个原子核被中子集中，且反应完成后没有剩余的完好原子核。问可能的反应情况数。

两种反应情况不同，当且仅当某个中子裂变时，对其他中子的影响不同(中子和光子)

答案对 998244353 取模，n\leq 2\times 10^5 ，时限\texttt{4s}.

当原子核 u 裂变时影响到的其他原子核向 u 连边，作为其儿子，构建一棵树。

这棵有标号有根树满足儿子的标号大于父亲。

而且，每个非叶节点有 2+c 个儿子,其中 c 个是叶子且满足 c∈S ，另外两个可以是也可以不是叶子。

设 f[i] 为 i 个点树的方案数。

考虑 DP ,枚举两个不一定是叶子的子树大小,并分配标号。

f[i]=\dfrac{1}{2}\sum\limits_{j,k}[i-1-j-k∈S]\dbinom{i-1}{j}\dbinom{i-1-j}{k}f[j]f[k]

由于两个子树之间无序，所以要 \frac{1}{2} 。

把组合数拆一下,能得到 :

f[i]=\dfrac{1}{2}\sum\limits_{j,k}[i-1-j-k∈S]\dfrac{(i-1)!(i-1-j)!}{j!k!(i-1-j)!(i-1-j-k)!}f[j]f[k] \dfrac{f[i]}{(i-1)!}=\dfrac{1}{2}\sum\limits_{j,k}\dfrac{[i-1-j-k∈S]}{(i-1-j-k)!}\dfrac{f[j]f[k]}{j!k!}

似乎长得挺像 EGF ? 实际上直接考虑组合意义估计也能得到吧。

设 S(x) 为 [i∈S] 的 EGF。 F(x) 为 f[i] 的 EGF。

则有 F'(x)=\dfrac{1}{2}S(x)F(x)^2+1 (考虑清楚常数项)

暂不会牛顿迭代解一阶微分方程……而且看起来也不太实用？

我们也可以将其转化成半在线卷积，就可以 O(n\log^2n) 分治 FFT 做了。更好写，效率也不差。

提取系数得 (n+1)F[n+1]=\dfrac{1}{2}\sum\limits_{i=0}^nS[i]F^2[n-i]

即 F[n]=\dfrac{1}{2n}\sum\limits_{i=1}^{n}S[i-1]F^2[n-i]。

具体怎么做可见 半在线卷积小记。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=1<<18|500;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn<<1]={1};
  for (int i=0;i<n;i++)f[i]=g[tr[i]];
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int F[Maxn],F2[Maxn],G[Maxn],s1[Maxn],s2[Maxn],s3[Maxn];
void cdq(int l,int r)
{
  if (r==1)return ;
  if (l+1==r){
    if (l^1)F[l]=powM(2*l)*F[l]%mod;
    return ;
  }int mid=(l+r)>>1,n=r-l;
  cdq(l,mid);
  if (!l){
    cpy(s1,F,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    px(s1,s1,n);NTT(s1,0,n);
    for (int i=n/2;i<n;i++)
      F2[i]=(F2[i]+s1[i])%mod;
    cpy(s1,F2,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,G,n/2);clr(s2+(n/2),n/2);NTT(s2,1,n);
    px(s1,s2,n);NTT(s1,0,n);
    for (int i=n/2;i<n;i++)
      F[i]=(F[i]+s1[i])%mod;
  }else {
    cpy(s1,F+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,F,n);NTT(s2,1,n);
    px(s1,s2,n);NTT(s1,0,n);
    for (int i=n/2;i<n;i++)
      F2[l+i]=(F2[l+i]+2ll*s1[i])%mod;
    cpy(s1,G+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,F2,n);NTT(s2,1,n);px(s1,s2,n);
    for (int i=0;i<n;i++)s3[i]=s1[i];
    cpy(s1,F2+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
    cpy(s2,G,n);NTT(s2,1,n);px(s1,s2,n);
    for (int i=0;i<n;i++)s3[i]=(s3[i]+s1[i])%mod;
    NTT(s3,0,n);
    for (int i=n/2;i<n;i++)
      F[l+i]=(F[l+i]+s3[i])%mod;
  }cdq(mid,r);
}
ll fac[Maxn],ifac[Maxn];
void Init(int n)
{
  fac[0]=1;
  for (int i=1;i<=n;i++)
    fac[i]=fac[i-1]*i%mod;
  ifac[n]=powM(fac[n]);
  for (int i=n;i;i--)
    ifac[i-1]=ifac[i]*i%mod;
}
int n,m;
char s[Maxn];
int main()
{
  scanf("%d%s",&m,s);
  Init(++m);F[1]=1;
  for (int i=0;i<m;i++)
    if (s[i]=='1')G[i+1]=ifac[i];
  for (n=1;n<m;n<<=1);
  cdq(0,n);
  for (int i=1;i<m;i++)
    printf("%lld\n",fac[i]*F[i]%mod);
  return 0;
}