为什么这段代码胜过 rint() 以及如何保护它不受 -ffast-math 和朋友的影响?

Why does this code beat rint() and how to I protect it from -ffast-math and friends?

我正在寻找一种方法来保护某些代码免受跨 C 编译器工作的 -ffast-math(或 msvc/icc 等价物等)的影响。

我的内部循环正在搜索数据以寻找接近整数值(例如在 ~0.1 以内)的数字。数据值是有符号的,通常少于几千,没有 inf/nan。我发现最快的版本使用了一个带有大幻数的技巧:

 remainder = h - ((h+MAGIC)-MAGIC) ;

有人有办法保持上面关键行括号的优先顺序吗?这似乎比 rint(x) 快了 3 倍,所以我很好奇为什么它仍然有效。会不会跟矢量化有关?

大多数编译器 "simplify" 使用 -ffast-math 或等价物时表达式停止工作。我想保持性能(3 倍相当多),但也保持它的便携性(考虑到 MAGIC 取决于拥有正确的 ieee)。如果我添加一个 volatile 然后它变慢但似乎用 fast-math 给出正确的答案但是比 rint 慢:

 volatile t = h+MAGIC; t-=MAGIC;
 remainder = h - t;

下面是一个完整的例子。我尝试了一些像 __attribute__((optimize("-fno-associative-math"))) 这样的 gcc 东西,但它似乎不是最终适用于 icc/gcc/msvc/clang 等的正确方法。相关的 C99 标准 pragma 似乎也没有广泛使用。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

/*  */
union i_cast {double d; int i[2];};
#define MAGIC 6755399441055744.0
/* x86_64 for me today */
#define ENDIANLOC 0
#define REMAINDER_LUA  \
    {volatile union i_cast u; u.d = (h) + MAGIC; r = h - u.i[ENDIANLOC]; }

#define REMAINDER_MAGIC  r=(h - ((h+MAGIC)-MAGIC));
#define REMAINDER_RINT   r=(h - rint(h));
#define REMAINDER_TRUNC  r=(h - ( (h>0) ? ((int)(h+0.5)) : ((int)(h-0.5))) );
#define REMAINDER_FLOOR  r=(h - floor(h+0.5));
#define REMAINDER_REMAIN r=(remainder(h, 1.0));
#define REMAINDER_ROUND  r=(h - round(h));
#define REMAINDER_NEARBY r=(h - nearbyint(h));

#define block(MACRO) {                                                  \
    for(i=0 ; i<3 ; i++){                                               \
      gettimeofday(&start, NULL);                                       \
      n = 0;                                                            \
      for (k = 0; k < ng; k++) {                                        \
        h = mul * gv[k];                                                \
        MACRO                                                           \
          if ( (r*r) < tol ) n++;                                       \
      }                                                                 \
      gettimeofday(&end, NULL);                                         \
      double dt = (double)(end.tv_sec - start.tv_sec);                  \
      dt += (1e-6)*(double)(end.tv_usec - start.tv_usec);               \
      if(i==2)                                                          \
      printf("%20s %d indexed in %lf s %f ns/value\n",#MACRO,           \
             n,dt,1e9*dt/ng);                                           \
    }                                                                   \
  } 

int main(){
  struct timeval start, end;
  // Make some test data
  double h, r, tol = 0.02, mul = 123.4;
  int i, n, k, ng = 1024*1024*32;
  srand(42);
  double *gv = (double *) malloc(ng*sizeof(double));
  for(int i=0;i<ng;i++) { gv[i] = ((double)rand())/RAND_MAX*2.-1.; }
  // Measure some timing
  block(REMAINDER_MAGIC);
  block(REMAINDER_LUA);
  block(REMAINDER_RINT);
  block(REMAINDER_FLOOR);
  block(REMAINDER_TRUNC);
  block(REMAINDER_ROUND);
  block(REMAINDER_REMAIN);
  block(REMAINDER_NEARBY);
  free( gv );
  return 0;
}

对我来说,今天 gcc -O3 的输出是这样的:

     REMAINDER_MAGIC 9489537 indexed in 0.017953 s 0.535041 ns/value
       REMAINDER_LUA 9489537 indexed in 0.048870 s 1.456439 ns/value
      REMAINDER_RINT 9489537 indexed in 0.050894 s 1.516759 ns/value
     REMAINDER_FLOOR 9489537 indexed in 0.086768 s 2.585888 ns/value
     REMAINDER_TRUNC 9489537 indexed in 0.162564 s 4.844785 ns/value
     REMAINDER_ROUND 9489537 indexed in 0.417856 s 12.453079 ns/value
    REMAINDER_REMAIN 9489537 indexed in 0.517612 s 15.426040 ns/value
    REMAINDER_NEARBY 9489537 indexed in 0.786896 s 23.451328 ns/value

也许其他一些语言 (rust/go/opencl/whatever) 在这里会比 C 做得更好?或者最好控制编译器标志并在代码中添加运行时测试以确保正确性?

没有标准的方法来控制这种非标准行为。每个带有 -ffast-math-style 选项的编译器都有属性和编译指示来控制它,但它们各不相同,选项的精确效果也是如此。就此而言,其中一些编译器的不同版本具有 显着 不同的快速数学行为,因此这不仅仅是一组适当的 pragma 的问题。获得标准行为的标准方法是让编译器遵循语言标准。

-ffast-math 及其同类产品主要面向那些不关心浮点数学细节,只希望他们的程序(有限且保守地使用 FP)的程序员 运行 快点。在任何情况下,-ffast-math 的大部分有用效果都可以通过精心编写的代码来复制。