「POJ2888」Magic Bracelet - Polya定理 + 欧拉函数 + 矩阵快速幂

对$n$个珠子构成的环染$m$种颜色, 并且规定一些颜色不能相邻. 旋转后相同算是同一种方案, 求本质不同的着色方案数

$n\le 10^9, m\le 10$

Link

POJ 2888

Solution

一道很好的Polya定理的题

1. Burnside引理：对于一个置换$f$,若一个染色方案$s$经过置换后不变，称$s$为$f$的不动点。将$f$的不动点数目记为$C(f)$，则可以证明等价类数目为所有$C(f)$的平均值

2. Polya定理：假设一个置换有$k$个循环，易知每个循环对应的所有位置颜色需一致，而任意两个循环之间选什么颜色互不影响。因此，如果有$m$种可选颜色，则该置换对应的不动点个数为$m^k$。所以等价类数目为：$\displaystyle \frac{\sum_{i=0}^{|F}m^{k_i}}{|F|}$

   其中$|F|$表示置换的数目，$k_i$表示第$i$个置换包含的循环个数

我们称把环旋转$\frac{k}{n}\cdot \pi$度的置换为$k$置换。显然$k\in[0, n)$，即总共只有$n$种置换，但本质不同的置换只有$n$的约数种

这是因为，$k$置换的循环长度为$\frac{n}{\gcd(k, n)}$，循环个数为$\gcd(k, n)$

想象一下在环上，一个点从$1$开始，每次往后跳$k$步，最终跳回$1$号点所需要的步数，就是循环长度

根据Polya定理化式子：

$\begin{aligned} ans&=\frac{1}{n}\cdot \sum_{i=1}^{n}m^{gcd(n,i)}\\ &=\frac{1}{n}\cdot \sum_{d|n}\cdot m^{d}\cdot\varphi(\frac{n}{d}) \end{aligned}$

枚举$\gcd(n, i) = d$，那么这样的$i$会有$\displaystyle \varphi(\frac{n}{d})$个（$i$要有$d$这个因子，且把$n$和$i$同除$d$后要互质）

式子中"$m^d$"部分即为，当$\gcd(n, i) = d$时不动点的个数。当没有限制时，这$d$个循环每个都可以任选一种颜色

而这道题有颜色相邻的限制，所以需要预处理出$f_{i, j}$表示填了$i$种颜色，第$i$种颜色为$j$的方案数。显然可以矩阵快速幂优化

一个小技巧：因为首尾也要满足条件，所以可以先多矩乘一次，然后把所有的$G_{i, i}$相加即为答案

复杂度$O(\sqrt n m^3\log n)$，跑不满

Code

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>

#define x first
#define y second
#define y1 Y1
#define y2 Y2
#define mp make_pair
#define pb push_back
#define DEBUG(x) cout << #x << " = " << x << endl;

using namespace std;

typedef long long LL;
typedef pair <int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }
template <typename T> inline T read ()
{
	T sum = 0, fl = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
	for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + ch - '0';
	return sum * fl;
}

const int Mod = 9973;
const int Maxn = 15;

int N, M, K;

namespace MATH
{
	inline void Add (int &a, int b) { if ((a += b) >= Mod) a -= Mod; }

	inline int Pow (int a, int b)
	{
		int ans = 1;
		for (int i = b; i; i >>= 1, a = (LL) a * a % Mod) if (i & 1) ans = (LL) ans * a % Mod;
		return ans;
	}

	const int maxn = 1e5 + 100;

	int prime[maxn], phi[maxn];
	bool vis[maxn];

	inline void init ()
	{
		int n = 1e5;
		phi[1] = 1;
		for (int i = 2; i <= n; ++i)
		{
			if (!vis[i]) prime[++prime[0]] = i, phi[i] = i - 1;
			for (int j = 1; j <= prime[0] && (LL) i * prime[j] <= n; ++j)
			{
				vis[i * prime[j]] = 1;
				if (!(i % prime[j]))
				{
					phi[i * prime[j]] = phi[i] * prime[j];
					break;
				}
				phi[i * prime[j]] = phi[i] * phi[prime[j]];
			}
		}
	}

	inline int get_phi (int n)
	{
		if (n <= 1e5) return phi[n];
		int ans = n;
		for (int i = 1; prime[i] * prime[i] <= n; ++i)
			if (!(n % prime[i]))
			{
				ans -= ans / prime[i];
				while (!(n % prime[i])) n /= prime[i];
			}
		if (n > 1) ans -= ans / n;
		return ans % Mod;
	}
}

using namespace MATH;

struct mat
{
	int A[Maxn][Maxn];
	inline mat () { memset (A, 0, sizeof A); }

	inline int* operator [] (const int &x) { return A[x]; }

	inline mat operator * (const mat &rhs) const
	{
		mat now;
		for (int i = 1; i <= M; ++i)
			for (int k = 1; k <= M; ++k)
				if (A[i][k])
					for (int j = 1; j <= M; ++j)
						Add (now.A[i][j], A[i][k] * rhs.A[k][j] % Mod);
		return now;
	}

	inline mat operator ^ (int b)
	{
		mat ans, a = *this;
		for (int i = 1; i <= M; ++i) ans.A[i][i] = 1;
		for (int i = b; i; i >>= 1, a = a * a) if (i & 1) ans = ans * a;
		return ans;
	}

	inline void init () 
	{
		for (int i = 1; i <= M; ++i) 
			for (int j = 1; j <= M; ++j) 
				A[i][j] = 1; 
	}

	inline void print ()
	{
		for (int i = 1; i <= M; ++i, puts(""))
			for (int j = 1; j <= M; ++j)
				cout << A[i][j] << " ";
	}
} trans;

int ans;

inline void Calc (int x)
{
	mat res = trans ^ x;
	int sum = 0;
	for (int i = 1; i <= M; ++i) Add (sum, res[i][i]);
	sum = (LL) sum * get_phi (N / x) % Mod;
	Add (ans, sum);
}

inline void Solve ()
{
	ans = 0;
	int m = sqrt (N);
	for (int i = 1; i <= m; ++i)
		if (!(N % i))
		{
			Calc (i);
			if (i * i != N) Calc (N / i);
		}
	ans = (LL) ans * Pow (N, Mod - 2) % Mod;
	printf("%d\n", ans);
}

inline void Input ()
{
	N = read<int>(), M = read<int>(), K = read<int>();
	trans.init ();
	for (int i = 1; i <= K; ++i)
	{
		int x = read<int>(), y = read<int>();
		trans[x][y] = trans[y][x] = 0;
	}
}

int main()
{

#ifdef hk_cnyali
	freopen("C.in", "r", stdin);
	freopen("C.out", "w", stdout);
#endif

	MATH :: init ();
	int T = read<int>();
	while (T--)
	{
		Input ();
		Solve ();
	}

	return 0;
}