我应该使用什么算法进行高性能大整数除法？

Question

我将大整数编码到size_t 的数组中。我已经进行了其他操作（加、减、乘）；以及除以一位数。但如果可能的话，我想匹配我的乘法算法的时间复杂度（目前是 Toom-Cook）。

我收集了一些线性时间算法来获取我的股息乘法逆的各种概念。这意味着理论上我可以在与乘法相同的时间复杂度下实现除法，因为无论如何线性时间运算相比之下都是“微不足道的”。

我的问题是，我该怎么做呢？哪种类型的乘法逆在实践中最好？取模64^digitcount?当我将乘法逆元乘以除数时，我可以避免计算由于整数截断而被丢弃的数据部分吗？任何人都可以提供 C 或 C++ 伪代码或准确解释应该如何完成吗？

或者有没有比基于逆的方法更好的专用除法算法？

编辑：我挖掘了上面提到的“逆”方法。在“计算机编程艺术，第 2 卷：半数值算法”的第 312 页上，Knuth 提供了“算法 R”，它是一种高精度倒数。他说它的时间复杂度低于乘法。然而，将其转换为 C 并对其进行测试并非易事，并且在我编写代码之前不清楚将消耗多少开销内存等，这需要一段时间。如果没有人超过我，我会发布它。

Answer 1

GMP 库通常是好的算法的一个很好的参考。他们的documented algorithms for division主要取决于选择一个非常大的基数，因此您将一个4位数字除以一个2位数字，然后通过长除法进行。

长除法需要计算 2 位乘 1 位的商；这可以递归完成，也可以通过预先计算逆并估计商来完成，就像使用 Barrett 约简一样。

将2n 位数除以n 位数时，递归版本的成本为O(M(n) log(n))，其中M(n) 是n 位数相乘的成本。

如果使用牛顿算法计算逆，使用 Barrett 约简的版本将花费 O(M(n))，但根据 GMP 的文档，隐藏常数要大得多，因此这种方法只适用于非常大的除法。

---

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

x = qy + r

但没有0 <= r < y 的限制。典型的循环是

• 估计q 的商x/y
• 计算对应的约简r = x - qy
• 可以选择调整商，以便减少 r 在某个所需的区间内
• 如果r 太大，则重复使用r 代替x。

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

例如，教科书长除法就是这种形式。例如第 3 步涵盖了您猜测的数字太大或太小的情况，您可以对其进行调整以获得正确的值。

分治法通过计算x'/y' 来估计x/y 的商，其中x' 和y' 是x 和y 的前导数字。通过调整它们的大小有很大的优化空间，但是如果x' 是y' 的数字的两倍，IIRC 你会得到最好的结果。

如果您坚持整数运算，IMO 乘以逆方法是最简单的方法。基本方法是

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

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

错误分析起来很痛苦，但是如果我回忆一下这样做的方法，您想选择i和j，以便x ~ 2^(i+j)由于错误的累积方式，您想选择@987654357 @ 以最小化整体工作。

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

Answer 2

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

我的 bigint 除法有点不同。

见bignum division。尤其是看看近似除法器和那里的 2 个链接。一个是我的定点除法器，其他的是带有测量值的快速乘法算法（如 NTT 上的 karatsuba、Schönhage-Strassen），以及指向我的 32 位 Base 的快速 NTT 实现的链接。

我不确定逆乘数是否可行。

它主要用于除数为常数的模运算。恐怕对于任意除法，获取 bigint inverse 所需的时间和运算可能比标准除法本身更大，但由于我不熟悉它我可能是错的。

我在实现中看到的最常用的除法器是 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 是用于存储 bigint 数的 32 位 DWORDs 的数量。

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

速度提升是因为它不使用乘法（如果你不计算位移）

如果您使用像 2^32（ALU 块）这样的大底数，那么您可以使用 32 位内置 ALU 操作以多项式风格重写整个数字。
这通常比二进制长除法更快，其想法是将每个 DWORD 处理为单个数字，或者递归地将使用的算术除以一半，直到达到 CPU 的能力。
见division by half-bitwidth arithmetics

最重要的是使用 bignums 进行计算

如果您已经优化了基本运算，那么复杂性可以进一步降低，因为随着迭代（改变基本运算的复杂性）子结果变得更小（改变基本运算的复杂性），基于 NTT 的乘法就是一个很好的例子。

开销可能会搞砸。

因此，运行时有时不会复制大的 O 复杂度，因此您应该始终测量阈值并使用更快的方法来计算已使用的位数，以获得最大性能并尽可能优化。