浮点数平方根倒数方法正确舍入

Question

我已经使用牛顿-拉夫逊方法 (汇编语言) 实现了一个 32 位 IEEE-754 浮点数平方根，这种方法是基于找到平方根的倒数。我使用的是四舍五入的舍入方法。我的平方根算法只接受规范化的值和零，不接受非规范化的值或特殊值 （NaN、Inf等）。

我想知道如何实现正确的舍入（使用类似指令的汇编语言），以便我的结果对所有输入都正确（符合 IEEE-754 标准）。基本上，我知道如何测试我的结果是否正确，但我想调整以下算法，以便我获得正确的舍入结果。我应该向算法中添加哪些指令？

有关更多信息，请参见：确定浮点数平方根 谢谢！

- Veridian

Answer 1

只有大约20亿个与您描述相符的float。尝试使用C库中的sqrtf对它们进行比较，并检查所有差异。如果您担心，可以使用C库中的sqrt或sqrtl获得更高精度的平方根。然而，通常的C库会正确地舍入sqrt、sqrtf和sqrtl，因此直接比较应该可以解决问题。

- tmyklebu

---

不确定为什么这个被踩了，但作为一个初始测试，这似乎很好。对于2^31个输入的详尽测试是相当可行的，事实上，在测试单输入单精度数学函数时很常见。虽然第三方库使用作为参考可能存在漏洞，但在实践中，像sqrt()这样的东西几乎没有风险。而tmyklebu建议检查任何不匹配，而不是简单地假设参考必须是正确的。类似的方法是针对相关的SSE内部函数进行测试。 - njuffa

实际上，如果你考虑到指数只有两个值是有意义的（输入中的4的因子在结果中变为2的因子，仅影响结果的指数，这应该不是问题），那么实际上小于20亿的float数需要检查的只有33554432个。 - Walter Tross

@tmyklebu，我想知道需要添加哪些指令到我发布的算法中，以便我可以获得正确舍入的结果。我已经知道如何将我的结果与其他算法进行比较，但我需要的是确保倒数平方根算法获得正确舍入结果所缺少的内容。我认为，基于Sterbenz定理，我可能会得到一些答案。 - Veridian

为了仅测试尾数部分（这是较难的部分），确实可以限制测试以覆盖两个连续二进制指数中所有可能的操作数，例如 1.0 <= x < 4.0，这需要 16.7M 个测试向量。但为了确保一切正常，全面的测试非常必要。 - njuffa

@njuffa，请看下面的帖子，我在评论中提到了这种方法。但是，你说得对，我想要的是“如何实现”。我更新了问题以反映这一点。 - Veridian

Answer 2

为什么不对结果进行平方，如果结果不等于输入，则根据差异的符号添加或减去一个最低有效位，平方并检查是否会得到更好的结果？这里的更好指的是绝对差异更小。唯一可能棘手的情况是当与尾数相交的√2，但这可以一次性检查。

编辑

我意识到上述答案是不充分的。在32位FP中简单平方并与输入进行比较并不能给出足够的信息。假设y = your_sqrt(x)。您将y的平方与x进行比较，发现y的平方>x，从y减去1个LSB并获得z（在您的注释中称为y1），然后将z的平方与x进行比较，并发现不仅z的平方

从你的评论中，我怀疑你使用的是严格的32位硬件，但是假设你有一个可用于64位结果的32位乘法器（如果没有，可以构造出来）。如果您将y的23位尾数作为整数，并在前面放置1，然后将其乘以自身，您将得到一个数字，除了可能需要额外移位1之外，您可以直接将其与以同样方式处理的x的尾数进行比较。通过这种方式，您可以使用所有48位进行比较，并且可以在不进行任何近似的情况下决定abs(y^2-x)≷abs(z^2-x)。

如果您不确定是否距离最终结果仅相差一个LSB（但是您确定不会比这更远），则应重复上述步骤，直到y^2-x改变符号或达到0。但要注意边缘情况，这些情况本质上应该是由于尾数越过2的幂而调整指数的情况。

还可以记住，正浮点数可以作为整数正确比较，至少在那些1.0F为0x3f800000的机器上。

- Walter Tross

---

您可能希望使用尝试结果的全宽乘积来计算残差，以获取尾随位。这就是为什么 FMA 如此有用的原因。在没有 FMA 的情况下，我已经成功地使用整数算术进行残差计算（请参见 CUDA 的 device_functions.h 文件中的 __fsqrt_rn()）。 - njuffa

@starbox：在没有溢出或下溢的情况下，两个不同的、相同符号的IEEE-754二进制浮点值不可能产生相同的平方。考虑一个正的s和s+u，其中u是s的ULP。那么(s+u)2就是s2 + 2su + uu。在这个区域内的ULP最多是2s乘以u，因此(s+u)2超过s超过一个ULP。因此，它不能舍入为与s2相同的浮点值。 - Eric Postpischil

@starbox，你问题的第一部分答案应该在我的回答中（本质上是循环，如果你知道迭代的最大次数，可以展开循环）。至于第二部分的答案……我真的不知道，但如果我是你，我会先通过tmyklebu的回答来探索情况。 - Walter Tross

tmyklebu的答案涉及采用公式(ss-foo)/s并进行最后一次牛顿迭代。我的方法涉及公式Xi(3/2-(1/2*Xi^2))。因此，我不确定如何将他的公式应用于这个倒数方法。 - Veridian

让我们在聊天中继续这个讨论：http://chat.stackoverflow.com/rooms/33781/discussion-between-starbox-and-walter-tross - Veridian