Точная сумма чисел с плавающей запятой

Ответ 1

алгоритм суммирования Kahan значительно более точен, чем простое суммирование, и он работает в O (n) (где-то между 1-4 раза медленнее, чем прямое суммирование в зависимости от того, насколько быстро с плавающей запятой сравнивается доступ к данным. Определенно меньше, чем в 4 раза медленнее на настольном оборудовании и без каких-либо перетасовки данных).

В качестве альтернативы, если вы используете обычное аппаратное обеспечение x86, и если ваш компилятор разрешает доступ к 80-разрядному типу long double, просто используйте простой алгоритм суммирования с аккумулятором типа long double. Только преобразовать результат в double в самом конце.

Если вам действительно нужна большая точность, вы можете объединить два вышеупомянутых решения, используя long double для переменных c, y, t, sum в алгоритме суммирования Каха.

Ответ 2

Если вы обеспокоены уменьшением числовой ошибки в суммировании, вы можете быть заинтересованы в алгоритме Кахана.

Ответ 3

Я предполагаю, что ваша двоичная декомпозиция будет работать почти так же, как суммирование Кахана.

Вот пример, чтобы проиллюстрировать это:

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

Я объявил, что мои операнды volatile и скомпилированы с -ffloat-store, чтобы избежать дополнительной точности архитектуры x86

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

и получить: (0.03125 - 1ULP)

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

Это заслуживает небольшого объяснения.

• Сначала я показываю наивное суммирование
• Тогда суммирование с двойной точностью (Kahan примерно эквивалентно такому)
• Первое приближение совпадает с вашим двоичным разложением. Кроме того, что я храню сумму в данных [0] и что я забочусь о хранении остатков. Таким образом, точная сумма данных до и после суммирования не изменяется
• Это позволяет мне аппроксимировать ошибку, суммируя остатки на 2-й итерации, чтобы исправить 1-ю итерацию (что эквивалентно применению Kahan при двоичном суммировании)
• Следуя дальнейшему дальнейшему уточнению результата, мы видим сходимость

Ответ 4

Элементы будут помещаться в кучу в порядке возрастания, поэтому вы можете использовать две очереди. Это дает O (n), если числа предварительно отсортированы.

Этот псевдокод дает те же результаты, что и ваш алгоритм, и работает в O(n), если вход предварительно отсортирован и алгоритм сортировки обнаруживает, что:

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