Are floating point contractions only allowed within one expression?

Question

This C/C++ reduced testcase:

int x = ~~~;
const double s = 1.0 / x;
const double r = 1.0 - x * s;

assert(r >= 0); // fail

isn't numerically stable, and hits the assert. The reason is that the last calculation can be done with FMA, which brings r into the negatives.

Clang has enabled FMA by default (since version 14), so it is causing some interesting regressions. Here's a running version: https://godbolt.org/z/avvnnKo5E

---

Interestingly enough, if one splits the last calculation in two, then the FMA is not emitted and the result is always non-negative:

int x = ~~~;
const double s = 1.0 / x;
const double tmp = x * s;
const double r = 1.0 - tmp;

assert(r >= 0); // OK

Is this a guaranteed behavior of IEEE754 / FP_CONTRACT, or is this playing with fire and one should just find a more numerically stable approach? I can't find any indication that fp contractions are only meant to happen "locally" (within one expression), and a simple split like to above is sufficient to prevent them.

(Of course, in due time, one could also think of replacing the algorithm with a more numerically stable one. Or add a clamp into the [0.0, 1.0] range, but that feels hackish.)

Answer 1

The C++ standard allows floating-point expressions to be evaluated with extra range and precision because C++ 2020 draft N4849 7.1 [expr.pre] 6 says:

The values of the floating-point operands and the results of floating-point expressions may be represented in greater precision and range than that required by the type; the types are not changed thereby.

However, note 51 tells us:

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.

The meaning of this is that an assignment or a cast must convert the value to the nominal type. So if extra range or precision has been used, that result must be converted to an actual double value when an assignment to a double is performed. (I expect, for this purpose, assignment includes initialization in a definition.)

So 1.0 - x * s may used a fused multiply-add, but const double tmp = x * s; const double r = 1.0 - tmp; must compute a double result for x * s, then subtract that double result from 1.0.

Note that this does not preclude const double tmp = x * s; from using extra precision to compute x * s and then rounding again to get a double result. In rare instances, this can produce a double-rounding error, where the result differs slightly from what one would get by rounding the real-number-arithmetic result of x•s directly to a double. That is unlikely to occur in practice; a C++ implementation would have no reason to compute x * s with extra precision and then round it to double.

Also note that C and C++ implementations do not necessarily conform to the C or C++ standard.