[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 ;
}