计算范围内的值时的舍入误差答案

【问题标题】：Round-off error when computing value inside range计算范围内的值时的舍入误差
【发布时间】：2021-12-26 23:00:03
【问题描述】：

对于有限值v0、v1 和r 在[0, 1] 范围内的值v，计算如下，是否始终属于[v0, v1] 范围，或者它可以（稍微）由于四舍五入错误而超出？

double v0;  // Finite
double v1;  // Finite
double r;   // In [0, 1]

double v = v0 * r + v1 * (1.0 - r);

if (v0 <= v1)
    assert(v0 <= v && v <= v1);
else
    assert(v1 <= v && v <= v0);

【问题讨论】：

当v0、v1和r都设置为1./3.（使用IEEE-754 binary64）时，断言失败。
@EricPostpischil 是的，我会对v0 != v1 和 FTZ（将 denorms 刷新为 0）模式感兴趣...
备注：如果还有疑问，可以使用std::clamp保证v在区间内
没有v0 等于v1 且远离次正规的情况是v0 = 7./13、v1 = 6./11 和r = 0x1.5p-50。
@user2052436 即使v0 != v1 和所有涉及的数字都正常，断言也会失败。

标签： c++ c floating-point rounding

【解决方案1】：

是的，可以。这是一个例子：

#include <assert.h>

int main() {

    double v0 = 2.670088631008241e-307;
    double v1 = 2.6700889402193536e-307;
    double  r = 0.9999999999232185;

    double v = v0 * r + v1 * (1.0 - r);

    if (v0 <= v1)
        assert(v0 <= v && v <= v1);
    else
        assert(v1 <= v && v <= v0);

    return 0;
}

这会产生：

Assertion failed: (v0 <= v && v <= v1), function main, file b.cpp, line 12.

在这种情况下计算的v 的值是：

2.67009e-307

【讨论】：

double v = v1 - (v1 - v0) * r; 也有同样的问题。我想知道是否有一个表达式（没有条件）总是会映射到 [v0, v1] ？
@user2052436 好吧，v0 或 v1 将适合（假设两者都可以表示为 doubles）
@user2052436 试试fma(v0, r, fma (v1, -r, v1))。确保使用fma() 的正确工作实现，最好是映射到硬件指令的实现（在现代 x86-64、ARm64、PowerPC 和 GPU 上可用，或例如）。