多精度算术：和算法
5.00/5 (1投票)
最简单的多精度加法算法。
本文是系列文章的一部分
多精度算术：一个有趣的爱好项目
引言
计算机以给定大小的数值进行运算。在撰写本文时，64 位算术在廉价硬件上很常见。但有一个固定的限制，如果你需要操作超过这个限制，你必须使用算法来指示你的计算机如何去做，这并不是“开箱即用”的。
在本系列文章中，我将分享我的知识以及一些代码。我的目标是编写一套最小化的、可工作的算术算法集，并尽可能地分享给其他程序员。我无意创建一套严肃的库，因为市面上已经有很多由专业数学家编写的优秀的 MP 库。
本系列分为四篇背景文章，最后一部分将解释如何使用我的代码，以防你想挑战我实现的简单算术运算的“参考”实现。
一个真正的 MP 库需要实现更多功能，需要处理负数和分数、多项式，而我将在这里涵盖最基本的四种算术运算。
长整数
不使用 阿拉伯数字 来传达“那个罐子里有 12334 块石头”的概念可能会很难。好吧，即使使用阿拉伯数字也很难，但至少你能对罐子里石头的数量级有一个概念，这比发送一张摆满石头的桌子图片要好，至少在我看来是这样，这就是阿拉伯数字如此优越的原因；那些伟大的古代数学家在没有阿拉伯数字的情况下发现了许多数学知识，他们有多么伟大。
顺便说一句，不用深入研究 位值制，我们可以看到“12334”本身不是一个数字，它是一串数字，以一种非常紧凑的形式表示一个数字，比一张摆满石头的桌子更容易让大脑处理。我们可以使用这个字符串来比较数字，并以比使用一堆石头（石头在古罗马拉丁语中翻译为 calculi）更简单的方式进行算术运算。
我们在童年时至少学习了四种处理长整数的算法：长加法、长减法、长乘法和长除法，虽然存在更好的算法，但我们将从这些开始。
简短说明
在本系列中，我将只关注无符号算术。
背景
加法是最基本的算术运算之一。要对一对数字进行加法，你只需要知道如何计数。计算 23+98 意味着在 23 之后计数 98。但当你需要计算 100 000 + 133 212 时，最好使用 十进制位值制 和一种需要计数到 10（以 10 为基数）的算法。我们将这种算法称为 **长加法算法**，我们使用它是因为它比计数所需的努力少。
现在，我们将定义一个基本操作，然后构建长加法算法。
十进制简单加法（我们的基本操作）
我将为两个数字 A 和 B 的和定义一个名称，其中 A 和 B 只有 1 位数字（数字可以有任何基数，例如 10、16、64、232）：我称之为 **简单加法**，在十进制中，你可以通过计数到 9 来计算它。
示例 1
3+3 = 从第一个数字 3 之后开始计数（是 4），数 3 次： 4, 5, 6：所以，结果是 6
示例 2：（如果结果需要多于 1 位数字会发生什么？）
溢出
8 + 4 = 从 9 开始计数 4 次： 9, 10（溢出：在纸上标记你 有进位 并从 1 继续计数），1, 2；所以，结果是 2 并且你有进位，进位（如果存在）总是 1。
简单加法的一些有趣性质
无论基数 **b** 是多少（例如 10），如果 A + B >= **b**，那么
- 我们有进位
- 结果有 2 位有效数字
- 最高有效数字为 1
- 最低有效数字小于任何一个操作数，换句话说，如果 A+B >= b，那么结果是两位数字 [1][R]，其中 R < A 且 R < B（如果你需要证明，请参阅附录 A1）。
对称地，如果 A + B < **b**，那么
- 我们没有进位
- 结果有 1 位数字
- 该数字将始终大于或等于任何一个操作数（附录 A1）
推论
- 长加法的结果比最长的加数多 0 或 1 位数字（便于为结果分配空间）。
- 保存结果的最高有效位数字将 **始终** 是 1 或 0，即使你将 1 加到简单加法结果的 2 位数字中（这意味着你总可以在不分配新数字的情况下处理前一项操作的进位），让我更好地解释一下
- 我们说过，如果 A+B 已经有进位，那么 A+B+1 的进位肯定为 0。例如：9 是十进制中最大的个位数，如果我们计算 9+9 得到
- 9+9 = 18,
- 我们发生溢出，因此最高有效数字是 1
- 最低有效数字是 8
- 你可以安全地在 8 上加 1 而无需再次溢出，因此最高位永远不会是 2。
- 我们说过，如果 A+B 已经有进位，那么 A+B+1 的进位肯定为 0。例如：9 是十进制中最大的个位数，如果我们计算 9+9 得到
注意：4+5=9，如果你加 1，就会产生进位，但结果仍然是 2 位数字，左边是 [1]，这意味着你可能需要将进位向左传播，这使得并行化加法算法变得困难。
- 重要：**如果最低有效数字小于任何一个操作数，那么我们就有进位**：如果我们对单个数字 [A] + [B] 进行加法，结果将是 [进位?] [R]，我们可以通过以下表达式推断进位的值
- 如果 R < A （R < B 也有效）则
carry = 1
else
carry = 0。
- 如果 R < A （R < B 也有效）则
这意味着我们可以忽略硬件进位标志，通过比较结果的低位数字和任何一个操作数来推断进位。
十进制长加法
为了对多于 1 位数字的数进行加法，我们将应用 **长加法算法**。
我将使用我们 6 岁时学会的算法，在 CPU 上实现长加法。
这是一个串行算法，我认为可以调整上述有趣的性质来创建一个并行算法，用于 SIMD 硬件，例如 x64 CPU 上的 MMX 或 多核 GPU，但这应该是一篇新文章，如果我以后写的话。
CPU 寄存器数字与十进制数字的类比
为了继续深入，我们需要将十进制数字与 CPU 寄存器数字（也称为 CPU 字）进行类比。从我们的角度来看，64 位字与十进制数字没有区别。
长加法算法对儿童和计算机来说是相同的，区别仅在于数字的大小和使用的基本简单加法，对我们人类来说是计数到十，而软件实现则使用 CPU 内置的加法指令。
数字计算机上的长加法（与十进制相同，但使用基数 2integer_register_size）
CPU 可以对寄存器大小的数字（字）进行加法。现代 CPU 可以在 1 个时钟周期内完成 32/64/… 位的加法。我们的方法是使用寄存器大小的数字和 CPU 进位标志，或进位检测算法。
算法
此处提出的算法是对实际算法的高层概述，想法是将 A+B 的长加法分解为一系列简单加法 An+Bn，从两个数字的最低有效数字开始，我们还需要在每次迭代中传播进位。可能存在的优化被完全忽略了。
进位检测变体
上述算法在“进位 = 进位标志”点有一个子算法，在这里我们检测前一次简单加法的进位。
CPU 具有检测和位标志来指示进位的硬件。
该 CPU 标志通常无法用于更高级的语言，但可以通过使用我们之前讨论过的属性来获得进位。
我知道有四种检测进位的方法。
- 使用硬件进位标志。
- 利用属性：(当进行 An+Bn 时) **当且仅当 (R < An 或 R < Bn) 时进位为 1**，这有点棘手，因为我们不再谈论 R = An+Bn，而是 R = An+Bn+previous_carry。因此，实际的测试将是：**carry_flag = if(previous_carry = 1) (R <= An ? 1 : 0), else if(previous_carry=0) (R < An) ? 1 : 0**，或者通过执行：**carry_flag = (R < An || R < Bn) ? 1 : 0**）这种选择似乎是最具可移植性的，也是我处理进位的首选方式，但从性能角度来看，我们应该衡量这两个子变体。一种 hacky 的方法（如果编译器支持）来避免分支可能是：**carry_flag = (R < An | R < Bn); **如果编译器生成预期的代码，这应该可以工作，单元测试（或查看反汇编代码）可以告诉你生成的代码是否实际有效。
- 使用一个数字系统，你保留最高位而不是进位（例如：在 32 位 CPU 上，将数字表示为 31 位字符串，这在语言不提供无符号类型时也可能有用），但这只是无用的复杂化。
- 如果高级语言提供至少比寄存器大小大 1 位的整数类型，我们可以使用该类型来存储 R，高级语言会自动将进位标志移动到高位字的最低位，我在 C 语言中举了两个使用第三种变例的例子。
/* First Example: Using union Using a UNION to avoid bitwise ops: we must pay attention to do right casts and to the endianness of the target hardware */ #include <stdio.h> /* this is not the complete algorithm, we are only showing the carry detection part in the evidenced lines suppose that word_t corresponds to the hardware CPU register_size and that dword_t corresponds to 2 times the register_size, NOTE THAT we just need a type having 1 bit more than the register_size for sum, but later on we might need 2 times the register_size to do multiplications */ typedef union _reg_long { dword_t DW; struct _reg_words { /*on big endian field order should be inverted I guess*/ word_t Lower; word_t Higher; } words; } reg_long; int main() { reg_log rl; word_t R; short prev_carry = 1; short carry; word_t An = ((word_t)-1); word_t Bn = 0; rl.DW = (dword_t)An + Bn + prev_carry; carry = rl.words.Higher; R = rl.words.Lower; printf("New carry: %d, Result:%x", carry, R); return 0; } /* Second Example: Using bitwise ops template for languages not having union types, but have some bitwise operations */ dword_t rl; word_t R; short prev_carry = 1; short carry; word_t An = ((word_t)-1); word_t Bn = 0; rl = (dword_t)An + Bn + prev_carry; carry = rl >> (8*sizeof(word_t)) ; R = rl & ((word_t)-1); printf("New carry: %d, Result:%x", carry, R);
使用 理解代码
为了完整性，我将提供该算法的 C 实现，请注意在不使用进位标志的情况下，如何使用高级语言检测溢出（进位）。接口比较粗糙，目的是提供一个用于构建更复杂事物的基本功能，它需要尽可能快，因此不浪费时间进行输入检查或创建漂亮的数据结构，这不是基本函数应该负责的。这可以被其他代码封装以生成更简洁的接口。
/*
Preconditions: Will not check array bounds on R, thus R must have space
to accommodate MAX(ASIZE, BSIZE)+1 unsigned reg_t.
This implementation requires the least significant digit is at index 0,
I did not optimize to keep code clear, but consider that second and third loop could be
divided into 2 loops: after the first time carry becomes zero
it will never be 1 again, thus we could eliminate the carry detection
and just memcpy the remaining bits of A into result.
*/
typedef unsigned long long reg_t;
reg_t LongSum(reg_t* A, reg_t ASize, reg_t * B, reg_t BSize, reg_t* R)
{
reg_t carry;
reg_t i;
reg_t min_a_b = ASize > BSize ? BSize : ASize;
carry = i = 0;
while (min_a_b > i)
{
R[i] = A[i] + B[i] + carry;
carry = R[i] < A[i] + carry ? 1 : 0;
++i;
}
while (ASize > i)
{
R[i] = A[i] + carry;
carry = R[i] == 0 ? 1 : 0;
++i;
}
while (BSize > i)
{
R[i] = B[i] + carry;
carry = R[i] == 0 ? 1 : 0;
++i;
}
if (carry)
R[i] = 1;
return i;
}
略微优化的“参考”版本
/* Will not check array bounds on R, thus must have space to
accommodate MAX(ASIZE,BSIZE)+1 unsigned ints
returns length of Result
*/
numsize_t LongSumWithCarryDetection
(reg_t* A, numsize_t ASize, reg_t * B, numsize_t BSize, reg_t* R)
{
register reg_t carry;
register numsize_t i;
reg_t min_a_b = ASize > BSize ? BSize : ASize;
carry = i = 0;
while (min_a_b > i)
{
R[i] = A[i] + B[i] + carry;
carry = R[i] < A[i] + carry ? _R(1) : _R(0);
++i;
}
while (ASize > i && carry == 1)
{
R[i] = A[i] + carry;
carry = R[i] == _R(0) ? _R(1) : _R(0);
++i;
}
while (ASize > i) /*just copy words*/
{
R[i] = A[i];
++i;
}
while (BSize > i && carry == 1)
{
R[i] = B[i] + carry;
carry = R[i] == _R(0) ? _R(1) : _R(0);
++i;
}
while (BSize > i)
{
R[i] = B[i];
++i;
}
if (carry)
{
R[i] = _R(1);
++i;
}
return (numsize_t)i;
}
我也尝试在 x64 汇编中做同样的事情，我不擅长汇编，但这是我的 x64 代码
.data
.code
;int LongSumAsm(
; qword * A, RCX
; int ASizeInQWords, RDX
; qword * B, R8
; int BSizeInQWords, R9
; qword * R ON STACK
; );
; RAX, RCX, RDX, R8, R9, R10, R11 VOLATILE
; RAX used to store cpu flags
LongSumAsm proc
;stack frame
push rbp
mov rbp, rsp
push rbx ;rbx is non volatile...
; the following three instructions move the minimum of ASize, BSize into R9
mov r11, rdx ; suppose minimum is ASize, thus Asize in R11
cmp rdx, r9 ; compare ASize, BSize
mov [rbp + 20o], rdx ;ASize in [rbp + 20o] that is second qword in shadow space
mov [rbp + 30o], r9 ;BSize in [rbp + 30o] that is third qword in shadow space
cmova r11, r9 ; if ASize > BSize , move BSize into R11
clc ;clear carry
lahf ;load flags in ah
; use shadow stack space to save inputs
mov r10, qword ptr [rbp + 60o] ; keep R in R10, R is in stack at rbp + 60o
xor r9, r9 ; zero in r9, from now on r9 will be
; the index on the arrays
lea r11, [r11*8]
A_And_B_Loop: ; while (ASize > i && BSize > i)
cmp r11, r9 ; compare MinSize (r11)
; and array index (r9)
jle A_And_B_Loop_Ends ; if index >= minsize quit loop
mov rbx, qword ptr [rcx + r9] ; A[index] in rbx
mov rdx, qword ptr [r8 + r9] ; B[index] in r11
sahf ; reload flags from rax
adc rbx, rdx ; add with carry A[index] + B[index]
lahf ; store flags into rax
mov qword ptr [r10+ r9], rbx ; save to R[index]
; the result of A[index] + b[index]
add r9, 8 ; increase index
jmp A_And_B_Loop ; continue loop
A_And_B_Loop_Ends:
xor r11, r11 ; zero into r11
mov rdx, [rbp + 20o] ; rdx = size A
lea rdx, [rdx*8]
A_Only_Loop: ; while (ASize > index)
cmp rdx, r9 ; compare ASize (rdx) and array index
jle A_Only_Loop_Ends ; if(index <= ASize) quit loop
mov rbx, qword ptr [rcx + r9] ; A[index] into rbx
sahf ; reload flags from rax into flags register
adc rbx, 0 ; add with carry A[index] + 0 (0 in r11)
lahf ; store flags into rax
jnc A_NoCarry_Loop ; carry is zero no need to
; do more additions
mov qword ptr [r10 + r9 ], rbx ; save to R[index]
; the result that's in rbx
add r9,8 ; increase index
jmp A_Only_Loop ; continue loop
A_NoCarry_Loop:
mov qword ptr [r10+ r9], rbx ; save to R[index]
; the result of A[index] + b[index]
add r9, 8 ; increase index
cmp rdx, r9 ; compare MinSize (rdx)
; and array index (r9)
mov rbx, qword ptr [rcx + r9] ; A[index] into rbx
jle A_Only_Loop_Ends ; if index >= minsize quit loop
jmp A_NoCarry_Loop ; continue loop
A_Only_Loop_Ends:
mov rdx, [rbp + 30o] ; move BSize into rdx
lea rdx, [rdx*8]
B_Only_Loop: ; while (BSize > index)
cmp rdx, r9 ; compare BSize (rdx) with array index
jle B_Only_Loop_Ends ; if index <= BSize quit loop
mov rbx, qword ptr [r8 + r9] ; B[index] into rbx
sahf ; reload carry flag into flag registers
adc rbx, r11 ; add with carry B[index] + 0 (0 in r11)
lahf ; store flags into rax
jnc B_NoCarry_Loop ; carry is zero no need to
; do more additions
mov qword ptr [r10 + r9], rbx ; save to R[index]
; the result that's in rbx
add r9, 8i ; increase index
jmp B_Only_Loop ; continue loop
B_NoCarry_Loop:
mov qword ptr [r10+ r9], rbx ; save to R[index] the result of B[index]
add r9, 8i ; increase index
cmp rdx, r9 ; compare BSIze (rdx) and array index (r9)
jle B_Only_Loop_Ends ; if index >= BSIZE quit loop
mov rbx, qword ptr [r8 + r9] ; rbx into result
jmp B_NoCarry_Loop ; continue loop
B_Only_Loop_Ends:
sahf ; reload flags into flag registers
jnc Prepare_to_return
mov r11, 1 ; add with carry 0 + 0 (last carry in r11)
mov qword ptr[r10 + r9], r11 ; store the last carry flag into R[index]
add r9, 8
; return value
Prepare_to_return:
mov rax, r9
shr rax, 3
; resume stack frame
pop rbx
mov rsp, rbp
pop rbp
ret
LongSumAsm endp
end
后果
这是汇编 x64 与 C 第二个“参考”版本的一些速度比较，我的 CPU 是 2011 年的四核 i7 3Ghz
C 版本以发布模式编译似乎相当不错，在小数字与大数字的加法中比汇编更快，但在所有其他测试中 ASM 比 C 快：这些数字必须谨慎看待，因为我不确定我是否能完全控制所有方面，也不能确定 C 算法是否与汇编完全相同，因为我是从头开始编写汇编的，而不是翻译 C 代码。ASM 比 C 快 50%（0.0125 秒 vs 0.018 秒），VC++2019 在 i7 (2011) 上。
但更换了编译器（VC++2022）和电脑（i9 2021），现在 ASM 在同一测试中比 C 慢。
在这里，它们几乎相等（0.64 秒 ASM vs 0.56 秒 C），但 C 版本在小数字对大数字时表现良好。
VC++2019 在 i7 (2011) 上
在 VC++2022 on i9 (2021) 上进行相同的测试，差异更小，它们处于同一水平，没有图片，但…相信我。
下面是最后一个测试，这是认真的，因为我们现在要对 2 个非常大的数字（每个 3300 万位）进行加法。
此测试中的 ASM 版本比 C 快 4 倍（VC++2019，i7 2011），也许我忘了在那次运行前启动优化器，我现在没有那台 CPU 了……我永远不会知道。
在 VC++2022 (i9 2021) 上进行相同的测试，由于新编译器的最新 CPU 利用率，C 版本有了巨大的改进。
关注点
有些事情值得注意，ASM 语言在使用高级语言（如 C）中不可用的硬件指令时似乎更快。这些优势很难衡量，因为 C 编译器通常很智能，能够生成机器代码来优化 CPU 的乱序执行等。对于涉及长加/减/乘/除的任意长整数，ASM 优于 C，但值得注意的是，C 也不算太差，表现相当好，我想说，有了新的 vc2022，编译器比 ASM 更快，考虑到我是在猜测进位而不是使用硬件进位标志。
第二点需要注意的是，我认为一个好的 C 编译器，对于给定的平台，可以提供映射到硬件指令（内在函数）的函数库，从而向开发人员公开与 ASM 相同的硬件，如果编译器足够智能，这就足以击败大多数手动编写的 ASM 实现，尤其是我自己写的 ASM 实现。
附录 A1
如果 A+B（A 和 B 是单个数字）发生溢出，在任何基数 b 下，结果都将小于任何一个操作数。
执行简单加法 **A**+**B** 会得到两个数字的结果：Result（最低有效数字）和 Carry（进位只能等于 1 或 0，它是最高有效数字）。
我们想证明，当 \(A+B \geq b\) 时，最低有效数字小于 A 和 B，我们可以换句话说：
直观的解释是，如果我们从 \(A\) 开始计数 \(B\) 次，当我们达到 \(b\) 时，我们至少消耗了 \(B\) 的 1 个单位，所以 \(B\) 将大于 \(Result\)，A 也一样，反之，如果 A + B < b，我们在不从 0 重新开始计数的情况下，将单位从 A 转移到 B，那么结果将大于两者。
我关于左箭头的直观想法是把数字想象成两个可以装 b-1 个水分子的小桶，如果你想把第 b 个分子倒入一个数字中，但这个数字满了，你必须倒空这个桶，剩余的分子必须添加到它左边的桶中。所以用这个比喻，试着想想当倒入 A 个分子 + B 个分子时，结果桶（最初是空的）发生了什么，因为结果没有被倒空（因为它能装下 A 和 B 的所有分子），所以结果中的水位比 A 和 B 中的水位高。
好的，我们来对右箭头进行严格证明，如果有人要求，我会提供左箭头的证明。
我首先展示以下等式
正如你通过仅仅看上面的表达式就可以证明其相等性一样，现在我们重新排列一下
好的，(2) 正确地告诉我们，当 \( A + B > b \) 时，结果是 2 位数字，我们可以重写 (2) 以体现其位值表示法，如下所示：
A + B = (1)\color{blue}{b^{1}} + \color{gray}{(B + A - b)}\color{blue}{b^{0}}
\end{equation}
请注意，(3) 中的灰色部分 \(\color{gray}{B + A - b} \geq 0\)，因为假设是 \(B + A \geq b\)，但我们也希望它小于 b，否则它将不适合单个数字。
所以我们想要 \(\color{gray}{B + A - b}< b\)（换句话说，我们希望 \(\color{gray}{B + A - b}\) 能够装入一个数字）；让我展示 \(0 \leq \color{gray}{B + A - b} < b\)
再次将 A 和 B 重写为与 b 的差值
这是正确的，现在灰色部分是一个正数（因为 A 和 B 都是单个数字，因此它们位于整数区间 \([0, b-1]\)）……简单来说，\(0 < i \leq b\) 且 \(0 < j \leq b\)
但我们也需要 \(i + j \leq b \)，幸运的是这是真的，因为
好的，现在所有这三个表达式为真 \(0 < i \leq b\) 且 \(0 < j \leq b\) 且 \( i + j \leq b \) 等价于
利用 (11)，接下来我们将证明
让我们通过用 (b-j) 和 (b-i) 替换 A 和 B 来做到这一点
由于 (11) 为真，那么 (16) 也为真，因此我们揭示了 (12) 为真。
好的，让我们从等式 (3) 继续
\(\color{gray}{(B + A - b)}\) 是 \(A + B\) 的最低有效数字，而 \(1\) 是最高有效数字：我们将最高有效数字重命名为 \(carry\)（我们现在也知道它的值总是 \(1\)），然后我们将 \(\color{gray}{(B + A - b)}\) 标记为 \(Result\)。
最初，我们想证明当假设成立时（假设是 \(A + B \geq b\)），\(Result\) 小于 A 和 B。
为了证明这个论点，只需要证明 (17) \(B + A - b < A\) 和 (18) \(B + A - b < B\)，方法是相同的，让我们看看情况 (18) \(B + A - b < B\)
由于 \(b > A\)，根据数字的定义，我们可以说 \(A - b < 0\)
因此，表达式 (18) 为真，因为 \(A - b < 0 \) **证明完毕**
现在我们必须证明 (17) \(B + A - b < A\)，让我重新排序一下 \(A + B - b < B\)
\(b - B\)，根据数字的定义，我们可以说 \(B - b < 0\)
因此，表达式 (17) \(B + A - b < A\) 为真，因为 \(B - b < 0\) **证明完毕**
历史
- 2017：开始写作
- 2019 年 2 月 4 日：几乎准备好发布
- 2022 年 6 月 11 日：在线发布