1 /* Boost interval/arith2.hpp template implementation file
3 * This header provides some auxiliary arithmetic
4 * functions: fmod, sqrt, square, pov, inverse and
5 * a multi-interval division.
7 * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
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)
14 #ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
15 #define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
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>
26 #include <boost/config/no_tr1/cmath.hpp>
31 template<class T, class Policies> inline
32 interval<T, Policies> fmod(const interval<T, Policies>& x,
33 const interval<T, Policies>& y)
35 if (interval_lib::detail::test_input(x, y))
36 return interval<T, Policies>::empty();
37 typename Policies::rounding rnd;
38 typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
39 T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
40 T n = rnd.int_down(rnd.div_down(x.lower(), yb));
41 return (const I&)x - n * (const I&)y;
44 template<class T, class Policies> inline
45 interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
47 if (interval_lib::detail::test_input(x, y))
48 return interval<T, Policies>::empty();
49 typename Policies::rounding rnd;
50 typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
51 T n = rnd.int_down(rnd.div_down(x.lower(), y));
52 return (const I&)x - n * I(y);
55 template<class T, class Policies> inline
56 interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
58 if (interval_lib::detail::test_input(x, y))
59 return interval<T, Policies>::empty();
60 typename Policies::rounding rnd;
61 typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
62 T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
63 T n = rnd.int_down(rnd.div_down(x, yb));
64 return x - n * (const I&)y;
67 namespace interval_lib {
69 template<class T, class Policies> inline
70 interval<T, Policies> division_part1(const interval<T, Policies>& x,
71 const interval<T, Policies>& y, bool& b)
73 typedef interval<T, Policies> I;
75 if (detail::test_input(x, y))
78 if (!user::is_zero(y.lower()))
79 if (!user::is_zero(y.upper()))
80 return detail::div_zero_part1(x, y, b);
82 return detail::div_negative(x, y.lower());
84 if (!user::is_zero(y.upper()))
85 return detail::div_positive(x, y.upper());
89 return detail::div_non_zero(x, y);
92 template<class T, class Policies> inline
93 interval<T, Policies> division_part2(const interval<T, Policies>& x,
94 const interval<T, Policies>& y, bool b = true)
96 if (!b) return interval<T, Policies>::empty();
97 return detail::div_zero_part2(x, y);
100 template<class T, class Policies> inline
101 interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
103 typedef interval<T, Policies> I;
104 if (detail::test_input(x))
106 T one = static_cast<T>(1);
107 typename Policies::rounding rnd;
109 typedef typename Policies::checking checking;
110 if (!user::is_zero(x.lower()))
111 if (!user::is_zero(x.upper()))
114 return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
116 if (!user::is_zero(x.upper()))
117 return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
121 return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
126 template<class T, class Rounding> inline
127 T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
130 T y = (pwr & 1) ? x_ : static_cast<T>(1);
133 x = rnd.mul_down(x, x);
134 if (pwr & 1) y = rnd.mul_down(x, y);
140 template<class T, class Rounding> inline
141 T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
144 T y = (pwr & 1) ? x_ : static_cast<T>(1);
147 x = rnd.mul_up(x, x);
148 if (pwr & 1) y = rnd.mul_up(x, y);
154 } // namespace detail
155 } // namespace interval_lib
157 template<class T, class Policies> inline
158 interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
160 BOOST_USING_STD_MAX();
161 using interval_lib::detail::pow_dn;
162 using interval_lib::detail::pow_up;
163 typedef interval<T, Policies> I;
165 if (interval_lib::detail::test_input(x))
169 if (interval_lib::user::is_zero(x.lower())
170 && interval_lib::user::is_zero(x.upper()))
173 return I(static_cast<T>(1));
175 return interval_lib::multiplicative_inverse(pow(x, -pwr));
177 typename Policies::rounding rnd;
179 if (interval_lib::user::is_neg(x.upper())) { // [-2,-1]
180 T yl = pow_dn(static_cast<T>(-x.upper()), pwr, rnd);
181 T yu = pow_up(static_cast<T>(-x.lower()), pwr, rnd);
182 if (pwr & 1) // [-2,-1]^1
183 return I(-yu, -yl, true);
185 return I(yl, yu, true);
186 } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
187 if (pwr & 1) { // [-1,1]^1
188 return I(-pow_up(static_cast<T>(-x.lower()), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
190 return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true);
193 return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true);
197 template<class T, class Policies> inline
198 interval<T, Policies> sqrt(const interval<T, Policies>& x)
200 typedef interval<T, Policies> I;
201 if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
203 typename Policies::rounding rnd;
204 T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
205 return I(l, rnd.sqrt_up(x.upper()), true);
208 template<class T, class Policies> inline
209 interval<T, Policies> square(const interval<T, Policies>& x)
211 typedef interval<T, Policies> I;
212 if (interval_lib::detail::test_input(x))
214 typename Policies::rounding rnd;
215 const T& xl = x.lower();
216 const T& xu = x.upper();
217 if (interval_lib::user::is_neg(xu))
218 return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
219 else if (interval_lib::user::is_pos(x.lower()))
220 return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
222 return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
225 namespace interval_lib {
228 template< class I > inline
229 I root_aux(typename I::base_type const &x, int k) // x and k are bigger than one
231 typedef typename I::base_type T;
233 I y(static_cast<T>(1), x, true);
236 I yy = intersect(y, y0 - (pow(I(y0, y0, true), k) - x) / (tk * pow(y, k - 1)));
237 if (equal(y, yy)) return y;
242 template< class I > inline // x is positive and k bigger than one
243 typename I::base_type root_aux_dn(typename I::base_type const &x, int k)
245 typedef typename I::base_type T;
246 typedef typename I::traits_type Policies;
247 typename Policies::rounding rnd;
249 if (x > one) return root_aux<I>(x, k).lower();
250 if (x == one) return one;
251 return rnd.div_down(one, root_aux<I>(rnd.div_up(one, x), k).upper());
254 template< class I > inline // x is positive and k bigger than one
255 typename I::base_type root_aux_up(typename I::base_type const &x, int k)
257 typedef typename I::base_type T;
258 typedef typename I::traits_type Policies;
259 typename Policies::rounding rnd;
261 if (x > one) return root_aux<I>(x, k).upper();
262 if (x == one) return one;
263 return rnd.div_up(one, root_aux<I>(rnd.div_down(one, x), k).lower());
266 } // namespace detail
267 } // namespace interval_lib
269 template< class T, class Policies > inline
270 interval<T, Policies> nth_root(interval<T, Policies> const &x, int k)
272 typedef interval<T, Policies> I;
273 if (interval_lib::detail::test_input(x)) return I::empty();
275 if (k == 1) return x;
276 typename Policies::rounding rnd;
277 typedef typename interval_lib::unprotect<I>::type R;
278 if (!interval_lib::user::is_pos(x.upper())) {
279 if (interval_lib::user::is_zero(x.upper())) {
281 if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,0]^/2 or [0,0]
282 return I(zero, zero, true);
284 return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), zero, true);
285 } else if (!(k & 1)) // [-2,-1]^/2
288 return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k),
289 -interval_lib::detail::root_aux_dn<R>(-x.upper(), k), true);
292 T u = interval_lib::detail::root_aux_up<R>(x.upper(), k);
293 if (!interval_lib::user::is_pos(x.lower()))
294 if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,1]^/2 or [0,1]
295 return I(static_cast<T>(0), u, true);
297 return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), u, true);
299 return I(interval_lib::detail::root_aux_dn<R>(x.lower(), k), u, true);
302 } // namespace numeric
305 #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP