¿En qué orden deben añadirse flotadores para obtener el resultado más preciso?

También hay un algoritmo diseñado para este tipo de operación de acumulación, llamado Kahan Sumation, que probablemente debería tener en cuenta.

Según Wikipedia,

El algoritmo de suma de Kahan (también conocido como suma compensada) reduce significativamente el error numérico en el total obtenido mediante la adición de una secuencia de números de coma flotante de precisión finita, en comparación con el enfoque obvio. Esto se hace manteniendo un separado compensación de ejecución (una variable para acumular pequeños errores).

En pseudocódigo, el algoritmo es:

function kahanSum(input)
 var sum = input[1]
 var c = 0.0          //A running compensation for lost low-order bits.
 for i = 2 to input.length
  y = input[i] - c    //So far, so good: c is zero.
  t = sum + y         //Alas, sum is big, y small, so low-order digits of y are lost.
  c = (t - sum) - y   //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
  sum = t             //Algebraically, c should always be zero. Beware eagerly optimising compilers!
 next i               //Next time around, the lost low part will be added to y in a fresh attempt.
return sum

Probé el ejemplo extremo en la respuesta proporcionada por Steve Jessop.

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    for (long i = 0; i < billion; ++i)
        sum += small;
    std::cout << std::scientific << std::setprecision(1) << big << " + " << billion << " * " << small << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    sum = 0;
    for (long i = 0; i < billion; ++i)
        sum += small;
    sum += big;
    std::cout  << std::scientific << std::setprecision(1) << billion << " * " << small << " + " << big << " = " <<
        std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Obtuve el siguiente resultado:

1.0e+00 + 1000000000 * 1.0e-09 = 2.000000082740371    (difference = 0.000000082740371)
1000000000 * 1.0e-09 + 1.0e+00 = 1.999999992539933    (difference = 0.000000007460067)

El error en la primera línea es más de diez veces mayor en la segunda.

Si cambio los double s a float s en el código anterior, obtengo:

1.0e+00 + 1000000000 * 1.0e-09 = 1.000000000000000    (difference = 1.000000000000000)
1000000000 * 1.0e-09 + 1.0e+00 = 1.031250000000000    (difference = 0.968750000000000)

Ninguna respuesta está ni siquiera cerca de 2.0 (pero la segunda está un poco más cerca).

Usando la suma de Kahan (con double s) como lo describe Daniel Pryden:

#include <iostream>
#include <iomanip>
#include <cmath>

int main()
{
    long billion = 1000000000;
    double big = 1.0;
    double small = 1e-9;
    double expected = 2.0;

    double sum = big;
    double c = 0.0;
    for (long i = 0; i < billion; ++i) {
        double y = small - c;
        double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }

    std::cout << "Kahan sum  = " << std::fixed << std::setprecision(15) << sum <<
        "    (difference = " << std::fabs(expected - sum) << ")" << std::endl;

    return 0;
}

Consigo exactamente 2.0:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

E incluso si cambio las doubles a float s en el código anterior, obtengo:

Kahan sum  = 2.000000000000000    (difference = 0.000000000000000)

Parece que Kahan es el camino a seguir!

Hay una clase de algoritmos que resuelven este problema exacto, sin la necesidad de ordenar o reordenar los datos.

En otras palabras, la suma se puede hacer en una pasada sobre los datos. Esto también hace que dichos algoritmos sean aplicables en situaciones en las que el conjunto de datos no se conoce de antemano, por ejemplo, si los datos llegan en tiempo real y la suma en ejecución debe mantenerse.

Aquí está el resumen de un artículo reciente:

Presentamos una novela, en línea algoritmo para la suma exacta de una secuencia de números de coma flotante. Por "en línea" queremos decir que el algoritmo necesita ver solo una entrada a la vez, y puede tomar una arbitraria longitud flujo de entrada de tales entradas mientras que requiere solo constante memoria. Por "exacto", queremos decir que la suma de la matriz interna de nuestro algoritmo es exactamente igual a la suma de todas las entradas, y el el resultado devuelto es la suma redondeada correctamente. La prueba de corrección es válido para todas las entradas (incluyendo números no normalizados pero módulo desbordamiento intermedio), y es independiente del número de sumandos o el número de condición de la suma. El algoritmo necesita asintóticamente solo 5 FLOPs por sumando, y debido al paralelismo a nivel de instrucción corre solo alrededor de 2 3 3 veces más lento que el obvio, rápido pero tonto bucle "suma recursiva ordinaria" cuando el número de sumandos es más de 10.000. Por lo tanto, a nuestro conocimiento, es el más rápido, más precisa, y la memoria más eficiente entre algoritmos conocidos. De hecho, es difícil ver cómo un algoritmo más rápido o uno que requiere significativamente menos fracasos podrían existir sin mejoras de hardware. Se proporciona una solicitud para un gran número de sumandos.

Fuente: Algoritmo 908: Suma Exacta en Línea de Flujos de Coma Flotante.

Basándome en la respuesta de Steve de ordenar primero los números en orden ascendente, presentaría dos ideas más:

1. Decidir sobre la diferencia en el exponente de dos números por encima de la cual podría decidir que perdería demasiada precisión.

2. Luego agregue los números en orden hasta que el exponente del acumulador sea demasiado grande para el siguiente número, luego coloque el acumulador en una cola temporal e inicie el acumulador con el siguiente número. Continúa hasta que agota la lista original.

Se repite el proceso con la cola temporal (habiéndola ordenado) y con una diferencia posiblemente mayor en el exponente.

Creo que esto será bastante lento si tienes que calcular exponentes todo el tiempo.

Tuve un rápido ir con un programa y el resultado fue 1.99903

Creo que puedes hacer algo mejor que ordenar los números antes de acumularlos, porque durante el proceso de acumulación, el acumulador se hace más y más grande. Si tiene una gran cantidad de números similares, comenzará a perder precisión rápidamente. Esto es lo que sugeriría en su lugar:

while the list has multiple elements
    remove the two smallest elements from the list
    add them and put the result back in
the single element in the list is the result

Por supuesto, este algoritmo será más eficiente con una cola de prioridad en lugar de una lista. Código C++:

template <typename Queue>
void reduce(Queue& queue)
{
    typedef typename Queue::value_type vt;
    while (queue.size() > 1)
    {
        vt x = queue.top();
        queue.pop();
        vt y = queue.top();
        queue.pop();
        queue.push(x + y);
    }
}

Controlador:

#include <iterator>
#include <queue>

template <typename Iterator>
typename std::iterator_traits<Iterator>::value_type
reduce(Iterator begin, Iterator end)
{
    typedef typename std::iterator_traits<Iterator>::value_type vt;
    std::priority_queue<vt> positive_queue;
    positive_queue.push(0);
    std::priority_queue<vt> negative_queue;
    negative_queue.push(0);
    for (; begin != end; ++begin)
    {
        vt x = *begin;
        if (x < 0)
        {
            negative_queue.push(x);
        }
        else
        {
            positive_queue.push(-x);
        }
    }
    reduce(positive_queue);
    reduce(negative_queue);
    return negative_queue.top() - positive_queue.top();
}

Los números en la cola son negativos porque topproduce el número más grande, pero queremos el más pequeño. Podría haber proporcionado más argumentos de plantilla a la cola, pero este enfoque parece más simple.

Esto no responde completamente a su pregunta, pero una cosa inteligente es ejecutar la suma dos veces, una vez con modo de redondeo "redondear hacia arriba" y una vez con "redondear hacia abajo". Compare las dos respuestas, y usted sabe / cómo / inexacto sus resultados son, y si por lo tanto necesita utilizar una estrategia de suma más inteligente. Desafortunadamente, la mayoría de los idiomas no hacen que cambiar el modo de redondeo de punto flotante sea tan fácil como debería ser, porque la gente no sabe que realmente es útil en todos los días calculo.

Echa un vistazo a Aritmética de intervalos donde haces todas las matemáticas como esta, manteniendo los valores más altos y más bajos a medida que avanzas. Conduce a algunos resultados interesantes y optimizaciones.

El más simple ordenar que mejora la precisión es ordenar por el valor absoluto ascendente. Eso permite que los valores de magnitud más pequeños tengan la oportunidad de acumularse o cancelarse antes de interactuar con valores de magnitud más grandes que desencadenarían una pérdida de precisión.

Dicho esto, puede hacerlo mejor rastreando múltiples sumas parciales no superpuestas. Aquí hay un documento que describe la técnica y presenta una prueba de precisión: www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps

Ese algoritmo y otros enfoques para la suma exacta de coma flotante se implementan en Python simple en: http://code.activestate.com/recipes/393090 / Al menos dos de ellos se pueden convertir trivialmente a C++.

Para IEEE 754 números de precisión simple o doble o de formato conocido, otra alternativa es usar una matriz de números (pasados por el llamante, o en una clase para C++) indexados por el exponente. Cuando se agregan números a la matriz, solo se agregan números con el mismo exponente (hasta que se encuentra un espacio vacío y se almacena el número). Cuando se pide una suma, la matriz se suma de menor a mayor para minimizar el truncamiento. Ejemplo de precisión simple:

/* clear array */
void clearsum(float asum[256])
{
size_t i;
    for(i = 0; i < 256; i++)
        asum[i] = 0.f;
}

/* add a number into array */
void addtosum(float f, float asum[256])
{
size_t i;
    while(1){
        /* i = exponent of f */
        i = ((size_t)((*(unsigned int *)&f)>>23))&0xff;
        if(i == 0xff){          /* max exponent, could be overflow */
            asum[i] += f;
            return;
        }
        if(asum[i] == 0.f){     /* if empty slot store f */
            asum[i] = f;
            return;
        }
        f += asum[i];           /* else add slot to f, clear slot */
        asum[i] = 0.f;          /* and continue until empty slot */
    }
}

/* return sum from array */
float returnsum(float asum[256])
{
float sum = 0.f;
size_t i;
    for(i = 0; i < 256; i++)
        sum += asum[i];
    return sum;
}

Doble precisión ejemplo:

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7ff){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

Sus flotadores deben añadirse con doble precisión. Eso le dará más precisión adicional que cualquier otra técnica puede. Para un poco más de precisión y significativamente más velocidad, puede crear digamos cuatro sumas, y sumarlas al final.

Si está agregando números de doble precisión, use long double para la suma; sin embargo, esto solo tendrá un efecto positivo en implementaciones donde long double realmente tiene más precisión que double (típicamente x86, PowerPC dependiendo de configuración del compilador).

Con respecto a la ordenación, me parece que si espera una cancelación, los números deben agregarse en orden descendente de magnitud, no ascendente. Por ejemplo:

((-1 + 1) + 1e-20) dará 1e-20

Pero

((1e-20 + 1) - 1) dará 0

En la primera ecuación que dos números grandes se cancelan, mientras que en la segunda el término 1e-20 se pierde cuando se agrega a 1, ya que no hay suficiente precisión para retenerlo.

También, la suma en pares es bastante decente para sumar muchos números.