Link
先考虑有多少种珠子，这个旋转/翻转同构的要求实际上就是无序。
所以我们可以先求出有序的方案数，然后再除以(3!)得到无序的方案数。
设(s1,s2,s3)分别表示无限制情况下(1,2,3)元组的方案数，那么这部分的答案就是(k=frac{s3+3s2+s1}{3!})。
不难得到(s1=1,s2=sumlimits_{d=1}^amu(d)lfloorfrac ad
floor^2,s3=sumlimits_{d=1}^amu(d)lfloorfrac ad
floor^3)，可以线性筛(mu)之后数论分块求。
然后考虑用Burnside引理求出环的方案数，设(f(i))表示旋转(i)个单位的不动点数，那么答案为(frac1n(f imesvarphi)(n))。
不难发现若(d|n)，则(f(d))就是长度为(d)的序列的个数，我们可以得到其递推式(f(n)=(m-2)f(n-1)+(m-1)f(n-2))，其中(f(1)=0,f(2)=m^2-m)。
不难解得(f(n)=(m-1)^n+(-1)^n(m-1))。
那么我们dfs枚举(n)的约数同时维护当前的(varphi)就可以(O(d(n)log n))地计算((f imesvarphi)(n))了。
最后要除以(n)，由于可能(P|n)因此我们在计算过程中改为对(P^2)取模最后特判一下即可。
#include<cstdio>
#include<vector>
#include<utility>
using i64=long long;
using pi=std::pair<i64,int>;
const int N=10000007;const i64 p=1000000007,P=p*p,i6=833333345000000041;
int tot,pr[N],mu[N];i64 n,m,ans;std::vector<pi>fac;
i64 read(){i64 x;scanf("%lld",&x);return x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
i64 mul(i64 a,i64 b){i64 c=0;for(;b;b>>=1,inc(a,a))if(b&1)inc(c,a);return c;}
i64 pow(i64 a,i64 b){i64 c=1;for(;b;b>>=1,a=mul(a,a))if(b&1)c=mul(c,a);return c;}
i64 Pow(i64 a,i64 b){i64 c=1;for(;b;b>>=1,a=a*a%p)if(b&1)c=c*a%p;return c;}
void Sieve(int n)
{
static int np[N];mu[1]=1;
for(int i=2;i<=n;++i)
{
if(!np[i]) pr[++tot]=i,mu[i]=-1;
for(int j=1;i*pr[j]<=n&&j<=tot;++j)
{
np[i*pr[j]]=1;
if(i%pr[j]) mu[i*pr[j]]=-mu[i]; else break;
}
}
for(int i=2;i<=n;++i) mu[i]+=mu[i-1];
}
i64 calc(i64 n)
{
i64 a=0,b=0;
for(i64 l=1,r;l<=n;l=r+1) r=n/(n/l),inc(a,mul(pow(n/l,3),P+mu[r]-mu[l-1])),inc(b,mul(pow(n/l,2),P+mu[r]-mu[l-1]));
return mul(i6,a+3*b+2);
}
void Factor(i64 x)
{
for(int i=1,cnt;i<=tot&&1ll*pr[i]*pr[i]<=x;++i)
{
for(cnt=0;!(x%pr[i]);x/=pr[i],++cnt);
if(cnt) fac.emplace_back(pr[i],cnt);
}
if(x^1) fac.emplace_back(x,1);
}
i64 f(i64 n){return ((n&1? P+1-m:m-1)+pow(m-1,n))%P;}
void dfs(int i,i64 now,i64 phi)
{
if(i==(int)fac.size()) return inc(ans,mul(phi,f(n/now)));
dfs(i+1,now,phi),now*=fac[i].first,phi*=fac[i].first-1,dfs(i+1,now,phi);
for(int j=2;j<=fac[i].second;++j) now*=fac[i].first,phi*=fac[i].first,dfs(i+1,now,phi);
}
void solve(){n=read(),m=calc(read()),fac.clear(),ans=0,Factor(n),dfs(0,1,1),printf("%lld
",n%p? ans%p*Pow(n%p,p-2)%p:ans/p*Pow(n/p,p-2)%p);}
int main()
{
Sieve(10000000);
for(int t=read();t;--t) solve();
}