[BZOJ1488][HNOI2009]图的同构（Polya定理+组合数）

考虑把这个图视为 n n 阶完全图，每一条边可以染成两种颜色00和 1 1 。那么可以根据Polya定理：

ans=1|G|∑g∈Gtc(g)ans=1|G|∑g∈Gtc(g)

来求解。 G G 为置换的集合，tt为颜色数量， c(g) c ( g ) 为置换 g g 下的循环节个数。

然而，此题的置换个数|G||G|为 n! n ! ，暴力枚举是过不去的。但可以发现 c(g) c ( g ) 的值总是在 [1,n(n−1)2] [ 1 , n ( n − 1 ) 2 ] 之间。

因此可以考虑对于任何一个 i∈[1,n(n−1)2] i ∈ [ 1 , n ( n − 1 ) 2 ] 有多少种置换 g g 使得c(g)=ic(g)=i。

从 1 1 到nn的排列的循环节上去考虑。如果一个 1 1 到nn的排列有 m m 个循环节，循环节的长度分别为：

a1,a2,a3,...,ama1,a2,a3,...,am

那么通过推理可以得出，这时候边的循环节个数为：

∑mi=1⌊ai2⌋+∑1≤i<j≤mgcd(ai,aj) ∑ i = 1 m ⌊ a i 2 ⌋ + ∑ 1 ≤ i < j ≤ m gcd ( a i , a j )

也就是所有循环节长度的一半的下取整之和，再加上循环节长度的两两的 gcd gcd 之和。

由于 n=60 n = 60 时 n n 的正整数拆分方法数只有966467966467，所以可以暴力搜索每一个拆分后求出边的循环节个数。暴力出奇迹啊！

最后的问题就是循环节长度分别为 a1,a2,...,am a 1 , a 2 , . . . , a m 的排列有多少个。可以将 a a <script type="math/tex" id="MathJax-Element-5407">a</script>的值分类后用组合数学的知识求出（具体见代码）。

记得写个逆元

代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = , MX = , M = N * N;
int n, m, a[N], fac[N], inv[N], res[M];
int qpow(int a, int b) {
    int res = ;
    while (b) {
        if (b & ) res = res * a % MX;
        a = a * a % MX;
        b >>= ;
    }
    return res;
}
int C(int n, int m) {
    return fac[n] * inv[m] * inv[n - m] % MX;
}
void dfs(int orz, int dep, int minv) {
    if (!orz) {
        dep--; int i, j, cnt = ; for (i = ; i <= dep; i++) cnt += a[i] >> ;
        for (i = ; i < dep; i++) for (j = i + ; j <= dep; j++)
            cnt += __gcd(a[i], a[j]); int co = , ti = , tn = n;
        for (i = ; i <= dep; i++) {
            ti++; if (i == dep || a[i] != a[i + ]) {
                int de = ; for (j = ; j <= ti; j++)
                    de = de * C(tn, a[i]) * fac[a[i] - ] % MX, tn -= a[i];
                de = de * inv[ti] % MX; co = co * de % MX; ti = ;
            }
        }
        res[cnt] = (res[cnt] + co) % MX; return;
    }
    int i; for (i = minv; (i << ) <= orz; i++)
        a[dep] = i, dfs(orz - i, dep + , i);
    if (minv <= orz) a[dep] = orz, dfs(, dep + , orz);
}
int main() {
    int i, ans = ; cin >> n; fac[] = ; m = n * (n - ) >> ;
    if (n == ) return puts("1"), ;
    for (i = ; i <= n; i++) fac[i] = fac[i - ] * i % MX;
    inv[n] = qpow(fac[n], MX - ); for (i = n - ; i >= ; i--)
        inv[i] = inv[i + ] * (i + ) % MX;
    dfs(n, , ); for (i = ; i <= m; i++)
        ans = (ans + qpow(, i) * res[i]) % MX;
    cout << ans * qpow(fac[n], MX - ) % MX << endl;
    return ;
}