From 9a687d7c474978504ec0513cc7d44bee71ec4f85 Mon Sep 17 00:00:00 2001 From: Sayan Date: Wed, 21 Oct 2026 18:59:38 -0400 Subject: [PATCH 1/3] Use x87 fsqrt for rounded square roots on 32-bit MSVC --- .../boost/numeric/interval/rounded_arith.hpp | 39 +++++++++-- test/CMakeLists.txt | 2 +- test/Jamfile.v2 | 1 + test/sqrt_rounding.cpp | 65 +++++++++++++++++++ 4 files changed, 102 insertions(+), 5 deletions(-) create mode 100644 test/sqrt_rounding.cpp diff --git a/include/boost/numeric/interval/rounded_arith.hpp b/include/boost/numeric/interval/rounded_arith.hpp index 4bb8691..cb15ae1 100644 --- a/include/boost/numeric/interval/rounded_arith.hpp +++ b/include/boost/numeric/interval/rounded_arith.hpp @@ -18,6 +18,37 @@ 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 x87 +// instruction there, as msvc_rounding_control.hpp does for rint. +template inline T sqrt_rounded(const T& x) +{ BOOST_NUMERIC_INTERVAL_using_math(sqrt); return sqrt(x); } + +#if defined(_MSC_VER) && defined(_M_IX86) +inline float sqrt_rounded(const float& x) +{ + float v = x; + __asm fld v + __asm fsqrt + __asm fstp v + return v; +} +inline double sqrt_rounded(const double& x) +{ + double v = x; + __asm fld v + __asm fsqrt + __asm fstp v + return v; +} +inline long double sqrt_rounded(const long double& x) +{ return sqrt_rounded(static_cast(x)); } // same format on MSVC +#endif + +} // namespace detail + /* * Three classes of rounding: exact, std, opp * See documentation for details. @@ -65,9 +96,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 +133,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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a9930f7..7f3630a 100644 --- a/test/CMakeLists.txt +++ b/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() diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index e2f8505..368e207 100644 --- a/test/Jamfile.v2 +++ b/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 : : no ] ] [ run cmp.cpp ] diff --git a/test/sqrt_rounding.cpp b/test/sqrt_rounding.cpp new file mode 100644 index 0000000..30ff6e3 --- /dev/null +++ b/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 +#include + +template +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 +void test() +{ + using namespace boost::numeric::interval_lib; + test_policy >, T>(); + test_policy >, T>(); +} + +int main() +{ + test(); + test(); + // 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(); + return boost::report_errors(); +} From fbfc287371786ba15056a7cd00919be1c4ca0a92 Mon Sep 17 00:00:00 2001 From: Sayan Date: Thu, 22 Oct 2026 18:29:54 -0400 Subject: [PATCH 2/3] Use SSE sqrt on 32-bit MSVC when the compiler targets SSE Select by _M_IX86_FP: _mm_sqrt_ss for float with /arch:SSE or higher, _mm_sqrt_sd for double with /arch:SSE2 or higher, and x87 fsqrt otherwise, so the square root runs on the same unit as the rest of the arithmetic. --- .github/workflows/ci.yml | 7 ++++++ .../boost/numeric/interval/rounded_arith.hpp | 24 +++++++++++++++++-- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 21a2735..9f535e8 100644 --- a/.github/workflows/ci.yml +++ b/.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, diff --git a/include/boost/numeric/interval/rounded_arith.hpp b/include/boost/numeric/interval/rounded_arith.hpp index cb15ae1..0288205 100644 --- a/include/boost/numeric/interval/rounded_arith.hpp +++ b/include/boost/numeric/interval/rounded_arith.hpp @@ -13,6 +13,9 @@ #include #include #include +#if defined(_MSC_VER) && defined(_M_IX86) && _M_IX86_FP >= 1 +# include +#endif namespace boost { namespace numeric { @@ -21,12 +24,23 @@ 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 x87 -// instruction there, as msvc_rounding_control.hpp does for rint. +// Microsoft runtime rounds to nearest whatever the mode, so use the square +// root instruction directly there. template 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; @@ -35,6 +49,11 @@ inline float sqrt_rounded(const float& x) __asm 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; @@ -43,6 +62,7 @@ inline double sqrt_rounded(const double& x) __asm fstp v return v; } +# endif inline long double sqrt_rounded(const long double& x) { return sqrt_rounded(static_cast(x)); } // same format on MSVC #endif From 83f68cd60a07a39073dbfe9ae7f1a929ab7a4078 Mon Sep 17 00:00:00 2001 From: Sayan Date: Fri, 23 Oct 2026 06:03:52 -0400 Subject: [PATCH 3/3] Use a single __asm block for the x87 square root on 32-bit MSVC --- include/boost/numeric/interval/rounded_arith.hpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/include/boost/numeric/interval/rounded_arith.hpp b/include/boost/numeric/interval/rounded_arith.hpp index 0288205..31e3e94 100644 --- a/include/boost/numeric/interval/rounded_arith.hpp +++ b/include/boost/numeric/interval/rounded_arith.hpp @@ -44,9 +44,11 @@ inline float sqrt_rounded(const float& x) inline float sqrt_rounded(const float& x) { float v = x; - __asm fld v - __asm fsqrt - __asm fstp v + __asm { + fld v + fsqrt + fstp v + } return v; } # endif @@ -57,9 +59,11 @@ inline double sqrt_rounded(const double& x) inline double sqrt_rounded(const double& x) { double v = x; - __asm fld v - __asm fsqrt - __asm fstp v + __asm { + fld v + fsqrt + fstp v + } return v; } # endif