迭代可汗求和的优化实现

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

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

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

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

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

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

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
extern "C" void csumasm(double&, double, double&);
__asm__(
"csumasm:n"
"movsd  (%rcx), %xmm0n" //xmm0 = a
"subsd  (%r8), %xmm1n"  //xmm1 - r8 (c) | y = b-c
"movapd %xmm0, %xmm2n"  
"addsd  %xmm1, %xmm2n"  //xmm2 + xmm1 (y) | b = a+y
"movapd %xmm2, %xmm3n" 
"subsd  %xmm0, %xmm3n"  //xmm3 - xmm0 (a) | b - a
"movapd %xmm3, %xmm0n"  
"subsd  %xmm1, %xmm0n"  //xmm0 - xmm1 (y) | - y
"movsd  %xmm0, (%r8)n"  //xmm0 to c
"movsd  %xmm2, (%rcx)n" //b to a
"retn"
);
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)nn",sum1,sum1-sum3);
printf("Kahan compensated summation [g++]:n %.16e (error: %.16e)nn",sum2,sum2-sum3);
printf("Kahan compensated summation [asm]:n %.16en",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"函数可以。
• 总和必须作为迭代过程计算(校正后的项必须应用于每个添加项)。这是因为预期的求和函数采用的输入取决于前一个总和。
• 预期的求和函数被无限期调用，每秒调用数亿次，这促使人们追求高性能的低级实现。
• 由于性能原因，更高精度的算法(如长双精度、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