「ZJOI2018」树
「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);
}
#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);
}
吐槽：实际上套了板子之后已经比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());
}
