有没有办法用i387fsqrt指令获得正确的舍入？

Question

有没有办法用i387fsqrt指令获得正确的舍入？...

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

我正在处理的问题如下:的x87 fsqrt 操作码以fpu寄存器的精度执行一个正确的四舍五入（每IEEE754）平方根操作，我将假设它是扩展（80位）精度。但是，我想用它来实现高效的单精度和双精度平方根函数，并将结果正确舍入（根据当前舍入模式）。由于结果具有多余的精度，第二步将结果再次转换为单精度或双精度轮，可能留下一个不正确的四舍五入的结果。

通过一些操作，可以通过偏见来解决这个问题。例如，我可以通过以2的幂的形式添加偏置来避免加法结果中的过量精度，该偏置迫使双精度值的52有效位进入63位扩展精度尾数的最后52位。但我没有看到任何明显的方法来用平方根做这样的伎俩。

有什么聪明的主意吗？

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

Answer 1

首先，让我们把显而易见的事情弄清楚:您应该使用SSE而不是x87。上证所 sqrtss 和 sqrtsd 指令完全符合您的要求，在所有现代x86系统上都得到支持，而且速度也要快得多。

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

现在坏消息: 80 < 2*53 + 2, ，所以双精度没有这样的运气。我可以建议几种解决方法;这是一个很好的简单的一个从我的头顶。

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 s.t. y1*y1 = a1 + b1.比较一下 r1 = x - a1 - b1 到 r, ，并返回任 y 或 y1, ，取决于哪个具有较小的残差（或低阶位为零的残差，如果残差大小相等）。
6. if (r < 0), ，做同样的事情 y1 = y - 1 ulp.

此proceure仅处理默认舍入模式;但是，在有向舍入模式中，简单地舍入到目标格式可以做正确的事情。

Answer 2

好吧，我想我有更好的解决方案:

1. 计算/计算 y=sqrt(x) 在扩展精度（fsqrt).
2. 如果最后11位不是 0x400, ，简单地转换为双精度并返回。
3. 添加 0x100-(fpu_status_word&0x200) 到扩展精度表示的低字。
4. 转换为双精度并返回。

步骤3是基于这样一个事实，即状态字的C1位（0x200）是1当且仅当 fsqrt结果被四舍五入。这是有效的，因为，由于在步骤2中的测试, x 不是一个完美的正方形;如果它是一个完美的正方形, y 不会有超过双精度的位。

使用条件浮点操作执行步骤3可能会更快，而不是处理位表示并重新加载。

这是代码（似乎在所有情况下都有效）:

sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret

Answer 3

它可能不是你想要的，因为它没有利用387 fsqrt 指令，但有一个令人惊讶的效率 sqrtf(float) 在 格里布克 用32位整数算术实现。它甚至可以正确处理NaNs，Infs，subnormals-可能可以用真正的x87指令/FP控制字标志来消除其中的一些检查。见: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

该 dbl-64/e_sqrt.c 代码并不那么友好。很难一目了然地说出正在做出的假设。奇怪的是，图书馆的i386 sqrt[f|l] 实现只需调用 fsqrt, ，但加载值不同。 flds 对于SP, fldl 为DP。