1 | // (C) Copyright John Maddock 2008. |
2 | // Use, modification and distribution are subject to the |
3 | // Boost Software License, Version 1.0. (See accompanying file |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
5 | |
6 | #ifndef BOOST_MATH_SPECIAL_NEXT_HPP |
7 | #define BOOST_MATH_SPECIAL_NEXT_HPP |
8 | |
9 | #ifdef _MSC_VER |
10 | #pragma once |
11 | #endif |
12 | |
13 | #include <boost/math/special_functions/math_fwd.hpp> |
14 | #include <boost/math/policies/error_handling.hpp> |
15 | #include <boost/math/special_functions/fpclassify.hpp> |
16 | #include <boost/math/special_functions/sign.hpp> |
17 | #include <boost/math/special_functions/trunc.hpp> |
18 | |
19 | #include <float.h> |
20 | |
21 | #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) |
22 | #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) |
23 | #include "xmmintrin.h" |
24 | #define BOOST_MATH_CHECK_SSE2 |
25 | #endif |
26 | #endif |
27 | |
28 | namespace boost{ namespace math{ |
29 | |
30 | namespace concepts { |
31 | |
32 | class real_concept; |
33 | class std_real_concept; |
34 | |
35 | } |
36 | |
37 | namespace detail{ |
38 | |
39 | template <class T> |
40 | struct has_hidden_guard_digits; |
41 | template <> |
42 | struct has_hidden_guard_digits<float> : public boost::false_type {}; |
43 | template <> |
44 | struct has_hidden_guard_digits<double> : public boost::false_type {}; |
45 | template <> |
46 | struct has_hidden_guard_digits<long double> : public boost::false_type {}; |
47 | #ifdef BOOST_HAS_FLOAT128 |
48 | template <> |
49 | struct has_hidden_guard_digits<__float128> : public boost::false_type {}; |
50 | #endif |
51 | template <> |
52 | struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public boost::false_type {}; |
53 | template <> |
54 | struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public boost::false_type {}; |
55 | |
56 | template <class T, bool b> |
57 | struct has_hidden_guard_digits_10 : public boost::false_type {}; |
58 | template <class T> |
59 | struct has_hidden_guard_digits_10<T, true> : public boost::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {}; |
60 | |
61 | template <class T> |
62 | struct has_hidden_guard_digits |
63 | : public has_hidden_guard_digits_10<T, |
64 | std::numeric_limits<T>::is_specialized |
65 | && (std::numeric_limits<T>::radix == 10) > |
66 | {}; |
67 | |
68 | template <class T> |
69 | inline const T& normalize_value(const T& val, const boost::false_type&) { return val; } |
70 | template <class T> |
71 | inline T normalize_value(const T& val, const boost::true_type&) |
72 | { |
73 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
74 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
75 | |
76 | boost::intmax_t shift = (boost::intmax_t)std::numeric_limits<T>::digits - (boost::intmax_t)ilogb(val) - 1; |
77 | T result = scalbn(val, shift); |
78 | result = round(result); |
79 | return scalbn(result, -shift); |
80 | } |
81 | |
82 | template <class T> |
83 | inline T get_smallest_value(boost::true_type const&) |
84 | { |
85 | // |
86 | // numeric_limits lies about denorms being present - particularly |
87 | // when this can be turned on or off at runtime, as is the case |
88 | // when using the SSE2 registers in DAZ or FTZ mode. |
89 | // |
90 | static const T m = std::numeric_limits<T>::denorm_min(); |
91 | #ifdef BOOST_MATH_CHECK_SSE2 |
92 | return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;; |
93 | #else |
94 | return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m; |
95 | #endif |
96 | } |
97 | |
98 | template <class T> |
99 | inline T get_smallest_value(boost::false_type const&) |
100 | { |
101 | return tools::min_value<T>(); |
102 | } |
103 | |
104 | template <class T> |
105 | inline T get_smallest_value() |
106 | { |
107 | #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310) |
108 | return get_smallest_value<T>(boost::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>()); |
109 | #else |
110 | return get_smallest_value<T>(boost::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>()); |
111 | #endif |
112 | } |
113 | |
114 | // |
115 | // Returns the smallest value that won't generate denorms when |
116 | // we calculate the value of the least-significant-bit: |
117 | // |
118 | template <class T> |
119 | T get_min_shift_value(); |
120 | |
121 | template <class T> |
122 | struct min_shift_initializer |
123 | { |
124 | struct init |
125 | { |
126 | init() |
127 | { |
128 | do_init(); |
129 | } |
130 | static void do_init() |
131 | { |
132 | get_min_shift_value<T>(); |
133 | } |
134 | void force_instantiate()const{} |
135 | }; |
136 | static const init initializer; |
137 | static void force_instantiate() |
138 | { |
139 | initializer.force_instantiate(); |
140 | } |
141 | }; |
142 | |
143 | template <class T> |
144 | const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer; |
145 | |
146 | template <class T> |
147 | inline T calc_min_shifted(const boost::true_type&) |
148 | { |
149 | BOOST_MATH_STD_USING |
150 | return ldexp(tools::min_value<T>(), tools::digits<T>() + 1); |
151 | } |
152 | template <class T> |
153 | inline T calc_min_shifted(const boost::false_type&) |
154 | { |
155 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
156 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
157 | |
158 | return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1); |
159 | } |
160 | |
161 | |
162 | template <class T> |
163 | inline T get_min_shift_value() |
164 | { |
165 | static const T val = calc_min_shifted<T>(boost::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>()); |
166 | min_shift_initializer<T>::force_instantiate(); |
167 | |
168 | return val; |
169 | } |
170 | |
171 | template <class T, class Policy> |
172 | T float_next_imp(const T& val, const boost::true_type&, const Policy& pol) |
173 | { |
174 | BOOST_MATH_STD_USING |
175 | int expon; |
176 | static const char* function = "float_next<%1%>(%1%)" ; |
177 | |
178 | int fpclass = (boost::math::fpclassify)(val); |
179 | |
180 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
181 | { |
182 | if(val < 0) |
183 | return -tools::max_value<T>(); |
184 | return policies::raise_domain_error<T>( |
185 | function, |
186 | "Argument must be finite, but got %1%" , val, pol); |
187 | } |
188 | |
189 | if(val >= tools::max_value<T>()) |
190 | return policies::raise_overflow_error<T>(function, 0, pol); |
191 | |
192 | if(val == 0) |
193 | return detail::get_smallest_value<T>(); |
194 | |
195 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>())) |
196 | { |
197 | // |
198 | // Special case: if the value of the least significant bit is a denorm, and the result |
199 | // would not be a denorm, then shift the input, increment, and shift back. |
200 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
201 | // |
202 | return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); |
203 | } |
204 | |
205 | if(-0.5f == frexp(val, &expon)) |
206 | --expon; // reduce exponent when val is a power of two, and negative. |
207 | T diff = ldexp(T(1), expon - tools::digits<T>()); |
208 | if(diff == 0) |
209 | diff = detail::get_smallest_value<T>(); |
210 | return val + diff; |
211 | } // float_next_imp |
212 | // |
213 | // Special version for some base other than 2: |
214 | // |
215 | template <class T, class Policy> |
216 | T float_next_imp(const T& val, const boost::false_type&, const Policy& pol) |
217 | { |
218 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
219 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
220 | |
221 | BOOST_MATH_STD_USING |
222 | boost::intmax_t expon; |
223 | static const char* function = "float_next<%1%>(%1%)" ; |
224 | |
225 | int fpclass = (boost::math::fpclassify)(val); |
226 | |
227 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
228 | { |
229 | if(val < 0) |
230 | return -tools::max_value<T>(); |
231 | return policies::raise_domain_error<T>( |
232 | function, |
233 | "Argument must be finite, but got %1%" , val, pol); |
234 | } |
235 | |
236 | if(val >= tools::max_value<T>()) |
237 | return policies::raise_overflow_error<T>(function, 0, pol); |
238 | |
239 | if(val == 0) |
240 | return detail::get_smallest_value<T>(); |
241 | |
242 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>())) |
243 | { |
244 | // |
245 | // Special case: if the value of the least significant bit is a denorm, and the result |
246 | // would not be a denorm, then shift the input, increment, and shift back. |
247 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
248 | // |
249 | return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits); |
250 | } |
251 | |
252 | expon = 1 + ilogb(val); |
253 | if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix) |
254 | --expon; // reduce exponent when val is a power of base, and negative. |
255 | T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits); |
256 | if(diff == 0) |
257 | diff = detail::get_smallest_value<T>(); |
258 | return val + diff; |
259 | } // float_next_imp |
260 | |
261 | } // namespace detail |
262 | |
263 | template <class T, class Policy> |
264 | inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol) |
265 | { |
266 | typedef typename tools::promote_args<T>::type result_type; |
267 | return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
268 | } |
269 | |
270 | #if 0 //def BOOST_MSVC |
271 | // |
272 | // We used to use ::_nextafter here, but doing so fails when using |
273 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own |
274 | // - albeit slower - code instead as at least that gives the correct answer. |
275 | // |
276 | template <class Policy> |
277 | inline double float_next(const double& val, const Policy& pol) |
278 | { |
279 | static const char* function = "float_next<%1%>(%1%)" ; |
280 | |
281 | if(!(boost::math::isfinite)(val) && (val > 0)) |
282 | return policies::raise_domain_error<double>( |
283 | function, |
284 | "Argument must be finite, but got %1%" , val, pol); |
285 | |
286 | if(val >= tools::max_value<double>()) |
287 | return policies::raise_overflow_error<double>(function, 0, pol); |
288 | |
289 | return ::_nextafter(val, tools::max_value<double>()); |
290 | } |
291 | #endif |
292 | |
293 | template <class T> |
294 | inline typename tools::promote_args<T>::type float_next(const T& val) |
295 | { |
296 | return float_next(val, policies::policy<>()); |
297 | } |
298 | |
299 | namespace detail{ |
300 | |
301 | template <class T, class Policy> |
302 | T float_prior_imp(const T& val, const boost::true_type&, const Policy& pol) |
303 | { |
304 | BOOST_MATH_STD_USING |
305 | int expon; |
306 | static const char* function = "float_prior<%1%>(%1%)" ; |
307 | |
308 | int fpclass = (boost::math::fpclassify)(val); |
309 | |
310 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
311 | { |
312 | if(val > 0) |
313 | return tools::max_value<T>(); |
314 | return policies::raise_domain_error<T>( |
315 | function, |
316 | "Argument must be finite, but got %1%" , val, pol); |
317 | } |
318 | |
319 | if(val <= -tools::max_value<T>()) |
320 | return -policies::raise_overflow_error<T>(function, 0, pol); |
321 | |
322 | if(val == 0) |
323 | return -detail::get_smallest_value<T>(); |
324 | |
325 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>())) |
326 | { |
327 | // |
328 | // Special case: if the value of the least significant bit is a denorm, and the result |
329 | // would not be a denorm, then shift the input, increment, and shift back. |
330 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
331 | // |
332 | return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); |
333 | } |
334 | |
335 | T remain = frexp(val, &expon); |
336 | if(remain == 0.5f) |
337 | --expon; // when val is a power of two we must reduce the exponent |
338 | T diff = ldexp(T(1), expon - tools::digits<T>()); |
339 | if(diff == 0) |
340 | diff = detail::get_smallest_value<T>(); |
341 | return val - diff; |
342 | } // float_prior_imp |
343 | // |
344 | // Special version for bases other than 2: |
345 | // |
346 | template <class T, class Policy> |
347 | T float_prior_imp(const T& val, const boost::false_type&, const Policy& pol) |
348 | { |
349 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
350 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
351 | |
352 | BOOST_MATH_STD_USING |
353 | boost::intmax_t expon; |
354 | static const char* function = "float_prior<%1%>(%1%)" ; |
355 | |
356 | int fpclass = (boost::math::fpclassify)(val); |
357 | |
358 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
359 | { |
360 | if(val > 0) |
361 | return tools::max_value<T>(); |
362 | return policies::raise_domain_error<T>( |
363 | function, |
364 | "Argument must be finite, but got %1%" , val, pol); |
365 | } |
366 | |
367 | if(val <= -tools::max_value<T>()) |
368 | return -policies::raise_overflow_error<T>(function, 0, pol); |
369 | |
370 | if(val == 0) |
371 | return -detail::get_smallest_value<T>(); |
372 | |
373 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>())) |
374 | { |
375 | // |
376 | // Special case: if the value of the least significant bit is a denorm, and the result |
377 | // would not be a denorm, then shift the input, increment, and shift back. |
378 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
379 | // |
380 | return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits); |
381 | } |
382 | |
383 | expon = 1 + ilogb(val); |
384 | T remain = scalbn(val, -expon); |
385 | if(remain * std::numeric_limits<T>::radix == 1) |
386 | --expon; // when val is a power of two we must reduce the exponent |
387 | T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits); |
388 | if(diff == 0) |
389 | diff = detail::get_smallest_value<T>(); |
390 | return val - diff; |
391 | } // float_prior_imp |
392 | |
393 | } // namespace detail |
394 | |
395 | template <class T, class Policy> |
396 | inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol) |
397 | { |
398 | typedef typename tools::promote_args<T>::type result_type; |
399 | return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
400 | } |
401 | |
402 | #if 0 //def BOOST_MSVC |
403 | // |
404 | // We used to use ::_nextafter here, but doing so fails when using |
405 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own |
406 | // - albeit slower - code instead as at least that gives the correct answer. |
407 | // |
408 | template <class Policy> |
409 | inline double float_prior(const double& val, const Policy& pol) |
410 | { |
411 | static const char* function = "float_prior<%1%>(%1%)" ; |
412 | |
413 | if(!(boost::math::isfinite)(val) && (val < 0)) |
414 | return policies::raise_domain_error<double>( |
415 | function, |
416 | "Argument must be finite, but got %1%" , val, pol); |
417 | |
418 | if(val <= -tools::max_value<double>()) |
419 | return -policies::raise_overflow_error<double>(function, 0, pol); |
420 | |
421 | return ::_nextafter(val, -tools::max_value<double>()); |
422 | } |
423 | #endif |
424 | |
425 | template <class T> |
426 | inline typename tools::promote_args<T>::type float_prior(const T& val) |
427 | { |
428 | return float_prior(val, policies::policy<>()); |
429 | } |
430 | |
431 | template <class T, class U, class Policy> |
432 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol) |
433 | { |
434 | typedef typename tools::promote_args<T, U>::type result_type; |
435 | return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol); |
436 | } |
437 | |
438 | template <class T, class U> |
439 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction) |
440 | { |
441 | return nextafter(val, direction, policies::policy<>()); |
442 | } |
443 | |
444 | namespace detail{ |
445 | |
446 | template <class T, class Policy> |
447 | T float_distance_imp(const T& a, const T& b, const boost::true_type&, const Policy& pol) |
448 | { |
449 | BOOST_MATH_STD_USING |
450 | // |
451 | // Error handling: |
452 | // |
453 | static const char* function = "float_distance<%1%>(%1%, %1%)" ; |
454 | if(!(boost::math::isfinite)(a)) |
455 | return policies::raise_domain_error<T>( |
456 | function, |
457 | "Argument a must be finite, but got %1%" , a, pol); |
458 | if(!(boost::math::isfinite)(b)) |
459 | return policies::raise_domain_error<T>( |
460 | function, |
461 | "Argument b must be finite, but got %1%" , b, pol); |
462 | // |
463 | // Special cases: |
464 | // |
465 | if(a > b) |
466 | return -float_distance(b, a, pol); |
467 | if(a == b) |
468 | return T(0); |
469 | if(a == 0) |
470 | return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)); |
471 | if(b == 0) |
472 | return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); |
473 | if(boost::math::sign(a) != boost::math::sign(b)) |
474 | return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)) |
475 | + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); |
476 | // |
477 | // By the time we get here, both a and b must have the same sign, we want |
478 | // b > a and both positive for the following logic: |
479 | // |
480 | if(a < 0) |
481 | return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol); |
482 | |
483 | BOOST_ASSERT(a >= 0); |
484 | BOOST_ASSERT(b >= a); |
485 | |
486 | int expon; |
487 | // |
488 | // Note that if a is a denorm then the usual formula fails |
489 | // because we actually have fewer than tools::digits<T>() |
490 | // significant bits in the representation: |
491 | // |
492 | (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon); |
493 | T upper = ldexp(T(1), expon); |
494 | T result = T(0); |
495 | // |
496 | // If b is greater than upper, then we *must* split the calculation |
497 | // as the size of the ULP changes with each order of magnitude change: |
498 | // |
499 | if(b > upper) |
500 | { |
501 | int expon2; |
502 | (void)frexp(b, &expon2); |
503 | T upper2 = ldexp(T(0.5), expon2); |
504 | result = float_distance(upper2, b); |
505 | result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1); |
506 | } |
507 | // |
508 | // Use compensated double-double addition to avoid rounding |
509 | // errors in the subtraction: |
510 | // |
511 | expon = tools::digits<T>() - expon; |
512 | T mb, x, y, z; |
513 | if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) |
514 | { |
515 | // |
516 | // Special case - either one end of the range is a denormal, or else the difference is. |
517 | // The regular code will fail if we're using the SSE2 registers on Intel and either |
518 | // the FTZ or DAZ flags are set. |
519 | // |
520 | T a2 = ldexp(a, tools::digits<T>()); |
521 | T b2 = ldexp(b, tools::digits<T>()); |
522 | mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2); |
523 | x = a2 + mb; |
524 | z = x - a2; |
525 | y = (a2 - (x - z)) + (mb - z); |
526 | |
527 | expon -= tools::digits<T>(); |
528 | } |
529 | else |
530 | { |
531 | mb = -(std::min)(upper, b); |
532 | x = a + mb; |
533 | z = x - a; |
534 | y = (a - (x - z)) + (mb - z); |
535 | } |
536 | if(x < 0) |
537 | { |
538 | x = -x; |
539 | y = -y; |
540 | } |
541 | result += ldexp(x, expon) + ldexp(y, expon); |
542 | // |
543 | // Result must be an integer: |
544 | // |
545 | BOOST_ASSERT(result == floor(result)); |
546 | return result; |
547 | } // float_distance_imp |
548 | // |
549 | // Special versions for bases other than 2: |
550 | // |
551 | template <class T, class Policy> |
552 | T float_distance_imp(const T& a, const T& b, const boost::false_type&, const Policy& pol) |
553 | { |
554 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
555 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
556 | |
557 | BOOST_MATH_STD_USING |
558 | // |
559 | // Error handling: |
560 | // |
561 | static const char* function = "float_distance<%1%>(%1%, %1%)" ; |
562 | if(!(boost::math::isfinite)(a)) |
563 | return policies::raise_domain_error<T>( |
564 | function, |
565 | "Argument a must be finite, but got %1%" , a, pol); |
566 | if(!(boost::math::isfinite)(b)) |
567 | return policies::raise_domain_error<T>( |
568 | function, |
569 | "Argument b must be finite, but got %1%" , b, pol); |
570 | // |
571 | // Special cases: |
572 | // |
573 | if(a > b) |
574 | return -float_distance(b, a, pol); |
575 | if(a == b) |
576 | return T(0); |
577 | if(a == 0) |
578 | return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)); |
579 | if(b == 0) |
580 | return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); |
581 | if(boost::math::sign(a) != boost::math::sign(b)) |
582 | return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)) |
583 | + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); |
584 | // |
585 | // By the time we get here, both a and b must have the same sign, we want |
586 | // b > a and both positive for the following logic: |
587 | // |
588 | if(a < 0) |
589 | return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol); |
590 | |
591 | BOOST_ASSERT(a >= 0); |
592 | BOOST_ASSERT(b >= a); |
593 | |
594 | boost::intmax_t expon; |
595 | // |
596 | // Note that if a is a denorm then the usual formula fails |
597 | // because we actually have fewer than tools::digits<T>() |
598 | // significant bits in the representation: |
599 | // |
600 | expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a); |
601 | T upper = scalbn(T(1), expon); |
602 | T result = T(0); |
603 | // |
604 | // If b is greater than upper, then we *must* split the calculation |
605 | // as the size of the ULP changes with each order of magnitude change: |
606 | // |
607 | if(b > upper) |
608 | { |
609 | boost::intmax_t expon2 = 1 + ilogb(b); |
610 | T upper2 = scalbn(T(1), expon2 - 1); |
611 | result = float_distance(upper2, b); |
612 | result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1); |
613 | } |
614 | // |
615 | // Use compensated double-double addition to avoid rounding |
616 | // errors in the subtraction: |
617 | // |
618 | expon = std::numeric_limits<T>::digits - expon; |
619 | T mb, x, y, z; |
620 | if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) |
621 | { |
622 | // |
623 | // Special case - either one end of the range is a denormal, or else the difference is. |
624 | // The regular code will fail if we're using the SSE2 registers on Intel and either |
625 | // the FTZ or DAZ flags are set. |
626 | // |
627 | T a2 = scalbn(a, std::numeric_limits<T>::digits); |
628 | T b2 = scalbn(b, std::numeric_limits<T>::digits); |
629 | mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2); |
630 | x = a2 + mb; |
631 | z = x - a2; |
632 | y = (a2 - (x - z)) + (mb - z); |
633 | |
634 | expon -= std::numeric_limits<T>::digits; |
635 | } |
636 | else |
637 | { |
638 | mb = -(std::min)(upper, b); |
639 | x = a + mb; |
640 | z = x - a; |
641 | y = (a - (x - z)) + (mb - z); |
642 | } |
643 | if(x < 0) |
644 | { |
645 | x = -x; |
646 | y = -y; |
647 | } |
648 | result += scalbn(x, expon) + scalbn(y, expon); |
649 | // |
650 | // Result must be an integer: |
651 | // |
652 | BOOST_ASSERT(result == floor(result)); |
653 | return result; |
654 | } // float_distance_imp |
655 | |
656 | } // namespace detail |
657 | |
658 | template <class T, class U, class Policy> |
659 | inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol) |
660 | { |
661 | typedef typename tools::promote_args<T, U>::type result_type; |
662 | return detail::float_distance_imp(detail::normalize_value(static_cast<result_type>(a), typename detail::has_hidden_guard_digits<result_type>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<result_type>::type()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
663 | } |
664 | |
665 | template <class T, class U> |
666 | typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b) |
667 | { |
668 | return boost::math::float_distance(a, b, policies::policy<>()); |
669 | } |
670 | |
671 | namespace detail{ |
672 | |
673 | template <class T, class Policy> |
674 | T float_advance_imp(T val, int distance, const boost::true_type&, const Policy& pol) |
675 | { |
676 | BOOST_MATH_STD_USING |
677 | // |
678 | // Error handling: |
679 | // |
680 | static const char* function = "float_advance<%1%>(%1%, int)" ; |
681 | |
682 | int fpclass = (boost::math::fpclassify)(val); |
683 | |
684 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
685 | return policies::raise_domain_error<T>( |
686 | function, |
687 | "Argument val must be finite, but got %1%" , val, pol); |
688 | |
689 | if(val < 0) |
690 | return -float_advance(-val, -distance, pol); |
691 | if(distance == 0) |
692 | return val; |
693 | if(distance == 1) |
694 | return float_next(val, pol); |
695 | if(distance == -1) |
696 | return float_prior(val, pol); |
697 | |
698 | if(fabs(val) < detail::get_min_shift_value<T>()) |
699 | { |
700 | // |
701 | // Special case: if the value of the least significant bit is a denorm, |
702 | // implement in terms of float_next/float_prior. |
703 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
704 | // |
705 | if(distance > 0) |
706 | { |
707 | do{ val = float_next(val, pol); } while(--distance); |
708 | } |
709 | else |
710 | { |
711 | do{ val = float_prior(val, pol); } while(++distance); |
712 | } |
713 | return val; |
714 | } |
715 | |
716 | int expon; |
717 | (void)frexp(val, &expon); |
718 | T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon); |
719 | if(val <= tools::min_value<T>()) |
720 | { |
721 | limit = sign(T(distance)) * tools::min_value<T>(); |
722 | } |
723 | T limit_distance = float_distance(val, limit); |
724 | while(fabs(limit_distance) < abs(distance)) |
725 | { |
726 | distance -= itrunc(limit_distance); |
727 | val = limit; |
728 | if(distance < 0) |
729 | { |
730 | limit /= 2; |
731 | expon--; |
732 | } |
733 | else |
734 | { |
735 | limit *= 2; |
736 | expon++; |
737 | } |
738 | limit_distance = float_distance(val, limit); |
739 | if(distance && (limit_distance == 0)) |
740 | { |
741 | return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode." , val, pol); |
742 | } |
743 | } |
744 | if((0.5f == frexp(val, &expon)) && (distance < 0)) |
745 | --expon; |
746 | T diff = 0; |
747 | if(val != 0) |
748 | diff = distance * ldexp(T(1), expon - tools::digits<T>()); |
749 | if(diff == 0) |
750 | diff = distance * detail::get_smallest_value<T>(); |
751 | return val += diff; |
752 | } // float_advance_imp |
753 | // |
754 | // Special version for bases other than 2: |
755 | // |
756 | template <class T, class Policy> |
757 | T float_advance_imp(T val, int distance, const boost::false_type&, const Policy& pol) |
758 | { |
759 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized); |
760 | BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2); |
761 | |
762 | BOOST_MATH_STD_USING |
763 | // |
764 | // Error handling: |
765 | // |
766 | static const char* function = "float_advance<%1%>(%1%, int)" ; |
767 | |
768 | int fpclass = (boost::math::fpclassify)(val); |
769 | |
770 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) |
771 | return policies::raise_domain_error<T>( |
772 | function, |
773 | "Argument val must be finite, but got %1%" , val, pol); |
774 | |
775 | if(val < 0) |
776 | return -float_advance(-val, -distance, pol); |
777 | if(distance == 0) |
778 | return val; |
779 | if(distance == 1) |
780 | return float_next(val, pol); |
781 | if(distance == -1) |
782 | return float_prior(val, pol); |
783 | |
784 | if(fabs(val) < detail::get_min_shift_value<T>()) |
785 | { |
786 | // |
787 | // Special case: if the value of the least significant bit is a denorm, |
788 | // implement in terms of float_next/float_prior. |
789 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. |
790 | // |
791 | if(distance > 0) |
792 | { |
793 | do{ val = float_next(val, pol); } while(--distance); |
794 | } |
795 | else |
796 | { |
797 | do{ val = float_prior(val, pol); } while(++distance); |
798 | } |
799 | return val; |
800 | } |
801 | |
802 | boost::intmax_t expon = 1 + ilogb(val); |
803 | T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon); |
804 | if(val <= tools::min_value<T>()) |
805 | { |
806 | limit = sign(T(distance)) * tools::min_value<T>(); |
807 | } |
808 | T limit_distance = float_distance(val, limit); |
809 | while(fabs(limit_distance) < abs(distance)) |
810 | { |
811 | distance -= itrunc(limit_distance); |
812 | val = limit; |
813 | if(distance < 0) |
814 | { |
815 | limit /= std::numeric_limits<T>::radix; |
816 | expon--; |
817 | } |
818 | else |
819 | { |
820 | limit *= std::numeric_limits<T>::radix; |
821 | expon++; |
822 | } |
823 | limit_distance = float_distance(val, limit); |
824 | if(distance && (limit_distance == 0)) |
825 | { |
826 | return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode." , val, pol); |
827 | } |
828 | } |
829 | /*expon = 1 + ilogb(val); |
830 | if((1 == scalbn(val, 1 + expon)) && (distance < 0)) |
831 | --expon;*/ |
832 | T diff = 0; |
833 | if(val != 0) |
834 | diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits); |
835 | if(diff == 0) |
836 | diff = distance * detail::get_smallest_value<T>(); |
837 | return val += diff; |
838 | } // float_advance_imp |
839 | |
840 | } // namespace detail |
841 | |
842 | template <class T, class Policy> |
843 | inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol) |
844 | { |
845 | typedef typename tools::promote_args<T>::type result_type; |
846 | return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
847 | } |
848 | |
849 | template <class T> |
850 | inline typename tools::promote_args<T>::type float_advance(const T& val, int distance) |
851 | { |
852 | return boost::math::float_advance(val, distance, policies::policy<>()); |
853 | } |
854 | |
855 | }} // boost math namespaces |
856 | |
857 | #endif // BOOST_MATH_SPECIAL_NEXT_HPP |
858 | |
859 | |