SDOI2013 项链

首先我们可以求出来满足条件的珠子的总数 $m$。有以下结论：

• 对于一个大小为 $n$ 的环，用 $m$ 种颜色给它染色，要求环上相邻两个点的颜色不同，那么方案数为：$(m-1)^n+(-1)^{n}(m-1)$。

我是硬推出来这个式子的2333

具体来说，设 $f(n)$ 表示大小为 $n$ 的环的答案，考虑容斥：

• 先不管第 $1$ 个元素和第 $n$ 个元素之间的限制，那么可以发现答案就是 $m\times (m-1)^{n-1}$。这是因为第 $1$ 个元素有 $m$ 种选择，后面的元素都只有 $m-1$ 种选择。
• 再去掉第一个元素和第 $n$ 个元素相同时的答案，此时可以把这两个元素看成一个，那么方案数就是 $f(n-1)$。
• 但需要特判 $n\le 2$ 的情况：$f(1)=m,f(2)=m\times(m-1)$。

我们有 $f(n)=m\times (m-1)^{n-1}-f(n-1)$，其中 $n\ge 3$。展开得到

$$f(n)=m\times \sum_{i=1}^{n-1}(-1)^{n+1-i}(m-1)^i=m\times\left(\sum_{i=0}^{n-1}(-1)^{n+1-i}(m-1)^i\right)-m$$

把后面个求和里面的 $(-1)^{n+1}$ 提出来得到

$$f(n)=m\times (-1)^{n+1}\times \sum_{i=0}^{n-1}(-1)^{-i}(m-1)^i=m\times (-1)^{n+1}\times \sum_{i=0}^{n-1}(1-m)^i-(-1)^{n+1}m$$

由等比数列求和公式有

$$\begin{aligned} f(n)&=m\times (-1)^{n+1}\times\dfrac{1-(1-m)^{n}}{1-(1-m)}=(-1)^{n+1}\times (1-(1-m)^n)-(-1)^{n+1}m\\ &=(-1)^{n+1}+(m-1)^n-(-1)^{n+1}m=(-1)^{n}(m-1)+(m-1)^n \end{aligned}$$

这就证明了结论。

考虑怎么算出 $m$。再用一次 $\text{Burnside}$ 引理，我们设

$$S_1=1\\ S_2=\sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1]\\ S_3=\sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1]$$

这些都可以 $O(\sqrt{a})$ 求解。那么 $m=\frac{1}{6}(2S_1+3S_2+S_3)$。

本题中置换群由 $n$ 种旋转构成。

如果转了 $k$ 格，那么会形成 $\gcd(k,n)$ 个置换环，答案就是 $(m-1)^{\gcd(k,n)}+(-1)^{\gcd(k,n)+1}$。显然可以 $O(\sqrt{n}\log n)$ 求出来。

看上去 $n=10^{14}$ 非常不可过，但是你发现这个算法的复杂度实际上是 $O(d(n)\log n)$，我们都知道 $d(n)\sim O(\sqrt[3]{n})$，于是这个算法也就可以通过本题了。

然后我们发现 $\text{Burnside}$ 引理算完结束之后需要除以 $n$，但这个题 $n$ 可能是 $p=10^9+7$ 的倍数导致逆元不存在

我们可以在 $\bmod p^2$ 意义下进行计算，这样一来如果我们算出来那个 $\sum_{g}D(g)=A$，那么如果 $p\mid n$ 我们就计算 $(A/p)\times (n/p)^{-1}\bmod p$ 就可以了。

#include<bits/stdc++.h>

#define int long long

using namespace std;

int read(){
	int x=0,f=1;char c=getchar();
	for(;(!(c>='0'&&c<='9'));c=getchar())if(c=='-')f=-1;
	for(;(c>='0'&&c<='9');c=getchar())x=x*10+c-'0';
	return x*f;
}

const int N=1e7;
int pri[N+5],cnt=0,mu[N+5];
bool vis[N+5];

void init(){
	for(int i=2;i<=N;i++){
		if(!vis[i])pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pri[j]<=N;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
			mu[i*pri[j]]=-mu[i];
		}
	}
	mu[1]=1;for(int i=1;i<=N;i++)mu[i]+=mu[i-1];
}

const int mod=1e9+7;

int mul(int x,int y,int p){
	__int128 xx=x,yy=y,pp=p;
	return (int)(xx*yy%pp);
}

const int MOD=mod*mod;

int ksm(int x,int y,int p=MOD){
	int ans=1;
	for(int i=y;i;i>>=1,x=mul(x,x,p))if(i&1)ans=mul(ans,x,p);
	return ans;
}

int inv(int x,int p=MOD){
	return ksm(x,p-2,p)%p;
}

const int iv=833333345000000041;

int calc(int a){
	int res=0;
	for(int l=1,r=0;l<=a;l=r+1){
		r=a/(a/l);int w=a/l;
//		cout<<l<<" "<<r<<endl;
		res=(res+mul((mul(w,mul(w,w,MOD),MOD)+mul(3,mul(w,w,MOD),MOD)),(mu[r]-mu[l-1]+MOD)%MOD,MOD))%MOD;
//		cout<<res<<endl;
	}
	res+=2;
//	cout<<"res = "<<res<<endl;
	return mul(res,iv,MOD);
}

int p[30],c[30],r[30],k,n,m,ans=0;
void dfs(int now){
	if(now==k+1){
		int phi=1,d=1;
		for(int i=1;i<=k;i++){
			if(r[i]>0)phi=phi*ksm(p[i],r[i]-1)*(p[i]-1);
			d=d*ksm(p[i],c[i]-r[i]);
		}
//		cout<<d<<" | "<<phi<<endl;
		ans=(ans+(mul((ksm(m-1,d)%MOD+mul(d%2==0?1:-1,m-1,MOD)+MOD)%MOD,phi%MOD,MOD)))%MOD;
		return ;
	}
	for(int i=0;i<=c[now];i++)r[now]=i,dfs(now+1),r[now]=0;
}

void solve(){
	k=0;ans=0;
	n=read();int a=read();
	m=calc(a);int div=n;
//	cout<<m<<endl;
	for(int i=2;i*i<=n;i++){
		if(n%i==0){
			p[++k]=i,c[k]=0;
			while(n%i==0)n/=i,c[k]++;
		}
	}
	if(n>1)p[++k]=n,c[k]=1;
//	for(int i=1;i<=k;i++)cout<<p[i]<<" ^ "<<c[i]<<endl;
	dfs(1);
	if(div%mod!=0)cout<<ans%mod*inv(div%mod,mod)%mod<<endl;
	else cout<<(ans/mod)*inv(div/mod,mod)%mod<<endl;
}

signed main(void){

#ifndef ONLINE_JUDGE
	freopen("in.in","r",stdin);
#endif

	init();
	int tt=read();
	while(tt--)solve();

	return 0;
}