1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/interval/arith2.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,227 @@
1.4 +/* Boost interval/arith2.hpp template implementation file
1.5 + *
1.6 + * This header provides some auxiliary arithmetic
1.7 + * functions: fmod, sqrt, square, pov, inverse and
1.8 + * a multi-interval division.
1.9 + *
1.10 + * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
1.11 + *
1.12 + * Distributed under the Boost Software License, Version 1.0.
1.13 + * (See accompanying file LICENSE_1_0.txt or
1.14 + * copy at http://www.boost.org/LICENSE_1_0.txt)
1.15 + */
1.16 +
1.17 +#ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
1.18 +#define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
1.19 +
1.20 +#include <boost/config.hpp>
1.21 +#include <boost/numeric/interval/detail/interval_prototype.hpp>
1.22 +#include <boost/numeric/interval/detail/test_input.hpp>
1.23 +#include <boost/numeric/interval/detail/bugs.hpp>
1.24 +#include <boost/numeric/interval/detail/division.hpp>
1.25 +#include <boost/numeric/interval/arith.hpp>
1.26 +#include <boost/numeric/interval/policies.hpp>
1.27 +#include <algorithm>
1.28 +#include <cmath>
1.29 +
1.30 +namespace boost {
1.31 +namespace numeric {
1.32 +
1.33 +template<class T, class Policies> inline
1.34 +interval<T, Policies> fmod(const interval<T, Policies>& x,
1.35 + const interval<T, Policies>& y)
1.36 +{
1.37 + if (interval_lib::detail::test_input(x, y))
1.38 + return interval<T, Policies>::empty();
1.39 + typename Policies::rounding rnd;
1.40 + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
1.41 + T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
1.42 + T n = rnd.int_down(rnd.div_down(x.lower(), yb));
1.43 + return (const I&)x - n * (const I&)y;
1.44 +}
1.45 +
1.46 +template<class T, class Policies> inline
1.47 +interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
1.48 +{
1.49 + if (interval_lib::detail::test_input(x, y))
1.50 + return interval<T, Policies>::empty();
1.51 + typename Policies::rounding rnd;
1.52 + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
1.53 + T n = rnd.int_down(rnd.div_down(x.lower(), y));
1.54 + return (const I&)x - n * I(y);
1.55 +}
1.56 +
1.57 +template<class T, class Policies> inline
1.58 +interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
1.59 +{
1.60 + if (interval_lib::detail::test_input(x, y))
1.61 + return interval<T, Policies>::empty();
1.62 + typename Policies::rounding rnd;
1.63 + typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
1.64 + T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
1.65 + T n = rnd.int_down(rnd.div_down(x, yb));
1.66 + return x - n * (const I&)y;
1.67 +}
1.68 +
1.69 +namespace interval_lib {
1.70 +
1.71 +template<class T, class Policies> inline
1.72 +interval<T, Policies> division_part1(const interval<T, Policies>& x,
1.73 + const interval<T, Policies>& y, bool& b)
1.74 +{
1.75 + typedef interval<T, Policies> I;
1.76 + b = false;
1.77 + if (detail::test_input(x, y))
1.78 + return I::empty();
1.79 + if (zero_in(y))
1.80 + if (!user::is_zero(y.lower()))
1.81 + if (!user::is_zero(y.upper()))
1.82 + return detail::div_zero_part1(x, y, b);
1.83 + else
1.84 + return detail::div_negative(x, y.lower());
1.85 + else
1.86 + if (!user::is_zero(y.upper()))
1.87 + return detail::div_positive(x, y.upper());
1.88 + else
1.89 + return I::empty();
1.90 + else
1.91 + return detail::div_non_zero(x, y);
1.92 +}
1.93 +
1.94 +template<class T, class Policies> inline
1.95 +interval<T, Policies> division_part2(const interval<T, Policies>& x,
1.96 + const interval<T, Policies>& y, bool b = true)
1.97 +{
1.98 + if (!b) return interval<T, Policies>::empty();
1.99 + return detail::div_zero_part2(x, y);
1.100 +}
1.101 +
1.102 +template<class T, class Policies> inline
1.103 +interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
1.104 +{
1.105 + typedef interval<T, Policies> I;
1.106 + if (detail::test_input(x))
1.107 + return I::empty();
1.108 + T one = static_cast<T>(1);
1.109 + typename Policies::rounding rnd;
1.110 + if (zero_in(x)) {
1.111 + typedef typename Policies::checking checking;
1.112 + if (!user::is_zero(x.lower()))
1.113 + if (!user::is_zero(x.upper()))
1.114 + return I::whole();
1.115 + else
1.116 + return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
1.117 + else
1.118 + if (!user::is_zero(x.upper()))
1.119 + return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
1.120 + else
1.121 + return I::empty();
1.122 + } else
1.123 + return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
1.124 +}
1.125 +
1.126 +namespace detail {
1.127 +
1.128 +template<class T, class Rounding> inline
1.129 +T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
1.130 +{
1.131 + T x = x_;
1.132 + T y = (pwr & 1) ? x_ : static_cast<T>(1);
1.133 + pwr >>= 1;
1.134 + while (pwr > 0) {
1.135 + x = rnd.mul_down(x, x);
1.136 + if (pwr & 1) y = rnd.mul_down(x, y);
1.137 + pwr >>= 1;
1.138 + }
1.139 + return y;
1.140 +}
1.141 +
1.142 +template<class T, class Rounding> inline
1.143 +T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
1.144 +{
1.145 + T x = x_;
1.146 + T y = (pwr & 1) ? x_ : static_cast<T>(1);
1.147 + pwr >>= 1;
1.148 + while (pwr > 0) {
1.149 + x = rnd.mul_up(x, x);
1.150 + if (pwr & 1) y = rnd.mul_up(x, y);
1.151 + pwr >>= 1;
1.152 + }
1.153 + return y;
1.154 +}
1.155 +
1.156 +} // namespace detail
1.157 +} // namespace interval_lib
1.158 +
1.159 +template<class T, class Policies> inline
1.160 +interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
1.161 +{
1.162 + BOOST_USING_STD_MAX();
1.163 + using interval_lib::detail::pow_dn;
1.164 + using interval_lib::detail::pow_up;
1.165 + typedef interval<T, Policies> I;
1.166 +
1.167 + if (interval_lib::detail::test_input(x))
1.168 + return I::empty();
1.169 +
1.170 + if (pwr == 0)
1.171 + if (interval_lib::user::is_zero(x.lower())
1.172 + && interval_lib::user::is_zero(x.upper()))
1.173 + return I::empty();
1.174 + else
1.175 + return I(static_cast<T>(1));
1.176 + else if (pwr < 0)
1.177 + return interval_lib::multiplicative_inverse(pow(x, -pwr));
1.178 +
1.179 + typename Policies::rounding rnd;
1.180 +
1.181 + if (interval_lib::user::is_neg(x.upper())) { // [-2,-1]
1.182 + T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);
1.183 + T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);
1.184 + if (pwr & 1) // [-2,-1]^1
1.185 + return I(-yu, -yl, true);
1.186 + else // [-2,-1]^2
1.187 + return I(yl, yu, true);
1.188 + } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
1.189 + if (pwr & 1) { // [-1,1]^1
1.190 + return I(-pow_up(-x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
1.191 + } else { // [-1,1]^2
1.192 + return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);
1.193 + }
1.194 + } else { // [1,2]
1.195 + return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
1.196 + }
1.197 +}
1.198 +
1.199 +template<class T, class Policies> inline
1.200 +interval<T, Policies> sqrt(const interval<T, Policies>& x)
1.201 +{
1.202 + typedef interval<T, Policies> I;
1.203 + if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
1.204 + return I::empty();
1.205 + typename Policies::rounding rnd;
1.206 + T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
1.207 + return I(l, rnd.sqrt_up(x.upper()), true);
1.208 +}
1.209 +
1.210 +template<class T, class Policies> inline
1.211 +interval<T, Policies> square(const interval<T, Policies>& x)
1.212 +{
1.213 + typedef interval<T, Policies> I;
1.214 + if (interval_lib::detail::test_input(x))
1.215 + return I::empty();
1.216 + typename Policies::rounding rnd;
1.217 + const T& xl = x.lower();
1.218 + const T& xu = x.upper();
1.219 + if (interval_lib::user::is_neg(xu))
1.220 + return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
1.221 + else if (interval_lib::user::is_pos(x.lower()))
1.222 + return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
1.223 + else
1.224 + return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
1.225 +}
1.226 +
1.227 +} // namespace numeric
1.228 +} // namespace boost
1.229 +
1.230 +#endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP