多精度算术： 乘法算法
1.00/5 (1投票)
乘法算法中最简单的一种，对于不太大的数来说足够快。
本文是系列文章的一部分
多精度算术：一个有趣的爱好项目
教会数字计算机长乘法算法
这是第三部分，我们将教会我们的计算机长乘法。在这篇文章中，我们将重点关注朴素的小学算术算法。
存在许多不同的长乘法算法，其中最简单的是小学算术算法。在计算机上正确实现它非常重要，因为尽管大O符号告诉我们这是一个 O(n2) 的算法，但它的每个基本部分都非常快，不需要额外的内存分配，因此对于较小的 n 值（假设小于一百个 64 位字），它比所有“最快”的算法都有优势。
存在更好的算法，但是小学算术算法在某个阈值之前比任何其他算法都快。换个角度来说：任何更快的算法在 n 低于给定阈值时都会比它慢；因此，将渐进速度更快的算法与“小学”算法混合使用是有益的，当 n 低于一个精确选择的阈值时；在源代码中，您会发现 Karatsuba 算法与这些算法结合使用时会更快。
所以，除了最基础之外，这个算法也值得最高级别的优化，您可能会在任何分治算法的叶子层中使用它，这意味着它将被大量使用。
我希望很快能写一些关于更快的算法的内容。
好了，说了这么多，让我们回顾一下如何手动进行，然后将其移植到计算机。本章分为三个部分：第一部分，我们将找到一个基本运算，然后分析算法，最后讨论一些关于硬件实现的内容。
1. 使用位值表示法进行简单乘法，基数 10
让我们从“简单乘法”（或基本乘法，随你怎么称呼）开始，即 A 和 B 两个数字的乘法，其中 A 和 B 都只有一个数字。这将是我们构建长乘法的基础。
要进行简单乘法，您可以重复求和，例如 3x4 可以计算为 3+3+3+3（4 次）。
……或者只是查找预先计算好的基数 b 的毕达哥拉斯表。
示例：使用毕达哥拉斯表进行 3x4 乘法，基数 10：搜索第 3 行和第 4 列的交叉点，结果是 12。
单个数字乘以单个数字可能得到两位数（绝不会超过两位数）。
现在我们有了一个基本运算来构建长乘法，我们称之为“简单乘法”。手动进行时我们将使用毕达哥拉斯表，在计算机硬件上则使用 CPU 的 MUL 指令（或同等指令）。
2. 长整数乘法
长乘法算法如下：你有一个数字字符串 A，一个数字字符串 B，对于字符串 A 和 B 的每个位置 i、j，将 A[i] * B[j] 相乘，然后将部分结果写在位置 i+j，最后将部分结果相加。
3. 长整数乘法与数字计算机
我们用纸和笔回顾了十进制数字的算法。CPU 版本完全相同。我们做了一个类比，CPU 寄存器就像一个十进制数字，并且我们使用等效的硬件指令（即乘法指令）而不是毕达哥拉斯表查找。
有使用基数 2 的实现（因为 1 位毕达哥拉斯表非常简单），但我的推测是这种方法效率低下（除非在硬件中实现）。
我现在尝试展示为什么使用最大的 CPU 寄存器比使用较小的寄存器或更差的单比特更好：例如，我将比较基数 4 和基数 16 的乘法。
考虑十六进制 E * D 的乘法（十进制结果是 182）。
在基数 16 下，我们进行一次简单乘法：ExD=十进制 182
在基数 4 下，相同的乘法是 32*31，因此我们必须计算：(1*2),(1*3),(3*2),(3*3)，得到它们各自的结果，将部分结果放入矩阵中，然后计算长加法。
推测：故事的寓意是：将比特数翻倍，我们可以将大数的算法工作量减少至少 4 倍。
推论：64 位 CPU 使我们的算法比 32 位 CPU（在处理大数时）的工作量减少至少 4 倍。
推论：即使每个 64 位乘法比 32 位乘法（并非如此，它消耗更多电能，但时间不会是 4 倍）成本高 4 倍，我们仍然会因使用 64 位而获得巨大收益，因为需要的长加法次数更少，内部循环迭代次数更少，更不用说基数 2 了。
推论：为了乘 32 位数，最好使用 32 位操作，因为单个 32 位操作的成本低于 64 位操作。
推论：任何基于基数 2 的实现都无法与任何 n 位软件长乘法竞争。
免责声明：上述推测并非基于证据。
算法
基本算法与我们在学校学到的相同，它既适用于计算机，也适用于纸笔计算。一如既往，我在此提供算法的高级概述。它比其他更快的算法有一个很大的优点：其内存需求易于预测，除了结果数组之外，不需要额外的分配。
在将其转化为软件之前的一些考虑
要将算法转化为程序，我们需要稍微思考一下如何有效地实现它。
在加法和减法方面很简单，但现在需要我们高级语言提供更多支持来有效地利用硬件。我的推测是，我们使用的寄存器越大，乘法就越快。为了证伪这个假设，我们需要检查代码的两个变体，一个利用 CPU 的全部能力，另一个利用我们使用的语言可以提供的功能。
让我稍微澄清一下：假设在一个 64 位计算机上工作：并非所有高级语言都有一个乘法函数可以从一次乘法中返回两个 64 位值（正如 64 位 CPU 通常所做的那样），也不是它们都有一个 128 位整数类型。
我们必须在两个极端之间选择一种策略
- 牺牲可移植性，编写特定于机器的代码，为每台特定机器使用最快的指令
- 牺牲性能，使用你的高级语言提供的任何功能
我的初衷是选择方案 2：我们只想了解它是如何工作的，如果您想要更高的吞吐量，市面上已经有真正的专家创建的 MP 库。在我决定使用一点汇编语言之后，这是因为我的编译器在 int 大小上有一些限制，我并没有真正使用“最快的机器指令”，但也不是纯粹的高级语言解决方案。
让我们从一个可移植的算法开始，我将使用高级语言（例如 C 语言）中易于获得的结构，并尽量保持算法的简单性。只有一小部分代码需要访问“秘密”硬件功能（将两个寄存器大小的 int 相乘并返回 2 个寄存器 int），因为我的编译器不允许我使用 128 位整数类型。
使用代码理解代码
致读者的说明：在这里，我解释了我如何决定实现这个东西，这不是唯一的方法。
如果我们有一个支持 128 位整数的编译器，我会这样做：我会使用语言提供的两种无符号整数数据类型：一种我们称之为“word 类型”，然后我们需要另一种无符号类型，它的大小是前者的两倍，我们称之为“double word”。理想情况是，根据我在前几段中做出的推测，word 类型具有最大的 CPU 寄存器大小。
一些语言提供了联合类型（union），使用联合类型可以使事情更简单，否则我们将使用位移来获取乘法结果的高位字。
在这篇文章中，我将使用 union，一个不错的起点可能是这样的（注意字节序）
/*
We must provide this header for each target platform, or use conservative settings,
the rest of the algorithm will not change.
A configuration script could figure out what native type to use for the target machine
*/
typedef uint32_t _word_t; /*Suppose hour CPU is 32 bit:
ensure we use CPU sized words, in case we have a 64 bit CPU
*/
typedef uint64_t _dword_t; /*machine double words*/
typedef union
{
_dword_t dword;
struct
{
_word_t L; /*little endian*/
_word_t H;
} Pair;
} _word_struct;
但是，我不得不与这样一个事实作斗争：我的 VC++ 编译器不支持（或曾经不支持）128 位类型，但我的机器是 x64 机器，这迫使我使用 32 位 word 类型（在乘法中），这并不理想。
作为折衷，我提出了这个解决方案：编写一个单一的算法，然后将毕达哥拉斯查找委托给另一个函数，该函数可以用 C 编写，或者也许用汇编语言编写（汇编例程很可能不会内联调用，但可能比降级的半字乘法更快）。
结果如下
#include "BigIntSimple.h"
#include <assert.h>
/* here is the pythagorean function callback declaration,
the portable algorithm will call this without knowing
which is the actual implementation
*/
typedef reg_t (*_mul_func)(reg_t A, reg_t B, reg_t* Hiword) ;
/* here i have a small x64 assembler routine implementation
of the above callback, which unfortunately is not going to be
inlined though, so you might replace it with platform specific intrinsics
*/
EXTERN reg_t cpu_multiply(reg_t A, reg_t B, reg_t * high);
#ifdef NO_DWORD_INTS
/*
this version is a C version for compilers with no dword ints, as mine,
so this will software emulate the 64 bit multiplication using 32 bit ints
(but please don't think "64 bit" and "32 bit" as absolute values, it depends
on underlying hardware, it could be "128 bit word" and "256 bit dword")
the definition / declaration are private for this module
*/
static reg_t c_multiply(reg_t A, reg_t B, reg_t * high){ ... }
#else
/* BEWARE: this version is for compilers with dword ints which I did not test...
below is the definition of a version working for compilers that support
a "dword" integer type (dword i mean an int type twice as big as the register size)
the definition / declaration are private for this module
*/
static reg_t c_multiply(reg_t A, reg_t B, reg_t * high){ ... }
#endif
/*
Here we have a private multiplication function that can receive an additional
parameter containing the callback function for the multiplication
*/
static numsize_t LongMultiplicationPortable
(reg_t *A, numsize_t m, reg_t *B, numsize_t n, reg_t * R, _mul_func mfunc)
{...}
/* these are the public interface implementation */
/* first we define a version using the assembler version
"pythagorean lookup" as callback
*/
numsize_t LongMultiplication(reg_t *A, numsize_t m, reg_t *B, numsize_t n, reg_t * R)
{
return LongMultiplicationPortable(A, m, B, n, R, cpu_multiply);
}
/* then we define a version using one of the 2 "C" version callbacks
*/
numsize_t LongMultiplicationNoAssembly
(reg_t *A, numsize_t m, reg_t *B, numsize_t n, reg_t * R)
{
return LongMultiplicationPortable(A, m, B, n, R, c_multiply);
}
这是毕达哥拉斯查找的小型 x64 汇编版本，它（我们稍后会看到）使 C 代码的运行速度提高了 2.5 倍。
.code
cpu_multiply proc
mov rax, rdx
mul rcx
mov qword ptr[r8], rdx
ret
cpu_multiply endp
end
这是 LongMultiplicationPortable 函数之一的完整 C 代码列表，在源代码中您会找到不同的版本。
static numsize_t LongMultiplicationPortableV2
(reg_t *A, numsize_t m, reg_t *B, numsize_t n, reg_t * R, _mul_func mfunc)
{
numsize_t outBuffSize = m + n, k;
reg_t hiword, loword;
int carry;
for (k = 0; k < outBuffSize; ++k)
R[k] = 0; /*reset output array*/
for (numsize_t j = 0; j < m; j++)/* read left number*/
for (numsize_t i = 0; i < n; i++) /* read right number*/
{
k = i + j;
loword = mfunc(A[j], B[i], &hiword);
R[k] += loword;
carry = R[k] < loword; /* detect the carry */
++k; /* next digit */
R[k] += hiword + carry;
carry = (R[k] < hiword) | (carry & (R[k] == hiword)); /* ugly carry
detection */
while (carry) /* that is unlikely to loop too much */
{
R[++k] += 1; /*sum the previous carry to Result[i+j+k] */
carry = R[k] == _R(0);
}
}
while (outBuffSize > 0 && R[outBuffSize - 1] == 0)
outBuffSize--;
return outBuffSize;
}
为了让你理解上述代码，还需要注意一些事项：A 和 B 是输入的长整数，reg_t 类型的大小与机器寄存器大小相同，m 是 A 的大小，n 是 B 的大小，我们还接收一个预先分配的 R（结果）数（调用者应分配它，考虑到它可能包含多达 n+m 位数）。没有输入检查，没有前置条件检查，没有干净的数据结构，纯粹而原始，这不应该是一个 API，这是一个 API 将会封装的基本部分，它也将被长除法算法使用，所以我们不想在该级别进行昂贵的操作。
请记住，A 和 B 数字的最低有效数字在索引 0 处，最高有效数字在索引 A[m-1] 和 B[n-1] 处。
假定 A 和 B 的最高有效数字不为 0，但算法仍然有效，只是浪费时间和内存。
我还在匆忙中创建了一个汇编 x64 版本，我已在一定程度上进行了测试，看起来它还在工作，然后我并不满意，于是创建了一个新版本，它稍微好一些，我将它粘贴在下面，任何 x64 专家都会做得比这更好，仅供参考。
.data
.code
;int LongMulAsm2(
; qword * A, RCX
; int ASizeInQWords, eDX
; qword * B, R8
; int BSizeInQWords, R9d
; qword * R ON STACK
; );
; RAX, RCX, RDX, R8, R9, R10, R11 VOLATILE
LongMulAsmVariant_1 proc
; mul by zero ret 0
xor rax, rax
test rdx, rdx
jz immediate_ret
test r9,r9
jz immediate_ret
push rbp
mov rbp, rsp
push r12
push r13
push r14
push r15
mov r10, rdx
mov r11, [rbp + 60o] ; r11 <- R
;reset R array
add rdx, r9
xor rax,rax
_reset_loop:
mov [r11+rdx*8], rax
sub edx, 1
jns _reset_loop
xor r12, r12
_outer_loop:
cmp r12, r10 ; r12 = j
je _outer_loop_end
; it is unlikely that one digit is zero so... don't optimize that case
xor r13,r13
; k= i + j
mov r15, [rcx + r12 * 8]
_inner_loop:
cmp r13, r9
je _inner_loop_end
lea r14, [r11 + r12 * 8]
; pointer to result
; multiply *A * *B
mov rax, [r8 + r13 * 8]
lea r14, [r14 + r13*8] ;randomly mixing instructions
;in the hope to achieve better results
mul r15
add rax, [r14] ; R[i+j] += loword (first part)
adc rdx, [r14+8] ; R[i+j+1] += highword + carry (first part)
mov [r14], rax ; R[i+j] += loword (second part)
mov [r14+8], rdx ; R[i+j+1] += highword + carry (second part)
jnc _inner_loop_continue
lea r14, [r14+8]
_carry_loop:
lea r14, [r14+8]
mov rax, [r14]
add rax, 1
mov [r14], rax
jnc _carry_loop_ends
jmp _carry_loop
_carry_loop_ends:
_inner_loop_continue:
inc r13
jmp _inner_loop
_inner_loop_end:
_outer_loop_continue:
add r12,1
jmp _outer_loop
_outer_loop_end:
;compute size of result, the number of iteration will be 1 or 2 if
;input numbers have no leading zeroes
lea rax, [r11 + r10 * 8]
lea rax, [rax + r9 * 8]
_compute_size_loop:
cmp r11, rax ; pointer to begin of array == pointer to end of array ?
jg _end_compute_size_loop
sub rax, 8 ; pointer to end --
mov r10, [rax]
test r10, r10 ; *pointer_to_end == 0 ?
jz _compute_size_loop
_end_compute_size_loop:
add rax, 8; need size not index
sub rax, r11 ; distance in bytes
shr eax, 3 ; distance in words , this is the return value
; resume stack frame & non volatile
pop r15
pop r14
pop r13
pop r12
mov rsp, rbp
pop rbp
immediate_ret:
ret
LongMulAsmVariant_1 endp
end
测试实现
在不过多深入细节的情况下，我们基本上需要一些测试来验证任何实现的正确性，在源代码中，您会找到其中一些测试。我在此列出一些可能性。
- 检查 A * B = B * A，交换律
- 检查 A * 1 = A。（1 是单位元）
- 检查 A * 0 = 0。（这很明显）
- 如果您已经实现了 SUM，请检查序列 An = (An-1+An-1) = A1*2n 。我提供了一个算法示例。
- 从 A1 = 5 开始，同时设置 An=5。
- 重复，比如说 7 次 An=An+An，An 的连续值是（5,10,20,40,80,160,320,640），最后，An 的值为 640。
- 然后执行 Temp=A1*27，即 5*128 = 640，最后检查 Temp= An。
- 检查结合律：(A*B)*C=A*(B*C)
- 检查 (A*B MOD m) = (A MOD m) * (B MOD m)，如果您在基数 10 中选择模数 9，这被称为“九余检验”，对于 32 位数，一个简单的选择是 232-1。九余检验（在我们的例子中是 232-1 检验）很容易做到，因为您只需要将所有单个数字相加来计算模数：以基数 10 为例。
53 * 31 = 1643 -> 1+6+4+4 = 14 -> 1+ 4 = 5, so 53*31 MOD 9 = 5. 53 -> 5 + 3 -> 8, so 53 MOD 9 = 8 31 -> 3 + 1 -> 4, so 31 MOD 9 = 4 8 * 4 = 32 -> 3 + 2 = 5, so (53 MOD 9) * (31 MOD 9) = 5
粗体数字的两个模数具有相同的值，到目前为止一切正常。
在源代码中，我提供了计算 A MOD 264-1 的算法，这个测试并非万无一失，但如果它说“您的算法是错误的”，那么该算法绝对是错误的。
- 将实现彼此进行比较，如果其中一个出错，您会注意到它；如果两者都出错，您很有可能注意到，除非它们具有相同的错误。
- 对结果已知的特殊输入数字进行测试。
- 还要对一些手动制作的数字进行测试，以便您确定能触发您代码的某些部分：例如，如果您针对数字为零的情况进行了优化，最好创建一个包含一些零的数字，否则您的分支在随机测试中很少会被执行。
- 我们应该有一些测试来处理实现者可能弄错进位检测算法的可能性，这很难做到，所以我没有这样做，只需极其小心，并在单独的项目中彻底测试进位检测部分，因为这是随机数字测试很难发现的一种错误。
通过让一个例程通过以上所有测试，我们可以对实现的正确性有信心。当然，我们必须对其进行精心设计的数字测试，然后对任何可能长度的随机选择的数字进行一些测试。
通过数学证明正确性
如果您想获得更高的信心，除了添加越来越多的测试之外，您还可以始终通过数学技术来证明您的实现，这超出了我能做到的范围，但大致来说，它是这样的：对于代码的每一行，都要构建一个证明，证明指令执行后程序的が状态是正确的（符合断言），例如。
carry = ((operand + operand2 + carry) < (operand + carry)
要证明该指令的正确性，首先您需要对程序的が状态有一些断言，例如：执行后，如果操作数 + 操作数 2 + 进位溢出，则进位值将为 1，然后您应该证明对于每个可能的输入が状态（结果、操作数、操作数2、进位）（您知道这一点，因为您对所有之前的指令都有证明），输出が状态符合您的断言。
例如，在这种情况下，假设您的计算机工作在基数 10，因此结果、操作数和操作数2是基数 10 的数字，进位只能是 1 或 0。
那么，如果我们有操作数 = 9，操作数2 = 0 且进位 = 1，我们可以说该指令不符合断言。
carry = (9 + 0 + 1) < (9 + 1)
因此
carry = 10 < 10
而且因为我们是模 10 进行计算
carry = 0 < 0
因此
carry = 0
但是我们发生了溢出，因为 9 + 0 + 1 不适合一个基数 10 的数字。
那么我证明了指令不正确吗？这取决于我们对输入が状态的证明，如果我们能证明操作数永远不会是 9，那么指令就是正确的……
这就是想法，事情更深奥，但也许您喜欢这个想法，并会深入研究。
后果
请注意，提出的算法比长加法需要更多的工作。
有三个嵌套循环，其中第二个循环使用乘法，这是一个相对昂贵的操作，但在现代硬件上，它在功耗和晶体管方面比时间更昂贵。
总的来说，乘法比加法昂贵很多倍，但比除法便宜……因为除法需要多次乘法，所以乘法算法的性能将直接影响除法性能。
除法将是该系列的最后一部分。由于算法复杂度为 O(N^2)，我进行的测试是在 4KB*2KB 数字上进行的（每次运行 10K 次迭代），如果我像处理加法那样保持 512KB 数字，它将运行得太久。
上述结果由 VS 2017 C++ 在 2011 年的 i7 处理器上生成。（该算法不是并行的，因此使用 1 核。）
总之，我认为我已经证明了该算法的 32 位版本（黄色条上的版本，对应于 LongMultiplicationNoAssembly func）运行速度比使用原生 64 位 CPU MUL 的版本慢两倍以上，正如理论上应该的那样，能够使用 64 位寄存器的良好 C 版本将会很不错。
我在删除跳过循环的所谓“优化”后重新运行了测试，即当其中一个操作数为零时。我开始考虑一个完整的 64 位操作数为零的概率……嗯，在十进制中，改进一个可能发生 1/10 次的情况是明智的，但在六十四进制中，编写一些代码来处理一个可能发生十亿分之一的机会是一个非常糟糕的主意，上面的灰色汇编过程已经忽略了这种“优化”，并且比其他过程快得多，让我们让另一个算法做同样的事情，以便我们可以进行公平的比较。
这是更新后的统计数据，改进虽然不大但可衡量。
上述结果由 VS 2017 C++ 在 2011 年的 i7 处理器上生成。（该算法不是并行的，因此使用 1 核。）
最后的想法是，我们应该对可以瞄准的最高速度极限有一个更好的想法，以便找到更好的实现。
在我最好的实现中，我每秒每核可以运行 1.8 亿次字乘法，如果我没算错的话，请查看我的计算。
我的 3GHZ CPU 花费了，如果 Windows 性能计数器正确的话，平均每个字乘法 16 个时钟周期，使用这个算法或许可以做得更好。
下面是 VS 2022 C++ 在 i9-10900 处理器上生成的一些新结果。（该算法不是并行的，因此使用 1 个核心）。
历史
- 2017：开始写作
- 2019 年 2 月 6 日：几乎准备发布
- 2019 年 3 月 20 日：添加了新的实现来尝试证明我的推测
- 2020 年 6 月 9 日：添加了关于测试实现正确性的建议
- 2022 年 6 月 11 日：上线