[UOJ50] 链式反应

题面

高质量好题。

写的时候脑子一片混乱,内容狗屁不通

算法一

首先根据题面,可以设 \(f_{i}\) 表示序列长度为 \(i\) 的时候的方案数,其中 \(f_{1}=1\) 。

那么有 \(f_{i}=\sum_{j=1}\sum_{k=1}[s_{i-j-k-1}=1]f_{j}f_{k}\binom{i-1}{i-j-k-1}\binom{j+k}{j}\)

两个组合数分别是选出 \(i-j-k-1\) 个点被删除和给两个剩下来的组重标号。

但是注意到这东西实际上会算重,因为 \((j,k)\) 和 \((k,j)\) 这个方案被计算了两次,而实际上它们是同一种方案,所以答案还要除以 \(2\) 。

虽然 \((j,j)\) 这个二元组只被算了一次,但是可以发现后面重标号部分使得这个地方被算了两次,刚好抵消了。所以直接除以 \(2\) 是没问题的。

于是我们就得到了一个 \(O(n^3)\) 的做法。

算法二

对上面的式子做一个变形,可以得到 :

\[\begin{aligned} \frac{2f_{i}}{(i-1)!}=\sum_{j=1}\sum_{k=1}[s_{i-j-k-1}=1]\frac{f_{j}}{j!}\frac{f_{k}}{k!}\frac{1}{(i-j-k-1)!} \end{aligned} \]

如果我们令 \(g_{t}=\sum_{i=1}\dfrac{f_{i}}{i!}\dfrac{f_{t-i}}{(t-i)!}\)

那么有

\[\begin{aligned} \frac{2f_{i}}{(i-1)!}=\sum^{i-1}_{t=0}[s_{i-t-1}=1]g_{t}\frac{1}{(i-t-1)!} \end{aligned} \]

就目前为止,已经可以做到 \(O(n^2)\) 了。

算法三

如果我们令 \(S_{i}=[s_{i}=1]\frac{1}{i!}\),那么有

\[\begin{aligned} \frac{2f_{i}}{(i-1)!}&=\sum^{i-1}_{t=0}S_{i-t-1}g_{t}\\ g_{t}&=\sum^{t}_{i=0}f_{i}f_{t-i} \end{aligned} \]

可以发现均为卷积的形式,这种半在线的卷积考虑用 CDQ 来优化。

为了方便,这里定义 \(F(x)_{l,r}\) 表示\(\sum^{r}_{i=l}f_{i}x^i\)

先看 $ F(x)$ 的部分,假设我们已经解决到了区间 \([l,r]\) 的部分,其中 \([l,mid]\) 已经解决,要计算 \([mid+1,r]\) 部分的答案:

\[\begin{aligned} S(x)G(x)_{l,r}&=S(x)G(x)_{l,mid}+S(x)G(x)_{mid+1,r}\\ \end{aligned} \]

看上去我们需要\([0,n-1]\) 的 \(S(x)\) ,但是实际上我们真正需要的部分是 \([mid+1,r]\) ,有很多地方是浪费的,尝试标一下上下界。

\[\begin{aligned} S(x)G(x)_{l,r}&=S(x)_{1,r-l}G(x)_{l,mid}+S(x)_{}G(x)_{mid+1,r}\\ \end{aligned} \]

左边部分可以 \(O(L\log L)\) 解决,右边部分是一个子问题,递归解决。

同样的,我们来看 \(G(x)\) 的部分:

\[\begin{aligned} F(x)F(x)_{l,r}&=F(x)F(x)_{l,mid}+F(x)F(x)_{mid+1,r}\\ &=F(x)_{1,r-l}F(x)_{l,mid}+F(x)F(x)_{mid+1,r}\\ \end{aligned} \]

这里似乎出现了问题, \(r-l\) 有可能比 \(mid\) 大,答案还没有计算出来。

可以发现此时的 \(l\) 一定等于 \(1\)。 (不难证明,自己可以推推看)

\[\begin{aligned} F(x)_{1,r-l}F(x)_{1,mid}&=F(x)_{1,mid}F(x)_{1,mid}+F(x)_{mid+1,r-1}F(x)_{1,mid}\\ \end{aligned} \]

前面一部分已经可以计算了,后面一部分把它丢到下一次递归里再计算。

也就是在 \(solve(mid+1,r)\) 的时候还需要计算 \(F(x)_{1,mid}F(x)_{mid+1,r-1}\)

这个东西就是在正常情况中 \(F(x)_{1,l-1}F(x)_{l,r-1}\) 。

继续把这个东西拆成\(F(x)_{1,l-1}F(x)_{l,mid}+F(x)_{1,l-1}F(x)_{mid+1,r-1}\)

左边这部分就是直接让 \(l\) 前面的系数乘二,右边那部分和\(F(x)_{1,mid}F(x)_{mid+1,r-1}\) 是等价的(因为多出来的部分不会被统计)。

这样子之后恰好就是每一个部分 \(l\) 前面的系数乘二。

Code
#include<bits/stdc++.h>
using namespace std;
#define il inline
#define ri register int
#define ll long long
#define ui unsigned int
il ll read(){
    bool f=true;ll x=0;
    register char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=false;ch=getchar();}
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(f) return x;
    return ~(--x);
}
il void write(const ll &x){if(x>9) write(x/10);putchar(x%10+'0');}
il void print(const ll &x) {x<0?putchar('-'),write(~(x-1)):write(x);putchar('\n');}
il ll max(const ll &a,const ll &b){return a>b?a:b;}
il ll min(const ll &a,const ll &b){return a<b?a:b;}
namespace Poly{
    const ll mod=998244353,inv2=(mod+1)/2;
    ll ksm(ll d,ll t){
        ll res=1;
        for(;t;t>>=1){
            if(t&1) res=res*d%mod;
            d=d*d%mod;
        }
        return res;
    }
    il ll dec(ll x,ll y) {return x>=y?x-y:x-y+mod;}
    il ll add(ll x,ll y) {return x+y<mod?x+y:x+y-mod;}
    const int N=262144*2;
    vector<ll> w[23];
    ll fac[N],ifac[N],inv[N];
    il ll C(ll x,ll y){
        if(x<y) return 0;
        return fac[x]*ifac[y]%mod*ifac[x-y]%mod;
    }
    void init(){
        int n=N-1;
        fac[0]=1;
        for(ri i=1;i<=n;++i) fac[i]=i*fac[i-1]%mod;
        ifac[n]=ksm(fac[n],mod-2);
        for(ri i=n-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%mod;
        for(ri i=n;i;--i) inv[i]=ifac[i]*fac[i-1]%mod;
        inv[0]=1;
        int d=log(N)/log(2)+0.5;
        for(ri i=1;i<=d;++i){
            w[i].resize(1<<i);
            w[i][0]=1,w[i][1]=ksm(3,(mod-1)>>i);
            for(ri j=2;j<(1<<i);++j) w[i][j]=w[i][j-1]*w[i][1]%mod;
        }
    }
    int r[N];
    void DFT(int limit,ll *a,int flag){
        for(ri i=1;i<limit;++i) if(r[i]<i) swap(a[i],a[r[i]]);
        for(ri l=1,t=1;l<limit;l<<=1,++t){
            for(ri i=0;i<limit;i+=l<<1){
                ll *W=&w[t][0];
                for(ri j=0;j<l;++j){
                    ll tmp=a[i+j+l]*(*W++)%mod;
                    a[i+j+l]=dec(a[i+j],tmp);
                    a[i+j]=add(a[i+j],tmp);
                }
            }
        }
        if(flag==-1){
            reverse(a+1,a+limit);
            ll inv=ksm(limit,mod-2);
            for(ri i=0;i<limit;++i) a[i]=a[i]*inv%mod;
        }
    }
    ll _F[N];
    void rev(int limit,int ws){
        for(ri i=0;i<limit;i+=2){
            r[i]=r[i>>1]>>1;
            r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
        }
    }
    void NTT(int n,ll *f,int m,ll *g,int lim=0){
        if(!lim) lim=n+m;
        int limit=1,ws=0;
        for(;limit<=n+m;++ws,limit<<=1);
        rev(limit,ws);
        for(ri i=0;i<=m;++i) _F[i]=g[i];
        for(ri i=m+1;i<limit;++i) _F[i]=0;
        for(ri i=n+1;i<limit;++i) f[i]=0;
        DFT(limit,f,1),DFT(limit,_F,1);
        for(ri i=0;i<limit;++i) f[i]=f[i]*_F[i]%mod;
        DFT(limit,f,-1);
        for(ri i=lim+1;i<limit;++i) f[i]=0;
    }
}
using namespace Poly;
ll f[N],g[N],s[N],F[N],G[N],S[N];
int n;
/*
0 l+1
mid-l mid+1
r-l-1 r
*/
void solve(int l,int r){
    if(l==r){
        if(l==1) f[1]=1;
        return;
    }
    int mid=l+r>>1;
    solve(l,mid);
    for(ri i=0;i<=(r-l)*2;++i) S[i]=G[i]=0;
    for(ri i=0;i<=r-l-1;++i) S[i]=s[i];
    for(ri i=0;i<=mid-l;++i) G[i]=g[l+i];
    NTT(mid-l,G,r-l-1,S);
    for(ri i=mid-l;i<r-l;++i) f[i+l+1]=add(f[i+l+1],G[i]*inv[i+l+1]%mod*inv2%mod);
    for(ri i=0;i<=(r-l)*2;++i) S[i]=G[i]=0;
    for(ri i=0;i<=r-l-1;++i){
        if(i+1<l) S[i]=add(f[i+1],f[i+1]);
        else if(i+1<=mid) S[i]=f[i+1];
    }
    for(ri i=0;i<=mid-l;++i) G[i]=f[l+i];
    NTT(mid-l,G,r-l+1,S);
    for(ri i=mid-l;i<r-l;++i) g[i+l+1]=add(g[i+l+1],G[i]);
    solve(mid+1,r);
}
char ch[N];
int main(){
    // freopen("ex_reaction3.in","r",stdin);
    // freopen("1.out","w",stdout);
    init();
    n=read();
    scanf("%s",ch);
    for(ri i=0;i<n;++i) 
        if(ch[i]=='1') s[i]=ifac[i];
    solve(1,n);
    for(ri i=1;i<=n;++i) print(f[i]*fac[i]%mod);
    return 0;
}
上一篇:hdu 4736 This Is The Job The Bear Finds(2013年成都ACM网络赛)


下一篇:TODO、FIXME和XXX转载