我正在用 x86-64 汇编语言编写代码库,为s0128
, s0256
, s0512
,s1024
有符号整数类型和f0128
, f0256
, f0512
,f1024
浮点类型提供所有常规的按位、移位、逻辑、比较、算术和数学函数。到目前为止,我正在研究有符号整数类型,因为浮点函数可能会调用一些为整数类型编写的内部例程。
到目前为止,我已经编写并测试了执行各种一元运算符、比较运算符以及加法、减法和乘法运算符的函数。
现在我正在尝试决定如何为除法运算符实现功能。
我的第一个想法是,“Newton-Raphson 一定是最好的方法”。为什么?因为给定一个好的种子(开始猜测)它会很快收敛,我想我应该能够弄清楚如何在操作数上执行本机 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% 的时间我们会发现我们需要从被除数中减去除数,所以无论如何我们都需要减去它们的整个宽度。在这种情况下,我们实际上执行了比传统方法更多的指令(我们首先将它们相减,然后测试标志以确定除数是否 <= 被除数)。因此,一半的时间我们实现储蓄,一半的时间我们实现亏损。在大型类型s1024
(如包含 -16-64 位组件)上,节省大量且损失很小,因此这种方法应该实现大量净节省。在像s0128
(包含 -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...)
注意:我的猜测是,实际的实现只会对无符号整数进行操作。在乘法的情况下,将负操作数转换为正数,然后执行无符号乘法,然后如果恰好一个原始操作数为负数,则将结果取反更容易和更有效(也许)。但是,如果它是有效的,我会考虑一个有符号整数算法。
注意:如果我的浮点类型 ( 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
这与我的代码库的“永不丢失精度”和“永不溢出”的策略一致(当答案因精度丢失或上溢/下溢而无效时返回错误)。但是,当调用双宽度返回值函数时,不会出现此类错误。