【ybt金牌导航8-5-7】魔法手镯（Burnside 引理）（矩阵乘法优化DP）

魔法手镯

题目链接：ybt金牌导航8-5-7

题目大意

给你一个长度为 n 的珠子环，然后规定一些珠子不能相邻。
然后有多少种本质不同的珠子环。
本质相同时指它可以通过旋转变成一样的。

思路

看到循环置换，和颜色，不难想到 Polya 定理。

但是它有限制条件，那我们就考虑用 Burnside 引理。
然后你发现它的限制条件你可以按长度 DP 下去：
\(f_{i,j}\) 为一开始的颜色为 \(i\)，最后的颜色为 \(j\) 的方案数。

然后你发现这个 DP 可以用矩阵乘法优化。

然后你就直接枚举 \(\gcd(i,n)\) 的结果，每个都求出来结果乘上 \(\phi\) 再最后除 \(n\) 就是答案了。

代码

#include<cstdio>
#define mo 9973

using namespace std;

struct matrix {
	int n, m;
	int a[11][11];
}a, one, A;
int T, n, m, k, x, y, ans;

int jia(int x, int y) {
	return x + y >= mo ? x + y - mo : x + y;
}

int mul(int x, int y) {
	int re = x * y;
	if (re > mo) re = re - re / mo * mo;
	return re;
}

matrix operator *(matrix x, matrix y) {//矩阵乘法
	matrix re;
	re.n = x.n; re.m = y.m;
	for (int i = 1; i <= re.n; i++)
		for (int j = 1; j <= re.m; j++)
			re.a[i][j] = 0;
	for (int k = 1; k <= x.m; k++)
		for (int i = 1; i <= re.n; i++)
			for (int j = 1; j <= re.m; j++)
				re.a[i][j] = jia(re.a[i][j], mul(x.a[i][k], y.a[k][j]));
	return re;
}

matrix jzksm(matrix x, int y) {
	matrix re = one;
	while (y) {
		if (y & 1) re = re * x;
		x = x * x;
		y >>= 1;
	} 
	return re;
}

int ksm(int x, int y) {
	int re = 1; x = mul(x, 1);
	while (y) {
		if (y & 1) re = mul(re, x);
		x = mul(x, x);
		y >>= 1;
	}
	return re;
}

int phi(int n) {//求 phi
	int re = n;
	for (int i = 2; i * i <= n; i++) {
		if (n % i == 0) {
			re = re / i * (i - 1);
			while (n % i == 0) n /= i;
		}
	}
	if (n > 1) re = re / n * (n - 1);
	return mul(re, 1); 
}

int work(int d) {
	A = jzksm(a, n / d - 1);
	int re = 0;
	for (int i = 1; i <= m; i++)//直接上 DP 的结果
		for (int j = 1; j <= m; j++)
			if (a.a[i][j])//记得首尾要连接，所以要首尾能在一起
				re = jia(re, A.a[i][j]);
	ans = jia(ans, mul(re, phi(d)));//直接乘上取这个值的可能性（phi）
}

int main() {
	for (int i = 1; i <= 10; i++) one.a[i][i] = 1;
	
	scanf("%d", &T);
	while (T--) {
		scanf("%d %d %d", &n, &m, &k);
		
		one.n = one.m = m;
		a.n = a.m = m;
		for (int i = 1; i <= m; i++)
			for (int j = 1; j <= m; j++)
				a.a[i][j] = 1;
		for (int i = 1; i <= k; i++) {
			scanf("%d %d", &x, &y);
			a.a[x][y] = a.a[y][x] = 0;
		}
		
		ans = 0;
		for (int i = 1; i * i <= n; i++)
			if (n % i == 0) {//枚举 gcd(i,n) 的取值
				work(i); if (n / i != i) work(n / i);
			}
		
		ans = ans * ksm(n, mo - 2) % mo;
		printf("%d\n", ans);
	}
	
	return 0;
}