减少浮点数平方和的舍入误差

Question

我尝试计算浮点数组的平方和。

如何减少舍入误差？

我试图在我的实际程序的内部循环中总结大约 5,000,000 个浮点数。

test.cpp:

#include <iostream>
#include <stdint.h>
template <typename Sum, typename Element>
Sum sum(const size_t st, const size_t en) {
    Sum s = 0;
    for (size_t i = st; i < en; ++ i) {
        s += Element(i)*Element(i); 
    }
    return s;
}
int main() {
    size_t size = 100000;
    std::cout << "double, float: " 
              << sum<double, float>(0,size) << "\n";
    std::cout << "int, int: " 
              << sum<int, int>(0,size) << "\n";
}

输出：

double, float: 3.33328e+14
int, int: 216474736

Answer 1

如果浮点数的格式已知，例如 IEEE，则​​可以使用以浮点数的指数为索引的数组来存储部分和，然后求和以产生总和。在数组更新期间，仅将具有相同指数的浮点数相加并存储到适当位置的数组中。最后的总和从小到大。对于 C++，数组和函数可以是类的成员。

将数组作为参数传递给函数的浮点示例：

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

/* add a number into array */
void addtosum(float f, float asum[256])
{
size_t i;
    while(1){
        /* i = exponent of f */
        i = ((size_t)((*(unsigned int *)&f)>>23))&0xff;
        if(i == 0xff){          /* max exponent, could be overflow */
            asum[i] += f;
            return;
        }
        if(asum[i] == 0.f){     /* if empty slot store f */
            asum[i] = f;
            return;
        }
        f += asum[i];           /* else add slot to f, clear slot */
        asum[i] = 0.f;          /* and continue until empty slot */
    }
}

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

Answer 2

float 有 24 个有效位，而 double 有 53 个有效位。所以你有 29 个保护位，大约是 5000000 的 100 倍。

因此，只有比最大值小 100 倍的值才会出现舍入误差。

另请注意，在 Intel 架构中，浮点寄存器实际上保存 80 位的扩展精度数，其中 63 位有效。

那么只有小于最大 100000 倍的数字才会被截断。

你真的应该担心吗？

Answer 3

如果您只想将连续值的平方相加，请使用公式n*(n+1)*(2n+1)/6计算从1 到n 的所有值的平方和。

只要您使用可以表示结果的类型，就可以消除大多数舍入的影响。例如；

 template<typename Sum> Sum sumsq(size_t n)
 {
     // calculates sum of squares from 1 to x
     //   assumes size_t can be promoted to a Sum

     Sum temp(n);     // force promotion to type Sum
     return temp * (temp + 1)* (2*temp + 1)/6;
 }

 template<typename Sum> Sum alternate_sum(size_t st, size_t en)
 {
        Sum result = sumsq(en - 1);
        if (st > 0) result -= sumsq(st-1);
        return result;
 }

 int main()
 {
     size_t size = 100000;
     std::cout << "double, float: " 
              << alternate_sum<double>(0,size) << "\n";
     std::cout << "int, int: " 
          << alternate_sum<long long>(0,size) << "\n";
 }

请注意，对于 size 等于 100000，使用 int 保存结果会产生未定义的行为（有符号整数类型的溢出）。

alternate_sum() 中的 -1s 反映了您的循环的形式为 for (size_t i = st; i < en; ++ i)

您可以取消使用 size_t 类型作为固定功能，但我将把它留作练习。

顺便说一句：既然你说这段代码在一个内部循环中，值得注意的是，这个公式将比你一直使用的循环快得多。

Answer 4

当您对Element 使用int 类型时，在std::sqrt(std::numeric_limits<int>::max()) 之后的每个i 都会在方块上溢出，在您的系统上可能是46341。 当达到std::numeric_limits<int>::max() 时，总和也会溢出。

您可以使用long 或long long 类型而不是int 来增加这个数字。

在乘法之前将第一个float 存储或转换为double 或long double 也是一个好主意，以减少浮点平方运算的错误。 对一组计算的最后一步进行四舍五入总是比对早期步骤进行四舍五入获得更好的结果，因为您可以避免在内部计算中传播（和增加）表示错误。

如果您真的想要精度，并且不想使用一些复杂的技术重新发明轮子，您可以使用像 GNU Multi-Precision Library 或 Boost Multiprecision 这样的多精度库： https://en.wikipedia.org/wiki/List_of_arbitrary-precision_arithmetic_software

它们比您系统的 long double 类型更精确