STDC中的(1.0e300 + pow(2.0, -30.0) > 1.0)具体是做什么的？

Question

我遇到了一个计算atan(x)的函数（源代码在这里）。简化问题并稍微重新格式化一下，他们有类似下面的代码：

static const double one   = 1.0,
                   huge   = 1.0e300;

double atan(double x)
{
  /* A lot of uninteresting stuff here */

  if (ix < 0x3fdc0000) {              /* |x| < 0.4375 */

    if (ix < 0x3e200000) {            /* |x| < 2^-29 */
      if ((huge + x) > one) return x; /* raise inexact */
    }

    id = -1;
  }

  /* A lot of more uninteresting stuff */
}

我非常感兴趣了解代码行 if ((huge + x) ... 的作用和工作原理。

根据注释，如果 x 的绝对值小于 2^-29，则表达式或比较会引发一个 inexact 错误。

我的第一个问题是我目前不明白为什么要这样做：如果使用该函数计算 arctan 时，当 x 的绝对值太小时会导致不精确的结果，为什么他们不只是使用类似于 if (fabs(x) < [some_value_here]) ... 的东西？我怀疑这只是因为在他们的硬件/库中以这种方式不会引发 inexact 警告，但我想确定一下。

假设我是正确的，我的第二个问题是我不明白为什么需要比较。我认为这里的关键点是将一个非常小的数字添加到一个非常大的数字中，以便此添加不会足够或根本不会改变大数字。因此，引发 inexact 警告的是加法，而不是比较。因此，我在想比较应该做什么。这只是为了强制编译器实际计算 (huge + x)，否则可能会被优化掉吗？

最后，如果有人能解释一下数学，我会很感激。选择 1.0e300 作为 huge 的值似乎是一个相当随意的选择。但这只是一个额外的问题，因为我承认我还没有完成我的作业中的 数学 部分（我对 double 值及其 IEEE754 表示并不陌生，但理解此代码的数学方面需要我一些时间，除非有人给出简短的解释）。 编辑1 无意中看到：

该函数的 float32 版本，包括上述奇怪的行，几乎字面上仍然在 glibc 2.19 中！由于 glibc 应该是可移植的，因此该代码也应该如此。它位于子目录 sysdeps\ieee754\flt-32 中，因此我认为这是 float32 函数的软件仿真，其中可移植性不是问题，因为硬件相关的怪异不会显示出来（我认为软件仿真会按照 IEEE754 中定义的完全引发这些异常）。

- Binarus

---

标题与代码不一致。标题中的 2^-30 暗示了对 2 和 -30 进行按位异或（尽管由于 C++ 的常常令人惊讶的运算符优先级，这并不是实际情况）。 :-) - Adrian McCarthy

你是完全正确的。我会更改标题。 - Binarus

Answer 1

if ((huge + x) > one) return x;的意图是生成一个浮点数不精确的异常，然后从例程中返回。

浮点异常不是陷阱或处理器异常。它只是意味着在浮点操作中发生了一些异常情况。然后取决于操作的情况会发生什么。特别地，浮点环境可能被设置为仅在特殊寄存器中引发标志并继续操作，提供数值结果的情况下，不精确的异常仅会引发标志。或者可能设置为不精确的异常会导致陷阱，并将程序控制重定向到陷阱处理程序。

这个实现atan的代码不知道浮点环境如何设置。也许它可以获取设置，但它不想麻烦。考虑到它已经决定无法精确计算反正切函数，最简单的方法是触发一个浮点不精确的异常，只需执行具有不精确结果的简单加法即可。这种不精确的加法将具有所需的不精确反正切函数行为-它将引发标志或者根据设置引起陷阱。

至于为什么要使用ix < 0x3e200000进行比较，原因不清楚。首先，ix已经调整为反映绝对值，而x没有，那么为什么不使用已准备好的ix而不是使用另一个操作来产生fabs(x)？此外，与浮点比较相比，整数比较通常需要更少的处理器资源，特别是在编写该代码时的处理器中。或者可能作者只是恰巧使用了一种方法而不是另一种方法，也许大部分代码都使用ix来操作浮点编码而不是x来操作浮点值，并且他们不想不必要地来回切换。这也可能是因为在可用十六进制浮点符号之前编写了该代码（因此我们可以写出x < 0x1p-29f），并且编译器不能很好地将十进制数字转换为浮点值，因此他们不想在源代码中写出浮点值。

这种代码类型具有问题并且高度依赖于其所编写的C实现。一般来说，可能没有保证从编译器那里得到(huge + x) > one会在程序执行期间被评估。编译器可能会在编译时评估它。虽然可以假设此代码是为特定的C实现编写的，他们知道编译器将在编译时评估它，或者将确保达到相同的结果，包括引发浮点不精确异常。 (huge + x) > one表面上似乎没有比huge + x本身更有用，但也许作者知道我们不知道的一些关于编译器的东西。

huge 并不需要是 1.0e300。任何一个大到使得 huge 和 x 的和无法精确表示的值都可以。

- Eric Postpischil

---

感谢迄今为止的解释！有一点：[...] 另外，整数比较通常比浮点数比较需要更少的处理器资源[...] 是的，我知道他们首先针对2^-29进行检查，并且我理解他们的方法。但是以下代码让我困惑。显然，只有当x < 2 ^ -29时，才会构建和比较huge + x（我终将在某一天理解），但是对于任何给定的x，huge的选择都会影响是否触发inexact（例如，1.0e300 vs 1.0e307）。那么，他们是如何得出值1.0e300的呢？ - Binarus

@Binarus：由于 x 小于 2^-29，它的最高位小于 2^-30 或更低。这段代码所编写的 double 格式只有 53 位有效数字。因此，如果一个数的最高位是 2^23，则其有效数字中最低的位最多只能是 2^(23-52) = 2^-29。因此，如果将任何大于等于 2^23 的数加到小于 2^-29 的正数上，结果就无法被准确表示，因此会出现不精确的异常。但是，在这一点上，x 可能是负数，在这种情况下，总和可能会从 2^23 稍微降低，并且还有另一个可用的位。使用 2^24 可以避免这种情况。 - Eric Postpischil

还有一件有趣的事情（我刚刚偶然看到）：那个函数的 float32 版本，包括我提出问题中讨论的技巧，几乎仍然字面上在 glibc 2.19 中！由于 glibc 应该是可移植的，所以代码也应该如此。它在子目录 sysdeps\ieee754\flt-32 中，因此我认为这是 float32 函数的软件仿真，其中可移植性不是问题（我认为软件仿真会像 IEEE754 中定义的那样引发这些异常）。我将把这一发现添加到我的问题中，以强调这种技巧仍在使用。 - Binarus

因此，将任何大于2^24的数字添加到x中都会导致不精确异常。它不需要是1e300。 - Eric Postpischil

你的评论解释得非常好，非常感谢！已接受和+1。 - Binarus