Boost Interval with Multiprecision

Question

I try to use the boost interval arithmetic library together with the boost multiprecision library. If I use standard double precision with the native double datatype, everything works fine.

With the multiprecision library however, it produces results that are actually more inaccurate. Here some code:

#include <boost\numeric\interval.hpp>
#include <boost\multiprecision\cpp_dec_float.hpp>
#include <iostream>

using namespace boost::numeric;
using namespace boost::numeric::interval_lib;
using namespace boost::multiprecision;

template <typename T>
using Interval = interval<T, policies<save_state<rounded_transc_exact<T>>, checking_base<T>>>;
using BigFloat = cpp_dec_float_100;

int main()
{
    std::cout << sin(Interval<double>(0.0, 0.1)).upper() << "\n"; // 0.0998334
    std::cout << sin(Interval<BigFloat>(0.0, 0.1)).upper() << "\n"; // 1
}

As can be seen, the double version produces a very accurate result. The BigFloat version should be even much more accurate, however it produces a very large bound - actually the maximum value of the sin function, so this bound is completly useless.

How can I fix this such that the interval library actually takes advantage of the higher precision and produces sharper bounds?

Answer 1

To get started, I tested with cos instead of sin.

The interval library implements sin(x) as in terms of cos(x-½π). This means that sin([0, 0.1]) is transformed into cos([-½π,-½π+0.1]) (which recurses into cos([½π,½π+0.1])).

In the case of BigFloat, and due to the library not know Pi constants (pi<BigFloat>(), pi_half<BigFloat>() or pi_twice<BigFloat>()) it represents them as integer intervals, e.g.: pi_half<BigFloat> is represented as [1,2]. OOPS. the cos intervals have become [-2,-0.9] (recursing into [0,3.1]¹).

Adding some tracing:

DOUBLE--------------------
pi/2: [1.570796326794896558,1.57079632679489678]
sin: [0,0.10000000000000000555]
cos: [-1.57079632679489678,-1.4707963267948964692]
cos: [1.5707963267948961139,1.6707963267948979791]
[-5.0532154980743028885e-16,0.099833416646829500896]
BigFloat--------------------
pi/2: [1,2]
sin: [0,0.10000000000000000555]
cos: [-2,-0.89999999999999999445]
cos: [0,3.1000000000000000056]
[-1,1]

Solutions?

The best solution I can think of involves using cos directly or specializing pi_half:

Use cos directly

This is NOT a solution because it will still use some of the broken pi_*<BigFloat>() constants internally:

static BigFloat bf_pi_half() { return bmp::default_ops::get_constant_pi<BigFloat::backend_type>() / BigFloat(2); }

Now you could write

Live On Coliru

std::cout << "BigFloat--------------------\n";
std::cout << cos(ival - bf_pi_half()) << "\n";

Printing

BigFloat--------------------
[-0.909297,0.818277]

As you can see, this is not the output desired.

Specialize Constants

In fact, you should specialize the underlying constants:

Live On Coliru

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

namespace bn  = boost::numeric;
namespace bni = bn::interval_lib;
namespace bmp = boost::multiprecision;

template <typename T>
using Interval = bn::interval<T, bni::policies<bni::save_state<bni::rounded_transc_exact<T>>, bni::checking_base<T>>>;
using BigFloat = bmp::cpp_dec_float_100; // bmp::number<bmp::backends::cpp_dec_float<100>, bmp::et_off>;

static BigFloat bf_pi() { return bmp::default_ops::get_constant_pi<BigFloat::backend_type>(); }

namespace boost { namespace numeric { namespace interval_lib { namespace constants {
    template<> inline BigFloat pi_lower<BigFloat>()       { return bf_pi(); }
    template<> inline BigFloat pi_upper<BigFloat>()       { return bf_pi(); }
    template<> inline BigFloat pi_twice_lower<BigFloat>() { return bf_pi() * 2; }
    template<> inline BigFloat pi_twice_upper<BigFloat>() { return bf_pi() * 2; }
    template<> inline BigFloat pi_half_lower<BigFloat>()  { return bf_pi() / 2; }
    template<> inline BigFloat pi_half_upper<BigFloat>()  { return bf_pi() / 2; }
} } } }

int main() {
    std::cout << sin(Interval<BigFloat>(0, 0.1)) << "\n";
}

This prints:

[0,0.0998334]

---

---

¹ that's actually using pi_twice<> constant