如何使用uint64_t数组进行1024位操作

我试图找到一种方法来计算类型为uint1024_t(无符号1024位整数)的值，通过定义5个基本操作:加，减，乘，除，模。

我可以这样做的方法是创建一个结构，它将具有以下原型:

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

现在，因为用uint64_t作为块大小来包装我的头是复杂的，我首先写了一些代码来操纵uint8_t。下面是我想到的:

#define UINT8_HI(x) (x >> 4)
#define UINT8_LO(x) (((1 << 4) - 1) & x)
void uint8_add(uint8_t a, uint8_t b, uint8_t *res, int i) {
uint8_t s0, s1, s2;
uint8_t x = UINT8_LO(a) + UINT8_LO(b);
s0 = UINT8_LO(x);
x = UINT8_HI(a) + UINT8_HI(b) + UINT8_HI(x);
s1 = UINT8_LO(x);
s2 = UINT8_HI(x);
uint8_t result = s0 + (s1 << 4);
uint8_t carry = s2;
res[1 + i] = result;
res[0 + i] = carry;
}
void uint8_multiply(uint8_t a, uint8_t b, uint8_t *res, int i) {
uint8_t s0, s1, s2, s3;
uint8_t x = UINT8_LO(a) * UINT8_LO(b);
s0 = UINT8_LO(x);
x = UINT8_HI(a) * UINT8_LO(b) + UINT8_HI(x);
s1 = UINT8_LO(x);
s2 = UINT8_HI(x);
x = s1 + UINT8_LO(a) * UINT8_HI(b);
s1 = UINT8_LO(x);
x = s2 + UINT8_HI(a) * UINT8_HI(b) + UINT8_HI(x);
s2 = UINT8_LO(x);
s3 = UINT8_HI(x);
uint8_t result = s1 << 4 | s0;
uint8_t carry = s3 << 4 | s2;
res[1 + i] = result;
res[0 + i] = carry;
}

它似乎工作得很好，但是我无法为除法，减法和模数定义相同的操作…

此外，我似乎无法看到如何实现相同的原则，我的自定义uint1024_t结构，即使它是几乎相同的几行代码来管理溢出。

如果你能帮助我实现我的结构的5个基本操作，我将非常感激。

编辑:

我在下面给出了解决这个问题的实现。

找到一种方法来计算…5种基本运算:加、减、乘、除、取模。

如果uint1024_t使用uint32_t，会更容易。

我建议1)最宽字体uintmax_t宽度的一半，或2)unsigned，以较小者为准。例如32位。

(也要考虑uintN_t以外的东西，以避免与c的未来版本冲突)

typedef struct {
uint32_t chunk[1024/32];
} u1024;

一些未经测试的代码示例，让OP了解如何使用uint32_t简化任务。

void u1024_mult(u1024 *product, const u1024 *a, const u1024 *b) {
memset(product, 0, sizeof product[0]);
unsigned n = sizeof product->chunk / sizeof product->chunk[0];
for (unsigned ai = 0; ai < n; ai++) {
uint64_t acc = 0;
uint32_t m = a->chunk[ai];
for (unsigned bi = 0; ai + bi < n; bi++) {
acc += (uint64_t) m * b->chunk[bi] + product->chunk[ai + bi];
product->chunk[ai + bi] = (uint32_t) acc;
acc >>= 32;
}
}
}

+,-与上述非常相似。

/,%可以合并成一个例程，同时计算商和余数。

在这里发布这些函数并不难，因为它确实与小学数学相同，但不是以10为底，而是以2为底³²．我反对张贴它，虽然这是一个有趣的锻炼自己。

我希望上面的*示例代码能够启发而不是回答。

uint8_t数组的实现存在一些问题:

您没有在展开中插入宏参数。这很容易出错，因为如果参数是表达式，可能会导致意外的操作符优先级问题。你应该写:
```
#define UINT8_HI(x) ((x) >> 4)
#define UINT8_LO(x) (((1 << 4) - 1) & (x))
```
首先存储最高有效部分的数组元素是违反直觉的。多精度算法通常将较大的值表示为数组，最不有效部分放在前面。

对于uint8_t这样的小类型，没有必要将其分成两半，因为有更大的类型可用。此外，必须传播前一个加法的进位。下面是一个简单得多的加法实现:

void uint8_add(uint8_t a, uint8_t b, uint8_t *res, int i) {
uint16_t result = a + b + res[i + 0]; // add previous carry
res[i + 0] = (uint8_t)result;
res[i + 1] = (uint8_t)(result >> 8); // assuming res has at least i+1 elements and is initialized to 0
}

对于乘法，您必须将每个数字的每个部分相乘的结果与结果数的适当选择的部分相加，将进位传播到更高的部分。

除法更难有效地实现。我建议你学习一个开源的多精度包，比如QuickJS的libbf.c.

要将其转换为uint64_t的数组，如果平台上可用，可以使用无符号128位整数类型(64位编译器gcc, clang和vsc都支持此类类型)。

下面是加法和乘法的简单实现:

#include <limits.h>
#include <stddef.h>
#include <stdint.h>
#define NB_CHUNK  16
typedef __uint128_t uint128_t;
typedef struct {
uint64_t chunk[NB_CHUNK];
} uint1024_t;
void uint0124_add(uint1024_t *dest, const uint1024_t *a, const uint1024_t *b) {
uint128_t result = 0;
for (size_t i = 0; i < NB_CHUNK; i++) {
result += (uint128_t)a->chunk[i] + b->chunk[i];
dest->chunk[i] = (uint64_t)result;
result >>= CHAR_BIT * sizeof(uint64_t);
}
}
void uint0124_multiply(uint1024_t *dest, const uint1024_t *a, const uint1024_t *b) {
for (size_t i = 0; i < NB_CHUNK; i++)
dest->chunk[i] = 0;
for (size_t i = 0; i < NB_CHUNK; i++) {
uint128_t result = 0;
for (size_t j = 0, k = i; k < NB_CHUNK; j++, k++) {
result += (uint128_t)a->chunk[i] * b->chunk[j] + dest->chunk[k];
dest->chunk[k] = (uint64_t)result;
result >>= CHAR_BIT * sizeof(uint64_t);
}
}
}

如果128位整数不可用，您的1024位类型可以实现为一个32位整数数组。下面是一个灵活的实现，为数组元素和中间结果提供了可选择的类型:

#include <limits.h>
#include <stddef.h>
#include <stdint.h>
#if 1  // if platform has 128 bit integers
typedef uint64_t type1;
typedef __uint128_t type2;
#else
typedef uint32_t type1;
typedef uint64_t type2;
#endif
#define TYPE1_BITS  (CHAR_BIT * sizeof(type1))
#define NB_CHUNK    (1024 / TYPE1_BITS)
typedef struct uint1024_t {
type1 chunk[NB_CHUNK];
} uint1024_t;
void uint0124_add(uint1024_t *dest, const uint1024_t *a, const uint1024_t *b) {
type2 result = 0;
for (size_t i = 0; i < NB_CHUNK; i++) {
result += (type2)a->chunk[i] + b->chunk[i];
dest->chunk[i] = (type1)result;
result >>= TYPE1_BITS;
}
}
void uint0124_multiply(uint1024_t *dest, const uint1024_t *a, const uint1024_t *b) {
for (size_t i = 0; i < NB_CHUNK; i++)
dest->chunk[i] = 0;
for (size_t i = 0; i < NB_CHUNK; i++) {
type2 result = 0;
for (size_t j = 0, k = i; k < NB_CHUNK; j++, k++) {
result += (type2)a->chunk[i] * b->chunk[j] + dest->chunk[k];
dest->chunk[k] = (type1)result;
result >>= TYPE1_BITS;
}
}
}

相关内容

最新更新

热门标签：