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