问题标签 [fast-math]

简介
Kahan 求和/补偿求和是解决编译器无法尊重数字的关联属性的技术。截断误差导致 (a+b)+c 不完全等于 a+(b+c)，从而在较长的和系列上累积了不希望的相对误差，这是科学计算中的常见障碍。

任务
我希望 Kahan 求和的最佳实施。我怀疑使用手工汇编代码可以实现最佳性能。

尝试
下面的代码使用三种方法计算 [0,1] 范围内的 1000 个随机数的总和。

1. 标准求和：天真的实现，它累积一个增长为 O(sqrt(N)) 的均方根相对误差

2. Kahan summation [g++]：使用 c/c++ 函数“csum”的补偿求和。评论中的解释。请注意，某些编译器可能具有使此实现无效的默认标志（请参见下面的输出）。

3. Kahan summation [asm]：补偿求和实现为“csumasm”，使用与“csum”相同的算法。评论中的神秘解释。

-O3 的输出是：

-O3 -ffast-math 的输出

很明显，-ffast-math 破坏了 Kahan 求和算法，这是不幸的，因为我的程序需要使用 -ffast-math。

问题

1. 是否可以为 Kahan 的补偿求和构造更好/更快的 asm x64 代码？也许有一种巧妙的方法可以跳过一些 movapd 指令？

2. 如果没有更好的 asm 代码，是否有一种 c++ 方法来实现 Kahan 求和，可以与 -ffast-math 一起使用，而无需转向天真的求和？也许 C++ 实现对于编译器的优化通常更灵活。

想法或建议表示赞赏。

更多信息

• “fun”的内容不能内联，但“csum”函数可以。
• 总和必须作为一个迭代过程计算（校正项必须应用于每一个加法）。这是因为预期的求和函数采用取决于先前总和的输入。
• 预期的求和函数被无限期地调用，每秒数亿次，这促使人们追求高性能的低级实现。
• 由于性能原因，诸如 long double、float128 或任意精度库之类的更高精度算术不被视为更高精度的解决方案。

编辑：内联 csum （没有完整代码没有多大意义，但仅供参考）

我正在编写一些 CUDA 代码，我希望它根据是否--use_fast_math设置而表现不同。而且 - 我想在编译时做出决定，而不是在运行时。

NVCC 似乎没有添加或更改预处理器定义--use_fast_math的时间。我通过比较以下输出来检查这一点：

与输出

它们完全一样；所以那条大道似乎被封锁了。现在，如果编译用户会调用 NVCC，--use_fast_math -DUSING_FAST_MATH那么我也可以检测到；但假设它是库代码，我们不能将这些限制强加给用户。

是否有其他方法可以让正在编译的代码注意到它--use_fast_math的开启？

注意：“注意”可能意味着使用预处理器#if或#ifdef指令、使用 SFINAE、使用编译器内置值或 constexpr 函数 - 编译时可用的任何内容。

我不明白 Cuda 中是否存在函数重载。我想在以下两个函数上解释我的问题，我希望能够在 GPU 和 CPU 上都使用它们，并且我不在乎精度：

1. abs分别是哪个函数 cos我应该打电话，为什么？
   • std::abs/ abs/ fabs/ fabsf/anythingelse
   • std::cos/ cos/ cosf/ __cosf/anythingelse

（由于__cosf是 Cuda-intrinsic 和std::abs/std::cos在 Cuda 中不可用，我假设我必须在我的函数中使用预处理器指令来进行这些选择。）

2. 我应该包括哪些标题？

3. 前两个问题的答案是否取决于我是否使用快速数学标志（例如-ffast-math）进行编译。

如果这对答案很重要，我将在Ubuntu 18.04.4下使用nvcc 10.2进行编译。，但我对独立于平台的答案很感兴趣。

我已经阅读了这篇文章和do-denormal-flags-like-denormals-are-zero-daz-affect-comparisons-for-equality，我了解 FTZ 和 DAZ 标志之间的用法和区别。

DAZ 适用于输入，FTZ 适用于 FP 操作的输出。

让我感到困惑的是，如果设置了 FTZ，那么汇编视图中的非规范值来自哪里。我认为它只能是常量值，既可以是直接操作数，也可以是节.rodata（通过 RIP 相对寻址访问）。

但是我在我的二进制文件中发现，这些地方没有异常值，但它仍然存在 FP-ASSIST 问题，导致性能不佳。

如果我同时设置 DAZ 和 FTZ，问题就会消失并且性能会变得更好。实际上，我什至在我的源代码中都没有找到任何非正规输入。我真的很困惑，异常值来自哪里？

另外一个问题，对于指令vmovsd 0x9498(%rip),%xmm0，假设0x9498(%rip)是一个非正规值，xmm0如果我们分别设置 FTZ 或 DAZ，这条指令执行后会发生什么？

据我了解，DAZ 会将其0x9498(%rip)设为零并0移至 xmm0；FTZ 将移动0x9498(%rip)到 xmm0 并发现它是非正规的，因此刷新xmm0为零。我不确定，是否正确？

是否可以修改下面的代码，使其即使在启用快速数学的 GCC 编译时也能正常工作？

注意：我在头文件中有它，但我没有设法关闭头文件的快速数学。请参阅奇怪的 while 循环行为和如何禁用头文件函数的快速数学运算

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

显然，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 编写了一个小测试程序：

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

有五列的 Godbolt 站点

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

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

在其他一些语言中，您可以允许编译器执行一些浮点优化，其中之一是将您的代码重写为代数等效表达式。以下是我的意思的一些例子：

可以优化为

...

可以优化为

显然，编译器默认情况下不会这样做，因为优化的表达式会产生不同的结果，我认为您可以将一个属性分配给方法或类来表示您允许进行这些优化，但我不能找到任何东西。

请参阅这个简单的最小搜索（Godbolt）：

gcc 和 clang 都不会使用-O3. 如果我使用-ffinite-math-only，仍然不会发生自动矢量化。我需要使用-ffinite-math-only和-fno-signed-zeros编译器自动向量化代码。为什么-fno-signed-zeros需要自动矢量化？