二进制 GCD 学习笔记

前言

欧几里得算法可以在 log 的时间复杂度内求出 个数的 GCD，但是这还是太慢了。

在一些题目中 ，欧几里得算法就会 TLE。

欧几里得算法

理论：\(\gcd(a,b) = \gcd(b, a \bmod b)\)

二进制 GCD

更相减损术

已知两个数 \(a\), \(b\), 求 \(\gcd(a,b)\)。

设 \(a \ge b\)，若 \(a=b\)，则 \(\gcd(a,b) = a = b\)，否则，对于 \(\forall d \mid a, d \mid b\)，可以证明 \(d \mid a - b\)。

\(\because d \mid a\) \(\therefore a = dk_1\)。同理，\(b = dk_2\)

\(\because a - b = dk_1 - dk_2 = d(k_1-k_2)\)

\(\therefore d \mid a-b\)

因此，\(a\) 和 \(b\) 的公因数都是 \(a-b\) 和 \(b\) 的。即：\(\gcd(a,b) = \gcd(a-b, b)\)

原理

假设我们要求的 \(\gcd(a,b)\) 中，满足 \(a, b > 0\)

则有以下 \(5\) 种情况：

第一种：\(a = b\)，则有 \(\gcd(a,b) = a = b\)

第二种：\(a = 0\) 或 \(b = 0\)，若 \(a=0\)，则 \(\gcd(a,b) = \gcd(0, b) = b\)，若 \(b=0\)，则 \(\gcd(a,b) = \gcd(a,0) = a\)

第三种：\(a\) 和 \(b\) 均是偶数。很明显，此时 \(a\) 和 \(b\) 有公因数 \(2\)，所以可以将 \(a,b\) 同时除以 \(2\)，再求 GCD，用公式来表达，就是 \(\gcd(a,b) = \gcd(\frac{a}{2}, \frac{b}{2})\)，等式右边乘以 \(2\)，是因为我们将 \(a,b\) 都除以了 \(2\)，而没有算入答案，所以要乘 \(2\)。

第四种：\(a\) 和 \(b\) 中有一个偶数。不妨设 \(a\) 就是那个偶数，显然，\(2\) 不是 \(a\) 和 \(b\) 的公因数，所以我们直接把 \(a\) 除以 \(2\) 就可以了。表达式为：\(\gcd(a,b) = \gcd(\frac{a}{2}, b)\)

第五种：\(a\) 和 \(b\) 都是奇数，设 \(a > b\), 则有 \(\gcd(a,b) = \gcd(\frac{a-b}{2}, b)\)。下面是证明。

证明

设 \(a = \gcd(a,b) \times k_1, b = \gcd(a,b) \times k_2\)。

\(\because k_1\) 表示的数字和 \(b\) 互质，\(k_2\) 表示的数字和 \(a\) 互质

\(\therefore \gcd(k_1, k_2) = 1\)

接下来证明 \(\gcd(a - b, b) = \gcd(a, b)\)（这里跟前面的证明方法不一样）

\(\because a = k_1 \times \gcd(a,b), b = k_2 \times \gcd(a, b)\)

\(\therefore a - b = (k_1-k_2) \times \gcd(a, b)\)，又 \(\because b = k_2 \times \gcd(a, b)\)

\(\therefore \gcd(a-b, b) = \gcd(k_1 - k_2, k_2) \times \gcd(a, b)\)

接下来需要证明 \(\gcd(k_1 - k_2, k_2) = 1\)，下面采用 反证法。

设 \(\gcd(k_1 - k_2, k_2) = m, (m \in \mathbb{N^*}, m > 1), k_1 - k_2 = t_1m, k2 = t_2m\)

则 \(k_1 = (t_1 + t_2)m\)

\(\therefore \gcd(k_1, k_2) = m \times \gcd(t_1 + t_2, t_2) > 1\)

与 \(\gcd(k_1, k_2) = 1\) 矛盾，所以 \(\gcd(k_1 - k_2, k_2) = 1\)

\(\because \gcd(k_1 - k_2, k_2) = 1\)

\(\therefore \gcd(a, b) = \gcd(a - b, b)\)

综上所述，所以 \(\gcd(a, b) = \gcd(a-b, b)\)，证毕。

于是，第 \(5\) 种情况可以转换为 \(\gcd(a, b) = \gcd(a - b, b)\)

又因为 \(a, b\) 均为奇数，所以 \(a - b\) 为偶数，于是，第 \(5\) 种情况就能被转化为第 \(4\) 种情况

所以 \(\gcd(a-b, b) = \gcd(\frac{a-b}{2}, b)\)，即 \(\gcd(a, b) = \gcd(\frac{a-b}{2}, b)\)

通过以上的推理，我们可以将 GCD 转化为

int gcd(int a, int b) {
	if (a == b) return a; // 1
	if (a == 0) return b; // 2
	if (b == 0) return a; // 2
	if (!(a & 1) && !(b & 1)) // 3
		return gcd(a >> 1, b >> 1) << 1;
	if (!(a & 1) && (b & 1)) // 4
		return gcd(a >> 1, b);
	if (!(b & 1)) return gcd(a, b >> 1); // 4
	if (a > b) return gcd((a - b) >> 1, b); // 5
	else return gcd((b - a) >> 1, a); // 5
}

可能这个代码看上去比欧几里得算法的时间复杂度高，但其实它比欧几里得算法还快上不少。

接下来，我们进行一个优化，将递归转化为递推，进一步减少时间复杂度（用到一些二进制的知识）。

我们先把两个数中的 \(2\) 全部取出，也就是把两个数的二进制末尾的 \(0\) 全部取出，保证其中至少有一个奇数，然后两个数再分别消掉末尾的 \(0\)，保证两个数都是奇数，最终使用更相减损术。

int gcd(int a, int b) {
	if (!a) return b; // a = 0;
	if (!b) return a; // b = 0;
	int shift;
	for (shift = 0; ~(a | b) & 1; shift++) 
		a >>= 1, b >>= 1;
	/*
	先 a | b，方便 2 个数一起消 0
	取反再与 1，判断末尾是不是 0 
	*/
	while (~a & 1) a >>= 1;
	// 判断 a 的末尾是不是 0
	do {
		while (!(b & 1)) b >>= 1; // 消 b 末尾的 0
		// 现在都是奇数 
		if (a > b) swap(a, b);
		b -= a; // 更相减损术
	} while (b); 
	return a << shift; // 公因数 2 ^ shift 
}

但是，这还是不够快！每次都用 while，效率低下，于是，我在这里隆重介绍我们的硬件函数__builtin，以这个名字开头的函数名都是运行速度极快！这里用到的是__builtin_ctz()，它用来处理一个数二进制末尾 \(0\) 的个数。于是，我们的代码就变成了这样：

int gcd(int a, int b) {
	int az = __builtin_ctz(a), bz = __builtin_ctz(b);
	int z = az > bz ? bz : az, diff; // z 是后面要恢复的 2^z 
	b >>= bz;
	while (a) {
		a >>= az; // 右移 
		diff = b - a; // 更相减损术 
		az = __built_ctz(diff);
		if (a < b) b = a;  // 更新 b 
		a = diff < 0 ? -diff : diff; // 更新 a, diff < 0取相反数 
	}
	return b << z; // 答案 
}

到这里就结束啦，你学废了吗？

参考链接：这里 和 这里