如何基于uint64_t数组实现uint1024_t的五大基础运算

我们基于以下结构体实现1024位无符号整数：

#include <stdint.h>
#include <string.h>

typedef struct {
    uint64_t chunk[16];
} uint1024_t;

1. 加法运算 #

按小端存储（chunk[0]为最低64位），从最低位chunk开始逐位相加，传递进位：

void uint1024_add(const uint1024_t* a, const uint1024_t* b, uint1024_t* res) {
    uint64_t carry = 0;
    for (int i = 0; i < 16; i++) {
        __uint128_t sum = (__uint128_t)a->chunk[i] + b->chunk[i] + carry;
        res->chunk[i] = (uint64_t)sum;
        carry = (uint64_t)(sum >> 64);
    }
    // 无符号加法溢出时，高位进位自动丢弃，符合无符号整数溢出规则
}

2. 减法运算 #

处理借位逻辑，若被减数小于减数，结果为模2^1024的差值（无符号溢出结果）：

void uint1024_sub(const uint1024_t* a, const uint1024_t* b, uint1024_t* res) {
    uint64_t borrow = 0;
    for (int i = 0; i < 16; i++) {
        uint64_t a_chunk = a->chunk[i];
        uint64_t b_chunk = b->chunk[i] + borrow;
        
        borrow = (a_chunk < b_chunk) ? 1 : 0;
        res->chunk[i] = a_chunk - b_chunk + (borrow ? (1ULL << 64) : 0);
    }
}

3. 乘法运算 #

模仿竖式乘法，用每个chunk两两相乘（结果为128位），累加到对应位置后处理进位：

void uint1024_mul(const uint1024_t* a, const uint1024_t* b, uint1024_t* res) {
    memset(res, 0, sizeof(uint1024_t));
    
    __uint128_t temp[31] = {0}; // 中间结果最多需要31个64位存储单元
    
    for (int i = 0; i < 16; i++) {
        if (a->chunk[i] == 0) continue;
        for (int j = 0; j < 16; j++) {
            temp[i + j] += (__uint128_t)a->chunk[i] * b->chunk[j];
        }
    }
    
    uint64_t carry = 0;
    for (int i = 0; i < 16; i++) {
        __uint128_t total = temp[i] + carry;
        res->chunk[i] = (uint64_t)total;
        carry = (uint64_t)(total >> 64);
    }
}

4. 辅助比较函数 #

用于除法和模运算中判断数值大小：

// 返回值：1(a > b), 0(a == b), -1(a < b)
int uint1024_cmp(const uint1024_t* a, const uint1024_t* b) {
    for (int i = 15; i >= 0; i--) {
        if (a->chunk[i] > b->chunk[i]) return 1;
        if (a->chunk[i] < b->chunk[i]) return -1;
    }
    return 0;
}

5. 辅助移位函数 #

用于除法中的移位操作：

void uint1024_shift_left(uint1024_t* num, int bits) {
    if (bits == 0) return;
    
    if (bits >= 64) {
        int shift_blocks = bits / 64;
        int shift_bits = bits % 64;
        
        // 整体移动chunk
        for (int i = 15; i >= shift_blocks; i--) {
            num->chunk[i] = num->chunk[i - shift_blocks];
        }
        memset(num->chunk, 0, shift_blocks * sizeof(uint64_t));
        
        if (shift_bits > 0) {
            uint64_t carry = 0;
            for (int i = 0; i < 16; i++) {
                uint64_t new_carry = num->chunk[i] >> (64 - shift_bits);
                num->chunk[i] = (num->chunk[i] << shift_bits) | carry;
                carry = new_carry;
            }
        }
        return;
    }
    
    // 小于64位的移位
    uint64_t carry = 0;
    for (int i = 0; i < 16; i++) {
        uint64_t new_carry = num->chunk[i] >> (64 - bits);
        num->chunk[i] = (num->chunk[i] << bits) | carry;
        carry = new_carry;
    }
}

6. 除法与模运算 #

采用移位减法试商法，从高位到低位确定商的每一位，最终得到商和余数：

void uint1024_div_mod(const uint1024_t* dividend, const uint1024_t* divisor, uint1024_t* quotient, uint1024_t* remainder) {
    // 除数为0的情况可根据需求添加断言或报错逻辑
    if (uint1024_cmp(divisor, &(uint1024_t){0}) == 0) return;
    
    memset(quotient, 0, sizeof(uint1024_t));
    memcpy(remainder, dividend, sizeof(uint1024_t));
    
    uint1024_t shifted_divisor;
    memcpy(&shifted_divisor, divisor, sizeof(uint1024_t));
    
    // 找到除数左移后不超过余数的最大移位次数
    int shift = 0;
    while (uint1024_cmp(&shifted_divisor, remainder) <= 0 && shift < 1024) {
        uint1024_shift_left(&shifted_divisor, 1);
        shift++;
    }
    
    // 从高位到低位试商
    while (shift > 0) {
        shift--;
        uint1024_shift_left(&shifted_divisor, -1); // 右移1位
        
        if (uint1024_cmp(remainder, &shifted_divisor) >= 0) {
            uint1024_sub(remainder, &shifted_divisor, remainder);
            // 商的对应位置1
            uint1024_t bit_set = {0};
            bit_set.chunk[shift / 64] = 1ULL << (shift % 64);
            uint1024_add(quotient, &bit_set, quotient);
        }
    }
}

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