二进制gcd（Binary gcd）

二进制gcd（Binary gcd）

前置知识：$\gcd$（最大公约数），位运算

介绍

众所周知利用欧几里得算法（辗转相除法）可以做到在 $\log$ 的时间复杂度内求出两个数的最大公约数

但是有时候这还是不够快，例如这道题——Luogu P5435

如果直接辗转相除硬求 $\gcd$ 就会 T 得飞起，但是题解里那种方法 $\Theta(值域)$ 预处理 $\Theta(1)$ 查询又好麻烦，那么有没有一种更快的直接求两个数 $\gcd$ 的方法呢？

二进制 $\gcd$ 就是这样一种算法，它利用计算机优秀的二进制处理能力来计算 $\gcd$

原理

假设现在我们需要求 $\gcd(a,b)$ 满足 $a,b>0$

我们将 $a,b$ 分为 $5$ 种情况讨论

第一种，$a=b$。显然 $\gcd(a,b)=a=b$

第二种，$a=0$ 或 $b=0$。根据最大公约数的定义当 $b=0$ 时 $\gcd(a,b)=\gcd(a,0)=a$，当 $a=0$ 时 $\gcd(a,b)=\gcd(b,0)=b$

第三种，$a,b$ 均为偶数。对于这种情况，显然 $2$ 是公约数之一，直接将两个数都除以 $2$，再递归下去求 $\gcd(\frac a2,\frac b2)$ 即可，$\gcd(a,b)=2\cdot\gcd(\frac a2,\frac b2)$

第四种，$a,b$ 中有且仅有一个偶数。不妨设 $a$ 是偶数，那么显然 $2$ 不是公约数之一，直接将 $a$ 除以 $2$ 求 $\gcd(\frac a2,b)$ 即可，$\gcd(a,b)=\gcd(\frac a2,b)$

第五种，$a,b$ 均为奇数。不妨设 $a>b$，那么有 $\gcd(a,b)=\gcd(\frac{a-b}2,b)$，为什么呢？下面来证明一下

设 $a=k_1\gcd(a,b),\ b=k_2\gcd(a,b)$ 且 $\gcd(k_1,k_2)=1$，先考虑证明 $\gcd(a-b,b)=\gcd(a,b)$（即“更相减损术”）

$a-b=(k_1-k_2)\gcd(a,b),\ b=k_2\gcd(a,b)$，那么有 $\gcd(a-b,b)=\gcd(k_1-k_2,k_2)\cdot\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,\ k_2=t_2m$

就有 $k_1=(t_1+t_2)m$，那么 $\gcd(k_1,k_2)=m\cdot\gcd(t_1+t_2,t_2)>1$ 与 $k_1,k_2$ 互质矛盾，得证 $\gcd(k_1-k_2,k_2)=1$

综上 $\gcd(a-b,b)=\gcd(a,b)$

由于 $a,b$ 均为奇数，$a-b$ 就是偶数，那么问题就变成了和情况四一样，$2$ 显然不是公约数之一，故 $\gcd(a,b)=\gcd(\frac{a-b}2,b)$

实现

直接根据上图模拟，我们可以得到以下代码（注意：~ 的运算优先级高于 | 和 &）

int gcd(int a,int b)
{
	if(a==b) return a;//a=b
	if(!a) return b;//a=0
	if(!b) return a;//b=0
	if(~a&1)//a是偶数
	{
		if(b&1)//b是奇数
			return gcd(a>>1,b);
		else//b是偶数
			return gcd(a>>1,b>>1)<<1;
	}
	if(~b&1)//a是奇数b是偶数
		return gcd(a,b>>1);
	//均为奇数
	if(a>b) return gcd((a-b)>>1,b);
	return gcd((b-a)>>1,a);
}

这样使用了大量的位运算，虽然看起来比辗转相除法递归层数更多，但是实际上大大优化了运算速度

当然我们可以把递归展开写成循环的形式，效率还能更高

首先把两个数公约数中的所有 $2$ 提出来，也就是把二进制末尾的 $0$ 全部提出来，然后两个数就至少有一个是奇数，把两个数末尾的 $0$ 再分别去掉就能保证一定是两个奇数，利用更相减损术循环求解就行，具体看代码

int gcd(int a,int b)
{
	if(!a)return b;//a=0
	if(!b)return a;//b=0
	int shift;//a,b二进制下末尾0的个数
	for(shift=0;~(a|b)&1;++shift)a>>=1,b>>=1;//把两个数末尾的shift个0都去掉
	while(~a&1)a>>=1;//现在两个数至少有一个奇数，2一定不是公约数
	//现在a一定是奇数
	do
	{
		while(!(b&1))b>>=1;
		//现在两个数一定都是奇数
		if(a>b)swap(a,b);
		b-=a;//更相减损
	}while(b);
	return a<<shift;//把最开始去掉的2^shift的公约数加上
}

注意到每次还是都使用 while 来去除二进制下末尾的 $0$，效率低下，可以利用硬件函数 __builtin_ctz() 来继续加速

同时只有直接输入 $0$ 时才可能 $a=0$ 或是 $b=0$，所以如果输入保证不为 $0$ 就可以去掉这个判断

再稍加卡常就得到了终极代码

inline int gcd(int a, int b)
{
	register int az=__builtin_ctz(a),bz=__builtin_ctz(b),z=az>bz?bz:az,diff;
	b>>=bz;
	while(a)
	{
		a>>=az;
		diff=b-a;
		az=__builtin_ctz(diff);
		if(a<b)b=a;
		a=diff<0?-diff:diff;
	}
	return b<<z;
}

这样我们就成功地利用 $\Theta(\log n)$ 求 $\gcd$ 切掉了 Luogu P5435 这道需要 $\Theta(\text{值域})$ 预处理 $\Theta(1)$ 求 $\gcd$ 的题目

---

博客园传送门

知乎传送门