Something went wrong. Try again.
This repository has no description
Something went wrong. Try again.
10 kB · 305 lines
C++
at develop
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306/* 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 GECODE_BOOST_NUMERIC_INTERVAL_ARITH2_HPP#define GECODE_BOOST_NUMERIC_INTERVAL_ARITH2_HPP
#include <gecode/third-party/boost/config.hpp>#include <gecode/third-party/boost/numeric/interval/detail/interval_prototype.hpp>#include <gecode/third-party/boost/numeric/interval/detail/test_input.hpp>#include <gecode/third-party/boost/numeric/interval/detail/bugs.hpp>#include <gecode/third-party/boost/numeric/interval/detail/division.hpp>#include <gecode/third-party/boost/numeric/interval/arith.hpp>#include <gecode/third-party/boost/numeric/interval/policies.hpp>#include <algorithm>#include <cassert>#include <gecode/third-party/boost/config/no_tr1/cmath.hpp>
namespace gecode_boost {namespace numeric {
template<class T, class Policies> inlineinterval<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> inlineinterval<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> inlineinterval<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> inlineinterval<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> inlineinterval<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> inlineinterval<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> inlineT 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> inlineT 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> inlineinterval<T, Policies> pow(const interval<T, Policies>& x, int pwr){ GECODE_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 GECODE_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> inlineinterval<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> inlineinterval<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 > inlineI 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 onetypename 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 onetypename 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 > inlineinterval<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 gecode_boost
#endif // GECODE_BOOST_NUMERIC_INTERVAL_ARITH2_HPP