17

我正在用 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

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

4

5 回答 5

6

您肯定知道现有的任意精度包(例如,http ://gmplib.org/ )以及它们是如何运行的吗?它们通常被设计为“尽可能快地”运行以实现任意精度。

如果您将它们专门用于固定大小(例如,应用 [手动]部分评估技术来折叠常量和展开循环),我希望您能够获得非常好的例程,用于您想要的特定固定大小精度。

另外,如果您还没有看过它,请查看 D. Knuth 的Seminumerical Algorithms和 oldie but really goodie,它提供了多精度算术的关键算法。(大多数软件包都基于这些想法,但 Knuth 有很好的解释,而且非常正确)。

关键思想是将多精度数字视为非常大的基数(例如,基数 2^64),并将标准的 3 级算术应用于“数字”(例如 64 位字)。除法包括“估计商(大基数)数字,将估计乘以除数,从被除数中减去,左移一位,重复”,直到你得到足够的数字来满足你。对于除法,是的,它都是无符号的(在包装器中进行符号处理)。基本技巧是很好地估计商数(使用处理器为您提供的单精度指令),并进行快速多精度乘以个位数。有关详细信息,请参阅高德纳。请参阅有关多精度算术的技术研究论文(您可以进行一些研究)以获得奇异的(“尽可能快的”)改进。

于 2014-01-02T04:19:59.053 回答
1

对于您提到的大量数据类型,“大基数”方法更有效,特别是如果您可以用汇编语言执行 128 位除以 64 位指令。

虽然 Newton-Raphson 迭代确实很快收敛,但每次迭代都需要大量的乘法和累加步骤。

于 2014-01-16T06:12:08.203 回答
1

对于乘法,请看这里:

http://www.math.niu.edu/~rusin/known-math/99/karatsuba http://web.archive.org/web/20141114071302/http://www.math.niu.edu/~rusin /known-math/99/karatsuba

基本上,它允许使用三个(而不是四个)512 x 512 位乘法来进行 1024 x 1024 乘法。或九个 256 x 256 位,或二十七个 128 x 128 位。即使对于 1024 x 1024,增加的复杂性也可能无法击败蛮力,但可能对于更大的产品。这是最简单的“快速”算法,使用 n ^ (log 3 / log 2) = n^1.585 乘法。

我建议不要使用汇编程序。使用内联汇编器实现 64 x 64 -> 128 位乘法,与 add-with-carry 相同(我认为 gcc 和 clang 现在可能有内置操作);那么您可以例如并行乘以 n 位 x 256 位(任意数量的字乘以 4 个字),避免乘法的所有延迟,而不会对汇编程序发疯。

于 2014-03-18T19:24:18.413 回答
0

对于大量位,我了解到最快的算法是这样的:不是除 x / y,而是计算 1 / y 并乘以 x。计算 1 / y:

1 / y is the solution t of (1 / ty) - 1 = 0.
Newton iteration: t' = t - f (t) / f' (t) 
  = t - (1 / ty - 1) / (-1 / t^2 / y)
  = t + (t - t^2 y)
  = 2t - t^2 y

牛顿迭代二次收敛。现在的诀窍:如果你想要 1024 位精度,你从 32 位开始,一个迭代步骤给出 64 位,下一个迭代步骤给出 128 位,然后是 256,然后是 512,然后是 1024。所以你做了很多次迭代,但只有最后一个一个使用全精度。所以总而言之,你做了一个 512 x 512-> 1024 乘积 (t^2),一个 1024 x 1024 -> 1024 乘积 (t^2 y = 1 / y),另一个 1024 x 1024 乘积 (x * ( 1 / 年))。

当然,您必须非常准确地弄清楚每次迭代后的错误是什么;您可能必须从 40 位开始,在每一步中损失一点精度,这样最后您就足够了。

我不知道这会比你在学校学到的简单的蛮力除法运行得更快。并且 y 的位数可能少于全部位数。

于 2014-03-18T19:04:28.983 回答
0

另一种选择是蛮力。您可以取 x 的最高 128 位,除以 y 的最高 64 位,得到商的最高 64 位 r,然后从 x 中减去 rxy。并根据需要重复,仔细检查错误有多大。

分区很慢。所以你计算一次 2^127 /(y 的最高 64 位)。然后估计下一个 64 位,将 x 的最高 64 位乘以这个数字,并将所有内容移到正确的位置。乘法比除法快得多。

接下来你会发现所有这些操作都有很长的延迟。例如,5 个周期得到一个结果,但您可以在每个周期进行一次乘法运算。所以:估计结果的 64 位。在 x 的高端开始减去 r * y,这样您就可以尽快估计下一个 64 位。然后你同时从 x 中减去两个或多个乘积,以避免延迟带来的损失。实现这一点很困难。即使对于 1024 位(仅 16 个 64 位整数),有些事情可能也不值得。

于 2014-03-18T19:14:11.510 回答