此C/C++简化测试用例：
int x = ~~~;
const double s = 1.0 / x;
const double r = 1.0 - x * s;
assert(r >= 0); // fail在数值上不稳定，并且命中Assert，原因是最后一个计算可以用FMA来完成，这会使r变为负值。
Clang默认启用了FMA（从版本14开始），所以它会引起一些有趣的回归。https://godbolt.org/z/avvnnKo5E
有趣的是，如果将最后一个计算一分为二，则不会发出FMA，结果总是非负的：
int x = ~~~;
const double s = 1.0 / x;
const double tmp = x * s;
const double r = 1.0 - tmp;
assert(r >= 0); // OK这是IEEE 754/ FP_CONTRACT的一个保证行为吗？或者这是在玩火，我们应该找到一个数值上更稳定的方法？我找不到任何迹象表明fp压缩只会“局部”发生（在 one expression内），像上面到这样的简单拆分就足以防止它们。
(Of当然，在适当的时候，人们也可以考虑用一个数值上更稳定的算法来代替这个算法，或者在[0.0，1.0]范围内添加一个箝位，但这让人感觉有点笨拙。）
1条答案
按热度按时间t30tvxxf1#
C标准允许以额外的范围和精度计算浮点表达式，因为C2020草案N4849 7.1 [expr.pre] 6规定：
浮点操作数的值和浮点表达式的结果可以用比类型所需的精度和范围更大的精度和范围来表示;因此不改变类型。
然而，注51告诉我们：
The cast and assignment operators must still perform their specific conversions as described in 7.6.1.3, 7.6.3, 7.6.1.8 and 7.6.19.
这意味着赋值或强制转换必须将值转换为名义类型。因此，如果使用了额外的范围或精度，则在执行对
double的赋值时，结果必须转换为实际的double值。（为此，我希望赋值在定义中包括初始化。）因此
1.0 - x * s可以使用融合乘加，但是const double tmp = x * s; const double r = 1.0 - tmp;必须为x * s计算double结果，然后从1.0中减去double结果。请注意，这并不排除
const double tmp = x * s;使用额外的精度来计算x * s，然后再次舍入以获得double结果。在极少数情况下，这可能会产生双舍入错误。其中结果与将x·s的实数算术结果直接舍入为double所得到的结果略有不同。这在实践中不太可能发生; C实现将没有理由以额外的精度计算x * s，然后将其四舍五入为double。还要注意，C和C实现不一定符合C或C++标准。