| 1 | /* Boost interval/arith2.hpp template implementation file
|
|---|
| 2 | *
|
|---|
| 3 | * This header provides some auxiliary arithmetic
|
|---|
| 4 | * functions: fmod, sqrt, square, pov, inverse and
|
|---|
| 5 | * a multi-interval division.
|
|---|
| 6 | *
|
|---|
| 7 | * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
|
|---|
| 8 | *
|
|---|
| 9 | * Distributed under the Boost Software License, Version 1.0.
|
|---|
| 10 | * (See accompanying file LICENSE_1_0.txt or
|
|---|
| 11 | * copy at http://www.boost.org/LICENSE_1_0.txt)
|
|---|
| 12 | */
|
|---|
| 13 |
|
|---|
| 14 | #ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
|
|---|
| 15 | #define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
|
|---|
| 16 |
|
|---|
| 17 | #include <boost/config.hpp>
|
|---|
| 18 | #include <boost/numeric/interval/detail/interval_prototype.hpp>
|
|---|
| 19 | #include <boost/numeric/interval/detail/test_input.hpp>
|
|---|
| 20 | #include <boost/numeric/interval/detail/bugs.hpp>
|
|---|
| 21 | #include <boost/numeric/interval/detail/division.hpp>
|
|---|
| 22 | #include <boost/numeric/interval/arith.hpp>
|
|---|
| 23 | #include <boost/numeric/interval/policies.hpp>
|
|---|
| 24 | #include <algorithm>
|
|---|
| 25 | #include <cmath>
|
|---|
| 26 |
|
|---|
| 27 | namespace boost {
|
|---|
| 28 | namespace numeric {
|
|---|
| 29 |
|
|---|
| 30 | template<class T, class Policies> inline
|
|---|
| 31 | interval<T, Policies> fmod(const interval<T, Policies>& x,
|
|---|
| 32 | const interval<T, Policies>& y)
|
|---|
| 33 | {
|
|---|
| 34 | if (interval_lib::detail::test_input(x, y))
|
|---|
| 35 | return interval<T, Policies>::empty();
|
|---|
| 36 | typename Policies::rounding rnd;
|
|---|
| 37 | typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
|
|---|
| 38 | T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
|
|---|
| 39 | T n = rnd.int_down(rnd.div_down(x.lower(), yb));
|
|---|
| 40 | return (const I&)x - n * (const I&)y;
|
|---|
| 41 | }
|
|---|
| 42 |
|
|---|
| 43 | template<class T, class Policies> inline
|
|---|
| 44 | interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
|
|---|
| 45 | {
|
|---|
| 46 | if (interval_lib::detail::test_input(x, y))
|
|---|
| 47 | return interval<T, Policies>::empty();
|
|---|
| 48 | typename Policies::rounding rnd;
|
|---|
| 49 | typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
|
|---|
| 50 | T n = rnd.int_down(rnd.div_down(x.lower(), y));
|
|---|
| 51 | return (const I&)x - n * I(y);
|
|---|
| 52 | }
|
|---|
| 53 |
|
|---|
| 54 | template<class T, class Policies> inline
|
|---|
| 55 | interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
|
|---|
| 56 | {
|
|---|
| 57 | if (interval_lib::detail::test_input(x, y))
|
|---|
| 58 | return interval<T, Policies>::empty();
|
|---|
| 59 | typename Policies::rounding rnd;
|
|---|
| 60 | typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
|
|---|
| 61 | T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
|
|---|
| 62 | T n = rnd.int_down(rnd.div_down(x, yb));
|
|---|
| 63 | return x - n * (const I&)y;
|
|---|
| 64 | }
|
|---|
| 65 |
|
|---|
| 66 | namespace interval_lib {
|
|---|
| 67 |
|
|---|
| 68 | template<class T, class Policies> inline
|
|---|
| 69 | interval<T, Policies> division_part1(const interval<T, Policies>& x,
|
|---|
| 70 | const interval<T, Policies>& y, bool& b)
|
|---|
| 71 | {
|
|---|
| 72 | typedef interval<T, Policies> I;
|
|---|
| 73 | b = false;
|
|---|
| 74 | if (detail::test_input(x, y))
|
|---|
| 75 | return I::empty();
|
|---|
| 76 | if (in_zero(y))
|
|---|
| 77 | if (!user::is_zero(y.lower()))
|
|---|
| 78 | if (!user::is_zero(y.upper()))
|
|---|
| 79 | return detail::div_zero_part1(x, y, b);
|
|---|
| 80 | else
|
|---|
| 81 | return detail::div_negative(x, y.lower());
|
|---|
| 82 | else
|
|---|
| 83 | if (!user::is_zero(y.upper()))
|
|---|
| 84 | return detail::div_positive(x, y.upper());
|
|---|
| 85 | else
|
|---|
| 86 | return I::empty();
|
|---|
| 87 | else
|
|---|
| 88 | return detail::div_non_zero(x, y);
|
|---|
| 89 | }
|
|---|
| 90 |
|
|---|
| 91 | template<class T, class Policies> inline
|
|---|
| 92 | interval<T, Policies> division_part2(const interval<T, Policies>& x,
|
|---|
| 93 | const interval<T, Policies>& y, bool b = true)
|
|---|
| 94 | {
|
|---|
| 95 | if (!b) return interval<T, Policies>::empty();
|
|---|
| 96 | return detail::div_zero_part2(x, y);
|
|---|
| 97 | }
|
|---|
| 98 |
|
|---|
| 99 | template<class T, class Policies> inline
|
|---|
| 100 | interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
|
|---|
| 101 | {
|
|---|
| 102 | typedef interval<T, Policies> I;
|
|---|
| 103 | if (detail::test_input(x))
|
|---|
| 104 | return I::empty();
|
|---|
| 105 | T one = static_cast<T>(1);
|
|---|
| 106 | typename Policies::rounding rnd;
|
|---|
| 107 | if (in_zero(x)) {
|
|---|
| 108 | typedef typename Policies::checking checking;
|
|---|
| 109 | if (!user::is_zero(x.lower()))
|
|---|
| 110 | if (!user::is_zero(x.upper()))
|
|---|
| 111 | return I::whole();
|
|---|
| 112 | else
|
|---|
| 113 | return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
|
|---|
| 114 | else
|
|---|
| 115 | if (!user::is_zero(x.upper()))
|
|---|
| 116 | return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
|
|---|
| 117 | else
|
|---|
| 118 | return I::empty();
|
|---|
| 119 | } else
|
|---|
| 120 | return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
|
|---|
| 121 | }
|
|---|
| 122 |
|
|---|
| 123 | namespace detail {
|
|---|
| 124 |
|
|---|
| 125 | template<class T, class Rounding> inline
|
|---|
| 126 | T pow_aux(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
|
|---|
| 127 | {
|
|---|
| 128 | T x = x_;
|
|---|
| 129 | T y = (pwr & 1) ? x_ : static_cast<T>(1);
|
|---|
| 130 | pwr >>= 1;
|
|---|
| 131 | while (pwr > 0) {
|
|---|
| 132 | x = rnd.mul_up(x, x);
|
|---|
| 133 | if (pwr & 1) y = rnd.mul_up(x, y);
|
|---|
| 134 | pwr >>= 1;
|
|---|
| 135 | }
|
|---|
| 136 | return y;
|
|---|
| 137 | }
|
|---|
| 138 |
|
|---|
| 139 | } // namespace detail
|
|---|
| 140 | } // namespace interval_lib
|
|---|
| 141 |
|
|---|
| 142 | template<class T, class Policies> inline
|
|---|
| 143 | interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
|
|---|
| 144 | {
|
|---|
| 145 | BOOST_USING_STD_MAX();
|
|---|
| 146 | using interval_lib::detail::pow_aux;
|
|---|
| 147 | typedef interval<T, Policies> I;
|
|---|
| 148 |
|
|---|
| 149 | if (interval_lib::detail::test_input(x))
|
|---|
| 150 | return I::empty();
|
|---|
| 151 |
|
|---|
| 152 | if (pwr == 0)
|
|---|
| 153 | if (interval_lib::user::is_zero(x.lower())
|
|---|
| 154 | && interval_lib::user::is_zero(x.upper()))
|
|---|
| 155 | return I::empty();
|
|---|
| 156 | else
|
|---|
| 157 | return I(static_cast<T>(1));
|
|---|
| 158 | else if (pwr < 0)
|
|---|
| 159 | return interval_lib::multiplicative_inverse(pow(x, -pwr));
|
|---|
| 160 |
|
|---|
| 161 | typename Policies::rounding rnd;
|
|---|
| 162 |
|
|---|
| 163 | if (interval_lib::user::is_neg(x.upper())) { // [-2,-1]
|
|---|
| 164 | T yl = pow_aux(-x.upper(), pwr, rnd);
|
|---|
| 165 | T yu = pow_aux(-x.lower(), pwr, rnd);
|
|---|
| 166 | if (pwr & 1) // [-2,-1]^1
|
|---|
| 167 | return I(-yu, -yl, true);
|
|---|
| 168 | else // [-2,-1]^2
|
|---|
| 169 | return I(yl, yu, true);
|
|---|
| 170 | } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
|
|---|
| 171 | if (pwr & 1) { // [-1,1]^1
|
|---|
| 172 | return I(-pow_aux(-x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true);
|
|---|
| 173 | } else { // [-1,1]^2
|
|---|
| 174 | return I(static_cast<T>(0), pow_aux(max BOOST_PREVENT_MACRO_SUBSTITUTION(-x.lower(), x.upper()), pwr, rnd), true);
|
|---|
| 175 | }
|
|---|
| 176 | } else { // [1,2]
|
|---|
| 177 | return I(pow_aux(x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true);
|
|---|
| 178 | }
|
|---|
| 179 | }
|
|---|
| 180 |
|
|---|
| 181 | template<class T, class Policies> inline
|
|---|
| 182 | interval<T, Policies> sqrt(const interval<T, Policies>& x)
|
|---|
| 183 | {
|
|---|
| 184 | typedef interval<T, Policies> I;
|
|---|
| 185 | if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
|
|---|
| 186 | return I::empty();
|
|---|
| 187 | typename Policies::rounding rnd;
|
|---|
| 188 | T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
|
|---|
| 189 | return I(l, rnd.sqrt_up(x.upper()), true);
|
|---|
| 190 | }
|
|---|
| 191 |
|
|---|
| 192 | template<class T, class Policies> inline
|
|---|
| 193 | interval<T, Policies> square(const interval<T, Policies>& x)
|
|---|
| 194 | {
|
|---|
| 195 | typedef interval<T, Policies> I;
|
|---|
| 196 | if (interval_lib::detail::test_input(x))
|
|---|
| 197 | return I::empty();
|
|---|
| 198 | typename Policies::rounding rnd;
|
|---|
| 199 | const T& xl = x.lower();
|
|---|
| 200 | const T& xu = x.upper();
|
|---|
| 201 | if (interval_lib::user::is_neg(xu))
|
|---|
| 202 | return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
|
|---|
| 203 | else if (interval_lib::user::is_pos(x.lower()))
|
|---|
| 204 | return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
|
|---|
| 205 | else
|
|---|
| 206 | return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | } // namespace numeric
|
|---|
| 210 | } // namespace boost
|
|---|
| 211 |
|
|---|
| 212 | #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/gametools/chrome/site/your_project_logo.png){kind=logo}
