我知道类似的问题 ,但我想请人的意见对我的算法来概括尽可能准确浮点数与实际成本。
这是我的第一个解决方案:
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)。 我觉得有可能是一个更好的办法。 你能推荐/批?
Kahan的的求和算法是显著比简单求和更精确的,并且它在O(n)的比取决于浮点有多快相比数据访问简单求和慢1-4倍之间运行(某处。在桌面上慢绝对小于4倍硬件,并且不围绕数据的任何改组)。
另外,如果您使用的是常用的x86硬件,如果你的编译器可以访问80位long double
类型,只需使用简单的求和算法类型的蓄电池long double
。 只有结果转换为double
在最后。
如果你真的需要大量的精密,您可以通过使用结合上述两种解决方案long double
变量c
, y
, t
, sum
在Kahan的的总和算法。
如果您担心降低你的总和数值的错误,那么你可能有兴趣在Kahan的算法 。
我的猜测是,你的二进制分解将工作差不多,也是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的)
- 通过进一步迭代我可以进一步细化的结果,我们看到了一个收敛
该元素将被放入堆递增的顺序,这样你就可以使用两个队列来代替。 这产生为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();