整数除法算法

Question

整数除法算法

c++algorithmperformanceinteger-divisionbigint

23

我在思考一个大整数除法的算法：用余数bigint C 除以 bigint D，其中我们知道C在基数b下的表示，而D是形如 b^k-1 的形式。最好通过示例来展示。让我们尝试用D=999将C=21979182173除以。

我们将数字写成三个数字的集合：21 979 182 173
从左侧开始，我们对连续的集合进行求和（模999）：21 001 183 356
我们在"超过999"的那些集合之前加1：22 001 183 356

确实，21979182173/999=22001183，余数为356。

我已经计算了复杂度，并且如果我没有弄错的话，该算法应该在O(n)内工作，n是C在基数b表示中的位数。我还编写了一个非常粗糙且未经优化的版本（仅针对b=10）的算法，在C ++中进行了测试，与GMP的普通整数除法算法进行了比较，它确实似乎比GMP更好。我在任何我查找的地方都找不到类似的实现，因此我不得不使用一般除法进行测试。

我找到了几篇涉及到很相似的问题的文章，但是没有一篇集中讨论实际实现，特别是在不同于2的基数下。我想这是因为数字的内部存储方式，尽管即使考虑到这一点，上述算法在b=10的情况下也很有用。我也尝试联系其他人，但是无济于事。

因此，我的问题是：是否有一篇文章或一本书或其他东西描述了上述算法，可能讨论了其实现方式？如果没有，那么对我来说，在C/C ++中尝试实现和测试这样的算法是否有意义？这个算法本质上是否有缺陷？

另外，我不是程序员，虽然我在编程方面还可以，但我承认我对计算机的“内部”知识并不了解。因此，请原谅我的无知-这篇文章中可能存在一个或多个非常愚蠢的问题。再次抱歉。

非常感谢！

对评论/答案中提出的观点进行进一步澄清：

非常感谢大家 - 由于我不想用同样的话在所有很棒的答案和建议上进行评论，所以我只想解决你们讨论的一个问题。

我完全意识到，在2^n进制下工作通常是做事情最有效的方式。几乎所有的大整数库都使用2^32或其他类似的进制。然而，如果（仅仅针对这个特定算法！）我们将大整数实现为进制为b的数字数组，那么怎么办呢？当然，这里要求b是“合理的”：b=10，也就是最自然的情况，似乎足够合理。我知道这在内存和时间上都相对低效，考虑到数字的内部存储方式，但是如果我的（基本且可能存在缺陷的）测试正确，我已经能够比GMP的普通除法更快地产生结果，这就给实现这样一个算法提供了意义。

Ninefingers指出，这种情况下我必须使用一种昂贵的取模运算。我希望不需要：只需通过查看old+new+1的位数就可以看出是否越过了999。如果它有4个数字，我们就搞定了。更重要的是，由于old<999且new<=999，所以我们知道如果old+new+1有4个数字（它不能有更多），那么(old+new)%999等于删除(old+new+1)的最左边的数字，我想我们可以便宜地做到这一点。

当然，我并不否认这种算法明显存在的局限性，也不声称它无法改进 - 它只能除以某些类别的数字，而且我们必须预先知道在基数b中被除数的表示形式。然而，对于b=10，后者似乎很自然。

现在，假设我们已经按照上面的方式实现了大整数。假设C=(a_1a_2...a_n)是在基数b下表示的，并且D=b^k-1。该算法（可能还可以更加优化）将像这样进行。希望没有太多的错别字。

如果 k>n，则显然我们已经完成了
在 C 的开头添加一个零（即 a_0=0）（以防万一，比如我们试图用 99 去除 9999）
l=n%k （对于“常规”整数的模运算-不应太费力）
old=(a_0...a_l) （第一组数字，可能少于 k 个数字）
for (i=l+1; i < n; i=i+k) （我们将进行 floor(n/k) 次迭代或更多次）
- new=(a_i...a_(i+k-1))
- new=new+old （这是大整数加法，因此 O(k)）
- aux=new+1 （同样是大整数加法 - O(k)，我对此不太满意）
- 如果 aux 有超过 k 个数字
  - 删除 aux 的第一个数字
  - old=old+1 （再次进行大整数加法）
  - 在开头填充零，使 old 具有应有的位数
  - (a_(i-k)...a_(i-1))=old （如果 i=l+1，则为(a _ 0...a _ l)=old）
  - new=aux
- 在开头填充零，使 new 具有应有的位数
- (a_i...a_(i+k-1)=new
quot=(a_0...a_(n-k+1))
rem=new

好的，感谢与我讨论此事 - 如我所说，如果没有人发现任何致命缺陷，这似乎是一个有趣的“特殊情况”算法可以尝试实现、测试和讨论。如果这是迄今为止不广泛讨论的问题，那就更好了。请告诉我你的想法。对于这篇长文章，我感到抱歉。

还有一些个人评论：

@Ninefingers：我对GMP的工作原理、功能和大数除法算法有一些（非常基础的）了解，所以我能够理解你的很多论点。我也知道GMP高度优化，并且在某种程度上为不同的平台定制自己，因此我肯定不会试图在普遍情况下“战胜”它——这似乎与用尖木棍攻击坦克一样毫无成果。然而，这不是这个算法的想法——它适用于非常特殊的情况（GMP似乎没有涵盖）。另外，您确定常规除法是O(n)吗？我见过的最多只有M(n)。（如果我理解正确），这可以在实践中（例如，Schönhage–Strassen等）达到不到O(n)。如果我没记错的话，弗雷尔(Fürer's)算法几乎纯属理论。

@Avi Berger：这似乎并不完全与“抛弃九”的概念相同，尽管思路类似。然而，前面提到的算法应该始终有效，如果我没记错的话。

- mornik

2

所以您是建议将所有整数存储为BCD，以便使除法更快？将二进制转换为十进制涉及整数除法，对吗？ :-) - Ken

4

有趣的算法，但实际应用可能受限。从技术上讲，选择不同的基数可以让你将其用于任何任意除数，但关键在于首先将它转换为该基数。 - Karl Bielefeldt

@Ninefingers：看起来有点长，不太像是编辑。我试图删除它，但似乎已经有管理员删除了它。无论如何，还是谢谢。 - mornik

关于使用基数b的位操作而不总是使用2^32: 如果你需要经常除以某个b，这是一个有效的选择。例如，在一个代码高尔夫挑战中打印Fibonacci(10^9)的前1000个数字（有性能要求），我使用了半蛮力方法，通过在数字过大时除以10^9来保留最重要的1009个十进制数字。基数为10^9的位（32位元素）使得这种方法非常高效，并且手动进行进位比加法快。[105字节的x86机器码] (https://codegolf.stackexchange.com/a/135618/30206) - Peter Cordes

将一个数除以2^n-1是非常简单的。这与除以1-2^(-n)相同，因为 1/(1-2^-n ) = 1+2^-n+2^-2n.....所以只需将n左移并重复加即可。 - Paul

3个回答

5

根据我的评论，我觉得有必要补充一下。这不是一个答案，而是背景的解释。

大数库使用所谓的limbs——在gmp源代码中搜索mp_limb_t，它们通常是固定大小的整数字段。

当您执行类似于加法的操作时，一种（虽然低效的）方法是这样做：

doublelimb r = limb_a + limb_b + carryfrompreviousiteration

这个双倍大小的limb在limb_a + limb_b的和大于limb大小的情况下捕获溢出。因此，如果总数大于2^32（如果我们使用uint32_t作为limb大小），则可以捕获溢出。

为什么需要这个？好吧，你通常会循环遍历所有limbs - 你自己在将整数分解并逐个遍历时已经这样做了 - 但我们首先执行LSL（因此最小的limb首先），就像你手算一样。

这可能看起来效率低下，但这只是C语言处理事物的方式。要真正发挥它的威力，x86有一个名为adc的指令 - 带进位加法。它对您的字段进行算术运算，并在算术超出寄存器大小时设置进位标志。下次进行add或adc操作时，处理器也会考虑进位标志。在减法中，它称为借位标志。

这也适用于移位操作。因此，处理器的这个特性对于使大数快速计算至关重要。因此，事实是，芯片中有电子电路来处理这些内容 - 在软件中处理始终会更慢。

不详细解释，操作是通过添加、移位、减去等能力构建的。它们至关重要。如果你做得对，每个limb都要使用处理器寄存器的全宽度。

第二点 - 进制之间的转换。你不能取一个数字中间的值并改变它的进制，因为你无法考虑到原始进制下其下面的数字溢出，而且那个数字也无法考虑到其下面的数字溢出...等等。简而言之，每当你想改变进制时，你需要将整个大数从原始进制转换为新进制，然后再转回来。所以你必须至少遍历大数（所有limbs）三次。或者，作为替代方案，在所有其他操作中昂贵地检测溢出...记住，现在你需要进行模运算来计算是否溢出，而以前处理器是为我们做这件事的。

我还想补充一点，尽管你现在的方法对于这个案例可能很快，但请记住，作为一个大数库，gmp会为你做很多工作，如内存管理。如果你使用mpz_，那么你正在使用我在这里所描述的抽象层之上，作为起点。最后，gmp针对几乎所有你听说过的平台以及更多的平台使用手动优化的汇编和展开循环。它能与Mathematica、Maple等软件一起发布，这是有很好的理由的。

现在，只是为了参考一些阅读材料。

Modern Computer Arithmetic是一个针对任意精度库的类似Knuth的作品。
Donald Knuth，《半数值算法》（计算机程序设计艺术第二卷）。
William Hart's 博客介绍了实现bsdnt算法的各种除法算法。如果你对大数库感兴趣，这是一个很好的资源。在我开始关注这些东西之前，我认为自己是一个好程序员...

总的来说，除法汇编指令难以理解，因此人们通常计算倒数然后进行乘法运算，就像在模运算中定义除法时所做的那样。现有的各种技术（参见MCA）大多是O(n)的。

编辑：好的，并非所有技术都是O(n)。大多数称为div1（除以不大于一个limb的东西）的技术都是O(n)。当你增加规模时，复杂度会变成O(n^2)；这很难避免。

现在，您能将bigint实现为数字数组吗？当然可以。但是，请考虑加法下面的想法。

/* 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)之间的差异，但我们确实在谈论每个limb的昂贵的除法指令或者每次想要除法时的昂贵的转换。请参见这里。

接下来，您如何扩展通用情况除法？我的意思是，当您想要除以上代码中的两个数字x和y时。答案是您不能，除非使用昂贵的大数基础设施。请参阅Knuth。对大于您大小的数字取模是行不通的。

让我来解释一下。试试21979182173 mod 1099。为了简单起见，我们假设我们可以拥有的最大字段大小是三位数。这只是一个人为的例子，但我所知道的最大字段大小使用gcc扩展名使用128位。无论如何，重点是，你：

21 979 182 173

将您的数字分成几个部分。然后取模并求和：

21 1000 1182 1355

它不起作用。这里Avi是正确的，因为这是一种排除九法或其改编形式，但是由于我们的字段已溢出 - 你正在使用模运算来确保每个字段都在其枝/字段大小内，所以这不起作用。

那么解决方案是什么？将您的数字分成一系列适当大小的大数？然后开始使用大数函数来计算您需要的一切？这比任何现有的直接操作字段的方式都要慢得多。

现在也许您只是提出了除以一个枝而不是大数的情况，在这种情况下，它可以工作，但Hensel除法和预先计算的倒数等也可以做到没有转换要求。我不知道这个算法是否比Hensel除法更快；这将是一个有趣的比较；问题出现在整个大数库中的常见表示上。现有的大数库中选择的表示方式是我已经扩展的原因 - 它在首次执行时在汇编级别上是有意义的。
作为一条附言，您不必使用uint32_t来表示您的数据。您可以使用理想情况下与系统寄存器相同大小的变量（例如uint64_t），以便您可以利用汇编优化版本。因此，在64位系统上，只有在结果超过2^64位时，adc rax，rbx才设置溢出(CF)。 tl;dr版本：问题不在于您的算法或想法；而是在于转换进制时的问题，因为您算法所需的表示方式并非加/减/乘等操作中最高效的方法。用knuth的话来说：“这显示了数学优雅和计算效率之间的差异。”

- user257111

谢谢，这是一个非常好的评论/答案！然而，在我的算法中有一些关于基本转换的问题我想要澄清，但我现在没有时间，所以我可能会在几个小时后回到家写点什么。 - mornik

@mornik 啊，好的。我想是这样，因为它变得非常困难，但我只是想为完整性而说。那么，问题就变成了基本转换和您的表示方式 - 它每个除法都添加了一些费用，而其他任何方法都不存在，即使本身纯除法更快。 - user257111

好的。听到这个消息很高兴。现在，关于您的加法代码-那基本上是我所想的。然而，既然（t/somebase）<2，我们不必使用显式除法或模数：m=t-somebase。如果m<0，则令m=0。我们完成了。当然，那样我们不能使用uint32_t，而要用int32_t。这会有所帮助还是带来更多麻烦呢？（此外，无论如何，请注意在for循环的每一步中我只有三次加法，其中两个几乎总是微不足道的：我只加1，因此，在除非非常罕见的情况下，加法将在第一个，也许是两个或三个limbs之后停止。） - mornik

还要注意的是，在实际使用中，某些基数通常小于100。（很难想象一个人在对用13211902为基数表示的数字进行除法感兴趣时，能够轻松地理解它。）同样，最有用的情况可能是某些基数等于10（也许加上一些其他小的进制数）。 - mornik

要在C中检测进位，使用unsigned int a, b，请使用unsigned sum = a + b; unsigned carry = sum < a;。一个好的编译器将使用由add指令设置的进位标志。当尝试使用该进位进行下一次加法时，gcc和clang表现不佳，例如sum = a + b + carry。clang有时会实际发出一个adc，但为了获得进位，它实际上会进行比较。 - Peter Cordes

显示剩余2条评论

0

如果你需要经常除以相同的除数，使用它（或它的幂）作为基数可以使除法像对于二进制整数进行位移一样廉价。

你可以使用基数999，没有什么特别的，使用以10为底的幂使得转换为十进制整数非常便宜。（你可以逐个四肢地工作而不必在整个整数上执行完整的除法。这就像将二进制整数转换为十进制与将每4位转换成一个十六进制数字之间的区别。）二进制 - > 十六进制可以从最高有效位开始，但是转换为非2的幂次的基数必须使用除法从最低有效位开始。

举个例子，为了回答一个Code Golf问题并且需要满足性能要求，计算斐波那契数列的第1000位小数（10⁹），我的105字节x86机器码答案使用了与这个Python答案相同的算法：通常的a+=b; b+=a斐波那契迭代方法，但每当a变得太大时，就会除以（幂次方）10。

斐波那契数列的增长速度比进位传播快，所以偶尔丢弃低位小数不会改变高位数字的长期值。（您可以保留一些超出所需精度的额外数字）。

除以2的幂次方不适用于此算法，除非您跟踪已经丢弃了多少次幂次方，因为最终二进制->十进制转换将取决于该数量。

因此，对于该算法，您必须执行扩展精度加法，并进行除以10（或任何幂次方的10）。

我将十进制的10⁹位存储在32位整数元素中。通过除以10⁹来实现轻松便捷：只需指针增量跳过低位即可。我没有实际执行memmove，而是仅偏移下一个加法迭代所使用的指针。

我认为除以10的幂次方（而不是10^9）会比较便宜，但需要对每个位进行实际除法，并将余数传递到下一个位。

这种方式的扩展精度加法比二进制位稍微昂贵一些，因为我必须手动使用比较生成进位： sum[i] = a[i] + b[i]; carry = sum < a;（无符号比较）。并且还要根据该比较手动包装到10的9次方上，使用条件移动指令。但我能够将该进位作为输入用于adc（x86带进位加法指令）。

您不需要完整的模运算来处理加法包装，因为您知道最多只包装了一次。

这会浪费每个32位limb的2个bit：10^9而不是2^32 = 4.29... * 10^9。每个字节存储一个十进制数字将显着降低空间效率，并且对性能非常不利，因为在现代64位CPU上，8位二进制加法的成本与64位二进制加法相同。

我旨在实现代码大小：对于纯性能，我将使用64位limbs来保存基于10^19的“digits”。(2^64 = 1.84... * 10^19，因此每64位浪费不到1个bit。)这样可以让您在每个硬件add指令中完成两倍的工作量。嗯，实际上这可能是个问题：两个limbs的总和可能会超过64位整数，因此仅检查> 10^19不再足够。您可以使用基于5*10^18或基于10^18的进制，或者进行更复杂的进位检测，检查二进制进位以及手动进位。

使用每个4位半字节存储打包的BCD将更加影响性能，因为在一个字节内没有硬件支持从一个半字节到下一个半字节的进位阻塞。

总的来说，我的版本在相同硬件上比Python扩展精度版本快了约10倍（但它仍有很大的优化空间，通过减少除法次数可以提高速度）。（70秒或80秒对比12分钟）

尽管如此，我认为对于这个特定的算法实现（我只需要加法和除法，并且每隔几次加法后进行一次除法），选择基于10^9位的limbs非常好。对于第N个斐波那契数，有更高效的算法，不需要进行10亿次扩展精度加法。

- Peter Cordes

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Avi Berger · Accepted Answer

你的算法是基于一个称为“去掉九”的十进制算法的变形。你的例子使用的是基数1000和“去掉”999（比基数少1）的算法。这曾经被教在小学里，作为手算的快速检查方法。我有一位高中数学老师惊讶地发现它不再被教授，所以向我们介绍了它。

在基数1000中，“去掉”999并不能作为通用除法算法。它会生成与实际商和余数同余模999的值，而不是实际值。你的算法有点不同，我还没有验证过它是否可行，但它的基础是有效地使用基数1000和除数比基数少1。如果你想尝试用它来除以47，你必须先将其转换为一个基数为48的数字系统。

搜索“去掉九”以获取更多信息。

编辑：我最初阅读你的帖子有些匆忙，你确实知道这是一个可行的算法。如@Ninefingers和@Karl Bielefeldt在他们的评论中更清楚地表述的那样，你没有包括转换为适合特定除数的基数在内的性能估计。