多精度算术： 除法算法

标准猜测

我们这样找到我们的猜测值 Q_i

• 边界情况 A_{i <} B
  • Q_i = 0
• 否则如果 A_n = B_m , 那么
  • 如果 A_i 的位数与 B 相同，那么 Q_i = 1
    • 示例：A_i = 213 : B = 200, Q_i = 1
    • 示例：A_i = 2137 : B = 2122, Q_i = 1
    • _{Ai < B 的情况已经被前一个分支处理了……}
  • 否则（_{当 Ai 比 B 长 1 位时，没有其他可能的情况}） Q_i = base-1 （对于十进制字符串是 9，在 CPU 上是 0xFFFF....F）。
    • 示例：A_i = 2137 : B = 220, Q_i = 9
    • _{情况 Ai = 2237 : B = 220 不可能存在，因为： Ai 是 A 中大于 B 的最短最左边的数字组，在这种情况下 Ai 会是 223，因此这种情况应该由上面的分支处理}
• else
  • Q_i = simple_division(A_guess / B_m)，_{其中如果我们手工计算，简单除法可以是逆向查找；在计算机上，则是 CPU 的 Div 指令。}

注意 A_guess 只能是 1 位或 2 位数字，当它是 2 位时，A_n < B_m：所以我们对 Q_i = A_guess / B_m 的标准猜测可以使用内置的 CPU 指令来计算。

改进的猜测：稍后，我们将使用一些预处理：我们将摆脱一些分支判断，因为我们将预处理 A 和 B ，使得 B_m >= base/2。

距离真实的 Q_I 有多远？

现在我们需要确保我们猜对了 Q_i ：为此，我们必须运行长乘法算法 B * Q_i ，例如 56 * 2 = 112，由于 112 <= 123，我们就算完成了。

等等：我们难道不应该检查我们的猜测是否太小吗？为什么我们不检查 56 * 3 < 123 呢？嗯，那是不可能的，我们的猜测总是正确或太大，绝不会太小。证明在附录 1。

有时我们的猜测太大了

如果我们不对数字进行预处理，我们可能会遇到麻烦。

6789:18

使用我们的标准猜测算法，A_guess 将是 6，6 : 1 = 6，但 6*18=108，这比 67 大得多。

在这种情况下，我们将尝试减少我们的猜测，直到找到正确的结果。我们将我们的 Q_i 猜测值减 1 并尝试：5*18，然后是 4*18，最终我们发现 3*18 < 67，所以 Q_i = 3。

与其将我们的猜测减一然后乘以 B，我们将使用我们的中间结果（108）并持续减去 B（18），直到 Q_i 变得小于 67，因为这将节省少量计算时间：（我的猜想，没有证据支持）。

达到正确的 Q_I 需要多大成本？

1. 108-18=90
2. 90-18=72
3. 72-18=54

我们不得不多做三次减法才找到正确的数字，我们最初的猜测相当错误，正确的数字是 108 的一半，而且那还是十进制。如果使用 2²⁵⁶ 进制，试错过程将变得不切实际。

归一化除法

幸运的是，代数书上说，如果 \(B_m \geq \frac{base}{2}\)，那么我们最多只需要做 2 次额外的减法。我不会为此提供证明，如果有人需要书名，我会推荐我的书。

假设我们没那么幸运，\(B_m <\frac{base}{2}\)：原始除法是

e1: A = B * Q + R

我们可以预处理这个方程，我们可以合法地在两边都乘以 k，目的是使 B_m >= base/2。

e2: kA = kB * Q + kR

比较 e1 和 e2 表达式：Q 没有改变，改变的是 R，我们在返回前需要将 R 除以 k。

在我们的实际算法中，我们进行左移操作，直到最左边数字的最左边的位变为 1，而不是进行昂贵的乘法或除法。

回到例子

6789 / 18

B_m 是 1，所以我们不断乘以 2（在 CPU 上是左移），直到 B_m 大于 5（因为这个例子是十进制，在 CPU 上我们需要让 B_m 的最左边的位变为 1）。

1. 18 << 1 = 36
2. 36 << 1 = 72

我们需要 2 次移位（相当于乘以 4）。

同时也要移位被除数

6789 << 2 = 27156

现在 B_m 是 7，即 >= (10/2)。

原始除法是 6789/18 (Q=377, R=3)，现在除法变成 6789*4/18*4 (Q=377, R=3*4)，可以看作是 27156/72 (Q=377, R=12)。

代数书称之为归一化除法，Q 不变，R 在我们进行新的除法时被乘以了 4，但现在我们确定我们只需要做很少的猜测。

CPU 的“逆向位扫描”（Bit Scan Reverse）等功能在低级语言中可用，可以省去左移循环：我没有证据表明 CPU 的 BSR 指令比位移循环更快。

回顾算法，增加猜测算法并移除限制

现在我们将显式地将 FindQn (猜测) 子程序添加到我们的算法中。如果 B 是归一化的，其内部循环的迭代次数最多为 2，但余数将被移位。由于输入是归一化的，我们也可以简化猜测例程。

修订后的算法处理了 A <= B 的情况。为简单起见，“Find Qn”子程序被呈现为一个函数，但可以内联开发。它不包括归一化过程，让我们假装这是核心算法，而代码的另一部分将通过函数组合来处理归一化。

//TODO：在该图表中，缺少长乘法和测试。

使用伪代码的最终示例。

现在我们将逐步执行修订后的算法：我们将接收归一化的数字作为输入，以便 B 的最左边数字 >= base/2

因此，我们不是除以 6789 / 18，而是将两边都乘以 4，得到一个新的除法 27156 / 72，我们的返回值将是 R / 4 和 Q。（译者注：原文为 R * 4，应为笔误，根据上下文意为处理后的 R，在最后需要除以 4 得到最终余数）

[Q]=empty,

//An = 271 because that is the first group of digits on 27156 that is > 72
[An] = Select as many digits from [A] such that [An] >= [B];

while([An] > [B])

FIRST ITERATION:

loop:
    Qn = 0;

    if(An < B) //271 < 72 so continue
      return zero;

    Bn = Leftmost of B; //7

    Aguess = Select as many digits from [An] such that Aguess >= Bn; //Aguess = 27
    Qn = Aguess / Bn; //Qn = 3

    Test = Qn\*B; //Test = 3*72 = 216

    while(Test > An) //216 > 271 is false so leave subroutine

    {

       Qn--;

       Test -= B;

    }

   return done;

   //end of subroutine:

   Append Qn to Q; //Q become 3

   Remainder = An = An - Test; //remainder become 271 – 216 = 55

   if(there are more digits in A) //there are more digits in A

   {

     Append next digit from A to An; //An become 555

      goto loop; //continue iteration

   }


SECOND ITERATION:

loop:
    Qn = 0;

    if(An < B) //555 < 72 false

      return zero;

    Bn = Leftmost of B; //7

    Aguess = Select as many digits from \[An\] such that Aguess >= Bn; //Aguess = 55
    Qn = Aguess / Bn; //Qn = 7

    Test = Qn\*B; //Test = 7*72 = 504

    while(Test > An) //504 > 555 is false so leave subroutine

    {

       Qn--;

       Test -= B;

    }

   return done;

   //end of subroutine:

   Append Qn to Q; //Q become 37

   Remainder = An = An - Test; //remainder become 555 – 504 = 51

   if(there are more digits in A) //there are more digits in A

   {

     Append next digit from A to An; //An become 516

      goto loop; //continue iteration

   }


THIRD ITERATION:

loop:
    Qn = 0;

    if(An < B) //516 < 72 false

      return zero;

    Bn = Leftmost of B; //7

    Aguess = Select as many digits from \[An\] such that Aguess >= Bn; //Aguess = 51
    Qn = Aguess / Bn; //Qn = 7

    Test = Qn\*B; //Test = 7*72 = 504

    while(Test > An) //504 > 516 is false so leave subroutine

    {

       Qn--;

       Test -= B;

    }

   return done;

   //end of subroutine:

   Append Qn to Q; //Q become 377

   Remainder = An = An - Test; //remainder become 516 – 504 = 12

   if(there are more digits in A) //there are no more digits in A so exit

   {

     Append next digit from A to An;

      goto loop; //continue iteration

   }

检查结果：377*72+12=27156

但我们想计算 6789 / 18，余数必须除以 4；所以检查结果

6789 = 18*377 + (12 / 4) = 18 * 377 + 3，这是正确的。

一个计算机算法

计算机算法可能会根据你决定如何表示一个长数而有所不同。

如果你选择使用无符号整数数组来存储一个大数，这里有一个该算法的 C 语言示例版本：很抱歉我没有直接粘贴在这里，因为它太大了。

让我对上面的代码稍作评论：困难的部分（关于除法的一切都很困难，但只要有耐心就可以解决），是简单除法算法的实现。由于我的编译器不支持 128 位整数，我不得不编写一个非常小的 x64 汇编例程来更好地利用硬件，这使得该函数无法移植到 x32，但由于该代码非常简单，你也可以为你的硬件编写该支持函数。

.data

.code
    ;qword cpu_divide(
    ;     qword  low,         RCX
    ;     qword  high,         RDX
    ;     qword  divisor,             R8
    ;     qword * R,            R9
    ; );            
; RAX, RCX, RDX, R8, R9, R10, R11  VOLATILE
; RAX used to store cpu flags
cpu_divide proc
    mov rax, rcx    
    div r8
    mov qword ptr[r9], rdx
    ret
cpu_divide endp
end

因为以我的 x64 技能编写一个完整的 x64 实现会太耗时，我将只发布带有那段代码的文章；实际上我创建了两个 C 版本，第一个几乎不可读，因为我试图做“基于感觉的优化”，第二个更具可读性。

后果

以下是我使用上述实现得到的一些计时数据。

算法复杂度

有趣的是，除法比乘法慢得多，其速度主要取决于内部乘法的速度，否则该算法看起来是线性时间的，但事实并非如此。所有这些乘法都是用 1 位数乘以某个数来完成的，因此乘法算法总是以线性时间执行，即使使用简单的小学乘法算法也足够了……可能比使用任何其他算法都快。值得编写一个专门针对固定 N*1 大小并尽可能优化的算法，除法的速度取决于它。

所以上图的钟形曲线显示，当被除数大小为 N 且除数大小为 N/2 时，性能最差。

这可以解释为

当用 N 位数除以 1 位数时，你会得到（忽略减法）

• N 次 CPU 除法
• N 次大小为 (1 * 1) 的乘法
  总计 N + N * 1 次操作 => O(N)

当用 N 位数除以 N 位数时，你会得到

• 1 次 CPU 除法
• 1 次大小为 1 * N = N 的乘法
  总计 1 + N => O(N)

当用 N 位数除以 M 位数时 (M < N)，你会得到

• N-M+1 次 CPU 除法
• N-M+1 次大小为 1*M=M 的乘法
  总计 (N-M+1) + (N-M+1) * M => O((N-M)*M)
  因为 M<=N，我们有当 M = N/2 时，(N-M)*M 达到最大值

因此我们可以说，最好情况的复杂度是 \(O \left( N \right)\)，最坏情况是 \(O \left( \left( \dfrac{N}{2} \right) ^2 \right){\implies}O \left( N^2 \right) \)

目前，我不知道是否存在更好的算法，我订购了一些书希望能找到答案。

附录 A

本文提出的算法所猜测的 Q_n 大于或等于正确的 Q_n

由于我们讨论的是部分除法，我们只有两种情况

1. A/B 的除法，其中 A 和 B 具有相同的长度
2. A/B 的除法，其中 A 比 B 多 1 位数，但其最左边的数字 An 小于 B 的最左边的数字 Bn

这里是这两种可能情况的图示

只展示情况 1，情况 2 的证明类似。

案例 1

假设

• H.1 \(A > B\)
• H.2 \(A, B\) 具有相同数量的位数
• H.3 \(A_{n}\) 是 \(A\) 在 \(b\) 进制下的最左边的数字
• H.4 \(B_{n}\) 是 \(B\) 在 \(b\) 进制下的最左边的数字
• H.5 \(n\) 是 \(A\) 和 \(B\) 的位数

记住：

• R.1 \(A_{n} \geq B_{n}\) 根据 \eqref{def1} 和 \eqref{def2}
• R.2 \(A_{n} < b\) 根据 \(b\) 进制下数字的定义
• R.3 \(A_{n}, B_{n}\) 是非零的（最高有效位）

让我们将 \(A\) 和 \(B\) 展开为其位值表示法

我们的第一个猜测 \(Q\) 是简单除法 \(A_{n} / B_{n}\) 的单位数商，使得 \(A_{n}= B_{n} Q + R\)

论点

我们想证明，对于任何满足初始假设 H 的数对 \((A, B)\)，如果我们使用提出的算法计算 \(Q\)，那么 \((Q+1)B>A\)。

证明

我们必须证明：\((Q+1) B > A\)。只考虑前一个不等式的左边部分，并使用位值表示法重写 B。

上面的表示法在方程的右侧显示了长数 B 的字符串表示乘以 (Q+1)，根据定义，(Q+1) 是一个单位数。

这等于（根据加法和乘法的分配律）

让我解释一下前面的表达式，这是一个为了简化符号的简化。我们得到了一个新的位值展开数字串，其中最右边最低位的数字是 \([B_{1} * (Q+1) )\ MOD\ b]\)，我们将其重命名为 \([K_{1}]\)。

在进行内部乘法时，我们有一个余数向左传播，传播在表达式的最高有效部分结束，成为数字 \(C\)，我们不知道它有多大，但我们知道它小于 \(b\)。

因此，(eq5) 的左边部分是一个要附加到长 \(K\) 字符串左边的两位数。乘法修改了所有的数字，我们不知道 \(K\) 的数字值，但重要的是我们有了一个新的数，它有一个新的展开，可以有 n 或 n+1 位。

现在我们得到了两个位值表示法字符串 (eq6) 和 (eq7)

将上面的 A 字符串与 B 字符串进行比较，为了证明我们的论点，我们需要证明表达式 (eq6) > (eq7)，对于 \(\forall (A, B),\ H(A,B) \) 成立。

因为 \(0 \leq C < b\)，我将忽略 \(C\)，以得到一个更简单的表达式

如果你看 (eq6)、(eq7) 和 (eq8)，你会注意到，如果我们证明 \((eq8) > (eq7) \Rightarrow (eq6) > (eq7) \)

因为 (eq8) 更容易证明，我们将证明它。

要证明 (eq8) \(>\) (eq7)，有两种可能的情况

情况 1：如果

那么只需证明

情况2：否则，我们需要证明 (eq8) 的最左边部分严格大于 (eq7) 的最左边部分。在这种情况下，我们可以忽略最低有效位，只需要证明 (simpl1)。

鉴于 (情况 2) 比 (情况 1) 更具限制性，并且 \((情况 2) \Rightarrow (情况 1)\)，只需证明 (情况 2) 就足够了。因此我们可以忽略 A 和 B 的最低有效位（在两个表达式中都移除最低有效位），我们只需证明 \((eq11)\)。

我们可以将 \((eq11)\) 重写为（根据分配律）

去掉一些括号，重新排序：\((eq13)\) 是我们想要证明的

现在比较这个不等式系统

\((eq14)\) 与 \((eq13)\) 相同，\((eq15)\) 是我们已知的，Q 是从除法 \(A_{n}\ /\ B_{n}\) 得到的，所以 \(A_{n} = Q \cdot B_{n} + R\)

但是，根据除法余数的定义，\(R < B_{n}\)，所以 (proof1) 是真的，证毕。

历史

• 2022 年 6 月 11 日：在线发布