# 洛谷P3768 简单的数学题 莫比乌斯反演+杜教筛

## 题意简述

$\sum_{i=1}^n\sum_{j=1}^n ij(i,j) \bmod p$

## 做法

$\begin{split} \sum_{i=1}^n\sum_{j=1}^n ij(i,j)&=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n ij[(i,j)=d]\\ &=\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} id\times jd[(i,j)=1]\\ &=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} ij[(i,j)=1]\\ &=\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor} ij \sum_{k\mid i,k\mid j} \mu(k)\\ &=\sum_{d=1}^nd^3\sum_{k=1}^{\lfloor \frac{n}{d}\rfloor} \mu(k) \sum_{i=1}^{\lfloor \frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{kd}\rfloor} ik\times jk \\ &=\sum_{d=1}^nd^3\sum_{k=1}^{\lfloor \frac{n}{d}\rfloor} \mu(k) k^2\sum_{i=1}^{\lfloor \frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{kd}\rfloor} ij \\ &=\sum_{T=1}^n \sum_{d\mid T}d^3 \mu\left({T \over d}\right) \left({T \over d}\right)^2\sum_{i=1}^{\lfloor \frac{n}{T}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{T}\rfloor} ij \\ &=\sum_{T=1}^n T^2 \sum_{d\mid T}d \mu\left({T \over d}\right) \sum_{i=1}^{\lfloor \frac{n}{T}\rfloor} i^3 \\ &=\sum_{T=1}^n T^2 \varphi (T) \sum_{i=1}^{\lfloor \frac{n}{T}\rfloor} i^3 \\ \end{split}$

## 代码实现

#include<bits/stdc++.h>
using namespace std;
#define re register int
#define db double
#define ll long long
#define ak *
#define in inline
in char getch()
{
static char buf[10000],*p1=buf,*p2=buf;
}
char qwq;
#define gc() getch()
{
ll cz=0,ioi=1;qwq=gc();
while(qwq<'0'||qwq>'9') ioi=qwq=='-'?~ioi+1:1,qwq=gc();
while(qwq>='0'&&qwq<='9') cz=(cz<<3)+(cz<<1)+(qwq^48),qwq=gc();
return cz ak ioi;
}
const int lim=5000000;
ll inv2,inv6,mod,a,b,c,d,k,t,mu[5000005],vis[5000005],prime[500005],phi[5000005];
unordered_map<ll,ll>smu,sphi;
in void get()
{
mu[1]=phi[1]=1;
for(re i=2;i<=lim;i++)
{
if(!vis[i]) prime[++prime[0]]=i,mu[i]=-1,phi[i]=i-1;
for(re j=1;j<=prime[0];j++)
{
if(i*prime[j]>lim) break;
vis[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=0;phi[i*prime[j]]=prime[j]*phi[i]%mod;break;}
else mu[i*prime[j]]=-mu[i],phi[i*prime[j]]=(prime[j]-1)*phi[i]%mod;
}
}
for(re i=2;i<=lim;i++) phi[i]=(phi[i]*1ll*i%mod*i)%mod;
for(re i=2;i<=lim;i++) phi[i]=(phi[i-1]+phi[i])%mod;
}
in ll get1(ll l,ll r)
{
l%=mod;r%=mod;
ll dy=(l+r)%mod*(r-l+1)%mod*inv2%mod;
return dy*dy%mod;
}
in ll get2(ll l,ll r)
{
l%=mod;r%=mod;
return (r%mod*(r+1)%mod*(2ll*r+1)%mod*inv6%mod-(l-1)%mod*(2ll*(l-1)+1)%mod*l%mod*inv6%mod+mod)%mod;
}
ll get_phi(ll n)
{
if(n<=lim) return phi[n];
if(sphi[n]) return sphi[n];
ll res=get1(1,n);
for(ll l=2,r;l<=n;l=r+1)
r=n/(n/l),res=(res-get2(l,r)*get_phi(n/l)%mod+mod)%mod;
return sphi[n]=res;
}
inline ll qpow(ll x,ll y,ll z=1)
{
for(;y;y>>=1,x=x*x%mod) z=(y&1)?x*z%mod:z;
return z;
}
int main()
{