numeric/interval/arith.hpp

0001 /* Boost interval/arith.hpp template implementation file
0002  *
0003  * Copyright 2000 Jens Maurer
0004  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
0005  *
0006  * Distributed under the Boost Software License, Version 1.0.
0007  * (See accompanying file LICENSE_1_0.txt or
0008  * copy at http://www.boost.org/LICENSE_1_0.txt)
0009  */
0010
0011 #ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP
0012 #define BOOST_NUMERIC_INTERVAL_ARITH_HPP
0013
0014 #include <boost/config.hpp>
0015 #include <boost/numeric/interval/interval.hpp>
0016 #include <boost/numeric/interval/detail/bugs.hpp>
0017 #include <boost/numeric/interval/detail/test_input.hpp>
0018 #include <boost/numeric/interval/detail/division.hpp>
0019 #include <algorithm>
0020
0021 namespace boost {
0022 namespace numeric {
0023
0024 /*
0025  * Basic arithmetic operators
0026  */
0027
0028 template<class T, class Policies> inline
0029 const interval<T, Policies>& operator+(const interval<T, Policies>& x)
0030 {
0031   return x;
0032 }
0033
0034 template<class T, class Policies> inline
0035 interval<T, Policies> operator-(const interval<T, Policies>& x)
0036 {
0037   if (interval_lib::detail::test_input(x))
0038     return interval<T, Policies>::empty();
0039   return interval<T, Policies>(-x.upper(), -x.lower(), true);
0040 }
0041
0042 template<class T, class Policies> inline
0043 interval<T, Policies>& interval<T, Policies>::operator+=(const interval<T, Policies>& r)
0044 {
0045   if (interval_lib::detail::test_input(*this, r))
0046     set_empty();
0047   else {
0048     typename Policies::rounding rnd;
0049     set(rnd.add_down(low, r.low), rnd.add_up(up, r.up));
0050   }
0051   return *this;
0052 }
0053
0054 template<class T, class Policies> inline
0055 interval<T, Policies>& interval<T, Policies>::operator+=(const T& r)
0056 {
0057   if (interval_lib::detail::test_input(*this, r))
0058     set_empty();
0059   else {
0060     typename Policies::rounding rnd;
0061     set(rnd.add_down(low, r), rnd.add_up(up, r));
0062   }
0063   return *this;
0064 }
0065
0066 template<class T, class Policies> inline
0067 interval<T, Policies>& interval<T, Policies>::operator-=(const interval<T, Policies>& r)
0068 {
0069   if (interval_lib::detail::test_input(*this, r))
0070     set_empty();
0071   else {
0072     typename Policies::rounding rnd;
0073     set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low));
0074   }
0075   return *this;
0076 }
0077
0078 template<class T, class Policies> inline
0079 interval<T, Policies>& interval<T, Policies>::operator-=(const T& r)
0080 {
0081   if (interval_lib::detail::test_input(*this, r))
0082     set_empty();
0083   else {
0084     typename Policies::rounding rnd;
0085     set(rnd.sub_down(low, r), rnd.sub_up(up, r));
0086   }
0087   return *this;
0088 }
0089
0090 template<class T, class Policies> inline
0091 interval<T, Policies>& interval<T, Policies>::operator*=(const interval<T, Policies>& r)
0092 {
0093   return *this = *this * r;
0094 }
0095
0096 template<class T, class Policies> inline
0097 interval<T, Policies>& interval<T, Policies>::operator*=(const T& r)
0098 {
0099   return *this = r * *this;
0100 }
0101
0102 template<class T, class Policies> inline
0103 interval<T, Policies>& interval<T, Policies>::operator/=(const interval<T, Policies>& r)
0104 {
0105   return *this = *this / r;
0106 }
0107
0108 template<class T, class Policies> inline
0109 interval<T, Policies>& interval<T, Policies>::operator/=(const T& r)
0110 {
0111   return *this = *this / r;
0112 }
0113
0114 template<class T, class Policies> inline
0115 interval<T, Policies> operator+(const interval<T, Policies>& x,
0116                                 const interval<T, Policies>& y)
0117 {
0118   if (interval_lib::detail::test_input(x, y))
0119     return interval<T, Policies>::empty();
0120   typename Policies::rounding rnd;
0121   return interval<T,Policies>(rnd.add_down(x.lower(), y.lower()),
0122                               rnd.add_up  (x.upper(), y.upper()), true);
0123 }
0124
0125 template<class T, class Policies> inline
0126 interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y)
0127 {
0128   if (interval_lib::detail::test_input(x, y))
0129     return interval<T, Policies>::empty();
0130   typename Policies::rounding rnd;
0131   return interval<T,Policies>(rnd.add_down(x, y.lower()),
0132                               rnd.add_up  (x, y.upper()), true);
0133 }
0134
0135 template<class T, class Policies> inline
0136 interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)
0137 { return y + x; }
0138
0139 template<class T, class Policies> inline
0140 interval<T, Policies> operator-(const interval<T, Policies>& x,
0141                                 const interval<T, Policies>& y)
0142 {
0143   if (interval_lib::detail::test_input(x, y))
0144     return interval<T, Policies>::empty();
0145   typename Policies::rounding rnd;
0146   return interval<T,Policies>(rnd.sub_down(x.lower(), y.upper()),
0147                               rnd.sub_up  (x.upper(), y.lower()), true);
0148 }
0149
0150 template<class T, class Policies> inline
0151 interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y)
0152 {
0153   if (interval_lib::detail::test_input(x, y))
0154     return interval<T, Policies>::empty();
0155   typename Policies::rounding rnd;
0156   return interval<T,Policies>(rnd.sub_down(x, y.upper()),
0157                               rnd.sub_up  (x, y.lower()), true);
0158 }
0159
0160 template<class T, class Policies> inline
0161 interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y)
0162 {
0163   if (interval_lib::detail::test_input(x, y))
0164     return interval<T, Policies>::empty();
0165   typename Policies::rounding rnd;
0166   return interval<T,Policies>(rnd.sub_down(x.lower(), y),
0167                               rnd.sub_up  (x.upper(), y), true);
0168 }
0169
0170 template<class T, class Policies> inline
0171 interval<T, Policies> operator*(const interval<T, Policies>& x,
0172                                 const interval<T, Policies>& y)
0173 {
0174   BOOST_USING_STD_MIN();
0175   BOOST_USING_STD_MAX();
0176   typedef interval<T, Policies> I;
0177   if (interval_lib::detail::test_input(x, y))
0178     return I::empty();
0179   typename Policies::rounding rnd;
0180   const T& xl = x.lower();
0181   const T& xu = x.upper();
0182   const T& yl = y.lower();
0183   const T& yu = y.upper();
0184
0185   if (interval_lib::user::is_neg(xl))
0186     if (interval_lib::user::is_pos(xu))
0187       if (interval_lib::user::is_neg(yl))
0188         if (interval_lib::user::is_pos(yu)) // M * M
0189           return I(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)),
0190                    max BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_up  (xl, yl), rnd.mul_up  (xu, yu)), true);
0191         else                    // M * N
0192           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
0193       else
0194         if (interval_lib::user::is_pos(yu)) // M * P
0195           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
0196         else                    // M * Z
0197           return I(static_cast<T>(0), static_cast<T>(0), true);
0198     else
0199       if (interval_lib::user::is_neg(yl))
0200         if (interval_lib::user::is_pos(yu)) // N * M
0201           return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
0202         else                    // N * N
0203           return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
0204       else
0205         if (interval_lib::user::is_pos(yu)) // N * P
0206           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
0207         else                    // N * Z
0208           return I(static_cast<T>(0), static_cast<T>(0), true);
0209   else
0210     if (interval_lib::user::is_pos(xu))
0211       if (interval_lib::user::is_neg(yl))
0212         if (interval_lib::user::is_pos(yu)) // P * M
0213           return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
0214         else                    // P * N
0215           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
0216       else
0217         if (interval_lib::user::is_pos(yu)) // P * P
0218           return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
0219         else                    // P * Z
0220           return I(static_cast<T>(0), static_cast<T>(0), true);
0221     else                        // Z * ?
0222       return I(static_cast<T>(0), static_cast<T>(0), true);
0223 }
0224
0225 template<class T, class Policies> inline
0226 interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y)
0227 {
0228   typedef interval<T, Policies> I;
0229   if (interval_lib::detail::test_input(x, y))
0230     return I::empty();
0231   typename Policies::rounding rnd;
0232   const T& yl = y.lower();
0233   const T& yu = y.upper();
0234   // x is supposed not to be infinite
0235   if (interval_lib::user::is_neg(x))
0236     return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true);
0237   else if (interval_lib::user::is_zero(x))
0238     return I(static_cast<T>(0), static_cast<T>(0), true);
0239   else
0240     return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true);
0241 }
0242
0243 template<class T, class Policies> inline
0244 interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)
0245 { return y * x; }
0246
0247 template<class T, class Policies> inline
0248 interval<T, Policies> operator/(const interval<T, Policies>& x,
0249                                 const interval<T, Policies>& y)
0250 {
0251   if (interval_lib::detail::test_input(x, y))
0252     return interval<T, Policies>::empty();
0253   if (zero_in(y))
0254     if (!interval_lib::user::is_zero(y.lower()))
0255       if (!interval_lib::user::is_zero(y.upper()))
0256         return interval_lib::detail::div_zero(x);
0257       else
0258         return interval_lib::detail::div_negative(x, y.lower());
0259     else
0260       if (!interval_lib::user::is_zero(y.upper()))
0261         return interval_lib::detail::div_positive(x, y.upper());
0262       else
0263         return interval<T, Policies>::empty();
0264   else
0265     return interval_lib::detail::div_non_zero(x, y);
0266 }
0267
0268 template<class T, class Policies> inline
0269 interval<T, Policies> operator/(const T& x, const interval<T, Policies>& y)
0270 {
0271   if (interval_lib::detail::test_input(x, y))
0272     return interval<T, Policies>::empty();
0273   if (zero_in(y))
0274     if (!interval_lib::user::is_zero(y.lower()))
0275       if (!interval_lib::user::is_zero(y.upper()))
0276         return interval_lib::detail::div_zero<T, Policies>(x);
0277       else
0278         return interval_lib::detail::div_negative<T, Policies>(x, y.lower());
0279     else
0280       if (!interval_lib::user::is_zero(y.upper()))
0281         return interval_lib::detail::div_positive<T, Policies>(x, y.upper());
0282       else
0283         return interval<T, Policies>::empty();
0284   else
0285     return interval_lib::detail::div_non_zero(x, y);
0286 }
0287
0288 template<class T, class Policies> inline
0289 interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y)
0290 {
0291   if (interval_lib::detail::test_input(x, y) || interval_lib::user::is_zero(y))
0292     return interval<T, Policies>::empty();
0293   typename Policies::rounding rnd;
0294   const T& xl = x.lower();
0295   const T& xu = x.upper();
0296   if (interval_lib::user::is_neg(y))
0297     return interval<T, Policies>(rnd.div_down(xu, y), rnd.div_up(xl, y), true);
0298   else
0299     return interval<T, Policies>(rnd.div_down(xl, y), rnd.div_up(xu, y), true);
0300 }
0301
0302 } // namespace numeric
0303 } // namespace boost
0304
0305 #endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP