【51NOD 1237】最大公约数之和 V3

最新推荐文章于 2020-08-17 18:01:44 发布

原创最新推荐文章于 2020-08-17 18:01:44 发布 · 965 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#莫比乌斯反演 #欧拉函数 #杜教筛

莫比乌斯反演同时被 2 个专栏收录

23 篇文章

订阅专栏

杜教筛

8 篇文章

订阅专栏

本文探讨了数论中一个经典的求和问题：求∑i=1n∑j=1ngcd(i,j)的值，并给出了详细的解决方案。利用反演技巧简化了问题，并通过分块和杜教筛算法优化计算效率。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

求

\sum i = 1 n \sum j = 1 n gcd (i, j)

$\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)$

Solution

用反演，过程省略，

A n s = \sum d = 1 n d * \sum i = 1 ⌊ n d ⌋ μ (i) ⌊ n i d ⌋ 2

$Ans=\sum_{d=1}^n d*\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)\lfloor\frac{n}{id}\rfloor^2$

设 $T=id$

A n s = \sum T = 1 n ⌊ n T ⌋ 2 \sum d | T μ (d) * T d

$Ans=\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor^2\sum_{d|T}\mu(d)*\frac{T}{d}$
我们发现反演出了

φ $\varphi$ （证明点这里）

A n s = \sum T = 1 n ⌊ n T ⌋ 2 φ (T)

$Ans=\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor^2\varphi(T)$

用一个分块加上杜教筛即可。
欧拉函数的杜教筛点这里

复杂度： $O(n^{\frac{2}{3}})$

Code

#include <cstdio>
#include <cstdlib>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define min(q,w) ((q)<(w)?(q):(w))
using namespace std;
typedef long long LL;
const int N=5000500,M=2544657,mo=1e9+7;
const LL eni=500000004;
LL n,ans;
bool prz[N+10];
int pr[N+10];
int phi[N+10];
int Hx[M+1][2];
int HX(LL q)
{
    int i=q%M;
    while(Hx[i][0]&&Hx[i][0]!=q)i=(i+1)%M;
    return i;
}
LL Gphi(LL q)
{
    if(q<=N)return phi[q];
    int t=HX(q);
    if(Hx[t][0])return Hx[t][1];
    Hx[t][0]=q;
    LL ans=0;
    for(LL i=2,nx;i<=q;i=nx+1)
    {
        nx=q/(q/i);
        ans=(ans+(nx-i+1)%mo*Gphi(q/i))%mo;
    }
    q%=mo;
    return Hx[t][1]=(q*(q+1)%mo*eni%mo-ans+mo)%mo;
}
int main()
{
    phi[1]=1;
    fo(i,2,N)
    {
        if(!prz[i])pr[++pr[0]]=i,phi[i]=i-1;
        fo(j,1,pr[0])
        {
            int t=pr[j]*i;
            if(t>N)break;
            prz[t]=1;
            phi[t]=phi[i]*pr[j];
            if(i%pr[j]==0)break;
            phi[t]=phi[i]*(pr[j]-1);
        }
    }
    fo(i,2,N)phi[i]=(phi[i]+phi[i-1])%mo;
    scanf("%lld",&n);
    for(LL i=1,nx;i<=n;i=nx+1)
    {
        nx=n/(n/i);
        ans=(ans+(n/i)%mo*((n/i)%mo)%mo*(Gphi(nx)-Gphi(i-1))%mo)%mo;
    }
    printf("%lld\n",(ans+mo)%mo);
    return 0;
}