计算总和时使用比需要更高的精度答案

【问题标题】：Use higher precision than needed when computing a sum计算总和时使用比需要更高的精度
【发布时间】：2015-06-03 13:38:29
【问题描述】：

在计算总和时使用较大的精度并在算法结束时降低精度是否是一种好习惯？喜欢

float average(const float* begin, const float* end)
    {
    double sum=0;
    size_t N=end-begin;
    while(begin!=end)
        {
        sum+=(double)(*begin);
        ++begin;
        }

    return (float)( sum/N); //Assume range is not empty
    }

可能是，因为积累的错误少。另一方面，在数据类型之间转换时可能会出错。

【问题讨论】：

在 C 中，所有浮点运算都在 double-precision 无论如何中。 float 仅用于存储。
@EOF 他的观点是中间和的存储精度更高。
@EOF：那是在第一个 C 版本中。在当前版本中，它取决于编译器，但通常涉及浮点数的算术将使用浮点精度完成。
你的方法很有趣，但是使用大精度变量显然会降低程序性能。无论如何，我对你的想法有满意的一面。
@AnishSharma 根据编译器和架构，存储为 double 实际上可能会提高性能。

标签： c floating-accuracy

【解决方案1】：

这取决于您要避免什么，但可能不是。

如果您想避免灾难性的取消（10^100 + 1 - 10^100 的结果是 0 而不是 1），使用更宽的 FP 类型会有所帮助，但作用不大。

如果这些数字在数量级上更接近，但您仍然担心 LSB 随着总和的增长而落到最后（例如1e-8 + 1e-8 + (1e8 copies) != 1），那么更宽的类型可以提供帮助，但同样，只是在一定程度上。

真正有帮助的是更聪明的浮点求和方法。最简单的方法称为“成对求和”，您可以将数字数组视为二叉树的叶子，然后递归地对它们进行求和，直到只剩下一个数字。对于在那里进行的迭代求和，您还可以先对数字进行排序，这往往会减少错误。并且还有更复杂、更精确的方法可用……谷歌“补偿总和”了解详情。

这就是说，如果您怀疑舍入错误会给您带来问题，double sum 会有所帮助，但可能还不够。

哦，关于“在数据类型之间转换时可能会出错”：可能出错（特别是双舍入错误），但您可能会出现不精确与执行求和本身的误差相比，从它们中看出并不重要。

【讨论】：

假设来自某个分布的随机数，量化“不是很多”会很有趣
@user877329 发生特定取消级别的指数差异由尾数大小给出。也就是说，对于float 中间结果，1 + ~6e-8 - 1 处于取消为零的边缘，而对于double 中间结果，1 + ~1e-16 - 1 是。

【解决方案2】：

不是好的一件事是降低最后的精度。

无论如何，您的代码都被零除，因为当您进行除法时，开始 == 结束。

【讨论】：

【解决方案3】：

我不会：做这种事情会进一步将您的实现与特定平台联系起来。在 C 中无法保证 float 的精度低于 double，并且最后降低精度不是好的做法，而且从计算上讲也不是特别便宜。

我会让编译器去做它的工作。

在浮点数中添加数字时，最好先累积较小的数量级数字。这样他们就有更好的机会为这笔款项做出贡献。还有更高级的浮点求和方法；你也应该考虑一下。

【讨论】：

在硬件和软件实现中，float 和 double 之间的转换是您可以执行的最便宜的 FP 操作之一。
马里波恩的房子比梅菲尔便宜，但我还是不会买。顺便说一句，喜欢你的回答 - 加一个。

【解决方案4】：

Sneftel 提到了求和的方法。下面是一组函数，可以处理 2048 个 IEEE 64 位双精度数组（由调用者传递）。（假设 unsigned long long 也是 64 位）。

/* clear array */
void clearsum(double asum[2048])
{
size_t i;
    for(i = 0; i < 2048; i++)
        asum[i] = 0.;
}

/* add a number into array */
void addtosum(double d, double asum[2048])
{
size_t i;
    while(1){
        /* i = exponent of d */
        i = ((size_t)((*(unsigned long long *)&d)>>52))&0x7ff;
        if(i == 0x7ff){         /* max exponent, could be overflow */
            asum[i] += d;
            return;
        }
        if(asum[i] == 0.){      /* if empty slot store d */
            asum[i] = d;
            return;
        }
        d += asum[i];           /* else add slot to d, clear slot */
        asum[i] = 0.;           /* and continue until empty slot */
    }
}

/* return sum from array */
double returnsum(double asum[2048])
{
double sum = 0.;
size_t i;
    for(i = 0; i < 2048; i++)
        sum += asum[i];
    return sum;
}

【讨论】：