SP19985 GCDEX2 - GCD Extreme (hard) 题解

题目链接：GCDEX2 - GCD Extreme (hard)

题意：

给定 \(n~(1\le n\le 235711131719)\) ，计算
\[ G(n)=\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j) \] 答案对 \(2^{64}\) 取模。有多组询问 \((T\le10^4)\)。本题时限 20s 。

貌似写过这题的弱化版 UVA11424 GCD - Extreme (I) ，是一道 7 倍经验题。

不过之前的做法都不行了，现在我们来重新考察这个式子怎么做。

推式子，启动！ \[ \begin{aligned} \mathrm{Ans}&= \sum_{i=1}^n \sum_{j=i+1}^n \operatorname{gcd}(i, j) \\[6pt]& =\sum_{i=1}^n \sum_{j=1}^{i-1} \operatorname{gcd}(i, j) \\[6pt]& =\sum_{d=1}^n d \sum_{i=1}^n \sum_{j=1}^{i-1}[\operatorname{gcd}(i, j)=d] \\[6pt]& =\sum_{d=1}^n d \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j=1}^{i-1}[\operatorname{gcd}(i, j)=1] \\[6pt]& =\sum_{d=1}^n d \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} \varphi(i) \end{aligned} \] 现在我们需要快速求出 \(\varphi\) 的前缀和，考虑杜教筛即可，然后数论分块一下就完事了

注意本题对 unsigned long long （自然溢出）取模，有些地方需要特判一下。

时间复杂度 \(\mathcal{O}(n^{\frac{2}{3}})\)

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
typedef unsigned long long ull;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)(5e6 + 15))

char ck[N];
int pcnt, p[N]; ull phi[N];
unordered_map<int, ull> sum_phi;
void init(const int _ = N - 5)
{
    phi[1] = 1;
    rep(i, 2, _)
    {
        if(!ck[i]) { p[++pcnt] = i; phi[i] = i - 1; }
        for(int j = 1; j <= pcnt && i * p[j] <= _; j++)
        {
            const int pos = i * p[j]; ck[pos] = true;
            if(i % p[j] == 0) { phi[pos] = phi[i] * p[j]; break; }
            else phi[pos] = phi[i] * phi[p[j]];
        }
    }
    rep(i, 1, _) phi[i] += phi[i - 1];
}
ull sum(int n)
{
    if(n & 1) return 1ull * (n + 1) / 2 * n;
    else return 1ull * n / 2 * (n + 1);
}
ull get_sum(int n)
{
    if(n <= N - 5) return phi[n];
    if(sum_phi[n]) return sum_phi[n];
    ull res = sum(n);
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        res -= 1ull * (r - l + 1) * get_sum(n / l);
    }
    return sum_phi[n] = res;
}
void solve()
{
    int n; cin >> n; ull res = 0;
    for(int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        res += (sum(r) - sum(l - 1)) * (get_sum(n / l) - 1);
    }
    cout << res << '\n';
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int qwq; cin >> qwq; init(); while(qwq--) solve();
    return 0;
}

---

参考文献：

[1] https://www.luogu.com.cn/article/x33ds58f