Meilleur moyen d'ajouter 3 numéros (ou 4, ou N) en Java - Kahan Sums?

Question

J'ai trouvé une réponse complètement différente à cette question, toute la question originale n'a plus de sens. Cependant, la manière de réponse utile, donc je modifier un peu ...

Je veux résumer trois double chiffres, disent a, b et c, dans le plus façon stable numériquement possible. Je pense que l'utilisation d'un Kahan Sum serait le chemin à parcourir.

Cependant, une étrange pensée est produite pour moi: Ne serait-il logique de:

1. somme Première place a, b et c et rappelez-vous la (valeur absolue de la) compensation.
2. résumé ensuite a, c, b
3. Si la (valeur absolue de la) compensation de la seconde somme est plus petite, utilisez cette somme à la place.
4. Procéder de manière similaire avec b, a, c et d'autres permutations des nombres.
5. Renvoie la somme avec la plus petite compensation absolue associée.

Est-ce que j'aurais une addition plus "stable" de trois nombres de cette façon? Ou l'ordre des nombres dans la somme n'a-t-il pas d'impact (utilisable) sur la compensation restante à la fin de la Récapitulation? Avec (utilisable), je veux demander si la valeur de compensation elle-même est assez stable pour contenir des informations que je peux utiliser?

(J'utilise le langage de programmation Java, bien que je pense que cela n'a pas d'importance ici.)

Un grand merci, Thomas.

Answer 1

Je pense avoir trouvé une manière beaucoup plus fiable pour résoudre le « Ajouter 3 » (ou « Ajouter 4 » ou « Ajouter N » problème de nombres.

Tout d'abord, je mis en œuvre mon idée du message original . il a donné lieu à un certain code vestimentaire grand qui semblait, au début, de travailler Toutefois, il a échoué dans le cas suivant:. ajouter Double.MAX_VALUE, 1 et -Double.MAX_VALUE le résultat a été 0

@ commentaires de njuffa m'a inspiré creuser un peu plus profond.. et au http://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/, j'ai trouvé que dans Python, ce problème a été résolu assez bien Pour voir le code complet, j'ai téléchargé la source Python (Python 3 .5.1rc1 - 23/11/2015) à partir https://www.python.org/getit/source/, où l'on peut trouver la méthode suivante (sous PYTHON SOFTWARE FOUNDATION LICENCE 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; 
} 

Cette méthode de sommation est différente de la méthode de Kahan, il utilise une variable nombre de variables de compensation. Lors de l'ajout du numéro i, au plus i des variables de compensation supplémentaires (stockées dans le tableau p) sont utilisées. Cela signifie que si je veux ajouter 3 nombres, j'ai besoin de 3 variables supplémentaires. Pour 4 nombres, j'ai besoin de 4 variables supplémentaires.Étant donné que le nombre de variables utilisées peut augmenter n-n+1 seulement après la n e summand est chargé, je peux traduire le code ci-dessus pour Java comme suit:

/** 
* 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; 
} 

Cette fonction prendra le tableau summands et ajouter la éléments tout en l'utilisant simultanément pour stocker les variables de compensation. Puisque nous chargeons le summand à l'index iavant l'élément de tableau au dit index peut être utilisé pour la compensation, ceci fonctionnera. Puisque le tableau sera petit si le nombre de variables à ajouter est petit et n'échappera pas à la portée de notre méthode, je pense qu'il y a une chance décente qu'il soit alloué directement sur la pile par le JIT, ce qui peut rendre le code assez rapide. Je reconnais que je ne comprenais pas entièrement pourquoi les auteurs du code original traitaient les infinis, les débordements et les NaN comme ils le faisaient. Ici, mon code s'écarte de l'original. (Je l'espère, je ne l'ai pas gâcher.)

De toute façon, je peux maintenant résumer 3, 4 ou ndouble chiffres en faisant:

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 }); 
} 

Si je veux résumer 3 ou 4 long numéros et obtenir le résultat précis comme double, je vais devoir faire face au fait que double s ne peut représenter long s en -9007199254740992..9007199254740992L. Mais cela peut se faire facilement en divisant chaque long en deux parties:

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)), // 
      }); 
} 

Je pense que cela devrait être à peu près juste. Au moins, je peux maintenant ajouter Double.MAX_VALUE, 1, et -Double.MAX_VALUE et obtenir 1 comme résultat.