在 X87 和 SSE FPU 中用户定义的点之后截断浮点数和双精度数

Posted 2023-02-16

【中文标题】在 X87 和 SSE FPU 中用户定义的点之后截断浮点数和双精度数

【英文标题】：Truncate Floats and Doubles after user defined points in X87 and SSE FPUs

【发布时间】：2016-07-20 18:03:43

【问题描述】：

我已经创建了一个函数g，它能够在一定程度上逼近一个函数，这个函数给出了精确到小数点后 5 位的结果（1,23456xxxxxxxxxxxx，其中 x 位置只是舍入错误/垃圾）。

为避免将错误传播到将使用g 结果的其他计算，我想将所有 x 位置设置为零，更好的是，只需将小数点后 5 位的所有内容设置为 0。

我在 X87 和 SSE 文献中没有找到任何可以让我以我想要的方式使用 IEEE 754 位或其表示的内容。

对 X87 的 FISTP 指令有一个旧的引用，它显然在 SSE 世界中与 FISTTP 相映成趣，其好处是 FISTTP 不需要修改控制字，因此速度更快。

我注意到FISTTP 被称为“切碎模式”，但现在在更现代的文献中只是“向零舍入”或“截断”，这让我感到困惑，因为“切碎”意味着完全删除“向零舍入”的东西零”对我来说不一定意味着同样的事情。

我不需要四舍五入到零，我只需要保留最多 5 位小数作为函数的最后一步，然后将结果存储到内存中；如何在 X87（标量 FPU）和 SSE（矢量 FPU）中执行此操作？

【问题讨论】：

等等...我很困惑。您提到 IEEE 和 bits，这意味着您了解 FP 表示。但同时你问如何截断 decimals 这在二进制浮点中是不可能的。我错过了什么吗？

@user354900 我不同意，尽可能精确地工作，并截断或舍入显示。通过在计算过程中截断，您将“传播”更大的错误，而不是更少。

@user354900 IEEE 浮点中的“斩波模式”只是截断二进制表示。除非您恰好四舍五入为整数，否则它绝不会映射到十进制数字。

浮点值（通常）没有有“十进制值”——这是一种人类可读格式的转换。如果您不需要它们，那么您就不会对这些值进行任何进一步的计算。因此，您不妨保留精度并处理输出截断。如果您要进行更多计算，那么您确实需要最大精度。

@user354900 如果你真的在问这个问题。那我想你不知道浮点数学。所以我推荐你这个：docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

【参考方案1】：

正如几个人评论的那样，more early rounding doesn't help the final result be more accurate。如果您想阅读更多关于浮点比较和怪异/陷阱的信息，我强烈推荐 Bruce Dawson 的浮点系列文章。这是the one with the index的引述

我们终于在这个系列中达到了我一直在等待的点 为了。在this post我将分享最关键的部分 我拥有的浮点数学知识。这里是：

[浮点] 数学很难。

你简直不会相信它是多么巨大、巨大、令人难以置信的困难。 我的意思是，你可能认为很难计算火车从 芝加哥和洛杉矶会发生碰撞，但这只是小菜一碟 浮点数学。

（如果您认识到最后一段是对有关空间的著名行的释义，则可以加分。）

---

如何实际实施你的坏主意：

没有任何机器指令或 C 标准库函数可以截断或舍入到整数以外的任何值。

请注意，有些机器指令（和 C 函数）会将 double 舍入为最接近的（可表示的）整数，而不会将其转换为 intmax_t 或任何东西，只需 double->double。所以不会通过固定宽度的 2 的补码整数往返。

所以要使用它们，您可以将浮点数放大某个因子，四舍五入到最接近的整数，然后再缩小。 （就像 chux 的基于 round() 的函数，但我建议使用 C99 double rint(double) 而不是 round()。round 具有奇怪的舍入语义，与 x86 上的任何可用舍入模式都不匹配，因此它编译得更糟代码。

---

你一直提到的 x86 asm 指令没什么特别的，不要做任何你不能要求编译器用纯 C 做的事情。

FISTP (Float Integer STore (and Pop the x87 stack) 是编译器或 asm 程序员实现long lrint(double) or (int)nearbyint(double). 的一种方法@ 一些编译器为其中一个或另一个生成更好的代码。它使用当前的 x87 舍入模式进行舍入（默认：四舍五入），即same semantics as those ISO C standard functions。

FISTTP (Float Integer STore with Truncation (and Pop the x87 stack) is part of SSE3, 即使它在 x87 堆栈上运行。它让编译器避免将舍入模式设置为截断（向零舍入）实现(long)x的C截断语义，然后恢复旧的舍入模式。

这就是“不修改控制字”的内容。两条指令都没有这样做，但要在没有 FISTTP 的情况下实现 (int)x，编译器必须使用 other 指令来修改和恢复围绕 FIST 指令的舍入模式。或者只是使用SSE2 CVTTSD2SI 将 xmm 寄存器中的双精度值转换为截断，而不是传统 x87 堆栈上的 FP 值。

由于 FISTTP 仅适用于 SSE3，因此您只能将其用于 long double，或者用于在 x87 寄存器中具有 FP 值的 32 位代码，因为旧的 32 位 ABI 返回 FP x87 堆栈上的值。

---

PS。如果你不认识 Bruce 的 HHGTG 参考，the original is：

空间很大。真的很大。你不会相信有多么巨大 大得令人难以置信。我的意思是你可能认为这是一个很长的路要走 通往化学家之路，但那只是通往太空的路。

【讨论】：

可以使用double HHGTG_round(double) return 42;

【参考方案2】：

如何在 X87（标量 FPU）和 SSE（向量 FPU）中执行此操作？

以下不使用 X87 和 SSE。我已将其作为通用代码的社区参考。如果有的话，它可以用来测试 X87 解决方案。

g() 的结果的任何“斩波”至少会略微增加错误，希望可以容忍，因为 OP 说“为了避免将错误传播到其他计算......”

不清楚 OP 是否希望“精确到小数点后 5 位的结果”来反映绝对精度 (+/- 0.000005) 或相对精度 (+/- 0.000005 * 结果)。将假定“绝对精度”。

由于 float、double 通常是 二进制 浮点数，因此任何“印章”都会将 FP 数 最接近 的倍数0.00001.

文字方法：

//       - x    xxx...xxx    . xxxxx \0
char buf[1+1+ DBL_MAX_10_EXP+1  +5   +1];
sprintf(buf, "%.5f", x);
x = atof(buf);

round() rint() 方法：

#define SCALE 100000.0
if (fabs(x) < DBL_MAX/SCALE) 
  x = x*SCALE;
  x = rint(x)/SCALE;

x 的直接位操作。只需将有效数字中的选择位置零即可。

TBD code.

【讨论】：

我推荐C99 double rint(double) 而不是round()。 round 具有奇怪的舍入语义，与 x86 上的任何可用舍入模式都不匹配，so it compiles to worse code for that architecture.。 round 实际上可以内联到 ARM 的单个指令，因为似乎 ARM 具有具有这些语义的指令。

@Peter Cordes 在您的comment 中有好的想法。我认为rint() 和nearbyint() 都是竞争者。谢谢你的提示。根据您的需要修改帖子 - 这是社区 wiki。

哦，对，是的，double nearbyint(double) 应该总是至少一样高效，并且可能更好，因为它不必设置 errno。除此之外，我认为它们具有相同的语义。