无更大容器时高精度计算128位整数(a*b)/c的方法

嘿，这个问题我之前在处理大整数运算的时候踩过坑，确实128位乘128位没法直接依赖原生256位类型，但有两种靠谱的方法能实现无精度损失的(a*b)/c计算，刚好你的场景里a、b≤c，还能简化一些步骤，给你详细拆解：

1. 拆分模拟256位乘法+长除法 #

既然没有原生256位类型，我们可以手动把128位整数拆成高64位和低64位，分步计算乘积，再组合成完整的256位结果，最后做256位除以128位的长除法。因为a、b≤c，所以a*b≤c²，而c是128位，c²刚好是256位，整个过程不会有溢出问题。

具体步骤： #

• 拆分128位整数：把每个128位数拆成高64位和低64位，比如 a = a_high << 64 | a_low，b 同理拆成 b_high 和 b_low。
• 计算四部分乘积：分别计算 a_low*b_low、a_low*b_high、a_high*b_low、a_high*b_high，这四个都是64位乘64位，结果是128位。
• 组合成256位乘积：把四个部分按位累加，得到256位的乘积（用两个128位变量分别存高128位和低128位）。
• 256位除以128位：用长除法逐步计算商，因为a、b≤c，最终结果不会超过c，计算过程可以简化。

伪代码实现（C风格）： #

typedef unsigned __int128 uint128_t;

// 辅助函数：拆分128位整数为高低64位
void split_uint128(uint128_t x, uint64_t* high, uint64_t* low) {
    *low = (uint64_t)x;
    *high = (uint64_t)(x >> 64);
}

// 核心计算函数：(a*b)/c，a、b、c均为128位无符号整数，且a,b<=c
uint128_t mul_div_u128(uint128_t a, uint128_t b, uint128_t c) {
    uint64_t a_high, a_low, b_high, b_low;
    split_uint128(a, &a_high, &a_low);
    split_uint128(b, &b_high, &b_low);
    
    // 计算四个64x64的128位乘积
    uint128_t p0 = (uint128_t)a_low * b_low;
    uint128_t p1 = (uint128_t)a_low * b_high;
    uint128_t p2 = (uint128_t)a_high * b_low;
    uint128_t p3 = (uint128_t)a_high * b_high;
    
    // 组合成256位乘积：high128 << 128 | low128
    uint128_t low128 = p0;
    uint128_t carry = (uint128_t)(uint64_t)p1 + (uint64_t)(p0 >> 64);
    low128 = (carry << 64) | (uint64_t)low128;
    carry = (p1 >> 64) + (uint64_t)p2 + (carry >> 64);
    uint128_t high128 = p3 + (p2 >> 64) + carry;
    
    // 256位除以128位的长除法
    uint128_t result = 0;
    uint128_t div_high = high128;
    uint128_t div_low = low128;
    
    // 先处理高位部分，如果高位大于等于除数，先做初步除法
    if (div_high >= c) {
        result += div_high / c;
        div_high = div_high % c;
    }
    
    // 二进制长除法：从最高位到最低位试商
    for (int i = 127; i >= 0; i--) {
        // 被除数左移一位
        div_high = (div_high << 1) | (div_low >> 127);
        div_low <<= 1;
        
        // 如果当前高位部分大于等于除数，减去除数并标记商位
        if (div_high >= c) {
            div_high -= c;
            result |= ((uint128_t)1 << i);
        }
    }
    
    return result;
}

2. 牛顿迭代法求逆元（高性能场景） #

如果需要频繁计算这类除法，牛顿迭代法求乘法逆元的速度会更快。核心思路是先求c的乘法逆元inv(c)，然后用(a*b)*inv(c)取高位得到商，不过需要先处理c为偶数的情况。

具体步骤： #

• 提取2的因子：把c拆成c = 2^s * d，其中d是奇数（s是c末尾的0的个数）。
• 调整a和b：把a和b也分别除以2的幂次，使得总除以的2的幂次等于s（比如a' = a >> s_a，b' = b >> s_b，s_a + s_b = s）。
• 求奇数d的逆元：用牛顿迭代法计算d的模2^128逆元inv_d，满足d * inv_d ≡ 1 mod 2^128。
• 计算结果：(a'*b') * inv_d的高128位就是(a'*b')/d，也就是最终的(a*b)/c。

逆元计算的牛顿迭代伪代码： #

// 计算奇数d的模2^128逆元
uint128_t newton_inv(uint128_t d) {
    uint128_t inv = 3; // 初始值，3是2^2的逆元，适合奇数
    // 迭代4次足够收敛到128位精度
    for (int i = 0; i < 4; i++) {
        inv = inv * 2 - d * inv * inv;
    }
    return inv;
}

优缺点对比： #

• 拆分法：逻辑直观，容易实现，不需要额外处理特殊情况，适合一次性计算。
• 牛顿迭代法：计算速度快（迭代几次就能得到结果），适合大量重复计算的场景，但需要处理偶数除数的情况。

验证小例子 #

比如取a=3（128位），b=4（128位），c=5（128位），(3*4)/5=2，用上面的两种方法计算都会得到正确的结果，不会有精度损失。

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