https://www.luogu.com.cn/problem/P3768
推式子+杜教筛+整除分块，感觉往 (varphi) 的方向推比往 (mu) 上推要好
[egin{aligned}sum_{i=1}^nsum_{j=1}^n ijgcd(i,j)
&=sum_{i=1}^nsum_{j=1}^n ijsum_{k|gcd(i,j)} varphi(k)\
&=sum_{k}varphi(k)sum_{k|i}^nsum_{k|j}^nij\
&=sum_{k}k^2varphi(k)sum_{i=1}^{lfloorfrac{n}{k}
floor}isum_{j=1}^{lfloorfrac{n}{j}
floor}j\
&=sum_{k}k^2varphi(k)left(sum_{i=1}^{lfloorfrac{n}{k}
floor}i
ight)^2\
end{aligned}
]
然后这个 (left(sum_{i=1}^{lfloorfrac{n}{k}
floor}i
ight)^2) 可以整除分块并 (O(1)) 求，外面再套一个 (k^2varphi(k)) 用杜教筛算就行了
至于怎么杜教筛，设 (f=varphi cdot operatorname{Id}^2,g=operatorname{Id}^2)，那么就有 (f*g=sum_{d|n} varphi(d)d^2left(frac{n}{d}
ight)^2=n^2sum_{d|n} varphi(d)=n^3)
再设 (s) 是 (f) 的前缀和，根据杜教筛的套路式子，就有：
[egin{aligned}f(1)s(n)
&=sum_{i=1}^ng(i)s(lfloorfrac{n}{i}
floor)-sum_{i=2}^ng(i)s(lfloorfrac{n}{i}
floor)\
&=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)s(lfloorfrac{n}{i}
floor)\
end{aligned}
]
由于 (g(i)) 以及 (sum_{i=1}^n(f*g)(i)) 都是很好求的，所以就能直接杜教筛了
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<map>
#define reg register
#define LL_INF (long long)(0x3f3f3f3f3f3f3f3f)
#define INT_INF (int)(0x3f3f3f3f)
inline long long read(){
register long long x=0;register int y=1;
register char c=std::getchar();
while(c<'0'||c>'9'){if(c=='-') y=0;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+(c^48);c=getchar();}
return y?x:-x;
}
#define maxn 5000000
long long n;
long long mod,inv4,inv6;
std::map<long long,long long>map;
int notprime[maxn+6],prime[maxn+6];
long long phi[maxn+6];
inline long long power(long long a,long long b){
long long ans=1;
while(b){
if(b&1) ans=ans*a%mod;
a=a*a%mod;b>>=1;
}
return ans;
}
inline void pre(int n){
phi[1]=1;
for(reg int i=2;i<=n;i++){
if(!notprime[i]) prime[++prime[0]]=i,phi[i]=i-1;
for(reg int x,j=1;i*prime[j]<=n&&j<=prime[0];j++){
x=i*prime[j];notprime[x]=1;
if(i%prime[j]) phi[x]=phi[i]*(prime[j]-1);
else{
phi[x]=phi[i]*prime[j];break;
}
}
phi[i]=phi[i]*i%mod*i%mod;
}
for(reg int i=2;i<=n;i++) phi[i]=(phi[i]+phi[i-1])%mod;
}
inline long long S2(long long n){
n%=mod;
return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
inline long long S3(long long n){
n%=mod;
return n*n%mod*(n+1)%mod*(n+1)%mod*inv4%mod;
}
long long get(long long n){
if(n<=maxn) return phi[n];
if(map.find(n)!=map.end()) return map[n];
long long ans=S3(n);
for(reg long long l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans=(ans-get(n/l)*(S2(r)-S2(l-1))%mod+mod)%mod;
}
map[n]=ans;
return ans;
}
inline long long work(long long n){
long long ans=0;
for(reg long long l=1,r;l<=n;l=r+1){
r=n/(n/l);
ans+=S3(n/l)*(get(r)-get(l-1)+mod)%mod;
ans%=mod;
}
return ans;
}
int main(){
mod=read();n=read();
pre(maxn);
inv4=power(4,mod-2);inv6=power(6,mod-2);
printf("%lld
",work(n));
return 0;
}