題目：https://www.lydsy.com/JudgeOnline/problem.php?id=3202

參考了部落格：

https://www.cnblogs.com/zhoushuyu/p/9657640.html

https://www.cnblogs.com/DUXT/p/5957944.html?utm_source=itdadao&utm_medium=referral

https://blog.csdn.net/Maxwei_wzj/article/details/83184110

https://blog.csdn.net/a_crazy_czy/article/details/50688526

據 Narh 的想法，其實算珠子個數也可以從置換的角度，三稜柱有6種置換，迴圈節個數為1的有2個，個數為2的有3個，個數為3的有1個；

然後一個迴圈節內數字相同，於是也是那樣算...

還不太懂 O(1) 快速乘...

注意模 P 和模 mod 。

程式碼如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=1e7+5,P=1e9+7;
ll const mod=(ll)P*P;int A,pri[xn],cnt,mu[xn],pk[xn],num;
ll n,m,ans,p[xn];
 
//p!!
bool vis[xn];
void init()
{
  mu[1]=1; int mx=1e7;
  for(int i=2;i<=mx;i++)
    {
      if(!vis[i])pri[++cnt]=i,mu[i]=-1;
      for(int j=1;j<=cnt&&(ll)i*pri[j]<=mx;j++)
    {
      vis[i*pri[j]]=1;
      if(i%pri[j]==0){mu[i*pri[j]]=0; break;}
      mu[i*pri[j]]=-mu[i];
    }
    }
  
 
for(int i=2;i<=mx;i++)mu[i]+=mu[i-1];
}
ll mul(ll a,ll b){return (a*b-(ll)(((long double)a*b+0.5)/(long double)mod)*mod+mod)%mod;}
ll upt(ll x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
ll pw(ll a,ll b)
{
  ll ret=1; a=a%mod;
  for(;b;b>>=1ll,a=mul(a,a))if(b&1)ret=mul(ret,a);
  return ret;
}
ll pw2(ll a,ll b)
{
  ll ret=1; a=a%P;
  for(;b;b>>=1ll,a=(a*a)%P)if(b&1)ret=(ret*a)%P;
  return ret;
}
void div(ll x)
{
  num=0;
  for(int i=1;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
    {
      if(x%pri[i])continue;
      p[++num]=pri[i]; pk[num]=0;
      while(x%pri[i]==0)pk[num]++,x/=pri[i];
    }
  if(x>1)p[++num]=x,pk[num]=1;
}
ll calf(ll x)
{
  ll tmp;
  if(x&1)tmp=upt(1-m); else tmp=upt(m-1);
  return upt(pw(upt(m-1),x)+tmp);
}
void dfs(int nw,ll d,ll phi)
{
  if(nw==num+1){ans=upt(ans+mul(calf(n/d),phi)%mod); return;}
  dfs(nw+1,d,phi);
  d*=p[nw]; phi*=p[nw]-1; dfs(nw+1,d,phi);
  for(int i=2;i<=pk[nw];i++)
    d*=p[nw],phi*=p[nw],dfs(nw+1,d,phi);
}
int main()
{
  int T; init();
  scanf("%d",&T);
  while(T--)
    {
      scanf("%lld%d",&n,&A);
      //if(n%P==0)mod=(ll)P*P; else mod=P;//??
      ll ans2=0,ans3=0;
      for(ll i=1,j;i<=A;i=j+1)
    {
      j=A/(A/i);
      ans2=upt(ans2+mul(mul(A/i,A/i),mu[j]-mu[i-1]+mod)%mod);
      ans3=upt(ans3+mul(mul(mul(A/i,A/i),A/i),mu[j]-mu[i-1]+mod)%mod);
    }
      m=upt(ans3+mul(ans2,3)); m=upt(m+2);
      m=mul(m,pw(6,(ll)P*(P-1)-1)%mod);//phi[mod]-1

      div(n); ans=0; dfs(1,1,1);
      
      if(n%P==0)ans=(ans/P*pw2(n/P,P-2))%P;//P-2
      else ans=(ans%P*pw2(n%P,P-2))%P;//pw,mul:%mod
      printf("%lld\n",ans);
    }
  return 0;
}