题目:link

---

原式
\[\sum_{i=1}^{n}\sum_{j=1}^n{ijgcd(i,j)}\]
枚举 \(gcd\) 得到
\[\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^3[gcd(i,j)=1]\]
把 \(d^3\) 提到前面
\[\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij[gcd(i,j)=1]\]
设函数
\[f(n)=\sum_{i=1}^{n}\sum_{j=1}^nij[gcd(i,j)=1]\]
\[s(n)=\sum_{i=1}^{n}\sum_{j=1}^nij\]
则原式等于
\[\sum_{d=1}^nd^3f(\lfloor\frac{n}{d}\rfloor)\]
可以看出是数论分块，其中立方部分和为
\[\sum_{i=1}^{n}i^3=(\sum_{i=1}^ni)^2=[\frac{n(n+1)}{2}]^2\]
求出前缀做差即可

---

考虑快速求 \(f(n)\)，非常显然
\[s(n)=\sum_{i=1}^{n}\sum_{j=1}^nij=[\frac{n(n+1)}{2}]^2\]
而对于 \(s(n)\) 同样枚举 \(gcd\) 可得
\[s(n)=\sum_{i=1}^{n}\sum_{j=1}^nij=\sum_{t=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{t}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{t}\rfloor}t^2ij[gcd(i,j)=1]\]
用 \(f(n)\) 替换上式得
\[s(n)=\sum_{t=1}^{n}t^2f(\lfloor\frac{n}{t}\rfloor)=f(n)+\sum_{t>1}t^2f(\lfloor\frac{n}{t}\rfloor)\]
所以有
\[f(n)=s(n)-\sum_{t>1}t^2f(\lfloor\frac{n}{t}\rfloor)=[\frac{n(n+1)}{2}]^2-\sum_{t>1}t^2f(\lfloor\frac{n}{t}\rfloor)\]
到了此刻已经非常容易求了，数论分块即可，其中
\[\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6}\]
做法同立方和，前缀做差

---

考虑预处理 \(f(n)\) 部分
\[f(n)=\sum_{i=1}^{n}\sum_{j=1}^nij[gcd(i,j)=1]=f(n-1)+2\sum_{i=1}^{n}in[gcd(i,n)=1]\]

那么需要快速求得
\[\sum_{i=1}^ni[gcd(i,n)=1]\]
发现这就是与 \(n\) 互质的数的和，就等于
\[n\frac{\varphi(n)}{2}\]
详细证明参考链接

带回原式得到
\[f(n)=f(n-1)+n^2\varphi(n)\]
线性筛预处理即可

---

Code：

#pragma GCC optimize(2)
#include<stdio.h>
#include<string.h>
#include<cassert>
#define Mod p
#define LL long long
#define N 5000000

template<class T>
inline void read(T &x){
    x=0;char c=getchar();T flag=1;
    while(c<'0'||c>'9'){if(c=='-')flag=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+c-48;c=getchar();}
    x*=flag;
}

LL n,p,inv,inv6;
int cnt,P[N];
bool mark[N];
LL phi[N],f[N],F[N];

inline LL qpow(LL a,LL b){
    if(!b) return 1LL;
    if(b==1LL) return a;
    LL t=qpow(a,b>>1LL);
    t=t*t%Mod;
    if(b&1LL) t=t*a%Mod;
    assert(t>=0);
    return t;
}
inline LL pow2(LL x){//平方前缀和
    return ((inv6*x)%Mod)*(((x+1LL)%Mod)*((2LL*x+1LL)%Mod)%Mod)%Mod;
}
inline LL pow2(LL l,LL r){平方区间和
    return ((pow2(r)-pow2(l-1LL))%Mod+Mod)%Mod;
}
inline LL pow3(LL x){立方前缀和
    return qpow(((x%Mod)*(x+1LL)%Mod)*inv%Mod,2LL);
}
inline LL pow3(LL l,LL r){立方区间和
    return ((pow3(r)-pow3(l-1LL))%Mod+Mod)%Mod;
}
inline void Prework(int x){//线性筛预处理 f(n)
    cnt=0;
    f[1]=phi[1]=1;
    for(int i=2;i<=x;++i){
        if(!mark[i]){
            P[cnt++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<cnt&&i*P[j]<=x;++j){
            mark[i*P[j]]=1;
            if(!(i%P[j])){
                phi[i*P[j]]=phi[i]*P[j];
                break;
            }else phi[i*P[j]]=phi[i]*(P[j]-1);
        }
        LL ret=qpow(i,2LL)*phi[i]%Mod;
        f[i]=(f[i-1]+ret)%Mod;
    }
}
inline LL solve(LL x){//杜教筛
    if(x<=N-3) return f[x];
    if(~F[n/x]) return F[n/x];
    F[n/x]=qpow(((x%Mod)*((x+1)%Mod)%Mod)*inv%Mod,2);
    for(LL l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        F[n/x]=((F[n/x]-(pow2(l,r)*solve(x/l)%Mod))%Mod+Mod)%Mod;
    }
    return F[n/x];
}
int main(){
    memset(F,-1,sizeof(F));
    read(p),read(n);
    inv=qpow(2LL,p-2);//2的逆元
    inv6=qpow(6LL,p-2);//6的逆元
    Prework(N-3);
    LL ans=0;
    for(LL l=1,r;l<=n;l=r+1){//原式
        r=n/(n/l);
        ans=(ans+(pow3(l,r)*solve(n/l))%Mod)%Mod;
    }
    printf("%lld",ans);
}
/*
998244353 2
*/