烷烃计数

正式而卤蛋 · 2022-01-21 10:24:33 · 个人记录

前置知识

~~有机化学~~，生成函数，多项式，牛顿迭代，群论基础。

引入

由于不考虑立体异构，问题可以抽象为统计无标号无根树每个点度数 \le4 的本质不同个数。

下文将由简到难逐渐解决该问题，讲的稍微细一点因为我是萌新。

约定

为了方便表述，下文规定用 F(x) 表示烷基的 \mathrm{OGF} ，G(x) 表示烯烃的 \mathrm{OGF} ，A(x) 表示烷烃的 \mathrm{OGF} 。

烷基计数

根据有机化学基础，烷基计数可以抽象为统计无标号有根树每个点儿子数 \le 3 的本质不同个数。

考虑树形 \mathrm{dp} ，f_i 表示 i 个点的有根树的方案数（默认 f_0=1 方便 \mathrm{dp} ）。

考虑怎么转移出 f_n ，发现需要去重计算本质不同的个数，考虑用 \mathrm{Burnside}引理，我们把三个子树能凑出树除根 n-1 个点的方案看作染色集，而去重即三个子树的交换看作置换群，根据公式 |X/G|=\frac1{|G|}\sum_{g\in G}|X^g| ，我们需要计算在 6 种置换下染色集的不动点的和，实际发现置换可以分为 3 类：

$$\sum_{i=0}^{n-1}\sum_{j=0}^{n-1-i}f_i\times f_j\times f_{n-1-i-j}$$
$$4\times\sum_{i=0}^{\lfloor\frac{n-1}2\rfloor}f_i\times f_{n-1-2\times i}$$
$$2\times[3\mid(n-1)]\sum_{i=0}^{\frac{n-1}3}f_i$$

背包 \mathrm{dp} 类似于卷积形式，不难想到生成函数，无标号的用 \mathrm{OGF} ，有标号的用 \mathrm{EGF} 。于是构造 f 的 \mathrm{OGF} 为 F ，根据上面不难得出：

F(x)=1+x\ast\frac{F^3(x)+4\times F(x^2)F(x)+2\times F(x^3)}6

考虑用牛顿迭代求 F(x) 多项式的前 n 项，牛顿迭代其实是一个分治算法。设（就是上面那式子）：

H(F(x))=-F(x)+1+x\ast\frac{F^3(x)+4\times F(x^2)\ast F(x)+2\times F(x^3)}6

我们就是要求 H 的零点，假设我们已经求出了 F(x) 的前 \lceil\frac n2\rceil 项 F_0(x) ，即 F_0(x)\equiv F(x)\pmod{\lceil\frac n2\rceil} ，考虑现在怎么求出前 n 项。

根据泰勒公式有 f(x)=f(x_0)+f'(x_0)(x-x_0)+\cdots ，不难发现对于取模多项式其实只有前两项有用，即：

H(F(x))\equiv H(F_0(x))+H'(F_0(x))(F(x)-F_0(x))\pmod n

我们带入 H(F(x))=0 ，有：

\therefore F(x)\equiv F_0(x)-\frac{H(F_0(x))}{H'(F_0(x))}

于是就可以进行牛顿迭代分治了，边界 F(x)\equiv f_0=1\pmod {x^1} ，注意 H'(F_0(x)) 是对 F_0(x) 求导，即把 x,F_0(x^2),F_0(x^3) 看作常数即可。复杂度 O(n\log^2n) 。

烯烃计数

其实是类似的，还是一个有根问题。把双链看成关键边，那两边是两个烷基，组合去重就是两个儿子的置换，有：

G(x)=\frac{F^2(x)+F(x^2)}2

烷烃计数

问题变成一个无根树。把无根树转换成有根树计数：钦定重心为根最后去除子树大小不合法和重心重复。

具体的，先统计以重心为根的方案，即根可以有 4 个儿子，而子树结点的儿子度数 \le3 ，即烷基的方案，最后合并即 4 个儿子的置换，类似的得出：

A'(x)=1+x\ast\frac{F^4(x)+6\,F(x^2)F^2(x)+3\,F^2(x^2)+8\,F(x^3)F(x)+6\,F(x^4)}{24}

当前的答案就是 a'_n（注意这里是还没有去除不合法和重心重复）。

然后考虑不合法，即子树大小 >\lfloor\frac n2\rfloor ，暴力枚举一下大小，即减去：\sum_{i=\lfloor\frac n2\rfloor+1}^{n-1}f_i\times f_{n-i} 。

还要考虑重心重复，不难发现重复当且仅当有两个重心且两个重心子树本质不同，所以再减去：[2\mid n]\,C_{f_\frac n2}^2 。

于是：

Ans=a'_n-\sum_{i=\lfloor\frac n2\rfloor+1}^{n-1}f_i\times f_{n-i}-[2\mid n]\,C_{f_\frac n2}^2

这样复杂度是预处理 O(n\log^2n) ，单次询问 O(n) 的，可以过掉询问较少的烷烃计数。

而会有更毒瘤的问法让你算 1-100000 的所有答案，即求出生成函数 A(x) 。

我们换一种统计答案的方法，对于任意一个无根树，我们定义无根树的点等价类数为以该树每一个点为根形成的本质不同的有根树的个数，记作 p ，边等价类树同理，记作 q ，再记 s 表示这棵树是否有两个等价的重心，于是有结论 p-q+s=1 。

证明: 我们不妨以重心为根（这样最大化保证了等价点与其对应的等价边的对称性），每个点的贡献与其上的边对应。发现如果有两个点属于同一个等价类，那么其上连的边也就属于一个等价类。而重心单独讨论。若 s=0 即只有一个重心或两个重心不等价，则重心没有对应的边对应，而 s=1 两个重心等价贡献 -1 ，相当于两个重心与其连边对应。

所以答案为：\sum p-q+s ，枚举的是本质不同的无根树，\sum 可以拆开，于是：

Ans=\sum p-\sum q+\sum s

不难发现两个本质不同的无根树的点等价类是没有交集的，即 p 不会算重，而这种枚举方式一定可以枚举到任意形态的有根树（感性一下），所以其 \mathrm{OGF} 本质就是 A'(x) 。 q 与 s 考虑同理，其 \mathrm{OGF} 分别是 G(x) ，F(x^2) ，于是：

\begin{aligned}A(x)&=P(x)-Q(x)+S(x)\\&=A'(x)-G(x)+F(x^2)\end{aligned}

时间复杂度 O(n\log^2n) 。

Code

贴一个当多项式板子吧，常数巨大。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
namespace Modular{
    const LL mod=998244353;
    inline LL add(const LL a,const LL b){
        return a<mod-b?a+b:a-mod+b;
    }
    inline LL dec(const LL a,const LL b){
        return a<b?a-b+mod:a-b;
    }
    inline LL mul(const LL a,const LL b){
        return a*b%mod;
    }
    inline void Add(LL &a,const LL b){
        a=a<mod-b?a+b:a-mod+b;
    }
    inline void Dec(LL &a,const LL b){
        a=a<b?a-b+mod:a-b;
    }
    inline void Mul(LL &a,const LL b){
        a=a*b%mod;
    }
    inline LL ksm(LL a,LL b){
        LL res=1; while(b) (b&1)&&(Mul(res,a),1),Mul(a,a),b>>=1; return res;
    }
    inline LL Inv(const LL a){
        return ksm(a,mod-2);
    }
};
using namespace Modular;
namespace Polynome{
    struct Poly{
        vector<LL> a;
        Poly(const int n=0,const LL v=0){ a.resize(n+1),a[n]=v; }
        inline LL& operator [] (const int i){ return a[i]; }
        inline const LL operator [] (const int i)const{ return a[i]; }
        inline Poly rsz(const int n){ (*this).a.resize(n+1);return *this; }
        inline int deg()const{ return a.size()-1; }
    };  
    inline void pri(const Poly &A){
        for(int i=0;i<=A.deg();i++) cout<<A[i]<<' '; puts("");
    }
    inline Poly operator + (const Poly &A,const Poly &B){
        Poly C(max(A.deg(),B.deg())); 
        for(int i=0;i<=A.deg();i++) Add(C[i],A[i]);
        for(int i=0;i<=B.deg();i++) Add(C[i],B[i]);
        return C;
    }
    inline Poly operator + (const Poly &A,const LL X){
        Poly C=A; Add(C[0],X); return C;
    }
    inline Poly operator - (const Poly &A,const Poly &B){
        Poly C(max(A.deg(),B.deg())); 
        for(int i=0;i<=A.deg();i++) Add(C[i],A[i]);
        for(int i=0;i<=B.deg();i++) Dec(C[i],B[i]);
        return C;
    }
    inline Poly operator - (const Poly &A,const LL X){
        Poly C=A; Dec(C[0],X); return C;
    }
    int lim; vector<int> pos;
    inline void init(const int n){
        lim=1; while(lim<=n) lim<<=1; pos.resize(lim);
        for(int i=0;i<lim;i++) pos[i]=(pos[i>>1]>>1)|((i&1)*(lim>>1));
    }
    inline void NTT(Poly &A,const int type){
        for(int i=0;i<lim;i++) if(i<pos[i]) swap(A[i],A[pos[i]]);
        for(int b=1;b<lim;b<<=1){ int len=b<<1;
            LL w=ksm(type==1?3:(mod+1)/3,(mod-1)/len);
            for(int k=0;k<lim;k+=len){ LL wk=1;
                for(int i=0;i<b;i++,Mul(wk,w)){
                    LL A0=A[k+i],A1=mul(wk,A[k+i+b]);
                    A[k+i]=add(A0,A1),A[k+i+b]=dec(A0,A1);
                }
            } 
        } 
        if(type==-1){ LL inv=Inv(lim); for(int i=0;i<lim;i++) Mul(A[i],inv); }
    }
    inline Poly operator * (const Poly &A,const Poly &B){
        int n=A.deg(),m=B.deg(); init(n+m); Poly A1(lim-1),B1(lim-1);
        for(int i=0;i<=n;i++) A1[i]=A[i]; for(int i=0;i<=m;i++) B1[i]=B[i]; NTT(A1,1),NTT(B1,1);
        for(int i=0;i<lim;i++) Mul(A1[i],B1[i]); NTT(A1,-1); return A1;
    }
    inline Poly operator * (const Poly &A,const LL X){
        Poly C=A; for(int i=0;i<=C.deg();i++) Mul(C[i],X); return C;
    }   
    inline Poly P_inv(Poly A,const int n){
        if(n==1) return Poly(0,Inv(A[0]));
        Poly B0=P_inv(A.rsz((n+1)>>1),(n+1)>>1);
        return (B0*2-(A*((B0*B0).rsz(n-1))).rsz(n-1)).rsz(n-1);
    }
    inline Poly P_X2(const Poly &A,const int n){
        Poly C(n-1); for(int i=0;i*2<=C.deg();i++) C[i*2]=A[i]; return C;
    } 
    inline Poly P_X3(const Poly &A,const int n){
        Poly C(n-1); for(int i=0;i*3<=C.deg();i++) C[i*3]=A[i]; return C;
    } 
    inline Poly P_X4(const Poly &A,const int n){
        Poly C(n-1); for(int i=0;i*4<=C.deg();i++) C[i*4]=A[i]; return C;
    } 
}
using namespace Polynome;

int n; LL inv2,inv6,inv24; Poly X(1);

inline Poly F(const Poly &A,const int n){
    return A*(-1)+X*((A*A*A+A*P_X2(A,n)*3+P_X3(A,n)*2)*inv6)+1;
}
inline Poly F_div(const Poly &A,const int n){
    return X*((A*A*3+P_X2(A,n)*3)*inv6)-1;
}
inline Poly solve(const int n){
    if(n==1) return Poly(0,1); Poly A0=solve((n+1)>>1);
    return (A0-F(A0,n)*P_inv(F_div(A0,n).rsz(n-1),n).rsz(n-1)).rsz(n-1);
}
inline void write(const LL x){
    if(x>=10) write(x/10); putchar(x%10+48); 
}
#define gc getchar()
#define in read()
inline int read(){
    int f=1,k=0; char cp=gc; while(cp!='-'&&(cp<'0'||cp>'9')) cp=gc; if(cp=='-') f=-1,cp=gc;
    while(cp>='0'&&cp<='9') k=(k<<3)+(k<<1)+cp-48,cp=gc; return f*k;
}
int main(){ inv2=Inv(2),inv6=Inv(6),inv24=Inv(24),X[1]=1;
    n=100000; 
    Poly F=solve(n+1),F2=P_X2(F,n+1),F3=P_X3(F,n+1),F4=P_X4(F,n+1);
    Poly P=1+X*(((F*F).rsz(n)*(F*F).rsz(n)).rsz(n)+((F*F).rsz(n)*F2).rsz(n)*6+(F*F3).rsz(n)*8+(F2*F2).rsz(n)*3+F4*6)*inv24;
    Poly Q=(((F-1)*(F-1)).rsz(n)+F2-1)*inv2;
    Poly S=P_X2(F,n+1);
    Poly A=P-Q+S; 
    for(int i=1;i<=n;i++) write(A[i]),puts("");

    return 0;
}