Keep the rounded primitives from being merged across rounding-mode changes

.github/workflows/ci.yml

@@ -212,7 +212,7 @@ jobs:
       - name: Setup multiarch
         if: matrix.multiarch
         run: |
-            sudo apt-get -o Acquire::Retries=$NET_RETRY_COUNT -y -q --no-install-suggests --no-install-recommends install binfmt-support qemu-user-static
+            sudo apt-get -o Acquire::Retries=$NET_RETRY_COUNT -y -q --no-install-suggests --no-install-recommends install binfmt-support
             git clone --depth=1 https://github.com/jeking3/bdde.git ../bdde
             sudo ../bdde/bin/linux/bdde-multiarch
             sudo ../bdde/bin/linux/bdde-multiarch --version

include/boost/numeric/interval/detail/x86_rounding_control.hpp

@@ -85,11 +85,20 @@ struct x86_rounding_control_long_double<false>: x86_rounding_control
   { volatile long double r_ = r; return r_; }
 };
 
+// An x87 long double needs no rounding to memory, but the result must still be
+// computed before the next change of rounding mode.  The empty asm takes the
+// value as an operand, which makes the computation a dependency of it, and the
+// "memory" clobber keeps it before the mode change.
 template<>
 struct x86_rounding_control_long_double<true>: x86_rounding_control
 {
+#if defined(__GNUC__)
+  static long double force_rounding(long double r)
+  { __asm__ __volatile__ ("" : "+t"(r) : : "memory"); return r; }
+#else
   static long double const &force_rounding(long double const &r)
   { return r; }
+#endif
 };
 
 } // namespace detail

include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp

@@ -40,7 +40,7 @@ struct x86_rounding
   static T to_int(T r)
   {
     T r_;
-    __asm__ ("frndint" : "=&t"(r_) : "0"(r));
+    __asm__ __volatile__ ("frndint" : "=&t"(r_) : "0"(r));
     return r_;
   }
 };

include/boost/numeric/interval/rounded_arith.hpp

@@ -23,11 +23,12 @@ namespace interval_lib {
 
 namespace detail {
 
-// Read a floating operand back through a volatile object. Two square roots
-// of the same operand taken under different rounding modes must not be
-// merged by the optimiser (GCC bug 34678); reading the operand afresh after
-// the mode change makes each evaluation depend on a value the compiler
-// cannot reuse. Used only inside the sqrt bodies below.
+// Read a floating operand back through a volatile object. The compiler does
+// not see a rounding-mode change as a dependency, so it may compute the same
+// expression once for both ends of an interval (GCC bug 34678); reading one
+// operand afresh after the mode change makes each evaluation depend on a
+// value the compiler cannot reuse. Used by the rounded_arith_std primitives
+// and the sqrt bodies below.
 template<class T> inline T reload_operand(const T& x) { return x; }
 inline float reload_operand(const float& x)
 { volatile float v = x; return v; }
@@ -123,23 +124,23 @@ struct rounded_arith_std: Rounding {
 # define BOOST_NR(EXPR) this->to_nearest(); return this->force_rounding(EXPR)
 # define BOOST_UP(EXPR) this->upward();     return this->force_rounding(EXPR)
   void init() { }
-  template<class U> T conv_down(U const &v) { BOOST_DN(v); }
-  template<class U> T conv_up  (U const &v) { BOOST_UP(v); }
-  T add_down(const T& x, const T& y) { BOOST_DN(x + y); }
-  T sub_down(const T& x, const T& y) { BOOST_DN(x - y); }
-  T mul_down(const T& x, const T& y) { BOOST_DN(x * y); }
-  T div_down(const T& x, const T& y) { BOOST_DN(x / y); }
-  T add_up  (const T& x, const T& y) { BOOST_UP(x + y); }
-  T sub_up  (const T& x, const T& y) { BOOST_UP(x - y); }
-  T mul_up  (const T& x, const T& y) { BOOST_UP(x * y); }
-  T div_up  (const T& x, const T& y) { BOOST_UP(x / y); }
+  template<class U> T conv_down(U const &v) { BOOST_DN(detail::reload_operand(v)); }
+  template<class U> T conv_up  (U const &v) { BOOST_UP(detail::reload_operand(v)); }
+  T add_down(const T& x, const T& y) { BOOST_DN(detail::reload_operand(x) + y); }
+  T sub_down(const T& x, const T& y) { BOOST_DN(detail::reload_operand(x) - y); }
+  T mul_down(const T& x, const T& y) { BOOST_DN(detail::reload_operand(x) * y); }
+  T div_down(const T& x, const T& y) { BOOST_DN(detail::reload_operand(x) / y); }
+  T add_up  (const T& x, const T& y) { BOOST_UP(detail::reload_operand(x) + y); }
+  T sub_up  (const T& x, const T& y) { BOOST_UP(detail::reload_operand(x) - y); }
+  T mul_up  (const T& x, const T& y) { BOOST_UP(detail::reload_operand(x) * y); }
+  T div_up  (const T& x, const T& y) { BOOST_UP(detail::reload_operand(x) / y); }
   T median(const T& x, const T& y) { BOOST_NR((x + y) / 2); }
   T sqrt_down(const T& x)
   { BOOST_DN(detail::sqrt_rounded(detail::reload_operand(x))); }
   T sqrt_up  (const T& x)
   { BOOST_UP(detail::sqrt_rounded(detail::reload_operand(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); }
+  T int_down(const T& x) { BOOST_DN(this->to_int(detail::reload_operand(x))); }
+  T int_up  (const T& x) { BOOST_UP(this->to_int(detail::reload_operand(x))); }
 # undef BOOST_DN
 # undef BOOST_NR
 # undef BOOST_UP
@@ -175,8 +176,8 @@ struct rounded_arith_opp: Rounding {
   { BOOST_DN(detail::sqrt_rounded(detail::reload_operand(x))); }
   T sqrt_up  (const T& x)
   { BOOST_UP(detail::sqrt_rounded(detail::reload_operand(x))); }
-  T int_down(const T& x) { return -this->to_int(-x); }
-  T int_up  (const T& x) { return  this->to_int(x); }
+  T int_down(const T& x) { BOOST_UP_NEG(this->to_int(-x)); }
+  T int_up  (const T& x) { BOOST_UP(this->to_int(x)); }
 # undef BOOST_DN
 # undef BOOST_NR
 # undef BOOST_UP

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 sqrt_rounding sqrt_point cmp cmp_exp cmp_lex cmp_set test_float)
+foreach(t add det fmod mul nth_root overflow pi pow rounding_merge sqrt_rounding sqrt_point 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 rounding_merge.cpp ]
     [ run sqrt_rounding.cpp ]
     [ run sqrt_point.cpp ]
     [ run rounding_mode.cpp : : : [ check-target-builds /boost/architecture//power : : <build>no ] ]

test/rounding_merge.cpp

@@ -0,0 +1,157 @@
+/* Boost test/rounding_merge.cpp
+ * the rounded primitives of the rounding policies follow the rounding mode
+ *
+ * Copyright 2026 Sayan Samanta
+ *
+ * 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 interval functions call the _down and the _up primitive of a policy on
+// the same operands, one after the other, with a rounding-mode change in
+// between.  GCC does not see the mode change as a dependency, so it may
+// evaluate x + y once and use the result for both ends (GCC bug 34678).  The
+// interval is then too narrow and may miss the true value.
+//
+// Each check below calls both primitives of the same policy object on the
+// same plain local operands, so that the compiler is free to merge them.
+// Only the input arrays are volatile, so that nothing is folded at build
+// time.  For an inexact operation the two results must differ; for an exact
+// one both must be the exact value.  The square root is checked in
+// sqrt_rounding.cpp and sqrt_point.cpp.
+//
+// Built by the test suite with -frounding-math (see Jamfile.v2 /
+// CMakeLists.txt).  Control: add -DCONTROL_EXACT_POLICY to run the same
+// checks on the exact (non-rounding) policy.  The control must FAIL; if it
+// passes, the test cannot see a merge.
+
+#include <boost/numeric/interval.hpp>
+#include <boost/core/lightweight_test.hpp>
+#include <boost/limits.hpp>
+#include <cmath>
+
+using namespace boost::numeric;
+using namespace boost::numeric::interval_lib;
+
+#define BINARY_TEST(OP, EXACT)                                         \
+template<class Rounding, class T>                                      \
+void test_##OP()                                                       \
+{                                                                      \
+  /* x = a / b and y = c / d, divided at run time: g++ 12.3 with      \
+     -frounding-math miscompiles a static initializer that has to     \
+     round, such as 0.1L for a float.  The quotients are stored to    \
+     volatile T, so that x87 excess precision is dropped. */          \
+  static volatile T inexact[][4] =                                     \
+    { { 3, 1, 1, 10 }, { 1, 10, 3, 1 }, { 7, 1, 7, 10 }, { 7, 10, 7, 1 } }; \
+  const int n = sizeof(inexact) / sizeof(inexact[0]);                  \
+  for (int i = 0; i < n; ++i) {                                        \
+    volatile T xv = T(inexact[i][0]) / T(inexact[i][1]);               \
+    volatile T yv = T(inexact[i][2]) / T(inexact[i][3]);               \
+    T x = xv, y = yv;                                                  \
+    T dn, up;                                                          \
+    {                                                                  \
+      Rounding rnd;                                                    \
+      dn = rnd.OP##_down(x, y);                                        \
+      up = rnd.OP##_up(x, y);                                          \
+    }                                                                  \
+    BOOST_TEST(dn < up);                                               \
+  }                                                                    \
+  static volatile T exact[2] = { 1, 2 };                               \
+  T x = exact[0], y = exact[1];                                        \
+  T dn, up;                                                            \
+  {                                                                    \
+    Rounding rnd;                                                      \
+    dn = rnd.OP##_down(x, y);                                          \
+    up = rnd.OP##_up(x, y);                                            \
+  }                                                                    \
+  BOOST_TEST(dn == up);                                                \
+  BOOST_TEST(dn == T(EXACT));                                          \
+}
+
+BINARY_TEST(add, 3)
+BINARY_TEST(sub, -1)
+BINARY_TEST(mul, 2)
+BINARY_TEST(div, 0.5)
+
+#undef BINARY_TEST
+
+template<class Rounding, class T>
+void test_int()
+{
+  using std::floor;
+  using std::ceil;
+  static volatile T values[] = { 2.5f, -2.5f };
+  const int n = sizeof(values) / sizeof(values[0]);
+  for (int i = 0; i < n; ++i) {
+    T x = values[i];
+    T dn, up;
+    {
+      Rounding rnd;
+      dn = rnd.int_down(x);
+      up = rnd.int_up(x);
+    }
+    BOOST_TEST(dn == floor(x));
+    BOOST_TEST(up == ceil(x));
+  }
+}
+
+// Only a conversion to float from double is inexact here.
+template<class Rounding>
+void test_conv()
+{
+  static volatile double values[] = { 0.1, -0.1 };
+  const int n = sizeof(values) / sizeof(values[0]);
+  for (int i = 0; i < n; ++i) {
+    double v = values[i];
+    float dn, up;
+    {
+      Rounding rnd;
+      dn = rnd.conv_down(v);
+      up = rnd.conv_up(v);
+    }
+    BOOST_TEST(dn < up);
+  }
+}
+
+template<class Rounding, class T>
+void test_policy()
+{
+  test_add<Rounding, T>();
+  test_sub<Rounding, T>();
+  test_mul<Rounding, T>();
+  test_div<Rounding, T>();
+  test_int<Rounding, T>();
+}
+
+#if defined(CONTROL_EXACT_POLICY)
+template<class T> struct std_policy
+{ typedef save_state<rounded_arith_exact<T> > type; };
+template<class T> struct opp_policy
+{ typedef save_state<rounded_arith_exact<T> > type; };
+#else
+template<class T> struct std_policy
+{ typedef save_state<rounded_arith_std<T> > type; };
+template<class T> struct opp_policy
+{ typedef save_state<rounded_arith_opp<T> > type; };
+#endif
+
+template<class T>
+void test()
+{
+  test_policy<typename std_policy<T>::type, T>();
+  test_policy<typename opp_policy<T>::type, T>();
+}
+
+int main()
+{
+  test<float>();
+  test<double>();
+  // IBM double-double long double (106 bits) does not follow the rounding
+  // mode: the libgcc division gives the same result in every mode.
+  if (std::numeric_limits<long double>::digits != 106)
+    test<long double>();
+  test_conv<std_policy<float>::type>();
+  test_conv<opp_policy<float>::type>();
+  return boost::report_errors();
+}

test/sqrt_rounding.cpp

@@ -10,12 +10,14 @@
 // 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).
+// and for an exact square both must be the exact root.  Both calls take the
+// same plain local operand, so the test also fails if the compiler merges
+// the two roots across the mode change (GCC bug 34678).  Only the input
+// arrays are volatile, so that nothing is folded at build time.
 
 #include <boost/numeric/interval.hpp>
 #include <boost/core/lightweight_test.hpp>
+#include <boost/limits.hpp>
 
 template<class Rounding, class T>
 void test_policy()
@@ -28,16 +30,14 @@ void test_policy()
     Rounding rnd;
     T a = non_squares[i];
     T dn = rnd.sqrt_down(a);
-    T b = non_squares[i];
-    T up = rnd.sqrt_up(b);
+    T up = rnd.sqrt_up(a);
     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 up = rnd.sqrt_up(a);
     T x = squares[i];
     BOOST_TEST(dn == up);
     BOOST_TEST(dn * dn == x);
@@ -56,10 +56,9 @@ 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))
+  // IBM double-double long double (106 bits) does not follow the rounding
+  // mode.
+  if (std::numeric_limits<long double>::digits != 106)
     test<long double>();
   return boost::report_errors();
 }