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 <cassert> |
---|
26 | #include <boost/config/no_tr1/cmath.hpp> |
---|
27 | |
---|
28 | namespace boost { |
---|
29 | namespace numeric { |
---|
30 | |
---|
31 | template<class T, class Policies> inline |
---|
32 | interval<T, Policies> fmod(const interval<T, Policies>& x, |
---|
33 | const interval<T, Policies>& y) |
---|
34 | { |
---|
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; |
---|
42 | } |
---|
43 | |
---|
44 | template<class T, class Policies> inline |
---|
45 | interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y) |
---|
46 | { |
---|
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); |
---|
53 | } |
---|
54 | |
---|
55 | template<class T, class Policies> inline |
---|
56 | interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y) |
---|
57 | { |
---|
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; |
---|
65 | } |
---|
66 | |
---|
67 | namespace interval_lib { |
---|
68 | |
---|
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) |
---|
72 | { |
---|
73 | typedef interval<T, Policies> I; |
---|
74 | b = false; |
---|
75 | if (detail::test_input(x, y)) |
---|
76 | return I::empty(); |
---|
77 | if (zero_in(y)) |
---|
78 | if (!user::is_zero(y.lower())) |
---|
79 | if (!user::is_zero(y.upper())) |
---|
80 | return detail::div_zero_part1(x, y, b); |
---|
81 | else |
---|
82 | return detail::div_negative(x, y.lower()); |
---|
83 | else |
---|
84 | if (!user::is_zero(y.upper())) |
---|
85 | return detail::div_positive(x, y.upper()); |
---|
86 | else |
---|
87 | return I::empty(); |
---|
88 | else |
---|
89 | return detail::div_non_zero(x, y); |
---|
90 | } |
---|
91 | |
---|
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) |
---|
95 | { |
---|
96 | if (!b) return interval<T, Policies>::empty(); |
---|
97 | return detail::div_zero_part2(x, y); |
---|
98 | } |
---|
99 | |
---|
100 | template<class T, class Policies> inline |
---|
101 | interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x) |
---|
102 | { |
---|
103 | typedef interval<T, Policies> I; |
---|
104 | if (detail::test_input(x)) |
---|
105 | return I::empty(); |
---|
106 | T one = static_cast<T>(1); |
---|
107 | typename Policies::rounding rnd; |
---|
108 | if (zero_in(x)) { |
---|
109 | typedef typename Policies::checking checking; |
---|
110 | if (!user::is_zero(x.lower())) |
---|
111 | if (!user::is_zero(x.upper())) |
---|
112 | return I::whole(); |
---|
113 | else |
---|
114 | return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true); |
---|
115 | else |
---|
116 | if (!user::is_zero(x.upper())) |
---|
117 | return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true); |
---|
118 | else |
---|
119 | return I::empty(); |
---|
120 | } else |
---|
121 | return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true); |
---|
122 | } |
---|
123 | |
---|
124 | namespace detail { |
---|
125 | |
---|
126 | template<class T, class Rounding> inline |
---|
127 | T pow_dn(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive |
---|
128 | { |
---|
129 | T x = x_; |
---|
130 | T y = (pwr & 1) ? x_ : static_cast<T>(1); |
---|
131 | pwr >>= 1; |
---|
132 | while (pwr > 0) { |
---|
133 | x = rnd.mul_down(x, x); |
---|
134 | if (pwr & 1) y = rnd.mul_down(x, y); |
---|
135 | pwr >>= 1; |
---|
136 | } |
---|
137 | return y; |
---|
138 | } |
---|
139 | |
---|
140 | template<class T, class Rounding> inline |
---|
141 | T pow_up(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive |
---|
142 | { |
---|
143 | T x = x_; |
---|
144 | T y = (pwr & 1) ? x_ : static_cast<T>(1); |
---|
145 | pwr >>= 1; |
---|
146 | while (pwr > 0) { |
---|
147 | x = rnd.mul_up(x, x); |
---|
148 | if (pwr & 1) y = rnd.mul_up(x, y); |
---|
149 | pwr >>= 1; |
---|
150 | } |
---|
151 | return y; |
---|
152 | } |
---|
153 | |
---|
154 | } // namespace detail |
---|
155 | } // namespace interval_lib |
---|
156 | |
---|
157 | template<class T, class Policies> inline |
---|
158 | interval<T, Policies> pow(const interval<T, Policies>& x, int pwr) |
---|
159 | { |
---|
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; |
---|
164 | |
---|
165 | if (interval_lib::detail::test_input(x)) |
---|
166 | return I::empty(); |
---|
167 | |
---|
168 | if (pwr == 0) |
---|
169 | if (interval_lib::user::is_zero(x.lower()) |
---|
170 | && interval_lib::user::is_zero(x.upper())) |
---|
171 | return I::empty(); |
---|
172 | else |
---|
173 | return I(static_cast<T>(1)); |
---|
174 | else if (pwr < 0) |
---|
175 | return interval_lib::multiplicative_inverse(pow(x, -pwr)); |
---|
176 | |
---|
177 | typename Policies::rounding rnd; |
---|
178 | |
---|
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); |
---|
184 | else // [-2,-1]^2 |
---|
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(-x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true); |
---|
189 | } else { // [-1,1]^2 |
---|
190 | return I(static_cast<T>(0), pow_up(max BOOST_PREVENT_MACRO_SUBSTITUTION(static_cast<T>(-x.lower()), x.upper()), pwr, rnd), true); |
---|
191 | } |
---|
192 | } else { // [1,2] |
---|
193 | return I(pow_dn(x.lower(), pwr, rnd), pow_up(x.upper(), pwr, rnd), true); |
---|
194 | } |
---|
195 | } |
---|
196 | |
---|
197 | template<class T, class Policies> inline |
---|
198 | interval<T, Policies> sqrt(const interval<T, Policies>& x) |
---|
199 | { |
---|
200 | typedef interval<T, Policies> I; |
---|
201 | if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper())) |
---|
202 | return I::empty(); |
---|
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); |
---|
206 | } |
---|
207 | |
---|
208 | template<class T, class Policies> inline |
---|
209 | interval<T, Policies> square(const interval<T, Policies>& x) |
---|
210 | { |
---|
211 | typedef interval<T, Policies> I; |
---|
212 | if (interval_lib::detail::test_input(x)) |
---|
213 | return I::empty(); |
---|
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); |
---|
221 | else |
---|
222 | return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true); |
---|
223 | } |
---|
224 | |
---|
225 | namespace interval_lib { |
---|
226 | namespace detail { |
---|
227 | |
---|
228 | template< class I > inline |
---|
229 | I root_aux(typename I::base_type const &x, int k) // x and k are bigger than one |
---|
230 | { |
---|
231 | typedef typename I::base_type T; |
---|
232 | T tk(k); |
---|
233 | I y(static_cast<T>(1), x, true); |
---|
234 | for(;;) { |
---|
235 | T y0 = median(y); |
---|
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; |
---|
238 | y = yy; |
---|
239 | } |
---|
240 | } |
---|
241 | |
---|
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) |
---|
244 | { |
---|
245 | typedef typename I::base_type T; |
---|
246 | typedef typename I::traits_type Policies; |
---|
247 | typename Policies::rounding rnd; |
---|
248 | T one(1); |
---|
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()); |
---|
252 | } |
---|
253 | |
---|
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) |
---|
256 | { |
---|
257 | typedef typename I::base_type T; |
---|
258 | typedef typename I::traits_type Policies; |
---|
259 | typename Policies::rounding rnd; |
---|
260 | T one(1); |
---|
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()); |
---|
264 | } |
---|
265 | |
---|
266 | } // namespace detail |
---|
267 | } // namespace interval_lib |
---|
268 | |
---|
269 | template< class T, class Policies > inline |
---|
270 | interval<T, Policies> nth_root(interval<T, Policies> const &x, int k) |
---|
271 | { |
---|
272 | typedef interval<T, Policies> I; |
---|
273 | if (interval_lib::detail::test_input(x)) return I::empty(); |
---|
274 | assert(k > 0); |
---|
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())) { |
---|
280 | T zero(0); |
---|
281 | if (!(k & 1) || interval_lib::user::is_zero(x.lower())) // [-1,0]^/2 or [0,0] |
---|
282 | return I(zero, zero, true); |
---|
283 | else // [-1,0]^/3 |
---|
284 | return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), zero, true); |
---|
285 | } else if (!(k & 1)) // [-2,-1]^/2 |
---|
286 | return I::empty(); |
---|
287 | else { // [-2,-1]^/3 |
---|
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); |
---|
290 | } |
---|
291 | } |
---|
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); |
---|
296 | else // [-1,1]^/3 |
---|
297 | return I(-interval_lib::detail::root_aux_up<R>(-x.lower(), k), u, true); |
---|
298 | else // [1,2] |
---|
299 | return I(interval_lib::detail::root_aux_dn<R>(x.lower(), k), u, true); |
---|
300 | } |
---|
301 | |
---|
302 | } // namespace numeric |
---|
303 | } // namespace boost |
---|
304 | |
---|
305 | #endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP |
---|