c - 有没有办法用 i387 fsqrt 指令获得正确的舍入？

Question

有什么方法可以使用 i387 fsqrt 指令获得正确的舍入吗？...

...除了更改 x87 控制字中的精度模式 - 我知道这是可能的，但这不是一个合理的解决方案，因为它有令人讨厌的重入类型问题，其中精度模式将是如果 sqrt 操作被中断，则错误。

我正在处理的问题如下:x87 fsqrt 操作码以 fpu 寄存器的精度执行正确舍入(根据 IEEE 754)平方根运算，我将假设是扩展(80 位)精度。但是，我想用它来实现高效的单精度和 double 平方根函数，结果正确舍入(根据当前舍入模式)。由于结果精度过高，第二步将结果再次转换为单精度或 double 舍入，可能会留下不正确舍入的结果。

通过一些操作，可以通过偏差来解决这个问题。例如，我可以通过添加 2 的幂形式的偏差来避免加法结果的精度过高，该偏差将 double 值的 52 位有效位强制为 63 位扩展精度尾数的最后 52 位.但我没有看到任何明显的方法可以用平方根来做这样的把戏。

有什么好主意吗？

(也标记为 C，因为预期的应用程序是 C sqrt 和 sqrtf 函数的实现。)

Answer 1

首先，让我们避开显而易见的问题:您应该使用 SSE 而不是 x87。上证所sqrtss和 sqrtsd所有现代 x86 系统都支持所有现代 x86 系统，并且速度也快得多。

现在，如果您坚持使用 x87，我将从好消息开始:您无需为 float 做任何事情。你需要2p + 2位以计算 p 位浮点格式的正确舍入平方根。因为80 > 2*24 + 2 ，对单精度的额外舍入将始终正确舍入，并且您有一个正确舍入的平方根。

现在是坏消息:80 < 2*53 + 2 ，所以 double 没有这样的运气。我可以建议几种解决方法；这是我想到的一个很好的简单方法。

1. 让y = round_to_double(x87_square_root(x));
2. 使用 Dekker(头尾)积计算 a和 b这样 y*y = a + b正是。
3. 计算残差r = x - a - b .
4. if (r == 0) return y
5. if (r > 0) , 让y1 = y + 1 ulp , 并计算 a1 , b1英石。 y1*y1 = a1 + b1 .比较 r1 = x - a1 - b1至 r ，并返回 y或 y1 ，取决于哪个具有较小的残差(或者如果残差大小相等，则低位为零的那个)。
6. if (r < 0) , 对 y1 = y - 1 ulp 做同样的事情.

这个过程只处理默认的舍入模式；然而，在定向舍入模式中，简单地舍入到目标格式就可以了。