我对这个问题找到了一个完全不同的答案，整个原始问题不再有意义。但是，答案方式很有用，所以我对其进行了一些修改...
我想以数字上最稳定的方式总结三个double数字，比如a、b和。我认为使用Kahan Sum是可行的方法。c
然而，我想到了一个奇怪的想法：这样做有意义吗：
a, b,c并记住补偿的（绝对值）。a，，，cbb、a和c其他数字排列。这样我会得到一个更“稳定”的三个数字相加吗？或者总和中的数字顺序对总和结束时留下的补偿没有（可用）影响？使用（可用）我的意思是问补偿值本身是否足够稳定以包含我可以使用的信息？
（我使用的是 Java 编程语言，尽管我认为这在这里无关紧要。）
非常感谢，托马斯。
我想我找到了一种更可靠的方法来解决“加 3”（或“加 4”或“加 N”数字问题。
首先，我从原始帖子中实现了我的想法。它产生了相当大的代码，最初似乎可以工作。但是，它在以下情况下失败：添加Double.MAX_VALUE、、1和-Double.MAX_VALUE。结果是0。
@njuffa 的评论启发了我深入挖掘，在http://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/，我发现在 Python 中，这个问题有已经很好地解决了。为了查看完整代码，我从https://www.python.org/getit/source/下载了 Python 源代码（Python 3.5.1rc1 - 2015-11-23），我们可以在其中找到以下方法（在 PYTHON SOFTWARE 下基金会许可版本 2)：
static PyObject*
math_fsum(PyObject *self, PyObject *seq)
{
PyObject *item, *iter, *sum = NULL;
Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
double x, y, t, ps[NUM_PARTIALS], *p = ps;
double xsave, special_sum = 0.0, inf_sum = 0.0;
volatile double hi, yr, lo;
iter = PyObject_GetIter(seq);
if (iter == NULL)
return NULL;
PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
for(;;) { /* for x in iterable */
assert(0 <= n && n <= m);
assert((m == NUM_PARTIALS && p == ps) ||
(m > NUM_PARTIALS && p != NULL));
item = PyIter_Next(iter);
if (item == NULL) {
if (PyErr_Occurred())
goto _fsum_error;
break;
}
x = PyFloat_AsDouble(item);
Py_DECREF(item);
if (PyErr_Occurred())
goto _fsum_error;
xsave = x;
for (i = j = 0; j < n; j++) { /* for y in partials */
y = p[j];
if (fabs(x) < fabs(y)) {
t = x; x = y; y = t;
}
hi = x + y;
yr = hi - x;
lo = y - yr;
if (lo != 0.0)
p[i++] = lo;
x = hi;
}
n = i; /* ps[i:] = [x] */
if (x != 0.0) {
if (! Py_IS_FINITE(x)) {
/* a nonfinite x could arise either as
a result of intermediate overflow, or
as a result of a nan or inf in the
summands */
if (Py_IS_FINITE(xsave)) {
PyErr_SetString(PyExc_OverflowError,
"intermediate overflow in fsum");
goto _fsum_error;
}
if (Py_IS_INFINITY(xsave))
inf_sum += xsave;
special_sum += xsave;
/* reset partials */
n = 0;
}
else if (n >= m && _fsum_realloc(&p, n, ps, &m))
goto _fsum_error;
else
p[n++] = x;
}
}
if (special_sum != 0.0) {
if (Py_IS_NAN(inf_sum))
PyErr_SetString(PyExc_ValueError,
"-inf + inf in fsum");
else
sum = PyFloat_FromDouble(special_sum);
goto _fsum_error;
}
hi = 0.0;
if (n > 0) {
hi = p[--n];
/* sum_exact(ps, hi) from the top, stop when the sum becomes
inexact. */
while (n > 0) {
x = hi;
y = p[--n];
assert(fabs(y) < fabs(x));
hi = x + y;
yr = hi - x;
lo = y - yr;
if (lo != 0.0)
break;
}
/* Make half-even rounding work across multiple partials.
Needed so that sum([1e-16, 1, 1e16]) will round-up the last
digit to two instead of down to zero (the 1e-16 makes the 1
slightly closer to two). With a potential 1 ULP rounding
error fixed-up, math.fsum() can guarantee commutativity. */
if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
(lo > 0.0 && p[n-1] > 0.0))) {
y = lo * 2.0;
x = hi + y;
yr = x - hi;
if (y == yr)
hi = x;
}
}
sum = PyFloat_FromDouble(hi);
_fsum_error:
PyFPE_END_PROTECT(hi)
Py_DECREF(iter);
if (p != ps)
PyMem_Free(p);
return sum;
}
这种求和方法不同于 Kahan 的方法，它使用可变数量的补偿变量。添加第ith 个数字时，最多使用i额外的补偿变量（存储在数组中p）。这意味着如果我想添加 3 个数字，我可能需要 3 个额外的变量。对于 4 个数字，我可能需要 4 个附加变量。由于仅在加载第 th 和之后使用的变量的数量可能会增加，因此我可以将上述代码转换为 Java，如下所示n：n+1n
/**
* Compute the exact sum of the values in the given array
* {@code summands} while destroying the contents of said array.
*
* @param summands
* the summand array – will be summed up and destroyed
* @return the accurate sum of the elements of {@code summands}
*/
private static final double __destructiveSum(final double[] summands) {
int i, j, n;
double x, y, t, xsave, hi, yr, lo;
boolean ninf, pinf;
n = 0;
lo = 0d;
ninf = pinf = false;
for (double summand : summands) {
xsave = summand;
for (i = j = 0; j < n; j++) {
y = summands[j];
if (Math.abs(summand) < Math.abs(y)) {
t = summand;
summand = y;
y = t;
}
hi = summand + y;
yr = hi - summand;
lo = y - yr;
if (lo != 0.0) {
summands[i++] = lo;
}
summand = hi;
}
n = i; /* ps[i:] = [summand] */
if (summand != 0d) {
if ((summand > Double.NEGATIVE_INFINITY)
&& (summand < Double.POSITIVE_INFINITY)) {
summands[n++] = summand;// all finite, good, continue
} else {
if (xsave <= Double.NEGATIVE_INFINITY) {
if (pinf) {
return Double.NaN;
}
ninf = true;
} else {
if (xsave >= Double.POSITIVE_INFINITY) {
if (ninf) {
return Double.NaN;
}
pinf = true;
} else {
return Double.NaN;
}
}
n = 0;
}
}
}
if (pinf) {
return Double.POSITIVE_INFINITY;
}
if (ninf) {
return Double.NEGATIVE_INFINITY;
}
hi = 0d;
if (n > 0) {
hi = summands[--n];
/*
* sum_exact(ps, hi) from the top, stop when the sum becomes inexact.
*/
while (n > 0) {
x = hi;
y = summands[--n];
hi = x + y;
yr = hi - x;
lo = y - yr;
if (lo != 0d) {
break;
}
}
/*
* Make half-even rounding work across multiple partials. Needed so
* that sum([1e-16, 1, 1e16]) will round-up the last digit to two
* instead of down to zero (the 1e-16 makes the 1 slightly closer to
* two). With a potential 1 ULP rounding error fixed-up, math.fsum()
* can guarantee commutativity.
*/
if ((n > 0) && (((lo < 0d) && (summands[n - 1] < 0d)) || //
((lo > 0d) && (summands[n - 1] > 0d)))) {
y = lo * 2d;
x = hi + y;
yr = x - hi;
if (y == yr) {
hi = x;
}
}
}
return hi;
}
该函数将获取数组summands并将元素相加，同时使用它来存储补偿变量。由于我们在所述索引处的数组元素可能用于补偿之前加载summandat 索引，因此这将起作用。i
由于如果要添加的变量数量很少并且不会超出我们方法的范围，那么数组会很小，我认为它很有可能会被 JIT 直接分配到堆栈上，这可能会使代码相当快。
我承认我并不完全理解为什么原始代码的作者会以他们的方式处理无穷大、溢出和 NaN。在这里，我的代码与原始代码不同。（我希望我没有搞砸。）
无论哪种方式，我现在都可以通过以下方式总结 3、4 或n double数字：
public static final double add3(final double x0, final double x1,
final double x2) {
return __destructiveSum(new double[] { x0, x1, x2 });
}
public static final double add4(final double x0, final double x1,
final double x2, final double x3) {
return __destructiveSum(new double[] { x0, x1, x2, x3 });
}
如果我想将 3 或 4 个long数字相加并获得精确的结果double，我将不得不处理doubles 只能表示longs in的事实-9007199254740992..9007199254740992L。但这可以通过将每个long部分分成两部分来轻松完成：
public static final long add3(final long x0, final long x1,
final long x2) {
double lx;
return __destructiveSum(new long[] {new double[] { //
lx = x0, //
(x0 - ((long) lx)), //
lx = x1, //
(x1 - ((long) lx)), //
lx = x2, //
(x2 - ((long) lx)), //
});
}
public static final long add4(final long x0, final long x1,
final long x2, final long x3) {
double lx;
return __destructiveSum(new long[] {new double[] { //
lx = x0, //
(x0 - ((long) lx)), //
lx = x1, //
(x1 - ((long) lx)), //
lx = x2, //
(x2 - ((long) lx)), //
lx = x3, //
(x3 - ((long) lx)), //
});
}
我认为这应该是正确的。至少我现在可以添加Double.MAX_VALUE,1和-Double.MAX_VALUE并得到1结果。