【洛谷3768】简单的数学题（莫比乌斯反演+杜教筛）

点此看题面

大致题意： 求\(\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)\%p\)。

前置技能

关于这道题目，我们首先需要了解以下知识：

• 莫比乌斯反演
• 杜教筛
• 狄利克雷卷积
• 除法分块

知道这些，你就可以对这题的做法有一定的理解了。

推式子

首先，按照常见的套路，我们可以枚举\(gcd(i,j)=d\)，得到这样一个式子：

\[\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^nijd[gcd(i,j)==d] \]

然后，我们可以进行一个比较简单的转化，将中括号里面的\(d\)提出：

\[\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ijd^3[gcd(i,j)==1] \]

不难发现\(d^3\)与\(i,j\)无关，于是我们可以将其提前：

\[\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij[gcd(i,j)==1] \]

然后我们可以枚举\(gcd(i,j)=p\)，并对原式做一遍莫比乌斯反演，可以得到这样一个式子：

\[\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac nd\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac n{d·p}\rfloor}\sum_{j=1}^{\lfloor\frac n{d·p}\rfloor}(i·p)(j·p) \]

然后将\(p\)提出，就变成了这样：

\[\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac nd\rfloor}\mu(p)p^2(\sum_{i=1}^{\lfloor\frac n{d·p}\rfloor}i)(\sum_{j=1}^{\lfloor\frac n{d·p}\rfloor}j) \]

对于式子中的\((\sum_{i=1}^{\lfloor\frac n{d·p}\rfloor}i)(\sum_{j=1}^{\lfloor\frac n{d·p}\rfloor}j)\)，可以化简得\((\frac{{\lfloor\frac n{d·p}\rfloor}·({\lfloor\frac n{d·p}\rfloor}+1)}2)^2\)。

那么原式就变成了这个样子：

\[\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac nd\rfloor}\mu(p)p^2(\frac{{\lfloor\frac n{d·p}\rfloor}·({\lfloor\frac n{d·p}\rfloor}+1)}2)^2 \]

如果设\(s=d·p\)，则原式就变成了这个样子：

\[\sum_{s=1}^n\sum_{d|s}d^3\mu(\frac sd)\frac{s^2}{d^2}(\frac{{\lfloor\frac ns\rfloor}·({\lfloor\frac ns\rfloor}+1)}2)^2 \]

注意到式子中的\(d^2\)似乎可以约去，于是我们可以得到下面这个式子：

\[\sum_{s=1}^ns^2(\frac{{\lfloor\frac ns\rfloor}·({\lfloor\frac ns\rfloor}+1)}2)^2\sum_{d|s}d\mu(\frac sd) \]

对于后面的一个式子\(\sum_{d|s}d\mu(\frac sd)\)，用狄利克雷卷积的形式来表示就是这个样子：

\[(\mu*id)(s) \]

而众所周知，\(id=I*\phi\)，所以我们可以考虑将其代入，得到：

\[(\mu*I*\phi)(s) \]

而\(\mu*I=e\)刚好约掉，所以这个式子的结果就是：

\[\phi(s) \]

再代回原式就化简出来一个比较简单的式子了：

\[\sum_{s=1}^n(s^2\phi(s))(\frac{{\lfloor\frac ns\rfloor}·({\lfloor\frac ns\rfloor}+1)}2)^2 \]

于是，我们可以考虑杜教筛筛出\(s^2\phi(s)\)，然后除法分块求出前面一部分的值；而后面的\((\frac{{\lfloor\frac ns\rfloor}·({\lfloor\frac ns\rfloor}+1)}2)^2\)一项是可以\(O(1)\)计算的。

于是这道题就这样解决了。

代码

#include<bits/stdc++.h>
#define LL long long
#define sum_sqr(x) ((((x)%MOD)*(((x)+1)%MOD))%MOD*((((x)<<1)+1)%MOD)%MOD*InvSix%MOD)//求1^2+2^2+...+x^2的结果
#define sum_cube(x) ((((x)%MOD)*(((x)+1)%MOD)%MOD*InvTwo%MOD)*(((x)%MOD)*(((x)+1)%MOD)%MOD*InvTwo%MOD)%MOD)//求1^3+2^3+...+x^3的结果
#define Inc(x,y) ((x+=(y))>=MOD&&(x-=MOD))
#define Dec(x,y) ((x-=(y))<0&&(x+=MOD))
using namespace std;
LL n,MOD,InvTwo,InvSix;//分别用InvTwo和InvSix存储模MOD意义下2和6的逆元
inline LL Sum(LL x,LL y) {return Inc(x,y),x;}
inline LL Sub(LL x,LL y) {return Dec(x,y),x;}
class Class_DuSieve//杜教筛
{
    private:
        #define Size 5000000
        int Prime_cnt,Prime[Size+5],phi[Size+5];bool IsNotPrime[Size+5];LL sum_phi[Size+5];
        map<LL,LL> MemoryPhi;
    public:
        inline void Init()
        {
            register LL i,j;
            for(phi[1]=1,i=2;i<=Size;++i)
            {
                if(!IsNotPrime[i]) phi[Prime[++Prime_cnt]=i]=i-1;
                for(j=1;j<=Prime_cnt&&i*Prime[j]<=Size;++j)
                    if(IsNotPrime[i*Prime[j]]=true,i%Prime[j]) phi[i*Prime[j]]=phi[i]*(Prime[j]-1);else {phi[i*Prime[j]]=phi[i]*Prime[j];break;}
            }
            for(i=1;i<=Size;++i) sum_phi[i]=Sum(sum_phi[i-1],i*i%MOD*phi[i]%MOD);
        }
        inline LL GetPhi(LL x)
        {
            if(x<=Size) return sum_phi[x];
            if(MemoryPhi[x]) return MemoryPhi[x];
            register LL l,r,res=sum_cube(x);
            for(l=2;l<=x;l=r+1) r=x/(x/l),Dec(res,Sub(sum_sqr(r),sum_sqr(l-1))*GetPhi(x/l)%MOD);//递归调用
            return MemoryPhi[x]=res;
        }
}DuSieve;
inline LL quick_pow(LL x,LL y,register LL res=1)//快速幂求逆元
{
    for(;y;(x*=x)%=MOD,y>>=1) if(y&1) (res*=x)%=MOD;
    return res;
}
int main()
{
    register LL l,r,ans=0;
    scanf("%lld%lld",&MOD,&n),DuSieve.Init(),InvTwo=quick_pow(2,MOD-2),InvSix=quick_pow(6,MOD-2);//预处理
    for(l=1;l<=n;l=r+1) r=n/(n/l),Inc(ans,Sub(DuSieve.GetPhi(r),DuSieve.GetPhi(l-1))*sum_cube(n/l)%MOD);//除法分块求答案
    return printf("%lld",ans),0;
}