洛谷

https://www.luogu.org/problemnew/show/P3768

(F(n)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}ijgcd(i,j))

首先加入方括号并枚举g，提gcd的g：
(sumlimits_{g=1}^{n}gsumlimits_{i=1}^{n}sumlimits_{j=1}^{n}ij[gcd(i,j)==g])

后面的方括号里的g也可以提出来，注意前面有两个id，所以：
(sumlimits_{g=1}^{n}g^3 sumlimits_{i=1}^{lfloorfrac{n}{g} floor}sumlimits_{j=1}^{lfloorfrac{n}{g} floor}ij[gcd(i,j)==1])

记 (G(n)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}ij[gcd(i,j)==1]) ，则 (F(n)=sumlimits_{g=1}^{n}g^3 G(lfloorfrac{n}{g} floor)) ，(F(n))可以一次对(G(lfloorfrac{n}{g} floor))分块求出来。

关注 (G(n)) 怎么求。

先把 (G(n)) 分成两半，设 (H_1(n)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{i}ij[gcd(i,j)==1]) ，显然 (G(n)=2*H_1(n)-1)

又设(H_2(n)=sumlimits_{i=1}^{n}in[gcd(i,n)==1]) ，则 $H_1(n)=sumlimits_{i=1}^{n}H_2(i) $

而(H_2(n)=nsumlimits_{i=1}^{n}i[gcd(i,n)==1])，这个不是在疯狂lcm里面见过吗？和一个数gcd为1的不比他大的数的和。

直接套上去：

(H_2(n)=nH(n)=frac{n^2}{2}([n==1]+varphi(n)))

所以：
(H_1(n)=sumlimits_{i=1}^{n}H_2(i) = sumlimits_{i=1}^{n} frac{i^2}{2}([i==1]+varphi(i)))

故：
(G(n)=sumlimits_{i=1}^{n} i^2 ([i==1]+varphi(i)) -1= sumlimits_{i=1}^{n} i^2 varphi(i))

感觉可以用杜教筛来求，先线性筛出1e7以内的 (G(n))。

---

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

ll n;
int mod;
int inv2;
const int MAXN=1e7;

int pri[MAXN+1];
int &pritop=pri[0];
int B[2*MAXN+1];
int pk[MAXN+1];

void sieve(int n=MAXN) {
    memset(B,-1,sizeof(B));
    B[0]=0;
    pk[1]=1;
    B[1]=1;
    for(int i=2; i<=n; i++) {
        if(!pri[i]) {
            pri[++pritop]=i;
            pk[i]=i;
            ll tmp=1ll*i*i;
            if(tmp>=mod)
                tmp%=mod;
            tmp*=(i-1);
            if(tmp>=mod)
                tmp%=mod;
            B[i]=tmp;
        }
        for(int j=1; j<=pritop; j++) {
            int &p=pri[j];
            int t=i*p;
            if(t>n)
                break;
            pri[t]=1;
            if(i%p) {
                pk[t]=p;
                ll tmp=1ll*B[i]*B[p];
                if(tmp>=mod)
                    tmp%=mod;
                B[t]=tmp;
            } else {
                pk[t]=pk[i]*p;
                if(pk[t]==t) {
                    ll tmp=1ll*t*t;
                    if(tmp>=mod)
                        tmp%=mod;
                    tmp*=(t-i);
                    if(tmp>=mod)
                        tmp%=mod;
                    B[t]=tmp;
                } else {
                    ll tmp=1ll*B[t/pk[t]]*B[pk[t]];
                    if(tmp>=mod)
                        tmp%=mod;
                    B[t]=tmp;
                }
                break;
            }
        }
    }
    for(int i=1; i<=n; i++) {
        ll tmp=(ll)B[i]+B[i-1];
        if(tmp>=mod)
            tmp-=mod;
        B[i]=tmp;
    }
}

inline ll qpow(ll x,ll n) {
    if(x>=mod)
        x%=mod;
    ll res=1;
    while(n) {
        if(n&1) {
            res*=x;
            if(res>=mod)
                res%=mod;
        }
        x*=x;
        if(x>=mod)
            x%=mod;
        n>>=1;
    }
    return res;
}

int inv4;
int inv6;

inline int s2(ll n) {
    if(n>=mod)
        n%=mod;
    ll tmp=n*(n+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=(n*2+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=inv6;
    if(tmp>=mod)
        tmp%=mod;
    return tmp;
}

inline int s3(ll n) {
    if(n>=mod)
        n%=mod;
    ll tmp=n*(n+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=tmp;
    if(tmp>=mod)
        tmp%=mod;
    tmp*=inv4;
    if(tmp>=mod)
        tmp%=mod;
    return tmp;
}

inline int id(ll x){
    if(x<=MAXN)
        return x;
    else
        return n/x+MAXN;
}

inline int S(ll n) {
    int idn=id(n);
    if(B[idn]!=-1)
        return B[idn];
    ll ret=s3(n);
    for(ll l=2,r; l<=n; l=r+1) {
        ll t=n/l;
        r=n/t;
        ll tmp=s2(r)-s2(l-1);
        if(tmp<0)
            tmp+=mod;
        tmp*=S(t);
        if(tmp>=mod)
            tmp%=mod;
        ret-=tmp;
    }
    ret%=mod;
    if(ret<0)
        ret+=mod;
    return B[idn]=ret;
}

inline int F(ll n){
    ll res=0;
    for(ll l=1,r;l<=n;l=r+1){
        ll t=n/l;
        r=n/t;
        ll tmp=s3(r)-s3(l-1);
        if(tmp<0)
            tmp+=mod;
        tmp*=S(t);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
    }
    if(res>=mod)
        res%=mod;
    return res;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    scanf("%d%lld",&mod,&n);
    inv2=qpow(2,mod-2);
    inv4=qpow(4,mod-2);
    inv6=qpow(6,mod-2);
    sieve(min(n,(ll)MAXN));
    printf("%d
",F(n));
    return 0;
}