「ZJOI2018」树

前言

置换同构计数真是令人头大，感觉依然不是特别懂

\[ \ \]

分析与初步构想

---

设按照题意生成的 \(n\) 个节点有标号有根树族为 \(\mathcal{T}_n\) ，对于某种树形 \(T\) 的生成方案数为 \(c(T)\)

则答案显然是 \(\cfrac{\sum _{T\in \mathcal{T}_n} c^k(T)}{(n-1)!^k}\)

\(k\) 的量级意味着我们无法完成有关 \(c^k(T)\) 的展开，因此要在最开始的计算中就将 \(k\) 这个幂次加入

为此定义 \(F_{i,n}=\sum _{T\in \mathcal{T}_n} c^{ik}(T)\)

则我们要计算的答案就是 \(\cfrac{F_{1,n}}{(n-1)!^k}\)

容易发现有意义的 \(F_{i,j}\) 规模为 \(O(\sum \sum [ij\leq n])=O(n\ln n)\)

关于为什么 \(F_{i,n}\) 会有第一维，会在下面出现

ps: 我的两个维度好像和别人是反的

计算过程分析

---

考虑递推计算 \(F_{i,n}\)

容易发现对于任意一个 \(F_{i,n}\) 的计算只需要令 \(k\rightarrow ik\)

为了简化描述，我们先以计算 \(F_{1,...}\) 为例

对于同构树计数，子树的叠加是无序的背包问题，树形之间存在着置换同构

而每棵子树的编号之间可以任意归并，因此需要 \(\text{EGF}\) 的背包合并每棵子树

显然同构仅出现在同种大小的子树中，因此可以对于不同大小的子树分离

假设对于大小为 \(n\) 的子树，出现了 \(m\) 种不同的树形 \(T_1,T_2,\ldots T_m\) ，每种出现了 \(a_i\) 个

则答案应该为 \(\displaystyle \sum_{T_i\ne T_j} \frac{1}{(n!)^m}\prod \frac{c^{a_ik}(T_i)}{(a_i!)^k}\)

其中 \((a_i!)^k\) 除去同构树之间无排列序的方案数

普通的背包计数难以处理 \(T_i\ne T_j\) 的限制

因此考虑用 \(\text{Burnside}\) 引理解决置换同构问题

---

同构出现于 \(m\) 棵树之间，因此我们置换群为对于 \(m\) 的排列置换群

显然任意一个排列置换的的结果是若干个置换环，环上的树树形相同

则对于一个置换 \(f\) ，设其生成了大小分别为 \(b_1,b_2,\ldots, b_m\) 的置换环

理想情况下其不动点的式子计算如下

\(\displaystyle \text{fix}(f)=\frac{(n\sum b_i)!}{(n!)^{\sum b_i}}\cdot \frac{F_{b_i,n}}{(b_i!)^k}\)

（这就是为什么要给 \(F\) 添加第一维）

其中 \(\cfrac{(n\sum b_i)!}{(n!)^{\sum b_i}}\) 处理了 \(\sum b_i\) 棵树的点编号，实际上这两个权值可以在计算的最后加入，因此下面忽略掉

而实际上，在不动点中直接加入 \(\cfrac{1}{(b_i!)^k}\) 是错的

原因在于让这样的 置换环权值 带入 \(\text{Burnside}\) 引理之后

最终的 同构类权值 并不是我们想要的 \(\cfrac{1}{(a_i!)^k}\)

---

考虑 待定求解 一个置换环系数的 \(\text{EGF}\) ，设其为 \(A(x)\)

相较于上面的枚举 \(f\) 的式子，这里的计算中系数还需要考虑对于每种 \(b_i\) ，等价的置换 \(f\) 个数

这同样需要对于每棵子树的 \(\text{EGF}\) 合并，系数为 \(\cfrac{1}{b_i!}\)

而一个环上的点存在一个环排列 \((b_i-1)!\) ，两者合并即 \(\cfrac{1}{b_i}\)

注意这一部分并未被加入 \(A(x)\) 中

然而这也恰好使得 \(\text{Burnside}\) 引理的系数 \(\frac{1}{m!}\) 和置换元素的 \(\text{EGF}\) 系数相抵消

不妨设添加环系数 \(\text{EGF}\) 的变换为 \(\hat A_i=\cfrac{A_i}{i}\)

所以，则最终的计算中 \(\text{Burnside}\) 引理的答案就是 \(\text{exp}(\hat A(x))\)

---

我们希望最终一个 合法的同构类 的权值为 \(\cfrac{1}{(a_i!)^k}\)

而容易发现实际上 最终的同构类 是由我们初始枚举的置换环 的 \(\text{exp}\)

因此 \(\cfrac{1}{a_i!}\) 应当是置换环系数经过环元素 \(\text{exp}\) 叠加的结果

设 \(B(x)=\cfrac{x^i}{(i!)^k},C_i=\cfrac{A_i}{F_{n,i}}\)

则 \(B(x)=\text{exp}(\hat C(x))\) ，对于 \(B(x)\) 取 \(\ln\) 得到 \(\hat C(x)\)

加入系数得到我们前面所待定的 \(A(x)\) ，即可进行最后的 \(\text{exp}(\hat A(x))\) 计算

ps:

你会发现可以直接忽略环 \(\text{EGF}\) 变换，全程只有 \(\hat C(x),\hat A(x)\) ，在过渡中这个变换的系数直接消失了

在这里提到这个系数是为了避免不必要的误解，同时也强调其他时候使用 \(\text{Burnside}\) 引理需要添加这个系数

算法实现简谈

---

倒序枚举 \(F_{i,j}\) 的 \(i\) ，正序枚举大小 \(j\) ，边界条件自然是 \(F_{i,1}=1\)

确定了一个 \(i\) 后，就可以预处理 \(\ln\) 求出置换环系数，这里有 \(\frac{n}{i}\) 个

按照每个 \(i\) ，上面式子中的 \(k\rightarrow ik\)

按照 \(j\) 从小到大计算 \(F_{i,j}\) ，每次得到 \(F_{i,j}\) 之后

计算关于大小为 \(j\) 的树的背包系数，这里系数的个数为 \(l=\frac{n}{ij}\)

将先前的系数补上 \(F_{d,j}\) ，再做 \(\text{exp}\) ，最后把前面扔掉的 \(\cfrac{1}{(n!)^{\sum a_i}}\) 补上，( \((n\sum a_i)!\) 直接作为后面 \(\text{EGF}\) 合并的系数）

然后将它补进前面累和的背包里，就能得到这一项的值

注意前面算的式子都是计算儿子的，最后还要加上自己的大小1

当然，计算 \(\text{exp},\ln\) 需要下面的帮助

\(\text{exp}\) 的 \(O(n^2)\) 方法

\(F(x)=\text{exp}(G(x))\)

\(F'(x)=\text{exp}(G(x))G'(x)\)

\(F'(x)=F(x)G'(x)\)

先计算出 \(G'(x)\) ，然后 \(O(n^2)\) 依次得到 \(F'(x)\) 的第 \(i\) 项，就能知道 \(F(x)\) 的第 \(i+1\) 项

\(\ln\) 的 \(O(n^2)\) 方法

\(F(x)=\ln G(x)\)

\(F'(x)G(x)=G'(x)\)

边界 \([x^0]G(x)=1,[x^0]F(x)=0\)

暴力推 \(F'(x)\) 的每一项即可

进一步优化

上面的计算时，每次求得一个 \(\hat A(x)\) ，都做一次 \(\text{exp}\) ，然后背包合并

但是实际上，我们可以先将 \(\hat A(x)\) 放在一起，然后一起做 \(\text{exp}\)

具体的，每次得到 \(F_{i,j}\) 之和，我们就可以确定 \(\sum \hat A(x)\) 的第 \(j\) 项

那么在维护 \(\sum \hat A(x)\) 的同时，也依次递推 \(\text{exp}(\sum \hat A(x))\) 的 \(j\) 项

这样不仅去掉的背包的过程，也少了很多次 \(\text{exp}\)

Montegomery

最后是喜闻乐见的套板子时间

\[ \ \]

复杂度分析

---

未优化

相对于 \(\text{exp}\) ，求 \(\ln\) 的复杂度可以忽略

而 \(\text{exp}\) 每次大小是 \(\frac{n}{ij}\) ，即 \(O(\sum \sum \cfrac{n^2}{i^2j^2})=O(n^2)\)

最后的复杂度反而在于背包合并 \(\text{EGF}\) ，为 \(O(n^2\ln n)\)

优化后

同步求 \(\text{exp}\) 的复杂度为 \(O((\cfrac{n}{i})^2)=O(n^2)\)

Code1:

const int N=2010;

int n,k,P;
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int F[N][N],T[N],H[N],C[N];
int I[N],J[N],Inv[N];
void Exp(int *a,int n){
	static int b[N];
	rep(i,1,n) b[i-1]=1ll*a[i]*i%P;
	rep(i,0,n) a[i]=0;
	a[0]=1;
	rep(i,0,n-1) {
		int s=0;
		rep(j,0,i) s=(s+1ll*a[i-j]*b[j])%P;
		a[i+1]=1ll*s*Inv[i+1]%P;
	}
}

void Ln(int *a,int n){
	static int b[N];
	rep(i,0,n) b[i]=a[i];
	rep(i,0,n-1) {
		int s=1ll*b[i+1]*(i+1)%P;
		rep(j,0,i-1) s=(s-1ll*a[j]*b[i-j])%P;
		Mod2(s),a[i]=s;
	}
	drep(i,n,1) a[i]=1ll*a[i-1]*Inv[i]%P;
	a[0]=0;
}

int IK[N],JK[N];

int main(){
	n=rd(),k=rd(),P=rd();
	Inv[0]=Inv[1]=1;
	rep(i,2,n) Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
	rep(i,*I=*J=1,n) I[i]=1ll*I[i-1]*qpow(Inv[i],k)%P,J[i]=1ll*J[i-1]*qpow(i,k)%P;
	
	drep(i,n,1) {
		int m=n/i;
		rep(j,*H=1,m) H[j]=0;
		rep(j,0,m) IK[j]=qpow(I[j],i),JK[j]=qpow(J[j],i);
		rep(j,*C=1,m) C[j]=IK[j];
		Ln(C,m);
		rep(j,1,m) {
			F[i][j]=1ll*H[j-1]*JK[j-1]%P;
			
			int u=m/j;
			T[0]=0;
			rep(x,1,u) T[x]=1ll*C[x]*F[i*x][j]%P;
			Exp(T,u);
			int t=1;
			rep(x,1,u) t=1ll*t*IK[j]%P,T[x]=1ll*T[x]*t%P;
			drep(x,m,1) rep(y,1,x/j) H[x]=(H[x]+1ll*H[x-y*j]*T[y])%P;
		}
	}
	int ans=1ll*F[1][n]*I[n-1]%P;
	printf("%d\n",ans);
}

Code2:

#include<cstdio>
using namespace std;
typedef long long ll;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
enum{N=2010};
int n,k,P;
ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}
int H[N],F[N][N],T[N],C[N];
int I[N],J[N],Inv[N],Fac[N];
int IK[N],JK[N];
void Ln(int *a,int n){
	static int b[N];
	rep(i,0,n) b[i]=a[i];
	rep(i,0,n-1) {
		int s=1ll*b[i+1]*(i+1)%P;
		rep(j,0,i-1) s=(s-1ll*a[j]*b[i-j])%P;
		Mod2(s),a[i]=s;
	}
	drep(i,n,1) a[i]=1ll*a[i-1]*Inv[i]%P;
	a[0]=0;
}


int main(){
	scanf("%d%d%d",&n,&k,&P);
	Inv[0]=Inv[1]=1;
	rep(i,2,n) Inv[i]=1ll*(P-P/i)*Inv[P%i]%P;
	rep(i,*I=*J=1,n) I[i]=1ll*I[i-1]*qpow(Inv[i],k)%P,J[i]=1ll*J[i-1]*qpow(i,k)%P;
	rep(i,*Fac=1,n) Fac[i]=1ll*Fac[i-1]*i%P;
	
	drep(i,n,1) {
		int m=n/i;
		rep(j,*H=1,m) H[j]=T[j]=0;
		rep(j,0,m) IK[j]=qpow(I[j],i),JK[j]=qpow(J[j],i);
		rep(j,*C=1,m) C[j]=IK[j];
		Ln(C,m);
		rep(j,1,m) {
			F[i][j]=1ll*H[j-1]*JK[j-1]%P;
			int u=m/j,t=1;
			T[0]=0;
			rep(x,1,u) t=1ll*t*IK[j]%P,T[x*j]=(T[x*j]+1ll*C[x]*F[i*x][j]%P*t%P*x*j)%P;
			rep(x,1,j) H[j]=(H[j]+1ll*H[j-x]*T[x])%P;
			H[j]=1ll*H[j]*Inv[j]%P;
		}
	}
	int ans=1ll*F[1][n]*I[n-1]%P;
	printf("%d\n",ans);
}

Code3:

Loj Submission

吐槽：实际上套了板子之后已经比loj上的所有人都快了

但是由于新旧评测机的问题~~~，总时间就显得慢了

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

enum{N=2010};
int n,k,P;
using u32=uint32_t;
using i32=int32_t;
using u64=uint64_t;
using i64=int64_t;

static u32 m,m2,inv,r2;
u32 getinv(){
    u32 inv=m;
    for(int i=0;i<4;++i) inv*=2-inv*m;
    return inv;
}
struct Mont{
private :
    u32 x;
public :
    static u32 reduce(u64 x){ 
        u32 y=(x+u64(u32(x)*inv)*m)>>32;
        return i32(y)<0?y+m:y;
    }
    Mont(){ ; }
    Mont(i32 x):x(reduce(u64(x)*r2)) { }
    Mont& operator += (const Mont &rhs) { return x+=rhs.x-m2,i32(x)<0&&(x+=m2),*this; }
    Mont& operator -= (const Mont &rhs) { return x-=rhs.x,i32(x)<0&&(x+=m2),*this; }
    Mont& operator *= (const Mont &rhs) { return x=reduce(u64(x)*rhs.x),*this; }
    friend Mont operator + (Mont x,const Mont &y) { return x+=y; }
    friend Mont operator - (Mont x,const Mont &y) { return x-=y; }
    friend Mont operator * (Mont x,const Mont &y) { return x*=y; }
    i32 get(){ 
        u32 res=reduce(x);
        return res>=m?res-m:res;
    }
} H[N],F[N][N],T[N],C[N],I[N],J[N],Inv[N],IK[N],JK[N];
Mont qpow(Mont x,ll k=P-2) {
	Mont res(1);
	for(;k;k>>=1,x*=x) if(k&1) res*=x;
	return res;
}
void Init(int m) { 
    ::m=m,m2=m*2;
    inv=-getinv();
    r2=-u64(m)%m;
}

void Ln(Mont *a,int n){
	static Mont b[N];
	rep(i,0,n) b[i]=a[i];
	rep(i,0,n-1) {
		Mont s=b[i+1]*(i+1);
		rep(j,0,i-1) s-=a[j]*b[i-j];
		a[i]=s;
	}
	drep(i,n,1) a[i]=a[i-1]*Inv[i];
	a[0]=0;
}

int main(){
	scanf("%d%d%d",&n,&k,&P),Init(P);
	Inv[0]=Inv[1]=1;
	rep(i,2,n) Inv[i]=(P-P/i)*Inv[P%i];
	I[0]=J[0]=1;
	rep(i,1,n) I[i]=I[i-1]*qpow(Inv[i],k),J[i]=J[i-1]*qpow(i,k);
	
	drep(i,n,1) {
		int m=n/i;
		rep(j,1,m) T[j]=0;
		rep(j,0,m) IK[j]=qpow(I[j],i),JK[j]=qpow(J[j],i);
		C[0]=H[0]=1;
		rep(j,1,m) C[j]=IK[j];
		Ln(C,m);
		rep(j,1,m) {
			F[i][j]=H[j-1]*JK[j-1];
			Mont t=1;
			rep(x,1,m/j) t=t*IK[j],T[x*j]+=C[x]*F[i*x][j]*t*(x*j);
			H[j]=0;
			rep(x,1,j) H[j]+=H[j-x]*T[x];
			H[j]=H[j]*Inv[j];
		}
	}
	Mont ans=F[1][n]*I[n-1];
	printf("%d\n",ans.get());
}