有没有办法使用 i387 fsqrt 指令获得正确的舍入?...
...除了改变精确模式在 x87 控制字中 - 我知道这是可能的,但这不是一个合理的解决方案,因为它存在令人讨厌的重入型问题,如果 sqrt 操作中断,精度模式将出错。
我正在处理的问题如下:x87fsqrt
操作码以 fpu 寄存器的精度执行正确舍入(根据 IEEE 754)平方根运算,我假设是扩展(80 位)精度。但是,我想用它来实现高效的单精度和双精度平方根函数,并且结果正确舍入(根据当前舍入模式)。由于结果精度过高,将结果转换为单精度或双精度的第二步再次舍入,可能会留下舍入不正确的结果。
通过一些操作,可以通过偏差来解决这个问题。例如,我可以通过添加 2 的幂形式的偏差来避免加法结果中的精度过高,该偏差强制双精度值的 52 个有效位进入 63 位扩展精度尾数的最后 52 位。但我没有看到任何明显的方法可以用平方根来完成这样的技巧。
有什么巧妙的想法吗?
(也标记为 C,因为预期的应用程序是 C 的实现sqrt
and sqrtf
功能。)
首先,让我们明确一点:您应该使用 SSE 而不是 x87。上交所sqrtss
and sqrtsd
指令完全按照您的要求执行,所有现代 x86 系统都支持,并且速度也显着加快。
现在,如果您坚持使用 x87,我将从好消息开始:您不需要为浮动做任何事情。你需要2p + 2
位以 p 位浮点格式计算正确舍入的平方根。因为80 > 2*24 + 2
,附加舍入到单精度将始终正确舍入,并且您有一个正确舍入的平方根。
现在坏消息是:80 < 2*53 + 2
,所以双精度就没有这样的运气了。我可以建议几种解决方法;这是我脑海中想到的一个简单的好方法。
- let
y = round_to_double(x87_square_root(x));
- 使用 Dekker(头尾)积来计算
a
and b
这样y*y = a + b
确切地。
- 计算残差
r = x - a - b
.
if (r == 0) return y
-
if (r > 0)
, let y1 = y + 1 ulp
,并计算a1
, b1
s.t. y1*y1 = a1 + b1
。比较r1 = x - a1 - b1
to r
,并返回y
or y1
,取决于哪一个具有较小的残差(或者低位为零的残差,如果残差大小相等)。
-
if (r < 0)
,做同样的事情y1 = y - 1 ulp
.
此过程仅处理默认舍入模式;然而,在定向舍入模式中,简单地舍入到目标格式就可以做到正确的事情。
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)