对称插值 & 编译器优化

Question

I had a function:

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。

问题：

1. 这种方法合理吗？
2. 我能相信编译器按照我的意愿执行吗？特别是在 Windows 上，它有时会为部分结果使用更高的精度，而且编译器可以进行一些代数运算；显然，1-(1-x)==x。

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

- Ben

---

什么编程语言？C、C++还是其他的？有特定的编译器吗？ - Eric Postpischil

如果在整个过程中使用float精度的IEEE-754算术，我相信您的附加属性就可以得到满足。C和C++允许灵活地评估w0*x0 + w1*x1。您可以通过编写t0 = w0*x0; t1 = w1*x1; return t0+t1;来有些抵消这种情况，因为赋值必须丢弃多余的精度。然后，语言规则仍然允许计算例如t0 = x0*x0的多余精度，然后缩小为float。但是，我相信有一篇论文表明，只要多余精度足够大，这等效于在float中进行计算。 - Eric Postpischil

你有相关的参考资料吗？这些规则在哪里可以找到？ - Ben

编程语言规则，如语言标准中所述。 - Eric Postpischil

你应该添加C++标签。这需要删除其他标签之一，因为限制是五个。优化标签可能是最不适用的。 - Eric Postpischil

Answer 1

如果在整个过程中使用IEEE-754二进制浮点数的一种格式（例如，基本的32位二进制浮点数，这是C++中常用的float格式），并且所有C++运算符以直接简单的方式映射到IEEE-754操作，则 lerp_symmetric(alpha, x0, x1)（以下简称为A）等于 lerp_symmetric(1-alpha, x1, x0)（B）

证明：

如果假设的范围[0, 1]中的alpha大于等于½，则根据Sterbenz引理，1-alpha是精确的（即计算得到的浮点结果等于数学结果，没有舍入误差）。因此，在计算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的w1）是精确的，因为其数学结果是beta，可以完全表示。

现在我们可以看到，在A中，w0等于B中的w1，而在A中，w1等于B中的w0。在上述任何一种情况下，让beta等于1-alpha，因此A和B返回(1-beta)*x0 + beta*x1和beta*x1 + (1-beta)*x0。IEEE-754加法是可交换的（除了NaN有效负载），因此A和B返回相同的结果。

回答问题：

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

2. 不，你不能信任你的编译器。

C++允许在计算浮点数运算时使用多余的精度。因此，即使所有操作数都是float类型，w0*x0 + w1*x1可能会使用double、long double或其他精度进行计算。

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”。这将防止缩减。

（通过包含并插入预处理器指令#pragma STDC FP_CONTRACT OFF，还可以禁用缩减。某些编译器可能不支持该功能。）

上述解决方法的一个问题是，值首先被舍入到评估精度，然后再舍入为float。对于某些数学值，例如某个值x，先将x舍入到double（或其他精度），然后再舍入为float与直接将x舍入到float产生的结果不同。Samuel A. Figueroa del Cid的论文“一个严格支持高级编程语言中IEEE浮点算术标准的框架”指出，在IEEE-754基本64位浮点数（通常用于double）执行单个乘法或加法操作，然后舍入为32位格式永远不会出现双舍入误差（因为这些操作针对32位格式的输入元素，永远不会生成上述麻烦的x值之一）。¹

如果我从记忆中对C++规范的报告正确，那么上述描述的解决方法应该是完整的，只要C++实现在评估浮点表达式时使用名义格式或足够宽的格式来满足Figueroa del Cid给出的要求。

脚注

¹ 根据Figueroa del Cid，如果x和y具有p位有效数字，并且精确计算x+y或x*y，然后四舍五入到q位小数，如果将结果再四舍五入到p位小数，将得到与直接将结果四舍五入到p位小数相同的答案，当且仅当p ≤ (q − 1)/2时。这适用于IEEE-754基本32位二进制浮点数（p = 24）和64位（q = 53）。这些格式通常用于float和double，在使用它们的C++实现中，上述解决方法应该足够。如果C++实现使用的精度不满足Figueroa del Cid给出的条件，则可能会发生双重舍入错误。

- Eric Postpischil

---

谢谢。现在我了解了Sterbenz引理，#pragma STDC FP_CONTRACT off和fma等其他内容。 - Ben