【题解】GCDEX 系列

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

$\mathrm{O}(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{array}{rl}\sum _{i=1}^n\sum _{j=i+1}^n\gcd (i,j)& =\frac{\sum _{i=1}^n\sum _{j=1}^n\gcd (i,j)-\sum _{i=1}^n\gcd (i,i)}{2}\\ & =\frac{\sum _{i=1}^n\sum _{j=1}^n\gcd (i,j)-\sum _{i=1}^{n}i}{2}\end{array}$$

运用你丰富的莫反技巧，得到
$$\sum _{i=1}^n\sum _{j=1}^n\gcd (i,j)=\sum _{d=1}^n\phi (d){\lfloor \frac{n}{d}\rfloor }^2$$
线性筛 + 整除分块。

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

会 $\mathrm{T}\mathrm{L}\mathrm{E}$。

---

换一种方法：

令 $f(n)={\sum }_{i=1}^n\gcd (i,n)$，那么
$$\begin{array}{rl}\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)-\frac{n\cdot (n+1)}{2}\end{array}$$
又有
$$\begin{array}{rl}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,\frac{n}{d}\right)=1\right]\\ & =\sum _{d\mid n}d\cdot \phi \left(\frac{n}{d}\right)\\ & =\sum _{d\mid n}\phi (d)\cdot \frac{n}{d}\\ & =(\phi \ast \mathrm{Id})(n)\end{array}$$
由于 $\phi ,\mathrm{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{array}{rl}f(p^k)& =\phi (1)\cdot \frac{p^k}{1}+\sum _{d\mid p^k,d>1}\phi (d)\cdot \frac{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{array}$$

  • $\gamma (i)\ne i$，此时 $p_j$ 是 $i$ 的最小质因数。
    $$f(i\cdot p_j)=f\left\{\frac{i}{\gamma (i)}\cdot [p_j\cdot \gamma (i)]\right\}$$
    此时 $i$ 的所有质因数 $p_j$ 都被除掉了，$\gcd \left[\frac{i}{\gamma (i)},p_j\cdot \gamma (i)\right]=1$，由 $f$ 是积性函数有
    $$f(i\cdot p_j)=f\left\{\frac{i}{\gamma (i)}\cdot [p_j\cdot \gamma (i)]\right\}=f\left[\frac{i}{\gamma (i)}\right]\cdot f[p_j\cdot \gamma (i)]$$

综上，预处理 $\mathrm{O}(n)$，询问 $\Theta (1)$，总时间复杂度为 $\mathrm{O}(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{array}{r}\left[\sum _{i=i1}^{i2}\sum _{j=j1}^{j2}\gcd (i,j)\right]\mathrm{m}\mathrm{o}\mathrm{d}(10^9+7)\end{array}$$

• t ≤ 5 × 1 0 4 , max ⁡ { i 2 } , max ⁡ { j 2 } ≤ 1 0 6 t\le 5\times 10^4, \max\{i2\}, \max\{j2\} \le 10^6 t≤5×104,max{i2},max{j2}≤106。

Solution

容斥。

设
$$\text{（不妨设 }n\le m\text{）}\begin{array}{rl}f(n,m)& =\sum _{i=1}^n\sum _{j=1}^m\gcd (i,j)\\ & =\sum _{i=1}^n\phi (i)\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor \end{array}$$

则
$$ans=f(i2,j2)-f(i2,j1-1)-f(i1-1,j2)+f(i1-1,j1-1)$$
时间复杂度为 $\mathrm{O}(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)}\phi (i)\left(\lfloor \frac{i2}{i}\rfloor -\lfloor \frac{i1-1}{i}\rfloor \right)\left(\lfloor \frac{j2}{i}\rfloor -\lfloor \frac{j1-1}{i}\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]\mathrm{m}\mathrm{o}\mathrm{d}{2}^{64}$$

• $t\le 10^4$，$1\le n\le 235711131719\approx 2\times 10^{11}$。

Solution

首先 $\mathrm{m}\mathrm{o}\mathrm{d}{2}^{64}$ 就是开 $\text{unsigned long long}$ 然后自然溢出。

Version 2 的方法要 $\mathrm{O}(n)$ 预处理，会废。

---

试着枚举 $\gcd$：
$$\begin{array}{rl}ans& =\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 _{d=1}^n\sum _{i=1}^n\sum _{j=1}^{i}d[\gcd (i,j)=d]-\sum _{i=1}^{n}i\\ & =\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^i[\gcd (i,j)=1]-\sum _{i=1}^{n}i\\ & =\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (i)-\sum _{i=1}^{n}i\end{array}$$
前面一项把杜教筛扔到整除分块里就 $\mathrm{O}(n^{\frac{2}{3}}+\sqrt{n})\simeq \mathrm{O}(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;
}