384位无符号整数Karatsuba乘法C实现错误排查及优化求助

问题描述 #

以下是实现384位无符号整数Karatsuba乘法的C代码：

typedef struct {
    uint64_t chunk[6];
} uint384_t;

void add_uint384(const uint384_t *a, const uint384_t *b, uint384_t *res) {
    uint64_t carry = 0;
    for (int i = 0; i < 6; i++) {
        uint64_t temp = a->chunk[i] + b->chunk[i] + carry;
        carry = (temp < a->chunk[i]) ? 1 : 0;
        res->chunk[i] = temp;
    }
}

void subtract_uint384(const uint384_t *a, const uint384_t *b, uint384_t *res) {
    uint64_t borrow = 0;
    for (int i = 0; i < 6; i++) {
        uint64_t temp = a->chunk[i] - b->chunk[i] - borrow;
        borrow = (a->chunk[i] < (b->chunk[i] + borrow)) ? 1 : 0;
        res->chunk[i] = temp;
    }
}

void multiply_uint192(const uint64_t *a, const uint64_t *b, uint64_t *res) {
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            __uint128_t temp = (__uint128_t)a[i] * b[j];
            uint64_t low = (uint64_t)temp;
            uint64_t high = (uint64_t)(temp >> 64);
            res[i + j] += low;
            if (res[i + j] < low) high++;
            res[i + j + 1] += high;
        }
    }
}

void multiply_two_variables_karatsuba(const uint384_t *a, const uint384_t *b, uint384_t *c) {
    uint64_t a_low[3], a_high[3], b_low[3], b_high[3];
    uint64_t z0[6] = {0}, z1[6] = {0}, z2[6] = {0};
    uint64_t a_sum[3] = {0}, b_sum[3] = {0};

    for (int i = 0; i < 3; i++) {
        a_low[i] = a->chunk[i];
        a_high[i] = a->chunk[i + 3];
        b_low[i] = b->chunk[i];
        b_high[i] = b->chunk[i + 3];
        a_sum[i] = a_low[i] + a_high[i];
        b_sum[i] = b_low[i] + b_high[i];
    }

    multiply_uint192(a_low, b_low, z0);
    multiply_uint192(a_high, b_high, z2);
    multiply_uint192(a_sum, b_sum, z1);

    for (int i = 0; i < 6; i++) {
        z1[i] -= z0[i];
        z1[i] -= z2[i];
    }

    for (int i = 0; i < 3; i++) {
        c->chunk[i] = z0[i];
        c->chunk[i + 3] = z0[i + 3] + z1[i];
        c->chunk[i + 6] = z2[i];
    }
}

测试时，将a和b初始化为：

0x00000000FFFFFFFF_00000000FFFFFFFF_00000000FFFFFFFF_00000000FFFFFFFF_00000000FFFFFFFF_00000000FFFFFFFF

实际输出结果：

0xfffffff40000000b_fffffff600000009_fffffff800000006_fffffffa00000004_fffffffc00000002_fffffffe00000001

预期结果：

0xfffffff40000000a_fffffff600000008_fffffff800000006_fffffffa00000004_fffffffc00000002_fffffffe00000001

仅前两个chunk值不符合预期（错误值比预期大1），需排查问题并做性能优化。

错误原因分析 #

1. 数组越界访问：uint384_t仅包含6个uint64_t chunk（索引0-5），但代码中c->chunk[i + 6] = z2[i]直接访问索引6-8，属于越界写操作，会破坏内存并导致结果异常。
2. 减法未处理借位：计算z1 = z1 - z0 - z2时，直接逐位执行无符号减法，未处理借位。当z1[i]小于z0[i]+z2[i]时，无符号溢出会导致结果错误。
3. 求和时进位丢失：a_sum[i] = a_low[i] + a_high[i]未处理uint64_t溢出后的进位，导致a_sum和b_sum的值不准确，进而影响z1的计算。
4. 结果合并未处理进位：z0[i + 3] + z1[i]可能产生溢出，但代码未处理该进位，会导致高位chunk计算错误。

修复方案 #

修复后的核心乘法函数 #

void multiply_two_variables_karatsuba(const uint384_t *a, const uint384_t *b, uint384_t *c) {
    uint64_t a_low[3], a_high[3], b_low[3], b_high[3];
    uint64_t z0[6] = {0}, z1[6] = {0}, z2[6] = {0};
    uint64_t a_sum[3] = {0}, b_sum[3] = {0};

    for (int i = 0; i < 3; i++) {
        a_low[i] = a->chunk[i];
        a_high[i] = a->chunk[i + 3];
        b_low[i] = b->chunk[i];
        b_high[i] = b->chunk[i + 3];
        
        // 处理求和进位，避免溢出丢失
        __uint128_t sum = (__uint128_t)a_low[i] + a_high[i];
        a_sum[i] = (uint64_t)sum;
        if (i < 2) a_sum[i+1] += (uint64_t)(sum >> 64);
        
        sum = (__uint128_t)b_low[i] + b_high[i];
        b_sum[i] = (uint64_t)sum;
        if (i < 2) b_sum[i+1] += (uint64_t)(sum >> 64);
    }

    multiply_uint192(a_low, b_low, z0);
    multiply_uint192(a_high, b_high, z2);
    multiply_uint192(a_sum, b_sum, z1);

    // 带借位执行z1 = z1 - z0 - z2
    uint64_t borrow = 0;
    for (int i = 0; i < 6; i++) {
        uint64_t val = z1[i] - z0[i] - borrow;
        borrow = (z1[i] < (z0[i] + borrow)) ? 1 : 0;
        z1[i] = val;
    }
    borrow = 0;
    for (int i = 0; i < 6; i++) {
        uint64_t val = z1[i] - z2[i] - borrow;
        borrow = (z1[i] < (z2[i] + borrow)) ? 1 : 0;
        z1[i] = val;
    }

    // 合并结果并处理进位（仅保留前384位）
    uint64_t carry = 0;
    // 低192位直接取自z0的低3个chunk
    for (int i = 0; i < 3; i++) {
        c->chunk[i] = z0[i];
    }
    // 高192位 = z0的高3个chunk + z1的低3个chunk + 进位
    for (int i = 0; i < 3; i++) {
        __uint128_t temp = (__uint128_t)z0[i + 3] + z1[i] + carry;
        c->chunk[i + 3] = (uint64_t)temp;
        carry = (uint64_t)(temp >> 64);
    }
    // 超出384位的进位直接忽略
}

性能优化建议 #

1. 内联小函数：给add_uint384、subtract_uint384、multiply_uint192添加static inline修饰，减少函数调用开销。
2. 用__uint128_t简化运算：在加法、乘法中直接使用__uint128_t类型计算，减少手动进位判断的代码复杂度与性能损耗。
3. 循环展开：对固定次数的循环（如3次、6次）手动展开，或开启编译器-O3优化自动展开循环。
4. 减少内存操作：避免不必要的数组复制，直接在uint384_t结构上操作，降低栈内存占用。
5. 编译器优化参数：编译时添加-O3 -march=native，让编译器针对目标平台生成最优机器码。
6. 优化分治策略：当前384位分拆为两个192位的策略已足够高效，若需支持更大整数可递归使用Karatsuba算法。

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