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