[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}\)

然后求由这些珍珠拼成的本质不同的项链个数
先不考虑本质不同

可以 \(dp\) 解决
设 \(F(n)\) 表示项链长度为 \(n\) 的答案
那么

\[F(n) = (c-2)F(n-1)+(c-1)F(n-2) \]

即考虑当前这个位置左右两边的珍珠种类
若不同，则这个位置有 \(c-2\) 种珍珠放法，转成子问题 \(F(n-1)\) （去除当前位，两边合并）
若相同，则这个位置有 \(c-1\) 种珍珠放法，转成子问题 \(f(n-2)\) （去除当前位，左右合并成相邻关系，但相邻珍珠相同不合法，去除一个）
初始值 \(F(0)=0，F(1)=c(c-1)\)

因为 \(n \le 10^{14}\) 很大，我们需要矩阵加速递推
却 \(TLE\) 了？！！
二阶常系数递推，当然用特征方程了！！
不难推出通项

\[F(i) = (c-1)^i + (-1)^i(c-1) \]

然后考虑本质不同（旋转后）
设旋转 \(i\) 位，\(x_j\) 位第 \(j\) 位的种类 \(x_j \equiv x_{j+i} \pmod n\)
可以视为 \(j\) 向 \(j+i \pmod n\) 连一条边，于是有了很多个环，环上珍珠种类相同
环的起点和终点相同，设环长为 \(k\)，则 \(x + ki \equiv x \pmod n\)
变形一下得

\[\frac{n}{\gcd(i,n)}|k \frac{i}{\gcd(i,n)} \]

得最小环长 \(k = \frac{n}{\gcd(i,n)}\)
所以环数量为

\[\frac{n}{k} = \gcd(i,n) \]

同环同种珍珠，不同环可不同，即对于旋转 \(i\) 的置换，有 \(F(\gcd(i,n))\) 个不动点
再次使用 \(Burnside\) 引理得

\[\begin{aligned} ans & = \frac{1}{n} \sum_{i=1}^n F(\gcd(i,n)) \\ & = \frac{1}{n} \sum_{d|n} \sum_{d|i}[\gcd(i,n)=d]F(d) \\ & = \frac{1}{n} \sum_{d|n} \sum_{d|i}[\gcd(\frac{i}{d},\frac{n}{d})=1]F(d) \\ & = \frac{1}{n} \sum_{d|n} \varphi(\frac{n}{d}) F(d) \end{aligned} \]

只要快速计算这个式子就好了
由于根号的乘上 \(\log\) 再乘上 \(T\) 的复杂度不可承受
我们需要用 \(dfs\) 的方式快速找出 \(n\) 的每个约数并根据欧拉函数的性质相应算出贡献
至此本题就差不多结束了

有一个很坑的情况就是 \(mod | n\)
这时需要特别处理，把 \(mod\) 平方作为这种情况的模数
对答案也要特别处理，详情见代码

\(\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);
	}
}