一、题目
二、解法
直接推式子，这几步你要有点莫比乌斯反演基础才行：
[sum_{i=1}^nsum_{j=1}^nijgcd(i,j)
]
[sum_{d=1}^nd^3sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)=1]
]
设 (sum(n)=sum_{i=1}^ni) ，我们使用莫比乌斯反演：
[sum_{d=1}^nd^3sum_{x=1}^{n/d}mu(x)x^2sum(frac{n}{dx})^2
]
如果我们推到这一步就开始做的话，后面的东西是不好处理的，令 (T=id) ，可以把 (sum) 换到前面去：
[sum_{T=1}^nsum(frac{n}{T})^2T^2sum_{d|T}mu(frac{T}{d})d
]
因为基本性质：(mu*id=varphi) （这个乘是迪利克雷卷积的意思），所以可以把后面那部分换掉
[sum_{T=1}^nsum(frac{n}{T})^2T^2varphi(T)
]
然后就可以整除分块了，但是我们要知道 (f(x)=x^2varphi(x)) 的前缀和怎么算，不难发现这是一个杜教筛的问题。令辅助函数 (g(x)=x^2)，那么 ((f*g)(n)=sum_{x|n}x^2varphi(x)(frac{n}{x})^2=n^2sum_{x|n}varphi(x)=n^3)
然后 (g) 的前缀和也是很好算的，所以直接上杜教筛就可以了，给出一些基本的柿子（其实可以直接插值）：
[sum_{i=1}^ni^2=frac{n(n+1)(2n+1)}{6}
]
[sum_{i=1}^ni^3=(frac{n(n+1)}{2})^2
]
#include <cstdio>
#include <iostream>
#include <map>
using namespace std;
const int N = 8000000;
const int M = N+5;
#define int long long
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int MOD,n,cnt,inv,ans,p[M],phi[M];map<int,int> sum;
void init(int n)
{
phi[1]=1;
for(int i=2;i<=n;i++)
{
if(!phi[i])
{
phi[i]=i-1;
p[++cnt]=i;
}
for(int j=1;j<=cnt && i*p[j]<=n;j++)
{
if(i%p[j]==0)
{
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
for(int i=2;i<=n;i++)
phi[i]=(phi[i]*i%MOD*i+phi[i-1])%MOD;
}
int zy(int x)
{
x%=MOD;
return x*(x+1)/2%MOD;
}
int jzm(int x)
{
x%=MOD;
return x*(x+1)%MOD*(2*x+1)%MOD*inv%MOD;
}
int get(int n)
{
if(n<=N) return phi[n];
if(sum[n]) return sum[n];
int ans=zy(n)*zy(n)%MOD;
for(int l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans-(jzm(r)-jzm(l-1))*get(n/l))%MOD;
}
return sum[n]=ans;
}
int qkpow(int a,int b)
{
int r=1;
while(b>0)
{
if(b&1) r=r*a%MOD;
a=a*a%MOD;
b>>=1;
}
return r;
}
signed main()
{
MOD=read();n=read();
init(N);inv=qkpow(6,MOD-2);
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+zy(n/l)*zy(n/l)%MOD*(get(r)-get(l-1)))%MOD;
}
printf("%lld
",(ans+MOD)%MOD);
}