分析
三倍经验题，本文以[BZOJ1478][SGU282]Isomorphism为例展开叙述，主体思路与另外两题大（wan）致（quan）相（yi）同（zhi）。
这可能是博主目前写过最长也是最认真的题解了。
题目中规定“若两个已染色的图，其中一个图可以通过结点重新编号而与另一个图完全相同， 就称这两个染色方案相同”，说明这个置换群是定义在点上的，而染色方案是定义在边上的。把边的染色方案转化为点的染色方案不太现实，所以说我们可以考虑如何将点的置换转化为边的置换。
一个显然的结论是点的置换和边的置换是一一对应的，但是我们不能暴力枚举所有的点的置换，因为可以想到点的置换又可以和(1 sim n)的排列一一对应，推出共有(n!)个点的置换。换一个角度考虑，一个点的置换是由多个循环节组成的，所以我们可以暴力枚举每个循环节的长度，直接由此转化为边的置换，计算出此时的边的置换的循环节个数后乘一个满足这种循环节长度组合的点的置换的个数，这样的复杂度是可以接受的。
假设我们枚举到这样一个数组(cnt)，(cnt[i])表示在这种点的置换中循环节长度为(i)的循环节有(cnt[i])个。
考虑一个长度为(l)的点的置换的循环节，现在有两个结点都在这个循环节上，这样的两个结点会形成(lfloor frac{l}{2} floor)个边的置换的循环节。
为什么？
1----6 考虑这样一个点的置换的循环节(1,2,3,4,5,6)（逆时针），
/ 会形成6/2=3个边的置换的循环节，分别为：
/ ((1,2),(2,3),(3,4),(4,5),(5,6),(6,1))
2 5 ((1,3),(2,4),(3,5),(4,6),(5,1),(6,2))
/ ((1,4),(2,5),(3,6))
/ 由于边是无向边，以上便是全部的3个边的置换的循环节。
3----4 点的置换的循环节长度为奇数时情况稍有不同。
如果两个结点中有一个在另一个长度为(l')的循环节上，那么这样的两个结点会形成(frac{l imes l'}{lcm(l,l')}=gcd(l,l'))个边的置换的循环节。
把所有以上两种边的置换的循环节的个数加起来就好了。
以上，我们已经得到了此时的边的置换的循环节个数，接下来我们想求出满足这种循环节长度组合的点的置换的个数。（又开始啰嗦了）
假设我们有一个(1 sim n)的排列（注意这里要和前文所提到的“点的置换可以和(1 sim n)的排列一一对应”进行区分，它们不一样！），可以从左到右在这个排列中截出(cnt[1])个(1)，(cnt[2])个(2)，......，以此类推。每一段截出来的数字我们可以认为这是一个点的置换的循环节。我们知道长度(1 sim n)的排列有(n!)个，又注意到循环节关心的是每个数字的下一个是什么，而不关心第几个数字是什么，所以要排列数要除以(prod i^{cnt[i]})。又双注意到对于长度相同的几个循环节，我们其实不会区分它们，更不会考虑它们的先后顺序，所以又要除以(prod (cnt[i]!))。综上所述，满足这种循环节长度组合的点的置换的个数为:
根据Polya定理，把求出来的两坨东西搞一搞，最后乘一个((n!)^{-1})即可。
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>
#include <algorithm>
#define rin(i,a,b) for(int i=(a);i<=(b);i++)
#define rec(i,a,b) for(int i=(a);i>=(b);i--)
#define trav(i,x) for(int i=head[(x)];i;i=e[i].nxt)
using std::cin;
using std::cout;
using std::endl;
typedef long long LL;
const int MAXN=55;
const int MAXM=1005;
int n,m;LL MOD,ans=0;
int cnt[MAXN];
LL fac[MAXN];
inline LL gcd(LL x,LL y){
if(!x||!y) return x+y;
while(y){
std::swap(x,y);
y%=x;
}
return x;
}
inline LL qpow(LL x,LL y){
LL ret=1,tt=x%MOD;
while(y){
if(y&1) ret=ret*tt%MOD;
tt=tt*tt%MOD;
y>>=1;
}
return ret;
}
void dfs(int rem,int now){
if(!rem){
LL temp=0,mot=1;
rin(i,1,now-1){
if(!cnt[i]) continue;
temp+=i/2*cnt[i];
temp+=i*cnt[i]*(cnt[i]-1)/2;
rin(j,i+1,now-1){
if(!cnt[j]) continue;
temp+=cnt[i]*cnt[j]*gcd(i,j);
}
}
rin(i,1,now-1){
if(!cnt[i]) continue;
mot=mot*qpow(i,cnt[i])%MOD*fac[cnt[i]]%MOD;
}
ans=(ans+qpow(m,temp)*fac[n]%MOD*qpow(mot,MOD-2))%MOD;
return;
}
if(now>rem) return;
for(int i=0;i*now<=rem;i++){
cnt[now]=i;
dfs(rem-i*now,now+1);
}
}
int main(){
scanf("%d%d%lld",&n,&m,&MOD);
fac[0]=1;
rin(i,1,n) fac[i]=fac[i-1]*i%MOD;
dfs(n,1);
printf("%lld
",ans*qpow(fac[n],MOD-2)%MOD);
return 0;
}