Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/interval/arith.hpp@bde4ae8d615e

     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  */

    11 #ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP

    12 #define BOOST_NUMERIC_INTERVAL_ARITH_HPP

    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>

    21 namespace boost {

    22 namespace numeric {

    24 /*

    25  * Basic arithmetic operators

    26  */

    28 template<class T, class Policies> inline

    29 const interval<T, Policies>& operator+(const interval<T, Policies>& x)

    30 {

    31   return x;

    32 }

    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 }

    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 }

    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 }

    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 }

    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 }

    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 }

    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 }

   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 }

   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 }

   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 }

   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 }

   135 template<class T, class Policies> inline

   136 interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)

   137 { return y + x; }

   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 }

   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 }

   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 }

   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();

   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 }

   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 }

   243 template<class T, class Policies> inline

   244 interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)

   245 { return y * x; }

   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 }

   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 }

   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 }

   302 } // namespace numeric

   303 } // namespace boost

   305 #endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP

author	sl@SLION-WIN7.fritz.box
	Fri, 15 Jun 2012 03:10:57 +0200
changeset 0	bde4ae8d615e
permissions	-rw-r--r--