对称 Lerp 和编译器优化

Question

我有一个功能：

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)，那么 1 - alpha == 1 等等 1 - (1-alpha) == 0。据我所知，1.0f - x == 1.0f - (1.0f - (1.0f - x)) 总是正确的。好像还有w0 + w1 == 1.0f.

的效果

问题：

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

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

Answer 1

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

证明：

如果我们假设在 [0, 1] 中的 alpha 大于或等于 ½，则 1-alpha 根据 Sterbenz 引理是准确的。（“精确”是指计算出的 floating-point 结果等于数学结果；没有舍入误差。）然后，在计算 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 引理，并且等于 A 的 w1，并且 w1（B) 是精确的，因为它的数学结果是 beta，可以精确表示。

现在我们可以看到 A 中的 w0 等于 B 中的 w1 并且 A 中的 w1 等于 w0 在 B。在上述任一情况下，让 beta 为 1-alpha，因此 A 和 B 分别为 return (1-beta)*x0 + beta*x1 和 beta*x1 + (1-beta)*x0。 IEEE-754 加法是可交换的（NaN 有效载荷除外），因此 A 和 B return 相同的结果。

回答问题：

我会说这是一个合理的方法。我不会断言没有进一步思考就可以做出改进。
不，你不能相信你的编译器：
- C++ 允许实现在评估 floating-point 算术时使用超额精度。因此 w0*x0 + w1*x1 可以使用 double、long double 或其他精度计算，即使所有操作数都是 float.
- 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++ 标准要求在赋值和强制转换中放弃过高的精度。这扩展到函数 returns。（我从记忆中报告了这个和其他 C++ 规范；应该检查标准。）因此，即使最初使用了额外的精度，以上每一个都会将其结果四舍五入到 float 。这将防止收缩。

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

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

如果我从记忆中报告的 C++ 规范是正确的，那么只要 C++ 实现使用标称格式或足够宽的格式评估 floating-point 表达式，上述解决方法就应该完成满足 Figueroa del Cid 给出的要求。

脚注

¹ Per Figueroa del Cid，如果 x 和 y 有 p 位有效数，并且x+y 或 x*y 被精确计算，然后四舍五入到 q 位，第二次四舍五入到 p 位将具有如果 p ≤ (q − 1)/2。这满足 IEEE-754 基本 32 位二进制 floating-point (p = 24) 和 64 位 (q = 53 ).这些格式通常用于 float 和 double，上述解决方法在使用它们的 C++ 实现中应该足够了。如果 C++ 实现使用不满足 Figueroa del Cid 给出的条件的精度评估 float，则double-rounding 可能会发生错误。

对称 Lerp 和编译器优化

Symmetrical Lerp & compiler optimizations

c++

floating-point

precision

numerical-stability

lerp

脚注