整数除法算法

Question

我在考虑一个大数除法的算法：用余数除以 bigint C 除以 bigint D，我们知道 C 在基数 b 中的表示，而 D 的形式为 b^k-1。在示例中显示它可能是最容易的。让我们尝试将 C=21979182173 除以 D=999。

• 我们将数字写成三位数：21 979 182 173
• 我们对连续集合求和（模 999），从左开始：21 001 183 356
• 我们在“超过 999”之前的那些集合中添加 1：22 001 183 356

确实，21979182173/999=22001183 和其余 356。

我已经计算了复杂度，如果我没记错的话，算法应该在 O(n) 中工作，n 是基 b 表示中 C 的位数。我还在 C++ 中做了一个非常粗略和未优化的算法版本（仅适用于 b=10），针对 GMP 的通用整数除法算法对其进行了测试，它确实似乎比 GMP 更好。我在任何地方都找不到类似的实现，所以我不得不求助于对一般除法进行测试。

我发现几篇文章讨论了似乎非常相似的问题，但没有一篇文章专注于实际实现，尤其是在不同于 2 的基数中。我想这是因为数字在内部存储的方式，尽管提到的算法似乎对例如 b=10 很有用，即使考虑到这一点。我也尝试联系其他人，但同样无济于事。

因此，我的问题是：是否有文章或书籍或其他东西描述了上述算法，可能讨论了实现？如果不是，我尝试在 C/C++ 中实现和测试这样的算法是否有意义，或者这个算法在某种程度上是天生不好的？

另外，我不是程序员，虽然我在编程方面还算可以，但我承认我对计算机“内部结构”知之甚少。因此，请原谅我的无知——这篇文章中很可能有一件或多件非常愚蠢的事情。再次抱歉。

非常感谢！

---

进一步澄清 cmets/answers 中提出的观点：

谢谢大家 - 因为我不想用同样的东西评论所有很棒的答案和建议，所以我只想谈谈你们很多人提到的一点。

我完全清楚，一般来说，以 2^n 为基数显然是最有效的做事方式。几乎所有 bigint 库都使用 2^32 或其他。但是，如果（而且，我强调，它只对这个特定的算法有用！）我们将 bigints 实现为以 b 为底的数字数组，该怎么办？当然，我们在这里要求 b 是“合理的”：b=10，最自然的情况，似乎足够合理。我知道考虑到内存和时间，考虑到数字是如何在内部存储的，我知道它或多或少是低效的，但我已经能够，如果我的（基本的和可能有某种缺陷的）测试是正确的，比 GMP 的一般部门更快地产生结果，这将有助于实现这样的算法。

Ninefingers 注意到在这种情况下我必须使用昂贵的模运算。我希望不会：我可以通过查看 old+new+1 的位数来查看 old+new 是否交叉，例如 999。如果它有 4 位数字，我们就完成了。更重要的是，由于old

当然，我不是在争论这个算法的明显局限性，也不是我声称它不能改进——它只能除以特定类别的数字，我们必须先验地知道以 b 为基数的被除数的表示.但是，例如，对于 b=10，后者似乎很自然。

现在，假设我们已经实现了我上面概述的 bignums。说以 b 为底的 C=(a_1a_2...a_n) 和 D=b^k-1。该算法（可能会更加优化）会像这样。希望不要有太多错别字。

• 如果 k>n，我们显然已经完成了
• 在 C 的开头添加一个零（即 a_0=0）（以防我们尝试将 9999 除以 99）
• l=n%k （“常规”整数的 mod - 不应该太贵）
• old=(a_0...a_l) （第一组数字，可能少于k个数字）
• for (i=l+1; i （我们将有 floor(n/k) 次左右的迭代）
  • new=(a_i...a_(i+k-1))
  • new=new+old （这是 bigint 加法，因此 O(k)）
  • aux=new+1 （再次，bigint 加法 - O(k) - 我不满意）
  • 如果 aux 有超过 k 位
    • 删除辅助的第一个数字
    • old=old+1 （再次添加 bigint）
    • old 在开头用零填充，这样它就有尽可能多的数字
    • (a_(ik)...a_(i-1))=old (如果 i=l+1, (a _ 0...a _ l)=旧)
    • 新=辅助
  • 在新的开头用零填充，这样它就有尽可能多的数字
  • (a_i...a_(i+k-1)=新
• quot=(a_0...a_(n-k+1))
• rem=新

感谢您与我讨论这个问题 - 正如我所说，如果没有人发现它有任何致命缺陷，我认为这确实是一个有趣的“特殊情况”算法，可以尝试实施、测试和讨论。如果这是迄今为止尚未广泛讨论的事情，那就更好了。请让我知道你在想什么。抱歉，帖子太长了。

另外，还有几个个人的 cmets：

@Ninefingers：我实际上对 GMP 的工作原理、它的作用以及一般的 bigint 除法算法有一些（非常基本的！）知识，所以我能够理解你的大部分论点。我也知道 GMP 是高度优化的，并且在某种程度上为不同的平台定制了自己，所以我当然不会试图“打败它”——这似乎与用尖棍攻击坦克一样富有成效。然而，这不是这个算法的想法——它适用于非常特殊的情况（GMP 似乎没有涵盖）。在不相关的说明中，您确定在 O(n) 中完成了一般划分吗？我见过的最多的是M（n）。 （如果我理解正确，那在实践中（Schönhage–Strassen 等）可能不会达到 O(n)。Fürer 的算法仍然没有达到 O(n)，如果我是正确的，几乎纯粹是理论。）

@Avi Berger：这实际上似乎与“抛出九”并不完全相同，尽管想法相似。但是，如果我没记错的话，上述算法应该一直有效。

Answer 1

您的算法是一种以 10 为基数的算法的变体，称为“投出 9”。您的示例使用基数 1000 并“淘汰” 999（比基数少一个）。这曾经在小学教过，作为快速检查手算的方法。我有一位高中数学老师，得知它不再被教了，吓坏了，并填补了我们的空缺。

在基数 1000 中排除 999 不能作为通用除法算法。它将生成与实际商和余数模 999 一致的值 - 而不是实际值。你的算法有点不同，我没有检查它是否有效，但它基于有效地使用基数 1000 并且除数比基数小 1。如果您想尝试除以 47，则必须先转换为以 48 为基数的数字系统。

谷歌“投出九”以获取更多信息。

编辑：我最初阅读您的帖子有点太快了，您确实知道这是一种有效的算法。正如@Ninefingers 和@Karl Bielefeldt 在他们的 cmets 中所说的比我更清楚，您在性能估计中没有包括转换为适合手头特定除数的基数。

Answer 2

我觉得有必要根据我的评论对此进行补充。这不是答案，而是对背景的解释。

bignum 库使用所谓的肢体 - 在 gmp 源中搜索 mp_limb_t，它通常是一个固定大小的整数字段。

当你做加法之类的事情时，一种方法（尽管效率低下）是这样做：

doublelimb r = limb_a + limb_b + carryfrompreviousiteration

在总和大于肢体大小的情况下，这个双倍大小的肢体会捕获肢体_a + 肢体_b 的溢出。因此，如果总数大于 2^32，如果我们使用 uint32_t 作为肢体大小，则可以捕获溢出。

为什么我们需要这个？好吧，你通常做的是循环遍历所有的肢体——你已经自己完成了将整数除以并遍历每个肢体——但我们首先做 LSL（所以首先是最小的肢体），就像你做算术一样手工。

这可能看起来效率低下，但这只是 C 的做事方式。要真正打破大炮，x86 有 adc 作为指令 - 添加进位。这是一个算术和你的字段，如果算术溢出寄存器的大小，则设置进位位。下次您执行add 或adc 时，处理器也会考虑进位位。在减法中，它被称为借用标志。

这也适用于班次操作。因此，处理器的这一特性对于使 bignums 快速运行至关重要。所以事实是，芯片中有电子电路来做这些事情 - 在软件中做这件事总是会更慢。

无需过多详细说明，操作就是通过这种加、移、减等能力建立起来的。它们至关重要。哦，如果你做得对的话，你会使用处理器寄存器的整个宽度。

第二点——碱基之间的转换。您不能在数字中间取一个值并更改它的基数，因为您无法解释原始基数中其下方数字的溢出，并且该数字无法解释从下面的数字溢出......等等。简而言之，每次你想改变基数时，你都需要将整个 bignum 从原来的基数重新转换回你的新基数。所以你必须至少走三遍（所有的四肢）。或者，在所有其他操作中昂贵地检测溢出......请记住，现在您需要进行模运算来计算是否溢出，而在处理器为我们执行此操作之前。

我还想补充一点，虽然对于这种情况，您所拥有的可能很快，但请记住，作为 bignum 库，gmp 为您做了很多工作，例如内存管理。如果您使用的是mpz_，那么对于初学者，您使用的是我在此处描述的抽象之上。最后，gmp 对您听说过的几乎所有平台以及更多平台都使用了带有展开循环的手动优化组装。它与 Mathematica、Maple 等人一起提供是有充分理由的。

现在，一些阅读材料仅供参考。

• Modern Computer Arithmetic 是针对任意精度库的类似 Knuth 的作品。
• Donald Knuth，半数值算法（计算机编程艺术第二卷）。
• William Hart's blog 为 bsdnt 实现算法，他在其中讨论了各种除法算法。如果您对 bignum 库感兴趣，这是一个很好的资源。在我开始关注这类东西之前，我认为自己是一名优秀的程序员......

总结一下：除法汇编指令很烂，所以人们通常计算逆和乘法，就像在模算术中定义除法时一样。现有的各种技术（参见 MCA）大多是 O(n)。

---

编辑：好的，并非所有技术都是 O(n)。大多数称为 div1 的技术（除以不大于肢体的东西是 O(n)。当你变大时，最终会得到 O(n^2) 复杂度；这很难避免。

现在，您可以将 bigints 实现为数字数组吗？嗯，是的，你当然可以。但是，考虑一下加法下的想法

/* you wouldn't do this just before add, it's just to 
   show you the declaration.
 */
uint32_t* x = malloc(num_limbs*sizeof(uint32_t));
uint32_t* y = malloc(num_limbs*sizeof(uint32_t));
uint32_t* a = malloc(num_limbs*sizeof(uint32_t));
uint32_t m;

for ( i = 0; i < num_limbs; i++ )
{
    m = 0;
    uint64_t t = x[i] + y[i] + m;
    /* now we need to work out if that overflowed at all */
    if ( (t/somebase) >= 1 ) /* expensive division */
    {
        m = t % somebase; /* get the overflow */
    }
}

/* frees somewhere */

这是您希望通过您的方案添加的内容的粗略草图。所以你必须在碱基之间运行转换。因此，您需要将基础转换为到您的表示，然后在完成后返回，因为这种形式在其他任何地方都非常慢。我们这里不是在讨论 O(n) 和 O(n^2) 之间的区别，而是在讨论一个昂贵的除法指令perlimb 或想分。 See this.

接下来，如何将您的部门扩展到一般案件部门？我的意思是当你想从上面的代码中除以这两个数字 x 和 y 时。答案是，如果不求助于昂贵的基于 bignum 的设施，你就无法做到。见克努特。取模数大于您的大小是行不通的。

让我解释一下。试试 21979182173 mod 1099。为了简单起见，我们在这里假设 我们可以拥有的最大尺寸字段是三位数。这是一个人为的例子，但我知道的最大字段大小是否使用 128 位使用 gcc 扩展。不管怎样，重点是，你：

21 979 182 173

把你的号码分成四部分。然后你取模和求和：

21 1000 1182 1355

它不起作用。这是 Avi 正确的地方，因为这是一种排除 9 的形式，或者它的改编形式，但它在这里不起作用，因为我们的字段一开始就溢出 - 您正在使用模数来确保每个字段都保持在范围内它的肢体/字段大小。

那么解决办法是什么？把你的号码分成一系列大小合适的大数字？并开始使用 bignum 函数来计算您需要的一切？这将比任何现有的直接操作字段的方式慢得多。

现在，也许您只是提出用肢体除法而不是大数来划分这种情况，在这种情况下它可以工作，但是亨塞尔除法和预先计算的逆等没有转换要求。我不知道这个算法是否会比亨塞尔除法更快；这将是一个有趣的比较； bignum 库中的一个通用表示带来了问题。在现有的 bignum 库中选择的表示是出于我已经扩展的原因 - 它在汇编级别是有意义的，它是第一次完成的。

作为旁注；您不必使用uint32_t 来代表您的四肢。您最好使用系统寄存器的大小（例如 uint64_t），以便您可以利用汇编优化版本。所以在 64 位系统上adc rax, rbx 仅在结果溢出 2^64 位时才设置溢出 (CF)。

tl;dr 版本：问题不在于您的算法或想法；这是在基数之间转换的问题，因为您的算法所需的表示并不是在 add/sub/mul 等中最有效的方法。套用 knuth 的话说：这向您展示了数学优雅和计算效率之间的区别。

Answer 3

如果您需要经常除以同一个除数，使用 it（或它的幂）作为基数可以使除法与位移位对于基数为 2 的二进制整数一样便宜。 p>

如果你愿意，你可以使用基数 999；使用 10 次方基数没有什么特别之处，只是它使转换为十进制整数非常便宜。 （您可以一次工作一个肢体，而不必对整个整数进行完全除法。这就像将二进制整数转换为十进制与将每 4 位转换为十六进制数字之间的区别。二进制 -> 十六进制可以开始具有最高有效位，但转换为非 2 次方基数必须使用除法以 LSB 优先。）

---

例如，要计算具有性能要求的代码高尔夫问题的斐波那契 (10⁹) 的前 1000 个十进制数字，my 105 bytes of x86 machine code answer 使用与 this Python answer 相同的算法：通常 a+=b; b+=a 斐波那契迭代，但每次 a 变得太大时除以 10。

斐波那契的增长速度比进位传播的速度快，因此偶尔丢弃低位小数位不会长期改变高位位。 （你保留一些超出你想要的精度的额外部分）。

除以 2 的幂 不起作用，除非您跟踪您丢弃了多少 2 的幂，因为最后的二进制 -> 十进制转换将取决于那个。

所以对于这个算法，你必须进行扩展精度加法，然后除以 10（或任何你想要的 10 的幂）。

---

我将 base-10⁹ 肢体存储在 32 位整数元素中。除以 10⁹ 非常便宜：只是一个指针增量来跳过下肢。我没有实际执行memmove，而是偏移了下一次添加迭代使用的指针。

我认为除 10^9 以外的 10 次幂会有点便宜，但需要在每个分支上进行实际除法，并将余数传播到下一个分支。

以这种方式扩展精度加法比使用二进制肢体要贵一些，因为我必须通过比较手动生成进位：sum[i] = a[i] + b[i];carry = sum < a;（无符号比较）。并且还使用条件移动指令根据该比较手动换行到 10^9。但我能够使用该进位作为adc 的输入（x86 add-with-carry 指令）。

你不需要一个完整的模来处理加法的包装，因为你知道你最多包装一次。

这会浪费每个 32 位肢体的 2 位多一点：10^9 而不是 2^32 = 4.29... * 10^9。每个字节存储一个以 10 为基的数字会显着降低空间效率，并且非常性能更差，因为 8 位二进制加法的成本与现代 64 上的 64 位二进制加法相同- 位 CPU。

我的目标是代码大小：为了纯粹的性能，我会使用 64 位肢体，持有 base-10^19“数字”。 （2^64 = 1.84... * 10^19，因此每 64 位浪费了不到 1 位。）这使您可以使用每个硬件 add 指令完成两倍的工作。嗯，实际上这可能是个问题：两个分支的总和可能包含 64 位整数，因此仅检查 > 10^19 已经不够了。您可以在 base 5*10^18 或 base 10^18 中工作，或者进行更复杂的进位检测，检查二进制进位和手动进位。

每 4 位半字节存储一个数字的压缩 BCD 对性能来说会更糟，因为没有硬件支持阻止一个字节内从一个半字节到下一个半字节的进位。

---

总体而言，我的版本在相同硬件上的运行速度比 Python 扩展精度版本快了大约 10 倍（但它有显着优化速度的空间，通过减少划分频率）。 （70 秒或 80 秒对比 12 分钟）

不过，我认为对于 that 算法的这个特定实现（我只需要加法和除法，并且每隔几次加法就会发生除法），base-10^9 肢体的选择非常好的。对于第 N 个斐波那契数，有更高效的算法，不需要进行 10 亿次扩展精度加法。