loj6482. LJJ 爱数数

题意

给出\(n \leq {10} ^ {12}\)，求

\[\sum_{a = 1} ^ n \sum_{b = 1} ^ n \sum_{c = 1} ^ n [\frac{1}{a} + \frac{1}{b} = \frac{1}{c}] [\gcd(a, b, c) = 1] \]

题解

\[\begin{aligned} = & \sum_{a = 1} ^ n \sum_{b = 1} ^ n \sum_{c = 1} ^ n [c(a + b) = ab] [(a, b, c) = 1] \\ = & \sum_{a = 1} ^ n \sum_{b = 1} ^ n [(a + b) | ab] [(a, b, \frac{ab}{a + b}) = 1] \\ \end{aligned} \]

设\(g = (a, b)\)，则一对\((a, b)\)合法的充要条件\(a + b = g ^ 2\)。
则原问题相当于求

\[\begin{aligned} \sum_{g = 1} ^ {\sqrt {2n}} \sum_i [(g ^ 2 - i g, i g) = g] = & \sum_{g = 1} ^ {\sqrt {2n}} \sum_i [(g, i) = 1] \\ \end{aligned} \]

考虑\(i\)的上下界

\[1 \leq i g \leq n \\ 1 \leq g ^ 2 - i g \leq n \\ \]

得

\[i \in \Z \cap \left[ \max(g - \lfloor \frac{n}{g} \rfloor,1), \min(\lfloor \frac{n}{g} \rfloor, g - 1) \right] \]

即求

\[\sum_{g = 1} ^ {\sqrt{2n}} \sum_{i = \text{lowerbound}} ^ {\text{upperbound}} [(g, i) = 1] \]

考虑对后面的东西做前缀和是

\[f(L) = \sum_{i = 1} ^ L [(g, i) = 1] \]

反演一下可以得到

\[f(L) = \sum_{d | g} \mu(d) \lfloor \frac{L}{d} \rfloor \]

则原式即求

\[\sum_{g = 1} ^ {\sqrt {2n}} f(\text{upperbound}) - f(\text{lowerbound} - 1) \]

复杂度\(\mathcal O(\sqrt n \log n)\)。

#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 2e6 + 5, M = 4e7 + 5;
ll n, m, ans;
int p[N], mu[N], vis[N], d[M], tmp[N], cnt[N];
int calc (int g, int n) {
	int ret = 0;
	for (int i = cnt[g - 1] + 1; i <= cnt[g]; ++i) {
		ret += mu[d[i]] * (n / d[i]);
	}
    return ret;
}
int main () {
	cin >> n, mu[1] = 1;
	for (int i = 2; i < N; ++i) {
		if (!vis[i]) {
			p[++p[0]] = i, mu[i] = -1;
		}
		for (int j = 1; j <= p[0] && i * p[j] < N; ++j) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) {
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
	m = sqrt(n * 2);
	for (int i = 1; i <= m; ++i) {
		if (mu[i]) {
			for (int j = i; j <= m; j += i) {
				++cnt[j];
			}
		}
		cnt[i] += cnt[i - 1];
	}
	for (int i = 1; i <= m; ++i) {
		if (mu[i]) {
			for (int j = i; j <= m; j += i) {
				d[cnt[j - 1] + (++tmp[j])] = i;
			}
		}
	}
	for (int g = 1; g <= m; ++g) {
		ans += calc(g, min(n / g, 0ll + g - 1)) - calc(g, max(0ll + g - n / g, 1ll) - 1);
	}
	cout << ans << endl;
	return 0;
}