0

考虑以下函数：

#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)。第二个表达式存在较大的舍入误差。

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

4

2 回答 2

1

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

  • 永远不要比较浮点值是否相等。我见过的一些替代比较包括：abs(ab) < delta，或 percent_difference (a,b) < delta 甚至 abs(a/b-1) < delta，其中 delta 是您确定有效的“适当小”值对于这个特定的测试。

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

于 2013-11-14T01:05:32.207 回答
0

在 a() 中，当您将 a0（恰好为 1）添加到 sqrt()*dx 的小而不精确的结果中时，您会失去精度。

使用提供的值，函数 b() 不会丢失任何精度。

当您在第二个输出中调用 b() 之前的 a() 时，您正在对一个已经不精确的数字进行数学运算，从而加剧了错误。

尝试构造数学运算，以便您首先执行不太可能产生浮点错误的操作，最后执行那些更有可能产生浮点错误的操作。

或者，在您的函数内部，确保它们在“long double”值上运行。例如，以下使用浮点提升在第一次数学运算期间将 double 提升为 long double（注意运算符优先级）：

template <typename Type>
inline Type a(const Type dx, const Type a0, const Type z0, const Type b1)
{
    return (std::sqrt(std::abs(2*static_cast<long double>(b1)-z0))*dx)+a0;
}
于 2013-11-14T00:53:29.147 回答