迭代卡汉求和的优化实现

Question

介绍
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]：使用与“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。

题

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

2. 如果没有更好的 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表示优化级编译指示和属性“不适合生产代码”。

更新：-O3显然问题的一部分是基于不安全的误解，还是其他什么？这是; ISO C++ 指定了类似于 GCC 的 FP 数学规则-fno-fast-math。使用 just 编译所有内容-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相同的性能（1 uop，3 到 5 个周期延迟，具体取决于微架构： https: //agner.org/optimize/）。addsdsubsd

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

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

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

定义一个实际上通过引用获取参数的独立函数看起来非常影响性能。（特别是当它甚至不返回其中任何一个作为返回值以避免存储/重新加载链之一时）。即使只是进行函数调用也会通过破坏许多寄存器来引入大量开销。（在 Windows x64 中并不像在 x86-64 System V 中那么糟糕，其中所有XMM 寄存器都被调用破坏，并且更多的整数寄存器。）

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

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

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

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

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

但最好通过将包含求和循环的整个函数放入单独的文件中来禁用快速数学。您可能会根据使用快速数学编译一些代码和其他不使用快速数学的代码来选择非内联函数调用边界需要在哪里。

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

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

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

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