浮点数与 GCC、Intel 和 MS C 的双倍性能

问：

我已经看到了之前古老的 x87 时代的线程，并认为是时候在 SSE2 和 AVX 的现代时代重新审视它了。我相当简单的 C 基准测试的结果大多符合我的预期，但也有一些惊喜。我很想知道在基于 Unix 的 GCC 上进行此基准测试的正确时间，而不是使用 MinGW FP 库。

我选择了一个相对密集的复杂“代表性”数字代码段（对不起，它的长度尽可能短）。它是求解开普勒方程的三次启动器，它恰好使用了大多数常用函数（可以修改为使用带有 DIY 立方根的 exp 和 log，但我已将其精简到最低限度，仍然适用于浮点和双精度）。

这些是每个编译器上双精度 64 位实数和 32 位浮点数的机器周期计时。编译器是 Intel 2023、MSC 2022 和 GCC 13.1.0（MinGW 端口），除了进行 x87 测试外，都是 x64。我的 GCC 结果是狡猾的，因为 Mingw FP 库对 sin、cos、exp、log 和 atan2 使用 80 位 x87 模式。我已经通过对这些 x87 指令的开销的最佳猜测估计来纠正它们。我也可能没有使用最佳编译器设置，因为我是 GCC 的新手。它们是在 i5-12600 CPU OS Win 11 pro 64 上完成的。

gcc -O3 -Ofast -march=native -mavx -mfpmath=sse benchbasic.cpp

令我惊讶的是它们的相似程度——只有在 MS 编译器上，性能才有明显的提高，从双倍变为浮点。第二个观察结果是，英特尔的结果快得令人怀疑。这是因为它已经找到了一种方法来完全矢量化问题，并一次将其作为 2 或 4 个连续的 M 值来执行。我知道 GCC 也可以从其他测试中做到这一点，但如果涉及 if 语句，则不会。MSC 将向量 FP 寄存器视为标量，因此 SSE2 可能是那里最好的选择。传统的 x87 浮点 32 位代码似乎在英特尔和 MS 上编译速度变慢 - 不知道为什么。违反直觉！

英特尔还通过在 x87 模式下使用更快的 SSE2 库代码进行 trig 作弊。它还将 pow（x，1.0/3） 换成 cbrt（x） 和其他可爱的转换 （-Ofast）。其他两个编译器都直接使用 x87 操作码。在 -O3 中，所有编译器都内联了大部分代码，因为它是一次性使用的，但如果我采取措施防止这种情况发生，结果几乎相同。x87 80 位浮点数（在 GCC 和 Intel 上设置了某些选项的长双精度）最好避免，但有时对于累积双精度的总和很有用。

在整理它时，我还发现了另一个好奇心，以便时间在 80 列屏幕上很好地排列。从 printf 中的字符串中删除一些空格从根本上改变了 MS 编译代码在 Microsoft VC 调试器中运行时的计时从 250 增加到 ~1000。我真的很想知道为什么会这样！它似乎在一段时间内是可重现的（但事实并非如此）——如果从命令行运行，它总是全速运行。

静态数据长度的执行时间存在较小的系统变化，我怀疑这是由于 MSC 将编译器生成的长 AVX 结构放在堆栈上，有时甚至跨越缓存行边界。

这是代码 - 它应该在大多数系统上未经修改地编译。它测量的是吞吐量，而不是当前配置的延迟（与我的问题相关的指标）。

// BenchBasic.cpp : This is a simple benchmark. float vs double
// define symbol ALLDOUBLE to switch to doubles

#include <stdio.h>
#define _USE_MATH_DEFINES // for C
#include <math.h>
#include <intrin.h>
#include <stdint.h>
#include <time.h> 

//#define ALLDOUBLE

#ifdef ALLDOUBLE
#define float double
#define sinf(x) sin(x)
#define sqrtf(x) sqrt(x)
#define cbrtf(x) cbrt(x)
#define atan2f(x,y) atan2(x,y)
#define fabsf(x) fabs(x)
#endif 


double QuickTime(double (*func)(double, double), const char* name)
{ 
  double dM, E, M;
  uint64_t start, end;
  unsigned int aux;
  time_t  start_t, end_t;
  int j, k, cycles;
  const int MDIV = 10000000;
  end = end_t = 0; // to stop compiler warnings
  E = M = 1e-10;   // some routines don't like 0 as an input
  dM = (M_PI - M) / (MDIV - 1);
  printf("\n%10s %i", name, 0); 
//  printf("\n%10s", name); // uncomment this line for 4x faster in MSVC debug!
  for (j = 0; j < 39; j++)
  {
    E = 0;
    end_t = clock();
    while (end_t == (start_t = clock()))
        ; // sync start tick
    start = __rdtscp(&aux);
    for (k = 0; k < MDIV; k++)
    {
        E += (*func)(M, 0.7);
        M += dM;
    }
    end = __rdtscp(&aux);
    if (E == 42) printf("\nI hate global optimisers\n"); 
    // side effect added to defeat global optimisers nulling out code
    end = end - start;  // CPU cycles approx 
    end_t = clock() - start_t;
    cycles = (int)((end + MDIV / 2) / MDIV);
    printf(" %5i %5.3f", cycles, (double)end_t / CLOCKS_PER_SEC);
  }
  return E;
}

double BasicF32(double ed, double Md)
{
// this 32 bit implementation uses cubic solver inline NB a,b divided by 3 at source
// Also has problems with some awkward values like e = 0.85, M = 0.70
// Cubic solver can overflow since there is not enough headroom on the 
// exponent max(tan(E)) ~ 10^7 and (10^7)^6 = 10^42 > 10^38
// Cubic solver needs additional work to protect against overflows

 float  a, b, c, d, e2, e, E, f3, M, q, r, r2, s, s2, t;
 if (Md == 0) return Md;
 M = (float)Md;
 e = (float)ed;
 e2 = e * e;
 s = sinf(M / 2);
 s2 = s * s;
 s2 += s2;
 s = sinf(M);
 f3 = 60 * (1 - e) + e2 * (3 + 7 * e) - (60 + 3 * e2) * s2;
 if (f3 == 0)   f3 = 1e-6f;  // defend against f3 == 0 when M~=pi/2-e
 a = -(20 + e2) * s / f3;    // a/3
 b = 20 * (1 - e - s2) / f3; // b/3
 c = -30 * s / f3;           // c/2
// start of cubic solver
 q = b - a * a;
 r = 1.5f * a * b - c - a * a * a;
 r2 = r * r;
 d = q * q * q + r2;
 if (d < 0) d = 0; // these problems should only have one real root.
 if (r2 < d * 0.0001f)
 { // avoid catastrophic canellation of roots
    s = 2 * r / (3 * cbrtf(d));
    t = r * r / (27 * d);
    t = s * (1 + t*(5 + 66 * t)) - a; // next highest order
 }
 else
 {
    d = sqrtf(d);
    if (r > 0)
        s = cbrtf(r + d);
    else
        s = -cbrtf(d - r);
    if (s != 0) t = s - q / s - a;
    else t = -a;  // special case arises when e=1 and M=0 or pi
 }
 if (t * M < 0) a = -1;  else a = 1;
 E = atan2f(a * t, a);
 if (fabsf(E) < fabsf(M)) E = M; // defensive fix for M~=pi rounding errors
 return E;
}

int main()
{
  QuickTime(BasicF32, "Basic");
}

我完全有可能在某个地方滑倒了，所以知道它在其他一些高端 CPU 上的表现会很有趣。浮点 32 位代码中存在一个已知问题，其中 tan 超过 3e6 时，影响不到 M 范围的 0.1%。我不认为它过度改变了时间（但我可能是错的）。

其他人对现代 SSE2/AVX 硬件上浮点数与双打数的优缺点有何体验？我原以为浮子会更快，做更低精度的三角函数等。因此，我对这些结果感到非常惊讶。