Chris@16: /* Boost interval/arith.hpp template implementation file
Chris@16:  *
Chris@16:  * Copyright 2000 Jens Maurer
Chris@16:  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
Chris@16:  *
Chris@16:  * Distributed under the Boost Software License, Version 1.0.
Chris@16:  * (See accompanying file LICENSE_1_0.txt or
Chris@16:  * copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  */
Chris@16: 
Chris@16: #ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP
Chris@16: #define BOOST_NUMERIC_INTERVAL_ARITH_HPP
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/numeric/interval/interval.hpp>
Chris@16: #include <boost/numeric/interval/detail/bugs.hpp>
Chris@16: #include <boost/numeric/interval/detail/test_input.hpp>
Chris@16: #include <boost/numeric/interval/detail/division.hpp>
Chris@16: #include <algorithm>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace numeric {
Chris@16: 
Chris@16: /*
Chris@16:  * Basic arithmetic operators
Chris@16:  */
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: const interval<T, Policies>& operator+(const interval<T, Policies>& x)
Chris@16: {
Chris@16:   return x;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator-(const interval<T, Policies>& x)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   return interval<T, Policies>(-x.upper(), -x.lower(), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator+=(const interval<T, Policies>& r)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(*this, r))
Chris@16:     set_empty();
Chris@16:   else {
Chris@16:     typename Policies::rounding rnd;
Chris@16:     set(rnd.add_down(low, r.low), rnd.add_up(up, r.up));
Chris@16:   }
Chris@16:   return *this;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator+=(const T& r)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(*this, r))
Chris@16:     set_empty();
Chris@16:   else {
Chris@16:     typename Policies::rounding rnd;
Chris@16:     set(rnd.add_down(low, r), rnd.add_up(up, r));
Chris@16:   }
Chris@16:   return *this;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator-=(const interval<T, Policies>& r)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(*this, r))
Chris@16:     set_empty();
Chris@16:   else {
Chris@16:     typename Policies::rounding rnd;
Chris@16:     set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low));
Chris@16:   }
Chris@16:   return *this;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator-=(const T& r)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(*this, r))
Chris@16:     set_empty();
Chris@16:   else {
Chris@16:     typename Policies::rounding rnd;
Chris@16:     set(rnd.sub_down(low, r), rnd.sub_up(up, r));
Chris@16:   }
Chris@16:   return *this;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator*=(const interval<T, Policies>& r)
Chris@16: {
Chris@16:   return *this = *this * r;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator*=(const T& r)
Chris@16: {
Chris@16:   return *this = r * *this;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator/=(const interval<T, Policies>& r)
Chris@16: {
Chris@16:   return *this = *this / r;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies>& interval<T, Policies>::operator/=(const T& r)
Chris@16: {
Chris@16:   return *this = *this / r;
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator+(const interval<T, Policies>& x,
Chris@16:                                 const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   return interval<T,Policies>(rnd.add_down(x.lower(), y.lower()),
Chris@16:                               rnd.add_up  (x.upper(), y.upper()), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   return interval<T,Policies>(rnd.add_down(x, y.lower()),
Chris@16:                               rnd.add_up  (x, y.upper()), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)
Chris@16: { return y + x; }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator-(const interval<T, Policies>& x,
Chris@16:                                 const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   return interval<T,Policies>(rnd.sub_down(x.lower(), y.upper()),
Chris@16:                               rnd.sub_up  (x.upper(), y.lower()), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   return interval<T,Policies>(rnd.sub_down(x, y.upper()),
Chris@16:                               rnd.sub_up  (x, y.lower()), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   return interval<T,Policies>(rnd.sub_down(x.lower(), y),
Chris@16:                               rnd.sub_up  (x.upper(), y), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator*(const interval<T, Policies>& x,
Chris@16:                                 const interval<T, Policies>& y)
Chris@16: {
Chris@16:   BOOST_USING_STD_MIN();
Chris@16:   BOOST_USING_STD_MAX();
Chris@16:   typedef interval<T, Policies> I;
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return I::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   const T& xl = x.lower();
Chris@16:   const T& xu = x.upper();
Chris@16:   const T& yl = y.lower();
Chris@16:   const T& yu = y.upper();
Chris@16: 
Chris@16:   if (interval_lib::user::is_neg(xl))
Chris@16:     if (interval_lib::user::is_pos(xu))
Chris@16:       if (interval_lib::user::is_neg(yl))
Chris@16:         if (interval_lib::user::is_pos(yu)) // M * M
Chris@16:           return I(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)),
Chris@16:                    max BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_up  (xl, yl), rnd.mul_up  (xu, yu)), true);
Chris@16:         else                    // M * N
Chris@16:           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
Chris@16:       else
Chris@16:         if (interval_lib::user::is_pos(yu)) // M * P
Chris@16:           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
Chris@16:         else                    // M * Z
Chris@16:           return I(static_cast<T>(0), static_cast<T>(0), true);
Chris@16:     else
Chris@16:       if (interval_lib::user::is_neg(yl))
Chris@16:         if (interval_lib::user::is_pos(yu)) // N * M
Chris@16:           return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
Chris@16:         else                    // N * N
Chris@16:           return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
Chris@16:       else
Chris@16:         if (interval_lib::user::is_pos(yu)) // N * P
Chris@16:           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
Chris@16:         else                    // N * Z
Chris@16:           return I(static_cast<T>(0), static_cast<T>(0), true);
Chris@16:   else
Chris@16:     if (interval_lib::user::is_pos(xu))
Chris@16:       if (interval_lib::user::is_neg(yl))
Chris@16:         if (interval_lib::user::is_pos(yu)) // P * M
Chris@16:           return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
Chris@16:         else                    // P * N
Chris@16:           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
Chris@16:       else
Chris@16:         if (interval_lib::user::is_pos(yu)) // P * P
Chris@16:           return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
Chris@16:         else                    // P * Z
Chris@16:           return I(static_cast<T>(0), static_cast<T>(0), true);
Chris@16:     else                        // Z * ?
Chris@16:       return I(static_cast<T>(0), static_cast<T>(0), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y)
Chris@16: { 
Chris@16:   typedef interval<T, Policies> I;
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return I::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   const T& yl = y.lower();
Chris@16:   const T& yu = y.upper();
Chris@16:   // x is supposed not to be infinite
Chris@16:   if (interval_lib::user::is_neg(x))
Chris@16:     return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true);
Chris@16:   else if (interval_lib::user::is_zero(x))
Chris@16:     return I(static_cast<T>(0), static_cast<T>(0), true);
Chris@16:   else
Chris@16:     return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)
Chris@16: { return y * x; }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator/(const interval<T, Policies>& x,
Chris@16:                                 const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   if (zero_in(y))
Chris@16:     if (!interval_lib::user::is_zero(y.lower()))
Chris@16:       if (!interval_lib::user::is_zero(y.upper()))
Chris@16:         return interval_lib::detail::div_zero(x);
Chris@16:       else
Chris@16:         return interval_lib::detail::div_negative(x, y.lower());
Chris@16:     else
Chris@16:       if (!interval_lib::user::is_zero(y.upper()))
Chris@16:         return interval_lib::detail::div_positive(x, y.upper());
Chris@16:       else
Chris@16:         return interval<T, Policies>::empty();
Chris@16:   else
Chris@16:     return interval_lib::detail::div_non_zero(x, y);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator/(const T& x, const interval<T, Policies>& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   if (zero_in(y))
Chris@16:     if (!interval_lib::user::is_zero(y.lower()))
Chris@16:       if (!interval_lib::user::is_zero(y.upper()))
Chris@16:         return interval_lib::detail::div_zero<T, Policies>(x);
Chris@16:       else
Chris@16:         return interval_lib::detail::div_negative<T, Policies>(x, y.lower());
Chris@16:     else
Chris@16:       if (!interval_lib::user::is_zero(y.upper()))
Chris@16:         return interval_lib::detail::div_positive<T, Policies>(x, y.upper());
Chris@16:       else
Chris@16:         return interval<T, Policies>::empty();
Chris@16:   else
Chris@16:     return interval_lib::detail::div_non_zero(x, y);
Chris@16: }
Chris@16: 
Chris@16: template<class T, class Policies> inline
Chris@16: interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y)
Chris@16: {
Chris@16:   if (interval_lib::detail::test_input(x, y) || interval_lib::user::is_zero(y))
Chris@16:     return interval<T, Policies>::empty();
Chris@16:   typename Policies::rounding rnd;
Chris@16:   const T& xl = x.lower();
Chris@16:   const T& xu = x.upper();
Chris@16:   if (interval_lib::user::is_neg(y))
Chris@16:     return interval<T, Policies>(rnd.div_down(xu, y), rnd.div_up(xl, y), true);
Chris@16:   else
Chris@16:     return interval<T, Policies>(rnd.div_down(xl, y), rnd.div_up(xu, y), true);
Chris@16: }
Chris@16: 
Chris@16: } // namespace numeric
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP