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

题意简述

求出这个式子

[sum_{i=1}^nsum_{j=1}^n ij(i,j) mod p ]

做法

先用莫比乌斯反演拆一下式子

[egin{split} sum_{i=1}^nsum_{j=1}^n ij(i,j)&=sum_{d=1}^ndsum_{i=1}^nsum_{j=1}^n ij[(i,j)=d]\ &=sum_{d=1}^ndsum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{n}{d} floor} id imes jd[(i,j)=1]\ &=sum_{d=1}^nd^3sum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{n}{d} floor} ij[(i,j)=1]\ &=sum_{d=1}^nd^3sum_{i=1}^{lfloor frac{n}{d} floor}sum_{j=1}^{lfloor frac{n}{d} floor} ij sum_{kmid i,kmid j} mu(k)\ &=sum_{d=1}^nd^3sum_{k=1}^{lfloor frac{n}{d} floor} mu(k) sum_{i=1}^{lfloor frac{n}{kd} floor}sum_{j=1}^{lfloor frac{n}{kd} floor} ik imes jk \ &=sum_{d=1}^nd^3sum_{k=1}^{lfloor frac{n}{d} floor} mu(k) k^2sum_{i=1}^{lfloor frac{n}{kd} floor}sum_{j=1}^{lfloor frac{n}{kd} floor} ij \ &=sum_{T=1}^n sum_{dmid T}d^3 muleft({T over d} ight) left({T over d} ight)^2sum_{i=1}^{lfloor frac{n}{T} floor}sum_{j=1}^{lfloor frac{n}{T} floor} ij \ &=sum_{T=1}^n T^2 sum_{dmid T}d muleft({T over d} ight) sum_{i=1}^{lfloor frac{n}{T} floor} i^3 \ &=sum_{T=1}^n T^2 varphi (T) sum_{i=1}^{lfloor frac{n}{T} floor} i^3 \ end{split} ]

所以杜教筛直接做(sum_i varphi(i))，此题结束。

代码实现

请使用C++11

#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;
	return p1==p2&&(p2=(p1=buf)+fread(buf,1,10000,stdin),p1==p2)?EOF:*p1++;
} 
char qwq;
#define gc() getch()
in ll read()
{
	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()
{
    mod=read();get();
    ll n=read(),ans=0ll;
    inv2=qpow(2,mod-2);inv6=qpow(6,mod-2);
    for(ll l=1,r,s;l<=n;l=r+1)
    {
        r=n/(n/(l));
        ans=(ans+(get_phi(r)-get_phi(l-1)+mod)%mod*get1(1,n/l)%mod)%mod;
        ans=(ans+mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}