boost/numeric/interval/arith.hpp

   1 /* Boost interval/arith.hpp template implementation file
   2  *
   3  * Copyright 2000 Jens Maurer
   4  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
   5  *
   6  * Distributed under the Boost Software License, Version 1.0.
   7  * (See accompanying file LICENSE_1_0.txt or
   8  * copy at http://www.boost.org/LICENSE_1_0.txt)
   9  */
  10
  11 #ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP
  12 #define BOOST_NUMERIC_INTERVAL_ARITH_HPP
  13
  14 #include <boost/config.hpp>
  15 #include <boost/numeric/interval/interval.hpp>
  16 #include <boost/numeric/interval/detail/bugs.hpp>
  17 #include <boost/numeric/interval/detail/test_input.hpp>
  18 #include <boost/numeric/interval/detail/division.hpp>
  19 #include <algorithm>
  20
  21 namespace boost {
  22 namespace numeric {
  23
  24 /*
  25  * Basic arithmetic operators
  26  */
  27
  28 template<class T, class Policies> inline
  29 const interval<T, Policies>& operator+(const interval<T, Policies>& x)
  30 {
  31   return x;
  32 }
  33
  34 template<class T, class Policies> inline
  35 interval<T, Policies> operator-(const interval<T, Policies>& x)
  36 {
  37   if (interval_lib::detail::test_input(x))
  38     return interval<T, Policies>::empty();
  39   return interval<T, Policies>(-x.upper(), -x.lower(), true);
  40 }
  41
  42 template<class T, class Policies> inline
  43 interval<T, Policies>& interval<T, Policies>::operator+=(const interval<T, Policies>& r)
  44 {
  45   if (interval_lib::detail::test_input(*this, r))
  46     set_empty();
  47   else {
  48     typename Policies::rounding rnd;
  49     set(rnd.add_down(low, r.low), rnd.add_up(up, r.up));
  50   }
  51   return *this;
  52 }
  53
  54 template<class T, class Policies> inline
  55 interval<T, Policies>& interval<T, Policies>::operator+=(const T& r)
  56 {
  57   if (interval_lib::detail::test_input(*this, r))
  58     set_empty();
  59   else {
  60     typename Policies::rounding rnd;
  61     set(rnd.add_down(low, r), rnd.add_up(up, r));
  62   }
  63   return *this;
  64 }
  65
  66 template<class T, class Policies> inline
  67 interval<T, Policies>& interval<T, Policies>::operator-=(const interval<T, Policies>& r)
  68 {
  69   if (interval_lib::detail::test_input(*this, r))
  70     set_empty();
  71   else {
  72     typename Policies::rounding rnd;
  73     set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low));
  74   }
  75   return *this;
  76 }
  77
  78 template<class T, class Policies> inline
  79 interval<T, Policies>& interval<T, Policies>::operator-=(const T& r)
  80 {
  81   if (interval_lib::detail::test_input(*this, r))
  82     set_empty();
  83   else {
  84     typename Policies::rounding rnd;
  85     set(rnd.sub_down(low, r), rnd.sub_up(up, r));
  86   }
  87   return *this;
  88 }
  89
  90 template<class T, class Policies> inline
  91 interval<T, Policies>& interval<T, Policies>::operator*=(const interval<T, Policies>& r)
  92 {
  93   return *this = *this * r;
  94 }
  95
  96 template<class T, class Policies> inline
  97 interval<T, Policies>& interval<T, Policies>::operator*=(const T& r)
  98 {
  99   return *this = r * *this;
 100 }
 101
 102 template<class T, class Policies> inline
 103 interval<T, Policies>& interval<T, Policies>::operator/=(const interval<T, Policies>& r)
 104 {
 105   return *this = *this / r;
 106 }
 107
 108 template<class T, class Policies> inline
 109 interval<T, Policies>& interval<T, Policies>::operator/=(const T& r)
 110 {
 111   return *this = *this / r;
 112 }
 113
 114 template<class T, class Policies> inline
 115 interval<T, Policies> operator+(const interval<T, Policies>& x,
 116                                 const interval<T, Policies>& y)
 117 {
 118   if (interval_lib::detail::test_input(x, y))
 119     return interval<T, Policies>::empty();
 120   typename Policies::rounding rnd;
 121   return interval<T,Policies>(rnd.add_down(x.lower(), y.lower()),
 122                               rnd.add_up  (x.upper(), y.upper()), true);
 123 }
 124
 125 template<class T, class Policies> inline
 126 interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y)
 127 {
 128   if (interval_lib::detail::test_input(x, y))
 129     return interval<T, Policies>::empty();
 130   typename Policies::rounding rnd;
 131   return interval<T,Policies>(rnd.add_down(x, y.lower()),
 132                               rnd.add_up  (x, y.upper()), true);
 133 }
 134
 135 template<class T, class Policies> inline
 136 interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)
 137 { return y + x; }
 138
 139 template<class T, class Policies> inline
 140 interval<T, Policies> operator-(const interval<T, Policies>& x,
 141                                 const interval<T, Policies>& y)
 142 {
 143   if (interval_lib::detail::test_input(x, y))
 144     return interval<T, Policies>::empty();
 145   typename Policies::rounding rnd;
 146   return interval<T,Policies>(rnd.sub_down(x.lower(), y.upper()),
 147                               rnd.sub_up  (x.upper(), y.lower()), true);
 148 }
 149
 150 template<class T, class Policies> inline
 151 interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y)
 152 {
 153   if (interval_lib::detail::test_input(x, y))
 154     return interval<T, Policies>::empty();
 155   typename Policies::rounding rnd;
 156   return interval<T,Policies>(rnd.sub_down(x, y.upper()),
 157                               rnd.sub_up  (x, y.lower()), true);
 158 }
 159
 160 template<class T, class Policies> inline
 161 interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y)
 162 {
 163   if (interval_lib::detail::test_input(x, y))
 164     return interval<T, Policies>::empty();
 165   typename Policies::rounding rnd;
 166   return interval<T,Policies>(rnd.sub_down(x.lower(), y),
 167                               rnd.sub_up  (x.upper(), y), true);
 168 }
 169
 170 template<class T, class Policies> inline
 171 interval<T, Policies> operator*(const interval<T, Policies>& x,
 172                                 const interval<T, Policies>& y)
 173 {
 174   BOOST_USING_STD_MIN();
 175   BOOST_USING_STD_MAX();
 176   typedef interval<T, Policies> I;
 177   if (interval_lib::detail::test_input(x, y))
 178     return I::empty();
 179   typename Policies::rounding rnd;
 180   const T& xl = x.lower();
 181   const T& xu = x.upper();
 182   const T& yl = y.lower();
 183   const T& yu = y.upper();
 184
 185   if (interval_lib::user::is_neg(xl))
 186     if (interval_lib::user::is_pos(xu))
 187       if (interval_lib::user::is_neg(yl))
 188         if (interval_lib::user::is_pos(yu)) // M * M
 189           return I(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)),
 190                    max BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_up  (xl, yl), rnd.mul_up  (xu, yu)), true);
 191         else                    // M * N
 192           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
 193       else
 194         if (interval_lib::user::is_pos(yu)) // M * P
 195           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
 196         else                    // M * Z
 197           return I(static_cast<T>(0), static_cast<T>(0), true);
 198     else
 199       if (interval_lib::user::is_neg(yl))
 200         if (interval_lib::user::is_pos(yu)) // N * M
 201           return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
 202         else                    // N * N
 203           return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
 204       else
 205         if (interval_lib::user::is_pos(yu)) // N * P
 206           return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
 207         else                    // N * Z
 208           return I(static_cast<T>(0), static_cast<T>(0), true);
 209   else
 210     if (interval_lib::user::is_pos(xu))
 211       if (interval_lib::user::is_neg(yl))
 212         if (interval_lib::user::is_pos(yu)) // P * M
 213           return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
 214         else                    // P * N
 215           return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
 216       else
 217         if (interval_lib::user::is_pos(yu)) // P * P
 218           return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
 219         else                    // P * Z
 220           return I(static_cast<T>(0), static_cast<T>(0), true);
 221     else                        // Z * ?
 222       return I(static_cast<T>(0), static_cast<T>(0), true);
 223 }
 224
 225 template<class T, class Policies> inline
 226 interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y)
 227 {
 228   typedef interval<T, Policies> I;
 229   if (interval_lib::detail::test_input(x, y))
 230     return I::empty();
 231   typename Policies::rounding rnd;
 232   const T& yl = y.lower();
 233   const T& yu = y.upper();
 234   // x is supposed not to be infinite
 235   if (interval_lib::user::is_neg(x))
 236     return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true);
 237   else if (interval_lib::user::is_zero(x))
 238     return I(static_cast<T>(0), static_cast<T>(0), true);
 239   else
 240     return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true);
 241 }
 242
 243 template<class T, class Policies> inline
 244 interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)
 245 { return y * x; }
 246
 247 template<class T, class Policies> inline
 248 interval<T, Policies> operator/(const interval<T, Policies>& x,
 249                                 const interval<T, Policies>& y)
 250 {
 251   if (interval_lib::detail::test_input(x, y))
 252     return interval<T, Policies>::empty();
 253   if (zero_in(y))
 254     if (!interval_lib::user::is_zero(y.lower()))
 255       if (!interval_lib::user::is_zero(y.upper()))
 256         return interval_lib::detail::div_zero(x);
 257       else
 258         return interval_lib::detail::div_negative(x, y.lower());
 259     else
 260       if (!interval_lib::user::is_zero(y.upper()))
 261         return interval_lib::detail::div_positive(x, y.upper());
 262       else
 263         return interval<T, Policies>::empty();
 264   else
 265     return interval_lib::detail::div_non_zero(x, y);
 266 }
 267
 268 template<class T, class Policies> inline
 269 interval<T, Policies> operator/(const T& x, const interval<T, Policies>& y)
 270 {
 271   if (interval_lib::detail::test_input(x, y))
 272     return interval<T, Policies>::empty();
 273   if (zero_in(y))
 274     if (!interval_lib::user::is_zero(y.lower()))
 275       if (!interval_lib::user::is_zero(y.upper()))
 276         return interval_lib::detail::div_zero<T, Policies>(x);
 277       else
 278         return interval_lib::detail::div_negative<T, Policies>(x, y.lower());
 279     else
 280       if (!interval_lib::user::is_zero(y.upper()))
 281         return interval_lib::detail::div_positive<T, Policies>(x, y.upper());
 282       else
 283         return interval<T, Policies>::empty();
 284   else
 285     return interval_lib::detail::div_non_zero(x, y);
 286 }
 287
 288 template<class T, class Policies> inline
 289 interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y)
 290 {
 291   if (interval_lib::detail::test_input(x, y) || interval_lib::user::is_zero(y))
 292     return interval<T, Policies>::empty();
 293   typename Policies::rounding rnd;
 294   const T& xl = x.lower();
 295   const T& xu = x.upper();
 296   if (interval_lib::user::is_neg(y))
 297     return interval<T, Policies>(rnd.div_down(xu, y), rnd.div_up(xl, y), true);
 298   else
 299     return interval<T, Policies>(rnd.div_down(xl, y), rnd.div_up(xu, y), true);
 300 }
 301
 302 } // namespace numeric
 303 } // namespace boost
 304
 305 #endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP