如何计算在 +、-、*、/ 运算中舍入误差后的浮点精度?

3
为了验证,我希望能够计算出在特定算术计算中舍入到可表示值所导致的累积误差的合理上限。
假设我们有一个名为foo()的函数声称执行某个特定算术计算。同时假设该涉及类型为floatdoublefoo()执行计算的被隐含保证对最大舍入误差(由此产生)的方式(或者说明)。
我想通过以一种追踪最坏情况下误差的方式进行计算来验证foo()的结果是否与特定输入值集的计算结果足够接近,然后检查这两个结果是否与最终最坏情况下误差要求足够接近。
我想可以通过引入一个新的算术类track_prec<T>,向基本浮点数类型之一添加精度跟踪,然后让该类的算术运算符的实现计算每个子表达式的最坏情况下误差。我的问题是我不知道如何在这些一般情况下计算这些最坏情况下的误差。
// T = float or double
template<class T> class track_prec {
public:
    T value;
    T ulp; // http://en.wikipedia.org/wiki/Unit_in_the_last_place

    track_prec& operator+=(const track_prec& v)
    {
        value += v.value;
        ulp = ???; // How to do this? And what about -=, *=, and /=?
    }

    friend bool operator==(T, const track_prec&)
    {
        // Exactly how should this comparison be done?
    }
};

假设,例如,foo() 是对一系列数字的简单求和。那么我们可以按照以下方式使用 track_prec<T>
std::vector<T> values { 0.4, -1.78, 1.3E4, -9.29E3, ... };
CHECK_EQUAL(std::accumulate(values.begin(), values.end(), track_prec<T>()),
            foo(values.begin(), values.end()));

当然,任何形式的帮助都受欢迎,但指向免费且有效代码的指针将非常好。
我在这些链接上找到了相关内容,但它们似乎并没有直接回答我的问题。

1
区间算术很有吸引力,但并不是非常流行。它的硬件支持很少;大多数机器要么不支持必要的舍入模式,要么在切换舍入模式时会产生高昂的执行时间成本。通常,在需要误差界限时,它们是由人类针对特定情况计算的,或者通过各种手段进行估计。 - Eric Postpischil
好的,在我的情况下,性能不是问题,因为只有验证代码才会进行精确跟踪。但是你知道有哪些可访问的C/C++实现可以自动计算误差界限吗? - Kristian Spangsege
2个回答

3

追踪浮点计算精度的最简单方法称为区间算术。它不需要IEEE 754算术,只需要将计算切换到向上或向下舍入,以便在每个步骤计算出的区间边界包含所有实数,这些实数可能是由相同的计算得出的,如果使用实际算术进行计算,则会得到相同的结果。

您应该能够找到许多现有的区间算术实现。

任何给定步骤的计算精度是在该步骤计算的区间宽度。

警惕常量:如果你希望近似于πx,你需要将包含π的浮点区间乘以x的区间。如果你将x的区间乘以双精度数3.1415926535897932,你将得到将x乘以该双精度数的误差(它既不等于π也不等于3.1415926535897932)。在你提出的代码中,常量0.4表示“最接近0.4的双精度数”。如果你需要有理数4/10,使用其边界分别为3.9999999999999997e-014.0000000000000002e-01的区间。

那么使用ULP的想法是个死胡同,对吗? - Kristian Spangsege
@KristianSpangsege:基本上是这样。假设您有两个正数,它们仅在ULP(最低有效位)上有所不同,并且已知它们的正确性在一个ULP之内。如果将它们相加,则结果在一个ULP之内是正确的(新的ULP的价值是旧的两倍)。但是如果你将它们相减呢? - rici
@KristianSpangsege 将错误表示为单个浮点计算结果周围的宽度是可能的,但如果您希望正确执行,则其中的数学并不愉快。假设您有实数n1和n2,表示为(d1 + e1)和(d2 + e2)。 d1和d2相乘的误差是d1 * d2的ULP的一半。但您不想将d1和d2相乘,而是想将n1和n2相乘。因此,您必须编写(d1 + e1)*(d2 + e2)= d1 * d2 + ...,并在新误差中包括由先前近似导致的复合误差。 - Pascal Cuoq
1
@KristianSpangsege,既然您询问了C++实现的问题,下面的文章涉及到了几个实现。其中至少有一个必须允许我推荐的“舍入区间算术”。http://www.ac.usc.es/arith19/sites/default/files/3670a231-spec-session-interval-paper1.pdf - Pascal Cuoq

2

1
这个答案把绝对误差和相对误差混淆了。舍入误差与人们可能用来描述具有绝对不确定度的测量的 x +/- d 表达式是不同的。机器 epsilon 代表相对误差,因为舍入会影响尾数。描述舍入误差的适当表达式是 rounded(x op y) = (x op y)*(1 + d),其中 |d| <= eps。 - Praxeolitic
@Praxeolitic 在我的回答中,ulp是绝对误差的估计。 - 0kcats
@Praxeolitic,您能具体说明原因吗?对于加法,绝对误差被累加,操作的绝对误差为eps * abs(result)。因此,总和考虑了这一点。 - 0kcats
为了保持累积的绝对舍入误差的运行上限,对于“a + b”,我们可以添加先前累积的a的误差、先前累积的b的误差和新引入的舍入误差。对于新的舍入误差,使用结果的ulp是有意义的。对于a和b中先前累积的误差,我们需要一直跟踪它们。没有办法仅从值(例如ulp(a))计算这些数字中的累积误差,因为舍入误差是由操作生成的,它不是浮点值固有的。 - Praxeolitic
这是有道理的。一个基本浮点运算结果的误差不会超过 abs(result) * eps,而且非常接近。 - 0kcats
显示剩余5条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接