如何仅用便携64位运算求(2^128-1)/x的商的低64位？

我需要计算(2^128 - 1) / x的值。除数x是无符号64位整数，被除数由两个无符号64位整数（高64位、低64位）组成，两个值均为UINT64_MAX。要求仅使用64位算术运算实现，且具备可移植性，不可使用GNU的__int128、MSVC的_udiv128、汇编或其他类似非标准扩展。我不需要商的高位部分，仅需获取商的低64位。
补充条件：x >= 3，且x不是2的幂。
补充说明：我已经实现了自己的解决方案，也欢迎各位提供性能更优的其他实现思路。

可移植基础实现 #

该方案完全基于标准C的64位无符号运算，无任何编译器扩展，可在所有支持C99及以上标准的平台运行：

#include <stdint.h>

// 内部函数：计算 (r << 64 | lo) / x，要求 r < x，返回商的低64位
static uint64_t div128_low_part(uint64_t r, uint64_t lo, uint64_t x) {
    const uint64_t HALF_BITS = 32;
    const uint64_t HALF_MASK = (1ULL << HALF_BITS) - 1;

    uint64_t x_high = x >> HALF_BITS;
    uint64_t x_low = x & HALF_MASK;

    // 计算近似商
    uint64_t approx_q;
    uint64_t r_shift32 = r << HALF_BITS;
    if ((r >> HALF_BITS) == x_high) {
        approx_q = UINT64_MAX;
    } else {
        approx_q = r_shift32 / x_high;
    }

    // 修正近似商
    uint64_t q_mul_xlow = approx_q * x_low;
    uint64_t temp_rem = r_shift32 + (lo >> HALF_BITS) - approx_q * x_high * (1ULL << HALF_BITS) - (q_mul_xlow >> HALF_BITS);
    uint64_t lo_low = lo & HALF_MASK;

    if (q_mul_xlow > (temp_rem << HALF_BITS) + lo_low) {
        approx_q--;
        temp_rem += x;
        if (temp_rem >= x && q_mul_xlow > (temp_rem << HALF_BITS) + lo_low) {
            approx_q--;
        }
    }

    return approx_q;
}

// 对外接口：计算(2^128 - 1)/x的低64位商
uint64_t compute_2pow128_minus1_div_x(uint64_t x) {
    const uint64_t UINT64_MAX_VAL = 0xFFFFFFFFFFFFFFFFULL;
    // 第一步：计算高64位除以x的余数
    uint64_t rem_high = UINT64_MAX_VAL % x;
    // 第二步：计算拼接后的值除以x的商，即为最终需要的低64位结果
    return div128_low_part(rem_high, UINT64_MAX_VAL, x);
}

原理解释 #

被除数2^128 - 1可以拆分为UINT64_MAX * 2^64 + UINT64_MAX，根据整数除法的拆分规则：
(a * 2^64 + b) / x = (a / x) * 2^64 + ((a % x) * 2^64 + b) / x
我们不需要商的高64位，所以只需要计算((a % x) * 2^64 + b) / x即可，其中a和b都是UINT64_MAX。
内部的div128_low_part函数模拟了手工长除法的逻辑，将64位除数拆分为高低32位，先通过高32位得到近似商，再经过最多两次修正得到精确结果，全程没有使用128位运算。

高性能优化思路（针对x为奇数的场景） #

因为x不是2的幂，所有x都包含奇数因子，如果x本身是奇数，可以用乘法逆元加速计算：

1. 先计算x模2^64的乘法逆元inv_x，满足x * inv_x ≡ 1 mod 2^64
2. 计算R = (2^128 - 1) % x，可以通过快速幂模运算快速得到
3. 最终结果等于( (UINT64_MAX * inv_x * 2) - R * inv_x ) & UINT64_MAX
   该方法运算量远小于长除法，适合需要多次计算不同x值的场景。

内容的提问来源于stack exchange，提问作者tgonzalez89