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

最新推荐文章于 2019-01-11 11:55:38 发布

原创最新推荐文章于 2019-01-11 11:55:38 发布 · 261 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#Polya计数

BZOJ、UOJ、LOJ 专栏收录该内容

341 篇文章

订阅专栏

本文介绍了一种利用Polya定理求解完全图中边的不同染色方案数量的方法，并详细解释了如何针对大规模问题进行优化，特别是当图的阶数较大时（如n=60），采用高效的算法策略。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

考虑把这个图视为 $n$ 阶完全图，每一条边可以染成两种颜色 $0$ 和 $1$ 。那么可以根据Polya定理：
$ans=\frac{1}{|G|}\sum_{g\in G}t^{c(g)}$
来求解。 $G$ 为置换的集合， $t$ 为颜色数量， $c(g)$ 为置换 $g$ 下的循环节个数。
然而，此题的置换个数 $|G|$ 为 $n!$ ，暴力枚举是过不去的。但可以发现 $c(g)$ 的值总是在 $[1,\frac{n(n-1)}{2}]$ 之间。
因此可以考虑对于任何一个 $i\in[1,\frac{n(n-1)}{2}]$ 有多少种置换 $g$ 使得 $c(g)=i$ 。
从 $1$ 到 $n$ 的排列的循环节上去考虑。如果一个 $1$ 到 $n$ 的排列有 $m$ 个循环节，循环节的长度分别为：
$a_1,a_2,a_3,...,a_m$
那么通过推理可以得出，这时候边的循环节个数为：
$\sum_{i=1}^m\lfloor\frac{a_i}{2}\rfloor+\sum_{1\leq i<j\leq m}\gcd(a_i,a_j)$
也就是所有循环节长度的一半的下取整之和，再加上循环节长度的两两的 $\gcd$ 之和。
由于 $n=60$ 时 $n$ 的正整数拆分方法数只有 $966467$ ，所以可以暴力搜索每一个拆分后求出边的循环节个数。~~暴力出奇迹啊！~~
最后的问题就是循环节长度分别为 $a_1,a_2,...,a_m$ 的排列有多少个。可以将的值分类后用组合数学的知识求出（具体见代码）。
~~记得写个逆元~~
代码：

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