Mercurial > hg > vamp-build-and-test
diff DEPENDENCIES/generic/include/boost/numeric/interval/arith2.hpp @ 16:2665513ce2d3
Add boost headers
author | Chris Cannam |
---|---|
date | Tue, 05 Aug 2014 11:11:38 +0100 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DEPENDENCIES/generic/include/boost/numeric/interval/arith2.hpp Tue Aug 05 11:11:38 2014 +0100 @@ -0,0 +1,305 @@ +/* Boost interval/arith2.hpp template implementation file + * + * This header provides some auxiliary arithmetic + * functions: fmod, sqrt, square, pov, inverse and + * a multi-interval division. + * + * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion + * + * Distributed under the Boost Software License, Version 1.0. + * (See accompanying file LICENSE_1_0.txt or + * copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP +#define BOOST_NUMERIC_INTERVAL_ARITH2_HPP + +#include <boost/config.hpp> +#include <boost/numeric/interval/detail/interval_prototype.hpp> +#include <boost/numeric/interval/detail/test_input.hpp> +#include <boost/numeric/interval/detail/bugs.hpp> +#include <boost/numeric/interval/detail/division.hpp> +#include <boost/numeric/interval/arith.hpp> +#include <boost/numeric/interval/policies.hpp> +#include <algorithm> +#include <cassert> +#include <boost/config/no_tr1/cmath.hpp> + +namespace boost { +namespace numeric { + +template<class T, class Policies> inline +interval<T, Policies> fmod(const interval<T, Policies>& x, + const interval<T, Policies>& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval<T, Policies>::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I; + T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper(); + T n = rnd.int_down(rnd.div_down(x.lower(), yb)); + return (const I&)x - n * (const I&)y; +} + +template<class T, class Policies> inline +interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval<T, Policies>::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I; + T n = rnd.int_down(rnd.div_down(x.lower(), y)); + return (const I&)x - n * I(y); +} + +template<class T, class Policies> inline +interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval<T, Policies>::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I; + T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper(); + T n = rnd.int_down(rnd.div_down(x, yb)); + return x - n * (const I&)y; +} + +namespace interval_lib { + +template<class T, class Policies> inline +interval<T, Policies> division_part1(const interval<T, Policies>& x, + const interval<T, Policies>& y, bool& b) +{ + typedef interval<T, Policies> I; + b = false; + if (detail::test_input(x, y)) + return I::empty(); + if (zero_in(y)) + if (!user::is_zero(y.lower())) + if (!user::is_zero(y.upper())) + return detail::div_zero_part1(x, y, b); + else + return detail::div_negative(x, y.lower()); + else + if (!user::is_zero(y.upper())) + return detail::div_positive(x, y.upper()); + else + return I::empty(); + else + return detail::div_non_zero(x, y); +} + +template<class T, class Policies> inline +interval<T, Policies> division_part2(const interval<T, Policies>& x, + const interval<T, Policies>& y, bool b = true) +{ + if (!b) return interval<T, Policies>::empty(); + return detail::div_zero_part2(x, y); +} + +template<class T, class Policies> inline +interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x) +{ + typedef interval<T, Policies> I; + if (detail::test_input(x)) + return I::empty(); + T one = static_cast<T>(1); + typename Policies::rounding rnd; + if (zero_in(x)) { + typedef typename Policies::checking checking; + if (!user::is_zero(x.lower())) + if (!user::is_zero(x.upper())) + return I::whole(); + else + return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true); + else + if (!user::is_zero(x.upper())) + return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true); + else + return I::empty(); + } else + return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true); +} + +namespace detail { + +template<class T, class Rounding> inline +T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive +{ + T x = x_; + T y = (pwr & 1) ? x_ : static_cast<T>(1); + pwr >>= 1; + while (pwr > 0) { + x = rnd.mul_down(x, x); + if (pwr & 1) y = rnd.mul_down(x, y); + pwr >>= 1; + } + return y; +} + +template<class T, class Rounding> inline +T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive +{ + T x = x_; + T y = (pwr & 1) ? x_ : static_cast<T>(1); + pwr >>= 1; + while (pwr > 0) { + x = rnd.mul_up(x, x); + if (pwr & 1) y = rnd.mul_up(x, y); + pwr >>= 1; + } + return y; +} + +} // namespace detail +} // namespace interval_lib + +template<class T, class Policies> inline +interval<T, Policies> pow(const interval<T, Policies>& x, int pwr) +{ + BOOST_USING_STD_MAX(); + using interval_lib::detail::pow_dn; + using interval_lib::detail::pow_up; + typedef interval<T, Policies> I; + + if (interval_lib::detail::test_input(x)) + return I::empty(); + + if (pwr == 0) + if (interval_lib::user::is_zero(x.lower()) + && interval_lib::user::is_zero(x.upper())) + return I::empty(); + else + return I(static_cast<T>(1)); + else if (pwr < 0) + return interval_lib::multiplicative_inverse(pow(x, -pwr)); + + typename Policies::rounding rnd; + + if (interval_lib::user::is_neg(x.upper())) { // [-2,-1] + T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd); + T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd); + if (pwr & 1) // [-2,-1]^1 + return I(-yu, -yl, true); + else // [-2,-1]^2 + return I(yl, yu, true); + } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1] + if (pwr & 1) { // [-1,1]^1 + return I(-pow_up(static_cast<T>(-x.lower()), pwr, rnd), pow_up(x.upper(), pwr, rnd), true); + } else { // [-1,1]^2 + return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true); + } + } else { // [1,2] + return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true); + } +} + +template<class T, class Policies> inline +interval<T, Policies> sqrt(const interval<T, Policies>& x) +{ + typedef interval<T, Policies> I; + if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper())) + return I::empty(); + typename Policies::rounding rnd; + T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower()); + return I(l, rnd.sqrt_up(x.upper()), true); +} + +template<class T, class Policies> inline +interval<T, Policies> square(const interval<T, Policies>& x) +{ + typedef interval<T, Policies> I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + const T& xl = x.lower(); + const T& xu = x.upper(); + if (interval_lib::user::is_neg(xu)) + return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true); + else if (interval_lib::user::is_pos(x.lower())) + return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true); + else + return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true); +} + +namespace interval_lib { +namespace detail { + +template< class I > inline +I root_aux(typename I::base_type const &x, int k) // x and k are bigger than one +{ + typedef typename I::base_type T; + T tk(k); + I y(static_cast<T>(1), x, true); + for(;;) { + T y0 = median(y); + I yy = intersect(y, y0 - (pow(I(y0, y0, true), k) - x) / (tk * pow(y, k - 1))); + if (equal(y, yy)) return y; + y = yy; + } +} + +template< class I > inline // x is positive and k bigger than one +typename I::base_type root_aux_dn(typename I::base_type const &x, int k) +{ + typedef typename I::base_type T; + typedef typename I::traits_type Policies; + typename Policies::rounding rnd; + T one(1); + if (x > one) return root_aux<I>(x, k).lower(); + if (x == one) return one; + return rnd.div_down(one, root_aux<I>(rnd.div_up(one, x), k).upper()); +} + +template< class I > inline // x is positive and k bigger than one +typename I::base_type root_aux_up(typename I::base_type const &x, int k) +{ + typedef typename I::base_type T; + typedef typename I::traits_type Policies; + typename Policies::rounding rnd; + T one(1); + if (x > one) return root_aux<I>(x, k).upper(); + if (x == one) return one; + return rnd.div_up(one, root_aux<I>(rnd.div_down(one, x), k).lower()); +} + +} // namespace detail +} // namespace interval_lib + +template< class T, class Policies > inline +interval<T, Policies> nth_root(interval<T, Policies> const &x, int k) +{ + typedef interval<T, Policies> I; + if (interval_lib::detail::test_input(x)) return I::empty(); + assert(k > 0); + if (k == 1) return x; + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect<I>::type R; + if (!interval_lib::user::is_pos(x.upper())) { + if (interval_lib::user::is_zero(x.upper())) { + T zero(0); + if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,0]^/2 or [0,0] + return I(zero, zero, true); + else // [-1,0]^/3 + return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), zero, true); + } else if (!(k & 1)) // [-2,-1]^/2 + return I::empty(); + else { // [-2,-1]^/3 + return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), + -interval_lib::detail::root_aux_dn<R>(-x.upper(), k), true); + } + } + T u = interval_lib::detail::root_aux_up<R>(x.upper(), k); + if (!interval_lib::user::is_pos(x.lower())) + if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,1]^/2 or [0,1] + return I(static_cast<T>(0), u, true); + else // [-1,1]^/3 + return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), u, true); + else // [1,2] + return I(interval_lib::detail::root_aux_dn<R>(x.lower(), k), u, true); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP