问题描述
clang和GCC具有int __builtin_ctz(unsigned)
函数。这将以整数形式计算尾随零。 Wikipedia article on this family of functions提到可以使用__builtin_ctz
加速二进制GCD算法,但我不知道如何。
二进制GCD的sample implementation看起来像这样:
unsigned int gcd(unsigned int u,unsigned int v)
{
// simple cases (termination)
if (u == v)
return u;
if (u == 0)
return v;
if (v == 0)
return u;
// look for factors of 2
if (~u & 1) // u is even
if (v & 1) // v is odd
return gcd(u >> 1,v);
else // both u and v are even
return gcd(u >> 1,v >> 1) << 1;
if (~v & 1) // u is odd,v is even
return gcd(u,v >> 1);
// reduce larger argument
if (u > v)
return gcd(u - v,v);
return gcd(v - u,u);
}
我怀疑我可以如下使用__builtin_ctz
:
constexpr unsigned int gcd(unsigned int u,unsigned int v)
{
// simplified first three ifs
if (u == v || u == 0 || v == 0)
return u | v;
unsigned ushift = __builtin_ctz(u);
u >>= ushift;
unsigned vshift = __builtin_ctz(v);
v >>= vshift;
// Note sure if max is the right approach here.
// In the if-else block you can see both arguments being rshifted
// and the result being leftshifted only once.
// I expected to recreate this behavior using max.
unsigned maxshift = std::max(ushift,vshift);
// The only case which was not handled in the if-else block before was
// the odd/odd case.
// We can detect this case using the maximum shift.
if (maxshift != 0) {
return gcd(u,v) << maxshift;
}
return (u > v) ? gcd(u - v,v) : gcd(v - u,u);
}
int main() {
constexpr unsigned result = gcd(5,3);
return result;
}
不幸的是,这还行不通。程序结果为4,应为1。那么,我在做什么错呢?如何在这里正确使用__builtin_ctz
? See my code so far on GodBolt。
解决方法
这是我从comments开始的迭代实现:
尽管尾部递归算法通常很优雅,但是迭代实现在实践中几乎总是更快。 (现代编译器实际上可以在非常简单的情况下执行此转换。)
unsigned ugcd (unsigned u,unsigned v)
{
unsigned t = u | v;
if (u == 0 || v == 0)
return t; /* return (v) or (u),resp. */
int g = __builtin_ctz(t);
while (u != 0)
{
u >>= __builtin_ctz(u);
v >>= __builtin_ctz(v);
if (u >= v)
u = (u - v) / 2;
else
v = (v - u) / 2;
}
return (v << g); /* scale by common factor. */
}
如上所述,|u - v| / 2
步骤通常实现为非常有效的无条件右移,例如 shr r32
,以除以(2)
- (u)
和(v)
都是奇数,因此|u - v|
必须是偶数。
这不是严格必要的,因为“整理”步骤:u >>= __builtin_clz(u);
将在下一次迭代中有效地执行此操作。
假设(u)
或(v)
具有“随机”比特分布,则(n)
通过 tzcnt
尾随零的概率为〜(1/(2^n))
。该说明是对 bsf
(IIRC的Haswell之前的__builtin_clz
的实现)的改进。
感谢有帮助的评论员,我发现了一个关键错误:我应该使用min
而不是max
这是最终解决方案:
#include <algorithm>
constexpr unsigned gcd(unsigned u,unsigned v)
{
if (u == v || u == 0 || v == 0)
return u | v;
// effectively compute min(ctz(u),ctz(v))
unsigned shift = __builtin_ctz(u | v);
u >>= __builtin_ctz(u);
v >>= __builtin_ctz(v);
const auto &[min,max] = std::minmax(u,v);
return gcd(max - min,min) << shift;
}
int main() {
constexpr unsigned g = gcd(25,15); // g = 5
return g;
}
此解决方案也很不错,nearly branch-free compile output。
以下是到目前为止所有答案中的some benchmark results(我们实际上击败了std::gcd
):