我有一个功能：

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 and 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 。

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

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

Proof:

• If alpha，我们假设在 [0, 1] 中，大于或等于 1/2，则1-alpha由 Sterbenz 引理所精确。 （“精确”是指计算出的浮点结果等于数学结果；不存在舍入误差。）然后，在计算中A, w0是准确的，因为它是1-alpha, and w1是精确的，因为它的数学结果是alpha，所以它是完全可以表示的。并且，在计算领域B, w0是精确的，因为它的数学结果是alpha, and w1是准确的，因为它又是1-alpha.
• If alpha小于 1/2，那么1-alpha可能有一些舍入误差。设结果为beta。然后，在A, w0 is beta。现在 ½ ≤beta，所以斯特本茨引理适用于评估w1 = 1.0f - w0, so w1是精确的（并且等于数学结果1-beta）。并且，在B, w0是精确的，同样由 Sterbenz 引理得出，并且等于w1 of A, and w1 (of B) 是精确的，因为它的数学结果是beta，这完全可以表示。

现在我们可以看到w0 in A equals w1 in B and w1 in A equals w0 in B。出租beta be 1-alpha在上述任何一种情况下，A and B因此返回(1-beta)*x0 + beta*x1 and beta*x1 + (1-beta)*x0， 分别。 IEEE-754 加法是可交换的（NaN 有效负载除外），因此A and B返回相同的结果。

回答问题：

1. 我想说这是一个合理的做法。如果没有进一步思考，我不会断言没有可以改进的地方。

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

• C++ 允许实现在评估浮点运算时使用超额精度。因此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++ 标准要求在赋值和转换中丢弃多余的精度。这扩展到函数返回。 （我从内存中报告了这个和其他 C++ 规范；应该检查标准。）因此上面的每一个都会将其结果舍入为float即使最初使用了额外的精度。这将防止收缩。

（人们还应该能够通过包括来禁用收缩<cmath>并插入预处理器指令#pragma STDC FP_CONTRACT OFF。有些编译器可能不支持。）

One problem with the workaround above is that values are first rounded to the evaluation precision and then rounded to float. There are mathematical values for which, for such a value x, rounding x first to double (or another precision) and then to float produces a different result than rounding x directly to float. The dissertation A Rigorous Framework for Fully Supporting the IEEE Standard for Floating-Point Arithmetic in High-Level Programming Languages by Samuel A. Figueroa del Cid establishes that evaluating a single operation of multiplication or addition in IEEE-754 basic 64-bit floating-point (commonly used for double) and then rounding to the 32-bit format never has a double-rounding error (because these operations, given inputs that are elements of the 32-bit format, can never produce one of the troublesome x values described above).¹

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

Footnote

¹ Per Figueroa del Cid, if x and y have p-bit significands, and x+y or x*y is computed exactly and then rounded to q places, a second rounding to p places will have the same answer as if the result were directly rounded to p places if p ≤ (q − 1)/2. This is satisfied for IEEE-754 basic 32-bit binary floating-point (p = 24) and 64-bit (q = 53). These formats are commonly used for float and double, and the workaround above should suffice in a C++ implementation that uses them. If a C++ implementation evaluated float using a precision that did not satisfy the condition Figueroa del Cid gives, then double-rounding errors could occur.