<cmath>: Restore std::pow(x, 2) accuracy by StephanTLavavej · Pull Request #5771 · microsoft/STL (original) (raw)

Fixes #5768, a regression that was introduced by #903 in VS 2019 16.8. Reported as DevCom-10968831 / internal VSO-2596129.

WG21-N5014 [cmath.syn]/3 requires std::pow(float, int) to effectively cast its arguments to double. #903 correctly fixed a bug by removing an overload float pow(float, int) that was performing the wrong computation with the wrong return type. See 5d4f4af for the merged commit.

But in the process of removing our handwritten overloads for pow(float, int), pow(double, int), and pow(long double, int), and using our _GENERIC_MATH2 machinery to always call UCRT ::pow with properly double-cast arguments, it revealed an underlying accuracy bug in the UCRT, which we had previously been avoiding.

Our int-power overloads had special cases for squaring, where we directly returned x * x. It turns out that the UCRT mishandles squaring in a small fraction of cases, returning answers that are off by 1 ULP (unit in the last place). This was not only an observable behavioral difference, but a correctness-damaging one for pow(double, int).

(pow(float, int) also had the special case for squaring, but because it was originally returning float instead of double, it was already wildly incorrect. Changing the computation and return type to be double changed the behavior, but towards increased correctness. In fact, so far I have been unable to find any float values that, when widened to double, are improperly squared by UCRT ::pow() compared to directly squaring the double, so I am not yet aware of examples where the special case is actually necessary for float.)

I've reported this to the UCRT maintainers who are investigating, but in the meantime we should fix this regression. My fix preserves the overload improvements of #903 while restoring the special case for squaring. Now that we have if constexpr, this is straightforward (although the analysis of the Standardese, UCRT behavior, and STL history to get here was quite involved). For non-bool integral exponents, we detect when they're 2. In that squaring case, we effectively cast the LHS to double (there's an existing comment in <cmath> explaining why this is always correct), and then directly square that. I have a TRANSITION comment for future maintenance.

(We need to move _Is_nonbool_integral up from <type_traits> to <xtr1common> for this. This avoids "warning C4806: '==': unsafe operation: no value of type '_Ty2' promoted to type 'int' can equal the given constant" emitted by a pathological libcxx test that's passing a bool as the second argument.)

This restores our original behavior for pow(double, int) and pow(long double, int). For pow(float, int) our original behavior was wildly wrong (too narrow), but now the behavior of pow(float, 2) is definitely correct. (As I mentioned, we might not need the special case to avoid misbehavior for float widened to double, but it's easier and safer to use it regardless of the floating-point type).

This PR does not attempt to extend the special case to our pow(float, float) and pow(long double, long double) overloads (that call UCRT powf and powl, respectively). They never had the special case, so there is no behavioral regression. I believe it's the case that pow(long double, 2.0L) will experience the underlying UCRT correctness bug, but attempting to patch it in the STL would lead to surprising behavior. We can't intercept pow(double, double), because that always selects the UCRT. It would be strange if we made pow(long double, long double) behave differently. (In practice it's unlikely to matter due to the low usage of long double.)

The resulting behavior is fairly simple to explain: calling std::pow(x, 2) results in correct answers (double-widened direct squaring). std::pow(x, 2.0) is still subject to the UCRT bug.

I'm adding test coverage with handwritten and randomized cases. This is not especially useful with the current fix (it just confirms that we are indeed directly squaring), but will be useful in the future if the UCRT is ever fixed and we can finally directly call it as #903 attempted to do.

The manually verified cases for float are not as interesting as they might seem, due to the double widening behavior (I wrote them before remembering how pow(float, int) would need to be tested), but I retained them for symmetry, and carefully commented why we need to test pow(float, int) differently.

Mister Pow, that's my name, that name again is Mister Pow.