Magic Bracelet（POJ 2888，Burnside 引理 + 矩阵快速幂优化 dp）

一、题目链接

Magic Bracelet

二、题目大意

给定 $m$ 种不同颜色的魔法珠子，每种颜色的珠子的个数都足够多。

现在要从中挑选 $n$ 个珠子，串成一个环形魔法手链。

魔法珠子之间存在 $k$ 对排斥关系，互相排斥的两种颜色的珠子不能相邻，否则会发生爆炸。（同一种颜色的珠子之间也可能存在排斥）

请问一共可以制作出多少种不同的手链。

注意，如果两个手链经旋转后能够完全重合在一起，对应位置的珠子颜色完全相同，则视为同一种手链。

答案对 $99739973$ 取模。

$1\le T\le 101\le n\le 10^{9}gcd(n,9973)=11\le m\le 100\le k\le \frac{m(m+1)}{2}1\le a,b\le m$.

三、分析

由于循环之间存在限制关系，故无法使用 $P\acute{o}lya$ 定理，而只能使用 $Burnside$ 引理.

旋转操作共有 $n$ 种置换，考虑一个置换 $P_i$：顺时针旋转 $i$ 个珠子（$1\le i\le n$）.

通过二元一次不定方程解得循环长度为 $\frac{n}{gcd(n,i)}$，于是循环个数为 $gcd(n,i)$.

因此置换 $P_i$ 的不动点个数等于长度为 $gcd(n,i)$ 的环的合法方案数，记为 $F(gcd(n,i))$.

根据 $Burnside$ 引理可得答案为 $\begin{array}{r}ANS=\frac{1}{n}\sum _{i=1}^{n}F(gcd(n,i))\end{array}$.

由于 $n$ 最大 $10^9$，直接枚举 $i$ 是无法接受的. 而 $gcd(n,i)$ 最多只有 $1600$ 种，所以我们去枚举 $gcd(n,i)$.

于是 $\begin{array}{r}ANS=\frac{1}{n}\sum _{d|n}\phi (\frac{n}{d})F(d)\end{array}$，其中 $\phi (\frac{n}{d})$ 可以先预处理出 $10^5$ 以内的质数再 $O(\frac{\sqrt{n/d}}{ln(n/d)})$ 求出.

考虑如何计算 $F(d)$，记 $c[i][j]$ 表示颜色 $i$ 和颜色 $j$ 的珠子能否相邻，$f_x[i][j]$ 为第一颗珠子的颜色为 $x$ 且 确定了前 $i$ 颗珠子的颜色 且 第 $i$ 颗珠子颜色为 $j$ 的方案数.

则有递推式 $\begin{array}{r}{f}_x[i][j]=\sum _{c[j][k]}{f}_x[i-1][k](2\le i\le n1\le j,k\le m)\end{array}$.

于是 $\begin{array}{r}F(d)=\sum _x\sum _{c[x][j]}{f}_x[d][j]\end{array}$

直接 $dp$ 的时间复杂度是 $O(n{m}^3)$ 的，不难发现可以使用矩阵快速幂加速 $dp$ 递推.

即

$\left[\begin{array}{cccc}{f}_1[i][1]& f_1[i][2]& \cdots & f_1[i][m]\\ f_2[i][1]& f_2[i][2]& \cdots & f_2[i][m]\\ ⋮& ⋮& & ⋮\\ f_m[i][1]& f_m[i][2]& \cdots & f_m[i][m]\end{array}\right]=\left[\begin{array}{cccc}{f}_1[i-1][1]& f_1[i-1][2]& \cdots & f_1[i-1][m]\\ f_2[i-1][1]& f_2[i-1][2]& \cdots & f_2[i-1][m]\\ ⋮& ⋮& & ⋮\\ f_m[i-1][1]& f_m[i-1][2]& \cdots & f_m[i-1][m]\end{array}\right]\left[\begin{array}{cccc}c[1][1]& c[1][2]& \cdots & c[1][m]\\ c[2][1]& c[2][2]& \cdots & c[2][m]\\ ⋮& ⋮& & ⋮\\ c[m][1]& c[m][2]& \cdots & c[m][m]\end{array}\right]$

于是

$\left[\begin{array}{cccc}{f}_1[d][1]& f_1[d][2]& \cdots & f_1[d][m]\\ f_2[d][1]& f_2[d][2]& \cdots & f_2[d][m]\\ ⋮& ⋮& & ⋮\\ f_m[d][1]& f_m[d][2]& \cdots & f_m[d][m]\end{array}\right]={\left[\begin{array}{cccc}c[1][1]& c[1][2]& \cdots & c[1][m]\\ c[2][1]& c[2][2]& \cdots & c[2][m]\\ ⋮& ⋮& & ⋮\\ c[m][1]& c[m][2]& \cdots & c[m][m]\end{array}\right]}^{d-1}$

四、代码实现

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;

const int M = (int)1e5;
const int N = (int)1e5;
const int inf = 0x3f3f3f3f;
const double eps = 1e-8;
const int mod = (int)9973;

int n, m, k;
bool c[10][10];

struct Matrix
{
    int D[10][10];

    Matrix operator*(const Matrix& B)const
    {
        Matrix C;
        for(int i = 0; i < m; ++i)
        {
            for(int j = 0; j < m; ++j)
            {
                C.D[i][j] = 0;
                for(int k = 0; k < m; ++k)
                {
                    C.D[i][j] = (C.D[i][j] + D[i][k] * B.D[k][j]) % mod;
                }
            }
        }
        return C;
    }
};

Matrix quick(Matrix A, int b)
{
    Matrix S;
    if(b == 0)
    {
        for(int i = 0; i < m; ++i)
            for(int j = 0; j < m; ++j) S.D[i][j] = (i == j);
    }
    else
    {
        S = A; --b;
        while(b)
        {
            if(b & 1) S = S * A;
            A = A * A;
            b >>= 1;
        }
    }
    return S;
}

bool is_prime[M + 5];
int prime[M + 5], cnt;

void init()
{
    memset(is_prime, 1, sizeof(is_prime));
    is_prime[0] = is_prime[1] = 0;
    for(int i = 2; i <= M; ++i)
    {
        if(is_prime[i]) prime[++cnt] = i;
        for(int j = 1; j <= cnt && i * prime[j] <= M; ++j)
        {
            is_prime[i * prime[j]] = 0;
            if(i % prime[j] == 0) break;
        }
    }
}

int F(int d)
{
    int ans = 0;
    Matrix A;

    for(int i = 0; i < m; ++i)
    {
        for(int j = 0; j < m; ++j)
        {
            A.D[i][j] = c[i][j];
        }
    }
    A = quick(A, d - 1);
    for(int i = 0; i < m; ++i)
    {
        for(int j = 0; j < m; ++j)
        {
            if(c[i][j]) ans = (ans + A.D[i][j]) % mod;
        }
    }
    return ans;
}

int phi(int n)
{
    int ans = n;
    for(int i = 1; i <= cnt && 1ll * prime[i] * prime[i] <= n; ++i)
    {
        if(n % prime[i] == 0)
        {
            ans = ans / prime[i] * (prime[i] - 1);
            while(n % prime[i] == 0) n /= prime[i];
        }
    }
    if(n > 1) ans = ans / n * (n - 1);
    return ans;
}

int quick(int a, int b, int p = mod)
{
    a %= p;
    int s = 1;
    while(b)
    {
        if(b & 1) s = s * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return s % p;
}

int inv(int n, int p = mod)
{
    return quick(n, p - 2, p);
}

void work()
{
    scanf("%d %d %d", &n, &m, &k);
    for(int i = 0; i < m; ++i) for(int j = 0; j < m; ++j) c[i][j] = 1;
    for(int i = 1, a, b; i <= k; ++i)
    {
        scanf("%d %d", &a, &b);
        --a, --b;
        c[a][b] = c[b][a] = 0;
    }

    int ans = 0;
    for(int d = 1; d * d <= n; ++d)
    {
        if(n % d) continue;
        ans = (ans + phi(n / d) % mod * F(d)) % mod;
        if(d != n / d) ans = (ans + phi(d) % mod * F(n / d)) % mod;
    }
    ans = ans * inv(n) % mod;
    printf("%d\n", ans);
}

int main()
{
//    freopen("input.txt", "r", stdin);
//    freopen("output.txt", "w", stdout);
    init();
    int T; scanf("%d", &T);
    while(T--) work();
//    work();
    return 0;
}