为什么这个使用 gcc、-mfpmath=387 和优化级别 -O2 或 -O3 编译的简单程序会产生 NaN 值？

Question

我有一个简短的程序，它执行数值计算，并在某些特定条件成立时获得不正确的 NaN 结果。我看不出这个 NaN 结果是如何产生的。请注意，我没有使用允许重新排序算术运算的编译器选项，例如-ffath-math.

\n

问题：我正在寻找 NaN 结果如何产生的解释。从数学上讲，计算中没有任何内容会导致除以零或类似的结果。我错过了一些明显的东西吗？

\n

请注意，我不是问如何解决问题\xe2\x80\x94，这很简单。我只是想了解 NaN 是如何出现的。

\n

最小的例子

\n

请注意，这个示例非常脆弱，甚至需要进行很小的修改，例如添加printf()在循环中添加调用以观察值）也会改变行为。这就是为什么我无法进一步减少它。

\n

// prog.c\n\n#include <stdio.h>\n#include <math.h>\n\ntypedef long long myint;\n\nvoid fun(const myint n, double *result) {\n    double z = -1.0;\n    double phi = 0.0;\n    for (myint i = 0; i < n; i++) {\n        double r = sqrt(1 - z*z);\n\n        /* avoids division by zero when r == 0 */\n        if (i != 0 && i != n-1) {\n            phi += 1.0 / r;\n        }\n\n        double x = r*cos(phi);\n        double y = r*sin(phi);\n\n        result[i + n*0] = x;\n        result[i + n*1] = y;\n        result[i + n*2] = z;\n\n        z += 2.0 / (n - 1);\n    }\n}\n\n#define N 11\n\nint main(void) {\n    // perform computation\n    double res[3*N];\n    fun(N, res);\n\n    // output result\n    for (int i=0; i < N; i++) {\n        printf("%g %g %g\\n", res[i+N*0], res[i+N*1], res[i+N*2]);\n    }\n\n    return 0;\n}\n

\n

编译：

\n

gcc -O3 -mfpmath=387 prog.c -o prog -lm\n

\n

输出的最后一行是：

\n

nan nan 1\n

\n

我期望一个接近于零的数字，而不是 NaN。

\n

该示例的关键特征

\n

要出现 NaN 输出，以下条件必须全部成立：

\n

\n

• 在 x86 平台上使用 GCC 进行编译。我能够在 macOS 10.14.6 上使用 GCC 12.2.0（来自 MacPorts）以及在 Linux（openSUSE Leap 15.3）上使用 GCC 版本 9.3.0、8.3.0 和 7.5.0 进行重现。

  \n

  我不能在 Linux 上使用 GCC 10.2.0 或更高版本，或在 macOS 上使用 GCC 11.3.0 重现它。

  \n
\n

• 选择使用 x87 指令-mfpmath=387，优化级别为-O2或-O3。

  \n
\n

• myint必须是签名的64 位类型。

  \n
\n

• 想着resultn×3 矩阵，它必须按列优先顺序存储。

  \n
\n

• printf()主循环中没有调用fun()调用。

  \n
\n

\n

如果没有这些功能，我确实会得到预期的输出，即类似1.77993e-08 -1.12816e-08 1或0 0 1作为最后一行的内容。

\n

程序说明

\n

尽管这对问题来说并不重要，但我还是对该程序的功能进行了简短的解释，以使其更容易理解。它计算x球体表面上特定排列的点的y、 、z三维坐标。值以相等的增量从 -1 到 1，但是，由于数值舍入误差，最后一个值不会精确为 1。坐标被写入一个×3 矩阵 ，以列优先顺序存储。和nznresultrphi是 (x, y) 平面中的极坐标。

\n

请注意，当zis-1或1thenr变为 0 时。这发生在第一个和最后一个迭代步骤中。这将导致表达式中除以 0 1.0 / r。但是，1.0 / r被排除在循环的第一次和最后一次迭代之外。

\n

Answer 1

这是由 x87 80 位内部精度、GCC 的不一致行为以及编译器版本之间不同的优化决策的相互作用引起的。

x87 仅支持 IEEE 二进制 32 和二进制 64 作为存储格式，在加载/存储时与其 80 位表示形式相互转换。为了使程序行为可预测，C 标准要求在赋值时放弃额外的精度，并允许通过宏检查中间精度FLT_EVAL_METHOD。其中-mfpmath=387,FLT_EVAL_METHOD为 2，因此您知道中间精度对应于类型long double。

不幸的是，GCC 不会降低赋值的额外精度，除非您通过-std=cNN（而不是-std=gnuNN）或显式传递 来请求更严格的一致性-fexcess-precision=standard。

在您的程序中，该z += 2.0 / (n - 1);语句应通过以下方式计算：

1. 2.0 / (n - 1)以中间 80 位精度进行计算。
2. 添加到之前的值z（仍为 80 位精度）。
3. 舍入到声明的类型z（即binary64）。

在以 NaN 结尾的版本中，GCC 执行以下操作：

1. 2.0 / (n - 1)在循环之前仅计算一次。
2. 将此分数从binary80 舍入到binary64 并存储在堆栈中。
3. 在循环中，它从堆栈中重新加载该值并添加到z。

这是不符合要求的，因为它2.0 / (n - 1)经历了两次舍入（首先到binary80，然后到binary64）。

上面解释了为什么您会看到不同的结果，具体取决于编译器版本和优化级别。然而，一般来说，您不能期望您的计算在最后一次迭代中不会产生 NaN。当n - 1不是 2 的幂时，2.0 / (n - 1)无法精确表示并且可以向上舍入。在这种情况下，“z”的增长速度可能比真实的 sum 快一点-1.0 + 2.0 / (n - 1) * i，并且最终可能会高于 1.0 i == n - 1，导致sqrt(1 - z*z)由于负参数而产生 NaN。

事实上，如果您在程序中更改#define N 11为#define N 12，您将确定性地获得具有 80 位和 64 位中间精度的 NaN。

Answer 2

... NaN 结果是如何产生的（？）

尽管更好地遵守 C 规范可能显然可以解决 OP 的直接问题，但我断言应该考虑其他预防实践。

sqrt(1 - z*z)当 时 是候选 NaN |z| > 1.0。

除以零的索引测试预防可能还不够，然后导致cos(INFINITE)，另一种 NaN 可能性。

// /* avoids division by zero when r == 0 */
//    if (i != 0 && i != n-1) {
//        phi += 1.0 / r;
//    }

为了避免这些问题，1）直接测试，2）使用更精确的方法。

if (r) {
  phi += 1.0 / r;
}

// double r = sqrt(1 - z*z);
double rr = (1-z)*(1+z);  // More precise than 1 - z*z
double r = rr < 0.0 ? 0.0 : sqrt(rr);