2

以下 3 行使用"gcc -Ofast -march=skylake"给出了不精确的结果：

int32_t  i = -5;
const double  sqr_N_min_1 = (double)i * i;
1. - ((double)i * i) / sqr_N_min_1

显然，sqr_N_min_1gets25.和在第 3 行(-5 * -5) / 25应该变成1.使得第 3 行的总体结果正好是0.。事实上，编译器选项"gcc -O3 -march=skylake"确实如此。

但是使用“-Ofast”，最后一行产生-2.081668e-17而不是0.和与其他i（-5例如6或7）它得到其他非常小的正或负随机偏差0.。我的问题是：这种不精确的根源究竟在哪里？

为了调查这一点，我用 C 编写了一个小测试程序：

#include <stdint.h>      /* int32_t */
#include <stdio.h>
#define MAX_SIZE 10

double W[MAX_SIZE];

int main( int argc, char *argv[] )
{
  volatile int32_t n = 6; /* try 6 7 or argv[1][0]-'0' */
  double           *w = W;
  int32_t          i = 1 - n;
  const int32_t    end = n - 1;
  const double     sqr_N_min_1 = (double)i * i;

  /* Here is the crucial part. The loop avoids the compiler replacing it with constants: */
  do {
    *w++ = 1. - ((double)i * i) / sqr_N_min_1;
  } while ( (i+=2) <= end );

  /* Then, show the results (only the 1st and last output line matters): */
  w = W;
  i = 1 - n;
  do {
    fprintf( stderr, "%e\n", *w++ );
  } while ( (i+=2) <= end );

  return( 0 );
}

Godbolt 向我展示了由"-Ofast -march=skylake"与"-O3 -march=skylake"选项的"x86-64 gcc9.3"生成的程序集。请检查网站的五个栏目（1.源代码，2.带有“-Ofast”的程序集，3.带有“-O3”的程序集，4.第一个程序集的输出，5.第二个程序集的输出）：

有五列的 Godbolt 站点

如您所见，程序集的差异很明显，但我无法弄清楚不精确的确切来源。那么，问题是，哪些汇编指令对此负责？

一个后续问题是：是否有可能通过重新编写 C 程序来避免使用“-Ofast -march=skylake”的这种不精确性？

4

2 回答 2

6

评论和另一个答案指出了您的情况正在发生的具体转变，即倒数和 FMA 而不是除法。

是否有可能通过重新编写 C 程序来避免“-Ofast -march=skylake”的这种不精确性？

一般不会。

-Ofast是（当前）的同义词-O3 -ffast-math。
见https://gcc.gnu.org/wiki/FloatingPointMath

-ffast-mathis的一部分-funsafe-math-optimizations，顾名思义，可以改变数值结果。（目标是允许更多优化，例如将 FP 数学视为关联以允许使用 SIMD 自动矢量化数组的总和，和/或展开多个累加器，或者甚至只是在一个表达式中重新排列一系列操作以组合两个单独的常数。）

这正是您通过使用该选项所要求的那种速度超过精度的优化。如果您不想这样，请不要启用所有-ffast-math子选项，只启用-fno-math-errno/之类的安全选项-fno-trapping-math。（请参阅如何强制 GCC 假设浮点表达式为非负数？）

没有办法制定您的来源来避免所有可能的问题。

可能您可以volatile在所有地方使用 tmp 变量来破坏语句之间的优化，但这会使您的代码比-O3使用 default 的常规代码慢-fno-fast-math。即使这样，对库函数的调用sin也log可能会解析为假设 args 是有限的版本，而不是 NaN 或无穷大，因为-ffinite-math-only.

-Ofast 的 GCC 问题？指出另一个效果：isnan()被优化为 compile-time 0。

于 2021-05-10T11:15:01.680 回答
3

从评论看来，对于-O3，编译器计算1. - ((double)i * i) / sqr_N_min_1：

  1. 转换i成双倍并平方。
  2. 除以sqr_N_min_1。
  3. 从 1 中减去它。

并且，对于-Ofast，计算它：

  1. 在循环之前，计算 的倒数sqr_N_min_1。
  2. 转换i成双倍并平方。
  3. 计算 1 的融合乘减法减去平方乘以倒数。

后者提高了速度，因为它只计算一次除法，并且乘法比目标处理器中的除法快得多。最重要的是，融合运算比单独的乘法和减法更快。

发生错误是因为倒数运算引入了原始表达式中不存在的舍入误差（1/25 不能完全以二进制格式表示，而 25/25 当然是）。这就是为什么编译器在尝试提供严格的浮点语义时不进行这种优化的原因。

此外，只需将倒数乘以 25 即可消除错误。（这有点“偶然”，因为舍入误差以复杂的方式变化。1./25*25产生 1，但1./49*49不是。）但是融合运算产生了更准确的结果（它产生的结果就好像产品是精确计算的一样，发生了舍入只有在减法之后），所以它保留了错误。

于 2021-05-10T11:04:09.540 回答