c++ - 舍入误差减少？

Question

考虑以下函数:

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>

template <typename Type>
inline Type a(const Type dx, const Type a0, const Type z0, const Type b1)
{
    return (std::sqrt(std::abs(2*b1-z0))*dx)+a0;
}

template <typename Type>
inline Type b(const Type dx, const Type a0, const Type z0, const Type a1)
{
    return (std::pow((a1-a0)/dx, 2)+ z0)/2;
}

int main(int argc, char* argv[])
{
    double dx = 1.E-6;
    double a0 = 1;
    double a1 = 2;
    double z0 = -1.E7;
    double b1 = -10;
    std::cout<<std::scientific;
    std::cout<<std::setprecision(std::numeric_limits<double>::digits10);
    std::cout<<a1-a(dx, a0, z0, b(dx, a0, z0, a1))<<std::endl;
    std::cout<<b1-b(dx, a0, z0, a(dx, a0, z0, b1))<<std::endl;
    return 0;
}

在我的机器上，它返回:

0.000000000000000e+00
-1.806765794754028e-07

而不是 (0, 0)。第二个表达式存在较大的舍入误差。

我的问题是:如何在不改变类型的情况下减少每个函数的舍入误差(我需要保留这 2 个函数声明(但公式可以重新排列):它们来自更大的程序)？

Answer 1

遗憾的是，所有浮点类型都因舍入错误而臭名昭著。没有它，它们甚至不能存储 0.1(您可以使用长除法手工证明这一点:二进制等价物是 0b0.0001100110011001100...)。您可能会尝试一些解决方法，例如将该 pow 扩展为硬编码的乘法，但您最终需要对程序进行编码以预测和最小化舍入误差的影响。这里有几个想法:

• 永远不要比较浮点值是否相等。我看到的一些替代比较包括:abs(a-b) < delta，或 percent_difference (a,b) < delta 或什至 abs(a/b-1) < delta，其中 delta 是您确定有效的“适当小”的值用于此特定测试。

• 避免将长数组添加到累加器中；随着累加器变大，数组的末尾可能会因舍入误差而完全丢失。在 Jason Sanders 和 Edward Kandrot 的“Cuda 实例”中，作者建议递归地分别添加每对元素，以便每一步生成的数组大小是前一步的一半，直到获得单元素数组。