任意精度（bignum）整数的乘法算法答案

【问题标题】：Multiplication algorithm for abritrary precision (bignum) integers任意精度（bignum）整数的乘法算法
【发布时间】：2010-05-02 21:33:20
【问题描述】：

我正在为一个家庭作业项目编写一个小型 bignum 库。我要实现 Karatsuba 乘法，但在此之前我想编写一个幼稚的乘法例程。

我正在关注 Paul Zimmerman 编写的名为“现代计算机算术”的指南，该指南是 freely available online。

在第 4 页，有一个名为 BasecaseMultiply 的算法的描述，它执行小学乘法。

我理解第 2、3 步，其中 B^j 是 1、j 次的数字移位。但我不明白第 1 步和第 3 步，我们有 A*b_j。如果 bignum 乘法还没有定义，这个乘法应该如何进行？

这个算法中的“*”操作会不会只是重复加法？

这是我到目前为止写的部分。我已经对它们进行了单元测试，因此它们在大多数情况下似乎是正确的：

我的 bignum 结构如下：

#define BIGNUM_DIGITS 2048
typedef uint32_t u_hw; // halfword
typedef uint64_t u_w; // word

typedef struct {
    unsigned int sign; // 0 or 1
    unsigned int n_digits;
    u_hw digits[BIGNUM_DIGITS];
} bn;

当前可用的例程：

bn *bn_add(bn *a, bn *b); // returns a+b as a newly allocated bn
void bn_lshift(bn *b, int d); // shifts d digits to the left, retains sign
int bn_cmp(bn *a, bn *b); // returns 1 if a>b, 0 if a=b, -1 if a<b

【问题讨论】：

您可以尝试使 bignum 乘法例程递归，因为您需要执行 A*b_j。我不确定它会如何定义，所以它几乎不是一个答案，但它至少是一个想法。我相信有一个更好的解决方案。 :)

标签： c algorithm math bignum multiplication

【解决方案1】：

我前段时间写了一个乘法算法，我在顶部有这个评论。如果你有两个大小相同的数字 x 和 y（相同的 n_digits），那么你可以像这样相乘得到 n，它的数字是数字的两倍。该算法的部分复杂性来自于计算出如果两个输入的 n_digits 不同时哪些位不相乘。

从右边开始，n0 是 x0*y0，你可以省去溢出。现在 n1 是 x1*y0 和 y1*x0 的总和，前一个溢出移动了您的数字大小。如果您在 64 位数学中使用 32 位数字，这意味着 n0 = low32(x0*y0) 并且您携带 high32(x0*y0) 作为溢出。您可以看到，如果您实际使用 32 位数字，则无法将中心列相加而不超过 64 位，因此您可能使用 30 或 31 位数字。

如果每个数字有 30 位，这意味着您可以将两个 8 位数字组合在一起。首先编写此算法以接受两个 n_digits 最多为 8 的小缓冲区，并使用本机数学进行算术。然后再次实现它，采用任意大小的 n_digits 并使用第一个版本以及您的 shift 和 add 方法，一次乘以 8x8 块的数字。

/*
    X*Y = N

                          x0     y3
                            \   /  
                             \ /   
                              X    
                      x1     /|\     y2
                        \   / | \   /  
                         \ /  |  \ /   
                          X   |   X    
                  x2     /|\  |  /|\     y1
                    \   / | \ | / | \   /  
                     \ /  |  \|/  |  \ /   
                      X   |   X   |   X    
              x3     /|\  |  /|\  |  /|\     y0
                \   / | \ | / | \ | / | \   /
                 \ /  |  \|/  |  \|/  |  \ /
                  V   |   X   |   X   |   V
                  |\  |  /|\  |  /|\  |  /|
                  | \ | / | \ | / | \ | / |
                  |  \|/  |  \|/  |  \|/  |
                  |   V   |   X   |   V   |
                  |   |\  |  /|\  |  /|   |
                  |   | \ | / | \ | / |   |
                  |   |  \|/  |  \|/  |   |
                  |   |   V   |   V   |   |
                  |   |   |\  |  /|   |   |
                  |   |   | \ | / |   |   |
                  |   |   |  \|/  |   |   |
                  |   |   |   V   |   |   |
                  |   |   |   |   |   |   |
              n7  n6  n5  n4  n3  n2  n1  n0
*/

【讨论】：

我不确定您所说的超过 64 位和 32 位数字是什么意思？如您所见，我有 32 位数字，并计划使用 64 位字长来管理整个范围。
顺便说一句 - 很酷的图表，乘法真的很方便。
basecase 乘法通过一次走一条对角线并构建结果来避免添加列，因此您可以忽略该问题。如果您将两个 32 位乘法的乘积相加，您将得到一个 65 位的值。你要么必须处理溢出（在 asm 中很容易，在 C 中很麻烦），要么你必须在每个数字中使用更少的位。
顺便说一句，执行 c+=ab 比仅执行 c=ab 更便宜，因为在后者中您必须将 c 初始化为零。
你的意思不是正好相反吗？即“在前者”中，您必须将 c 初始化为零，因为在后者 ('=') 中，您只需完全覆盖 c。

【解决方案2】：

要做A*b_j，你需要做一个大数与一位数的小学乘法。你最终不得不将一堆两位数的产品加在一起：

bn *R = ZERO;
for(int i = 0; i < n; i++) {
  bn S = {0, 2};
  S.digits[0] = a[i] * b_j;
  S.digits[1] = (((u_w)a[i]) * b_j) >> 32;  // order depends on endianness
  bn_lshift(S, i);
  R = bn_add(R, S);
}

当然，这是非常低效的。

【讨论】：