[SDOI2013]项链

( ext{Problem})

非常经典的题

( ext{Analysis})

显然将思考过程分为两步

( ext{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} floor) 种颜色的涂法（(i) 的倍数）
考虑旋转和翻折的六个置换
由 (Burnside) 引理得

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

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

( ext{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) 引理得

[egin{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) 平方作为这种情况的模数
对答案也要特别处理，详情见代码

( ext{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
", ans);
	}
}