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

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++，数组和函数可以是 class.

的成员

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

/* 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

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

时，总和也会溢出

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

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

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

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

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() 中的 -1 反映了您的循环的形式 for (size_t i = st; i < en; ++ i)

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

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

Answer 4

浮点数有 24 个有效位，而双精度数有 53 个有效位。所以你有 29 个保护位，大约是 5000000 的 100 倍。

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

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

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

你真的应该担心吗？

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

Reduce rounding error of sum of square of floats

c++

numerical-methods

numerical