Beating C++ std::gcd

["Ruud H.G. van Tol" <rvtol@isolution.nl>]

In a project I'm working on, I measured that about 25% of its runtime 
was spent in gcd.
The program runs with 10 threads, on a MacBook M1 Max.

On my clang, std::gcd came out as the best suited implementation.
Then I started looking for more speed, and tried a lookup table for 
smaller values.
And that made it 20% faster!

  gcd(a, b); 0 <= a <= 511; 0 <= b <= 1023.

Implementation attached.

-- Greetings, Ruud

using i128 = __int128_t;
using u128 = __uint128_t;

// Small-value GCD table: gcd(i,j) for i in [0,GCD_TABLE_NA), j in
// [0,GCD_TABLE_NB). Checked inside igcd's loop below; once both values drop
// under their threshold, the remaining reduction is a single lookup instead
// of further iterations. Asymmetric on purpose.
// NA=512,NB=1024 measured within 0.2 percentage points of a
// symmetric 1024x1024 table's speedup, at half the memory. Smaller on
// either axis loses more (down to 256/256: -69%, vs -92% here) -- range
// dominates over per-lookup cache cost, but only on the axis whose real
// distribution actually needs it. Populated once at program startup (not
// compile-time constexpr): building it at compile time needs a raised,
// toolchain-specific step limit (GCC's -fconstexpr-ops-limit, Clang's
// -fconstexpr-steps, neither portable to a fixed value across
// compilers/versions) -- not acceptable for a file meant to build with a
// plain g++ invocation anywhere.
constexpr int GCD_TABLE_NA = 512, GCD_TABLE_NB = 1024;
std::array<std::array<uint16_t, GCD_TABLE_NB>, GCD_TABLE_NA> make_gcd_table() {
    std::array<std::array<uint16_t, GCD_TABLE_NB>, GCD_TABLE_NA> t{};
    for (int i = 0; i < GCD_TABLE_NA; ++i)
        for (int j = 0; j < GCD_TABLE_NB; ++j)
            t[i][j] = (uint16_t)std::gcd(i, j);
    return t;
}
const auto GCD_TABLE = make_gcd_table();

constexpr int ctz128(const u128 x) { // trailing-zero count, x!=0
    uint64_t lo = (uint64_t)x;
    if (lo != 0) return __builtin_ctzll(lo);
    return 64 + __builtin_ctzll((uint64_t)(x >> 64));
}
// Hybrid Euclidean/binary GCD, u128 internally. One Euclidean mod reduces a
// large size mismatch between numerator and denominator in a single step,
// before switching to Stein's binary algorithm for the rest -- matches
// libc++'s std::gcd (itself based on lemire.me/blog/2024/04/13/), confirmed
// faster than plain Stein's alone on exactly this operand shape.
constexpr i128 igcd(i128 a, i128 b) {
    u128 ua = (u128)iabs128(a), ub = (u128)iabs128(b);
    if (ua == 0) return (i128)ub;
    if (ub == 0) return (i128)ua;

    if (ua > ub) { u128 t = ua; ua = ub; ub = t; }  // now ua < ub
    ub %= ua;
    if (ub == 0) return (i128)ua;

    int shift = ctz128(ua | ub);
    ua >>= ctz128(ua);
    bool use_table = !std::is_constant_evaluated();
    do {
        u128 t = ub >> ctz128(ub);
        if (use_table && ua < GCD_TABLE_NA && t < GCD_TABLE_NB)
            return (i128)((u128)GCD_TABLE[(size_t)ua][(size_t)t] << shift);
        if (ua > t) { ub = ua - t; ua = t; }
        else        { ub = t - ua; }
    } while (ub != 0);
    return (i128)(ua << shift);
}