18

我知道一个类似的问题，但我想征求人们对我的算法的意见，以尽可能准确地将浮点数与实际成本相加。

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

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)。我觉得可能有更好的方法。你能推荐/批评吗？

4

4 回答 4

21

Kahan 的求和算法比直接求和要精确得多，它的运行时间为 O(n)（比直接求和慢 1-4 倍，具体取决于浮点与数据访问相比的速度。在桌面上肯定慢不到 4 倍硬件，并且没有任何数据混洗）。

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

如果你真的需要很高的精度，你可以通过在 Kahan 求和算法中使用long doublefor variables c, y, t,来组合上述两种解决方案。sum

于 2012-11-16T13:40:32.667 回答
9

如果您担心减少求和中的数值误差，那么您可能会对Kahan 算法感兴趣。

于 2012-11-16T13:40:07.237 回答
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 应用于二进制求和）
  • 通过进一步迭代，我可以进一步细化结果，我们看到了收敛
于 2012-11-18T01:16:41.463 回答
1

元素将按递增顺序放入堆中，因此您可以使用两个队列。如果数字是预先排序的，这将产生 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();
于 2012-11-16T14:06:42.823 回答