<cmath>: Restore std::pow(x, 2) accuracy
#5771
Merged
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
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 todouble. #903 correctly fixed a bug by removing an overloadfloat 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), andpow(long double, int), and using our_GENERIC_MATH2machinery to always call UCRT::powwith properlydouble-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 returnedx * 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 forpow(double, int).(
pow(float, int)also had the special case for squaring, but because it was originally returningfloatinstead ofdouble, it was already wildly incorrect. Changing the computation and return type to bedoublechanged the behavior, but towards increased correctness. In fact, so far I have been unable to find anyfloatvalues that, when widened todouble, are improperly squared by UCRT::pow()compared to directly squaring thedouble, so I am not yet aware of examples where the special case is actually necessary forfloat.)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-boolintegral exponents, we detect when they're2. In that squaring case, we effectively cast the LHS todouble(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_integralup 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 aboolas the second argument.)This restores our original behavior for
pow(double, int)andpow(long double, int). Forpow(float, int)our original behavior was wildly wrong (too narrow), but now the behavior ofpow(float, 2)is definitely correct. (As I mentioned, we might not need the special case to avoid misbehavior forfloatwidened todouble, 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)andpow(long double, long double)overloads (that call UCRTpowfandpowl, respectively). They never had the special case, so there is no behavioral regression. I believe it's the case thatpow(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 interceptpow(double, double), because that always selects the UCRT. It would be strange if we madepow(long double, long double)behave differently. (In practice it's unlikely to matter due to the low usage oflong 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
floatare not as interesting as they might seem, due to thedoublewidening behavior (I wrote them before remembering howpow(float, int)would need to be tested), but I retained them for symmetry, and carefully commented why we need to testpow(float, int)differently.Mister Pow, that's my name, that name again is Mister Pow.