Use x87 fsqrt for rounded square roots on 32-bit MSVC

.github/workflows/ci.yml

@@ -434,7 +434,14 @@ jobs:
               toolset: msvc-14.3, cxxstd: "latest",          address-model: "64",    os: windows-2022 }
           - { name: "Collect coverage w/ VS 2026", coverage: "yes",
               toolset: msvc-14.5, cxxstd: "latest",          address-model: "64",    os: windows-2025-vs2026 }
+          # 32-bit code for the x87 unit only, and for SSE2 (the default)
+          - { name: "msvc-14.3, /arch:IA32",
+              toolset: msvc-14.3, cxxstd: "14,latest",       address-model: "32",    os: windows-2022, b2_flags: "cxxflags=/arch:IA32" }
+          - { name: "msvc-14.3, /arch:SSE2",
+              toolset: msvc-14.3, cxxstd: "14,latest",       address-model: "32",    os: windows-2022, b2_flags: "cxxflags=/arch:SSE2" }
           - { toolset: clang-win, cxxstd: "14,17,latest",    address-model: "32,64", os: windows-2025, b2_flags: "embed-manifest-via=linker" }
+          - { name: "clang-win, /arch:IA32",
+              toolset: clang-win, cxxstd: "14,latest",       address-model: "32",    os: windows-2025, b2_flags: "embed-manifest-via=linker cxxflags=/arch:IA32" }
           - { toolset: gcc,       cxxstd: "11,14,17,2a",     address-model: "64",    os: windows-2022 }
           # 32-bit x87 code from a current GCC (#41)
           - { name: "gcc, MinGW32", msystem: MINGW32, msys2-package: mingw-w64-i686-gcc,

include/boost/numeric/interval/rounded_arith.hpp

@@ -13,11 +13,66 @@
 #include <boost/numeric/interval/rounding.hpp>
 #include <boost/numeric/interval/detail/bugs.hpp>
 #include <boost/config/no_tr1/cmath.hpp>
+#if defined(_MSC_VER) && defined(_M_IX86) && _M_IX86_FP >= 1
+#  include <emmintrin.h>
+#endif
 
 namespace boost {
 namespace numeric {
 namespace interval_lib {
 
+namespace detail {
+
+// Square root in the current rounding mode.  The sqrt of the 32-bit
+// Microsoft runtime rounds to nearest whatever the mode, so use the square
+// root instruction directly there.
+template<class T> inline T sqrt_rounded(const T& x)
+{ BOOST_NUMERIC_INTERVAL_using_math(sqrt); return sqrt(x); }
+
+#if defined(_MSC_VER) && defined(_M_IX86)
+// Use the unit the compiler uses for the rest of the arithmetic on the type,
+// which /arch selects and _M_IX86_FP reports:
+//   /arch:IA32  (0): float and double on x87
+//   /arch:SSE   (1): float on SSE, double on x87
+//   /arch:SSE2  (2): float and double on SSE2 (the default)
+// This avoids moving the value between units, and relies only on the
+// rounding mode of the unit the type already uses.
+# if _M_IX86_FP >= 1
+inline float sqrt_rounded(const float& x)
+{ return _mm_cvtss_f32(_mm_sqrt_ss(_mm_set_ss(x))); }
+# else
+inline float sqrt_rounded(const float& x)
+{
+  float v = x;
+  __asm {
+    fld v
+    fsqrt
+    fstp v
+  }
+  return v;
+}
+# endif
+# if _M_IX86_FP >= 2
+inline double sqrt_rounded(const double& x)
+{ return _mm_cvtsd_f64(_mm_sqrt_sd(_mm_setzero_pd(), _mm_set_sd(x))); }
+# else
+inline double sqrt_rounded(const double& x)
+{
+  double v = x;
+  __asm {
+    fld v
+    fsqrt
+    fstp v
+  }
+  return v;
+}
+# endif
+inline long double sqrt_rounded(const long double& x)
+{ return sqrt_rounded(static_cast<double>(x)); }  // same format on MSVC
+#endif
+
+} // namespace detail
+
 /*
  * Three classes of rounding: exact, std, opp
  * See documentation for details.
@@ -65,9 +120,9 @@ struct rounded_arith_std: Rounding {
   T div_up  (const T& x, const T& y) { BOOST_UP(x / y); }
   T median(const T& x, const T& y) { BOOST_NR((x + y) / 2); }
   T sqrt_down(const T& x)
-  { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_DN(sqrt(x)); }
+  { BOOST_DN(detail::sqrt_rounded(x)); }
   T sqrt_up  (const T& x)
-  { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_UP(sqrt(x)); }
+  { BOOST_UP(detail::sqrt_rounded(x)); }
   T int_down(const T& x) { this->downward(); return this->to_int(x); }
   T int_up  (const T& x) { this->upward();   return this->to_int(x); }
 # undef BOOST_DN
@@ -102,9 +157,9 @@ struct rounded_arith_opp: Rounding {
   T div_up  (const T& x, const T& y) { BOOST_UP(x / y); }
   T median  (const T& x, const T& y) { BOOST_NR((x + y) / 2); }
   T sqrt_down(const T& x)
-  { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_DN(sqrt(x)); }
+  { BOOST_DN(detail::sqrt_rounded(x)); }
   T sqrt_up  (const T& x)
-  { BOOST_NUMERIC_INTERVAL_using_math(sqrt); BOOST_UP(sqrt(x)); }
+  { BOOST_UP(detail::sqrt_rounded(x)); }
   T int_down(const T& x) { return -this->to_int(-x); }
   T int_up  (const T& x) { return  this->to_int(x); }
 # undef BOOST_DN

test/CMakeLists.txt

@@ -47,7 +47,7 @@ if(BOOST_NUMERIC_INTERVAL_TEST_SELF_CONTAINED_HEADERS)
   endforeach()
 endif()
 
-foreach(t add det fmod mul nth_root overflow pi pow cmp cmp_exp cmp_lex cmp_set test_float)
+foreach(t add det fmod mul nth_root overflow pi pow sqrt_rounding cmp cmp_exp cmp_lex cmp_set test_float)
   boost_test(TYPE run SOURCES ${t}.cpp)
 endforeach()
 

test/Jamfile.v2

@@ -74,6 +74,7 @@ if [ os.environ BOOST_NUMERIC_INTERVAL_TEST_SELF_CONTAINED_HEADERS ]
     [ run overflow.cpp ]
     [ run pi.cpp       ]
     [ run pow.cpp      ]
+    [ run sqrt_rounding.cpp ]
     [ run rounding_mode.cpp : : : [ check-target-builds /boost/architecture//power : : <build>no ] ]
 
     [ run cmp.cpp ]

test/sqrt_rounding.cpp

@@ -0,0 +1,65 @@
+/* Boost test/sqrt_rounding.cpp
+ * the rounded square roots of the rounding policies follow the rounding mode
+ *
+ * Distributed under the Boost Software License, Version 1.0.
+ * (See accompanying file LICENSE_1_0.txt or
+ * copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+// The sqrt of the 32-bit Microsoft runtime rounds to nearest whatever the
+// rounding mode, so sqrt_down and sqrt_up returned the same value there.
+//
+// sqrt is correctly rounded, so for a non-square x the two ends must differ,
+// and for an exact square both must be the exact root.  Each call reads its
+// operand from volatile storage, so that GCC cannot merge the two roots of
+// the same operand (a separate bug, GCC bug 34678).
+
+#include <boost/numeric/interval.hpp>
+#include <boost/core/lightweight_test.hpp>
+
+template<class Rounding, class T>
+void test_policy()
+{
+  static volatile T non_squares[] = { 2, 3, 5, 7, 10, 1.5f, 0.75f };
+  static volatile T squares[] = { 1, 4, 9, 2.25f, 0.25f };
+  const int n_non = sizeof(non_squares) / sizeof(non_squares[0]);
+  const int n_sq = sizeof(squares) / sizeof(squares[0]);
+  for (int i = 0; i < n_non; ++i) {
+    Rounding rnd;
+    T a = non_squares[i];
+    T dn = rnd.sqrt_down(a);
+    T b = non_squares[i];
+    T up = rnd.sqrt_up(b);
+    BOOST_TEST(dn < up);
+  }
+  for (int i = 0; i < n_sq; ++i) {
+    Rounding rnd;
+    T a = squares[i];
+    T dn = rnd.sqrt_down(a);
+    T b = squares[i];
+    T up = rnd.sqrt_up(b);
+    T x = squares[i];
+    BOOST_TEST(dn == up);
+    BOOST_TEST(dn * dn == x);
+  }
+}
+
+template<class T>
+void test()
+{
+  using namespace boost::numeric::interval_lib;
+  test_policy<save_state<rounded_arith_opp<T> >, T>();
+  test_policy<save_state<rounded_arith_std<T> >, T>();
+}
+
+int main()
+{
+  test<float>();
+  test<double>();
+  // Only where long double has the format of double, as with MSVC.  An x87
+  // long double result is not yet ordered before the next mode change, which
+  // is a separate bug.
+  if (sizeof(long double) == sizeof(double))
+    test<long double>();
+  return boost::report_errors();
+}