有没有办法让这个函数更快？（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里面也有另一个解决这个问题的实现，但是速度更慢。

你能帮我吗？

- Jonathan Sánchez

---

除法（和取模）与加法（和减法）相比非常低效。if (a+b > 9) carry = 1; else carry = 0; - pmg

一个明显可以提高性能的点就是使用数组元素值范围更多的值。目前，您从可能的256个值范围中只使用了10个值。我建议您还应该使用不同的数据类型，例如 uint_fast32_t，因为这样您就需要进行更少的“真正的”加法来执行加法操作。你也许想要查看一下在 .net core 中它是如何实现的，可以查看。 - Ackdari

你显然需要将你的十进制数转换为一百二十八进制数，但是如果你对同一数据进行多次加法运算，这样可以至少稍微加快速度。 - Ackdari

@JonathanSánchez，这样做会更好，因为你需要做更少的for循环迭代。这也是为什么在纸上进行十进制数字加法比二进制数字手算更容易的原因。 - Ackdari

@JonathanSánchez 你在这里做的是所谓的“二进制编码十进制”（BCD）。它是一种数据格式，已经存在了很长时间，因为它允许对例如财务数据进行精确计算。想象一下，在基于8位或16位指令的没有浮点处理的CPU上，这将是多么有用。正如其他人所暗示的那样，有很多方法可以将多个BCD十进制打包成更大的数据元素，例如uint32_t。 - Andrew Henle

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 的最低位未被设置，且 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

然后我们移除掉已经进位的所有10（即十位数）。

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

我们还跟踪进位，它可能会在两个位置发生。

- Veedrac

---

@JonathanSánchez 已更新 - Veedrac

这段求和代码中，如何将 uint8_t* 的两个数组改写为 uint32_t* 类型？@Veedrac - Marc

@Marc，这应该不会有太大的区别，只需要使用u32 SIMD数组并调整比较即可。 - Veedrac

@Marc 我不知道，也许你忘记改变移位量了。你是那个编写代码的人。 - Veedrac

谢谢，但是如果我使用uint32_t*，我需要更改我的向量类型，对吗？@Veedrac - Marc

Answer 2

加速优化的候选方案：

优化

确保您已经使用编译器启用了其速度优化设置。

restrict

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

使用restrict指示对Vout []的写入不会影响Vin1 []，Vin2 []。

// 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 是在整数计算中使用的“常用”类型。与使用 unsigned char CARRY 或 char CARRY 相比，使用 unsigned 或者来自<inttypes.h>中的 uint_fast8_t 更好。

% alternative

sum = a+b+carry; if (sum >= 10) { sum -= 10; carry = 1; } else carry = 0; 可以使用像@pmg这样的 % 方案。

---

注意: 我希望LongNumAddition1()返回最终进位。

- chux - Reinstate Monica

---

将 restrict 添加到 Vin1 和 Vin2 不会破坏将数字加到自身的操作吗？ - 0x5453

@0x5453 是的，restrict 会破坏它 - 很好的观察。但是如果 OP 通常正在寻找最佳速度，那么最好也提供一个 += 函数来处理该特殊情况。 - chux - Reinstate Monica

事实上，我知道restrict和const，但我没有意识到我可以使用它们，真的非常感谢！另外，使用uint_fast8_t会加快速度，这是我需要的。为什么会加速呢？我应该只更改进位到这个新结构体吗？还是我也可以更改向量以获得更多的优化？ - Jonathan Sánchez

好的，那我会尝试一下对进位进行处理。只要我的数组超过10000个位置，也许改变它并不是最好的选择。另外，在你上面的注释中，restrict是什么意思？如果我添加了restrict关键字，那么我就不能做Vin1+Vin2吗？还是说我不能做Vin1[i]+=1？ 除此之外，为了澄清你在答案中的注释，实际上原始版本确实返回了它，但只要我不用它做任何事情，我就避免使用它了。 - Jonathan Sánchez

迄今为止最大的优化是将这两个循环合并为一个。我也尝试了一些限定词和整数大小等，但对机器代码的帮助不大。 - Lundin

Answer 3

为了提高大数加法的速度，您应该将更多的十进制数字打包到数组元素中。例如：您可以使用 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的编译器探索器上查看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;
}

- chqrlie

---

carry = (r >= 10) 是一个测试。它可能有助于编译器避免发出分支指令，而选择条件移动。它是否比 OP 中的 carry = r/10 更好取决于确切的 CPU 类型。在一些具有硬件除法但没有条件移动的 CPU 上，除法会更快，在其他一些具有条件移动指令但没有硬件除法的 CPU 上，比较会更快。 - cmaster - reinstate monica

我真的很喜欢这个答案。无论如何，我已经尝试了这段代码，正如你可以在这里看到的那样，它有一个非常奇怪的输出。两个操作数各有25个位置，我想当进行LongNumConvert时，它们会转换为100。之后进行加法运算，结果是一个230个位置的数字，我需要转换这个数字还是发生了什么？ - Jonathan Sánchez

@cmaster-reinstatemonica: (c >= 10) 是一个比较操作。根据CPU的不同，它可能会生成分支，但在当前架构上并不会。r/10通常不会编译成硬件除法，而是编译成乘法和移位操作。但正如你所说，这取决于CPU。 - chqrlie

哇，太厉害了，现在它确实起作用了！ 我想问你一个最后的问题，如果A和B都是长度为N，那么LongNumConvert需要做出任何更改吗？ - Jonathan Sánchez

再次您好，我进行了更多的测试，发现它以不同的顺序执行加法。例如：567854513218631683125312535 + 1 应该是 667854513218631683125312535。但实际输出是：567854513218631683225312535。可以看到，它将数字添加到了最后一组9中：12变成了22。如何解决这个问题呢？ 谢谢！ - Jonathan Sánchez

Answer 4

在没有特定系统的情况下讨论手动优化通常是没有意义的。假设您有一些主流的32位带有数据缓存、指令缓存和分支预测的处理器，则：

- 避免多次循环。您应该将它们合并成一个循环，从而获得重大性能提升。这样您就不必多次触及同一内存区域，也可以减少总分支数。程序必须检查每个`i

- Lundin

Answer 5

第一个循环

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;
}

这段代码包含一个循环依赖，它实际上使整个循环串行化（考虑将1加到99999999999999999，只能逐位计算）。循环依赖是现代计算机科学中最大的难题之一。

因此，第一个版本更快 - 它部分矢量化了。其他版本并非如此。

如何避免循环依赖？

由于计算机是基于2进制的设备，因此在10进制算术方面表现不佳。它不仅浪费空间，还会在每个数字之间创建人为的进位依赖关系。

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

不幸的是，C编译器仍然很难生成具有进位传播的有效循环。因此，例如libgmp不是使用C而是使用汇编语言实现大数加法，使用ADC（带进位加法）指令。顺便说一下，libgmp可以直接替换项目中许多大数算术函数。

- rustyx

---

正如你所说，我尝试通过将代码分成两个循环来进行向量化...事实上，现在我发现我可以将其分成三个循环，并向量化第三个循环。然后更改是将最后的Vout [i] = Vout [i]% 10移动到另一个for中，而不是第二个for中。我的主要问题现在仍然存在于第二个循环中，在那里我添加进位并计算下一个进位。我已经阅读了你链接的SSE问题帖子，但我无法提取如何执行部分字算术的信息。此外，我还看到了libgmp代码用于[加法进位]（https://github.com/haiku/buildtools/blob/master/gcc/gmp/mpn/cray/add_n.c）。 - Jonathan Sánchez

由于我在C语言方面是新手，因此我没有像他们一样深入理解低级操作以至于能够用它们代替我的代码...你能否帮我解释一下关于2进制表示的段落？我已经明白了我必须将我的数字从10进制转换为2进制，但我不理解uint64_t部分...你是指将8个“数字”组合成一个uint64吗？ 谢谢。 - Jonathan Sánchez

那是Cray版本。它只是简单地传递第64位，而不会在具有硬件进位标志的x86上执行。x86版本在这里。我没有看到任何可以有意义地改进您的第二个循环的方法，除非更改数据表示。您可以尝试展开，即在单次迭代中计算4..8个数字，但我想那就是全部了...通过数据表示，我指的是二进制表示，例如123表示为1*10^2 + 2*10 + 3 = 二进制01111011。 - rustyx

我的处理器具有x86-64 ISA，因此我将使用您在此回复中发布的新链接。但无论如何，我真正不理解的是我需要在我的数据表示中更改什么...正如您所说，123将以二进制形式表示为01111110，我可以在我的数组中获取3或4个位置并将它们打包成二进制并保存在uint64变量中吗？这就是我不明白的地方...对此我很抱歉。 - Jonathan Sánchez

是的，您可以在32位中打包9个数字或在64位中打包19个数字，并且每次迭代轻松添加9或19个数字，但这不会真正成为基于2进制。在基于2进制的情况下，您需要将所有数字打包到所需的位数中。libgmp有一个名为mpn_set_str的函数可用于此。 - rustyx