题意
用 (m) 种颜色，给 (n) 个点的无向完全图的 (frac{n(n-1)}{2}) 条边染色，两种方案相同当且仅当一种方案交换一些点的编号后可以变成另一种方案。问有多少本质不同的染色方案。
(nle 53, mle 1000, n<modle 10^9) 且 (mod) 为质数。
分析
- 考虑 (Polya) 定理。
- 假设已经枚举了一个点置换（对应唯一一种边置换），能否快速求出对应边的置换的循环个数？
- 对于两个点的循环（设长度分别为 (l_1,l_2)），它们之间的边构成了 (frac{l_1l_2}{lcm(l_1,l_2)}=gcd(l_1,l_2)) 个边的循环。
- 对于一个点循环内部的边：假设初始在序列上选定的两个点的位置相差了 (x) ，在右端点转回序列左端后差值会变成 (n-x) 。可以得到，当且仅当 (2x=n) 时，循环长度为 (frac{n}{2})（无向图） ，其余时刻为 (n) 。所以当 (n) 为偶数时的边循环个数还要加 1。由于总边数为 (frac{n(n-1)}{2}) ，所以边循环个数总可以写成 (lfloorfrac{n}{2} floor) 的形式。
- 那么对于一种点置换，假设其所有循环满足 (l_1le l_2le cdotsle l_k) ，这样的置换个数为 (frac{n!}{l_1l_2cdots l_kS_1!S_2!cdots S_{max}!}) ，其中 (S_i) 表示长度为 (i) 的循环个数。为什么是 (l_i) 而不是 (l_i!) 的原因是每次我们找 (l_i) 个位置作为一个循环之后，他们都连了一条指向其他位置的边，构成了一个环。环排列的个数是 (frac{n!}{n}=(n-1)!) 。
- 综上，对于一种点置换，对应边置换的循环个数为 (sum_{i=1}^klfloorfrac{l_i}{2} floor+sum_{i=1}^{k-1}sum_{j=i+1}^kgcd(l_i,l_j)) 。
- 由于 (nle 53) ，拆分的方案数在一个可接受的范围内。所以爆搜循环的拆分即可。
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
#define go(u) for(int i = head[u], v = e[i].to; i; i=e[i].lst, v=e[i].to)
#define rep(i, a, b) for(int i = a; i <= b; ++i)
#define pb push_back
#define re(x) memset(x, 0, sizeof x)
inline int gi() {
int x = 0,f = 1;
char ch = getchar();
while(!isdigit(ch)) { if(ch == '-') f = -1; ch = getchar();}
while(isdigit(ch)) { x = (x << 3) + (x << 1) + ch - 48; ch = getchar();}
return x * f;
}
template <typename T> inline bool Max(T &a, T b){return a < b ? a = b, 1 : 0;}
template <typename T> inline bool Min(T &a, T b){return a > b ? a = b, 1 : 0;}
const int N = 60;
int n, m, mod;
LL fac[N], invfac[N], inv[N], ans;
int g[N][N], stk[N], num[N];
LL Pow(LL a, LL b) {
LL res = 1ll;
for(; b; b >>= 1, a = a * a % mod) if(b & 1) res = res * a % mod;
return res;
}
void add(LL &a, LL b) {
a += b;
if(a >= mod) a -= mod;
}
void solve( int tp) {
re(num);
rep(i, 1, tp) num[stk[i]]++;
LL cnt = fac[n];
rep(i, 1, tp) cnt = cnt * inv[stk[i]] % mod;
rep(i, 1, n) cnt = cnt * invfac[num[i]] % mod;
LL c = 0;
rep(i, 1, tp) c = (c + stk[i] / 2) % (mod - 1);
rep(i, 1, tp - 1)
rep(j, i + 1, tp) c = (c + g[stk[i]][stk[j]]) % (mod - 1);
add(ans, cnt * Pow(m, c) % mod);
}
void dfs(int dep, int rest, int lst) {
if(!rest) {
solve(dep - 1);
return;
}
for(int i = lst; i <= rest; ++i) stk[dep] = i, dfs(dep + 1, rest - i, i);
}
int gcd(int a, int b) {
return !b ? a : gcd(b, a % b);
}
int main() {
n = gi(), m = gi(), mod = gi();
rep(i, 0, n)
rep(j, 0, n) g[i][j] = gcd(i, j);
inv[1] = fac[0] = invfac[0] = 1;
rep(i, 1, n) {
if(i ^ 1) inv[i] = (LL) (mod - mod / i) *inv[mod % i] % mod;
fac[i] = (LL) fac[i - 1] * i % mod;
invfac[i] = (LL) invfac[i - 1] * inv[i] % mod;
}
inv[0] = 1;
dfs(1, n, 1);
printf("%lld
", ans * invfac[n] % mod);
return 0;
}