POJ 2888 Magic Bracelet(Polya计数+dp+矩阵快速幂+欧拉函数+乘法逆元)

题目链接：
POJ 2888 Magic Bracelet
题意:
有一串$n$个珠子的项链，用$m$种颜色来染，有$k$个限制条件:$a[i]和b[i]$不能相邻。问本质不同的项链有多少种？(考虑旋转，答案对$9973$取模，且$gcd(n,9973)=1$）。数据范围:$n \leq 10^9,1\leq m \leq 10,0 \leq k \leq \frac{m*(m-1)}{2}$
分析：
因为需要考虑旋转，所以很容易想到$Polya计数和Burnside引理。$
假设有$k$个循环，用$link[i][j]$表示第$i$种颜色能否和第$j$种颜色相邻:当$link[i][j]=1$，表示$i$和$j$不能相邻，否则可以相邻。无向边:$link[i][j]=link[j][i]$ 。
用$dp[k][i][j]$表示经过$k$个循环从第$i$种颜色转移到第$j$种颜色的方案数：

$$dp[k][i][j]=\sum_{p\leq m}(1-link[p][j])*dp[k-1][i][p]$$
初始化：$dp[1][i][j] = 1-link[i][j]$
初始矩阵:$A[i][j]=1-link[i][j]$来表示颜色间的连通性
由上面的递推式，$A^k[i][j]$代表的是:从第$i$种颜色经过$k$个循环后变为第$j$种颜色的方法数。考虑循环，需要知道从$i$ 出发回到$i$ 的方案数即：$A^{k}[i][i]$。$solve(k)$表示循环个数为$k$时的方案数:
$$solve(k)=\sum_{i\leq m}A^{k}[i][i]$$

离散数学里有：如果用0,1矩阵A来表示无向图的连通情况的话，A^k代表的就是一个点经过k条路后能到达的地方的方法数。

假设对于每个循环的步长为$i$，也就是$0,i,2i,...$构成一个循环。这个循环的周期为$\frac{i*n}{gcd(i,n)}$，所以这个循环有$\frac{n}{gcd(i,n)}$ 个元素,共有$gcd(i,n)$个循环。所以枚举的循环个数一定是$n$的因子,即：$k \ |\ n$。
满足循环个数为$k$的置换的旋转步长$i$满足$gcd(i,n)=k$，此种置换的个数也就是$gcd(\frac{i}{k},\frac{n}{k})=1$的$i$的个数，即:$\phi(\frac{n}{k})$.
综上对于每个满足$n \% k = 0$的$k$可以得到的方案数是

$$solve(k)*\phi (\frac{n}{k})$$
枚举每个$k$然后求和最后还要除以$n$(因为总的置换数为$n$),又因为有模数且模数为素数，那么就相当于乘以$n^{mod-2}$(费马小定理)。

$$Ans = \{ \sum_{i \mid n}solve(i)*\phi(\frac{n}{i})\} * n^{p - 2} \ \% p$$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const ll mod = 9973;

int link[15][15];

struct Matrix{
    int row, col;
    ll data[15][15];

    Matrix operator * (const Matrix& rhs) const { //矩阵相乘条件:col = rhs.row
        Matrix res;
        res.row = row, res.col = rhs.col;
        for(int i = 1; i <= res.row; ++i) {
            for(int j = 1; j <= res.col; ++j) {
                res.data[i][j] = 0;
                for(int k = 1; k <= row; ++k) {
                    res.data[i][j] += data[i][k] * rhs.data[k][j];
                } 
                res.data[i][j] %= mod;
                //要是写成: 
                //res.data[i][j] = (res.data[i][j] + data[i][k] * rhs.data[k][j] % mod) % mod,
                //会TLE，辣鸡
            }
        }
        return res;
    }

    Matrix operator ^ (const int m) const {  //矩阵快速幂
        Matrix res, tmp;
        res.row = row, res.col = col; //row == col
        memset(res.data, 0, sizeof(res.data));

        tmp.row = row, tmp.col = col;
        memcpy(tmp.data, data, sizeof(data));
        for(int i = 1; i <= res.row; ++i) { res.data[i][i] = 1; }
        int mm = m;
        while(mm) {
            if(mm & 1) res = res * tmp;
            tmp = tmp * tmp;
            mm >>= 1;
        }
        return res;
    }
};

inline ll quick_pow(ll a, ll b, ll p)
{
    ll res = 1, tmp = a;
    while(b) {
        if(b & 1) res = res * tmp % p;
        tmp = tmp * tmp % p;
        b >>= 1;
    }
    return res;
}

inline ll phi(int x)
{
    ll res = x;
    for(int i = 2; i * i <= x; ++i) {
        if(x % i == 0) {
            res = res / i * (i - 1);
            while(x % i == 0) { x /= i; }
        }
    }
    if(x > 1) res = res / x * (x - 1);
    return res;
}

inline ll solve(int len, int m)
{
    Matrix tmp;
    tmp.row = tmp.col = m;
    for(int i = 1; i <= m; ++i) {
        for (int j = 1; j <= m; ++j) {
            tmp.data[i][j] = 1 - link[i][j];
        }
    }
    tmp = tmp ^ len;

    ll ans = 0;
    for (int i = 1; i <= m; ++i) {
        ans = (ans + tmp.data[i][i]) % mod;
    }
    return ans;
}

int main()
{
    int T, n, m, t;
    scanf("%d", &T);
    while (T--) {
        memset(link, 0, sizeof(link));
        scanf("%d%d%d", &n, &m, &t);
        for(int i = 0; i < t; ++i) {
            int former, later;
            scanf("%d%d", &former, &later);
            link[former][later] = link[later][former] = 1;
        }

        ll ans = 0;
        for(int i = 1; i * i <= n; ++i) {
            if(n % i) continue;
            ans = (ans + phi(n / i) * solve(i, m)) % mod;
            if(n / i == i) continue;
            ans = (ans + phi(i) * solve(n / i, m)) % mod;
        }
        printf("%lld\n", ans * quick_pow(n % mod, mod - 2, mod) % mod);
    }
    return 0;
}