Лучший способ добавить 3 числа (или 4, или N) в Java - Kahan Sums?

Question

Я нашел совершенно другой ответ на этот вопрос, весь оригинальный вопрос больше не имеет смысла. Однако способ ответа будет полезен, поэтому я немного его модифицирую...

Я хочу подвести итог три double числа, скажем a, b, а также c наиболее устойчивым способом. Я думаю, что использование Kahan Sum было бы хорошим способом.

Однако мне пришла в голову странная мысль: имеет ли смысл:

1. Первая сумма a, b, а также c и помните (абсолютное значение) компенсации.
2. Тогда подведите итог a, c, b
3. Если (абсолютное значение) компенсации второй суммы меньше, используйте эту сумму вместо.
4. Действовать аналогично b, a, c и другие перестановки чисел.
5. Вернуть сумму с наименьшей связанной абсолютной компенсацией.

Получу ли я таким образом более "стабильное" добавление трех чисел? Или порядок чисел в сумме не оказывает (пригодного для использования) влияния на компенсацию, оставленную в конце Суммирования? Имея (пригодный для использования), я хочу спросить, достаточно ли само значение компенсации достаточно, чтобы содержать информацию, которую я могу использовать?

(Я использую язык программирования Java, хотя я думаю, что здесь это не имеет значения.)

Большое спасибо, Томас.

Answer 1

Я думаю, что нашел гораздо более надежный способ решения проблемы "Добавить 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 эта проблема имеет было решено довольно приятно. Чтобы увидеть полный код, я загрузил исходный код Python (Python 3.5.1rc1 - 2015-11-23) с https://www.python.org/getit/source/, где мы можем найти следующий метод (в разделе 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;
}

Этот метод суммирования отличается от метода Кахана, он использует переменное количество переменных компенсации. При добавлении iне более i дополнительные переменные компенсации (хранятся в массиве p) привыкнуть. Это означает, что если я хочу добавить 3 числа, мне могут понадобиться 3 дополнительные переменные. Для 4 чисел мне могут понадобиться 4 дополнительные переменные. Поскольку число используемых переменных может увеличиться с n в n+1 только после nЭто слагаемое загружено, я могу перевести приведенный выше код в 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 – 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 и сложите элементы, одновременно используя его для хранения переменных компенсации. Так как мы загружаем summand по указателю i прежде чем элемент массива с указанным индексом может быть использован для компенсации, это будет работать.

Поскольку массив будет небольшим, если количество добавляемых переменных невелико и не выйдет за рамки нашего метода, я думаю, что есть неплохой шанс, что он будет размещен непосредственно в стеке JIT, что может сделать код довольно быстро.

Я признаю, что не до конца понимал, почему авторы исходного кода обрабатывали бесконечности, переполнения и NaN так, как они это делали. Здесь мой код отклоняется от оригинала. (Надеюсь, я не испортил это.)

В любом случае, теперь я могу подвести итог 3, 4 или ndouble номера, выполнив:

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 может представлять только longв -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 в результате.