#ifndef BOOST_NUMERIC_INTERVAL_HPP #define BOOST_NUMERIC_INTERVAL_HPP // MS compatible compilers support #pragma once #if defined(_MSC_VER) && (_MSC_VER >= 1020) # pragma once #endif // Copyright (c) 2012 Robert Ramey // // 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) #include #include #include #include #include #include #include #include "utility.hpp" // log #include "checked_result.hpp" #include "checked.hpp" // from stack overflow // http://stackoverflow.com/questions/23815138/implementing-variadic-min-max-functions namespace boost { namespace numeric { template struct interval { static_assert( std::is_literal_type< checked_result >::value, "is literal type" ); const R l; const R u; template constexpr interval(const T & lower, const T & upper) : l(lower), u(upper) {} template constexpr interval(const std::pair & p) : l(p.first), u(p.second) {} template constexpr interval(const interval & rhs) : l(rhs.l), u(rhs.u) {} constexpr interval(); // return true if this interval contains every point found in some // other inteval t template constexpr bool includes(const interval & t) const { return safe_compare::greater_than_equal(t.l, l) && safe_compare::less_than_equal(t.u, u) ; } template constexpr bool includes(const T & t) const { return safe_compare::greater_than_equal(t, l) && safe_compare::less_than_equal(t, u) ; } }; template constexpr interval make_interval(const R & r){ return interval(); } template constexpr interval::interval() : l(std::numeric_limits::lowest()), u(std::numeric_limits::max()) {} // account for the fact that for floats and doubles // the most negative value is called "lowest" rather // than min template<> constexpr interval::interval() : l(std::numeric_limits::lowest()), u(std::numeric_limits::max()) {} template<> constexpr interval::interval() : l(std::numeric_limits::lowest()), u(std::numeric_limits::max()) {} namespace detail { template constexpr interval minmax(const std::initializer_list & l){ using namespace boost::numeric; T minimum{exception_type::positive_overflow_error, ""}; T maximum{exception_type::negative_overflow_error, ""}; // note: we can't use for_each and a lambda because neither of these are // constexpr for( typename std::initializer_list::const_iterator i = l.begin(); i != l.end(); ++i ){ // std::cout << *i << ','; // std::cout.flush(); if(minimum != exception_type::negative_overflow_error){ // if it corresponds to the lowest value if(*i == exception_type::negative_overflow_error){ // initialize the minimum minimum = *i; } // skip any other exceptions else if(! (*i).exception()){ if(minimum.exception() || (*i).m_r < minimum.m_r){ minimum = *i; } } } if(maximum != exception_type::positive_overflow_error){ // if it corresponds to the highest value if(*i == exception_type::positive_overflow_error){ // initialize the maximum maximum = *i; } // skip any other exceptions else if(! (*i).exception()){ if(maximum.exception() || (*i).m_r > maximum.m_r){ maximum = *i; } } } } return interval{minimum, maximum}; } } // detail template constexpr interval> add( const interval & t, const interval & u ){ // adapted from https://en.wikipedia.org/wiki/Interval_arithmetic const checked_result lower = checked::add(t.l, u.l); const checked_result upper = checked::add(t.u, u.u); return interval>(lower, upper); } template constexpr interval> subtract( const interval & t, const interval & u ){ // adapted from https://en.wikipedia.org/wiki/Interval_arithmetic const checked_result lower = checked::subtract(t.l, u.u); const checked_result upper = checked::subtract(t.u, u.l); return interval>(lower, upper); } template constexpr interval> multiply( const interval & t, const interval & u ){ // adapted from https://en.wikipedia.org/wiki/Interval_arithmetic return detail::minmax>( std::initializer_list> { checked::multiply(t.l, u.l), checked::multiply(t.l, u.u), checked::multiply(t.u, u.l), checked::multiply(t.u, u.u) } ); } // divide two intervals. This will give a new range of the quotient. But // it doesn't consider the possibility that the divisor is zero. This will // have to be considered separately. template constexpr interval> divide( const interval & t, const interval & u ){ // assert(u.l != 0 || u.u != 0); return detail::minmax>( (u.l == 0) ? std::initializer_list> { checked::divide(t.l, +1), checked::divide(t.l, u.u), checked::divide(t.u, +1), checked::divide(t.u, u.u) } : (0 == u.u) ? std::initializer_list> { checked::divide(t.l, u.l), checked::divide(t.l, -1), checked::divide(t.u, u.l), checked::divide(t.u, -1) } : (u.u < 0 || 0 < u.l) // divisor range excludes 0 ? std::initializer_list> { checked::divide(t.l, u.l), checked::divide(t.l, u.u), checked::divide(t.u, u.l), checked::divide(t.u, u.u) } : std::initializer_list> { // divisor includes 0 checked::divide(t.l, u.l), checked::divide(t.l, -1), checked::divide(t.l, +1), checked::divide(t.l, u.u), checked::divide(t.u, u.l), checked::divide(t.u, -1), checked::divide(t.u, +1), checked::divide(t.u, u.u) } ); } // modulus of two intervals. This will give a new range of for the modulus. But // it doesn't consider the possibility that the radix is zero. This will // have to be considered separately. template constexpr interval> modulus( const interval & t, const interval & u ){ /* Turns out that getting the exact range is a suprisingly difficult * problem. But we can get a reasonable guaranteed range with the following * simple formula (due to Skona Brittain): * * [-m + 1, m - 1] where m is max(|c|, |d|). * * where the operation is [a,b] % [c,d] * * But a naive implementation of this algorithm doesn't work. We have a * problem with (8 bit integers for illustration) abs(-128) -> -128! So we * can't get the absolute value of the minimum value ! * * so we impement the above in a slightly different way. */ // const R xa = (t.l > 0) ? t.l - 1 : t.l + 1; // a // const R xb = (t.u > 0) ? t.u - 1 : t.u + 1; // b const checked_result xc = (u.l > 0) ? u.l - 1 : (u.l < 0) ? -(u.l + 1) : 0; // c const checked_result xd = (u.u > 0) ? u.u - 1 : (u.u < 0) ? -(u.u + 1) : 0; // d // special care to void problem when inverting -128 const checked_result mxc = checked::subtract(0, xc); const checked_result mxd = checked::subtract(0, xd); return detail::minmax>( std::initializer_list> { xc, mxc, xd, mxd, checked_result(0) } ); } template constexpr interval> left_shift( const interval & t, const interval & u ){ return interval>{ checked::left_shift(t.l, u.l), checked::left_shift(t.u, u.u), }; } template constexpr interval> right_shift( const interval & t, const interval & u ){ return interval>{ checked::right_shift(t.l, u.u), checked::right_shift(t.u, u.l) }; } template constexpr interval> bitwise_or( const interval & t, const interval & u ){ return detail::minmax>( std::initializer_list> { checked::bitwise_or(t.l, u.l), checked::bitwise_or(t.l, u.u), checked::bitwise_or(t.u, u.l), checked::bitwise_or(t.u, u.u) } ); } template constexpr interval> bitwise_and( const interval & t, const interval & u ){ return detail::minmax>( std::initializer_list> { checked::bitwise_and(t.l, u.l), checked::bitwise_and(t.l, u.u), checked::bitwise_and(t.u, u.l), checked::bitwise_and(t.u, u.u) } ); } template constexpr interval> bitwise_xor( const interval & t, const interval & u ){ return detail::minmax>( std::initializer_list> { checked::bitwise_xor(t.l, u.l), checked::bitwise_xor(t.l, u.u), checked::bitwise_xor(t.u, u.l), checked::bitwise_xor(t.u, u.u) } ); } template constexpr interval> intersection( const interval & t, const interval & u ){ const checked_result rl = checked::cast(safe_compare::greater_than(t.l, u.l) ? t.l : u.l); const checked_result ru = checked::cast(safe_compare::less_than(t, u) ? t.u : u.u); if(rl > ru){ return checked_result>( exception_type::uninitialized, "null intersection" ); } return interval(rl, ru); } template constexpr interval> union_interval( const interval & t, const interval & u ){ const checked_result rl = checked::cast(safe_compare::less_than(t.l, u.l) ? t.l : u.l); const checked_result ru = checked::cast(safe_compare::greater_than(t, u) ? t.u : u.u); return interval(rl, ru); } template constexpr boost::logic::tribool operator<( const interval & t, const interval & u ){ return // if every element in t is less than every element in u safe_compare::less_than(t.u, u.l) ? boost::logic::tribool(true) : // if every element in t is greater than every element in u safe_compare::greater_than(t.l, u.u) ? boost::logic::tribool(false) : // otherwise some element(s) in t are greater than some element in u boost::logic::indeterminate ; } template constexpr boost::logic::tribool operator>( const interval & t, const interval & u ){ return // if every element in t is greater than every element in u safe_compare::greater_than(t.l, u.u) ? boost::logic::tribool(true) : // if every element in t is less than every element in u safe_compare::less_than(t.u, u.l) ? boost::logic::tribool(false) : // otherwise some element(s) in t are greater than some element in u boost::logic::indeterminate ; } template constexpr bool operator==( const interval & t, const interval & u ){ // intervals have the same limits return safe_compare::equal(t.l, u.l) && safe_compare::equal(t.u, u.u) ; } template constexpr bool operator!=( const interval & t, const interval & u ){ return ! (t == u); } template constexpr boost::logic::tribool operator<=( const interval & t, const interval & u ){ return ! (t > u); } template constexpr boost::logic::tribool operator>=( const interval & t, const interval & u ){ return ! (t < u); } } // numeric } // boost #include namespace std { template inline std::basic_ostream & operator<<( std::basic_ostream & os, const boost::numeric::interval & i ){ return os << '[' << i.l << ',' << i.u << ']'; } template inline std::basic_ostream & operator<<( std::basic_ostream & os, const boost::numeric::interval & i ){ os << "[" << (unsigned)i.l << "," << (unsigned)i.u << "]"; return os; } template inline std::basic_ostream & operator<<( std::basic_ostream & os, const boost::numeric::interval & i ){ os << "[" << (int)i.l << "," << (int)i.u << "]"; return os; } } // std #endif // BOOST_NUMERIC_INTERVAL_HPP