9

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

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

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

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

有什么聪明的主意吗？

（也标记为 C，因为预期的应用是 Csqrt和sqrtf函数的实现。）

4

3 回答 3

15

首先，让我们明确一点：您应该使用 SSE 而不是 x87。SSEsqrtss和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等。by*y = a + b
  3. 计算残差r = x - a - b。
  4. if (r == 0) return y
  5. if (r > 0), 让y1 = y + 1 ulp, 和 计算a1, b1st y1*y1 = a1 + b1。比较r1 = x - a1 - b1，r并返回y或y1，这取决于哪个具有较小的残差（如果残差大小相等，则返回低位为零的那个）。
  6. if (r < 0), 做同样的事情y1 = y - 1 ulp。

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

于 2012-03-13T18:17:35.457 回答
3

好的，我想我有一个更好的解决方案：

  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
于 2012-03-15T04:00:00.663 回答
0

它可能不是您想要的，因为它没有利用 387fsqrt指令，但是在使用 32 位整数运算实现sqrtf(float)的glibc中具有惊人的效率。它甚至可以正确处理 NaN、Infs、次规范 - 可以使用真正的 x87 指令/FP 控制字标志来消除其中一些检查。看：glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

dbl-64/e_sqrt.c代码不是那么友好。很难一眼看出做出了什么假设。奇怪的是，该库的 i386sqrt[f|l]实现只是调用fsqrt，但加载值的方式不同。flds对于 SP，fldl对于 DP。

于 2012-03-13T17:36:52.230 回答