algorithm - 什么是最快的x86-64汇编语言除法算法？

Question

我正在用x86-64汇编语言编写一个代码库，为s0128，s0256，s0512，s1024有符号整数类型和f0128，f0256，f0512浮点类型提供所有常规的按位、移位、逻辑、比较、算术和数学函数。到目前为止，我正在研究有符号整数类型，因为浮点函数可能会调用一些为整数类型编写的内部例程。
到目前为止，我已经编写并测试了执行各种一元运算符、比较运算符以及加、减和乘运算符的函数。
现在我要决定如何实现除法运算符的函数。
我的第一个想法是，“牛顿-拉斐逊一定是最好的方法”。为什么？因为给定一个好的种子（开始猜测），它会很快收敛，我想我应该能够找出如何在操作数上执行本机64位除法指令以获得一个好的种子值。事实上，如果种子值精确到64位，要获得正确的答案只需：

`s0128` : 1~2 iterations : (or 1 iteration  plus 1~2 "test subtracts")
`s0256` : 2~3 iterations : (or 2 iterations plus 1~2 "test subtracts")
`s0512` : 3~4 iterations : (or 3 iterations plus 1~2 "test subtracts")
`s1024` : 4~5 iterations : (or 4 iterations plus 1~2 "test subtracts")

然而，对这个问题再多考虑一点，我就想知道了。例如，我记得我编写的核心例程，它对所有大整数类型执行乘法操作：

s0128 :   4 iterations ==   4 (128-bit = 64-bit * 64-bit) multiplies +  12 adds
s0256 :  16 iterations ==  16 (128-bit = 64-bit * 64-bit) multiplies +  48 adds
s0512 :  64 iterations ==  64 (128-bit = 64-bit * 64-bit) multiplies + 192 adds
s1024 : 256 iterations == 256 (128-bit = 64-bit * 64-bit) multiplies + 768 adds

对于更广泛的数据类型，操作的增长是相当可观的，尽管循环相当短且高效（包括缓存方面）。此循环只写入结果的每个64位部分一次，并且从不读取结果的任何部分以进行进一步处理。
这让我思考更传统的移位和减法除法算法是否更快，尤其是对于较大的类型。
我的第一个想法是：

result = dividend / divisor                  // if I remember my terminology
remainder = dividend - (result * divisor)    // or something along these lines

#1:在计算每一位时，如果除数小于或等于被除数，通常从被除数中减去除数。好吧，通常我们可以通过检查它们最重要的64位部分来确定除数肯定小于或肯定大于被除数。只有当ms64位部分相等时，例程才必须检查下一个较低的64位部分，只有当它们相等时，我们才必须检查更低的部分，以此类推。因此，在几乎每次迭代（计算结果的每一位）中，我们可以大大减少计算此测试所执行的指令。
#2:但是……平均来说，大约50%的时间我们需要从红利中减去除数，所以我们无论如何都需要减去它们的整个宽度。在这种情况下，我们实际执行的指令比我们在传统方法中执行的指令要多（我们首先减去它们，然后测试标志以确定除数<=被除数）。因此，一半的时间我们意识到节约，一半的时间我们意识到损失。在像 f1024（包含-16-64位组件）这样的大型类型上，节省的空间很大，损失也很小，因此这种方法应该实现较大的净节省。对于像 s1024（包含-2-64位组件）这样的小型类型，节省的空间很小，损失也很大，但不是很大。
所以，我的问题是，“什么是最有效的除法算法”，给出：

#1: modern x86-64 CPUs like FX-8350
#2: executing in 64-bit mode only (no 32-bit)
#3: implementation entirely in assembly-language
#4: 128-bit to 1024-bit integer operands (nominally signed, but...)

注意：我的猜测是，实际实现只对无符号整数进行操作。在乘法的情况下，将负操作数转换为正操作数，然后执行无符号乘法，如果原来只有一个操作数为负，则对结果求反（可能）会更容易、更有效。但是，如果有符号整数算法是有效的，我将考虑它。
注意：如果我的浮点类型（ s0128， f0128， f0256， f0512， f1024， ）的最佳答案不同，请解释原因。
注意：我的内部核心无符号乘法例程对所有这些整数数据类型执行乘法操作，产生双宽度结果。换句话说：

u0256 = u0128 * u0128     // cannot overflow
u0512 = u0256 * u0256     // cannot overflow
u1024 = u0512 * u0512     // cannot overflow
u2048 = u1024 * u1024     // cannot overflow

我的代码库为每种有符号整数数据类型提供了两种乘法版本：

s0128 = s0128 * s0128     // can overflow (result not fit in s0128)
s0256 = s0256 * s0256     // can overflow (result not fit in s0256)
s0512 = s0512 * s0512     // can overflow (result not fit in s0512)
s1024 = s1024 * s1024     // can overflow (result not fit in s1024)

s0256 = s0128 * s0128     // cannot overflow
s0512 = s0256 * s0256     // cannot overflow
s1024 = s0512 * s0512     // cannot overflow
s2048 = s1024 * s1024     // cannot overflow

这与我的代码库“永不丢失精度”和“永不溢出”的策略一致（当答案因精度丢失或溢出/下溢而无效时，将返回错误）。但是，当调用双宽度返回值函数时，不会发生此类错误。

Answer 1

你肯定知道现有的任意精确包（例如，http://gmplib.org/），它们是如何操作的？它们通常被设计为“尽可能快”地运行以获得任意精度。
如果您将它们专门用于固定大小（例如，应用[手动]partial evaluation技术来折叠常数和展开循环），我希望您能够获得所需类型的特定固定大小精度的非常好的例程。
如果你还没看过，看看d.knuth的Seminumerical Algorithms，和老的但是很好的，它提供了多精度算法的关键算法。（大多数软件包都基于这些想法，但是knuth有很好的解释，而且非常正确）。
其关键思想是将多精度数字视为非常大的基数（例如基数2^64），并对“数字”（例如64位字）应用标准的三级算法。除法由“估计商（大基数）数字、估计乘除数、从被除数中减去、左移一位数、重复”组成，直到你得到足够的数字来满足你。对于除法，是的，它都是无符号的（在包装器中进行符号处理）。最基本的技巧是很好地估计商位数（使用处理器提供给您的单精度指令），并用单位数快速进行多精度乘法运算。详情见Knuth。请参阅有关多精度算法的技术研究论文（您需要做一些研究）以了解奇异的（“最快的”）改进。