标准对 std::pow、std::log 等 cmath 函数有啥看法？

Posted 2023-02-22

【中文标题】标准对 std::pow、std::log 等 cmath 函数有啥看法？

【英文标题】：What does standard say about cmath functions like std::pow, std::log etc?标准对 std::pow、std::log 等 cmath 函数有什么看法？

【发布时间】：2015-10-18 18:17:47

【问题描述】：

标准是否保证函数在所有实现中返回完全相同的结果？

以pow(float,float) 为例，用于 32 位 IEEE 浮点数。如果传入相同的两个浮点数，所有实现的结果是否相同？

或者根据用于实现pow 的算法，标准是否允许一些关于微小差异的灵活性？

【问题讨论】：

标准不强制 float 和 double 为 IEEE-754 类型。它们可以是特定于机器的类型，具有不同的范围和精度。

如果您有两台机器对floats 和doubles 使用相同的表示，并且在舍入方面具有相同的设置，我不知道答案。在所有其他情况下，答案是否定的，它们会给出不同的结果。

在遵循 ieee754 的同一台机器上还有两个不同的编译器。如果 std 说是的答案应该是相同的，那么它必须定义正确的答案是什么。

相关：Does any floating point-intensive code produce bit-exact results in any x86-based architecture?：摘要：C 中没有，但 asm 中是。

【参考方案1】：

不，C++ 标准不要求 cmath 函数的结果在所有实现中都相同。对于初学者，您可能无法获得 IEEE-754/IEC 60559 浮点算法。

也就是说，如果实现确实使用 IEC 60559 并定义 __STDC_IEC_559__，那么它必须遵守 C 标准的附件 F（是的，您的问题是关于 C++，但 C++ 标准遵循 C 标头的 C 标准，例如 math.h)。附件 F 指出：

float 类型与 IEC 60559 单一格式匹配。 double 类型与 IEC 60559 双格式匹配。 long double 类型匹配 IEC 60559 扩展格式，否则 非 IEC 60559 扩展格式，否则为 IEC 60559 double 格式。

此外，它说正常的算术必须遵循 IEC 60559 标准：

+、−、* 和 / 运算符提供 IEC 60559 加、减、乘和除运算。

进一步要求sqrt 遵循IEC 60559：

<math.h> 中的 sqrt 函数提供 IEC 60559 平方根运算。

然后继续描述其他几个浮点函数的行为，您可能对这个问题不感兴趣。

最后，它到达math.h 标头，并指定各种数学函数（即sin、cos、atan2、exp 等）应如何处理特殊情况（即@987654341 @ 返回 ±0、atanh(x) 返回 NaN 并引发 |x| > 1 的“无效”浮点异常等）。但它从未确定正常输入的精确计算，这意味着您不能依赖所有实现产生精确相同的计算。

所以不，它不需要这些函数在所有实现中表现相同，即使所有实现都定义了__STDC_IEC_559__。

---

这都是从理论的角度来看的。在实践中，情况甚至更糟。 CPU 通常实现 IEC 60559 算术，但它可以有不同的舍入模式（因此结果会因计算机而异），并且编译器（取决于优化标志）可能会做出一些不符合严格标准的假设浮点运算。

因此，在实践中，它甚至比理论上更不严格，而且您很可能会看到两台计算机在某些时候产生略微不同的结果。

---

glibc 是一个真实的例子，它是 GNU C 库实现。 They have a table of known error limits for their math functions 跨不同的 CPU。如果所有 C 数学函数都是位精确的，那么这些表都将显示 0 个错误 ULP。但他们没有。表格显示他们的 C 数学函数确实存在不同程度的错误。我觉得这句话是最有趣的总结：

除了某些函数，如sqrt、fma 和rint，其结果完全通过引用相应的 IEEE 754 浮点运算以及字符串和浮点之间的转换来指定，GNU C 库不旨在为数学库中的函数正确舍入结果[...]

glibc 中唯一位精确的东西是 C 标准的附件 F 要求位精确的东西。正如您在他们的表格中看到的那样，大多数情况并非如此。

【讨论】：

感谢您的详细解答。因此，为了保证相同的结果，需要编写自己的 pow 函数。

差不多。事实上，对于位精确性很重要的情况（例如在音频和视频编解码器中），实现自己的定点算法并不少见。

你的句子中的 正常输入 是什么意思“但它从来没有确定正常输入的确切计算”？

@Zboson：我的意思是不应该产生任何错误或需要任何特殊处理的输入。例如，对于函数acos，0.5 将被视为正常输入，因为该输入值的反余弦定义良好。但是输入1 需要特殊处理，因为它导致值为零，并且IEC 60559 有两个零：+0 和-0。附件 F 声明 acos(1) 应该产生 +0，而不是 -0。但附件 F 没有提及acos(0.5)，因为那个“正常”输入值 不需要任何特殊处理。这样就清楚了吗？

@Zboson：确实，您不能依赖 acos(0.5) 在不同的编译器中是完全相同的值。哎呀，你甚至不能依赖它与不同架构上的相同编译器相同。 glibc has a nice table listing the ULP error of their math functions，如您所见，这些数学函数并非在所有架构中都是位精确的。