3

我有一个功能： 

float lerp(float alpha, float x0, float x1) {
    return (1.0f - alpha) * x0 + alpha * x1;
}

对于那些没有看过的人来说，这比后者更可取，x0 + (x1-x0) * alpha因为后者并不能保证lerp(1.0f, x0, x1) == x1.

现在，我希望我的lerp函数有一个额外的属性：我想要lerp(alpha, x0, x1) == lerp(1-alpha, x1, x0). （至于为什么：这是一个更复杂功能的玩具示例。）我想出的似乎可行的解决方案是 

float lerp_symmetric(float alpha, float x0, float x1) {
    float w0 = 1.0f - alpha;
    float w1 = 1.0f - w0;
    return w0 * x0 + w1 * x1;
}

这种双重减法具有接近零和接近一的舍入效果，因此如果alpha = std::nextafter(0)(1.4012985e-45), then 1 - alpha == 1and so 1 - (1-alpha) == 0。据我所知，这始终是正确的1.0f - x == 1.0f - (1.0f - (1.0f - x))。似乎也有这样的效果w0 + w1 == 1.0f。

问题：

  1. 这是一个合理的方法吗？
  2. 我可以相信我的编译器会做我想做的事吗？特别是，我知道在 Windows 上它有时对部分结果使用更高的精度，而且我知道编译器可以做一些代数；显然 1-(1-x)==x 代数。

这是在 C++11 中使用 Clang、VisualStudio 和 gcc。

4

1 回答 1

1

如果始终使用一种 IEEE-754 二进制浮点格式（例如，基本的 32 位二进制，C++ 常用的格式float），所有 C++ 运算符都以直接和简单的方式映射到 IEEE-754 运算，则lerp_symmetric(alpha, x0, x1)(以下简称A) 等于lerp_symmetric(1-alpha, x1, x0)( B)

证明：

  • 如果alpha我们假设在 [0, 1] 中的 大于或等于 ½，则1-alpha根据 Sterbenz 引理是精确的。（“精确”是指计算的浮点结果等于数学结果；没有舍入误差。）然后，在计算A时，w0因为它是1-alpha，所以w1是精确的，因为它的数学结果是alpha，所以它是精确的有代表性的。并且，在计算中B，w0因为它的数学结果是精确的alpha，并且w1因为它再次是精确的，所以它是精确的1-alpha。
  • 如果alpha小于 ½，则1-alpha可能有一些舍入误差。让结果是beta。那么，在 中A，w0是beta。现在 ½ ≤ beta，所以 Sterbenz 引理适用于 的评估w1 = 1.0f - w0，所以w1是精确的（并且等于 的数学结果1-beta）。并且，在 中B，w0是精确的，同样由 Sterbenz 引理，并且等于w1，A并且w1(of B) 是精确的，因为它的数学结果是beta，它是完全可表示的。

现在我们可以看到w0inA等于w1inB和w1inA等于w0in B。在上述任何一种情况下都设beta为，因此分别返回和。IEEE-754 加法是可交换的（NaN 有效载荷除外），因此返回相同的结果。1-alphaAB(1-beta)*x0 + beta*x1beta*x1 + (1-beta)*x0AB

回答问题：

  1. 我会说这是一个合理的方法。我不会断言没有进一步思考就无法做出改进。

  2. 不，你不能相信你的编译器：

    • C++ 允许实现在评估浮点运算时使用超额精度。因此，即使所有操作数都是 ，w0*x0 + w1*x1也可以使用double,或其他精度进行评估。long doublefloat
    • C++ 允许收缩，除非禁用，因此w0*x0 + w1*x1可以评估为fmaf(w0, x0, w1*x1)，因此对其中一个乘法使用精确算术，但对另一个不使用。

您可以使用以下方法部分解决此问题：

float w0 = 1.0f - alpha;
float w1 = 1.0f - w0;
float t0 = w0*x0;
float t1 = w1*x1;
return t0+t1;

C++ 标准要求在赋值和强制转换中丢弃多余的精度。这扩展到函数返回。（我从内存中报告了这个和其他 C++ 规范；应该检查标准。）因此，float即使最初使用了额外的精度，上述每一项都会将其结果四舍五入。这将防止收缩。

（还应该能够通过包含<cmath>和插入预处理器指令来禁用收缩#pragma STDC FP_CONTRACT off。一些编译器可能不支持。）

上述解决方法的一个问题是，值首先四舍五入到评估精度，然后四舍五入到float. 有一些数学值，对于这样的值x ，首先将 x 舍入到（或double另一个精度）然后产生与将x直接float舍入到 不同的结果。Samuel A. Figueroa del Cid的论文A Rigrous Framework for Fully Supporting the IEEE Standard for Floating-Point Arithmetic in High-Level Programming Languages确定了评估 IEEE-754 基本 64 位浮点中的单次乘法或加法运算（通常用于floatdouble) 然后舍入到 32 位格式永远不会出现双舍入错误（因为这些操作，给定的输入是 32 位格式的元素，永远不会产生上述麻烦的x值之一）。1

如果我对从内存报告的 C++ 规范是正确的，那么只要 C++ 实现使用标称格式或足够宽的格式来评估浮点表达式以满足 Figueroa del Cid 给出的要求，那么上述解决方法应该是完整的.

脚注

1根据 Figueroa del Cid，如果x和y具有p位有效数，并且x+yorx*y被精确计算然后四舍五入到q位，则第二次四舍五入到p位将得到相同的答案，就好像结果直接四舍五入到p位如果p ≤ ( q - 1 )/2。这对于 IEEE-754 基本 32 位二进制浮点 ( p = 24) 和 64 位 ( q = 53) 是满足的。这些格式通常用于floatand double，上面的解决方法在使用它们的 C++ 实现中应该足够了。如果评估了 C++ 实现float使用不满足 Figueroa del Cid 给出的条件的精度，则可能会出现双舍入误差。

于 2018-03-23T00:52:38.723 回答