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$ 种颜色的方案数：

d p [k] [i] [j] = \sum p \leq m (1 ? l i n k [p] [j]) ? d p [k ? 1] [i] [p]

$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] $dp[1][i][j] = 1-link[i][j]$
初始矩阵:

A[i][j]=1?link[i][j] $A[i][j]=1-link[i][j]$ 来表示颜色间的连通性
由上面的递推式，

Ak[i][j] $A^k[i][j]$ 代表的是:从第

i $i$ 种颜色经过

k $k$ 个循环后变为第

j $j$ 种颜色的方法数。考虑循环，需要知道从

i $i$ 出发回到

i $i$ 的方案数即：

Ak[i][i] $A^{k}[i][i]$ 。

solve(k) $solve(k)$ 表示循环个数为

k $k$ 时的方案数:

s o l v e (k) = \sum i \leq m A k [i] [i]

$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$ 可以得到的方案数是

s o l v e (k) ? ? (n k)

$solve(k)*\phi (\frac{n}{k})$
枚举每个

k $k$ 然后求和最后还要除以

n $n$ (因为总的置换数为

n $n$ ),又因为有模数且模数为素数，那么就相当于乘以

nmod?2 $n^{mod-2}$ (费马小定理)。

A n s = {\sum i ∣ n s o l v e (i) ? ? (n i)} ? n p ? 2 % p

$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.rowMatrix 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 == colmemset(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;
}