LOJ#6538. 烷基计数 加强版 加强版/LG P6598 烷烃计数

Description

$n$ 个碳原子的烷烃共有多少种同分异构体？

提示：如果你不知道什么是烷烃，那么你可以认为这个问题等价于求 $n$ 个点的无标号无根树并满足每个点的度数 $\leq 4$ 的树的个数。

Solution

先对烷基计数，即求$n$个碳原子的烷基共有多少种

先忽视烷基中单独的那个电子，那么烷基可以视为三叉树，如果有不饱和的C，可以看做其有一个大小为0的子树

根据burnside引理，设$n$个C的烷基有$f_n$个，那么：

$$f_n=\frac{\sum _{i+j+k=n-1}f_if_jf_k+3\sum _{2i+j=n-1}f_if_j+2\sum _{3i=n-1}f_i}{6} $$

设$A(x)$是其生成函数，那么：

$$A(x)=\frac{A^3(x)+3A(x^2)A(x)+2A(x^3)}{6}x+1 $$

递推$A$可以用分治FFT或者牛顿迭代解决

再尝试证明结论：

在一个已知的烷烃（一个无根树）中，设$p$为将一个烷烃中的C特殊标记后，有多少不同的烷烃，$q$为将一个烷烃中的碳碳键特殊标记后，有多少不同的烷烃，$b$为该烷烃是否有两个重心且两重心两侧同构的布尔变量，那么：

$$p-q+b=1$$

分类讨论，当$b=false$时，任选一个重心，该重心被标记后必不会与其它点标记的情况同构，如果有两个C被标记后同构，那么这两个C向父亲方向（以重心为根）的碳碳键被标记后一定同构

当$b=true$时，除重心外的C同上，两个重心只产生一个同构，两重心之间连边产生一个同构，所以$p=q$

那么枚举所有的无根树，就有$\sum p-\sum q+\sum b=\sum 1$，其实题中就是求$\sum 1$

所以现在要计算$\sum p$，$\sum q$和$\sum b$

计算$\sum p$：

　　就是计算有根树的个数（将被标记的C看作根），用与计算烷基相似的方法计算有根树的个数，最终得到OGF$P(x)$为：

·　　$$P(x)=\frac{A^4(x)+6A(x^2)A^2(x)+3A^2(x^2)+8A(x^3)A(x)+6A(x^4)}{24}x+1 $$

计算$\sum q$：

　　可以看做让两个烷基之间成键，最终得到OGF$Q(x)$为：

　　$$Q(x)=\frac{(A^2(x)-1)^2+A(x^2)-1}{2}+1 =\frac{A^2(x)+A(x^2)-2A^2(x)}{2}+1 $$

　　在$A(x)$平方时将每个情况计算了两次，但是两个烷基相同的情况只计算了一次，所以写成这样

计算$\sum b$

　　发现就是将两个相同的烷基之间成键，OGF为$A(x^2)$

答案OGF为$P(x)-Q(x)+A(x^2)$

#include<iostream>
#include<cstdio>
using namespace std;
int s,tot,rev[800005],n,T;
const int mod=998244353;
long long A[800005],inv[800005],A2[800005],A3[800005],ans[800005],C[800005];
inline long long read(){
    long long f=1,w=0;
    char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')w=(w<<1)+(w<<3)+ch-'0',ch=getchar();
    return f*w;
}
long long ksm(long long a,long long p){
    long long ret=1ll;
    while(p){
        if(p&1)(ret*=a)%=mod;
        (a*=a)%=mod,p>>=1;
    }
    return ret;
}
void ntt(long long *a,int N,int in){
    for(int i=0;i<N;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=1;i<N;i<<=1){
        long long wn=ksm(3,(mod-1)/i/2);
        if(in==-1)wn=ksm(wn,mod-2);
        for(int j=0;j<N;j+=i<<1){
            long long w=1;
            for(int k=j;k<i+j;k++){
                long long x=a[k],y=w*a[k+i]%mod;
                a[k]=(x+y)%mod,a[k+i]=(mod+x-y)%mod,(w*=wn)%=mod;
            }
        }
    }
    if(in==-1)for(int i=0;i<N;i++)(a[i]*=inv[N])%=mod;
}
void solve(int l,int r){
    if(l==r){(A[l]+=(l%3==1)*A[l/3]*inv[3]%mod)%=mod;return;}
    int mid=(l+r)>>1;
    solve(l,mid),s=2,tot=1;
    static long long a1[800005],a2[800005],b[800005],c[800005];
    while(s<=(2*(r-l+1)))s<<=1,++tot;
    for(int i=0;i<s;i++)a1[i]=i<=mid-l?A[l+i]:0,a2[i]=i<=r-l?A[i]:0,b[i]=(i<=r-l&&!(i&1))?A[i/2]:0;
    for(int i=0;i<s;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    ntt(a1,s,1),ntt(a2,s,1),ntt(b,s,1);
    for(int i=0;i<s;i++)c[i]=(a1[i]*a2[i]%mod*a2[i]%mod*(l?inv[2]:inv[6])%mod+a1[i]*b[i]%mod*inv[2]%mod)%mod;
    ntt(c,s,-1);
    for(int i=mid-l;i<=r-l-1;i++)(A[l+i+1]+=c[i])%=mod;
    solve(mid+1,r);
}
int main(){
    inv[0]=inv[1]=1,T=read();
    for(int i=2;i<=800000;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    A[0]=1,solve(0,100000),s=2,tot=1;
    while(s<=300000)s<<=1,++tot;
    for(int i=0;i<=100000;i++)A2[i]=i&1?0:A[i/2],A3[i]=i%3?0:A[i/3];
    for(int i=0;i<=100000;i++)if(i%4==1)ans[i]=A[i/4]*inv[4]%mod;
    for(int i=0;i<=100000;i++)if(!(i&1))ans[i]=A[i/2];
    for(int i=0;i<s;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    ntt(A,s,1),ntt(A2,s,1),ntt(A3,s,1);
    for(int i=0;i<s;i++)C[i]=(A2[i]*A2[i]%mod*inv[8]%mod+A[i]*A[i]%mod*A2[i]%mod*inv[4]%mod+A[i]*A3[i]%mod*inv[3]%mod+A[i]*A[i]%mod*A[i]%mod*A[i]%mod*inv[24]%mod)%mod;
    ntt(C,s,-1);
    for(int i=0;i<=100000;i++)(ans[i+1]+=C[i])%=mod;
    for(int i=0;i<s;i++)C[i]=(A[i]*A[i]%mod*inv[2]%mod+A2[i]*inv[2]%mod-A[i]+mod)%mod;
    ntt(C,s,-1);
    for(int i=0;i<=100000;i++)(ans[i]+=mod-C[i])%=mod;
    for(;T;T--)n=read(),printf("%lld\n",ans[n]);
    return 0;
}