numeric/interval/arith2.hpp

0001 /* Boost interval/arith2.hpp template implementation file
0002  *
0003  * This header provides some auxiliary arithmetic
0004  * functions: fmod, sqrt, square, pov, inverse and
0005  * a multi-interval division.
0006  *
0007  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
0008  *
0009  * Distributed under the Boost Software License, Version 1.0.
0010  * (See accompanying file LICENSE_1_0.txt or
0011  * copy at http://www.boost.org/LICENSE_1_0.txt)
0012  */
0013
0014 #ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
0015 #define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
0016
0017 #include <boost/config.hpp>
0018 #include <boost/numeric/interval/detail/interval_prototype.hpp>
0019 #include <boost/numeric/interval/detail/test_input.hpp>
0020 #include <boost/numeric/interval/detail/bugs.hpp>
0021 #include <boost/numeric/interval/detail/division.hpp>
0022 #include <boost/numeric/interval/arith.hpp>
0023 #include <boost/numeric/interval/policies.hpp>
0024 #include <algorithm>
0025 #include <cassert>
0026 #include <boost/config/no_tr1/cmath.hpp>
0027
0028 namespace boost {
0029 namespace numeric {
0030
0031 template<class T, class Policies> inline
0032 interval<T, Policies> fmod(const interval<T, Policies>& x,
0033                            const interval<T, Policies>& y)
0034 {
0035   if (interval_lib::detail::test_input(x, y))
0036     return interval<T, Policies>::empty();
0037   typename Policies::rounding rnd;
0038   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
0039   T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
0040   T n = rnd.int_down(rnd.div_down(x.lower(), yb));
0041   return (const I&)x - n * (const I&)y;
0042 }
0043
0044 template<class T, class Policies> inline
0045 interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
0046 {
0047   if (interval_lib::detail::test_input(x, y))
0048     return interval<T, Policies>::empty();
0049   typename Policies::rounding rnd;
0050   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
0051   T n = rnd.int_down(rnd.div_down(x.lower(), y));
0052   return (const I&)x - n * I(y);
0053 }
0054
0055 template<class T, class Policies> inline
0056 interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
0057 {
0058   if (interval_lib::detail::test_input(x, y))
0059     return interval<T, Policies>::empty();
0060   typename Policies::rounding rnd;
0061   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
0062   T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
0063   T n = rnd.int_down(rnd.div_down(x, yb));
0064   return x - n * (const I&)y;
0065 }
0066
0067 namespace interval_lib {
0068
0069 template<class T, class Policies> inline
0070 interval<T, Policies> division_part1(const interval<T, Policies>& x,
0071                                      const interval<T, Policies>& y, bool& b)
0072 {
0073   typedef interval<T, Policies> I;
0074   b = false;
0075   if (detail::test_input(x, y))
0076     return I::empty();
0077   if (zero_in(y))
0078     if (!user::is_zero(y.lower()))
0079       if (!user::is_zero(y.upper()))
0080         return detail::div_zero_part1(x, y, b);
0081       else
0082         return detail::div_negative(x, y.lower());
0083     else
0084       if (!user::is_zero(y.upper()))
0085         return detail::div_positive(x, y.upper());
0086       else
0087         return I::empty();
0088   else
0089     return detail::div_non_zero(x, y);
0090 }
0091
0092 template<class T, class Policies> inline
0093 interval<T, Policies> division_part2(const interval<T, Policies>& x,
0094                                      const interval<T, Policies>& y, bool b = true)
0095 {
0096   if (!b) return interval<T, Policies>::empty();
0097   return detail::div_zero_part2(x, y);
0098 }
0099
0100 template<class T, class Policies> inline
0101 interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
0102 {
0103   typedef interval<T, Policies> I;
0104   if (detail::test_input(x))
0105     return I::empty();
0106   T one = static_cast<T>(1);
0107   typename Policies::rounding rnd;
0108   if (zero_in(x)) {
0109     typedef typename Policies::checking checking;
0110     if (!user::is_zero(x.lower()))
0111       if (!user::is_zero(x.upper()))
0112         return I::whole();
0113       else
0114         return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
0115     else
0116       if (!user::is_zero(x.upper()))
0117         return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
0118       else
0119         return I::empty();
0120   } else
0121     return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
0122 }
0123
0124 namespace detail {
0125
0126 template<class T, class Rounding> inline
0127 T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
0128 {
0129   T x = x_;
0130   T y = (pwr & 1) ? x_ : static_cast<T>(1);
0131   pwr >>= 1;
0132   while (pwr > 0) {
0133     x = rnd.mul_down(x, x);
0134     if (pwr & 1) y = rnd.mul_down(x, y);
0135     pwr >>= 1;
0136   }
0137   return y;
0138 }
0139
0140 template<class T, class Rounding> inline
0141 T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
0142 {
0143   T x = x_;
0144   T y = (pwr & 1) ? x_ : static_cast<T>(1);
0145   pwr >>= 1;
0146   while (pwr > 0) {
0147     x = rnd.mul_up(x, x);
0148     if (pwr & 1) y = rnd.mul_up(x, y);
0149     pwr >>= 1;
0150   }
0151   return y;
0152 }
0153
0154 } // namespace detail
0155 } // namespace interval_lib
0156
0157 template<class T, class Policies> inline
0158 interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
0159 {
0160   BOOST_USING_STD_MAX();
0161   using interval_lib::detail::pow_dn;
0162   using interval_lib::detail::pow_up;
0163   typedef interval<T, Policies> I;
0164
0165   if (interval_lib::detail::test_input(x))
0166     return I::empty();
0167
0168   if (pwr == 0)
0169     if (interval_lib::user::is_zero(x.lower())
0170         && interval_lib::user::is_zero(x.upper()))
0171       return I::empty();
0172     else
0173       return I(static_cast<T>(1));
0174   else if (pwr < 0)
0175     return interval_lib::multiplicative_inverse(pow(x, -pwr));
0176
0177   typename Policies::rounding rnd;
0178
0179   if (interval_lib::user::is_neg(x.upper())) {        // [-2,-1]
0180     T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);
0181     T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);
0182     if (pwr & 1)     // [-2,-1]^1
0183       return I(-yu, -yl, true);
0184     else             // [-2,-1]^2
0185       return I(yl, yu, true);
0186   } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
0187     if (pwr & 1) {   // [-1,1]^1
0188       return I(-pow_up(static_cast<T>(-x.lower()), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
0189     } else {         // [-1,1]^2
0190       return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);
0191     }
0192   } else {                                // [1,2]
0193     return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
0194   }
0195 }
0196
0197 template<class T, class Policies> inline
0198 interval<T, Policies> sqrt(const interval<T, Policies>& x)
0199 {
0200   typedef interval<T, Policies> I;
0201   if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
0202     return I::empty();
0203   typename Policies::rounding rnd;
0204   T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
0205   return I(l, rnd.sqrt_up(x.upper()), true);
0206 }
0207
0208 template<class T, class Policies> inline
0209 interval<T, Policies> square(const interval<T, Policies>& x)
0210 {
0211   typedef interval<T, Policies> I;
0212   if (interval_lib::detail::test_input(x))
0213     return I::empty();
0214   typename Policies::rounding rnd;
0215   const T& xl = x.lower();
0216   const T& xu = x.upper();
0217   if (interval_lib::user::is_neg(xu))
0218     return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
0219   else if (interval_lib::user::is_pos(x.lower()))
0220     return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
0221   else
0222     return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
0223 }
0224
0225 namespace interval_lib {
0226 namespace detail {
0227
0228 template< class I > inline
0229 I root_aux(typename I::base_type const &x, int k) // x and k are bigger than one
0230 {
0231   typedef typename I::base_type T;
0232   T tk(k);
0233   I y(static_cast<T>(1), x, true);
0234   for(;;) {
0235     T y0 = median(y);
0236     I yy = intersect(y, y0 - (pow(I(y0, y0, true), k) - x) / (tk * pow(y, k - 1)));
0237     if (equal(y, yy)) return y;
0238     y = yy;
0239   }
0240 }
0241
0242 template< class I > inline // x is positive and k bigger than one
0243 typename I::base_type root_aux_dn(typename I::base_type const &x, int k)
0244 {
0245   typedef typename I::base_type T;
0246   typedef typename I::traits_type Policies;
0247   typename Policies::rounding rnd;
0248   T one(1);
0249   if (x > one) return root_aux<I>(x, k).lower();
0250   if (x == one) return one;
0251   return rnd.div_down(one, root_aux<I>(rnd.div_up(one, x), k).upper());
0252 }
0253
0254 template< class I > inline // x is positive and k bigger than one
0255 typename I::base_type root_aux_up(typename I::base_type const &x, int k)
0256 {
0257   typedef typename I::base_type T;
0258   typedef typename I::traits_type Policies;
0259   typename Policies::rounding rnd;
0260   T one(1);
0261   if (x > one) return root_aux<I>(x, k).upper();
0262   if (x == one) return one;
0263   return rnd.div_up(one, root_aux<I>(rnd.div_down(one, x), k).lower());
0264 }
0265
0266 } // namespace detail
0267 } // namespace interval_lib
0268
0269 template< class T, class Policies > inline
0270 interval<T, Policies> nth_root(interval<T, Policies> const &x, int k)
0271 {
0272   typedef interval<T, Policies> I;
0273   if (interval_lib::detail::test_input(x)) return I::empty();
0274   assert(k > 0);
0275   if (k == 1) return x;
0276   typedef typename interval_lib::unprotect<I>::type R;
0277   if (!interval_lib::user::is_pos(x.upper())) {
0278     if (interval_lib::user::is_zero(x.upper())) {
0279       T zero(0);
0280       if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,0]^/2 or [0,0]
0281         return I(zero, zero, true);
0282       else               // [-1,0]^/3
0283         return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), zero, true);
0284     } else if (!(k & 1)) // [-2,-1]^/2
0285       return I::empty();
0286     else {               // [-2,-1]^/3
0287       return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k),
0288                -interval_lib::detail::root_aux_dn<R>(-x.upper(), k), true);
0289     }
0290   }
0291   T u = interval_lib::detail::root_aux_up<R>(x.upper(), k);
0292   if (!interval_lib::user::is_pos(x.lower()))
0293     if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,1]^/2 or [0,1]
0294       return I(static_cast<T>(0), u, true);
0295     else                 // [-1,1]^/3
0296       return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), u, true);
0297   else                   // [1,2]
0298     return I(interval_lib::detail::root_aux_dn<R>(x.lower(), k), u, true);
0299 }
0300
0301 } // namespace numeric
0302 } // namespace boost
0303
0304 #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP