破烂衣裳

题目链接：YBT2023寒假Day15 A

题目大意

有一个 n 个点的环，有 k 种颜色，一开始都没有颜色。
每次你可以选择一个位置，把它染成 k 种颜色的其中一个，但是相邻的两个位置会变回没有颜色。
然后问你会有多少种不同的环。

思路

看到数环的个数不难想到 Burnside 引理和 Polya 定理。
发现这个对于不同位置之间颜色会相互制约，所以我们考虑用 Burnside 引理。
那置换方式就是循环同构，要求的 \(ans=\dfrac{1}{n}\sum\limits_{i=1}^nc(r_i)\)

那注意到这个今天置换会让数组出现循环节，第 \(i\) 个置换的循环节就是 \(\gcd(i,n)\)。
那我们就可以枚举作为循环节的长度 \(x(x|n)\)，那是这个循环节的个数就是 \(\varphi(n/x)\)。

那考虑对于每个循环节怎么求答案。
首先一个显然的 DP，这个一个位置可以有颜色（可以有 \(k\) 种），它的旁边就不能有颜色。
直接 \(f_{0/1,i,0/1}\) 表示一开始是否是有颜色，链的长度为 \(i\)，最后是否有颜色，直接转移即可。
最后把环给闭上。

那注意到 \(n\) 很大，循环节很大，但是转移肉眼可见的简单。
容易想到用矩阵乘法优化，然后就没了。

代码

#include<cstdio>
#define ll long long 
#define mo 1000000007

using namespace std;

const int N = 2e6 + 100;
const int sqN = 35001;
int n, k, prime[N];
ll f[2][N][2], ans;
bool np[N];

ll ksm(ll x, ll y) {
	ll re = 1;
	while (y) {
		if (y & 1) re = re * x % mo;
		x = x * x % mo; y >>= 1;
	}
	return re;
}

int gcd(int x, int y) {
	if (!y) return x;
	return gcd(y, x % y);
}

struct matrix {
	int n, m;
	ll f[2][2];
}g[51], tmp;

matrix operator *(matrix x, matrix y) {
	matrix re; for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) re.f[i][j] = 0;
	re.n = x.n; re.m = y.m;
	for (int i = 0; i < re.n; i++)
		for (int j = 0; j < re.m; j++)
			for (int k = 0; k < x.m; k++)
				(re.f[i][j] += x.f[i][k] * y.f[k][j] % mo) %= mo;
	return re;
}

ll phi(int now) {
	ll re = now;
	for (int i = 1; prime[i] * prime[i] <= now; i++) {
		if (now % prime[i] == 0) {
			re = re / prime[i] * (prime[i] - 1);
			while (now % prime[i] == 0) now /= prime[i]; 
		}
	}
	if (now != 1) re = re / now * (now - 1);
	return re;
}

ll slove(ll now) {
	if (!now) return 1;
	ll re = 0;
	tmp.n = 1; tmp.m = 2; tmp.f[0][0] = 1; tmp.f[0][1] = 0;
	for (int i = 50; i >= 0; i--) {
		if ((now >> i) & 1) tmp = tmp * g[i];
	}
	(re += tmp.f[0][1] + tmp.f[0][0]) %= mo;
	tmp.n = 1; tmp.m = 2; tmp.f[0][0] = 0; tmp.f[0][1] = k;
	for (int i = 50; i >= 0; i--) {
		if ((now >> i) & 1) tmp = tmp * g[i];
	}
	(re += tmp.f[0][0]) %= mo;
	return re;
}

ll work(int now) {
	ll re = phi(n / now);
	re = re * slove(now - 1) % mo;
	return re;
}

int main() {
	freopen("iwill.in", "r", stdin);
	freopen("iwill.out", "w", stdout);
	
	for (int i = 2; i < sqN; i++) {
		if (!np[i]) {
			prime[++prime[0]] = i;
		}
		for (int j = 1; j <= prime[0] && i * prime[j] < sqN; j++) {
			np[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
		}
	}
	
	scanf("%d %d", &n, &k);
	
	if (n <= 2000000) {
		f[0][1][0] = 1; f[1][1][1] = k;
		for (int i = 2; i <= n; i++) {
			for (int j = 0; j <= 1; j++) {
				f[j][i][0] = (f[j][i - 1][0] + f[j][i - 1][1]) % mo;
				f[j][i][1] = f[j][i - 1][0] * k % mo;
			}
		}
		
		ll ans = 0;
		for (int i = 1; i <= n; i++) {
			int tmp = gcd(i, n);
			(ans += f[0][tmp][1] + f[1][tmp][0] + f[0][tmp][0]) %= mo;
		}
		ans = ans * ksm(n, mo - 2) % mo;
		printf("%lld", ans);
	}
	else {
		g[0].n = g[0].m = 2;
		g[0].f[0][0] = 1; g[0].f[1][0] = 1;
		g[0].f[0][1] = k;
		for (int i = 1; i <= 50; i++) g[i] = g[i - 1] * g[i - 1];
		
		for (int i = 1; i * i <= n; i++)
			if (n % i == 0) {
				(ans += work(i)) %= mo; if (n / i != i) (ans += work(n / i)) %= mo;
			}
		printf("%lld", ans * ksm(n, mo - 2) % mo);
	}
	
	return 0;
}