浮点数的精确总和(Precise sum of floating point numbers)

2019-07-03 21:34发布

我知道类似的问题 ,但我想请人的意见对我的算法来概括尽可能准确浮点数与实际成本。

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

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,而不是正常的O(N)(LOGN 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的的求和算法是显著比简单求和更精确的,并且它在O(n)的比取决于浮点有多快相比数据访问简单求和慢1-4倍之间运行(某处。在桌面上慢绝对小于4倍硬件,并且不围绕数据的任何改组)。


另外,如果您使用的是常用的x86硬件,如果你的编译器可以访问80位long double类型,只需使用简单的求和算法类型的蓄电池long double 。 只有结果转换为double在最后。


如果你真的需要大量的精密,您可以通过使用结合上述两种解决方案long double变量cytsum在Kahan的的总和算法。



Answer 2:

如果您担心降低你的总和数值的错误,那么你可能有兴趣在Kahan的算法 。



Answer 3:

我的猜测是,你的二进制分解将工作差不多,也是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;
}

我宣布我的操作数易挥发,并与-ffloat店内编译,以避免在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的大致等同于)
  • 第1近似相同的二进制分解。 不同之处在于我存储在[0]数据的总和,并且我关心存储的残基。 通过这种方式,数据的准确和前后总和不变
  • 这使我在为了修正第一次迭代的第二次迭代的残留物相加逼近误差(相当于二进制求和应用Kahan的)
  • 通过进一步迭代我可以进一步细化的结果,我们看到了一个收敛


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


文章来源: Precise sum of floating point numbers