algorithm - programacion - punto flotante informatica

Suma precisa de números de punto flotante

Question

Soy consciente de una pregunta similar , pero quiero pedir la opinión de la gente sobre mi algoritmo para sumar los números de punto flotante con la mayor precisión posible con costos prácticos.

Aquí está mi primera solución:

put all numbers into a min-absolute-heap. // EDIT as told by comments below
pop the 2 smallest ones.
add them.
put the result back into the heap.
continue until there is only 1 number in the heap.

Este tomaría O (n * logn) en lugar de O (n) normal. ¿Realmente vale la pena?

La segunda solución proviene de la característica de los datos en los que estoy trabajando. Es una enorme lista de números positivos con similar orden de magnitud .

a[size]; // contains numbers, start at index 0
for(step = 1; step < size; step<<=1)
    for(i = step-1; i+step<size; i+=2*step)
        a[i+step] += a[i];
    if(i < size-1)
        a[size-1] += a[i];

La idea básica es hacer la suma en forma de ''árbol binario''.

Nota: es un código pseudo C. step<<=1 significa multiplicar paso por 2. Este tomaría O (n). Siento que podría haber un mejor enfoque. ¿Puedes recomendar / criticar?

Answer 1

Los elementos se colocarán en el montón en orden creciente, por lo que puede usar dos colas en su lugar. Esto produce O (n) si los números están pre-ordenados.

Este pseudocódigo produce los mismos resultados que su algoritmo y se ejecuta en O(n) si la entrada se clasifica previamente y el algoritmo de clasificación detecta que:

Queue<float> leaves = sort(arguments[0]).toQueue();
Queue<float> nodes = new Queue();

popAny = #(){
       if(leaves.length == 0) return nodes.pop();
  else if(nodes.length == 0) return leaves.pop();
  else if(leaves.top() > nodes.top()) return nodes.pop();
  else return leaves.pop();
}

while(leaves.length>0 || nodes.length>1) nodes.push(popAny()+popAny());

return nodes.pop();

Answer 2

Mi conjetura es que su descomposición binaria funcionará casi tan bien como la suma de Kahan.

Aquí hay un ejemplo para ilustrarlo:

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>

void sumpair( float *a, float *b)
{
    volatile float sum = *a + *b;
    volatile float small = sum - std::max(*a,*b);
    volatile float residue = std::min(*a,*b) - small;
    *a = sum;
    *b = residue;
}

void sumpairs( float *a,size_t size, size_t stride)
{
    if (size <= stride*2 ) {
        if( stride<size )
            sumpair(a+i,a+i+stride);
    } else {
        size_t half = 1;
        while(half*2 < size) half*=2;;
        sumpairs( a , half , stride );
        sumpairs( a+half , size-half , stride );
    }
}

void sumpairwise( float *a,size_t size )
{
    for(size_t stride=1;stride<size;stride*=2)
        sumpairs(a,size,stride);
}

int main()
{
    float data[10000000];
    size_t size= sizeof data/sizeof data[0];
    for(size_t i=0;i<size;i++) data[i]=((1<<30)*-1.0+random())/(1.0+random());

    float naive=0;
    for(size_t i=0;i<size;i++) naive+=data[i];
    printf("naive      sum=%.8g/n",naive);

    double dprec=0;
    for(size_t i=0;i<size;i++) dprec+=data[i];
    printf("dble prec  sum=%.8g/n",(float)dprec);

    sumpairwise( data , size );
    printf("1st approx sum=%.8g/n",data[0]);
    sumpairwise( data+1 , size-1);
    sumpairwise( data , 2 );
    printf("2nd approx sum=%.8g/n",data[0]);
    sumpairwise( data+2 , size-2);
    sumpairwise( data+1 , 2 );
    sumpairwise( data , 2 );
    printf("3rd approx sum=%.8g/n",data[0]);
    return 0;
}

Declaré mis operandos volátiles y compilé con -float-store para evitar una mayor precisión en la arquitectura x86

g++  -ffloat-store  -Wl,-stack_size,0x20000000 test_sum.c

y obtener: (0.03125 es 1ULP)

naive      sum=-373226.25
dble prec  sum=-373223.03
1st approx sum=-373223
2nd approx sum=-373223.06
3rd approx sum=-373223.06

Esto merece una pequeña explicación.

• Primero muestro sumas ingenuas
• Luego la suma de precisión doble (Kahan es aproximadamente equivalente a eso)
• La primera aproximación es la misma que tu descomposición binaria. Excepto que almaceno la suma en datos [0] y que me importa almacenar los residuos. De esta manera, la suma exacta de los datos antes y después de la suma no se modifica.
• Esto me permite aproximar el error sumando los residuos en la segunda iteración para corregir la primera iteración (equivalente a aplicar Kahan en la suma binaria)
• Al iterar más puedo refinar aún más el resultado y vemos una convergencia

Answer 3

Si le preocupa reducir el error numérico en su resumen, entonces puede estar interesado en el algoritmo de Kahan .

Answer 4

El algoritmo de suma de Kahan es significativamente más preciso que la suma directa, y se ejecuta en O (n) (en algún lugar entre 1 y 4 veces más lento que la suma directa según la velocidad del punto flotante en comparación con el acceso a los datos. Definitivamente, menos de 4 veces más lento en el escritorio hardware, y sin ningún tipo de mezcla de datos).

Alternativamente, si está utilizando el hardware x86 habitual, y si su compilador permite el acceso al tipo long double 80 bits de long double , simplemente use el algoritmo de resumen directo con el acumulador de tipo long double . Solo convierte el resultado al double al final.

Si realmente necesita mucha precisión, puede combinar las dos soluciones anteriores utilizando el long double para las variables c , y , t , sum en el algoritmo de sum de Kahan.