diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9f535e8..81a25db 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/include/boost/numeric/interval/detail/x86_rounding_control.hpp b/include/boost/numeric/interval/detail/x86_rounding_control.hpp index 9078692..05bdcbb 100644 --- a/include/boost/numeric/interval/detail/x86_rounding_control.hpp +++ b/include/boost/numeric/interval/detail/x86_rounding_control.hpp @@ -85,11 +85,20 @@ struct x86_rounding_control_long_double: 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: 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 diff --git a/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp b/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp index 07af267..b62139c 100644 --- a/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp +++ b/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_; } }; diff --git a/include/boost/numeric/interval/rounded_arith.hpp b/include/boost/numeric/interval/rounded_arith.hpp index 2eac74e..79bca32 100644 --- a/include/boost/numeric/interval/rounded_arith.hpp +++ b/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 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 T conv_down(U const &v) { BOOST_DN(v); } - template 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 T conv_down(U const &v) { BOOST_DN(detail::reload_operand(v)); } + template 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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ff38b12..19e29b2 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 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() diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 3c24f39..6a52ff9 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 rounding_merge.cpp ] [ run sqrt_rounding.cpp ] [ run sqrt_point.cpp ] [ run rounding_mode.cpp : : : [ check-target-builds /boost/architecture//power : : no ] ] diff --git a/test/rounding_merge.cpp b/test/rounding_merge.cpp new file mode 100644 index 0000000..18dcffc --- /dev/null +++ b/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 +#include +#include +#include + +using namespace boost::numeric; +using namespace boost::numeric::interval_lib; + +#define BINARY_TEST(OP, EXACT) \ +template \ +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 +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 +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 +void test_policy() +{ + test_add(); + test_sub(); + test_mul(); + test_div(); + test_int(); +} + +#if defined(CONTROL_EXACT_POLICY) +template struct std_policy +{ typedef save_state > type; }; +template struct opp_policy +{ typedef save_state > type; }; +#else +template struct std_policy +{ typedef save_state > type; }; +template struct opp_policy +{ typedef save_state > type; }; +#endif + +template +void test() +{ + test_policy::type, T>(); + test_policy::type, T>(); +} + +int main() +{ + test(); + test(); + // 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::digits != 106) + test(); + test_conv::type>(); + test_conv::type>(); + return boost::report_errors(); +} diff --git a/test/sqrt_rounding.cpp b/test/sqrt_rounding.cpp index 30ff6e3..728357c 100644 --- a/test/sqrt_rounding.cpp +++ b/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 #include +#include template 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(); 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)) + // IBM double-double long double (106 bits) does not follow the rounding + // mode. + if (std::numeric_limits::digits != 106) test(); return boost::report_errors(); }