「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());
}

orangejuice's blog

按ctrl会跳搜索，右上角选择day/night mode 似乎修好了搜索

「ZJOI2018」树

「ZJOI2018」树

前言

分析与初步构想

计算过程分析

算法实现简谈

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

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

进一步优化

Montegomery

复杂度分析

未优化

优化后