题目：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+,P=1e9+;
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[]=; int mx=1e7;
for(int i=;i<=mx;i++)
{
if(!vis[i])pri[++cnt]=i,mu[i]=-;
for(int j=;j<=cnt&&(ll)i*pri[j]<=mx;j++)
{
vis[i*pri[j]]=;
if(i%pri[j]==){mu[i*pri[j]]=; break;}
mu[i*pri[j]]=-mu[i];
}
}
for(int i=;i<=mx;i++)mu[i]+=mu[i-];
}
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<)x+=mod; return x;}
ll pw(ll a,ll b)
{
ll ret=; a=a%mod;
for(;b;b>>=1ll,a=mul(a,a))if(b&)ret=mul(ret,a);
return ret;
}
ll pw2(ll a,ll b)
{
ll ret=; a=a%P;
for(;b;b>>=1ll,a=(a*a)%P)if(b&)ret=(ret*a)%P;
return ret;
}
void div(ll x)
{
num=;
for(int i=;i<=cnt&&(ll)pri[i]*pri[i]<=x;i++)
{
if(x%pri[i])continue;
p[++num]=pri[i]; pk[num]=;
while(x%pri[i]==)pk[num]++,x/=pri[i];
}
if(x>)p[++num]=x,pk[num]=;
}
ll calf(ll x)
{
ll tmp;
if(x&)tmp=upt(-m); else tmp=upt(m-);
return upt(pw(upt(m-),x)+tmp);
}
void dfs(int nw,ll d,ll phi)
{
if(nw==num+){ans=upt(ans+mul(calf(n/d),phi)%mod); return;}
dfs(nw+,d,phi);
d*=p[nw]; phi*=p[nw]-; dfs(nw+,d,phi);
for(int i=;i<=pk[nw];i++)
d*=p[nw],phi*=p[nw],dfs(nw+,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=,ans3=;
for(ll i=,j;i<=A;i=j+)
{
j=A/(A/i);
ans2=upt(ans2+mul(mul(A/i,A/i),mu[j]-mu[i-]+mod)%mod);
ans3=upt(ans3+mul(mul(mul(A/i,A/i),A/i),mu[j]-mu[i-]+mod)%mod);
}
m=upt(ans3+mul(ans2,)); m=upt(m+);
m=mul(m,pw(,(ll)P*(P-)-)%mod);//phi[mod]-1 div(n); ans=; dfs(,,); if(n%P==)ans=(ans/P*pw2(n/P,P-))%P;//P-2
else ans=(ans%P*pw2(n%P,P-))%P;//pw,mul:%mod
printf("%lld\n",ans);
}
return ;
}