<cmath>: Restore std::pow(x, 2) accuracy (#5771)

Author: StephanTLavavej

stl/inc/cmath

@@ -567,6 +567,17 @@ double frexp(_Ty _Value, _Out_ int* const _Exp) noexcept /* strengthened */ {
     return _CSTD frexp(static_cast<double>(_Value), _Exp);
 }
 
+template <class _Ty1, class _Ty2, _STD enable_if_t<_STD is_arithmetic_v<_Ty1> && _STD is_arithmetic_v<_Ty2>, int> = 0>
+_NODISCARD _STD _Common_float_type_t<_Ty1, _Ty2> pow(_Ty1 _Left, _Ty2 _Right) noexcept /* strengthened */ {
+    if constexpr (_STD _Is_nonbool_integral<_Ty2>) {
+        if (_Right == 2) {
+            return static_cast<double>(_Left) * static_cast<double>(_Left); // TRANSITION, see GH-5768
+        }
+    }
+
+    return _CSTD pow(static_cast<double>(_Left), static_cast<double>(_Right));
+}
+
 template <class _Ty1, class _Ty2, class _Ty3,
     _STD enable_if_t<_STD is_arithmetic_v<_Ty1> && _STD is_arithmetic_v<_Ty2> && _STD is_arithmetic_v<_Ty3>, int> = 0>
 _NODISCARD _STD _Common_float_type_t<_Ty1, _STD _Common_float_type_t<_Ty2, _Ty3>> fma(
@@ -861,7 +872,7 @@ _GENERIC_MATH1(cbrt)
 _GENERIC_MATH1(fabs)
 _GENERIC_MATH2(hypot)
 // 3-arg hypot() is hand-crafted
-_GENERIC_MATH2(pow)
+// pow() is hand-crafted
 _GENERIC_MATH1(sqrt)
 _GENERIC_MATH1(erf)
 _GENERIC_MATH1(erfc)

stl/inc/type_traits

@@ -1116,9 +1116,6 @@ struct _NO_SPECIALIZATIONS_OF_TYPE_TRAITS is_unsigned : bool_constant<_Sign_base
 _EXPORT_STD template <class _Ty>
 _NO_SPECIALIZATIONS_OF_TYPE_TRAITS constexpr bool is_unsigned_v = _Sign_base<_Ty>::_Unsigned;
 
-template <class _Ty>
-constexpr bool _Is_nonbool_integral = is_integral_v<_Ty> && !is_same_v<remove_cv_t<_Ty>, bool>;
-
 template <bool>
 struct _Select { // Select between aliases that extract either their first or second parameter
     template <class _Ty1, class>

stl/inc/xtr1common

@@ -197,6 +197,9 @@ _NO_SPECIALIZATIONS_OF_TYPE_TRAITS constexpr bool is_integral_v = _Is_any_of_v<r
 _EXPORT_STD template <class _Ty>
 struct _NO_SPECIALIZATIONS_OF_TYPE_TRAITS is_integral : bool_constant<is_integral_v<_Ty>> {};
 
+template <class _Ty>
+constexpr bool _Is_nonbool_integral = is_integral_v<_Ty> && !is_same_v<remove_cv_t<_Ty>, bool>;
+
 _EXPORT_STD template <class _Ty>
 _NO_SPECIALIZATIONS_OF_TYPE_TRAITS constexpr bool is_floating_point_v =
     _Is_any_of_v<remove_cv_t<_Ty>, float, double, long double>;

tests/std/test.lst

@@ -270,6 +270,7 @@ tests\GH_005421_vector_algorithms_integer_class_type_iterator
 tests\GH_005472_do_not_overlap
 tests\GH_005546_containers_size_type_cast
 tests\GH_005553_regex_character_translation
+tests\GH_005768_pow_accuracy
 tests\LWG2381_num_get_floating_point
 tests\LWG2510_tag_classes
 tests\LWG2597_complex_branch_cut

tests/std/tests/GH_005768_pow_accuracy/env.lst

@@ -0,0 +1,4 @@
+# Copyright (c) Microsoft Corporation.
+# SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+
+RUNALL_INCLUDE ..\usual_matrix.lst

tests/std/tests/GH_005768_pow_accuracy/test.cpp

@@ -0,0 +1,125 @@
+// Copyright (c) Microsoft Corporation.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <cstdio>
+#include <functional>
+#include <random>
+#include <vector>
+
+// INTENTIONALLY AVOIDED: using namespace std;
+
+void initialize_randomness(std::mt19937_64& gen) {
+    constexpr std::size_t n = std::mt19937_64::state_size;
+    constexpr std::size_t w = std::mt19937_64::word_size;
+    static_assert(w % 32 == 0, "w should be evenly divisible by 32");
+    constexpr std::size_t k = w / 32;
+
+    std::vector<std::uint32_t> vec(n * k);
+
+    std::random_device rd;
+    std::generate(vec.begin(), vec.end(), std::ref(rd));
+
+    std::printf("This is a randomized test.\n");
+    std::printf("DO NOT IGNORE/RERUN ANY FAILURES.\n");
+    std::printf("You must report them to the STL maintainers.\n\n");
+
+    std::seed_seq seq(vec.cbegin(), vec.cend());
+    gen.seed(seq);
+}
+
+void check_equal(const float val, const float actual, const float squared) {
+    if (actual != squared) {
+        std::printf("val: %.6a; actual: %.6a; squared: %.6a\n", val, actual, squared);
+    }
+    assert(actual == squared);
+}
+
+void check_equal(const double val, const double actual, const double squared) {
+    if (actual != squared) {
+        std::printf("val: %a; actual: %a; squared: %a\n", val, actual, squared);
+    }
+    assert(actual == squared);
+}
+
+void check_equal(const long double val, const long double actual, const long double squared) {
+    if (actual != squared) {
+        std::printf("val: %La; actual: %La; squared: %La\n", val, actual, squared);
+    }
+    assert(actual == squared);
+}
+
+void test_square_flt(const float val, const float squared) {
+    // float * float is float, N5014 [expr.arith.conv]/1.4.1.
+    check_equal(val, val * val, squared);
+
+    // For std::pow(float, int), the arguments are effectively cast to double, N5014 [cmath.syn]/3.
+    const double as_dbl = static_cast<double>(val);
+    check_equal(as_dbl, std::pow(val, 2), as_dbl * as_dbl);
+}
+
+void test_square_dbl(const double val, const double squared) {
+    check_equal(val, val * val, squared);
+    check_equal(val, std::pow(val, 2), squared);
+}
+
+void test_square_ldbl(const long double val, const long double squared) {
+    check_equal(val, val * val, squared);
+    check_equal(val, std::pow(val, 2), squared);
+}
+
+// GH-5768 <cmath>: Calling UCRT ::pow(x, 2) is less accurate than x * x
+void test_gh_5768_manually_verified() {
+    // Manually verified to be correctly rounded:
+    test_square_flt(0x1.bffff4p-1f, 0x1.87ffecp-1f);
+    test_square_flt(0x1.bc0040p-2f, 0x1.810870p-3f);
+    test_square_flt(0x1.2b7bc2p-1f, 0x1.5e5a52p-2f);
+
+    test_square_dbl(0x1.ec9a50154a6f9p-1, 0x1.d9f0c06b2463ep-1);
+    test_square_dbl(0x1.12814d2dd432cp-1, 0x1.26590a84f9b12p-2);
+    test_square_dbl(0x1.33994b0b751ccp-3, 0x1.719905c84494ap-6);
+
+    test_square_ldbl(0x1.ec9a50154a6f9p-1L, 0x1.d9f0c06b2463ep-1L);
+    test_square_ldbl(0x1.12814d2dd432cp-1L, 0x1.26590a84f9b12p-2L);
+    test_square_ldbl(0x1.33994b0b751ccp-3L, 0x1.719905c84494ap-6L);
+}
+
+void test_gh_5768_randomized(std::mt19937_64& gen) {
+    const int Trials = 1'000'000;
+
+    {
+        std::uniform_real_distribution<float> dist{-10.0f, 10.0f};
+        for (int i = 0; i < Trials; ++i) {
+            const auto val = dist(gen);
+            test_square_flt(val, val * val);
+        }
+    }
+
+    {
+        std::uniform_real_distribution<double> dist{-10.0, 10.0};
+        for (int i = 0; i < Trials; ++i) {
+            const auto val = dist(gen);
+            test_square_dbl(val, val * val);
+        }
+    }
+
+    {
+        std::uniform_real_distribution<long double> dist{-10.0L, 10.0L};
+        for (int i = 0; i < Trials; ++i) {
+            const auto val = dist(gen);
+            test_square_ldbl(val, val * val);
+        }
+    }
+}
+
+int main() {
+    test_gh_5768_manually_verified();
+
+    std::mt19937_64 gen;
+    initialize_randomness(gen);
+    test_gh_5768_randomized(gen);
+}