恶心的取模调了一下午。。。
求$\sum\limits_{i=1}^n\sum\limits_{j=1}^nijgcd(i,j)$
枚举$gcd$：
$=\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^nijd[gcd(i,j)=d]$
枚举$d$的倍数：
$=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^3[gcd(i,j)=1]$
套路反演：
$=\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum\limits_{k|gcd(i,j)}\mu(k)$
$=\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[k|i][k|j]ij\mu(k)$
枚举$k$的倍数：
$=\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{kd}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{n}{kd}\rfloor}j$
到这里复杂度还是不够优秀，继续往下推：
设$g(n)=\sum\limits_{i=1}^ni$，即自然数前$n$项和。
令$T=kd$，则有$d|T$。
枚举$T$：
$=\sum\limits_{d=1}^nd^3\sum\limits_{T=1}^n[d|T]\mu(\dfrac{T}{d})\left(\dfrac{T}{d}\right)^2g^2(\dfrac{n}{T})$
变换求和顺序：
$=\sum\limits_{T=1}^ng^2(\dfrac{n}{T})T^2\sum\limits_{d|T}\mu(\dfrac{T}{d})d$
后面那个$\sum$其实是个卷积：
$\sum\limits_{d|T}\mu(\dfrac{T}{d})d$
$=\mu*id$
$=\mu*\varphi*1$
$=\varphi$
我们要求的就成了
$\sum\limits_{T=1}^ng^2(\dfrac{n}{T})T^2\varphi(T)$
考虑杜教筛$\sum\limits_{i=1}^ni^2\varphi(i)$。
给$i^2\varphi(i)$卷个$id^2$：
$\sum\limits_{d|n}d^2\varphi(d)\dfrac{n^2}{d^2}=n^3$
然后就能$O(n^{\frac{2}{3}})$做了。
通过查阅题解还发现了一种更$NB$的推法：
根据$n=\sum\limits_{d|n}\varphi(d)$
原式直接化为$\sum\limits_{i=1}^n\sum\limits_{j=1}^nij\sum\limits_{d|gcd(i,j)}\varphi(d)$
剩下的就都是套路了。
代码：
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#define maxn 10000001
#define inf 0x3f3f3f3f
using namespace std;
inline int read(){
int x=0,y=0;
char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')y=1;ch=getchar();}
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
return y?-x:x;
}
const int Mod = 5e6 + 7;
int mod,phi[maxn],prime[maxn>>2],cnt,inv6;
bool is[maxn];
struct hash{
int h[Mod],d[Mod],nex[Mod],cnt;
long long s[Mod];
void insert(long long x,int y){
s[++cnt]=x,x%=Mod,nex[cnt]=h[x],d[h[x]=cnt]=y;
}
int find(long long x){
for(register int i=h[x%Mod];i;i=nex[i])
if(s[i]==x)return d[i];
return mod;
}
}mphi;
int INV(int x){return x==1?1:1ll*(mod-mod/x)*INV(mod%x)%mod;}
int pow2(long long n){
if(n>=mod)n%=mod;
return n*(n+1)%mod*((n<<1)+1)%mod*inv6%mod;
}
int pow3(long long n){
if(n>=mod)n%=mod;
int ans=(n*(n+1)>>1)%mod;
return 1ll*ans*ans%mod;
}
int sphi(long long n){
if(n<maxn)return phi[n];
int ans=mphi.find(n);
if(ans!=mod)return ans;
ans=pow3(n);
for(long long l=2,r;l<=n;l=r+1){
r=n/(n/l);
(ans-=1ll*(pow2(r)-pow2(l-1))*sphi(n/l)%mod)%=mod;
}
mphi.insert(n,ans);
return ans;
}
int main(){
mod=read(),inv6=INV(6);
is[1]=phi[1]=1;
for(register int i=2;i<maxn;++i){
if(!is[i])prime[++cnt]=i,phi[i]=i-1;
for(register int j=1;j<=cnt&&i*prime[j]<maxn;++j){
is[i*prime[j]]=1;
if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
for(register int i=2;i<maxn;++i)phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i%mod)%mod;
long long n=read<long long>(),p;
int ans=0;
for(long long l=1,r;l<=n;l=r+1){
p=n/l,r=n/p,p%=mod,p=(p*(p+1)>>1)%mod;
(ans+=p*p%mod*(sphi(r)-sphi(l-1))%mod)%=mod;
}
printf("%d\n",(ans+mod)%mod);
}
