我应该使用哪种算法来实现高性能的大整数除法？

Question

我应该使用哪种算法来实现高性能的大整数除法？

8

我正在将大整数编码成一个size_t数组。我已经完成了其他操作（例如加、减、乘），以及除以单个数字。但如果可能的话，我希望与我的乘法算法相匹配的时间复杂度。目前使用Toom-Cook算法。

据我所知，有线性时间算法可以获取被除数的各种乘法逆元。这意味着，从理论上讲，我可以在与乘法相同的时间复杂度下进行除法，因为相比之下，线性时间操作“微不足道”。

我的问题是，我应该如何做到这一点？实践中最好使用哪种类型的乘法逆元？模64^digitcount？当我将乘法逆元乘以我的除数时，由于整数截断，我能否缩小计算数据的部分而不必计算即将被丢弃的部分？是否有专门的除法算法比逆元方法更好？

编辑：我找到了上面提到“逆元”方法的来源。在《计算机程序设计艺术》第2卷《半数值算法》的第312页上，Knuth提供了“算法R”，这是一种高精度的倒数。他说它的时间复杂度小于乘法。然而，将其转换为C并测试它并不容易，并且不清楚会消耗多少开销的内存等等，直到我编写这个过程需要一段时间。如果没有人赶上我，那么我会发布它。

- VoidStar

你知道那些方法的渐近复杂度吗？以传入函数的数字计数为单位？与台式机乘法的O(n^2)进行比较等。 - VoidStar

“O(n*log(n))”听起来太快了，比最快的乘法还要快。我怀疑由于某些原因它可能会稍微慢一些，但如果我能找出原因，我会再回复你的。 - VoidStar

将注释移到答案中，添加了一些信息的二进制长除法示例... - Spektre

2个回答

3

我不知道乘法逆算法，但这听起来像是Montgomery Reduction或巴雷特补偿的修改。

我使用了不同的大整数除法方法。

请参见bignum division。特别是查看近似除法器和其中的2个链接。其中一个是我的定点除法器，另一个是快速乘法算法（如Karatsuba、NTT上的Schönhage-Strassen）和我的非常快速的32位基数NTT实现的度量，并带有一个链接。

我不确定逆乘因子是否是正确的方法。

它主要用于取模操作，其中除数是固定的。恐怕对于任意除法，获取bigint逆所需的时间和操作可能比标准除法本身更大，但由于我不熟悉它，我可能错了。

在实现中使用最常见的除法器是Newton-Raphson除法器，它与上面链接中的近似除法器非常相似。

近似/迭代除法器通常使用定义其速度的乘法。

如果数字足够小，则通常是长二进制除法和32/64位数码基数除法足够快：通常它们具有较小的开销，并且让n成为处理的最大值（不是数字的位数！）

二进制除法示例：

是O(log32(n).log2(n)) = O(log^2(n))。
它循环遍历所有重要的位。在每次迭代中，您需要执行compare, sub, add, bitshift。这些操作中的每一个都可以在log32(n)中完成，而log2(n)是位数。

这里是我一个bigint模板（C++）的二进制除法示例：

template <DWORD N> void uint<N>::div(uint &c,uint &d,uint a,uint b)
    {
    int i,j,sh;
    sh=0; c=DWORD(0); d=1;
    sh=a.bits()-b.bits();
    if (sh<0) sh=0; else { b<<=sh; d<<=sh; }
    for (;;)
        {
        j=geq(a,b);
        if (j)
            {
            c+=d;
            sub(a,a,b);
            if (j==2) break;
            }
        if (!sh) break;
        b>>=1; d>>=1; sh--;
        }
    d=a;
    }

N是用于存储大整数的32位DWORD数的数量。

c = a / b
d = a % b
qeq(a,b)是一个比较：如果a>=b，则为大于或等于（在log32(n)=N中完成）
它返回a<b时的0，返回a>b时的1，返回a == b时的2
sub(c,a,b) 是 c = a - b

加速的提升源于不使用乘法（如果不计算位移的话）

如果您使用像2^32（ALU块）这样的大基底数字，则可以使用32位内置ALU操作将整个重写为类似于多项式的样式。这通常甚至比二进制长除法更快，其思想是将每个DWORD处理为单个数字，或者将已使用的算术逐渐减半直到达到CPU能力。
请参阅通过半位宽算术进行除法运算

当与大整数计算时

如果您优化了基本操作，则复杂度甚至可以进一步降低，因为随着迭代，子结果会变得更小（从而改变基本操作的复杂度）。 NTT基础乘法的一个很好的例子。

开销可能会搞砸事情。

由于这个原因，运行时间有时不会复制大O复杂度，因此您应该始终测量阈值并针对所使用的位数使用更快的方法来获得最大性能并进行优化。

- Spektre

在大O符号中，您应该始终剥离标量常数。O(log32(n))=O(log(N))，因为它们与描述增长率无关。其次，大O符号最有用且最常用的表述方式是基于输入位数而非可处理值的大小。因此，您应该以数字计数为基础进行计算，而不是值的大小。您展示的是一个O(n^2)算法，这是可以通过Knuth的高速倒数结合快速乘法来加速的（对于极大的输入，您的算法适用于中等规模的数据）。 - VoidStar

@VoidStar 在这种情况下，二进制长除法的结果是 O(n^2)。 - Spektre

1

@VoidStar 出于好奇，您所说的“非常大”和“中等大小”是指多少位数字？ - Fabio says Reinstate Monica

@FabioTurati那取决于实现的方式...例如看一下fast bignum sqr, 我基于NTT的Sqr阈值是操作数的310*32=9920位(结果为19840位)，而NTT mul则有1396*32=44672位的结果，这些数字非常巨大...当你改变实现（优化或其他），阈值也会发生改变，同样适用于计算平台的变化。 - Spektre

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

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

- user1084944 · Accepted Answer

GMP库通常是好算法的良好参考。他们在除法的文档算法中，主要依赖于选择非常大的基数，这样您就可以通过长除法将4位数字除以2位数字并进行处理。长除法需要计算1位乘以2位商数；这可以通过递归或预先计算逆和估计商数来完成巴雷特规约。

当对一个2n位数进行n位数的除法时，递归版本的成本为O(M(n) log(n))，其中M(n)是乘以n位数的成本。

使用巴雷特规约的版本将花费O(M(n))，如果使用牛顿算法计算逆，则之后的隐藏常量会更大，因此仅适用于非常大的除法。

更详细地说，大多数除法算法背后的核心算法是“带有规约的估计商数”计算，计算出(q,r)，使得：

x = qy + r

但不限制 0 <= r < y。典型的循环如下：

估算商 q = x/y
计算相应的余数 r = x - qy
可选择调整商，使余数 r 在某个期望区间内
如果 r 太大，则使用 r 替换 x 重复上述过程。

x/y 的商将是所有生成的 q 的总和，而最终值的 r 将是真正的余数。

例如，学校长除法就是这种形式。例如，步骤3涵盖了那些你猜测的数字太大或太小的情况，并且你需要调整它以获得正确的值。

分治方法通过计算 x 和 y 的前导数字来估算 x/y 的商。可以通过调整它们的大小来进行优化，但我记得如果 x' 是 y' 的两倍数字位数，则可以获得最佳结果。

在我看来，如果你坚持使用整数算术，那么乘以倒数的方法是最简单的。基本方法如下：

使用 m = floor(2^k / y) 估算 y 的倒数
使用 q = 2^(i+j-k) floor(floor(x / 2^i) m / 2^j) 估算 x/y

事实上，如果这意味着您可以使用更快的倒数实现，则实际实现可以容忍额外的 m 中的误差。

误差很难分析，但如果我记得正确的话，您需要选择 i 和 j，以便由于误差如何积累，x ~ 2^(i+j)，并且您需要选择 x / 2^i ~ m^2 来最小化整体工作量。

随后的余数将具有 r ~ max(x/m, y)，因此这给出了选择 k 的经验法则：您希望 m 的大小大约为每次迭代计算的商位数，或者等效地，您希望从 x 中删除的位数。