洛谷 - P3768 - 简单的数学题 - 欧拉函数 - 莫比乌斯反演
https://www.luogu.org/problemnew/show/P3768
\(F(n)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ijgcd(i,j)\)
首先加入方括号并枚举g，提gcd的g：
\(\sum\limits_{g=1}^{n}g\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij[gcd(i,j)==g]\)
后面的方括号里的g也可以提出来，注意前面有两个id，所以：
\(\sum\limits_{g=1}^{n}g^3 \sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{g}\rfloor}ij[gcd(i,j)==1]\)
记 \(G(n)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}ij[gcd(i,j)==1]\) ，则 \(F(n)=\sum\limits_{g=1}^{n}g^3 G(\lfloor\frac{n}{g}\rfloor)\) ，\(F(n)\)可以一次对\(G(\lfloor\frac{n}{g}\rfloor)\)分块求出来。
关注 \(G(n)\) 怎么求。
先把 \(G(n)\) 分成两半，设 \(H_1(n)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i}ij[gcd(i,j)==1]\) ，显然 \(G(n)=2*H_1(n)-1\)
又设\(H_2(n)=\sum\limits_{i=1}^{n}in[gcd(i,n)==1]\) ，则 $H_1(n)=\sum\limits_{i=1}^{n}H_2(i) $
而\(H_2(n)=n\sum\limits_{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)=\sum\limits_{i=1}^{n}H_2(i) = \sum\limits_{i=1}^{n} \frac{i^2}{2}([i==1]+\varphi(i))\)
故：
\(G(n)=\sum\limits_{i=1}^{n} i^2 ([i==1]+\varphi(i)) -1= \sum\limits_{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\n",F(n));
return 0;
}
