简介
Kahan 求和/补偿求和是解决编译器无法尊重数字的关联属性的技术。截断误差导致 (a+b)+c 不完全等于 a+(b+c)，从而在较长的和系列上累积了不希望的相对误差，这是科学计算中的常见障碍。
任务
我希望 Kahan 求和的最佳实施。我怀疑使用手工汇编代码可以实现最佳性能。
尝试
下面的代码使用三种方法计算 [0,1] 范围内的 1000 个随机数的总和。
标准求和：天真的实现，它累积一个增长为 O(sqrt(N)) 的均方根相对误差
Kahan summation [g++]：使用 c/c++ 函数“csum”的补偿求和。评论中的解释。请注意，某些编译器可能具有使此实现无效的默认标志（请参见下面的输出）。
Kahan summation [asm]：补偿求和实现为“csumasm”，使用与“csum”相同的算法。评论中的神秘解释。
#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