POJ 2888 Burnside 引理 + 矩阵快速幂 + DFS

题意

传送门 POJ 2888 Magic Bracelet

题解

根据 Burnside 引理，令 $C(i)$ 代表旋转 $i$ 个单位的置换下的不动点，则答案为
$$\frac{1}{n}\sum _0^{n-1}C(i)$$ 旋转 $i$ 个单位的置换下，循环有 $gcd(i,n)$ 个，每个循环的元素在项链上距离为 $gcd(i,n)$，项链上相邻的 $gcd(i,n)$ 个珠子属于不同的置换。那么求解 $C(i)$ 转化为求解大小为 $gcd(i,n)$ 的环上，合法的染色方案数。

任取环上一点为 $DP$ 起点，枚举其前一个节点的颜色，利用冲突关系构造的矩阵进行矩阵快速幂，统计答案即可。

满足 $gcd(i,n)=d$ 的 $i$ 共有 $\varphi (n/d)$ 个，合并相同的 $C(i)$ 即可，答案为
$$\frac{1}{n}\sum _{d|n}\varphi (d)C(n/d)$$ 预处理出质因子，$DFS$ 求解欧拉函数。总时间复杂度 $O(\sqrt{n}+{\sum }_{d|n}{m}^3\log n)$。

#include <algorithm>
#include <cstring>
#include <iostream>
using namespace std;
#define rep(i, l, r) for (int i = l, _ = r; i < _; ++i)
const int MAXLG = 35, MAXM = 11, MOD = 9973;
int T, N, M, K;
int pn, prime[MAXLG], num[MAXLG], Res;
struct matrix
{
    int a[MAXM][MAXM], n;
    matrix() {}
    matrix(int n) : n(n) { memset(a, 0, sizeof(a)); }
    void operator*=(const matrix &o)
    {
        matrix res(n);
        rep(i, 0, n) rep(k, 0, n) rep(j, 0, n) res.a[i][j] = (res.a[i][j] + a[i][k] * o.a[k][j]) % MOD;
        memcpy(a, res.a, sizeof(a));
    }
} tran;

matrix mpow(matrix a, int n)
{
    matrix res(a.n);
    rep(i, 0, res.n) res.a[i][i] = 1;
    while (n)
    {
        if (n & 1)
            res *= a;
        a *= a, n >>= 1;
    }
    return res;
}

int pow_mod(int x, int n, int mod)
{
    int res = 1;
    x %= mod;
    while (n)
    {
        if (n & 1)
            res = res * x % mod;
        x = x * x % mod, n >>= 1;
    }
    return res;
}

void dfs(int k, int n, int ph, int d)
{
    if (k == pn)
    {
        int t = N / d;
        matrix a = mpow(tran, t);
        int res = 0;
        rep(i, 0, M)  res = (res + a.a[i][i]) % MOD;
        Res = (Res + ph % MOD * res % MOD) % MOD;
        return;
    }
    int p = prime[k];
    if (n < num[k])
        dfs(k, n + 1, ph * (n == 0 ? p - 1 : p), d * p);
    dfs(k + 1, 0, ph, d);
}

int main()
{
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> T;
    while (T--)
    {
        cin >> N >> M >> K;
        tran.n = M;
        rep(i, 0, M) rep(j, 0, M) tran.a[i][j] = 1;
        rep(i, 0, K)
        {
            int u, v;
            cin >> u >> v;
            --u, --v;
            tran.a[u][v] = tran.a[v][u] = 0;
        }
        pn = Res = 0;
        int t = N;
        for (int i = 2; i * i <= t; ++i)
            if (t % i == 0)
            {
                prime[pn] = i, num[pn] = 0;
                while (t % i == 0)
                    t /= i, ++num[pn];
                ++pn;
            }
        if (t != 1)
            prime[pn] = t, num[pn++] = 1;
        dfs(0, 0, 1, 1);
        cout << Res * pow_mod(N, MOD - 2, MOD) % MOD << '\n';
    }
    return 0;
}