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

Posted 2023-02-22

【中文标题】如何使用 __builtin_ctz 加速二进制 GCD 算法？

【英文标题】：How can I speed up the binary GCD algorithm using __builtin_ctz?

【发布时间】：2020-08-26 19:57:25

【问题描述】：

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.

【问题讨论】：

二进制 GCD 的另一个样本基本上有 min，而你有 max，但它的整体工作方式有点不同

与std::gcd相比如何？这应该更快吗？

@TedLyngmo 不确定。 I tried to benchmark it 但我的基准测试出现了段错误。请注意，该链接指向 Godbolt，因为 quick-bench 不会让您分享失败的基准测试。你知道这个基准有什么问题吗？无论如何，它并不真正应该更快，问题在于了解如何使用 CTZ。

这是你gcd 中的一个***。它进入了一个非常深的递归@u=3508125240，v=2952784951。这里是max、min和shiftvalues。

@TedLyngmo 感谢您的澄清。我的实现比std::gcd 快，BrettHale 的甚至比这更快。 （有关基准结果，请参阅我的答案）

【参考方案1】：

这是我在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)，@ 987654330@ 是奇数，因此|u - v| 必须是偶数。

这不是严格必要的，因为“修改”步骤：u >>= __builtin_clz(u); 将在下一次迭代中有效地执行此操作。

假设(u) 或(v) 具有“随机”位分布，则(n) 通过tzcnt 出现尾随零的概率为~(1/(2^n))。此指令是对 bsf 的改进，后者是 Haswell, IIRC 之前的 __builtin_clz 的实现。

【讨论】：

就像这个算法的 ARM 输出一样，它可以使用 rev/clz 来处理 ctz。

quick-bench: clang, gcc - 看起来真不错。

64bit version 比 std::gcd 的优势更大。尽管在 Intel Ice Lake 上结果可能有所不同，但它提高了除法速度。

tzcnt 有一个错误的依赖关系——实际上列为勘误表——经常插入xor r, r 来破坏它们。我不确定最近的微架构是否解决了这个问题。 bsf 得到了很多的改进，但是很多编译器都在使用它。

@harold - 那么我们就可以使用：return ((v == 0) ? u : u64_gcd(v, u % v)); ...但我已经读到 Icelake 的部门得到了很大的推动。

【参考方案2】：

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

【讨论】：

顺便说一句，您可以使用__builtin_ctz(u | v)获取最小尾随零计数

使用 max 而不是 min 是一个错误。但是可以使用单独的班次。因为您从其中一个中删除的 2 的额外因数不会改变最大的 COMMON 除数。

如果我取 2 个奇数并乘以 2 的幂，我没有改变 gcd。

只有当您不包括递归调用开销的成本时，它才是“无分支”的。作为尾递归算法，转换成iterative implementation比较容易。

更喜欢@BrettHale 的回答。编译器更擅长改进这一点，递归很少是最快的答案。