sl@0: /* Boost interval/arith2.hpp template implementation file
sl@0:  *
sl@0:  * This header provides some auxiliary arithmetic
sl@0:  * functions: fmod, sqrt, square, pov, inverse and
sl@0:  * a multi-interval division.
sl@0:  *
sl@0:  * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
sl@0:  *
sl@0:  * Distributed under the Boost Software License, Version 1.0.
sl@0:  * (See accompanying file LICENSE_1_0.txt or
sl@0:  * copy at http://www.boost.org/LICENSE_1_0.txt)
sl@0:  */
sl@0: 
sl@0: #ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
sl@0: #define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
sl@0: 
sl@0: #include <boost/config.hpp>
sl@0: #include <boost/numeric/interval/detail/interval_prototype.hpp>
sl@0: #include <boost/numeric/interval/detail/test_input.hpp>
sl@0: #include <boost/numeric/interval/detail/bugs.hpp>
sl@0: #include <boost/numeric/interval/detail/division.hpp>
sl@0: #include <boost/numeric/interval/arith.hpp>
sl@0: #include <boost/numeric/interval/policies.hpp>
sl@0: #include <algorithm>
sl@0: #include <cmath>
sl@0: 
sl@0: namespace boost {
sl@0: namespace numeric {
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> fmod(const interval<T, Policies>& x,
sl@0:                            const interval<T, Policies>& y)
sl@0: {
sl@0:   if (interval_lib::detail::test_input(x, y))
sl@0:     return interval<T, Policies>::empty();
sl@0:   typename Policies::rounding rnd;
sl@0:   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
sl@0:   T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
sl@0:   T n = rnd.int_down(rnd.div_down(x.lower(), yb));
sl@0:   return (const I&)x - n * (const I&)y;
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
sl@0: {
sl@0:   if (interval_lib::detail::test_input(x, y))
sl@0:     return interval<T, Policies>::empty();
sl@0:   typename Policies::rounding rnd;
sl@0:   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
sl@0:   T n = rnd.int_down(rnd.div_down(x.lower(), y));
sl@0:   return (const I&)x - n * I(y);
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
sl@0: {
sl@0:   if (interval_lib::detail::test_input(x, y))
sl@0:     return interval<T, Policies>::empty();
sl@0:   typename Policies::rounding rnd;
sl@0:   typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
sl@0:   T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
sl@0:   T n = rnd.int_down(rnd.div_down(x, yb));
sl@0:   return x - n * (const I&)y;
sl@0: }
sl@0: 
sl@0: namespace interval_lib {
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> division_part1(const interval<T, Policies>& x,
sl@0:                                      const interval<T, Policies>& y, bool& b)
sl@0: {
sl@0:   typedef interval<T, Policies> I;
sl@0:   b = false;
sl@0:   if (detail::test_input(x, y))
sl@0:     return I::empty();
sl@0:   if (zero_in(y))
sl@0:     if (!user::is_zero(y.lower()))
sl@0:       if (!user::is_zero(y.upper()))
sl@0:         return detail::div_zero_part1(x, y, b);
sl@0:       else
sl@0:         return detail::div_negative(x, y.lower());
sl@0:     else
sl@0:       if (!user::is_zero(y.upper()))
sl@0:         return detail::div_positive(x, y.upper());
sl@0:       else
sl@0:         return I::empty();
sl@0:   else
sl@0:     return detail::div_non_zero(x, y);
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> division_part2(const interval<T, Policies>& x,
sl@0:                                      const interval<T, Policies>& y, bool b = true)
sl@0: {
sl@0:   if (!b) return interval<T, Policies>::empty();
sl@0:   return detail::div_zero_part2(x, y);
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
sl@0: {
sl@0:   typedef interval<T, Policies> I;
sl@0:   if (detail::test_input(x))
sl@0:     return I::empty();
sl@0:   T one = static_cast<T>(1);
sl@0:   typename Policies::rounding rnd;
sl@0:   if (zero_in(x)) {
sl@0:     typedef typename Policies::checking checking;
sl@0:     if (!user::is_zero(x.lower()))
sl@0:       if (!user::is_zero(x.upper()))
sl@0:         return I::whole();
sl@0:       else
sl@0:         return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
sl@0:     else
sl@0:       if (!user::is_zero(x.upper()))
sl@0:         return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
sl@0:       else
sl@0:         return I::empty();
sl@0:   } else 
sl@0:     return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
sl@0: }
sl@0: 
sl@0: namespace detail {
sl@0: 
sl@0: template<class T, class Rounding> inline
sl@0: T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
sl@0: {
sl@0:   T x = x_;
sl@0:   T y = (pwr & 1) ? x_ : static_cast<T>(1);
sl@0:   pwr >>= 1;
sl@0:   while (pwr > 0) {
sl@0:     x = rnd.mul_down(x, x);
sl@0:     if (pwr & 1) y = rnd.mul_down(x, y);
sl@0:     pwr >>= 1;
sl@0:   }
sl@0:   return y;
sl@0: }
sl@0: 
sl@0: template<class T, class Rounding> inline
sl@0: T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
sl@0: {
sl@0:   T x = x_;
sl@0:   T y = (pwr & 1) ? x_ : static_cast<T>(1);
sl@0:   pwr >>= 1;
sl@0:   while (pwr > 0) {
sl@0:     x = rnd.mul_up(x, x);
sl@0:     if (pwr & 1) y = rnd.mul_up(x, y);
sl@0:     pwr >>= 1;
sl@0:   }
sl@0:   return y;
sl@0: }
sl@0: 
sl@0: } // namespace detail
sl@0: } // namespace interval_lib
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
sl@0: {
sl@0:   BOOST_USING_STD_MAX();
sl@0:   using interval_lib::detail::pow_dn;
sl@0:   using interval_lib::detail::pow_up;
sl@0:   typedef interval<T, Policies> I;
sl@0: 
sl@0:   if (interval_lib::detail::test_input(x))
sl@0:     return I::empty();
sl@0: 
sl@0:   if (pwr == 0)
sl@0:     if (interval_lib::user::is_zero(x.lower())
sl@0:         && interval_lib::user::is_zero(x.upper()))
sl@0:       return I::empty();
sl@0:     else
sl@0:       return I(static_cast<T>(1));
sl@0:   else if (pwr < 0)
sl@0:     return interval_lib::multiplicative_inverse(pow(x, -pwr));
sl@0: 
sl@0:   typename Policies::rounding rnd;
sl@0:   
sl@0:   if (interval_lib::user::is_neg(x.upper())) {        // [-2,-1]
sl@0:     T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);
sl@0:     T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);
sl@0:     if (pwr & 1)     // [-2,-1]^1
sl@0:       return I(-yu, -yl, true);
sl@0:     else             // [-2,-1]^2
sl@0:       return I(yl, yu, true);
sl@0:   } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
sl@0:     if (pwr & 1) {   // [-1,1]^1
sl@0:       return I(-pow_up(-x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
sl@0:     } else {         // [-1,1]^2
sl@0:       return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);
sl@0:     }
sl@0:   } else {                                // [1,2]
sl@0:     return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
sl@0:   }
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> sqrt(const interval<T, Policies>& x)
sl@0: {
sl@0:   typedef interval<T, Policies> I;
sl@0:   if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
sl@0:     return I::empty();
sl@0:   typename Policies::rounding rnd;
sl@0:   T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
sl@0:   return I(l, rnd.sqrt_up(x.upper()), true);
sl@0: }
sl@0: 
sl@0: template<class T, class Policies> inline
sl@0: interval<T, Policies> square(const interval<T, Policies>& x)
sl@0: {
sl@0:   typedef interval<T, Policies> I;
sl@0:   if (interval_lib::detail::test_input(x))
sl@0:     return I::empty();
sl@0:   typename Policies::rounding rnd;
sl@0:   const T& xl = x.lower();
sl@0:   const T& xu = x.upper();
sl@0:   if (interval_lib::user::is_neg(xu))
sl@0:     return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
sl@0:   else if (interval_lib::user::is_pos(x.lower()))
sl@0:     return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
sl@0:   else
sl@0:     return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
sl@0: }
sl@0: 
sl@0: } // namespace numeric
sl@0: } // namespace boost
sl@0: 
sl@0: #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP