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