在 Java 中添加 3 个数字（或 4，或 N）的最佳方法 - Kahan 总和？

Question

我找到了这个问题的完全不同的答案，整个原始问题不再有意义了。不过回答的方式还是有用的，所以我稍微修改了一下...

我想总结三个double数，比如a、b和c，在最数值稳定的方式成为可能。我认为使用 Kahan Sum 是可行的方法。

然而，我突然想到一个奇怪的想法：这样做是否有意义：

先总结一下a、b、c，记住补偿的（绝对值）
然后总结a,c,b
如果第二个金额的补偿（的绝对值）较小，则用这个金额代替。
与b、a、c和其他数字排列类似。
Return 具有最小相关绝对补偿的总和。

这样我会得到更多"stable"三个数字的加法吗？还是求和中数字的顺序对求和结束时留下的补偿没有（可用的）影响？ With (use-able) 我的意思是问补偿值本身是否足够稳定，可以包含我可以使用的信息？

（我使用的是Java编程语言，虽然我认为这在这里并不重要。）

非常感谢，托马斯.

Answer 1

我想我找到了一种更可靠的方法来解决 "Add 3"（或 "Add 4" 或 "Add N" 数字问题。

首先，我实现了我原来的想法post。它产生了相当大的代码，最初似乎可以工作。但是，它在以下情况下失败了：添加 Double.MAX_VALUE、1 和 -Double.MAX_VALUE。结果是 0.

@njuffa 的评论启发了我更深入地挖掘 http://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/, I found that in Python, this problem has been solved quite nicely. To see the full code, I downloaded the Python source (Python 3.5.1rc1 - 2015-11-23) from https://www.python.org/getit/source/，在那里我们可以找到以下方法（在 PYTHON SOFTWARE FOUNDATION LICENSE VERSION 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;
}

这种求和方法与卡汉的方法不同，它使用可变数量的补偿变量。添加第 i 个数字时，最多使用 i 个额外的补偿变量（存储在数组 p 中）。这意味着如果我想添加 3 个数字，我可能需要 3 个额外的变量。对于 4 个数字，我可能需要 4 个额外的变量。由于仅在加载第 n 个被加数后，使用的变量数量可能会从 n 增加到 n+1，因此我可以将上面的代码转换为 Java 如下：

/**
 * 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 &ndash; 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 并将元素相加，同时使用它来存储补偿变量。由于我们在索引 i 之前加载 summand，因此在 处的数组元素可能会用于补偿，这将起作用。

由于如果要添加的变量数量少并且不会超出我们方法的范围，数组就会很小，我认为它很有可能被 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 只能表示的事实-9007199254740992..9007199254740992L 中的 long 秒。但这可以很容易地通过将每个 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 作为结果。

在 Java 中添加 3 个数字（或 4，或 N）的最佳方法 - Kahan 总和？

Best Way to Add 3 Numbers (or 4, or N) in Java - Kahan Sums?

precision

double

numerical

numerical-methods

numerical-stability