Why does boost::numeric::interval::widen not behave the same as manually applying the same logic

Root Cause

The problem is floating point inexact representation. If we simplify the test program and use a decimal representation, there is no issue:

Live On Coliru

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/numeric/interval.hpp>
#include <boost/numeric/interval/io.hpp>
#include <iostream>
#include <limits>

int main() {
    using T = boost::multiprecision::cpp_dec_float_50;
    using Interval = boost::numeric::interval<T>;
    using LIM      = std::numeric_limits<T>;

    std::cout << "Precision: " << LIM::max_digits10 << "\n"
              << std::setprecision(LIM::max_digits10);

    const T tol("0.1"); // 1e-12;

    for (int j = 0; j < 100; j++) {
        T lower = -j;
        T upper = j;

        Interval const manual(lower - tol, upper + tol);
        Interval       library{lower, upper};
        library = widen(library, tol);

        using namespace boost::numeric::interval_lib::compare::lexicographic;

        if (library != manual)
            std::cout << " error " << manual << " vs " << library << "\n";
    }
    std::cout << "Done\n";
}

Prints

Precision: 80
Done

Interval Policies & Rounding Modes

It looks like the naive manual implementation results in boundaries that are closest to the intended representation if the representation cannot be exact.

By contrast it looks like widen makes sure that if the result cannot be represented exactly, the widening is guaranteed to be at least the requested amount, never accidentally slightly less due to representation issues. This means that in these situations

• the widened lower bound may be lower than the manually calculated one
• the widened upper bound may be higher than the manually calculated one

The behaviour can be adjusted, because it comes from the rounding policy which supplies sub_down and add_up primitives.

I would suggest not to change it, though:

if you need precise computations with the float or double types, use the default rounded_math;

rounded_math<T> is already the default.

If you insist on having identical results as the manual method (even if they are inferior):

using T        = double; // boost::multiprecision::cpp_dec_float_50;

namespace I    = boost::numeric::interval_lib;
using Interval = I::change_rounding<boost::numeric::interval<T>,
                                    I::rounded_arith_exact<T>>::type;

Prints Live On Coliru

Precision: 17
Done

Further Caveats

There's a documentation warning:

Warning! Guaranteed interval arithmetic for native floating-point format is not supported on every combination of processor, operating system, and compiler.

There is a list of specific compiler quirks/flags you have to keep in mind.

TL;DR Boost is complicated and somewhere it sets the mode of the floating-point unit (FPU) which changes rounding behavior.

The answer by sehe is spot on but in case you are interested in a more technical answer then read on.

The lines of interest in your example are

  using Interval = boost::numeric::interval<double>;
  Interval i{lower, upper};
  Interval i_ = widen(i, -tol);

Let's go down the rabbit hole and try to understand how boost implements these lines.

The main class is

// boost/numeric/interval/interval.hpp:l.37
template<class T, class Policies>
class interval;

The template argument Policies gets its default value defined at

// boost/numeric/interval/detail/interval_prototype.hpp:l.27
template<class T>
struct default_policies
{
  typedef policies<rounded_math<T>, checking_strict<T> > type;
};
// <truncated>
template<class T, class Policies = typename interval_lib::default_policies<T>::type >
class interval;

The general rounded_math template is defined at

// boost/numeric/interval/rounding.hpp:l.93
template<class T>
struct rounded_math: save_state_nothing<rounded_arith_exact<T> >
{};

but we actually need to look at the explicit template specialization

// boost/numeric/interval/hw_rounding.hpp:l.59
template<>
struct rounded_math<double>
  : save_state<rounded_arith_opp<double> >
{};

The save_state template is a derived class from its argument, i.e.

// boost/numeric/interval/rounding.hpp:l.75
template<class Rounding>
struct save_state: Rounding
{
  typename Rounding::rounding_mode mode;
  save_state() {
    this->get_rounding_mode(mode);
    this->init();
  }
  ~save_state() { this->set_rounding_mode(mode); }
  typedef detail::save_state_unprotected<Rounding> unprotected_rounding;
};

Note here that Rounding::init() will be called (as well as Rounding::set_rounding_mode(mode) at destruction)!

Then lets have a look at rounded_arith_opp

// boost/numeric/interval/rounded_arith.hpp:l.78
template<class T, class Rounding>
struct rounded_arith_opp: Rounding {
  void init() { this->upward(); }
  // <truncated>
# define BOOST_UP(EXPR) return this->force_rounding(EXPR)
# define BOOST_UP_NEG(EXPR) return -this->force_rounding(EXPR)
  T sub_down(const T& x, const T& y) { BOOST_UP_NEG(y - x); }
  T add_up  (const T& x, const T& y) { BOOST_UP(x + y); }
  // <truncated>
};

Its default Rounding argument is defined at

// boost/numeric/interval/rounding.hpp:l.46
template<class T, class Rounding = rounding_control<T> >
struct rounded_arith_opp;

The rounding control is defined at

// boost/numeric/interval/rounding.hpp:l.20
template<class T>
struct rounding_control
{
  // <truncated>
  static const T& force_rounding(const T& x) { return x; }
};

It also has explicit template specializations which vary by platforms, e.g.

// boost/numeric/interval/detail/c99_rounding_control.hpp:l.28
template<>
struct rounding_control<double>:
  detail::c99_rounding_control
{
  static double force_rounding(double const &r)
  { volatile double r_ = r; return r_; }
};

The c99_rounding_control class implements the change of rounding modes, i.e.

// boost/numeric/interval/detail/c99sub_rounding_control.hpp:l.23
struct c99_rounding_control
{
  static void set_rounding_mode(rounding_mode mode)  { fesetround(mode); }
  static void get_rounding_mode(rounding_mode &mode) { mode = fegetround(); }
  static void upward()      { set_rounding_mode(FE_UPWARD);     }
  // <truncated>
};

Now we understand the interval class and can have a look at the widen function

// boost/numeric/interval/detail/c99_rounding_control.hpp:l.82
template<class T, class Policies> inline
interval<T, Policies> widen(const interval<T, Policies>& x, const T& v)
{
  if (interval_lib::detail::test_input(x))
    return interval<T, Policies>::empty();
  typename Policies::rounding rnd;
  return interval<T, Policies>(rnd.sub_down(x.lower(), v),
                               rnd.add_up  (x.upper(), v), true);
}

It creates the rnd object which will call in the end c99_rounding_control::upward and together with rounded_arith_opp::{add_up,sub_down} one gets the observed results.

I will supply a simple example to show how the rounding takes place (godbolt-link)

#include <boost/numeric/interval.hpp>
#include <iomanip>
#include <iostream>
#include <limits>

using namespace std;

double force_rounding(const double& r) {
  volatile double r_ = r;
  return r_;
}

int main() {
  cout << setprecision(numeric_limits<double>::max_digits10) << endl;

  auto tol = 0.1;
  const auto lower = -2;
  const auto upper = 3;

  cout << "default rounding:\t" << lower + tol << ", " << upper - tol << '\n';

  {
    using Interval = boost::numeric::interval<double>;
    const auto i = widen(Interval{lower, upper}, -tol);
    cout << "widen:\t\t\t" << i.lower() << ", " << i.upper() << '\n';
  }

  {
    boost::numeric::interval_lib::detail::c99_rounding_control::upward();
    const auto lower_ = -force_rounding(-lower - tol);
    const auto upper_ = force_rounding(upper - tol);
    cout << "set rounding:\t\t" << lower_ << ", " << upper_ << '\n';
  }

  return 0;
}

Output:

default rounding:       -1.8999999999999999, 2.8999999999999999
widen:                  -1.9000000000000001, 2.9000000000000004
set rounding:           -1.9000000000000001, 2.9000000000000004

See also

• boost doc about interval's rounding: boost.org
• boost version: 1.79.0 (for line numbers)