为什么此代码击败rint（），以及如何保护它免受-ffast-math和朋友的侵害？

Question

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

我的内部循环正在搜索数据以寻找接近整数值（例如，在〜0.1范围内）的数字。数据值是有符号的，通常少于几千，没有inf / nan。我发现最快的版本使用了一个带有大幻数的把戏：

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

有人想出一种方法来保持上方关键行括号的优先顺序吗？这似乎比rint(x)高出3倍，所以我对它为什么仍然起作用感到好奇。可能与向量化有关吗？

大多数编译器在使用-ffast-math或同等功能时会“简化”该表达式，并且它将停止工作。我想保持性能（3倍很多），但又要隐约地保持便携性（考虑到MAGIC取决于拥有正确的ieee）。如果我添加了挥发物，那么它的速度会变慢，但似乎可以通过快速方法给出正确的答案，但是比rint慢：

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

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

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

/* https://stackoverflow.com/questions/17035464/a-fast-method-to-round-a-double-to-a-32-bit-int-explained */
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 /其他）在这里会比C更好？还是更好地控制编译器标志并在代码中添加运行时测试以确保正确性？

Answer 1

没有控制这种非标准行为的标准方法。每个带有-ffast-math-style选项的编译器都有属性和编译指示来控制它们，但是这些和选项的精确效果也有所不同。因此，其中某些编译器的不同版本具有[[非常明显]]快速运算的行为，因此，这不仅仅是适当的实用程序集合的问题。获得标准行为的标准方法是使编译器遵循语言标准。-ffast-math及其类似对象主要是为那些不关心浮点数学的细节而只希望其程序（有限且保守地使用FP）运行速度更快的程序员提供的。在任何情况下，-ffast-math的大多数有用作用都可以用精心编写的代码来复制。