「ZJOI2018」树
前言
置换同构计数真是令人头大，感觉依然不是特别懂
分析与初步构想
设按照题意生成的(n)个节点有标号有根树族为(mathcal{T}_n)，对于某种树形(T)的生成方案数为(c(T))
则答案显然是(cfrac{sum _{Tin mathcal{T}_n} c^k(T)}{(n-1)!^k})
(k)的量级意味着我们无法完成有关(c^k(T))的展开，因此要在最开始的计算中就将(k)这个幂次加入
为此定义(F_{i,n}=sum _{Tin mathcal{T}_n} c^{ik}(T))
则我们要计算的答案就是(cfrac{F_{1,n}}{(n-1)!^k})
容易发现有意义的(F_{i,j})规模为(O(sum sum [ijleq n])=O(nln n))
关于为什么(F_{i,n})会有第一维，会在下面出现
ps: 我的两个维度好像和别人是反的
计算过程分析
考虑递推计算(F_{i,n})
容易发现对于任意一个(F_{i,n})的计算只需要令(k ightarrow ik)
为了简化描述，我们先以计算(F_{1,...})为例
对于同构树计数，子树的叠加是无序的背包问题，树形之间存在着置换同构
而每棵子树的编号之间可以任意归并，因此需要( ext{EGF})的背包合并每棵子树
显然同构仅出现在同种大小的子树中，因此可以对于不同大小的子树分离
假设对于大小为(n)的子树，出现了(m)种不同的树形(T_1,T_2,ldots T_m)，每种出现了(a_i)个
则答案应该为(displaystyle sum_{T_i e T_j} frac{1}{(n!)^m}prod frac{c^{a_ik}(T_i)}{(a_i!)^k})
其中((a_i!)^k)除去同构树之间无排列序的方案数
普通的背包计数难以处理(T_i e T_j)的限制
因此考虑用( ext{Burnside})引理解决置换同构问题
同构出现于(m)棵树之间，因此我们置换群为对于(m)的排列置换群
显然任意一个排列置换的的结果是若干个置换环，环上的树树形相同
则对于一个置换(f)，设其生成了大小分别为(b_1,b_2,ldots, b_m)的置换环
理想情况下其不动点的式子计算如下
(displaystyle ext{fix}(f)=frac{(nsum b_i)!}{(n!)^{sum b_i}}cdot frac{F_{b_i,n}}{(b_i!)^k})
（这就是为什么要给(F)添加第一维）
其中(cfrac{(nsum b_i)!}{(n!)^{sum b_i}})处理了(sum b_i)棵树的点编号，实际上这两个权值可以在计算的最后加入，因此下面忽略掉
而实际上，在不动点中直接加入(cfrac{1}{(b_i!)^k})是错的
原因在于让这样的 置换环权值 带入( ext{Burnside})引理之后
最终的 同构类权值 并不是我们想要的(cfrac{1}{(a_i!)^k})
考虑 待定求解 一个置换环系数的( ext{EGF})，设其为(A(x))
相较于上面的枚举(f)的式子，这里的计算中系数还需要考虑对于每种(b_i)，等价的置换(f)个数
这同样需要对于每棵子树的( ext{EGF})合并，系数为(cfrac{1}{b_i!})
而一个环上的点存在一个环排列((b_i-1)!)，两者合并即(cfrac{1}{b_i})
注意这一部分并未被加入(A(x))中
然而这也恰好使得( ext{Burnside})引理的系数(frac{1}{m!})和置换元素的( ext{EGF})系数相抵消
不妨设添加环系数( ext{EGF})的变换为(hat A_i=cfrac{A_i}{i})
所以，则最终的计算中( ext{Burnside})引理的答案就是( ext{exp}(hat A(x)))
我们希望最终一个 合法的同构类 的权值为(cfrac{1}{(a_i!)^k})
而容易发现实际上 最终的同构类 是由我们初始枚举的置换环 的( ext{exp})
因此(cfrac{1}{a_i!})应当是置换环系数经过环元素( ext{exp})叠加的结果
设(B(x)=cfrac{x^i}{(i!)^k},C_i=cfrac{A_i}{F_{n,i}})
则(B(x)= ext{exp}(hat C(x)))，对于(B(x))取(ln)得到(hat C(x))
加入系数得到我们前面所待定的(A(x))，即可进行最后的( ext{exp}(hat A(x)))计算
ps:
你会发现可以直接忽略环( ext{EGF})变换，全程只有(hat C(x),hat A(x))，在过渡中这个变换的系数直接消失了
在这里提到这个系数是为了避免不必要的误解，同时也强调其他时候使用( ext{Burnside})引理需要添加这个系数
算法实现简谈
倒序枚举(F_{i,j})的(i)，正序枚举大小(j)，边界条件自然是(F_{i,1}=1)
确定了一个(i)后，就可以预处理(ln) 求出置换环系数，这里有(frac{n}{i})个
按照每个(i)，上面式子中的(k ightarrow ik)
按照(j)从小到大计算(F_{i,j})，每次得到(F_{i,j})之后
计算关于大小为(j)的树的背包系数，这里系数的个数为(l=frac{n}{ij})
将先前的系数补上(F_{d,j})，再做( ext{exp})，最后把前面扔掉的(cfrac{1}{(n!)^{sum a_i}})补上，( ((nsum a_i)!)直接作为后面( ext{EGF})合并的系数）
然后将它补进前面累和的背包里，就能得到这一项的值
注意前面算的式子都是计算儿子的，最后还要加上自己的大小1
当然，计算( ext{exp},ln)需要下面的帮助
( ext{exp})的(O(n^2))方法
(F(x)= ext{exp}(G(x)))
(F'(x)= ext{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))，都做一次( ext{exp})，然后背包合并
但是实际上，我们可以先将(hat A(x))放在一起，然后一起做( ext{exp})
具体的，每次得到(F_{i,j})之和，我们就可以确定(sum hat A(x))的第(j)项
那么在维护(sum hat A(x))的同时，也依次递推( ext{exp}(sum hat A(x)))的(j)项
这样不仅去掉的背包的过程，也少了很多次( ext{exp})
Montegomery
复杂度分析
未优化
相对于( ext{exp})，求(ln)的复杂度可以忽略
而( ext{exp})每次大小是(frac{n}{ij})，即(O(sum sum cfrac{n^2}{i^2j^2})=O(n^2))
最后的复杂度反而在于背包合并( ext{EGF})，为(O(n^2ln n))
优化后
同步求( ext{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
",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
",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
",ans.get());
}