(ans=(sum_{i=1}^nsum_{j=1}^nijgcd(i,j))\%p)
预处理(n^{frac 2 3})后，时间复杂度为(O(n^{frac23}))
#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read(){
int x=0,f=1;char c=getchar();
while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
return f==1?x:-x;
}
const int N=5e6+4;
int mod,tot,inv6,ans,s[N],pri[N];
bitset<N>b;
unordered_map<int,int>ms;
inline int ksm(int x,int r){
int ret=1;
for(int i=0;(1ll<<i)<=r;i++){
if((r>>i)&1)ret=ret*x%mod;
x=x*x%mod;
}
return ret;
}
inline void pre(int n){
s[1]=1;
for(int i=2;i<=n;i++){
if(!b[i]){pri[++tot]=i;s[i]=i-1;}
for(int j=1;j<=tot&&pri[j]*i<=n;j++){
b[i*pri[j]]=1;
if(i%pri[j])s[i*pri[j]]=s[i]*(pri[j]-1)%mod;
else{
s[i*pri[j]]=s[i]*pri[j]%mod;
break;
}
}
}
for(int i=2;i<=n;i++)s[i]=(s[i-1]+s[i]*i%mod*i)%mod;
}
inline int sum(int x){
x%=mod;
return x*(x+1)/2%mod;
}
inline int sum2(int x){
x%=mod;
return x*(x+1)%mod*(x*2+1)%mod*inv6%mod;
}
inline int S(int n){
if(n<N)return s[n];
if(ms[n])return ms[n];
int ret=sum(n)*sum(n)%mod;
for(int l=2,r,x;l<=n;l=r+1){
r=n/(n/l);
x=sum2(r)-sum2(l-1);
ret=(ret-x*S(n/l))%mod;
}
ms[n]=ret;
return ret;
}
signed main(){
mod=read();int n=read();
pre(N-1);
inv6=ksm(6,mod-2);
for(int l=1,r,pr=0,nw;l<=n;l=r+1){
r=n/(n/l);
nw=S(r);
ans=(ans+sum(n/l)*sum(n/l)%mod*(nw-pr))%mod;
pr=nw;
}
cout<<(ans+mod)%mod;
return (0-0);
}