浮点数的精确总和

Question

我知道a similar question，但我想征求人们对我的算法的意见，以尽可能准确地将浮点数与实际成本相加。

这是我的第一个解决方案：

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.

这个需要 O(n*logn) 而不是正常的 O(n)。真的值得吗？

第二种解决方案来自我正在处理的数据的特征。 这是一个巨大的正数列表，数量级相似。

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];

基本思想是以“二叉树”方式求和。

注意：这是一个伪 C 代码。 step<<=1 表示乘以 2。 这个需要 O(n)。 我觉得可能有更好的方法。你能推荐/批评吗？

Answer 1

Kahan's summation algorithm 比直接求和要精确得多，它的运行时间为 O(n)（比直接求和慢 1-4 倍，具体取决于浮点与数据访问相比的速度。绝对小于 4 倍在桌面硬件上速度较慢，并且没有任何数据移动）。

---

或者，如果您使用通常的 x86 硬件，并且您的编译器允许访问 80 位 long double 类型，则只需使用带有 long double 类型累加器的简单求和算法即可。仅在最后将结果转换为double。

---

如果你真的需要很高的精度，你可以结合以上两种解决方案，在 Kahan 的求和算法中对变量 c、y、t、sum 使用 sum。

Answer 2

我的猜测是你的二元分解几乎和 Kahan 求和一样有效。

这里有一个例子来说明它：

#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 大致相当于那个）
• 第一个近似值与您的二进制分解相同。除了我将总和存储在 data[0] 中并且我关心存储残差。这样求和前后数据的准确总和不变
• 这使我能够通过在第 2 次迭代中对残差求和来近似误差，以纠正第 1 次迭代（相当于对二进制求和应用 Kahan）
• 通过进一步迭代，我可以进一步细化结果，我们会看到收敛

Answer 3

元素将按升序放入堆中，因此您可以使用两个队列。如果数字是预先排序的，这将产生 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();

Answer 4

如果您担心减少求和中的数字误差，那么您可能会对Kahan's algorithm 感兴趣。