Mercurial > hg > vamp-build-and-test
comparison DEPENDENCIES/generic/include/boost/numeric/interval/arith2.hpp @ 16:2665513ce2d3
Add boost headers
author | Chris Cannam |
---|---|
date | Tue, 05 Aug 2014 11:11:38 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
15:663ca0da4350 | 16:2665513ce2d3 |
---|---|
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(static_cast<T>(-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 |