【代码超详解】UVA 11426 GCD Extreme (II)（欧拉函数 + 思维）

一、题目描述

二、算法分析说明与代码编写指导

N ≤ 4000000，如果直接执行 N(N - 1) / 2 次 GCD，时间复杂度会达到 O(n^2 log n)，必 T 无疑。
注意到：N = 1 时 G = 0，N = 2 时 G = GCD(1, 2)，N = 3 时 G = GCD(1, 2) + GCD(1, 3) + GCD(2, 3)，N = 4 时 G = GCD(1, 2) + GCD(1, 3) + GCD(1, 4) + GCD(2, 3) + GCD(2, 4) + GCD(3, 4)，……。
因此写出递推关系：G(N) = G(N - 1) + f(N)，其中 f(N) = GCD(1, N) + GCD(2, N) + … + GCD(N - 1, N)。
于是问题转化为求解 f(N)。设 X < N。
注意到 GCD(X, N) = i，i | N。因为 N 较大时的因数个数远小于 N，所以 GCD(X, N) = i 肯定会出现当 X 取多个不同的值的时候 i 却相同的情况。那么对于每个 i，需要求出符合 GCD(X, N) = i 的 X 的数量，才可以计算 f(N)。如何求得呢？
GCD(X, N) = i，则 GCD(X / i，N / i) = 1，符合该等式的 X 数量不变。那么与 N / i 互质的 X 共有多少个呢？正是 φ(N / i) 个。不同的 i 的个数则是 N 的因数个数，所以可以得出
再根据 G(N) = G(N - 1) + f(N)，求得最后的答案为 G = G(N)。
在具体实现上，提供 2 种方法。
法一：对 N = 1，2，……，4000000，依次分解质因数并获得 N 的全部因子，然后对每个 N 遍历其每个因子，计算 f(N) 并累加。
当然，分解质因数的耗时还是有点大。好在这题时限为 10 s，这个方法足够通过。
法二：在枚举欧拉函数的值期间逐步求出 f(N)。


欧拉函数的模板如下：

template<class _Pty> inline void gen_phi(_Pty* const phi, const _Pty& n) {
	_Pty m = n / 2; phi[1] = 1;
	for (_Pty i = 2; i <= m; ++i) {
		if (!phi[i]) {
			for (_Pty j = i; j <= n; j += i) {
				if (!phi[j])phi[j] = j;
				phi[j] = phi[j] / i * (i - 1);
			}
		}
	}
	for (_Pty i = m + 1; i <= n; ++i)
		if (!phi[i]) { phi[i] = i - 1; }
}

我们在第一个 for 循环的 if 语句后面添加如下语句：

for (_Pty j = 1; j <= L; ++j)f[i * j] += j * phi[i];

然后在第二个 for 循环的 if 语句后面添加如下语句：

f[i] += phi[i];

注意新添加的语句不要写在 if 条件句内，否则会导致部分 i 对应的 f[i] 未进行计算，从而 WA。
我们比照
发现新添加的语句其实就是：f(ij) += jφ(i)，也就是说正是上面这条公式的一个计算步骤。对每一个 n = ij，i 总是从小到大的：i = 1，……，n；于是有对应的 j = n，……，1。也就是说在外层循环遍历到不同的 i 时，只要 i 是 n 的因数（ij = n），那么就相当于为计算 f(N) 进行了一次取了不同的因子 i 的累加。如此进行下去，直到最后，每个 f(N) 都会被计算完毕。
f[i] += φ(i) 其实就是取 j = 1 的特殊情形。毕竟这时候 i > n / 2，2i > n，j > 1 时一定有 i × j > n，越界。

三、AC 代码

法一：（4190 ms）

#include<cstdio>
#include<bitset>
#include<cmath>
#include<algorithm>
#include<vector>
#pragma warning(disable:4996)
using namespace std;
unsigned long long prime[300000], phi[4000001], _PrimeTy, MaxPrime, * prime_end = prime; bitset<4000001> notprime;
template<class _Pty> inline void gen_prime_phi(_Pty* const phi, const _Pty& n) {
	_Pty m = n / 2, L; notprime[0] = notprime[1] = true, phi[1] = 1;
	for (_Pty i = 2; i <= m; ++i) {
		if (!notprime[i]) { *prime_end = i, ++prime_end, phi[i] = i - 1; }
		L = n / i;
		for (auto j = prime; j != prime_end && *j <= L; ++j) {
			notprime[i * *j] = true;
			if (i % *j == 0) { phi[i * *j] = *j * phi[i]; break; }
			phi[i * *j] = (*j - 1) * phi[i];
		}
	}
	for (_Pty i = m + 1; i <= n; ++i)
		if (!notprime[i]) { *prime_end = i, ++prime_end, phi[i] = i - 1; }
	MaxPrime = *(prime_end - 1);
}
template<class _Ty> inline _Ty Power(_Ty radix, _Ty exp) {
	_Ty ans = 1;
	while (exp) {
		if (exp & 1)ans *= radix;
		exp >>= 1, radix *= radix;
	}
	return ans;
}
template<class _Ty> inline bool factor_enum(_Ty X, vector<_Ty>& dest) {//返回false代表数X太大，分解结果可能不正确。打表的最大质数的平方不能超过X的类型_Ty能够存储的最大值，否则计算出的分解上限L会溢出，导致分解成功与否判断出错。
	static _Ty L = MaxPrime * MaxPrime; static vector<_Ty> f, e, c; bool reliable;
	_Ty t = min((_Ty)sqrt(X), MaxPrime), D, q = 1; f.clear(), e.clear(); dest.clear();
	for (auto d = prime; *d <= t; ++d) {
		if (X % *d == 0) {
			X /= *d, f.emplace_back(*d), e.emplace_back(1);
			while (X % *d == 0) { X /= *d, ++* e.rbegin(); }
		}
		t = min((_Ty)sqrt(X), MaxPrime);
		if (X == 1) { reliable = true; break; }
	}
	if (X > 1 && X < L) { f.emplace_back(X), e.emplace_back(1), reliable = true; }
	else reliable = false;
	for (size_t i = 0; i < e.size(); ++i)q *= e[i] + 1;
	c.resize(e.size(), 0);
	for (_Ty h = 0; h != q; ++h) {
		D = 1;
		for (size_t i = 0; i != f.size(); ++i) { D *= Power(f[i], c[i]); }
		dest.emplace_back(D);
		for (size_t i = c.size() - 1; i != SIZE_MAX; --i) {
			++c[i];
			if (c[i] > e[i]) { c.pop_back(); continue; }
			else break;
		}
		while (c.size() != e.size())c.emplace_back(0);
	}
	return reliable;
}
unsigned long long n, s[4000001]; vector<unsigned long long> F;
inline unsigned long long f(const unsigned long long& i) {
	unsigned long long r = 0;
	for (size_t j = 0; j < F.size(); ++j) {
		r += F[j] * phi[i / F[j]];
	}
	return r - i;
}
int main() {
	gen_prime_phi(phi, 4000000ull);
	for (unsigned long long i = 2; i <= 4000000; ++i) { factor_enum(i, F); s[i] = s[i - 1] + f(i); }
	for (;;) {
		scanf("%llu", &n); if (n == 0)return 0;
		printf("%llu\n", s[n]);
	}
}

法二：（810 ms）

#include<cstdio>
#pragma warning(disable:4996)
unsigned long long phi[4000001]; unsigned long long n/*, s[4000001]*/, f[4000001];
template<class _Pty> inline void gen_phi(_Pty* const phi, const _Pty& n) {
	_Pty m = n / 2, L; phi[1] = 1;
	for (_Pty i = 2; i <= m; ++i) {
		if (!phi[i]) {
			for (_Pty j = i; j <= n; j += i) {
				if (!phi[j])phi[j] = j;
				phi[j] = phi[j] / i * (i - 1);
			}
		}
		L = n / i;
		for (_Pty j = 1; j <= L; ++j)f[i * j] += j * phi[i];
	}
	for (_Pty i = m + 1; i <= n; ++i) {
		if (!phi[i]) { phi[i] = i - 1; }
		f[i] += phi[i];
	}
}
int main() {
	gen_phi(phi, 4000000llu);
	//for (unsigned long long i = 2; i <= 4000000; ++i)s[i] = s[i - 1] + f[i];
	for (unsigned long long i = 2; i <= 4000000; ++i)f[i] += f[i - 1];
	for (;;) {
		scanf("%llu", &n); if (n == 0)return 0;
		//printf("%llu\n", s[n]);
		printf("%llu\n", f[n]);
	}
}