$$\begin{eqnarray}&\sum_{i=1}^{n}\sum_{j=1}^{n}ij\gcd(i,j)\\&\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}ij[\gcd(i,j)=d] \\&\sum_{d=1}^{n}d^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij[\gcd(i,j)=1] \\&\sum_{d=1}^{n}d^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}ij\sum_{x|\gcd(i,j)}\mu(x) \\&\sum_{d=1}^{n}d^3\sum_{x=1}^{\frac{n}{d}}x^2\mu(x)\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{n}{dx}}ij \\&\sum_{d=1}^{n}d^3\sum_{x=1}^{\frac{n}{d}}x^2\mu(x)(1+2+3+…\lfloor \frac{n}{xd} \rfloor )^2 \\&令s(x)=(1+x)*x/2 \\&\sum_{d=1}^{n}d^3\sum_{x=1}^{\frac{n}{d}}x^2\mu(x)s(\lfloor\frac{n}{xd}\rfloor)^2\\&令T=dx \\&\sum_{T=1}s(\lfloor\frac{n}{T}\rfloor)^2\sum_{d|T}d^3\frac{T}{d}^2\mu(\frac{T}{d})\\&\sum_{T=1}s(\lfloor\frac{n}{T}\rfloor)^2T^2\sum_{d|T}d\mu(\frac{T}{d})\\&\sum_{T=1}s(\lfloor\frac{n}{T}\rfloor)^2T^2\varphi(T)\\&令f(x)=x^2\varphi(x)\\&sum(n)=\sum_{i=1}^{n}(g*f)(i)-\sum_{i=2}^{n}g(i)sum(n/i)(杜教筛式子)\\&(g*f)(i)=i*i\sum_{d|i}\varphi(d)=i^3\\&sum(n)=\sum_{i=1}^{n} i^{3}-\sum_{i=2}^{n} i^{2} sum\left(\frac{n}{i}\right)\\&ans=\sum_{T=1}^{n} \operatorname{sum}\left(\frac{n}{T}\right)^{2} T^{2} \sum_{d|T} d \mu\left(\frac{T}{d}\right)\\&\end{eqnarray}$$
#include<bits/stdc++.h> #define int long long using namespace std; const long long INF=1LL<<31; int N=8000000; int cnt,n; long long p[8001000],inv2,inv6,ans,zhi[8001000],mod; bool he[8001000]; map<long long,long long>M; long long S(long long x){x%=mod;return x*(x+1)%mod*inv2%mod;} long long Sump(long long x){x%=mod;return x*(x+1)%mod*(x+x+1)%mod*inv6%mod;} void xxs() { he[1]=p[1]=1; for(int i=2;i<=N;i++) { if(he[i]==0) { p[i]=(i-1)%mod; zhi[++cnt]=i; } for(int j=1;j<=cnt&&i*zhi[j]<=N;j++) { he[i*zhi[j]]=true; if(i%zhi[j]==0) { p[i*zhi[j]]=1LL*p[i]*zhi[j]%mod; break; } else { p[i*zhi[j]]=1LL*p[i]*(zhi[j]-1)%mod; } } } for(int i=1;i<=N;i++)p[i]=(p[i-1]+1ll*p[i]*i%mod*i%mod)%mod; } long long SF(long long x) { if(x<=N) return p[x]; if(M.find(x)!=M.end()) return M[x]; long long ret=S(x); ret=ret*ret%mod; for(long long i=2,r;i<=x;i=r+1) { r=x/(x/i); long long tt=(Sump(r)-Sump(i-1))%mod; ret-=SF(x/i)*tt%mod; ret%=mod; } return M[x]=(ret+mod)%mod; } long long quick_pow(long long a,long long b) { long long res=1; while(b>0) { if(b&1) { res*=a; res%=mod; } a*=a; a%=mod; b>>=1; } return res; } signed main() { scanf("%lld%lld",&mod,&n); inv2=quick_pow(2,mod-2); inv6=quick_pow(6,mod-2); xxs(); for(long long i=1,r;i<=n;i=r+1) { r=n/(n/i); long long tt=S(n/i); tt=tt*tt%mod; long long gg=(SF(r)-SF(i-1))%mod; ans+=gg*tt%mod; ans%=mod; } printf("%lld\n",(ans+mod)%mod); return 0; }