sqrt of a point interval can lose its upper end with g++ -O2 -frounding-math on x86-64 (1.92.0 and develop) #60

With the default policy, sqrt(interval<double>(x)) for a point interval can return an interval whose upper end is below the true square root. Seen with g++ 13.3.0 at -O2 -frounding-math (also -O3) on x86-64 Linux, with Boost 1.92.0 and with develop at ef305d5.

Reproducer

#include <boost/numeric/interval.hpp>
#include <cmath>
#include <cstdio>

typedef boost::numeric::interval<double> I;
volatile double inputs[2] = {3.0, 2.0};

int main() {
  int lost = 0;
  for (int i = 0; i < 2; ++i) {
    double x = inputs[i];
    I y = sqrt(I(x));
    bool ok = std::fma(y.upper(), y.upper(), -x) >= 0;  // exact sign of u*u - x
    std::printf("sqrt(%g): [%a, %a] %s\n", x, y.lower(), y.upper(),
                ok ? "encloses" : "UPPER END BELOW THE TRUE ROOT");
    if (!ok) ++lost;
  }
  return lost;
}

g++ -std=c++17 -O2 -frounding-math prints:

sqrt(3): [0x1.bb67ae8584caap+0, 0x1.bb67ae8584caap+0] UPPER END BELOW THE TRUE ROOT
sqrt(2): [0x1.6a09e667f3bccp+0, 0x1.6a09e667f3bccp+0] UPPER END BELOW THE TRUE ROOT

The expected result for 3 is [0x1.bb67ae8584caap+0, 0x1.bb67ae8584cabp+0]. That is what g++ -O0 and clang 18.1.3 at -O0, -O2 and -O3 print. Code shape matters. A single straight-line call did not miss at plain -O2 -frounding-math, because one of the two roots went through the C library sqrt call (the errno path). Adding -fno-math-errno (which -ffast-math also implies) makes that single call miss too.

Cause

rounded_arith_opp::sqrt_down sets downward rounding and evaluates sqrt(x). sqrt_up evaluates sqrt(x) again, relying on upward mode. For a point interval both receive the same value. Inside a loop, g++ emits one sqrtsd under downward mode and uses its result for both ends. The mode-change calls still happen, but the second root is never computed. The volatile store in force_rounding protects the result, not the operand. rounded_arith_std has the same shape and fails the same way. This is GCC bug 34678 (a value reused across a rounding-mode change), reached through this library.

Wider test. 96 point inputs (14 fixed values, including four exact squares, and 82 seeded values in [1, 4)), with each end checked exactly:

The existing tests mul.cpp and test_float.cpp pass in the failing builds. Their square-root checks use wide intervals, where the two ends have different operands.

At g++ -O2 and -O3, rounded_arith_std also loses point add, subtract, multiply and divide: 44, 10, 60 and 62 misses out of 64 pairs. The default policy computes the lower end with a negated expression and has 0 misses there.

Proposed fix (one store and one load per root). After setting the mode, read the operand back through a volatile local, so that the two roots cannot be merged:

namespace detail {
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; }
inline double reload_operand(const double& x) { volatile double v = x; return v; }
inline long double reload_operand(const long double& x) { volatile long double v = x; return v; }
}
// rounded_arith_opp
T sqrt_down(const T& x)
{ BOOST_NUMERIC_INTERVAL_using_math(sqrt);
  this->downward(); T y = detail::reload_operand(x);
  T r = this->force_rounding(sqrt(y)); this->upward(); return r; }
T sqrt_up(const T& x)
{ BOOST_NUMERIC_INTERVAL_using_math(sqrt);
  T y = detail::reload_operand(x); return this->force_rounding(sqrt(y)); }

rounded_arith_std gets the same change. A unified diff applies cleanly to 1.92.0, to develop and to PR 58. With it, there are 0 misses in every build above, on 1.92.0 and on develop. A second, optional diff does the same for the eight add, subtract, multiply and divide functions of rounded_arith_std. With both diffs, there are 0 misses at g++ and clang -O0, -O2 and -O3.

Limits

• One machine, two compilers, 96 inputs.
• A volatile read stops this merge, but it is not a language guarantee. The real fix is compiler support for the dynamic rounding mode.
• The released documentation lists x86 support with the x87 unit, so SSE could be read as outside what was promised. As I read it, the open PR 58 would extend support to SSE; if it is merged, this case falls inside it.
• I searched open and closed issues and pull requests for "sqrt". The closest is test_sqrt function in mul.cpp got failed. Need to change rounding mode for MXCSR. #23 (test_sqrt failing because MXCSR was not set), which Set the x86 rounding mode with fesetround (fixes #23, #35, #41) #58 addresses with fesetround. This is a different cause: here the mode is set, but the optimiser reuses one root for both ends. With Set the x86 rounding mode with fesetround (fixes #23, #35, #41) #58 applied, the same 92 of 96 misses remain at g++ -O2. nth_root does not restore the rounding mode and can return an upper bound below the true root #54 (nth_root returning an upper bound below the true root) is the same symptom in another function.

What I can send: the two diffs, a regression test using lightweight_test.hpp (it includes an exact-policy control that must fail, and does), the probe for the four operations, and the full run log. A pull request with the fix and the test follows and links here.

Related: #23 and #54; GCC bug 34678; Trac ticket 12018 and pull request 4 (the same kind of merge in rounded_arith_opp add, multiply and divide under MSVC and clang); PR 58.

Sonny Guardian