Version 1

UVA11417 GCD

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

UVA11424 GCD - Extreme (I)

SP3871 GCDEX - GCD Extreme

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

\[\begin{aligned} \sum_{i = 1}^n \sum_{j = i + 1}^n \gcd(i, j) & = \dfrac{\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n \gcd(i, j) - \sum\limits_{i = 1}^n \gcd(i, i)}{2} \\ & = \dfrac{\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n \gcd(i, j) - \sum\limits_{i = 1}^n i}{2} \end{aligned} \]

运用你丰富的莫反技巧，得到

\[\sum_{i = 1}^n \sum_{j = 1}^n \gcd(i, j) = \sum_{d = 1}^n \varphi(d) \left\lfloor\dfrac{n}{d}\right\rfloor^2 \]

线性筛 + 整除分块。

时间复杂度为 \(\Omicron(n + t \sqrt{n})\)。

会 \(\rm TLE\)。

---

换一种方法：

令 \(f(n) = \sum_{i = 1}^n \gcd(i, n)\)，那么

\[\begin{aligned} \sum_{i = 1}^n \sum_{j = i + 1}^n \gcd(i, j) & = \sum_{j = 1}^n \sum_{i = 1}^{j - 1} \gcd(i, j) \\ & = \sum_{i = 1}^n \sum_{j = 1}^i \gcd(i, j) - \sum_{i = 1}^n i \\ & = \sum_{i = 1}^n f(i) - \dfrac{n \cdot (n +1)}{2} \end{aligned} \]

又有

\[\begin{aligned} f(n) & = \sum_{i = 1}^n \gcd(i, n) \\ & = \sum_{d \mid n} \sum_{i = 1}^n d[\gcd(i, n) = d] \\ & = \sum_{d \mid n} d \sum_{i = 1}^{\frac{n}{d}} \left[\gcd\left(i, \dfrac{n}{d} \right) = 1 \right] \\ & = \sum_{d \mid n} d \cdot \varphi\left(\dfrac{n}{d} \right) \\ & = \sum_{d \mid n} \varphi(d) \cdot \dfrac{n}{d} \\ & = (\varphi * \operatorname{Id})(n) \end{aligned} \]

由于 \(\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

容斥。

设

\[\textsf{（不妨设 } n \le m \textsf{）} \begin{aligned} f(n, m) & = \sum_{i = 1}^n \sum_{j = 1}^m \gcd(i, j) \\ & = \sum_{i = 1}^n \varphi(i) \left\lfloor\dfrac{n}{i}\right\rfloor \left\lfloor\dfrac{m}{i}\right\rfloor \end{aligned} \]

则

\[ans = f(i2, j2) - f(i2, j1 - 1) - f(i1 - 1, j2) + f(i1 - 1, j1 - 1) \]

时间复杂度为 \(\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 = 1 + 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;
}