大致题意： 求(sum_{i=1}^nsum_{j=1}^nijgcd(i,j)\%p)。
前置技能
关于这道题目，我们首先需要了解以下知识：
知道这些，你就可以对这题的做法有一定的理解了。
推式子
首先，按照常见的套路，我们可以枚举(gcd(i,j)=d)，得到这样一个式子：
然后，我们可以进行一个比较简单的转化，将中括号里面的(d)提出：
不难发现(d^3)与(i,j)无关，于是我们可以将其提前：
然后我们可以枚举(gcd(i,j)=p)，并对原式做一遍莫比乌斯反演，可以得到这样一个式子：
然后将(p)提出，就变成了这样：
对于式子中的((sum_{i=1}^{lfloorfrac n{d·p} floor}i)(sum_{j=1}^{lfloorfrac n{d·p} floor}j))，可以化简得((frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+1)}2)^2)。
那么原式就变成了这个样子：
如果设(s=d·p)，则原式就变成了这个样子：
注意到式子中的(d^2)似乎可以约去，于是我们可以得到下面这个式子：
对于后面的一个式子(sum_{d|s}dmu(frac sd))，用狄利克雷卷积的形式来表示就是这个样子：
而众所周知，(id=I*phi)，所以我们可以考虑将其代入，得到：
而(mu*I=e)刚好约掉，所以这个式子的结果就是：
再代回原式就化简出来一个比较简单的式子了：
于是，我们可以考虑杜教筛筛出(s^2phi(s))，然后除法分块求出前面一部分的值；而后面的((frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+1)}2)^2)一项是可以(O(1))计算的。
于是这道题就这样解决了。
代码
#include<bits/stdc++.h>
#define LL long long
#define sum_sqr(x) ((((x)%MOD)*(((x)+1)%MOD))%MOD*((((x)<<1)+1)%MOD)%MOD*InvSix%MOD)//求1^2+2^2+...+x^2的结果
#define sum_cube(x) ((((x)%MOD)*(((x)+1)%MOD)%MOD*InvTwo%MOD)*(((x)%MOD)*(((x)+1)%MOD)%MOD*InvTwo%MOD)%MOD)//求1^3+2^3+...+x^3的结果
#define Inc(x,y) ((x+=(y))>=MOD&&(x-=MOD))
#define Dec(x,y) ((x-=(y))<0&&(x+=MOD))
using namespace std;
LL n,MOD,InvTwo,InvSix;//分别用InvTwo和InvSix存储模MOD意义下2和6的逆元
inline LL Sum(LL x,LL y) {return Inc(x,y),x;}
inline LL Sub(LL x,LL y) {return Dec(x,y),x;}
class Class_DuSieve//杜教筛
{
private:
#define Size 5000000
int Prime_cnt,Prime[Size+5],phi[Size+5];bool IsNotPrime[Size+5];LL sum_phi[Size+5];
map<LL,LL> MemoryPhi;
public:
inline void Init()
{
register LL i,j;
for(phi[1]=1,i=2;i<=Size;++i)
{
if(!IsNotPrime[i]) phi[Prime[++Prime_cnt]=i]=i-1;
for(j=1;j<=Prime_cnt&&i*Prime[j]<=Size;++j)
if(IsNotPrime[i*Prime[j]]=true,i%Prime[j]) phi[i*Prime[j]]=phi[i]*(Prime[j]-1);else {phi[i*Prime[j]]=phi[i]*Prime[j];break;}
}
for(i=1;i<=Size;++i) sum_phi[i]=Sum(sum_phi[i-1],i*i%MOD*phi[i]%MOD);
}
inline LL GetPhi(LL x)
{
if(x<=Size) return sum_phi[x];
if(MemoryPhi[x]) return MemoryPhi[x];
register LL l,r,res=sum_cube(x);
for(l=2;l<=x;l=r+1) r=x/(x/l),Dec(res,Sub(sum_sqr(r),sum_sqr(l-1))*GetPhi(x/l)%MOD);//递归调用
return MemoryPhi[x]=res;
}
}DuSieve;
inline LL quick_pow(LL x,LL y,register LL res=1)//快速幂求逆元
{
for(;y;(x*=x)%=MOD,y>>=1) if(y&1) (res*=x)%=MOD;
return res;
}
int main()
{
register LL l,r,ans=0;
scanf("%lld%lld",&MOD,&n),DuSieve.Init(),InvTwo=quick_pow(2,MOD-2),InvSix=quick_pow(6,MOD-2);//预处理
for(l=1;l<=n;l=r+1) r=n/(n/l),Inc(ans,Sub(DuSieve.GetPhi(r),DuSieve.GetPhi(l-1))*sum_cube(n/l)%MOD);//除法分块求答案
return printf("%lld",ans),0;
}