有没有办法让这个功能更快？ （C）

Question

我有一个 C 代码，它以与人类相同的方式进行添加，因此，例如，如果我有两个数组 A[0..n-1] 和 B[0..n-1]，则该方法将执行 C[0]=A[0]+B[0]、C[1]=A[1]+B[1].. .

即使解决方案使用内在函数，我也需要帮助来加快此功能。

我的主要问题是我有一个非常大的依赖性问题，因为迭代 i+1 取决于迭代 i 的进位，只要我使用基数 10。所以如果 A[0]=6 和 B[0]=5 , C[0] 必须是 1 并且我有一个 1 的进位用于下一次添加。

我能做的更快的代码是这个：

void LongNumAddition1(unsigned char *Vin1, unsigned char *Vin2,
                      unsigned char *Vout, unsigned N) {
    for (int i = 0; i < N; i++) {
        Vout[i] = Vin1[i] + Vin2[i];
    } 

    unsigned char carry = 0;

    for (int i = 0; i < N; i++) {
        Vout[i] += carry;
        carry = Vout[i] / 10;
        Vout[i] = Vout[i] % 10;
    }
}

但我也尝试了这些方法，结果速度较慢：

void LongNumAddition1(unsigned char *Vin1, unsigned char *Vin2,
                      unsigned char *Vout, unsigned N) {
    unsigned char CARRY = 0;
    for (int i = 0; i < N; i++) {
        unsigned char R = Vin1[i] + Vin2[i] + CARRY;
        Vout[i] = R % 10; CARRY = R / 10;
    }
}

void LongNumAddition1(char *Vin1, char *Vin2, char *Vout, unsigned N) {
    char CARRY = 0;
    for (int i = 0; i < N; i++) {
        char R = Vin1[i] + Vin2[i] + CARRY;
        if (R <= 9) {
            Vout[i] = R;
            CARRY = 0;
        } else {
            Vout[i] = R - 10;
            CARRY = 1;
        }
    }
}

我一直在谷歌研究，发现了一些与我实现的类似的伪代码，在 GeeksforGeeks 中也有另一个实现这个问题的方法，但速度也较慢。

你能帮帮我吗？

Answer 1

如果不想改变数据的格式，可以试试SIMD。

typedef uint8_t u8x16 __attribute__((vector_size(16)));

void add_digits(uint8_t *const lhs, uint8_t *const rhs, uint8_t *out, size_t n) {
    uint8_t carry = 0;
    for (size_t i = 0; i + 15 < n; i += 16) {
        u8x16 digits = *(u8x16 *)&lhs[i] + *(u8x16 *)&rhs[i] + (u8x16){carry};

        // Get carries and almost-carries
        u8x16 carries = digits >= 10; // true is -1
        u8x16 full = digits == 9;

        // Shift carries
        carry = carries[15] & 1;
        __uint128_t carries_i = ((__uint128_t)carries) << 8;
        carry |= __builtin_add_overflow((__uint128_t)full, carries_i, &carries_i);

        // Add to carry chains and wrap
        digits += (((u8x16)carries_i) ^ full) & 1;
        // faster: digits = (u8x16)_mm_min_epu8((__m128i)digits, (__m128i)(digits - 10));
        digits -= (digits >= 10) & 10;

        *(u8x16 *)&out[i] = digits;
    }
}

这是每个数字约 2 条指令。您需要添加代码来处理尾端。

---

这是算法的运行过程。

首先，我们将我们的数字与上次迭代的进位相加：

lhs           7   3   5   9   9   2
rhs           2   4   4   9   9   7
carry                             1
         + -------------------------
digits        9   7   9  18  18  10

我们计算哪些数字会产生进位（≥10），哪些数字会传播（=9）。无论出于何种原因，对于 SIMD，true 为 -1。

carries       0   0   0  -1  -1  -1
full         -1   0  -1   0   0   0

我们将carries 转换为整数并将其移位，还将full 转换为整数。

              _   _   _   _   _   _
carries_i  000000001111111111110000
full       111100001111000000000000

现在我们可以将这些加在一起来传播进位。注意只有最低位是正确的。

              _   _   _   _   _   _
carries_i  111100011110111111110000
(relevant) ___1___1___0___1___1___0

有两个指标需要注意：

1. carries_i 设置了最低位，digit ≠ 9。这个方格有进位。

2. carries_i 的最低位 un 已设置，digit = 9。有一个进位 这个方格，重置位。

我们用(((u8x16)carries_i) ^ full) & 1 计算，然后添加到digits。

(c^f) & 1     0   1   1   1   1   0
digits        9   7   9  18  18  10
         + -------------------------
digits        9   8  10  19  19  10

然后我们把10s去掉，已经全部携带了。

digits        9   8  10  19  19  10
(d≥10)&10     0   0  10  10  10  10
         - -------------------------
digits        9   8   0   9   9   0

我们还跟踪执行，这可能发生在两个地方。

Answer 2

提高速度的候选人：

优化

确保您已启用编译器的速度优化设置。

restrict

编译器不知道更改Vout[] 不会影响Vin1[], Vin2[]，因此在某些优化中受到限制。

使用restrict 表示Vin1[], Vin2[] 不受写入Vout[] 的影响。

// void LongNumAddition1(unsigned char  *Vin1, unsigned char *Vin2, unsigned char *Vout, unsigned N)
void LongNumAddition1(unsigned char * restrict Vin1, unsigned char * restrict Vin2,
   unsigned char * restrict Vout, unsigned N)

注意：这会限制调用者使用与Vin1, Vin2 重叠的Vout 调用函数。

const

还可以使用const 来帮助优化。 const 还允许将 const 数组作为 Vin1, Vin2 传递。

// void LongNumAddition1(unsigned char * restrict Vin1, unsigned char * restrict Vin2,
   unsigned char * restrict Vout, unsigned N)
void LongNumAddition1(const unsigned char * restrict Vin1, 
   const unsigned char * restrict Vin2, 
   unsigned char * restrict Vout, 
   unsigned N)

unsigned

unsigned/int 是用于整数数学的“goto”类型。不要使用unsigned char CARRY 或char CARRY，而是使用unsigned 或uint_fast8_t 来自<inttypes.h>。

% 替代

sum = a+b+carry; if (sum >= 10) { sum -= 10; carry = 1; } else carry = 0;@pmg 之类的。

---

注意：我希望 LongNumAddition1() 返回最后的进位。

Answer 3

在没有考虑特定系统的情况下讨论手动优化总是毫无意义的。如果我们假设你有某种具有数据缓存、指令缓存和分支预测的主流 32 位，那么：

• 避免多个循环。您应该能够将它们合并为一个，从而获得重大的性能提升。这样您就不必多次触摸同一个内存区域，并且您将减少分支的总数。每个i < N 都必须由程序检查，因此减少检查量应该会提供更好的性能。此外，这可以提高数据缓存的可能性。

• 对支持的最大对齐字长执行所有操作。如果你有 32 个苦味，你应该能够让这个算法一次处理 4 个字节，而不是一个字节一个字节。这意味着以某种方式逐个字节地换出 memcpy 的分配，一次执行 4 个字节。这就是库质量代码的作用。

• 正确限定参数。你真的应该熟悉const 正确性这个术语。 Vin1 和 Vin2 没有改变，所以这些应该是 const，不仅仅是为了性能，而是为了程序的安全性和可读性/可维护性。

• 同样，如果您可以保证参数没有指向重叠的内存区域，您可以restrict 限定所有指针。

• 除法在许多 CPU 上是一项昂贵的操作，因此如果可以更改算法以摆脱 / 和 %，那么就这样做。如果算法是逐字节完成的，那么你可以牺牲 256 字节的内存来保存一个查找表。

  （这假设您可以在 ROM 中分配这样一个查找表而不引入等待状态依赖等）。

• 将进位更改为 32 位类型可能会在某些系统上提供更好的代码，而在其他系统上则更差。当我在 x86_64 上尝试这个时，它通过一条指令给出了稍差的代码（非常小的差异）。

Answer 4

第一个循环

for (int i = 0; i < N; i++) {
    Vout[i] = Vin1[i] + Vin2[i];
} 

由编译器自动矢量化。但是下一个循环

for (int i = 0; i < N; i++) {
    Vout[i] += carry;
    carry = Vout[i] / 10;
    Vout[i] = Vout[i] % 10;
}

包含一个loop-carried dependence，它本质上是对整个循环进行序列化（考虑将 1 加到 99999999999999999 - 它只能按顺序计算，一次一位数）。循环携带依赖是现代计算机科学中最令人头疼的问题之一。

这就是第一个版本更快的原因 - 它是部分矢量化的。任何其他版本都不是这种情况。

如何避免循环携带的依赖？

作为 base-2 设备的计算机在使用 base-10 算术方面是出了名的糟糕。不仅浪费空间，还会在每个数字之间造成人为的进位依赖。

如果您可以将数据从 base-10 转换为 base-2 表示，那么机器添加两个数组将变得更容易，因为机器可以轻松地在一次迭代中执行多个位的二进制加法。例如，对于 64 位机器，性能良好的表示可能是 uint64_t。请注意，对于 SSE，带有进位的流式添加仍然存在问题，但那里也存在一些选项。

不幸的是，C 编译器仍然很难通过进位传播生成有效的循环。出于这个原因，例如libgmp 不是在 C 中而是在汇编语言中使用 ADC（带进位相加）指令实现 bignum 加法。顺便说一句，libgmp 可以直接替代您项目中的许多 bignum 算术函数。

Answer 5

为了提高 bignum 加法的速度，您应该将更多的十进制数字打包到数组元素中。例如：您可以使用uint32_t 代替unsigned char 并一次存储9 个数字。

另一个提高性能的技巧是避免出现分支。

这是未经测试的代码修改版本：

void LongNumAddition1(const char *Vin1, const char *Vin2, char *Vout, unsigned N) {
    char carry = 0;
    for (int i = 0; i < N; i++) {
        char r = Vin1[i] + Vin2[i] + CARRY;
        carry = (r >= 10);
        Vout[i] = r - carry * 10;
    }
}

这是一次处理 9 位数字的修改版本：

#include <stdint.h>

void LongNumAddition1(const uint32_t *Vin1, const uint32_t *Vin2, uint32_t *Vout, unsigned N) {
    uint32_t carry = 0;
    for (int i = 0; i < N; i++) {
        uint32_t r = Vin1[i] + Vin2[i] + CARRY;
        carry = (r >= 1000000000);
        Vout[i] = r - carry * 1000000000;
    }
}

你可以在GodBolt's Compiler Explorer上查看gcc和clang生成的代码。

这是一个小测试程序：

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

int LongNumConvert(const char *s, uint32_t *Vout, unsigned N) {
    unsigned i, len = strlen(s);
    uint32_t num = 0;
    if (len > N * 9)
        return -1;
    while (N * 9 > len + 8)
        Vout[--N] = 0;
    for (i = 0; i < len; i++) {
        num = num * 10 + (s[i] - '0');
        if ((len - i) % 9 == 1) {
            Vout[--N] = num;
            num = 0;
        }
    }
    return 0;
}

int LongNumPrint(FILE *fp, const uint32_t *Vout, unsigned N, const char *suff) {
    int len;
    while (N > 1 && Vout[N - 1] == 0)
        N--;
    len = fprintf(fp, "%"PRIu32"", Vout[--N]);
    while (N > 0)
        len += fprintf(fp, "%09"PRIu32"", Vout[--N]);
    if (suff)
        len += fprintf(fp, "%s", suff);
    return len;
}

void LongNumAddition(const uint32_t *Vin1, const uint32_t *Vin2,
                     uint32_t *Vout, unsigned N) {
    uint32_t carry = 0;
    for (unsigned i = 0; i < N; i++) {
        uint32_t r = Vin1[i] + Vin2[i] + carry;
        carry = (r >= 1000000000);
        Vout[i] = r - carry * 1000000000;
    }
}

int main(int argc, char *argv[]) {
    const char *sa = argc > 1 ? argv[1] : "123456890123456890123456890";
    const char *sb = argc > 2 ? argv[2] : "2035864230956204598237409822324";
#define NUMSIZE  111  // handle up to 999 digits
    uint32_t a[NUMSIZE], b[NUMSIZE], c[NUMSIZE];
    LongNumConvert(sa, a, NUMSIZE);
    LongNumConvert(sb, b, NUMSIZE);
    LongNumAddition(a, b, c, NUMSIZE);
    LongNumPrint(stdout, a, NUMSIZE, " + ");
    LongNumPrint(stdout, b, NUMSIZE, " = ");
    LongNumPrint(stdout, c, NUMSIZE, "\n");
    return 0;
}