如何使用 __builtin_ctz 加速二进制 GCD 算法？

Question

clang 和 GCC 有一个 int __builtin_ctz(unsigned) 函数。这将计算整数中的尾随零。 Wikipedia article on this family of functions 提到二进制 GCD 算法可以使用 __builtin_ctz 加速，但我不明白如何。

二进制 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.

Answer 1

感谢有帮助的评论员，我发现了关键错误：我应该使用 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）：

Answer 2

这是我在 comments 中的迭代实现：

虽然 tail-recursive 算法通常很优雅，但迭代实现在实践中几乎总是更快。（现代编译器实际上可以在非常简单的情况下执行此转换。）

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| 必须是偶数。

严格来说 没有必要，因为 'oddifying' 步骤：u >>= __builtin_clz(u); 将在下一次迭代中有效地执行此操作。

假设 (u) 或 (v) 具有 'random' 位分布，(n) 尾随零的概率，通过 tzcnt ，是~(1/(2^n))。该指令是对 bsf 的改进，__builtin_clz 在 Haswell，IIRC 之前的实现。

如何使用 __builtin_ctz 加速二进制 GCD 算法？

How can I speed up the binary GCD algorithm using __builtin_ctz?

c++

algorithm

bit-manipulation

built-in

greatest-common-divisor