c++ - 迭代Kahan求和的最佳实现

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

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

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

标准求和:天真的实现，累积了均方根相对误差，误差随着O(sqrt(N))的增长而增加

Kahan求和[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的补偿总和构造更好/更快的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

最佳答案

您可以将不需要使用-ffast-math的函数(如csum循环)放在一个没有-ffast-math进行编译的单独文件中。

可能您也可以使用__attribute__((optimize("no-fast-math")))，但是https://gcc.gnu.org/onlinedocs/gcc/Common-Function-Attributes.html表示不幸的是，优化级别的编译指示和属性“不适合生产代码”。

更新:显然是问题的一部分，原因是人们对-O3不安全或存在某种误解。它是; ISO C++指定了类似于GCC的-fno-fast-math的FP数学规则。仅使用-O3编译所有内容显然可以使OP的代码快速安全地运行。有关诸如OpenMP之类的变通方法，请参见此答案的底部，以在不实际启用-ffast-math的情况下获得代码某些部分的快速计算的好处。

ICC默认为快速路径，因此您必须专门启用FP = strict才能使用-O3使其安全，但是gcc/clang默认为完全严格的FP，而不考虑其他优化设置。 (-Ofast = -O3 -ffast-math除外)

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

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

因此，您实际上是在并行执行8个补偿求和，每个求和得到第8个输入元素。

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

如果要使用内联汇编，最好以内联方式实际使用它，这样您就可以在带有扩展汇编的XMM寄存器中获得多个输入和多个输出，而不是通过内存进行存储/重新加载。

定义一个实际上通过引用来接受args的独立函数看起来会严重影响性能。 (特别是当它甚至没有返回任何一个作为返回值时，避免使用存储/重载链之一)。即使只是进行函数调用，也会浪费大量寄存器，从而造成大量开销。 (在Windows x64中，不如在x86-64 System V中糟糕透了，所有XMM reg都被调用，以及更多的整数reg。)

同样，您的独立函数特定于Windows x64调用约定，因此它比函数内部的内联汇编的移植性差。

顺便说一句， clang设法仅使用两个csum(double&, double, double&):指令而不是您的asm中的3个(我假设您从GCC的asm输出中复制了)来实现movapd。 https://godbolt.org/z/lw6tug。如果可以假定AVX可用，则可以避免任何情况。

顺便说一句，movaps小1字节，应改为使用。没有CPU拥有double和float的单独数据域/转发网络，只有vec-FP与vec-int(相对于GP整数)

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

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

但是最好将包含求和循环的整个函数禁用快速运算，方法是将其放入单独的文件中。您可能会因为编译一些带有快速运算符的代码而没有编译一些带有快速运算符的代码，而难以选择非内联函数调用边界的位置。

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

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

或者也许为您的Kahan求和文件尝试-ftrapping-math，但这是自动向量化涉及归约法的FP循环所需要的主要方法，并且在其他情况下会有所帮助。但是，您仍然可以获得其他一些有值(value)的优化。

另一种获得通常需要快速运算的优化的方法是-ffast-math -fno-associative-math ，即使在未经自动矢量化编译的文件中也可以使用OpenMP进行自动矢量化。您可以声明一个用于减少的累加器变量，以使gcc对其上的操作进行重新排序，就好像它们是关联的一样。

关于c++ - 迭代Kahan求和的最佳实现，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/57301596/