题目：
参考了博客：
据 Narh 的想法，其实算珠子个数也可以从置换的角度，三棱柱有6种置换，循环节个数为1的有2个，个数为2的有3个，个数为3的有1个；
然后一个循环节内数字相同，于是也是那样算...
还不太懂 O(1) 快速乘...
注意模 P 和模 mod 。
代码如下：
#include#include #include #include 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;}