Luogu

分析

$$
\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\
=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{d|i,d|j}\varphi(d)\\
=&\sum_{d=1}^nd^2\varphi(d)\left(\sum_{i=1}^{n/d}i\right)^2
\end{aligned}\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\
=&\sum_{i=1}^n\sum_{j=1}^nij\sum_{d|i,d|j}\varphi(d)\\
=&\sum_{d=1}^nd^2\varphi(d)\left(\sum_{i=1}^{n/d}i\right)^2
\end{aligned}
$$

显然可以数论分块，我们只需要求 $\mathbf{f}=\mathbf{id}^2\cdot\mathbf{\varphi}$ 的前缀和即可。

考虑杜教筛，根据一些奇妙的方法可以找到 $\mathbf{g}=\mathbf{id}^2$，满足
$$
\mathbf{f}*\mathbf{g}=\mathbf{id^2}\cdot(\mathbf{\varphi*\mathbf{I}})=\mathbf{id}^3
$$
这样就很好杜教筛了。

代码

// ====================================
//   author: M_sea
//   website: https://m-sea-blog.com/
// ====================================
#include <bits/stdc++.h>
#define file(x) freopen(#x".in","r",stdin); freopen(#x".out","w",stdout)
using namespace std;
typedef long long ll;

ll read() {
    ll X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=10000000+10;

ll n; int mod,inv2,inv6;
int qpow(int a,int b) { int c=1;
    for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) c=1ll*c*a%mod;
    return c;
}

int p[N],v[N],cnt=0,f[N];
void sieve(int n) {
    f[1]=1;
    for (int i=2;i<=n;++i) {
        if (!v[i]) p[++cnt]=i,f[i]=1ll*i*i%mod*(i-1)%mod;
        for (int j=1;j<=cnt&&i*p[j]<=n;++j) {
            v[i*p[j]]=1;
            if (i%p[j]) f[i*p[j]]=1ll*f[i]*f[p[j]]%mod;
            else { f[i*p[j]]=1ll*f[i]*p[j]%mod*p[j]%mod*p[j]%mod; break; }
        }
    }
    for (int i=1;i<=n;++i) f[i]=(f[i-1]+f[i])%mod;
}

int S1(ll n) { n%=mod; return n*(n+1)%mod*inv2%mod;}
int S2(ll n) { n%=mod; return n*(n+1)%mod*(n+n+1)%mod*inv6%mod; }
int S3(ll n) { return 1ll*S1(n)*S1(n)%mod; }

unordered_map<ll,int> M;
int S(ll n) {
    if (n<=10000000) return f[n];
    if (M.count(n)) return M[n];
    int res=S3(n);
    for (ll l=2,r;l<=n;l=r+1) {
        r=n/(n/l);
        res=(res-1ll*S(n/l)*(S2(r)-S2(l-1)+mod)%mod+mod)%mod;
    }
    return M[n]=res;
}

int main() {
    mod=read(),inv2=qpow(2,mod-2),inv6=qpow(6,mod-2);
    n=read(); sieve(min(n,10000000ll));
    int ans=0;
    for (ll l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        ans=(ans+1ll*S1(n/l)*S1(n/l)%mod*(S(r)-S(l-1)+mod))%mod;
    }
    printf("%d\n",ans);
    return 0;
}