Luogu 1390

Luogu 1390 - 公约数的和

UVa 11426 - GCD - Extreme (II)

题意

给定数字(n)，计算下式结果

[sum_{i=1}^nsum_{j=i+1}^ngcd(i,j) ]

限制

(Case=1, nleq2 imes10^6, 1000ms) （洛谷）

(Caseleq100, nleq4 imes10^6, 10000ms) （UVa）

思路一（欧拉函数）

刘汝佳蓝书 P124-P125

设(f(n)=gcd(1,n)+gcd(2,n)+...+gcd(n-1,n)=sum_{i=1}^{n-1}gcd(i,n))

则所求答案即为(sum_{i=1}^nsum_{j=i+1}^ngcd(i,j)=f(2)+f(3)+...+f(n)=sum_{i=2}^nf(i))

洛谷仅单测试点，故最后直接(O(n))累加所有(f(i))输出即可

UVa有多测试点，需要预处理前缀和后再输出

---

设(g(n,i))为满足(gcd(n,j)=i, j<n)的所有满足条件的(j)的个数

即(g(n,i)=sum_{j=1}^{n-1}[gcd(n,j)=i])

又因为(gcd(n,j)=iiff gcd(frac n i,frac j i)=1)

即(g(n,i))表示所有小于(n/i)的并与(n/i)互质的数的个数

根据欧拉函数可得，所有满足(gcd(frac n i,frac j i)=1)的(frac j i)个数即为(varphi(frac n i))

所以(g(n,i)=varphi(frac n i))，可以根据欧拉筛法预处理求出欧拉函数先

然后再考虑(f(n)=sum i imes g(n,i), i|n)

要想得到(f(n))，首先就得找到所有小于(n)的因子（不包括(n)）

所以可以根据埃氏筛法，从小到大枚举因子，将答案加给因子的倍数即可

代码一 - Luogu 1390

Case Max (291ms/1000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;

ll eular[N+50],prim[N+50];
ll f[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            eular[i]=i-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                eular[k]=eular[i]*prim[j];
                break;
            }
            eular[k]=eular[i]*eular[prim[j]];
        }
    }
}

int main()
{
    int n;
    scanf("%d",&n);
    init(n);
    for(int i=1;i<=n;i++)
    {
        for(int j=i*2;j<=n;j+=i)
            f[j]+=eular[j/i]*i;
    }
    ll ans=0;
    for(int i=2;i<=n;i++)
        ans+=f[i];
    printf("%lld
",ans);
    
    return 0;
}

代码一 - UVa 11426

(630ms/10000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;

ll eular[N+50],prim[N+50];
ll f[N+50],ans[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            eular[i]=i-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                eular[k]=eular[i]*prim[j];
                break;
            }
            eular[k]=eular[i]*eular[prim[j]];
        }
    }
}

int main()
{
    init(N);
    for(int i=1;i<=N;i++)
    {
        for(int j=i*2;j<=N;j+=i)
            f[j]+=eular[j/i]*i;
    }
    for(int i=2;i<=N;i++)
        ans[i]=ans[i-1]+f[i];
    
    int n;
    while(scanf("%d",&n)!=EOF&&n)
    {
        printf("%lld
",ans[n]);
    }
    
    return 0;
}

思路二（莫比乌斯反演）（简单）

类似于思路一，我们将所有不同的(gcd(i,j)=k)分开来考虑，以“数量*乘积”的方式表示答案

则原式则会变成

[sum_{i=1}^nsum_{j=i+1}^ngcd(i,j)=sum_{k=1}^nksum_{i=1}^nsum_{j=i+1}^n[gcd(i,j)=k] ]

这里我们需要枚举(k)来求和所有(gcd(i,j)=k)的情况

因为(gcd(a,b)=ciff gcd(frac{a}{c},frac{b}{c})=1)

我们还能继续化简上式，得到

[sum_{k=1}^nksum_{i=1}^nsum_{j=i+1}^n[gcd(i,j)=k]=sum_{k=1}^nksum_{i=1}^{lfloorfrac n k floor}sum_{j=i+1}^{lfloorfrac n k floor}[gcd(i,j)=1]=sum_{k=1}^nksum_{i=1}^{lfloorfrac n k floor}mu(x)lfloor frac {lfloorfrac n k floor} x floor lfloor frac {lfloorfrac n k floor} x floor ]

明显发现，该式子后半部分即为莫比乌斯反演的模板题——求两区间内互质数的对数

直接套板子即可

代码二 - Luogu 1390

Case Max(139ms/1000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;

ll mu[N+50],prim[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
}

int main()
{
    int n;
    scanf("%d",&n);
    init(n);
    ll ans=0;
    for(int d=1;d<=n;d++)
    {
        ll ansd=0;
        int m=n/d;
        for(int i=1;i<=m;i++)
            ansd+=mu[i]*(m/i)*(m/i);
        ans+=(ansd>>1)*d;
    }
    printf("%lld
",ans);
    
    return 0;
}

代码二 - UVa 11426

(2310ms/10000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;

ll mu[N+50],prim[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
}

int main()
{
    init(N);
    int n;
    while(scanf("%d",&n)&&n)
    {
        ll ans=0;
        for(int d=1;d<=n;d++)
        {
            ll ansd=0;
            int m=n/d;
            for(int i=1;i<=m;i++)
                ansd+=mu[i]*(m/i)*(m/i);
            ans+=(ansd>>1)*d;
        }
        printf("%lld
",ans);
    }
    
    return 0;
}

思路三（莫比乌斯反演、分块）（进阶）

在思路二中，我们将式子化简到了直接套板子的情况

但是会发现，我们每次枚举(d)时，对于(gcd=d)的种类数就必须计算(frac n d)次才能出结果

总时间复杂度严格为(O(nsqrt n))

有没有办法再继续降低复杂度呢？当然，我们可以引入分块的思想

可以发现，虽然枚举的(d)在数值小的时候(frac n d)变化很大，以至于我们每次都需要重新计算一次

但是如果(d)逐渐变大，由于整除的关系，会有一段连续的(d)使得(frac n d)相同

那我们就可以根据这个性质，使得这一整段只计算一次就得出结果

假设当前(d)枚举到某个分块的左边界，那么(d_r=lfloor frac n {lfloor frac n d floor} floor)即为这一分块的右边界，即满足(lfloorfrac n d floor=lfloorfrac n {d_r} floor)的最大的(d_r)

代码三 - Luogu 1390

Case Max(79ms/1000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2000000;

ll mu[N+50],prim[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
}

int main()
{
    int n;
    scanf("%d",&n);
    init(N);
    ll ans=0;
    for(int d=1;d<=n;)
    {
        int m=n/d;
        ll ansd=0;
        for(int i=1;i<=m;i++)
            ansd+=mu[i]*(m/i)*(m/i);
        ansd>>=1;
        int nxt=n/(n/d); //计算这一块的右边界
        ans+=ansd*d;
        for(int i=d+1;i<=nxt;i++)
            ans+=ansd*i; //需要遍历乘上权值
        d=nxt+1; //d跳到下一分块的左边界
    }
    printf("%lld
",ans);
    
    return 0;
}

代码三 - UVa 11426

(1650ms/10000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4000000;

ll mu[N+50],prim[N+50];
bool vis[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    int p=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[p++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<p;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
}

int main()
{
    init(N);
    int n;
    while(scanf("%d",&n)&&n)
    {
        ll ans=0;
        for(int d=1;d<=n;)
        {
            int m=n/d;
            ll ansd=0;
            for(int i=1;i<=m;i++)
                ansd+=mu[i]*(m/i)*(m/i);
            ansd>>=1;
            int nxt=n/(n/d); //计算这一块的右边界
            ans+=ansd*d;
            for(int i=d+1;i<=nxt;i++)
                ans+=ansd*i; //需要遍历乘上权值
            d=nxt+1; //d跳到下一分块的左边界
        }
        printf("%lld
",ans);
    }
    
    return 0;
}