[SDOI2013]项链

\(\text{Problem}\)

非常经典的题

\(\text{Analysis}\)

显然将思考过程分为两步

\(\text{Part1}\)

求合法珍珠的种类数
设 \(f(i)\) 表示 \(\gcd = i\) 的本质不同的珍珠种类数
\(g(i)\) 表示 \(i | \gcd\) 的本质不同的珍珠种类数
那么

\[g(i) = \sum_{i|d} f(d) \]

莫比乌斯反演一下

\[f(i) = \sum_{i|d} \mu(\frac d i) g(i) \]

答案就是

\[f(1) = \sum_{i=1}^a \mu(i) g(i) \]

考虑求 \(g(i)\)
显然一个面有 \(\lfloor \frac{a}{i} \rfloor\) 种颜色的涂法（\(i\) 的倍数）
考虑旋转和翻折的六个置换
由 \(Burnside\) 引理得

\[g(i) = \frac{\lfloor \frac{a}{i} \rfloor^3+3 \lfloor \frac{a}{i} \rfloor^2 + 2 \lfloor \frac{a}{i} \rfloor}{6} \]

为保证复杂度，需要预处理 \(\mu\) 函数前缀和，并数论分块
记 \(c = f(1)\)

\(\text{Part2}\)

\(\text{Code}\)

#include<cstdio>
#define LL long long
using namespace std;

const int N = 1e7 + 5;
int flag;
LL P = 1e9 + 7, mu[N], p[105][2], totp, cnt;
LL n, a, inv6, inv61, inv62, c, ans;

int vis[N], prime[N];
inline void Euler()
{
	vis[1] = 1, mu[1] = 1;
	for(register int i = 2; i <= 1e7; i++)
	{
		if (!vis[i]) prime[++totp] = i, mu[i] = -1;
		for(register int j = 1; j <= totp && prime[j] * i <= 1e7; j++)
		{
			vis[prime[j] * i] = 1;
			if (i % prime[j] == 0) break;
			mu[prime[j] * i] = -mu[i];
		}
	}
	for(register int i = 1; i <= 1e7; i++) mu[i] += mu[i - 1];
}

inline LL fmul(LL x, LL y, LL P){return (x * y - (LL)((long double)x / P * y) * P + P) % P;}
inline LL fpow(LL x, LL y, LL P)
{
	LL res = 1;
	for(; y; y >>= 1)
	{
		if (y & 1) res = fmul(res, x, P);
		x = fmul(x, x, P);
	}
	return res;
}

inline LL calc(LL m, LL P){return (fmul(m * m % P, m, P) + m * m % P * 3 + m * 2 % P) % P;}
inline void Part1()
{
	int j;
	for(register int i = 1; i <= a; i = j + 1)
	{
		j = a / (a / i);
		c = (c + fmul(calc(a / i, P), (mu[j] - mu[i - 1] + P) % P, P)) % P;
	}
	c = fmul(c, inv6, P);
}

inline LL F(LL x){return (fpow(c - 1, x, P) + (x & 1 ? -1 : 1) * (c - 1) + P) % P;}
void dfs(int o, LL x, LL phi)
{
	if (o > cnt) return void(ans = (ans + fmul(phi, F(n / x), P)) % P);
	dfs(o + 1, x, phi);
	for(int i = 1; i < p[o][1]; i++)
	{
		x /= p[o][0], phi /= p[o][0];
		dfs(o + 1, x, phi);
	}
	dfs(o + 1, x / p[o][0], phi / (p[o][0] - 1));
}
inline void Part2()
{
	LL phi = n, x = n; cnt = 0;
	for(register int i = 1; i <= totp && prime[i] <= x; i++)
	if (x % prime[i] == 0)
	{
		p[++cnt][0] = prime[i], p[cnt][1] = 0, phi = phi - phi / prime[i];
		while (x % prime[i] == 0) x /= prime[i], ++p[cnt][1];
	}
	if (x > 1) p[++cnt][0] = x, p[cnt][1] = 1, phi = phi - phi / x;
	dfs(1, n, phi);
	if (!flag) ans = ans * fpow(n % P, P - 2, P) % P;
	else P = 1e9 + 7, ans = ans / P * fpow(n / P, P - 2, P) % P;
}

int main()
{
	int T; scanf("%d", &T);
	Euler(), inv61 = fpow(6, P - 2, P), inv62 = 833333345000000041;
	for(; T; --T)
	{
		scanf("%lld%lld", &n, &a); flag = 0;
		if (n % P == 0) flag = 1, P = P * P, inv6 = inv62;
		else inv6 = inv61;
		ans = 0, c = 0, Part1(), Part2();
		printf("%lld\n", ans);
	}
}