偶然发现oeis首页的数列是点数为 nn 的无标号无根树的数量！于是心血来潮地去学了。

无标号有根树计数

之前做过烷基计数，而它是对节点度数有限制的树的计数，直接套用burnside引理即可。
而这里我们大概也要考虑使用类似的思路。仍然考虑burnside。
而这是遇到的麻烦就是置换群太庞大了。不可能枚举的。那么考虑置换是由好多轮换生成的，用生成函数代替枚举置换群的过程即可。
这里先写出最后的式子：

A(x)=x⋅exp(∑i>0A(xi)i)A(x)=x\cdot exp(\sum_{i>0} \frac{A(x^i)}{i})

因为我们使用了burnside引理，这样使得无标号计数转化为有标号计数，使得我们可以使用指数生成函数。
exp里面的东西就是轮换的指数生成函数，而把 xix^i 换成了 A(xi)A(x^i)。因为同一个轮换中的元素是被钦定相同的。这样再用 expexp 组合一下就变成了置换。因为是指数生成函数，所以最后选择了 kk 个轮换的部分需要乘 k!k!。而正好burnside引理又需要除掉 k!k!，所以就变成了这个样子啦！

（其中 exp(∑i>0A(xi)i)=MESTA=∏i(1−xi)−aiexp(\sum_{i>0} \frac{A(x^i)}{i})=MEST A=\prod_{i}(1-x^i)^{-a_i}，被称作 Polya指数 here，它在无标号题目中扮演了有标号中exp的位置，在链接中的题目可以体现）

不过看起来它仍然不能快速计算啊其实是可以快速计算的...还是考虑倍增地求 AA，那么在一次从膜 xnx^n 到膜 x2nx^{2n} 的扩展中，A(xi),i>1A(x^i),i>1 都是已知的。然后发现这部分有用的项是调和级数的，把与它们有关的部分算出之后，就成了一个与 A(x)A(x) 有关的方程，用牛顿迭代的方法扩展即可。复杂度仍然是 T(n)=(T/2)+O(nlog⁡n)=O(nlog⁡n)T(n)=(T/2)+O(n\log n)=O(n\log n)，不过常数大得离谱就是了。不过EI的做法用了解微分方程，我也不知道为什么。。。它微分方程的式子是原式两边除 BB 再取 ln⁡\ln 再求导得到的。这样就不用处理exp套起来的东西了...

于是我们换一个角度。如果不从polya定理考虑，我们会发现没有度数限制反而方便了我们。

发现一个点的所有儿子的子树相当于是把所有类型的子树做了一个完全背包！即

A(x)=x∏i(∑jxij)aiA(x)=x\prod_i(\sum_jx^{ij})^{a_i}

A(x)=x∏i(11−xi)aiA(x)=x\prod_i(\frac{1}{1-x^i})^{a_i}

ln⁡A(x)=ln⁡x−∑iailn⁡(1−xi)\ln A(x)=\ln x-\sum_i a_i\ln(1-x^i)

因为我们不知道 aia_i，所以貌似不能像平常完全背包一样去拆 ln⁡\ln 的展开再贡献，最后exp。而且还有一个碍事的 ln⁡x\ln x 呢。所以继续推。两边求导。

A′(x)A(x)=1x+∑iaiixi−11−xi\frac{A'(x)}{A(x)}=\frac{1}{x}+\sum_i a_ii\frac{x^{i-1}}{1-x^i}

xA′(x)=A(x)+A(x)∑iaiixi1−xixA'(x)=A(x)+A(x)\sum_i a_ii\frac{x^i}{1-x^i}

xA′(x)=A(x)+A(x)∑iaii(∑j>0xij)xA'(x)=A(x)+A(x)\sum_i a_ii(\sum_{j>0}x^{ij})

然后比较第 nn 项系数

nan=an+∑iai∑jajj[j∣n−i]na_n=a_n+\sum_i a_i\sum_j a_jj[j|n-i]

an=∑iai∑j∣n−iajjn−1a_n=\frac{\sum_i a_i\sum_{j|n-i} a_jj}{n-1}

然后根据 A0=0,A1=1A_0=0,A_1=1 就可以递推了。可以设 bi=∑j∣iajjb_i=\sum_{j|i}a_jj，然后大概就是一个 AA 和 BB 卷积的形式。算出每一个 aa 的时候往它倍数上贡献一下即可算出 bb。
这样暴力是 O(n2)O(n^2)。而它计算是一个卷积的形式，不过很难直接写成生成函数方程的形式，于是不能牛顿迭代，可以用分治NTT做到 O(nlog⁡2n)O(n\log^2n)。类似 LOJ575 不等关系。这个推导也告诉我们完全背包和01背包有一种不用exp的分治NTT做法。

还有一个比较有意思的事情就是 lim⁡n→∞an+1an=2.955765...\lim_{n\to \infty}\frac{a_{n+1}}{a_n}=2.955765... 不太知道这个是咋算的hhh。渐进表达式是 an=c∗dn∗n−1.5a_n=c*d^n*n^{-1.5},c =0.439924...,d =2.955765...。

无标号无根树计数

先给出式子吧。设无根树的生成函数为 T(x)T(x)，有根树生成函数为 A(x)A(x)，则

T(x)=1+A(x)−A2(x)2+A(x2)2T(x)=1+A(x)-\frac{A^2(x)}{2}+\frac{A(x^2)}{2}

这个式子显然算出有根树之后再卷积一次就可以计算。
关于这个式子的理解有两种方式。

第一种也是比较简单的一种：关于重心。
因为一棵树只会有一到两个重心，所以我们考虑用重心去计数。用有根树的数量减去当前根不是重心的数量即可。而当前根不是重心则有且仅有一个大小超过一半的子树，于是

tn=an−∑i>n/2aian−it_n=a_n-\sum_{i>n/2}a_ia_{n-i}

然后把这个上三角样子的卷积乘 22 得到一个正常的卷积。但是要单独考虑对角线上的元素。这里对角线上的元素就是枚举的这两棵树一样。所以出现了 A(x2)A(x^2)。然后细节推一推就会发现就是上面这个式子。

第二种：构造关于等价类的计数。
这篇关于烷烃计数的题解讲得很明白。而烷烃计数也可以用上面的这个重心的思路来理解。

貌似撒花了？

#include<bits/stdc++.h>
using namespace std;
const int N=4e5+50,mod=998244353;
int n,r[N],rt[N],irt[N],lo[N],inv[N],ans=1;
int Glim(int n){int z=1;while(z<n)z<<=1;return z;}
void getr(int n){for(int i=1;i<n;i++)r[i]=r[i>>1]>>1^(i&1)<<lo[n];}
int power(int x,int y){
    int z=1;
    for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
    return z;
}
void init(int n){
    int lim=Glim(n);inv[1]=1;
    for(int i=2;i<=lim;i++)inv[i]=mod-1ll*mod/i*inv[mod%i]%mod;
    for(int i=2,j=0;i<=lim;i<<=1,j++)lo[i]=j;
    for(int i=1;i<n;i<<=1){
        int Wn=power(3,(mod-1)/(i<<1));
        for(int k=0,w=1;k<i;k++,w=1ll*w*Wn%mod)rt[i+k]=w;
    }
    for(int i=1;i<lim;i++)irt[i]=power(rt[i],mod-2);
}
void NTT(int *a,int n,int op){
    for(int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
    for(int i=1,*w=op>0?rt:irt;i<n;i<<=1)
        for(int j=0;j<n;j+=i<<1)
            for(int k=0;k<i;k++){
                int x=a[j+k],y=1ll*a[i+j+k]*w[i+k]%mod;
                a[j+k]=(x+y)%mod;a[i+j+k]=(x-y)%mod;
            }
    if(op<0)for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv[n]%mod;
}
int a[N],b[N],A[N],B[N],C[N],D[N],E[N],F[N],G[N];
void Inv(int *a,int *b,int n){
    if(n==1){b[0]=power(a[0],mod-2);return;}
    int m=(n+1)>>1,lim=Glim(n+m);Inv(a,b,m);
    for(int i=0;i<lim;i++)E[i]=i<n?a[i]:0,F[i]=i<m?b[i]:0;
    getr(lim);NTT(E,lim,1);NTT(F,lim,1);
    for(int i=0;i<lim;i++)F[i]=1ll*E[i]*F[i]%mod*F[i]%mod;
    NTT(F,lim,-1);for(int i=m;i<n;i++)b[i]=-F[i];
}
void dao(int *a,int *b,int n){for(int i=0;i<n-1;i++)b[i]=1ll*a[i+1]*(i+1)%mod;b[n-1]=0;}
void ji(int *a,int *b,int n){for(int i=0;i<n-1;i++)b[i+1]=1ll*a[i]*inv[i+1]%mod;b[0]=0;}
void ln(int *a,int *b,int n){
    dao(a,A,n);Inv(a,B,n);int lim=Glim(n*2);
    for(int i=n;i<lim;i++)A[i]=B[i]=0;
    getr(lim);NTT(A,lim,1);NTT(B,lim,1);
    for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
    NTT(A,lim,-1);ji(A,b,n);
}
void exp(int *a,int *b,int n){
    if(n==1){b[0]=1;return;}
    int m=(n+1)>>1,lim=Glim(n);exp(a,b,m);
    for(int i=0;i<lim;i++)C[i]=i<m?b[i]:0;
    ln(C,G,n);for(int i=0;i<lim;i++)G[i]=i<n?a[i]-G[i]:0;
    getr(lim);NTT(C,lim,1);NTT(G,lim,1);
    for(int i=0;i<lim;i++)C[i]=1ll*C[i]*G[i]%mod;
    NTT(C,lim,-1);for(int i=m;i<n;i++)b[i]=C[i];
}
int P[N],Q[N];
void solve1(int *a,int n){//nlogn
    if(n==1){a[0]=0;return;}
    int m=n+1>>1,lim=Glim(n*2);
    solve1(a,m);for(int i=m;i<n;i++)a[i]=0;
    for(int i=0;i<n;i++)P[i]=0;
    for(int i=2;i<n;i++)
        for(int j=1;i*j<n;j++)
            P[i*j]=(P[i*j]+1ll*a[j]*inv[i])%mod;
    exp(P,b,n);exp(a,P,n);
    for(int i=0;i<lim;i++)Q[i]=i<n?b[i]:P[i]=0;
    getr(lim);NTT(P,lim,1);NTT(Q,lim,1);
    for(int i=0;i<lim;i++)P[i]=1ll*P[i]*Q[i]%mod;
    NTT(P,lim,-1);
    for(int i=n-1;~i;i--)P[i]=i?P[i-1]:-1;
    Inv(P,Q,n);P[0]=0;
    for(int i=0;i<n;i++)(P[i]-=a[i])%=mod;
    for(int i=n;i<lim;i++)P[i]=Q[i]=0;
    getr(lim);NTT(P,lim,1);NTT(Q,lim,1);
    for(int i=0;i<lim;i++)P[i]=1ll*P[i]*Q[i]%mod;
    NTT(P,lim,-1);
    for(int i=m;i<n;i++)a[i]=-P[i];
}
void solve2(int l,int r){//nlog^2n
	if(l==r){if(l^1)a[l]=1ll*a[l]*inv[l-1]%mod;for(int i=l;i<=n;i+=l)b[i]=(b[i]+1ll*a[l]*l)%mod;return;}
	int mid=(l+r)>>1;solve2(l,mid);
	if(l==1){
		int lim=Glim(2*mid+1);getr(lim);
		for(int i=0;i<lim;i++)A[i]=i<=mid?a[i]:0,B[i]=i<=mid?b[i]:0;
		NTT(A,lim,1);NTT(B,lim,1);
		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
		NTT(A,lim,-1);
		for(int i=mid+1;i<=r;i++)(a[i]+=A[i])%=mod;
	}
	else{
		int lim=Glim(r+mid-2*l+1);getr(lim);
		for(int i=0;i<lim;i++)A[i]=i+l<=mid?a[i+l]:0,B[i]=i+l<=r?b[i]:0;
		NTT(A,lim,1);NTT(B,lim,1);
		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
		NTT(A,lim,-1);
		for(int i=mid+1;i<=r;i++)(a[i]+=A[i-l])%=mod;
		for(int i=0;i<lim;i++)A[i]=i+l<=mid?b[i+l]:0,B[i]=i+l<=r?a[i]:0;
		NTT(A,lim,1);NTT(B,lim,1);
		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
		NTT(A,lim,-1);
		for(int i=mid+1;i<=r;i++)(a[i]+=A[i-l])%=mod;
	}
	solve2(mid+1,r);
}
double c=0.439924,d=2.955765;
int main(){
	scanf("%d",&n);init((n+1)*2);
	solve1(a,n+1);/* a[1]=1;solve2(1,n); */
    b[0]=1;for(int i=1;i<=n;i++)b[i]=a[i];
	for(int i=1;i<=n;i++)if(i%2==0)b[i]=(b[i]+1ll*a[i/2]*inv[2])%mod;
	int lim=Glim(n*2+1);getr(lim);
	for(int i=0;i<lim;i++)A[i]=i<=n?a[i]:0;
	NTT(A,lim,1);
	for(int i=0;i<lim;i++)A[i]=1ll*A[i]*A[i]%mod;
	NTT(A,lim,-1);
	for(int i=1;i<=n;i++)b[i]=(b[i]-1ll*A[i]*inv[2])%mod;
	for(int i=1;i<=n;i++)(a[i]+=mod)%=mod,(b[i]+=mod)%=mod;
	for(int i=0;i<=n;i++)printf("%d %.4lf\n",a[i],c*pow(d,i)*pow(i,-1.5));//后面这个是渐进表达式，可以看出它非常准确。
	for(int i=0;i<=n;i++)printf("%d ",b[i]);
    return 0;
}

10510^5 的数据 O(nlog⁡n)O(n\log n) 跑了1326ms，O(nlog⁡2n)O(n\log^2n)跑了650ms。单log又一次被双log爆踩。

改了vector的板子之后写了一遍，vector真香啊，是真香，不过超级吃O2...

#include<bits/stdc++.h>
#define vec vector<z>
using namespace std;
const int mod=998244353,N=4e5+50;
int r[N],lo[N];
struct z{
    int x;
    z(int x=0):x(x){}
    friend z operator +(z a,z b){return (a.x+=b.x)>=mod?a.x-mod:a.x;}
    friend z operator -(z a,z b){return (a.x-=b.x)<0?a.x+mod:a.x;}
    friend z operator *(z a,z b){return 1ll*a.x*b.x%mod;}
}w[N],inv[N];
z power(z x,int y){z t=1;for(;y;y>>=1,x=x*x)if(y&1)t=t*x;return t;}
int Glim(int n){int t=1;while(t<n)t<<=1;return t;}
void init(int n){
    int lim=Glim(n);inv[1]=1;
    for(int i=2;i<=lim;i++)inv[i]=mod-mod/i*inv[mod%i];
    for(int i=2,j=0;i<=lim;i<<=1,j++)lo[i]=j;
    for(int i=1;i<lim;i<<=1){
        z Wn=power(3,(mod-1)/(i<<1)),W=1;
        for(int j=0;j<i;j++,W=W*Wn)w[i+j]=W;
    }
}
void dft(vec &a,int lim){
    a.resize(lim);
    for(int i=0;i<lim;i++)r[i]=r[i>>1]>>1^(i&1)<<lo[lim];
    for(int i=0;i<lim;i++)if(i>r[i])swap(a[i],a[r[i]]);
    for(int i=1;i<lim;i<<=1)
        for(int j=0;j<lim;j+=i<<1)
            for(int k=0;k<i;k++){
                z x=a[j+k],y=a[i+j+k]*w[i+k];
                a[j+k]=x+y;a[i+j+k]=x-y;
            }
}
void idft(vec &a,int lim){
    dft(a,lim);reverse(&a[1],&a[lim]);
    for(int i=0;i<lim;i++)a[i]=a[i]*inv[lim];
}
vec operator *(vec a,vec b){
    int n=a.size()+b.size()-1,lim=Glim(n);
    dft(a,lim);dft(b,lim);
    for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
    idft(a,lim);a.resize(n);return a;
}
vec Inv(const vec &a,int n=0){//注意这里n不能超过a.size()
    if((n=n?n:a.size())==1)return vec(1,power(a[0],mod-2));
    int m=n+1>>1,lim=Glim(n+m);
    vec b=Inv(a,m),d=b,c(&a[0],&a[n]);
    dft(b,lim);dft(c,lim);
    for(int i=0;i<lim;i++)b[i]=mod-b[i]*b[i]*c[i];
    idft(b,lim);for(int i=0;i<m;i++)b[i]=d[i];
    b.resize(n);return b;
}
vec dao(vec a){for(int i=0;i<a.size()-1;i++)a[i]=(i+1)*a[i+1];*a.rbegin()=0;return a;}
vec ji(vec a){for(int i=a.size()-1;i;i--)a[i]=inv[i]*a[i-1];a[0]=0;return a;}
vec ln(const vec &a){vec b=dao(a)*Inv(a);b.resize(a.size());return ji(b);}
vec exp(const vec &a,int n=0){
    if((n=n?n:a.size())==1)return vec(1,1);
    int m=n+1>>1,lim=Glim(n);
    vec b=exp(a,m),d=b;b.resize(n);vec e=ln(b);
    for(int i=0;i<n;i++)e[i]=a[i]-e[i];
    dft(b,lim);dft(e,lim);
    for(int i=0;i<lim;i++)b[i]=b[i]*e[i];
    idft(b,lim);for(int i=0;i<m;i++)b[i]=d[i];
    b.resize(n);return b;
}
vec solve1(int n){
    if(n==1)return vec(1,0);
    int m=n+1>>1,lim=Glim(n);
    vec a=solve1(m),d=a,b=vec(n,0);
    for(int i=2;i<n;i++)for(int j=1;i*j<n;j++)b[i*j]=b[i*j]+a[j]*inv[i];
    a.resize(n);b=exp(b)*exp(a);b.resize(n);a.resize(m);
    for(int i=n-1;~i;i--)b[i]=i?b[i-1]:mod-1;a[0]=mod-1;
    b=Inv(b);dft(a,lim);dft(b,lim);
    for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
    idft(a,lim);for(int i=0;i<m;i++)a[i]=d[i];
    a.resize(n);return a;
}
vec a,b,c;int n;
void solve2(int l,int r){
    if(l==r){a[l]=l^1?a[l]*inv[l-1]:1;for(int i=l;i<n;i+=l)b[i]=b[i]+a[l]*l;return;}
    int mid=(l+r)>>1;solve2(l,mid);
    if(l==1){
        vec c=vec(&a[0],&a[mid+1])*vec(&b[0],&b[mid+1]);
        for(int i=mid+1;i<=r;i++)a[i]=a[i]+c[i];
    }
    else{
        vec c=vec(&a[0],&a[r-l+1])*vec(&b[l],&b[mid+1]);//一定注意是左闭右开（这样写分治fft也太简单了吧QAQ）
        vec d=vec(&b[0],&b[r-l+1])*vec(&a[l],&a[mid+1]);
        for(int i=mid+1;i<=r;i++)a[i]=a[i]+c[i-l]+d[i-l];
    }
    solve2(mid+1,r);
}
int main(){
    scanf("%d",&n);init(2*n);
    a=solve1(n);
    // a.resize(n);b.resize(n);solve2(1,n-1);
    c=a*a;c.resize(n);
    for(int i=0;i<n;i++)c[i]=!i+a[i]+((i&1?0:a[i/2])-c[i])*inv[2];
    for(int i=0;i<n;i++)printf("%d ",c[i]);
    return 0;
}

开O2这甚至跑得比数组快，不开O2慢十倍，10510^5 要跑10秒。。。。而数组版不开O2只会慢两倍多

有标号有向无环图

这个好像跑题了，不过没找到合适的地方放它233
我们设 g(n,i)g(n,i) 表示 nn 个点，至少 ii 个点入度为 00 的方案数，f(n,i)f(n,i) 表示 nn 个点，恰好 ii 个点入度为 00 的方案数。这两个数组满足二项式反演关系。注意到 g(n,0)g(n,0) 就是 nn 个点时的答案。（而不是 g(n,1)g(n,1)，因为 g(n,1)=∑ii⋅f(n,i)g(n,1)=\sum_ii\cdot f(n,i) ）而我们有这样的式子：g(n,i)=(ni)2i(n−i)g(n−i,0)g(n,i)={n\choose i}2^{i(n-i)}g(n-i,0)，然后推一下式子发现可以递推的。

g(n,0)=∑if(n,i)g(n,0)=\sum_if(n,i)

=∑i∑j(−1)j−i(ji)g(n,j)=\sum_i\sum_j(-1)^{j-i}{j\choose i}g(n,j)

=∑i∑j(−1)j−i(ji)(nj)2j(n−j)g(n−j,0)=\sum_i\sum_j(-1)^{j-i}{j\choose i}{n\choose j}2^{j(n-j)}g(n-j,0)

=∑j(nj)2j(n−j)g(n−j,0)∑i(−1)j−i(ji)=\sum_j{n\choose j}2^{j(n-j)}g(n-j,0)\sum_i(-1)^{j-i}{j\choose i}

=∑j(nj)2j(n−j)g(n−j,0)[j=0]=\sum_j{n\choose j}2^{j(n-j)}g(n-j,0)[j=0]

推了个寂寞...要加上 i>0i>0 的限制。然后推得

g(n,0)=∑i>0(ni)(−1)i−12i(n−i)g(n−i,0)g(n,0)=\sum_{i>0}{n\choose i}(-1)^{i-1}2^{i(n-i)}g(n-i,0)

于是可以 n2n^2 递推。然后发现这是卷积，可以分治FFT做到 nlog⁡2nn\log^2n...

7.30upd。。
上面的这个推导太丑了。。。贴上一个2019冬令营课件上的推导，非常清晰。

这个DAG计数还是可以优化的。发现这个式子好像很难看，也不是OGF也不是EGF。实际上我们有这样的trick：ij=(i+j2)−(i2)−(j2)ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}，这个式子只需要考虑组合意义即可。这个东西就是Bluestein算法的那个trick。这个trick经常用来把指数上的乘积项拆掉，还是很好用的。
值得一提的是这个东西还有另一种转化，即 ij=12((i+j)2−i2−j2)ij=\frac{1}{2}((i+j)^2-i^2-j^2)，这个就是完全平方公式。然而有时会出现没有逆元的问题。组合数优秀的地方在于它不需要求逆元。

根据这个，我们可以重新定义一种生成函数：F(x)=∑i≥0fixii!2(i2)F(x)=\sum\limits_{i\geq 0}\frac{f_ix^i}{i!2^\binom{i}{2}}。冬令营课件上把它叫作组合生成函数SGF...（貌似只在有标号有向图计数上有应用，因为 2j(i−j)2^{j(i-j)} 这个形式在有向图中经常出现）
然后设 G(x)G(x) 是 gig_i 的组合生成函数，H(x)H(x) 是 (−1)i−1(-1)^{i-1} 的组合生成函数，则递推式可以写成 G(x)=G(x)H(x)+1G(x)=G(x)H(x)+1。然后只需求逆即可。

关于强连通图有向的计数课件中也有提到，与DAG计数的思想比较相似（主要是一般有向图考虑缩点后变成DAG，然后建立关系），也可以做到 O(nlog⁡n)O(n\log n)。

限制点双形态的无向图计数

这个是看了rqy课件...

就是说我们有一个集合 SS，要求合法图的所有点双都在集合 SS 里（同构意义的），计合法图的数量。
注意到树计数，仙人掌计都可以规约到这个问题：树的集合 SS 是一条边，仙人掌的集合 SS 是一条边和所有环。（这里点双集合改成边双集合大概也可以，主要是它得是图中可分割可表示的元素）

还可以扩展：定义仙人图为每条边最多在两个环上的图，求仙人图个数。这个相当于是在仙人掌的基础上给集合 SS 加了恰好两个点度数为 33，即两个环有公共链的情况。

我们考虑解决一般的问题...

首先是有标号：我们先求有根，再除以 nn 就是无根了。考虑我们现在 SS 里的图都是有标号的。若根只在一个点双里，那么我们枚举这个点双里其他点伸出去的有根点双，然后再把根插进这个序列（乘 ∣G∣|G|）。然后把这个序列对应到枚举的 SS 中的点双 GG。但是这样会算重，因为 GG 有可能自同构（即编号变换之后得到的图是一样的），所以还要除以图的自同构数。

这里的 A(x)A(x) 是有标号合法有根图的指数生成函数，B(x)B(x) 上面求的那个东西，Aut(G)Aut(G) 是 GG 的自同构群。那么还是很好理解的。如果根在多个点双，那么就是 B(x)B(x) 构成的集合，exp一下就行。
所以有 A(x)=xexp⁡B(x)A(x)=x\exp B(x)，再把 BB 替换一下就可以牛顿迭代了。

无标号：这时要用到polya定理。我们考虑枚举根在点双 GG 中的位置，然后其他点在组合的时候要用polya定理。在 GG 中固定了根 vv 之后的同构群其实就是 Aut(G)Aut(G) 中满足 g(v)=vg(v)=v 的置换，设为 HvH_v。（在仙人掌计数中体现就是我们枚举了根之后满足 g(v)=vg(v)=v 的置换只有翻转，那么只需考虑序列的翻转。）
但是枚举根也会重，如果两个点在图中的同一等价类，那么就只能枚举一次。解决办法：
这里的 tt 是
这个最后的结果十分漂亮，不过计算量会很大。不过我们在实际用的时候，在一些可以分析的图上面，大概只需要注意一下每一个等价类只枚举一次就行了。（比如仙人掌，它环只有一个等价类，就没有枚举了）
对很多边双的组合只需要再做一次欧拉变换（polya指数）

不过这个只有理论价值，毕竟仙人掌计数牛顿迭代的式子就已经很不能看了/cy

如果无标号无根仙人掌的话，考虑它的圆方树。一棵仙人掌与圆方树对应的话，是要求方点的儿子是有顺序的，但是可翻转，而圆点的儿子没有顺序。这时我们转而对圆方树计数。而圆方树的重心也有很好的性质，所以与无标号无根树类似，是可以数无标号无根仙人掌的。

点双、边双计数
我的理解是还是用上面限制点双形态的那个计数技巧，得到点双、边双生成函数与一般联通图生成函数的关系，这时再来套拉格朗日反演就可以 O(nlog⁡n)O(n\log n) 求单点了。