UVA11424 GCD - Extreme (I) 题解

题目链接：UVA11424 GCD - Extreme (I)

题意：

给定 $n$ ，求
$\sum\limits_{i=1}^n\sum\limits_{j=i+1}^n\gcd(i,j)$
$1 < n \le 2\times 10^5$ ，最多 $2 \times 10^4$ 组数据。

据说是「7倍经验」，那我干脆一并写了吧。

这道题有很多种解法，常见的解法是欧拉函数和莫比乌斯反演。

解法1：欧拉函数

正解解法。记 $f(n) = \sum_{i = 1} ^ {n} \gcd(i,n)$ ，则答案为 $\sum_{i = 1}^{n}f(i)$

考虑怎么快速计算 $f$ 。设 $g(i,n)$ 为对于所有的 $x(1\le x < n)$ ， $\gcd(x,n) = i$ 的 $x$ 的个数，则

$f(n) = \sum_{d \mid n}g(d,n)$

注意到 $\gcd(x,n)=i$ ，意味着 $\gcd(x/i, n/i) = 1$ ，也就是 $x/i$ 与 $n/i$ 互质，故

$f(n) = \sum_{d \mid n}\varphi(n / i)$

~~感觉有点对劲，但好像又不是很对劲~~。直接枚举因数会比较慢，因此我们可以反过来计算贡献。

时间复杂度 $\mathcal{O}(n \log n)$

代码：

cpp

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(2e5 + 15))


bitset<N> ck; int pcnt,phi[N],f[N],sum[N],prime[N];
void Euler(int n = N - 10)
{
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!ck[i]) { prime[++pcnt] = i; phi[i] = i - 1; }
        for(int j = 1; j <= pcnt && i * prime[j] <= n; j++)
        {
            int pos = i * prime[j]; ck[pos] = true;
            if(i % prime[j]) { 
                phi[pos] = phi[i] * phi[prime[j]];
            } else { phi[pos] = phi[i] * prime[j]; break; }
        }
    }
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    Euler();
    for(int i = 1; i <= N - 10; i++)
        for(int j = i * 2; j <= N - 10; j += i) f[j] += i * phi[j / i];
    sum[2] = f[2];
    for(int i = 3; i <= N - 10; i++) sum[i] = sum[i - 1] + f[i];
    for(int n; cin >> n && n != 0; ) cout << sum[n] << '\n';
    return 0;
}

解法2：莫比乌斯反演

这个也是我自己想到的解法（不是正解），不过在这道题就需要卡一卡常。

我们先来推推式子：（这个式子的结果 $x$ 与原答案 $x_0$ 的关系为 $x=2x_0+\frac{n(n+1)}{2}$ ）

$\begin{aligned} \sum_{i=1}^n \sum_{j=1}^n \operatorname{gcd}(i, j) &= \sum_{k=1}^{n}k\sum_{i=1}^n \sum_{j=1}^{n}[\gcd(i,j)=k] \\[6pt]&= \sum_{k=1}^{n}k\sum_{i=1}^{n / k} \sum_{j=1}^{n / k}[\gcd(i,j)=1] \\[6pt]&= \sum_{k=1}^{n}k\sum_{i=1}^{n / k} \sum_{j=1}^{n / k}\sum_{d\mid \gcd(i,j)}\mu(d) \\[6pt]&= \sum_{k=1}^{n}k\sum_{d=1}^{n / k} \mu(d) \sum_{i=1}^{n /k}[d \mid i] \sum_{j=1}^{n / k}[d \mid j] \\[6pt]&= \sum_{k=1}^{n}k\sum_{d=1}^{n / k}\mu(d)\left\lfloor\frac{n}{k d}\right\rfloor^2 \end{aligned}$

然而直接搞是过不了的，因为数据组数太多：

$\mathcal{O}\left(\sqrt{\frac{n}{1}} + \sqrt{\frac{n}{2}}+\cdots+\sqrt{\frac{n}{n}} \right) = \mathcal{O}\left(\sqrt{n}\left(1 + \frac{1}{\sqrt{2}}+\cdots+\frac{1}{\sqrt{n}}\right)\right) = \mathcal{O}\left(n\right)$

考虑设 $m = n / k$ ，则

$\sum_{k=1}^{n}k\sum_{d=1}^{m}\mu(d)\left\lfloor\frac{m}{d}\right\rfloor^2$

后面的部分在 $m$ 固定时与 $k$ 无关，考虑预处理为 $f(m)$ ，则

$\sum_{k=1}^{n}k\cdot f\left(n/k\right)$

时间复杂度 $\mathcal{O}(N\sqrt{N})$ ，但是需要卡常，原理参考link。

代码：（用时 1.17s ）

cpp

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
namespace FastIO
{
    #define gc() readchar()
    #define pc(a) putchar(a)
    #define SIZ (int)(1e6+15)
    char buf1[SIZ],*p1,*p2;
    char readchar()
    {
        if(p1==p2)p1=buf1,p2=buf1+fread(buf1,1,SIZ,stdin);
        return p1==p2?EOF:*p1++;
    }
    template<typename T>void read(T &k)
    {
        char ch=gc();T x=0,f=1;
        while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
        while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();}
        k=x*f;
    }
    template<typename A, typename ...B>
        void read(A &x, B &...y) { return read(x), read(y...); }
    template<typename T>void write(T k)
    {
        if(k<0){k=-k;pc('-');}
        static T stk[66];T top=0;
        do{stk[top++]=k%10,k/=10;}while(k);
        while(top){pc(stk[--top]+'0');}
    }
}using namespace FastIO;
#define N ((int)(2e5 + 15))


bitset<N> ck; int pcnt,mu[N],f[N],v[N],prime[N];
void Mobius(int n = N - 10)
{
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!ck[i]) { prime[++pcnt] = i; mu[i] = -1; }
        for(int j = 1; j <= pcnt && i * prime[j] <= n; j++)
        {
            int pos = i * prime[j]; ck[pos] = true;
            if(i % prime[j]) { mu[pos] = -mu[i]; } else { mu[pos] = 0; break; }
        }
    }
    for(int i = 1; i <= n; i++) mu[i] += mu[i - 1];
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    Mobius();
    for(int i = 1; i <= N - 10; i++)
    {
        int sq = sqrt(i), cnt = 2 * sq;
        for(int j = 1; j <= sq; j++) { v[j] = j; v[cnt - j + 1] = i / j; }
        for(int k = 1, j = cnt; k <= cnt; k++,j--) f[i] += v[k] * v[k] * (mu[v[j]] - mu[v[j - 1]]);
    }
    for(int n; read(n), n != 0; )
    {
        int res = 0;
        for(int l = 1, r; l <= n; l = r + 1)
        {
            r = n / (n / l);
            res += (f[n / l] - 1) * (l + r) * (r - l + 1) / 2;
        }
        write(res / 2); pc('\n');
    }
    return 0;
}