boost/numeric/interval/arith2.hpp

   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 <cassert>
  26 #include <boost/config/no_tr1/cmath.hpp>
  27
  28 namespace boost {
  29 namespace numeric {
  30
  31 template<class T, class Policies> inline
  32 interval<T, Policies> fmod(const interval<T, Policies>& x,
  33                            const interval<T, Policies>& y)
  34 {
  35   if (interval_lib::detail::test_input(x, y))
  36     return interval<T, Policies>::empty();
  37   typename Policies::rounding rnd;
  38   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
  39   T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
  40   T n = rnd.int_down(rnd.div_down(x.lower(), yb));
  41   return (const I&)x - n * (const I&)y;
  42 }
  43
  44 template<class T, class Policies> inline
  45 interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
  46 {
  47   if (interval_lib::detail::test_input(x, y))
  48     return interval<T, Policies>::empty();
  49   typename Policies::rounding rnd;
  50   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
  51   T n = rnd.int_down(rnd.div_down(x.lower(), y));
  52   return (const I&)x - n * I(y);
  53 }
  54
  55 template<class T, class Policies> inline
  56 interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
  57 {
  58   if (interval_lib::detail::test_input(x, y))
  59     return interval<T, Policies>::empty();
  60   typename Policies::rounding rnd;
  61   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
  62   T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
  63   T n = rnd.int_down(rnd.div_down(x, yb));
  64   return x - n * (const I&)y;
  65 }
  66
  67 namespace interval_lib {
  68
  69 template<class T, class Policies> inline
  70 interval<T, Policies> division_part1(const interval<T, Policies>& x,
  71                                      const interval<T, Policies>& y, bool& b)
  72 {
  73   typedef interval<T, Policies> I;
  74   b = false;
  75   if (detail::test_input(x, y))
  76     return I::empty();
  77   if (zero_in(y))
  78     if (!user::is_zero(y.lower()))
  79       if (!user::is_zero(y.upper()))
  80         return detail::div_zero_part1(x, y, b);
  81       else
  82         return detail::div_negative(x, y.lower());
  83     else
  84       if (!user::is_zero(y.upper()))
  85         return detail::div_positive(x, y.upper());
  86       else
  87         return I::empty();
  88   else
  89     return detail::div_non_zero(x, y);
  90 }
  91
  92 template<class T, class Policies> inline
  93 interval<T, Policies> division_part2(const interval<T, Policies>& x,
  94                                      const interval<T, Policies>& y, bool b = true)
  95 {
  96   if (!b) return interval<T, Policies>::empty();
  97   return detail::div_zero_part2(x, y);
  98 }
  99
 100 template<class T, class Policies> inline
 101 interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
 102 {
 103   typedef interval<T, Policies> I;
 104   if (detail::test_input(x))
 105     return I::empty();
 106   T one = static_cast<T>(1);
 107   typename Policies::rounding rnd;
 108   if (zero_in(x)) {
 109     typedef typename Policies::checking checking;
 110     if (!user::is_zero(x.lower()))
 111       if (!user::is_zero(x.upper()))
 112         return I::whole();
 113       else
 114         return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
 115     else
 116       if (!user::is_zero(x.upper()))
 117         return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
 118       else
 119         return I::empty();
 120   } else
 121     return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
 122 }
 123
 124 namespace detail {
 125
 126 template<class T, class Rounding> inline
 127 T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
 128 {
 129   T x = x_;
 130   T y = (pwr & 1) ? x_ : static_cast<T>(1);
 131   pwr >>= 1;
 132   while (pwr > 0) {
 133     x = rnd.mul_down(x, x);
 134     if (pwr & 1) y = rnd.mul_down(x, y);
 135     pwr >>= 1;
 136   }
 137   return y;
 138 }
 139
 140 template<class T, class Rounding> inline
 141 T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
 142 {
 143   T x = x_;
 144   T y = (pwr & 1) ? x_ : static_cast<T>(1);
 145   pwr >>= 1;
 146   while (pwr > 0) {
 147     x = rnd.mul_up(x, x);
 148     if (pwr & 1) y = rnd.mul_up(x, y);
 149     pwr >>= 1;
 150   }
 151   return y;
 152 }
 153
 154 } // namespace detail
 155 } // namespace interval_lib
 156
 157 template<class T, class Policies> inline
 158 interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
 159 {
 160   BOOST_USING_STD_MAX();
 161   using interval_lib::detail::pow_dn;
 162   using interval_lib::detail::pow_up;
 163   typedef interval<T, Policies> I;
 164
 165   if (interval_lib::detail::test_input(x))
 166     return I::empty();
 167
 168   if (pwr == 0)
 169     if (interval_lib::user::is_zero(x.lower())
 170         && interval_lib::user::is_zero(x.upper()))
 171       return I::empty();
 172     else
 173       return I(static_cast<T>(1));
 174   else if (pwr < 0)
 175     return interval_lib::multiplicative_inverse(pow(x, -pwr));
 176
 177   typename Policies::rounding rnd;
 178
 179   if (interval_lib::user::is_neg(x.upper())) {        // [-2,-1]
 180     T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);
 181     T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);
 182     if (pwr & 1)     // [-2,-1]^1
 183       return I(-yu, -yl, true);
 184     else             // [-2,-1]^2
 185       return I(yl, yu, true);
 186   } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
 187     if (pwr & 1) {   // [-1,1]^1
 188       return I(-pow_up(static_cast<T>(-x.lower()), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
 189     } else {         // [-1,1]^2
 190       return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);
 191     }
 192   } else {                                // [1,2]
 193     return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
 194   }
 195 }
 196
 197 template<class T, class Policies> inline
 198 interval<T, Policies> sqrt(const interval<T, Policies>& x)
 199 {
 200   typedef interval<T, Policies> I;
 201   if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
 202     return I::empty();
 203   typename Policies::rounding rnd;
 204   T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
 205   return I(l, rnd.sqrt_up(x.upper()), true);
 206 }
 207
 208 template<class T, class Policies> inline
 209 interval<T, Policies> square(const interval<T, Policies>& x)
 210 {
 211   typedef interval<T, Policies> I;
 212   if (interval_lib::detail::test_input(x))
 213     return I::empty();
 214   typename Policies::rounding rnd;
 215   const T& xl = x.lower();
 216   const T& xu = x.upper();
 217   if (interval_lib::user::is_neg(xu))
 218     return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
 219   else if (interval_lib::user::is_pos(x.lower()))
 220     return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
 221   else
 222     return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
 223 }
 224
 225 namespace interval_lib {
 226 namespace detail {
 227
 228 template< class I > inline
 229 I root_aux(typename I::base_type const &x, int k) // x and k are bigger than one
 230 {
 231   typedef typename I::base_type T;
 232   T tk(k);
 233   I y(static_cast<T>(1), x, true);
 234   for(;;) {
 235     T y0 = median(y);
 236     I yy = intersect(y, y0 - (pow(I(y0, y0, true), k) - x) / (tk * pow(y, k - 1)));
 237     if (equal(y, yy)) return y;
 238     y = yy;
 239   }
 240 }
 241
 242 template< class I > inline // x is positive and k bigger than one
 243 typename I::base_type root_aux_dn(typename I::base_type const &x, int k)
 244 {
 245   typedef typename I::base_type T;
 246   typedef typename I::traits_type Policies;
 247   typename Policies::rounding rnd;
 248   T one(1);
 249   if (x > one) return root_aux<I>(x, k).lower();
 250   if (x == one) return one;
 251   return rnd.div_down(one, root_aux<I>(rnd.div_up(one, x), k).upper());
 252 }
 253
 254 template< class I > inline // x is positive and k bigger than one
 255 typename I::base_type root_aux_up(typename I::base_type const &x, int k)
 256 {
 257   typedef typename I::base_type T;
 258   typedef typename I::traits_type Policies;
 259   typename Policies::rounding rnd;
 260   T one(1);
 261   if (x > one) return root_aux<I>(x, k).upper();
 262   if (x == one) return one;
 263   return rnd.div_up(one, root_aux<I>(rnd.div_down(one, x), k).lower());
 264 }
 265
 266 } // namespace detail
 267 } // namespace interval_lib
 268
 269 template< class T, class Policies > inline
 270 interval<T, Policies> nth_root(interval<T, Policies> const &x, int k)
 271 {
 272   typedef interval<T, Policies> I;
 273   if (interval_lib::detail::test_input(x)) return I::empty();
 274   assert(k > 0);
 275   if (k == 1) return x;
 276   typename Policies::rounding rnd;
 277   typedef typename interval_lib::unprotect<I>::type R;
 278   if (!interval_lib::user::is_pos(x.upper())) {
 279     if (interval_lib::user::is_zero(x.upper())) {
 280       T zero(0);
 281       if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,0]^/2 or [0,0]
 282         return I(zero, zero, true);
 283       else               // [-1,0]^/3
 284         return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), zero, true);
 285     } else if (!(k & 1)) // [-2,-1]^/2
 286       return I::empty();
 287     else {               // [-2,-1]^/3
 288       return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k),
 289                -interval_lib::detail::root_aux_dn<R>(-x.upper(), k), true);
 290     }
 291   }
 292   T u = interval_lib::detail::root_aux_up<R>(x.upper(), k);
 293   if (!interval_lib::user::is_pos(x.lower()))
 294     if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,1]^/2 or [0,1]
 295       return I(static_cast<T>(0), u, true);
 296     else                 // [-1,1]^/3
 297       return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), u, true);
 298   else                   // [1,2]
 299     return I(interval_lib::detail::root_aux_dn<R>(x.lower(), k), u, true);
 300 }
 301
 302 } // namespace numeric
 303 } // namespace boost
 304
 305 #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP