无标号荒漠计数

比较综合的一道图计数题目。

本文中，以大写字母表示对应小写字母数列的普通生成函数（OGF）。

前置知识

多项式全家桶（乘法，求逆，ln，exp）

欧拉变换

\operatorname{Euler}[F(x)]=\prod_{i\geq 1}\left(\frac{1}{1-x^i}\right)^{f_i}=\exp\sum_{i\geq 1}\frac{F(x^i)}{i}

可以这样理解其组合意义：假设有很多种数量无穷多的物品，每种物品有一个大小，大小为 i 的物品共 f_i 种，求选出的物品大小之和为 n 的方案数。实际上就是一个背包。

Burnside 引理 / Polya 定理

l=\frac{1}{|G|}\sum_{g\in G}m^{\sigma_g}

即对于置换群 G，等价类的个数等于每种置换下不动点个数的平均值。而对于染色问题，因为一个环中染色必须相同，所以置换 g 下的不动点个数为 m^{\sigma_g}，其中 m 为染色方案，\sigma_g 表示置换 g 的环个数。

无标号圆根圆方树（有根仙人掌）计数

与无标号树计数相同，我们先考虑有根的情况，稍后再转化为无根。

我们知道每个仙人掌对应着一棵圆方树。而圆方树与普通树有一些共性，也有一些不同的地方，比如：

圆点和方点本质不同，且只有圆点占据体积。
圆点与圆点，方点与方点不能相邻。
方点的度数至少为 2。
圆点的儿子是无序的，而方点的儿子构成的是一个环，因此是有序的。但是对称同构将视为同一种情况，因为环可以翻转。

由后两条性质可以看出圆点为根与方点为根是不一样的。因为圆点为根时不需要考虑环的问题，所以我们考虑先求出这种情况的答案。

设 f_n 表示 n 个圆点，以圆点为根的圆方树个数；再设 g_n 表示 n 个圆点，以方点为根的有父亲圆方树个数。注意这里的 g_n 是与普通方根圆方树不同的，因为父亲的方向就是根的方向，它是「独特的」，所以只需考虑对称同构（即翻转）而无需考虑循环同构。为了方便，我们钦定 f_0=g_0=0。

为了求出 g_n，我们定义 h_n 表示将 g_n 去掉对称同构限制的结果。发现 h 数组就是对 f 数组的一个背包，即

H(x)=F(x)+F(x)^2+\cdots=\frac{F(x)}{1-F(x)}

容易看出 h_n 的常数项为 0，也就是说满足了方点至少有一棵子树的性质，这也就是之前钦定 f_0=0 的目的。

而 g 数组考虑了对称同构，我们知道对称同构等价类数等于所有方案数与对称方案数之和除以二，即

g_n=\frac{1}{2}\left(h_n+\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\max(h_i,1)\max(f_{n-2i},1)\right)

前面的 h_n 是总方案数，这不难理解；而后面枚举 i 表示枚举两侧对称部分的大小，f_{n-2i} 是中心子树的方案数。取 \max 是为了考虑对称部分为空或中心子树为空的情况，即把 f_0 与 h_0 补成 1。写成生成函数运算，得

G(x)=\dfrac{H(x)+[H(x^2)+1][F(x)+1]}{2}

（我第一次推的时候没有把 f_0 补成 1，导致后面的推导全部有误，只能重写 /dk）

那么圆点为根的情况与普通树计数相同，是一堆有父亲方根子树的无序排列，即

F(x)=x\operatorname{Euler}[G(x)]=x\exp\sum_{i\geq 1}\frac{G(x^i)}{i}

而 G(x) 可以展开成由 F(x) 表示的形式，因此上式是一个与 F(x) 有关的方程，牛顿迭代即可求解。

然而牛顿迭代里面套 exp 的常数比较大，可以考虑对两边同时取 ln，得

\ln \frac{F(x)}{x}-\sum_{i\geq 1}\frac{G(x^i)}{i}=0

对这个式子求导时注意到 i\geq 2 时 \dfrac{G(x^i)}{i} 只与 F(x^i) 有关，可以由上一层迭代直接得到，因此可将其视为常数。这样我们能在 O(n\log n) 的时间内解决这个问题了。

无标号方根圆方树计数

方点做根的情况还要难办一些。前面已经提到了有父亲和无父亲情况的不同，因为此时并没有一个独特的子树，所以它的所有圆根子树构成的环既可以旋转也可以翻转。这是一个无标号环同构问题，首先要考虑的就是 Burnside 引理。

假设根节点有 k 个圆根子树，那么置换一共有 2k 种，其中 k 种是旋转 i\in[0,k-1] 次，另外 k 种是旋转了 i 次以后再翻转。枚举 k,i，可以写出生成函数式子：

\begin{aligned} 2F^\square(x)=&\sum_{k\geq 2}\frac{1}{k}\sum_{i=0}^{k-1}\left(F\left(x^{k/\gcd(k,i)}\right)^{\gcd(k,i)}+\left\{\begin{matrix}F(x)F(x^2)^\frac{k-1}{2}&2\nmid k\\F(x^2)^{k/2}&2\mid k,2\nmid i\\F(x)^2F(x^2)^{k/2-1}&2\mid k,2\mid i\end{matrix}\right.\right)\\ =&\sum_{k\geq 2}\frac{1}{k}\sum_{d\mid k}\varphi(d)F(x^d)^{k/d}+\sum_{i\geq 1}F(x)F(x^2)^i+\frac{1}{2}\sum_{i\geq 1}F(x^2)^i+\frac{1}{2}\sum_{i\geq 0}F(x)^2F(x^2)^i\\ =&\sum_{d\geq 1}\frac{\varphi(d)}{d}\sum_{i\geq 1}\frac{F(x^d)^i}{i}+\left(F(x)+\frac{F(x)^2+1}{2}\right)\sum_{i\geq 0}F(x^2)^i-F(x)+c\\ =&-\sum_{d\geq 1}\frac{\varphi(d)}{d}\ln [1-F(x^d)]+\dfrac{[F(x)+1]^2}{2-2F(x^2)}-F(x)+c \end{aligned}

如果看不太明白也没有关系，我们来分析一下这个式子。

第一行：为了方便我将 \dfrac{1}{2} 移项，所以后面的 \dfrac{1}{2k} 都变成了 \dfrac{1}{k}。旋转 i 次的置换会产生 \gcd(k,i) 个环，每个环的大小为 \dfrac{k}{\gcd(k,i)}，环里面的元素需要相同，就相当于占了多倍的空间。而旋转 i 次再翻转的置换会产生大小为 1 或 2 的一些环，可以按 k,i 的奇偶性讨论分为三类。

第二行：枚举最大公因数为 k/d，根据一些基本的数论套路，可以得出第一项和式。后面的三项分别为三种情况的贡献之和。

第三行：前面交换和式，后面可以提出一个公比为 F(x^2) 的等比数列来。注意，因为改变了某些求和下界，这里会多出一个 -F(x)+c 的余项出来，其中 c 是常数项，不过并没有什么意义，可以忽略。

第四行：我们知道 -\ln(1-x) 是幂级数 \sum\limits_{i\geq 1}\dfrac{x^i}{i} 的封闭形式，代入原式可得前半部分。后半部分则利用了等比数列求和公式。

因为 \ln[1-F(x^d)] 只与 x^d 有关，所以可以只计算出 \ln[1-F(x)] 然后扩倍即可。这部分的时间复杂度也是 O(n\log n)。

无标号无根仙人掌计数

有了有根情况的答案，我们接下来仍可以参考无标号树计数的套路，选出一个具有重心性质的点来代表这个仙人掌，并用总方案数减去不合法方案。那么这个特殊的点是哪一个呢？

关注圆方树的一个性质：对于任意一棵圆方树，有且仅有一个节点，使得如果它是圆点则所有子树中圆点个数都 \leq \left\lfloor\dfrac{n-1}{2}\right\rfloor，如果它是方点则所有子树中圆点个数都 \leq \left\lfloor\dfrac{n}{2}\right\rfloor。

（这个性质是我玩样例的时候乱猜的，实际上是对的，但我现在还不会证 qwq）

那么我们枚举每个节点，用所有以它为根的方案减去它不是重心的方案，得

f^*_n=f_n-\sum_{i=\lfloor\frac{n+1}{2}\rfloor}^{n-1}g_if_{n-i}+f^\square_n-\sum_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}f_ig_{n-i} F^*(x)=F(x)+F^\square(x)-F(x)G(x)

两部分巧妙的拼合在了一起，变成了只需要一个卷积就能处理的问题。所以复杂度也是 O(n\log n)。

一个无标号荒漠由若干个无标号仙人掌无序组合构成，也就是对 F^*(x) 做欧拉变换。需要做一次多项式 exp，复杂度依然是 O(n\log n)，与前三部分相同，因此总时间复杂度为 O(n\log n)。

下面是我的代码，但看懂上面的推理后就不建议看代码了，式子写得很乱。

#include <cstdio>
#include <algorithm>
using namespace std;
const int MAXN=1<<19;
const int P=998244353;
const int G=3, IG=332748118;
const int I2=499122177;
int qpow(int n, int k)
{
    int r=1;
    while (k)
    {
        if (k&1) r=1ll*r*n%P;
        n=1ll*n*n%P, k>>=1;
    }
    return r;
}
/* poly templates begin */
int l, r[MAXN];
void getl(int);
void NTT(int*, int);
void poly_inv(int*, int*, int);
void poly_diff(int*, int*, int);
void poly_int(int*, int*, int);
void poly_ln(int*, int*, int);
void poly_exp(int*, int*, int);
/* poly templates end */
int f[MAXN], g[MAXN], h[MAXN], t1[MAXN], t2[MAXN];
int a[MAXN], b[MAXN];
void solve(int n)
{
    if (n==2) return a[0]=0, a[1]=1, void();
    solve(n/2+1);
    for (int i=0; i<n; i++) t1[i]=a[i+1];
    poly_ln(t1, h, n);
    for (int i=0; i<n; i++)
        f[i]=a[i], g[i]=i&1?0:a[i>>1];
    getl(5*n);
    for (int i=n; i<l; i++) f[i]=g[i]=0;
    NTT(f, 1), NTT(g, 1);
    for (int i=0; i<l; i++)
        t1[i]=2ll*(1-f[i]+P)*(g[i]-1+P)%P;
    NTT(t1, -1);
    poly_inv(t1, t2, n);
    getl(5*n);
    for (int i=n; i<l; i++) h[i]=t2[i]=0;
    NTT(t2, 1);
    for (int i=0; i<l; i++)
        b[i]=(((1ll*f[i]+2*g[i]-2)%P*f[i]-g[i])%P*t2[i]%P+P)%P;
    NTT(b, -1);
    for (int i=1; i<n; i++)
    {
        int ii=(P-qpow(i, P-2))%P;
        for (int j=1; i*j<n; j++)
            h[i*j]=(h[i*j]+1ll*b[j]*ii)%P;
    }
    NTT(h, 1);
    for (int i=0; i<l; i++)
        t1[i]=((((1ll*f[i]+2*g[i]-4)%P*f[i]-5ll*g[i]+6)%P*f[i]+2*g[i]-2)%P+P)%P;
    NTT(t1, -1);
    poly_inv(t1, t2, n);
    getl(5*n);
    for (int i=n; i<l; i++) t2[i]=0;
    NTT(t2, 1);
    for (int i=0; i<l; i++)
        f[i]=2ll*f[i]%P*(f[i]-1)%P*(f[i]-1)%P*(g[i]-1+P)%P*h[i]%P*t2[i]%P;
    NTT(f, -1);
    for (int i=0; i<n; i++) a[i]=(a[i]-f[i]+P)%P;
}
int pri[MAXN], phi[MAXN];
bool mark[MAXN];
int main()
{
//  freopen("cactus.in", "r", stdin);
//  freopen("cactus.out", "w", stdout);
    int n;
    scanf("%d", &n);
    phi[1]=1;
    for (int i=2, cnt=0; i<=n; i++)
    {
        if (!mark[i]) pri[++cnt]=i, phi[i]=i-1;
        for (int j=1; j<=cnt&&i*pri[j]<=n; j++)
        {
            mark[i*pri[j]]=1;
            if (i%pri[j]==0)
            {
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    }
    solve(++n);
    for (int i=0; i<n; i++)
        f[i]=a[i], g[i]=i&1?0:a[i>>1];
    getl(3*n);
    for (int i=n; i<l; i++) f[i]=g[i]=0;
    NTT(f, 1), NTT(g, 1);
    for (int i=0; i<l; i++)
        t1[i]=2ll*(1-f[i]+P)*(g[i]-1+P)%P;
    NTT(t1, -1);
    poly_inv(t1, t2, n);
    getl(3*n);
    for (int i=n; i<l; i++) h[i]=t2[i]=0;
    NTT(t2, 1);
    for (int i=0; i<l; i++)
        b[i]=(((1ll*f[i]+2*g[i]-2)%P*f[i]-g[i])%P*t2[i]%P+P)%P;
    NTT(b, -1);
    for (int i=n; i<l; i++) b[i]=0;
    for (int i=0; i<n; i++)
        t1[i]=i&1?0:2*(P-a[i>>1])%P;
    t1[0]=(t1[0]+2)%P;
    poly_inv(t1, t2, n);
    getl(3*n);
    for (int i=0; i<n; i++) t1[i]=a[i];
    for (int i=n; i<l; i++) t1[i]=t2[i]=0;
    NTT(t1, 1), NTT(t2, 1);
    for (int i=0; i<l; i++)
        f[i]=1ll*(t1[i]+1)*(t1[i]+1)%P*t2[i]%P;
    NTT(f, -1);
    f[0]=0;
    for (int i=0; i<n; i++) t1[i]=(P-a[i])%P;
    t1[0]=(t1[0]+1)%P;
    poly_ln(t1, t2, n);
    for (int i=1; i<n; i++)
    {
        int ii=1ll*(P-phi[i])*qpow(i, P-2)%P;
        for (int j=1; i*j<n; j++)
            f[i*j]=(f[i*j]+1ll*t2[j]*ii)%P;
    }
    getl(2*n);
    for (int i=n; i<l; i++) b[i]=0;
    NTT(a, 1), NTT(b, 1);
    for (int i=0; i<l; i++) a[i]=1ll*a[i]*b[i]%P;
    NTT(a, -1);
    for (int i=0; i<n; i++) a[i]=(1ll*f[i]*I2-a[i]+P)%P;
    for (int i=0; i<n; i++) b[i]=0;
    for (int i=1; i<n; i++)
    {
        int ii=qpow(i, P-2);
        for (int j=1; i*j<n; j++)
            b[i*j]=(b[i*j]+1ll*a[j]*ii)%P;
    }
    poly_exp(b, a, n+1);
    for (int i=1; i<n; i++) printf("%d\n", a[i]);
    return 0;
}

冷滟泽的个人博客