在Mac和Linux上，exp函数的结果略有不同。

Question

下面的C程序在我的Mac和Linux上产生了不同的结果。 我很惊讶，因为我认为libm的实现在某种程度上是标准化的。

#include<math.h>
#include<stdio.h>

int main()
{
  double x18=-6.899495205106946e+01;
  double x19 = exp(-x18);
  printf("x19     = %.15e\n", x19);
  printf("x19 hex = %llx\n", *((unsigned long long *)(&x19)));
}

在Mac上的输出是：

x19     = 9.207186811339878e+29
x19 hex = 46273e0149095886

还有在Linux平台上

x19     = 9.207186811339876e+29
x19 hex = 46273e0149095885

以下两个都是没有任何优化标志编译的：

gcc -lm ....

我知道不应该将浮点数进行严格的比较。

在调试过程中，这个问题出现了。遗憾的是，使用这个计算的算法证明是数值不稳定的，而这种微小的差异会导致最终结果的显著偏差。但这是一个不同的问题。

我只是惊讶于像exp这样的基本操作没有像IEEE 754规定的基本代数运算那样标准化。

对于不同机器或不同版本的libm实现，我可以依赖哪些精度假设？

---

由于下面的讨论，我使用mpmath以高于机器精度的方式计算值，并且我得到了两个更多数字的结果9.2071868113398768244，因此对于我的两个结果，最后一位已经是错误的。Linux上的结果可以通过向下舍入此值来解释，如果计算机使用向上舍入，则Mac上的结果也会出错。

- rocksportrocker

---

好的...标准操作已列在https://en.wikipedia.org/wiki/IEEE_754-1985#Standard_operations。不过你在问什么？ - Eugene Sh.

顺便提一下，您可以使用“％a”格式说明符以十六进制打印“精确”的浮点值。 - Jens Gustedt

啊，谢谢！我以前从没见过 %a。 - rocksportrocker

另外，还请查看红帽博客中Jonathan Wakely的文章《为什么<cstdlib>比您想象的更复杂》。Wakely是GCC的C ++标准库维护者之一。我认为<math.h>与<cmath>之间的区别是一个更有趣的案例研究。 - jww

我认为您应该明确指定使用的编译器、编译器标志和CPU架构。在许多设置中，最简单的math.h函数最终会被编译器替换为内部函数，直接映射到CPU指令。这可能意味着在内部执行80位x87操作，然后将其舍入为double类型，例如。 - cnettel

Answer 1

C99规范（其他版本应该类似）中指出：

J.3 实现定义的行为

1. 在此子条款列出的每个领域中，符合规范的实现都需要记录其选择的行为。以下是实现定义的内容：

...

J.3.6 浮点数

- 浮点操作和返回浮点结果的和库函数的精度（5.2.4.2.2）。

这意味着GNU libm和BSD libm可以具有不同级别的精度。很可能发生的情况是，在OSX上的BSD实现会四舍五入到最近的单位（在最后一位），而GNU实现则向下截取到下一个ULP。

- dlasalle

---

IEEE-754标准仅要求加法、减法、乘法、除法和平方根的正确舍入（误差<= 0.5 ULP）。其他操作不需要正确舍入。这两个结果在1 ULP范围内，因此误差在<= 1 ULP。 - casevh

我使用mpmath扩展了我的问题，并给出了更精确的结果，表明在Mac上误差超过1 ULP。 - rocksportrocker

@rocksportrocker 请看下面我的回答。它有更多的细节。 - casevh

Answer 2

IEEE-754的行为在二进制级别上进行规定。在使用Linux时，我发现Python的本地math库、mpmath和通过gmpy2实现的MPFR得到的值是相同的。然而，将其转换为十进制则会因这三种方法而异。

>>> import mpmath, gmpy2
>>> import mpmath, gmpy2, math
>>> x18=68.99495205106946
>>> x19=math.exp(x18)
>>> mp18=mpmath.mpf("68.99495205106946")
>>> mp19=mpmath.exp(mp18)
>>> gp18=gmpy2.mpfr("68.99495205106946")
>>> gp19=gmpy2.exp(gp18)
>>> x18 == mp18
True
>>> x18 == gp18
True
>>> x19 == mp19
True
>>> x19 == gp19
True
>>> print(x18, mp18, gp18)
68.99495205106946 68.9949520510695 68.994952051069461
>>> print(x19, mp19, gp19)
9.207186811339876e+29 9.20718681133988e+29 9.2071868113398761e+29

在转换为Python的任意精度整数形式后，这三个结果也显示为精确值。

>>> hex(int(x19))
'0xb9f00a484ac42800000000000'
>>> hex(int(mp19))
'0xb9f00a484ac42800000000000'
>>> hex(int(gp19))
'0xb9f00a484ac42800000000000'

所以（至少有一个）Linux数学库，mpmath和gmpy2.mpfr是一致的。

免责声明：我维护gmpy2并曾经为mpmath做出过贡献。

- casevh

---

值得注意的是，IEEE-754不是C标准（或C ++）的一部分，Linux也不需要它（考虑到OSX可以合法运行的硬件集合受限，可以说它的实现将使用IEEE-754）。 - dlasalle

@dlasalle 确实。虽然IEEE-754标准不是必需的，但GNU libc试图遵守IEEE-754-2008标准，因此许多/大多数Linux系统将符合该标准。但这并不是一个安全的假设。如果您的应用程序依赖于可重复的结果，包括超越函数，则应使用单独的库，例如MPFR。 - casevh