【题解】GCDEX 系列
Version 1
Description
-
多测，\(t\) 组数据。
-
给定整数 \(n\)，求
\[\sum_{i = 1}^n \sum_{j = i + 1}^n \gcd(i, j) \] -
\(1\le t\le 100, 1\le n\le 501\)。
Solution
\(\Omicron(n^2)\) 暴力。
Code
略。
Version 2
UVA11426 拿行李（极限版） GCD - Extreme (II)
Description
-
多测，\(t\) 组数据。
-
给定整数 \(n\)，求
\[\sum_{i = 1}^n \sum_{j = i +1}^n \gcd(i, j) \] -
\(1\le t\le 2\times 10^4, 1\le n\le 4\times 10^6\)。
Solution
运用你丰富的莫反技巧，得到
线性筛 + 整除分块。
时间复杂度为 \(\Omicron(n + t \sqrt{n})\)。
会 \(\rm TLE\)。
换一种方法：
令 \(f(n) = \sum_{i = 1}^n \gcd(i, n)\)，那么
又有
由于 \(\varphi, \operatorname{Id}\) 均为积性函数，因此 \(f\) 也为积性函数，可以线性筛（当然你可以 \(\Theta(n\ln n)\) 算但是这不够好
将 \(n\) 分解质因数：\(n = \prod_{i = 1}^k p_i^{\alpha_i}\)，记录 \(\beta(n) = \alpha_1, \gamma(n) = p_1^{\alpha_1}\)。
-
\(i\in \mathbb{P}\)：\(f(i) = i + (i - 1) = 2i - 1, \beta(i) = 1, \gamma(i) = i\)
-
\(p_j\nmid i\)：\(f(i\cdot p_j) = f(i)\cdot f(p_j), \beta(i\cdot p_j) = 1, \gamma(i\cdot p_j) = p_j\)
-
\(p_j\mid i\)：\(\beta(i\cdot p_j) = \beta(i) + 1, \gamma(i\cdot p_j) = \gamma(i) \cdot p_j\)
-
\(\gamma(i) = i\)，即 \(i\cdot p_j = p^k(p\in \mathbb{P})\)，其中 \(k = \beta(i \cdot p_j)\)，暴力算
\[\begin{aligned} f(p^k) & = \varphi(1) \cdot \dfrac{p^k}{1} + \sum_{d\mid p^k, d > 1} \varphi(d) \cdot \dfrac{p^k}{d} \\ & = p^k \cdot \sum_{i = 1}^k [p^{i - 1}(p - 1)] \cdot p^{k - i} \\ & = p^k \cdot \sum_{i = 1}^k (p^i - p^{i - 1}) \cdot p^{k - i} \\ & = p^k \cdot \sum_{i = 1}^k p^k - p^{k - 1} \\ & = (k + 1) p^k - k \cdot p^{k - 1} \\ & = [\beta(i\cdot p_j) + 1](i\cdot p_j) - \beta(i\cdot p_j)\cdot i \end{aligned} \] -
\(\gamma(i) \ne i\)，此时 \(p_j\) 是 \(i\) 的最小质因数。
\[f(i \cdot p_j) = f\left\{\dfrac{i}{\gamma(i)} \cdot [p_j \cdot \gamma(i)] \right\} \]此时 \(i\) 的所有质因数 \(p_j\) 都被除掉了，\(\gcd\left[\dfrac{i}{\gamma(i)}, p_j \cdot \gamma(i) \right] = 1\)，由 \(f\) 是积性函数有
\[f(i \cdot p_j) = f\left\{\dfrac{i}{\gamma(i)} \cdot [p_j \cdot \gamma(i)] \right\} = f\left[\dfrac{i}{\gamma(i)} \right] \cdot f[p_j \cdot \gamma(i)] \]
-
综上，预处理 \(\Omicron(n)\)，询问 \(\Theta(1)\)，总时间复杂度为 \(\Omicron(n + t)\)。
Code
// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
const int MAXN = 1e6 + 5;
const int N = 1e6;
int p[MAXN], beta[MAXN], gamma[MAXN];
ll f[MAXN], sum[MAXN];
bool vis[MAXN];
void pre()
{
f[1] = 1;
for (int i = 2; i <= N; i++)
{
if (!vis[i])
{
p[++p[0]] = i;
f[i] = 2 * i - 1;
beta[i] = 1;
gamma[i] = i;
}
for (int j = 1; j <= p[0] && i * p[j] <= N; j++)
{
vis[i * p[j]] = true;
if (i % p[j] == 0)
{
beta[i * p[j]] = beta[i] + 1;
gamma[i * p[j]] = gamma[i] * p[j];
if (gamma[i] == i)
{
f[i * p[j]] = (ll)(beta[i * p[j]] + 1) * (i * p[j]) - beta[i * p[j]] * i;
}
else
{
f[i * p[j]] = f[i / gamma[i]] * f[p[j] * gamma[i]];
}
break;
}
f[i * p[j]] = f[i] * f[p[j]];
beta[i * p[j]] = 1;
gamma[i * p[j]] = p[j];
}
}
for (int i = 1; i <= N; i++)
{
sum[i] = sum[i - 1] + f[i];
}
}
int main()
{
pre();
int n;
while (scanf("%d", &n), n)
{
printf("%lld\n", sum[n] - (ll)n * (n + 1) / 2);
}
return 0;
}
Version 3
SP26045 GCDMAT2 - GCD OF MATRIX (hard)
Description
-
多测，\(t\) 组数据。
-
给定正整数 \(i1, i2, j1, j2\)，求
\[\begin{aligned} \left[\sum_{i = i1}^{i2} \sum_{j = j1}^{j2} \gcd(i, j) \right] \bmod (10^9 +7) \end{aligned} \] -
\(t\le 5\times 10^4, \max\{i2\}, \max\{j2\} \le 10^6\)。
Solution
容斥。
设
则
时间复杂度为 \(\Omicron(n + t\sqrt{n})\)。
你算了算 \(10^6 + 5\times 10^4 \times 2\sqrt{10^6} \approx 10^8\)，而时限 \(2.5s\)，好像能过。
问题来了，你的整除分块做除法至少有 \(4\) 倍常数，而一次容斥又要做 \(4\) 次整除分块，总共是 \(16\) 倍大常数。
\(10^6 + 5\times 10^4 \times 2\sqrt{10^6} \times 16 \approx 1.6\times 10^9\)。
额，终于知道为什么评黑了。
开始卡常。
-
\(\text{Step } 1\)：不用容斥，直接算
\[ans = \sum_{i = 1}^{\min(i2, j2)} \varphi(i) \left(\left\lfloor\dfrac{i2}{i}\right\rfloor - \left\lfloor\dfrac{i1 - 1}{i}\right\rfloor \right) \left(\left\lfloor\dfrac{j2}{i}\right\rfloor - \left\lfloor\dfrac{j1 - 1}{i}\right\rfloor \right) \]这样整除分块变成 \(8\) 倍常数，但只用算一次。
-
\(\text{Step } 2\)：预处理 \(1\sim 10^6\) 的倒数，用
double存应该够了。乘法比除法快。 -
\(\text{Step } 3\)：根号分治。
Code
代码中有一步 inv[i] = (1.0 / i) * EXP，其中 \(EXP = 10^{-10}\)。
原因是像 \(i = 3\) 这样的会得到 \(inv[3] = 0.333\cdots\)，然后 \(3 \times inv[3] = 0.999\cdots\)，下取整成了 \(0\)，所以要乘上 \(EXP\) 防止精度误差。
//18 = 9 + 9 = 18.
#pragma GCC optimize("Ofast")
#include <iostream>
#include <cstdio>
#include <bitset>
#include <cmath>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
#define re register
#define il inline
namespace Fread
{
}
using Fread::read;
namespace Fwrite
{
}
using Fwrite::write;
const int MAXN = 1e6 + 5;
const int MOD = 1e9 + 7;
const double EXP = 1 + 1e-10;
int p[MAXN / 10], phi[MAXN], sum[MAXN];
bitset<MAXN> vis;
double inv[MAXN];
void pre(int n)
{
phi[1] = 1;
for (int i = 2; i <= n; i++)
{
if (!vis[i])
{
p[++p[0]] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= p[0] && i * p[j] <= n; j++)
{
vis[i * p[j]] = true;
if (i % p[j] == 0)
{
phi[i * p[j]] = phi[i] * p[j];
break;
}
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
for (int i = 1; i <= n; i++)
{
sum[i] = (sum[i - 1] + phi[i]) % MOD;
}
for (int i = 1; i <= n; i++)
{
inv[i] = (1.0 / i) * EXP;
}
}
il int GetSum(int l, int r)
{
return (sum[r] - sum[l - 1] + MOD) % MOD;
}
il int min2(int x, int y)
{
return x < y ? x : y;
}
il int div2(int a, int b)
{
if (b)
{
return (int)a * inv[b];
}
return 0x3f3f3f3f;
}
il void swap2(int &x, int &y)
{
x ^= y ^= x ^= y;
}
int block(int i1, int i2, int j1, int j2)
{
i1--, j1--;
if (i2 > j2)
{
swap2(i1, j1);
swap2(i2, j2);
}
int mid = sqrt(i2), res = 0;
for (int i = 1; i <= mid; i++)
{
res = (res + (ll)phi[i] * (div2(i2, i) - div2(i1, i)) % MOD * (div2(j2, i) - div2(j1, i)) % MOD) % MOD;
}
for (int l = mid + 1, r; l <= i2; l = r + 1)
{
int k1 = div2(i1, l), k2 = div2(i2, l), k3 = div2(j1, l), k4 = div2(j2, l);
r = min2(div2(i1, k1), min2(div2(i2, k2), min2(div2(j1, k3), div2(j2, k4))));
res = (res + (ll)GetSum(l, r) * (k2 - k1) % MOD * (k4 - k3) % MOD) % MOD;
}
return res;
}
int main()
{
int t = read(), n = read(), m = read();
pre(max(n, m));
while (t--)
{
int i1 = read(), j1 = read(), i2 = read(), j2 = read();
write(block(i1, i2, j1, j2), '\n');
}
return 0;
}
Version4
SP19985 GCDEX2 - GCD Extreme (hard)
Description
-
多测，\(t\) 组数据。
-
给定整数 \(n\)，求出
\[\left[\sum_{i = 1}^n \sum_{j = i + 1}^n \gcd(i, j) \right] \bmod 2^{64} \] -
\(t\le 10^4\)，\(1\le n\le 235711131719 \approx 2\times 10^{11}\)。
Solution
首先 \(\bmod 2^{64}\) 就是开 \(\text{unsigned long long}\) 然后自然溢出。
Version 2 的方法要 \(\Omicron(n)\) 预处理，会废。
试着枚举 \(\gcd\)：
前面一项把杜教筛扔到整除分块里就 \(\Omicron(n^{\frac{2}{3}} + \sqrt{n}) \simeq \Omicron(n^{\frac{2}{3}})\) 做完了。
别紧张时限有 \(20s\)。
注意在 Version 2 中不要用杜教筛因为 \(n\) 太小然后杜教筛常数太大。
Code
// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <unordered_map>
#define Debug(x) cout << #x << "=" << x << endl
#define int unsigned long long
using namespace std;
const int MAXN = 21544346 + 5;
const int N = 21544346;
int p[MAXN], phi[MAXN], sum[MAXN];
bool vis[MAXN];
void pre()
{
phi[1] = sum[1] = 1;
for (int i = 2; i <= N; i++)
{
if (!vis[i])
{
p[++p[0]] = i;
phi[i] = i - 1;
}
sum[i] = sum[i - 1] + phi[i];
for (int j = 1; j <= p[0] && i * p[j] <= N; j++)
{
vis[i * p[j]] = true;
if (i % p[j] == 0)
{
phi[i * p[j]] = phi[i] * p[j];
break;
}
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
}
int S(int n)
{
if (n % 2 == 0)
{
return n / 2 * (n + 1);
}
return (n + 1) / 2 * n;
}
unordered_map<int, int> dp;
int sublinear(int n)
{
if (n <= N)
{
return sum[n];
}
if (dp.count(n))
{
return dp[n];
}
int res = S(n);
for (int l = 2, r; l <= n; l = r + 1)
{
int k = n / l;
r = n / k;
res -= (r - l + 1) * sublinear(k);
}
return dp[n] = res;
}
int GetSum(int l, int r)
{
return S(r) - S(l - 1);
}
int block(int n)
{
int res = 0;
for (int l = 1, r; l <= n; l = r + 1)
{
int k = n / l;
r = n / k;
res += GetSum(l, r) * sublinear(k);
}
return res;
}
signed main()
{
pre();
int t;
scanf("%llu", &t);
while (t--)
{
int n;
scanf("%llu", &n);
printf("%llu\n", block(n) - S(n));
}
return 0;
}
