题目链接：Click here
Solution:
看到gcd，我们先尝试常规转化一下
[sum_{i=1}^n sum_{j=1} ^n ij \,gcd(i,j)\
sum_{d=1}^n d^3 sum_{i=1}^{lfloor{n over d}
floor} sum_{j=1} ^{lfloor{n over d}
floor} ij [gcd(i,j)=1]\
sum_{d=1}^n d^3 sum_{i=1}^{lfloor{n over d}
floor} sum_{j=1} ^{lfloor{n over d}
floor} ij sum_{t|i,t|j} mu(t)\
s(n)=sum_{i=1}^n i\
sum_{d=1}^n d^3 sum_{t=1} ^{nover d} mu(t) t^2s(lfloor frac{n}{dt}
floor )^2\
T=dt\
sum_{T=1} ^n s(lfloor frac{n}{T}
floor) ^2sum_{d|T} d^3mu({Tover d}) {T over d} ^2\
sum_{T=1} ^n s(lfloor frac{n}{T}
floor) ^2T^2sum_{d|T} dmu({Tover d}) \
]
做到这里，我们发现有点麻烦，前面可以数论分块，但后面部分线筛的时间我们是接受不了的
我们把后面那部分拿出来，看看怎么来处理
[f(T)=T^2sum_{d|T} dmu({Tover d}) \
]
后面是一个卷积的形式，于是我们有
[(id*mu)=varphi\
f(n)=n^2varphi(n)\
]
现在我们要求(f(n))的前缀和，考虑杜教筛
[S(n)=sum_{i=1}^n f(i)\
g(1)S(n)=sum_{i=1}^n (g*f)(i)-sum_{i=2}^n g(i)S(lfloor{nover i}
floor)\
]
我们让(g(n)=n^2)
[(g*f)(n)=sum_{d|n}f(d)g({nover d})\
(g*f)(n)=sum_{d|n}d^2 varphi(d){n^2 over d^2}\
(g*f)(n)=n^2sum_{d|n}varphi(d)\
]
我们知道(varphi(n))有很多美妙的性质
[sum_{d|n} varphi(d)=n\
(g*f)(n)=n^3\
g(1)S(n)=sum_{i=1}^n i^3-sum_{i=2}^n i^2S(lfloor{nover i}
floor)\
]
对于幂和，我们很多结论
[sum_{i=1}^n i^3=s(n)^2\
sum_{i=1}^n i^2={n(n+1)(2n+1)over 6}\
]
那么我们就可以利用杜教筛得到(S(n))
[sum_{T=1} ^n s(lfloor frac{n}{T}
floor) ^2f(T)\
]
那么我们只需要对上面的式子数论分块即可
Code:
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=4e6+11;
bool vis[N];
int un,cnt,mod,phi[N],p[N];
int inv2,inv6,ans;
map<int,int> vphi;
int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}
return x*f;
}
int qpow(int a,int b){
int re=1;
while(b){
if(b&1) re=re*a%mod;
b>>=1;a=a*a%mod;
}return re;
}
void prepare(){phi[1]=1;
for(int i=2;i<N;i++){
if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt&&i*p[j]<N;j++){
vis[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j]%mod;
break;
}phi[i*p[j]]=phi[i]*phi[p[j]]%mod;
}
}
for(int i=1;i<N;i++){
phi[i]=(i*i%mod)*phi[i]%mod;
phi[i]=(phi[i]+phi[i-1])%mod;
}
}
int sum(int x){x%=mod;return x*(x+1)%mod*inv2%mod;}
int sum2(int x){x%=mod;return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;}
int sum3(int x){x%=mod;return sum(x)*sum(x)%mod;}
int getS(int n){
if(n<N) return phi[n];
if(vphi[n]) return vphi[n];
int re=sum3(n);
for(int i=2,j;i<=n;i=j+1){
j=n/(n/i);
int tmp=(sum2(j)-sum2(i-1))%mod;
re=re-(tmp*getS(n/i)%mod);re=re%mod;
}return vphi[n]=(re+mod)%mod;
}
signed main(){
mod=read(),un=read();
prepare();
inv2=qpow(2,mod-2);
inv6=qpow(6,mod-2);
for(int i=1,j;i<=un;i=j+1){
j=un/(un/i);
int tmp=getS(j)-getS(i-1);
int uv=sum(un/i);uv=uv*uv%mod;
tmp=tmp%mod;ans=ans+(uv*tmp%mod);
ans=ans%mod;
}printf("%lld
",(ans+mod)%mod);
return 0;
}