设 (m) 为有多少种不同的珠子，先用 (Burnside) 引理求 (m)，发现对于珠子的置换有旋转和翻转，得其置换群有 (|G|=6)，对于有 (3) 个不同数字的珠子，得其稳定化子 (|G_x|=1)，对于有 (2) 个不同数字的珠子，得其稳定化子 (|G_x|=2)，对于数字都相同的珠子，得其稳定化子 (|G_x|=6)。设这三类珠子的个数为 (S_3,S_2,S_1)，得：
[large m=frac{1}{|G|}sum_{fin G}c(f)=frac{1}{6}left( S_3+2S_2+6S_1
ight)
]
注意到：
[largeegin{aligned}
sum_{i=1}^asum_{j=1}^asum_{k=1}^a[gcd(i,j,k)=1]&=S_3+S_2+S_1\
sum_{i=1}^asum_{j=1}^a[gcd(i,j)=1]&=frac{1}{3}S_2+S_1\
S_1&=1
end{aligned}
]
因此得：
[large m=frac{1}{6}left( sum_{i=1}^aleftlfloorfrac{a}{i}
ight
floor^3mu(i)+3sum_{i=1}^aleftlfloorfrac{a}{i}
ight
floor^2mu(i)+2
ight)
]
考虑求环的方案数，设 (f_n) 表示 (n) 个点的环染色的方案数，即长为 (n) 的序列染色且首尾不同的方案数，环的置换只有旋转，由 (Burnside) 引理得方案数为：
[large frac{1}{n}sum_{i=1}^nf_{gcd(i,n)}=frac{1}{n}sum_{dmid n}f_dsum_{i=1}^n[gcd(i,n)=d]=frac{1}{n}sum_{dmid n}f_dvarphileft(frac{n}{d}
ight)
]
不难得到 (f_n) 的递推式：
[large f_n=(m-2)f_{n-1}+(m-1)f_{n-2},f_{1}=0,f_{2}=m(m-1)
]
考虑第 (1) 个位置和第 (n-1) 个位置颜色是否相同即可。用特征根求通项公式得：
[largeegin{aligned}
x^2&-(m-2)x-(m-1)=0\
x_1&=m-1,x_2=-1\
f_n&=c_1x_1^n+c_2x_2^n\
f_1&=c_1(m-1)-c_2=0\
f_2&=c_1(m-1)^2+c_2=m(m-1)\
f_n&=(m-1)^n+(m-1)(-1)^n\
end{aligned}
]
然后就能快速幂计算 (f_n) 了。
因为 (n) 达到 (10^{14})，所以求 (varphi(n)) 的时候应该先质因数分解 (n)，然后 (dfs) 每一个因数，搜索的过程中求出 (varphi(n))。
同样是因为 (n) 很大，可能会出现 (p mid n) 的情况，这样的话，最后除 (n) 时可能不存在逆元。于是统计答案的时候对 (p^2) 取模，因为除 (n) 前答案一定是 (n) 的倍数，那么其也一定是 (p) 的倍数，于是就直接除 (p)，然后再除 (frac{n}{p})，因为 (n< p^2)，所以 (frac{n}{p}) 在模 (p) 意义下是存在逆元的。
#include<bits/stdc++.h>
#define maxn 10000010
#define all 10000000
#define p 1000000007
#define inv6 (ll)833333345000000041
#define P ((ll)1000000007*1000000007)
using namespace std;
typedef long long ll;
template<typename T> inline void read(T &x)
{
x=0;char c=getchar();bool flag=false;
while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
if(flag)x=-x;
}
ll T,n,m,tot,cnt,ans;
ll pri[maxn],mu[maxn],v[maxn],k[maxn];
bool tag[maxn];
void init()
{
mu[1]=1;
for(int i=2;i<=all;++i)
{
if(!tag[i]) pri[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot;++j)
{
int k=i*pri[j];
if(k>all) break;
tag[k]=true;
if(i%pri[j]) mu[k]=-mu[i];
else break;
}
}
for(int i=1;i<=all;++i) mu[i]+=mu[i-1];
}
ll mul(ll x,ll y)
{
ll c=(long double)x*y/P+0.5;
c=x*y-c*P;
return c<0?c+P:c;
}
ll qp(ll x,ll y)
{
ll v=1;
while(y)
{
if(y&1) v=mul(v,x);
x=mul(x,x),y>>=1;
}
return v;
}
ll inv(ll x)
{
ll v=1,y=p-2;
while(y)
{
if(y&1) v=v*x%p;
x=x*x%p,y>>=1;
}
return v;
}
void get(ll n)
{
cnt=0;
for(int i=1;i<=tot&&(ll)pri[i]*pri[i]<=n;++i)
{
if(n%pri[i]) continue;
v[++cnt]=pri[i],k[cnt]=0;
while(n%pri[i]==0) n/=pri[i],k[cnt]++;
}
if(n!=1) v[++cnt]=n,k[cnt]=1;
}
void work(ll n)
{
ll v1=0,v2=0;
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
v1=(v1+mul(mul(mul(n/l,n/l),n/l),(mu[r]-mu[l-1]+P)%P))%P;
v2=(v2+mul(mul(n/l,n/l),(mu[r]-mu[l-1]+P)%P))%P;
}
m=mul((v1+mul(v2,3)+2)%P,inv6);
}
ll calc(ll n)
{
return (qp(m-1+P,n)+((n&1)?1-m+P:m-1+P)%P)%P;
}
void dfs(int x,ll d,ll phi)
{
if(x==cnt+1)
{
ans=(ans+mul(phi,calc(n/d)))%P;
return;
}
dfs(x+1,d,phi),d*=v[x],phi*=v[x]-1,dfs(x+1,d,phi);
for(int i=2;i<=k[x];++i) d*=v[x],phi*=v[x],dfs(x+1,d,phi);
}
int main()
{
init(),read(T);
while(T--)
{
read(n),read(m),get(n),work(m),ans=0,dfs(1,1,1);
if(n%p) ans=ans%p*inv(n%p)%p;
else ans=ans/p*inv(n/p)%p;
printf("%lld
",ans);
}
return 0;
}