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

点此看题面

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

前置技能

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

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

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

推式子

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

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

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

[sum_{d=1}^nsum_{i=1}^{lfloorfrac nd floor}sum_{j=1}^{lfloorfrac nd floor}ijd^3[gcd(i,j)==1] ]

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

[sum_{d=1}^nd^3sum_{i=1}^{lfloorfrac nd floor}sum_{j=1}^{lfloorfrac nd floor}ij[gcd(i,j)==1] ]

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

[sum_{d=1}^nd^3sum_{p=1}^{lfloorfrac nd floor}mu(p)sum_{i=1}^{lfloorfrac n{d·p} floor}sum_{j=1}^{lfloorfrac n{d·p} floor}(i·p)(j·p) ]

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

[sum_{d=1}^nd^3sum_{p=1}^{lfloorfrac nd floor}mu(p)p^2(sum_{i=1}^{lfloorfrac n{d·p} floor}i)(sum_{j=1}^{lfloorfrac n{d·p} floor}j) ]

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

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

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

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

[sum_{s=1}^nsum_{d|s}d^3mu(frac sd)frac{s^2}{d^2}(frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+1)}2)^2 ]

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

[sum_{s=1}^ns^2(frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+1)}2)^2sum_{d|s}dmu(frac sd) ]

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

[(mu*id)(s) ]

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

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

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

[phi(s) ]

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

[sum_{s=1}^n(s^2phi(s))(frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+1)}2)^2 ]

于是，我们可以考虑杜教筛筛出(s^2phi(s))，然后除法分块求出前面一部分的值；而后面的((frac{{lfloorfrac n{d·p} floor}·({lfloorfrac n{d·p} floor}+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;
}