二进制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\) 的题目
