迭代Kahan求和的优化实现

Question

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

任务
我希望 Kahan 求和的最佳实现。我怀疑手工编写的汇编代码可能会达到最佳性能。

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

标准求和：朴素的实现累积了一个增长为 O(sqrt(N))
[= 的均方根相对误差71=]
卡汉求和[g++]：使用c/c++函数"csum"进行补偿求和。评论中的解释。请注意，某些编译器可能具有使此实现无效的默认标志（请参阅下面的输出）。
Kahan 求和 [asm]：使用与 "csum" 相同的算法实现为 "csumasm" 的补偿求和。评论中的神秘解释。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

extern "C" void csumasm(double&, double, double&);
__asm__(
    "csumasm:\n"
    "movsd  (%rcx), %xmm0\n" //xmm0 = a
    "subsd  (%r8), %xmm1\n"  //xmm1 - r8 (c) | y = b-c
    "movapd %xmm0, %xmm2\n"  
    "addsd  %xmm1, %xmm2\n"  //xmm2 + xmm1 (y) | b = a+y
    "movapd %xmm2, %xmm3\n" 
    "subsd  %xmm0, %xmm3\n"  //xmm3 - xmm0 (a) | b - a
    "movapd %xmm3, %xmm0\n"  
    "subsd  %xmm1, %xmm0\n"  //xmm0 - xmm1 (y) | - y
    "movsd  %xmm0, (%r8)\n"  //xmm0 to c
    "movsd  %xmm2, (%rcx)\n" //b to a
    "ret\n"
);

void csum(double &a,double b,double &c) { //this function adds a and b, and passes c as a compensation term
    double y = b-c; //y is the correction of b argument
    b = a+y; //add corrected b argument to a argument. The output of the current summation
    c = (b-a)-y; //find new error to be passed as a compensation term
    a = b;
}

double fun(double fMin, double fMax){
    double f = (double)rand()/RAND_MAX;
    return fMin + f*(fMax - fMin); //returns random value
}

int main(int argc, char** argv) {
    int N = 1000;

    srand(0); //use 0 seed for each method
    double sum1 = 0;
    for (int n = 0; n < N; ++n)
        sum1 += fun(0,1);

    srand(0);
    double sum2 = 0;
    double c = 0; //compensation term
    for (int n = 0; n < N; ++n)
        csum(sum2,fun(0,1),c);

    srand(0);
    double sum3 = 0;
    c = 0;
    for (int n = 0; n < N; ++n)
        csumasm(sum3,fun(0,1),c);

    printf("Standard summation:\n %.16e (error: %.16e)\n\n",sum1,sum1-sum3);
    printf("Kahan compensated summation [g++]:\n %.16e (error: %.16e)\n\n",sum2,sum2-sum3);
    printf("Kahan compensated summation [asm]:\n %.16e\n",sum3);
    return 0;
}

-O3 的输出是：

Standard summation:
 5.1991955320902093e+002 (error: -3.4106051316484809e-013)

Kahan compensated summation [g++]:
 5.1991955320902127e+002 (error: 0.0000000000000000e+000)

Kahan compensated summation [asm]:
 5.1991955320902127e+002

带-O3 -ffast-math的输出

Standard summation:
 5.1991955320902093e+002 (error: -3.4106051316484809e-013)

Kahan compensated summation [g++]:
 5.1991955320902093e+002 (error: -3.4106051316484809e-013)

Kahan compensated summation [asm]:
 5.1991955320902127e+002

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

问题

是否可以为 Kahan 的补偿求和构建 better/faster asm x64 代码？也许有一种聪明的方法可以跳过一些 movapd 指令？
如果没有更好的 asm 代码，是否有一种 c++ 方法来实现 Kahan 求和，可以与 -ffast-math 一起使用，而无需转变成简单的求和？也许 C++ 实现通常更灵活，以便编译器进行优化。

欢迎提出想法或建议。

更多信息

无法内联 "fun" 的内容，但可以内联 "csum" 函数。
总和必须作为迭代过程计算（更正后的项必须应用于每次加法）。这是因为预期的求和函数采用取决于先前总和的输入。
预期的求和函数被无限期调用，每秒调用数亿次，这激发了对高性能低级实现的追求。
由于性能原因，long double、float128 或任意精度库等更高精度的算法不被视为更高精度的解决方案。

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

        subsd   xmm0, QWORD PTR [rsp+32]
        movapd  xmm1, xmm3
        addsd   xmm3, xmm0
        movsd   QWORD PTR [rsp+16], xmm3
        subsd   xmm3, xmm1
        movapd  xmm1, xmm3
        subsd   xmm1, xmm0
        movsd   QWORD PTR [rsp+32], xmm1

Answer 1

您可以将需要而非使用 -ffast-math 的函数（如 csum 循环）放在单独的文件中，该文件无需 -ffast-math 即可编译。

也许您也可以使用 __attribute__((optimize("no-fast-math")))，但不幸的是 https://gcc.gnu.org/onlinedocs/gcc/Common-Function-Attributes.html 说优化级别的编译指示和属性不是 "suitable in production code"。

更新： 显然问题的一部分是基于对 -O3 不安全的误解，还是什么？这是; ISO C++ 指定类似于 GCC 的 -fno-fast-math 的 FP 数学规则。仅使用 -O3 编译所有内容显然可以使 OP 的代码运行快速且安全。请参阅此答案的底部，了解 OpenMP 等解决方法，以便在不实际启用 -ffast-math.

的情况下为代码的某些部分获得快速数学的一些好处

ICC 默认为快速路径，因此您必须特别启用 FP=strict 以使其与 -O3 一起安全，但 gcc/clang 默认为完全严格的 FP，而不管其他优化设置如何。（除了 -Ofast = -O3 -ffast-math）

您应该能够通过保持一个（或四个）总计向量和相等数量的补偿向量来向量化 Kahan 求和。您可以使用内在函数来做到这一点（只要您不为该文件启用快速数学）。

例如使用 SSE2 __m128d 为每条指令进行 2 次打包加法。或 AVX __m256d。在现代 x86 上，addpd / subpd 具有与 addsd 和 subsd 相同的性能（1 uop，3 到 5 个周期延迟，具体取决于微体系结构：https://agner.org/optimize/） .

所以你实际上是在并行地进行 8 次补偿求和，每次求和得到每 8 个输入元素。

使用 fun() 即时生成随机数比从内存中读取随机数要慢得多。如果您的正常用例在内存中有数据，您应该对其进行基准测试。否则我想标量很有趣。

如果您打算使用内联 asm，实际内联使用它会好得多，这样您就可以通过扩展 asm 在 XMM 寄存器中获得多个输入和多个输出，而不是 stored/reloaded 通过内存。

定义一个实际上通过引用获取 args 的独立函数看起来很不利于性能。（特别是当它甚至没有 return 它们中的任何一个作为 return 值以避免 store/reload 链之一时）。即使只是进行函数调用也会破坏许多寄存器，从而引入大量开销。（在 Windows x64 中并不像在 x86-64 系统 V 中那么糟糕，其中所有 XMM regs 被调用破坏，以及更多的整数 regs。）

此外，您的独立函数特定于 Windows x64 调用约定，因此它的可移植性不如函数内部的内联汇编。

顺便说一句，clang 仅用两条 movapd 指令 就成功实现了 csum(double&, double, double&):，而不是你的 asm 中的 3 条指令（这我假设您是从 GCC 的 asm 输出中复制的）。 https://godbolt.org/z/lw6tug。如果您可以假设 AVX 可用，则可以避免任何。

顺便说一句，movaps 小了 1 个字节，应该改用它。 double vs. float 没有 CPU 有单独的数据域/转发网络，只有 vec-FP vs. vec-int（vs. GP 整数）

但实际上到目前为止，您的赌注是让 GCC 在没有 -ffast-math 的情况下编译文件或函数。 https://gcc.gnu.org/wiki/DontUseInlineAsm。这让编译器在 AVX 可用时避免 movaps 指令，除了让它在展开时更好地优化。

如果您愿意为每个元素接受函数调用的开销，您不妨让编译器通过将 csum 放在一个单独的文件中来生成该 asm。（希望 link-time 优化尊重 -fno-fast-math 一个文件，也许不内联该函数。）

但是，通过将 it 放在一个单独的文件中，为包含求和循环的整个函数禁用快速数学会好得多。您可能会在选择需要非内联函数调用边界的位置时遇到困难，这取决于使用快速数学编译一些代码而其他代码没有。

理想情况下使用 -O3 -march=native 和配置文件引导优化编译您的所有代码。还有 -flto link-时间优化以启用跨文件内联。

-ffast-math 打破 Kahan 求和并不奇怪： 将 FP 数学视为关联是使用快速数学的主要原因之一。如果您需要 -ffast-math 的其他部分，例如 -fno-math-errno 和 -fno-trapping-math，以便数学函数可以更好地内联，请手动启用它们。这些基本上总是安全的，也是个好主意；没有人在调用 sqrt 后检查 errno ，因此为某些输入设置 errno 的要求只是 C 的一个可怕的错误设计，它给实现带来了不必要的负担。 GCC 的 -ftrapping-math 默认情况下处于打开状态，即使它已损坏（它并不总是准确地重现如果你揭露任何掩码你会得到的 FP 异常的数量）所以 it should really be off by default。关闭它不会启用任何会破坏 NaN 传播的优化，它只会告诉 GCC 异常数量不是可见的副作用。

或者尝试 -ffast-math -fno-associative-math 用于您的 Kahan 求和文件，但这是自动矢量化涉及归约的 FP 循环所需的主要文件，并且在其他情况下也有帮助。但是，您仍然可以获得其他一些有价值的优化。

获得通常需要快速数学的优化的另一种方法是 #pragma omp simd 以启用使用 OpenMP 的自动矢量化，即使在编译的文件中没有自动矢量化。您可以为归约声明一个累加器变量，让 gcc 对其操作重新排序，就好像它们是关联的一样。

迭代Kahan求和的优化实现

Optimal implementation of iterative Kahan summation

c++

sse

x86-64

inline-assembly

fast-math