MercurialSymaptic / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | raw | help
os/ossrv/ossrv_pub/boost_apis/boost/numeric/interval/arith2.hpp
author sl
Tue, 10 Jun 2014 14:32:02 +0200
changeset 1 260cb5ec6c19
permissions -rw-r--r--
Update contrib. 
     1 /* Boost interval/arith2.hpp template implementation file

     2  *

     3  * This header provides some auxiliary arithmetic

     4  * functions: fmod, sqrt, square, pov, inverse and

     5  * a multi-interval division.

     6  *

     7  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion

     8  *

     9  * Distributed under the Boost Software License, Version 1.0.

    10  * (See accompanying file LICENSE_1_0.txt or

    11  * copy at http://www.boost.org/LICENSE_1_0.txt)

    12  */

    13 

    14 #ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP

    15 #define BOOST_NUMERIC_INTERVAL_ARITH2_HPP

    16 

    17 #include <boost/config.hpp>

    18 #include <boost/numeric/interval/detail/interval_prototype.hpp>

    19 #include <boost/numeric/interval/detail/test_input.hpp>

    20 #include <boost/numeric/interval/detail/bugs.hpp>

    21 #include <boost/numeric/interval/detail/division.hpp>

    22 #include <boost/numeric/interval/arith.hpp>

    23 #include <boost/numeric/interval/policies.hpp>

    24 #include <algorithm>

    25 #include <cmath>

    26 

    27 namespace boost {

    28 namespace numeric {

    29 

    30 template<class T, class Policies> inline

    31 interval<T, Policies> fmod(const interval<T, Policies>& x,

    32                            const interval<T, Policies>& y)

    33 {

    34   if (interval_lib::detail::test_input(x, y))

    35     return interval<T, Policies>::empty();

    36   typename Policies::rounding rnd;

    37   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;

    38   T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();

    39   T n = rnd.int_down(rnd.div_down(x.lower(), yb));

    40   return (const I&)x - n * (const I&)y;

    41 }

    42 

    43 template<class T, class Policies> inline

    44 interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)

    45 {

    46   if (interval_lib::detail::test_input(x, y))

    47     return interval<T, Policies>::empty();

    48   typename Policies::rounding rnd;

    49   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;

    50   T n = rnd.int_down(rnd.div_down(x.lower(), y));

    51   return (const I&)x - n * I(y);

    52 }

    53 

    54 template<class T, class Policies> inline

    55 interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)

    56 {

    57   if (interval_lib::detail::test_input(x, y))

    58     return interval<T, Policies>::empty();

    59   typename Policies::rounding rnd;

    60   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;

    61   T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();

    62   T n = rnd.int_down(rnd.div_down(x, yb));

    63   return x - n * (const I&)y;

    64 }

    65 

    66 namespace interval_lib {

    67 

    68 template<class T, class Policies> inline

    69 interval<T, Policies> division_part1(const interval<T, Policies>& x,

    70                                      const interval<T, Policies>& y, bool& b)

    71 {

    72   typedef interval<T, Policies> I;

    73   b = false;

    74   if (detail::test_input(x, y))

    75     return I::empty();

    76   if (zero_in(y))

    77     if (!user::is_zero(y.lower()))

    78       if (!user::is_zero(y.upper()))

    79         return detail::div_zero_part1(x, y, b);

    80       else

    81         return detail::div_negative(x, y.lower());

    82     else

    83       if (!user::is_zero(y.upper()))

    84         return detail::div_positive(x, y.upper());

    85       else

    86         return I::empty();

    87   else

    88     return detail::div_non_zero(x, y);

    89 }

    90 

    91 template<class T, class Policies> inline

    92 interval<T, Policies> division_part2(const interval<T, Policies>& x,

    93                                      const interval<T, Policies>& y, bool b = true)

    94 {

    95   if (!b) return interval<T, Policies>::empty();

    96   return detail::div_zero_part2(x, y);

    97 }

    98 

    99 template<class T, class Policies> inline

   100 interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)

   101 {

   102   typedef interval<T, Policies> I;

   103   if (detail::test_input(x))

   104     return I::empty();

   105   T one = static_cast<T>(1);

   106   typename Policies::rounding rnd;

   107   if (zero_in(x)) {

   108     typedef typename Policies::checking checking;

   109     if (!user::is_zero(x.lower()))

   110       if (!user::is_zero(x.upper()))

   111         return I::whole();

   112       else

   113         return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);

   114     else

   115       if (!user::is_zero(x.upper()))

   116         return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);

   117       else

   118         return I::empty();

   119   } else 

   120     return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);

   121 }

   122 

   123 namespace detail {

   124 

   125 template<class T, class Rounding> inline

   126 T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive

   127 {

   128   T x = x_;

   129   T y = (pwr & 1) ? x_ : static_cast<T>(1);

   130   pwr >>= 1;

   131   while (pwr > 0) {

   132     x = rnd.mul_down(x, x);

   133     if (pwr & 1) y = rnd.mul_down(x, y);

   134     pwr >>= 1;

   135   }

   136   return y;

   137 }

   138 

   139 template<class T, class Rounding> inline

   140 T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive

   141 {

   142   T x = x_;

   143   T y = (pwr & 1) ? x_ : static_cast<T>(1);

   144   pwr >>= 1;

   145   while (pwr > 0) {

   146     x = rnd.mul_up(x, x);

   147     if (pwr & 1) y = rnd.mul_up(x, y);

   148     pwr >>= 1;

   149   }

   150   return y;

   151 }

   152 

   153 } // namespace detail

   154 } // namespace interval_lib

   155 

   156 template<class T, class Policies> inline

   157 interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)

   158 {

   159   BOOST_USING_STD_MAX();

   160   using interval_lib::detail::pow_dn;

   161   using interval_lib::detail::pow_up;

   162   typedef interval<T, Policies> I;

   163 

   164   if (interval_lib::detail::test_input(x))

   165     return I::empty();

   166 

   167   if (pwr == 0)

   168     if (interval_lib::user::is_zero(x.lower())

   169         && interval_lib::user::is_zero(x.upper()))

   170       return I::empty();

   171     else

   172       return I(static_cast<T>(1));

   173   else if (pwr < 0)

   174     return interval_lib::multiplicative_inverse(pow(x, -pwr));

   175 

   176   typename Policies::rounding rnd;

   177   

   178   if (interval_lib::user::is_neg(x.upper())) {        // [-2,-1]

   179     T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);

   180     T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);

   181     if (pwr & 1)     // [-2,-1]^1

   182       return I(-yu, -yl, true);

   183     else             // [-2,-1]^2

   184       return I(yl, yu, true);

   185   } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]

   186     if (pwr & 1) {   // [-1,1]^1

   187       return I(-pow_up(-x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);

   188     } else {         // [-1,1]^2

   189       return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);

   190     }

   191   } else {                                // [1,2]

   192     return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);

   193   }

   194 }

   195 

   196 template<class T, class Policies> inline

   197 interval<T, Policies> sqrt(const interval<T, Policies>& x)

   198 {

   199   typedef interval<T, Policies> I;

   200   if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))

   201     return I::empty();

   202   typename Policies::rounding rnd;

   203   T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());

   204   return I(l, rnd.sqrt_up(x.upper()), true);

   205 }

   206 

   207 template<class T, class Policies> inline

   208 interval<T, Policies> square(const interval<T, Policies>& x)

   209 {

   210   typedef interval<T, Policies> I;

   211   if (interval_lib::detail::test_input(x))

   212     return I::empty();

   213   typename Policies::rounding rnd;

   214   const T& xl = x.lower();

   215   const T& xu = x.upper();

   216   if (interval_lib::user::is_neg(xu))

   217     return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);

   218   else if (interval_lib::user::is_pos(x.lower()))

   219     return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);

   220   else

   221     return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);

   222 }

   223 

   224 } // namespace numeric

   225 } // namespace boost

   226 

   227 #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP
Symaptic