题意:求$\sum\limits_{i=1}^n\sum\limits_{j=1}^i\sum\limits_{k=1}^i(i,j,k)$
先令$f(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^n(n,i,j)$
$\begin{align*}f(n)&=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sum\limits_{d=1}^n\left[(n,i,j)=d\right]d\\&=\sum\limits_{d|n}\sum\limits_{i=1}^{\frac nd}\sum\limits_{j=1}^{\frac nd}\left[\left(\dfrac nd,i,j\right)=1\right]d\\&=\sum\limits_{d|n}\sum\limits_{i=1}^{\frac nd}\sum\limits_{j=1}^{\frac nd}\sum\limits_{\substack{e|\frac nd\\e|i\\e|j}}d\mu(e)\\&=\sum\limits_{d|n}\sum\limits_{de|n}d\mu(e)\left(\dfrac n{de}\right)^2\end{align*}$
令$T=de$,枚举$T$并将原来的$e$代换
$\begin{align*}f(n)&=\sum\limits_{d|n}\sum\limits_{T|n}d\mu\left(\dfrac Td\right)\left(\dfrac nT\right)^2\\&=\sum\limits_{T|n}\sum\limits_{d|n}d\mu\left(\dfrac Td\right)\left(\dfrac nT\right)^2\end{align*}$
因为式中有$\mu\left(\dfrac Td\right)$所以必须满足$d|T$,也就是改写第二个sigma
$\begin{align*}f(n)&=\sum\limits_{T|n}\sum\limits_{d|T}d\mu\left(\dfrac Td\right)\left(\dfrac nT\right)^2\\&=\sum\limits_{T|n}\varphi(T)\left(\dfrac nT\right)^2\\&=\sum\limits_{T|n}\varphi\left(\dfrac nT\right)T^2\end{align*}$
$\begin{align*}\sum\limits_{i=1}^nf(i)&=\sum\limits_{i=1}^n\sum\limits_{T|i}\varphi\left(\dfrac iT\right)T^2\\&=\sum\limits_{T=1}^nT^2\sum\limits_{j=1}^{\left\lfloor\frac nT\right\rfloor}\varphi(j)\end{align*}$
除法下取整$O\left(\sqrt n\right)$枚举除数,欧拉函数求和直接上杜教筛即可
#include<stdio.h>
#include<map>
using namespace std;
typedef long long ll;
const int T=10000000;
int mod,inv6,inv2;
int mul(int a,int b){return a*(ll)b%mod;}
int ad(int a,int b){return(a+b)%mod;}
int de(int a,int b){return(a-b)%mod;}
int pr[10000010],phi[10000010];
bool np[10000010];
void sieve(){
int i,j,m;
np[1]=1;
phi[1]=1;
m=0;
for(i=2;i<=T;i++){
if(!np[i]){
m++;
pr[m]=i;
phi[i]=i-1;
}
for(j=1;j<=m;j++){
if(pr[j]*(ll)i>T)break;
np[i*pr[j]]=1;
if(i%pr[j]==0){
phi[i*pr[j]]=phi[i]*pr[j];
break;
}
phi[i*pr[j]]=phi[i]*(pr[j]-1);
}
}
for(i=2;i<=T;i++)phi[i]=ad(phi[i],phi[i-1]);
}
int pow(int a,int b){
int s=1;
while(b){
if(b&1)s=mul(s,a);
a=mul(a,a);
b>>=1;
}
return s;
}
int s(int n){return mul(mul(n,n+1),mul(n<<1|1,inv6));}
map<int,int>res;
map<int,int>::iterator it;
int dj(int n){
if(n<=T)return phi[n];
it=res.find(n);
if(it!=res.end())return it->second;
int s,i,nex;
s=mul(mul(n,n+1),inv2);
for(i=2;i<=n;i=nex+1){
nex=n/(n/i);
s=de(s,mul(nex-i+1,dj(n/i)));
}
return res[n]=s;
}
int main(){
int n,i,nex,res;
scanf("%d%d",&n,&mod);
sieve();
inv6=pow(6,mod-2);
inv2=pow(2,mod-2);
res=0;
for(i=1;i<=n;i=nex+1){
nex=n/(n/i);
res=ad(res,mul(s(nex)-s(i-1),dj(n/i)));
}
printf("%d",(res+mod)%mod);
}