HDU7201 Yet Another Easy Function Sum Problem

题意：求

$$\sum_{i = 1} ^ n \sum_{j = 1} ^ n H(ij) [\gcd(i, j) = 1]$$

其中 $H(x)$ 表示 $x$ 所有质因子的乘积（去重后）。$n\leq 10 ^ 9$，$T(T\leq 5)$ 组数据，满足 $\sum n\leq 2\times 10 ^ 9$。20s。

类似 $\varphi$ 的，我们有一个这样的结论：

$$H(ij) = \dfrac {H(i) H(j)}{H(\gcd(i, j))}$$

对每个质因数分开看就好了。

首先看到 $\gcd$ 套路反演一下（所有除法默认下取整）：

$$\begin{aligned} & \sum_{i = 1} ^ n \sum_{j = 1} ^ n H(ij) [\gcd(i, j) = 1] \\ =& \sum_{i = 1} ^ n \sum_{j = 1} ^ n H(i)H(j) \sum_{d | i, d | j} \mu(d)\\ =& \sum_{d = 1} ^ n \mu(d)(\sum_{i = 1} ^ {\frac nd} H(id)) ^ 2 \end{aligned}$$

这个看样子就不能很好做，至少都得枚举 $d$，这样就不可能低于 $O(n)$。考虑对 $d$ 分治，钦定阈值 $B$，$d\leq B$ 就暴力枚举，否则我们换另外的方法。

1. $d\leq B$

记 $F_1(n, d) = \sum_{i = 1} ^ n H(id)$，那么原式就可以表示成 $\sum_{d = 1} ^ B \mu(d)F_1 ^ 2(\dfrac nd, d)$。考虑如何计算 $F_1(n, d)$。

注意到有一个 $H(id)$ 很烦，使用上面的套路，我们开始展开：

$$\begin{aligned} F_1(n, d) =& \sum_{i = 1} ^ n H(id)\\ =& \sum_{k | d} \sum_{i = 1} ^ n [\gcd(i, d) = k] \dfrac {H(i)H(d)}{H(k)}\\ =& \sum_{k | d} \dfrac 1{H(k)} \sum_{i = 1} ^ {\frac nk} [\gcd(i, \dfrac dk) = 1] H(ik)H(d)\\ =& \sum_{k | d} \dfrac 1{H(k)}\sum_{i = 1} ^ {\frac nk} H(ik)H(d) \sum_{t | i, tk | d} \mu(t)\\ =& \sum_{T | d}(\sum_{k | T} \dfrac {\mu(\frac Tk)}{H(k)}) (\sum_{i = 1} ^ {\frac nT} H(iT))\\ =& \sum_{T | d} H_1(T) F_1(\dfrac nT, T) \end{aligned}$$

这里令 $H_1 = H \times \mu$，表示狄利克雷卷积。容易发现这是一个简单的递归函数，直接递归并记忆化即可。

注意边界是 $n = 0$（答案显然是 0）或者是 $d = 1$，此时我们相当于是求 $\sum_{i = 1} ^ n H(i)$。这个东西直接使用 Min_25 筛去处理，时间复杂度 $O(\dfrac {n ^ {\frac 34}}{\log n})$。

不考虑 Min_25 筛的复杂度，我们考虑所有可能的 $F_1(n, d)$ 和可能的转移数。容易发现一个 $F_1(n, d)$ 可能的转移数为 $O(\sigma(d))$。对于从外部调用的 $F_1(n, d)$，就是 $O(\sum_{i = 1} ^ B \sigma(i)) = O(B\log B)$。对于自调用的的 $F_1(m, t)$，一定满足 $m = \dfrac n{t ^ 2x}$，因为前面一个假设不来自外部调用的话，那么上一级的 $d$ 需要满足 $t|d$，我们可以分析得到 $O(\sum_{d = 1} ^ B \dfrac{\sqrt n}d) = O(\sqrt n\log n)$。类似的分析转移数可以得到 $O(\sqrt n\log ^ 2 n)$。总复杂度 $O(B\log n + \sqrt n\log ^ 2n)$，后半部分不管了（

2. $d > B$

再看 $d$ 显然是不明智的选择，我们考虑把平方看作枚举 $i, j$，计算 $\sum_{d = 1} [id\leq n\land jd\leq n]H(id)H(jd)$。令 $G(m, i, j) = \sum_{d = 1} ^ m \mu(d) H(id)H(jd)$。那么就可以得到原式为：

$$\sum_{i = 1} ^ {\frac nB} \sum_{j = 1} ^ {\frac nB} (G(\dfrac n{\max(i, j)}, i, j) - G(B, i, j))$$

考虑如何计算 $G(m, i, j)$。仍然套路展开 $H(id)$ 和 $H(jd)$，大力计算：

$$\begin{aligned} & \sum_{d = 1} ^ m \mu(d) H(id) H(jd)\\ =& H(i)H(j)\sum_{d = 1} ^ m \mu(d) H ^ 2(d) \sum_{k_1 | i} [\gcd(i, d) = k_1] \sum_{k_2 | j} [\gcd(j, d) = k_2] \dfrac 1{H(k_1)H(k_2)}\\ =& H(i) H(j) \sum_{k_1 | i} \sum_{k_2 | j} \dfrac 1{H(k_1)H(k_2)} \sum_{d = 1} ^ {\frac mk} \mu(dk) H ^ 2(dk) [\gcd(\dfrac i{k_1}, d) = 1][\gcd(\dfrac i{k_2}, d) = 1]\\ =& H(i) H(j) \sum_{k_1 | i} \sum_{k_2 | j} \dfrac 1{H(k_1) H(k_2)} \sum_{d = 1} ^ {\frac mk} \mu(dk) H ^ 2(dk) \sum_{t_1 | d, t_1k_1 | i} \mu(t_1) \sum_{t_2 | d, t_2k_2 | j} \mu(t_2) \\ =& H(i)H(j) \sum_{T_1 | i} \sum_{T_2 | j} (\sum_{k_1 | T_1} \dfrac {\mu(\frac {T_1}{k_1})}{H(k_1)})(\sum_{k_2 | T_2} \dfrac {\mu(\frac {T_2}{k_2})}{H(k_2)}) \sum_{d = 1} ^ {\frac mT} \mu(dT) H ^ 2(dT) \end{aligned}$$

上述式子中，$T = \text{lcm} (T_1, T_2)$，$k = \text{lcm} (k_1, k_2)$，最后一个仍然很难算，我们考虑再设一个函数 $F_2$ 来表示。令 $F_2(n, d)$ 为 $\sum_{i = 1} ^ n \mu(di) H ^ 2(di)$。这个和 $F_1$ 是比较类似的，我们就不做每一步的推导了，直接给结论：

$$F_2(n, d) = \mu(d) H ^ 2(T) \sum_{T | d} \mu(k) F_2(\dfrac nT, T)$$

同样，$F_2$ 的复杂度和 $F_1$ 分析的类似。不是很会分析了 /kk，据说是 $\tilde O(\dfrac {n}{\sqrt B})$。

平衡规划与代码

两边平衡一下，取 $B = n ^ {\frac 23}$，可以得到复杂度约为 $O(\dfrac {n ^ {\frac 34}}{\log n} + n ^ {\frac 23} \log n)$。卡卡常大概就可以过了（？）。

说几点代码细节：

1. 容易发现递归 $F_2$ 的时候 $\mu(t) = 0$ 可以直接返回。
2. $B$ 必须是 $n$ 的约数才能过，原因好像是交界处可能寄。
3. Min_25 可以写迭代，比较少见呢。

void sieve()
{
	H[1] = mu[1] = 1;
	for (int i = 2; i < N; ++ i)
	{
		if (!st[i]) prime[cnt ++] = i, H[i] = i, mu[i] = -1;
		for (int j = 0; j < cnt && i * prime[j] < N; ++ j) {
			H[i * prime[j]] = H[i];
			st[i * prime[j]] = true;
			if (i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i], H[i * prime[j]] *= prime[j];
		}
	}
	for (int i = 1; i < N; ++ i) iH[i] = qpow(H[i]);
	for (int i = 1; i < N; ++ i) {
		if (!mu[i]) continue;
		if (mu[i] > 0) {
			for (int j = 1; i * j < N; ++ j)
				adj(H1[i * j] += iH[j] - Mod);
		} else {
			for (int j = 1; i * j < N; ++ j)
				adj(H1[i * j] -= iH[j]);
		}
	}
	for (int i = 1; i < N; ++ i)
		for (int j = 1; i * j < N; ++ j)
			if (mu[i * j]) fac[i * j].push_back(i);
}

inline int getid(int x) { return x <= sq ? id1[x] : id2[n / x]; }
inline int loc(int x) { return std::lower_bound(a + 1, a + tot + 1, x, std::greater<int>()) - a; }
inline int sum1(int x) { return x * (x + 1LL) / 2 % Mod; }
inline int sum2(int x) { return x * (x + 1LL) % Mod * (x << 1 | 1) % Mod * (Mod + 1) / 6 % Mod; }

void Min25_init()
{
	sq = std::sqrt(n), tot = 0;
	for (int l = 1, r, t; l <= n; l = r + 1)
	{
		t = n / l, r = n / t;
		a[++ tot] = t;
		if (t <= sq) id1[t] = tot;
		else id2[r] = tot;
	}
	
	for (int i = 0; i < cnt; ++ i)
		adj(sp1[i] = sp1[i - 1] + prime[i] - Mod), sp2[i] = (sp2[i - 1] + (LL) prime[i] * prime[i]) % Mod;
	for (int i = 1; i <= tot; ++ i) psum1[i] = sum1(a[i]) - 1, psum2[i] = sum2(a[i]) - 1;
	for (int i = 0; i < cnt; ++ i) {
		LL le = (LL) prime[i] * prime[i];
		for (int j = 1, p = prime[i]; j <= tot && a[j] >= le; ++ j)
			psum1[j] = (psum1[j] + (LL) (sp1[i - 1] + Mod - psum1[getid(a[j] / p)]) * p) % Mod,
			psum2[j] = (psum2[j] + (LL) (sp2[i - 1] + Mod - psum2[getid(a[j] / p)]) * p * p) % Mod;
	}
	
	for (int j = 1, x; j <= tot; ++ j)
		adj(f1[j] = psum1[j] - psum1[x = loc(std::sqrt(a[j]))]), adj(f2[j] = psum2[x] - psum2[j]);
	for (int i = cnt - 1; ~i; -- i)
		for (int j = 1, p = prime[i]; j <= tot && a[j] >= (LL) p * p; ++ j) {
			LL cur = p;
			int v, nxt;
			for (int e = 1; cur * p <= a[j]; ++ e, cur *= p) {
				int nxt = a[j] / cur, v;
				if (nxt < p * p) adj(v = psum1[getid(nxt)] - sp1[i]);
				else v = f1[getid(nxt)];
				f1[j] = (f1[j] + (v + 1LL) * p) % Mod;
			}
			if ((nxt = a[j] / p) < p * p) adj(v = sp2[i] - psum2[getid(nxt)]);
			else v = f2[getid(nxt)];
			adj(f1[j] += p - Mod);
			f2[j] = (f2[j] + (Mod - p * p) * (v + 1LL)) % Mod;
		}
	for (int j = 1; j <= tot; ++ j) f1[j] ++, f2[j] ++;
}

inline int getF1(int n, int d)
{
	if (!n) return 0;
	if (n == 1) return H[d];
	if (d == 1) return f1[getid(n)];
	if (mf1[d].count(n)) return mf1[d][n];
	int res = 0;
	for (int T : fac[d])
		res = (res + (LL) H1[T] * getF1(n / T, T)) % Mod;
	return mf1[d][n] = (LL) res * H[d] % Mod;
}

inline int getF2(int n, int d)
{
	if (!n) return 0;
	if (d == 1) return f2[getid(n)];
	if (!mu[d]) return 0;
	if (mf2[d].count(n)) return mf2[d][n];
	int res = 0;
	for (int k : fac[d])
		if (mu[k] > 0) adj(res += getF2(n / k, k) - Mod);
		else if (mu[k]) adj(res -= getF2(n / k, k));
	res = (LL) res * H[d] % Mod * H[d] % Mod;
	if (mu[d] < 0) adj(res = -res);
	return mf2[d][n] = res;
}

void work()
{
	std::cin >> n;
	Min25_init();
	B = n / (int) (std::pow(n, 1. / 3) + 1), B = std::max(B, 1);
	int res = 0;
	for (int i = 1, v; i <= B; ++ i)
		if (mu[i] > 0) v = getF1(n / i, i), res = (res + (LL) v * v) % Mod;
		else if (mu[i]) v = getF1(n / i, i), res = (res + (LL) (Mod - v) * v) % Mod;
	int lim = n / B;
	for (int T1 = 1; T1 <= lim; ++ T1)
		for (int T2 = 1; T2 <= lim; ++ T2) {
			LL Lcm = (LL) T1 / Gcd(T1, T2) * T2, mul = (LL) H1[T1] * H1[T2] % Mod, tmp;
			if (!mu[Lcm]) continue;
			for (int i = T1, ed = lim; i <= ed; i += T1)
				for (int j = T2; j <= ed; j += T2)
				{
					tmp = mul * H[i] % Mod * H[j] % Mod;
					res = (res + tmp * getF2(n / std::max(i, j) / Lcm, Lcm)) % Mod;
					if (B >= Lcm)
						res = (res + tmp * (Mod - getF2(B / Lcm, Lcm))) % Mod;
				}
		}
	std::cout << res << '\n';
}