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
30namespace boost{ namespace math{
31
32 namespace concepts {
33
34 class real_concept;
35 class std_real_concept;
36
37 }
38
39namespace detail{
40
41template <class T>
42struct has_hidden_guard_digits;
43template <>
44struct has_hidden_guard_digits<float> : public std::false_type {};
45template <>
46struct has_hidden_guard_digits<double> : public std::false_type {};
47template <>
48struct has_hidden_guard_digits<long double> : public std::false_type {};
49#ifdef BOOST_HAS_FLOAT128
50template <>
51struct has_hidden_guard_digits<__float128> : public std::false_type {};
52#endif
53template <>
54struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public std::false_type {};
55template <>
56struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public std::false_type {};
57
58template <class T, bool b>
59struct has_hidden_guard_digits_10 : public std::false_type {};
60template <class T>
61struct has_hidden_guard_digits_10<T, true> : public std::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
62
63template <class T>
64struct 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
70template <class T>
71inline const T& normalize_value(const T& val, const std::false_type&) { return val; }
72template <class T>
73inline 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
84template <class T>
85inline T get_smallest_value(std::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
100template <class T>
101inline T get_smallest_value(std::false_type const&)
102{
103 return tools::min_value<T>();
104}
105
106template <class T>
107inline T get_smallest_value()
108{
109#if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
110 return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
111#else
112 return get_smallest_value<T>(std::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//
120template <class T>
121T get_min_shift_value();
122
123template <class T>
124struct min_shift_initializer
125{
126 struct init
127 {
128 init()
129 {
130 do_init();
131 }
132 static void do_init()
133 {
134 get_min_shift_value<T>();
135 }
136 void force_instantiate()const{}
137 };
138 static const init initializer;
139 static void force_instantiate()
140 {
141 initializer.force_instantiate();
142 }
143};
144
145template <class T>
146const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
147
148template <class T>
149inline T calc_min_shifted(const std::true_type&)
150{
151 BOOST_MATH_STD_USING
152 return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
153}
154template <class T>
155inline T calc_min_shifted(const std::false_type&)
156{
157 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
158 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
159
160 return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
161}
162
163
164template <class T>
165inline T get_min_shift_value()
166{
167 static const T val = calc_min_shifted<T>(std::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
173template <class T, bool b = boost::math::tools::detail::has_backend_type<T>::value>
174struct exponent_type
175{
176 typedef int type;
177};
178
179template <class T>
180struct exponent_type<T, true>
181{
182 typedef typename T::backend_type::exponent_type type;
183};
184
185template <class T, class Policy>
186T float_next_imp(const T& val, const std::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, nullptr, 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//
231template <class T, class Policy>
232T float_next_imp(const T& val, const std::false_type&, const Policy& pol)
233{
234 typedef typename exponent_type<T>::type exponent_type;
235
236 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
237 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
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, nullptr, 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
281template <class T, class Policy>
282inline 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()), std::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//
294template <class Policy>
295inline 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, nullptr, pol);
306
307 return ::_nextafter(val, tools::max_value<double>());
308}
309#endif
310
311template <class T>
312inline typename tools::promote_args<T>::type float_next(const T& val)
313{
314 return float_next(val, policies::policy<>());
315}
316
317namespace detail{
318
319template <class T, class Policy>
320T float_prior_imp(const T& val, const std::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, nullptr, 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//
366template <class T, class Policy>
367T float_prior_imp(const T& val, const std::false_type&, const Policy& pol)
368{
369 typedef typename exponent_type<T>::type exponent_type;
370
371 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
372 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
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, nullptr, 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
417template <class T, class Policy>
418inline 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()), std::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//
430template <class Policy>
431inline 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, nullptr, pol);
442
443 return ::_nextafter(val, -tools::max_value<double>());
444}
445#endif
446
447template <class T>
448inline typename tools::promote_args<T>::type float_prior(const T& val)
449{
450 return float_prior(val, policies::policy<>());
451}
452
453template <class T, class U, class Policy>
454inline 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
460template <class T, class U>
461inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
462{
463 return nextafter(val, direction, policies::policy<>());
464}
465
466namespace detail{
467
468template <class T, class Policy>
469T float_distance_imp(const T& a, const T& b, const std::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_MATH_ASSERT(a >= 0);
506 BOOST_MATH_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_MATH_ASSERT(result == floor(result));
568 return result;
569} // float_distance_imp
570//
571// Special versions for bases other than 2:
572//
573template <class T, class Policy>
574T float_distance_imp(const T& a, const T& b, const std::false_type&, const Policy& pol)
575{
576 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
577 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
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_MATH_ASSERT(a >= 0);
614 BOOST_MATH_ASSERT(b >= a);
615
616 std::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 std::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_MATH_ASSERT(result == floor(result));
675 return result;
676} // float_distance_imp
677
678} // namespace detail
679
680template <class T, class U, class Policy>
681inline 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 static_assert(
687 (std::is_same<T, U>::value
688 || (std::is_integral<T>::value && !std::is_integral<U>::value)
689 || (!std::is_integral<T>::value && std::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 (!std::is_same<T, U>::value)
697 {
698 BOOST_IF_CONSTEXPR(std::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()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
711 }
712}
713
714template <class T, class U>
715typename 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
720namespace detail{
721
722template <class T, class Policy>
723T float_advance_imp(T val, int distance, const std::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//
805template <class T, class Policy>
806T float_advance_imp(T val, int distance, const std::false_type&, const Policy& pol)
807{
808 static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
809 static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized.");
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 std::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
891template <class T, class Policy>
892inline 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, std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
896}
897
898template <class T>
899inline 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

source code of include/boost/math/special_functions/next.hpp