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;
}