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)
6 #ifndef BOOST_MATH_SPECIAL_NEXT_HPP
7 #define BOOST_MATH_SPECIAL_NEXT_HPP
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>
21 #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3)))
22 #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__)
23 #include "xmmintrin.h"
24 #define BOOST_MATH_CHECK_SSE2
28 namespace boost{ namespace math{
33 struct std_real_concept;
40 struct has_hidden_guard_digits;
42 struct has_hidden_guard_digits<float> : public mpl::false_ {};
44 struct has_hidden_guard_digits<double> : public mpl::false_ {};
46 struct has_hidden_guard_digits<long double> : public mpl::false_ {};
47 #ifdef BOOST_HAS_FLOAT128
49 struct has_hidden_guard_digits<__float128> : public mpl::false_ {};
52 struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public mpl::false_ {};
54 struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public mpl::false_ {};
56 template <class T, bool b>
57 struct has_hidden_guard_digits_10 : public mpl::false_ {};
59 struct has_hidden_guard_digits_10<T, true> : public mpl::bool_<(std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
62 struct has_hidden_guard_digits
63 : public has_hidden_guard_digits_10<T,
64 std::numeric_limits<T>::is_specialized
65 && (std::numeric_limits<T>::radix == 10) >
69 inline const T& normalize_value(const T& val, const mpl::false_&) { return val; }
71 inline T normalize_value(const T& val, const mpl::true_&)
73 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
74 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
76 boost::intmax_t shift = std::numeric_limits<T>::digits - ilogb(val) - 1;
77 T result = scalbn(val, shift);
78 result = round(result);
79 return scalbn(result, -shift);
83 inline T get_smallest_value(mpl::true_ const&)
86 // numeric_limits lies about denorms being present - particularly
87 // when this can be turned on or off at runtime, as is the case
88 // when using the SSE2 registers in DAZ or FTZ mode.
90 static const T m = std::numeric_limits<T>::denorm_min();
91 #ifdef BOOST_MATH_CHECK_SSE2
92 return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;;
94 return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
99 inline T get_smallest_value(mpl::false_ const&)
101 return tools::min_value<T>();
105 inline T get_smallest_value()
107 #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
108 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
110 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
115 // Returns the smallest value that won't generate denorms when
116 // we calculate the value of the least-significant-bit:
119 T get_min_shift_value();
122 struct min_shift_initializer
130 static void do_init()
132 get_min_shift_value<T>();
134 void force_instantiate()const{}
136 static const init initializer;
137 static void force_instantiate()
139 initializer.force_instantiate();
144 const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
147 inline T calc_min_shifted(const mpl::true_&)
150 return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
153 inline T calc_min_shifted(const mpl::false_&)
155 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
156 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
158 return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
163 inline T get_min_shift_value()
165 static const T val = calc_min_shifted<T>(mpl::bool_<!std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
166 min_shift_initializer<T>::force_instantiate();
171 template <class T, class Policy>
172 T float_next_imp(const T& val, const mpl::true_&, const Policy& pol)
176 static const char* function = "float_next<%1%>(%1%)";
178 int fpclass = (boost::math::fpclassify)(val);
180 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
183 return -tools::max_value<T>();
184 return policies::raise_domain_error<T>(
186 "Argument must be finite, but got %1%", val, pol);
189 if(val >= tools::max_value<T>())
190 return policies::raise_overflow_error<T>(function, 0, pol);
193 return detail::get_smallest_value<T>();
195 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
198 // Special case: if the value of the least significant bit is a denorm, and the result
199 // would not be a denorm, then shift the input, increment, and shift back.
200 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
202 return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
205 if(-0.5f == frexp(val, &expon))
206 --expon; // reduce exponent when val is a power of two, and negative.
207 T diff = ldexp(T(1), expon - tools::digits<T>());
209 diff = detail::get_smallest_value<T>();
213 // Special version for some base other than 2:
215 template <class T, class Policy>
216 T float_next_imp(const T& val, const mpl::false_&, const Policy& pol)
218 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
219 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
222 boost::intmax_t expon;
223 static const char* function = "float_next<%1%>(%1%)";
225 int fpclass = (boost::math::fpclassify)(val);
227 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
230 return -tools::max_value<T>();
231 return policies::raise_domain_error<T>(
233 "Argument must be finite, but got %1%", val, pol);
236 if(val >= tools::max_value<T>())
237 return policies::raise_overflow_error<T>(function, 0, pol);
240 return detail::get_smallest_value<T>();
242 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
245 // Special case: if the value of the least significant bit is a denorm, and the result
246 // would not be a denorm, then shift the input, increment, and shift back.
247 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
249 return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
252 expon = 1 + ilogb(val);
253 if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix)
254 --expon; // reduce exponent when val is a power of base, and negative.
255 T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
257 diff = detail::get_smallest_value<T>();
261 } // namespace detail
263 template <class T, class Policy>
264 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
266 typedef typename tools::promote_args<T>::type result_type;
267 return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
270 #if 0 //def BOOST_MSVC
272 // We used to use ::_nextafter here, but doing so fails when using
273 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
274 // - albeit slower - code instead as at least that gives the correct answer.
276 template <class Policy>
277 inline double float_next(const double& val, const Policy& pol)
279 static const char* function = "float_next<%1%>(%1%)";
281 if(!(boost::math::isfinite)(val) && (val > 0))
282 return policies::raise_domain_error<double>(
284 "Argument must be finite, but got %1%", val, pol);
286 if(val >= tools::max_value<double>())
287 return policies::raise_overflow_error<double>(function, 0, pol);
289 return ::_nextafter(val, tools::max_value<double>());
294 inline typename tools::promote_args<T>::type float_next(const T& val)
296 return float_next(val, policies::policy<>());
301 template <class T, class Policy>
302 T float_prior_imp(const T& val, const mpl::true_&, const Policy& pol)
306 static const char* function = "float_prior<%1%>(%1%)";
308 int fpclass = (boost::math::fpclassify)(val);
310 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
313 return tools::max_value<T>();
314 return policies::raise_domain_error<T>(
316 "Argument must be finite, but got %1%", val, pol);
319 if(val <= -tools::max_value<T>())
320 return -policies::raise_overflow_error<T>(function, 0, pol);
323 return -detail::get_smallest_value<T>();
325 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
328 // Special case: if the value of the least significant bit is a denorm, and the result
329 // would not be a denorm, then shift the input, increment, and shift back.
330 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
332 return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
335 T remain = frexp(val, &expon);
337 --expon; // when val is a power of two we must reduce the exponent
338 T diff = ldexp(T(1), expon - tools::digits<T>());
340 diff = detail::get_smallest_value<T>();
344 // Special version for bases other than 2:
346 template <class T, class Policy>
347 T float_prior_imp(const T& val, const mpl::false_&, const Policy& pol)
349 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
350 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
353 boost::intmax_t expon;
354 static const char* function = "float_prior<%1%>(%1%)";
356 int fpclass = (boost::math::fpclassify)(val);
358 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
361 return tools::max_value<T>();
362 return policies::raise_domain_error<T>(
364 "Argument must be finite, but got %1%", val, pol);
367 if(val <= -tools::max_value<T>())
368 return -policies::raise_overflow_error<T>(function, 0, pol);
371 return -detail::get_smallest_value<T>();
373 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
376 // Special case: if the value of the least significant bit is a denorm, and the result
377 // would not be a denorm, then shift the input, increment, and shift back.
378 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
380 return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
383 expon = 1 + ilogb(val);
384 T remain = scalbn(val, -expon);
385 if(remain * std::numeric_limits<T>::radix == 1)
386 --expon; // when val is a power of two we must reduce the exponent
387 T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
389 diff = detail::get_smallest_value<T>();
393 } // namespace detail
395 template <class T, class Policy>
396 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
398 typedef typename tools::promote_args<T>::type result_type;
399 return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
402 #if 0 //def BOOST_MSVC
404 // We used to use ::_nextafter here, but doing so fails when using
405 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
406 // - albeit slower - code instead as at least that gives the correct answer.
408 template <class Policy>
409 inline double float_prior(const double& val, const Policy& pol)
411 static const char* function = "float_prior<%1%>(%1%)";
413 if(!(boost::math::isfinite)(val) && (val < 0))
414 return policies::raise_domain_error<double>(
416 "Argument must be finite, but got %1%", val, pol);
418 if(val <= -tools::max_value<double>())
419 return -policies::raise_overflow_error<double>(function, 0, pol);
421 return ::_nextafter(val, -tools::max_value<double>());
426 inline typename tools::promote_args<T>::type float_prior(const T& val)
428 return float_prior(val, policies::policy<>());
431 template <class T, class U, class Policy>
432 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
434 typedef typename tools::promote_args<T, U>::type result_type;
435 return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
438 template <class T, class U>
439 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
441 return nextafter(val, direction, policies::policy<>());
446 template <class T, class Policy>
447 T float_distance_imp(const T& a, const T& b, const mpl::true_&, const Policy& pol)
453 static const char* function = "float_distance<%1%>(%1%, %1%)";
454 if(!(boost::math::isfinite)(a))
455 return policies::raise_domain_error<T>(
457 "Argument a must be finite, but got %1%", a, pol);
458 if(!(boost::math::isfinite)(b))
459 return policies::raise_domain_error<T>(
461 "Argument b must be finite, but got %1%", b, pol);
466 return -float_distance(b, a, pol);
470 return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
472 return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
473 if(boost::math::sign(a) != boost::math::sign(b))
474 return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
475 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
477 // By the time we get here, both a and b must have the same sign, we want
478 // b > a and both postive for the following logic:
481 return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
483 BOOST_ASSERT(a >= 0);
484 BOOST_ASSERT(b >= a);
488 // Note that if a is a denorm then the usual formula fails
489 // because we actually have fewer than tools::digits<T>()
490 // significant bits in the representation:
492 frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
493 T upper = ldexp(T(1), expon);
496 // If b is greater than upper, then we *must* split the calculation
497 // as the size of the ULP changes with each order of magnitude change:
503 T upper2 = ldexp(T(0.5), expon2);
504 result = float_distance(upper2, b);
505 result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1);
508 // Use compensated double-double addition to avoid rounding
509 // errors in the subtraction:
511 expon = tools::digits<T>() - expon;
513 if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
516 // Special case - either one end of the range is a denormal, or else the difference is.
517 // The regular code will fail if we're using the SSE2 registers on Intel and either
518 // the FTZ or DAZ flags are set.
520 T a2 = ldexp(a, tools::digits<T>());
521 T b2 = ldexp(b, tools::digits<T>());
522 mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
525 y = (a2 - (x - z)) + (mb - z);
527 expon -= tools::digits<T>();
531 mb = -(std::min)(upper, b);
534 y = (a - (x - z)) + (mb - z);
541 result += ldexp(x, expon) + ldexp(y, expon);
543 // Result must be an integer:
545 BOOST_ASSERT(result == floor(result));
547 } // float_distance_imp
549 // Special versions for bases other than 2:
551 template <class T, class Policy>
552 T float_distance_imp(const T& a, const T& b, const mpl::false_&, const Policy& pol)
554 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
555 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
561 static const char* function = "float_distance<%1%>(%1%, %1%)";
562 if(!(boost::math::isfinite)(a))
563 return policies::raise_domain_error<T>(
565 "Argument a must be finite, but got %1%", a, pol);
566 if(!(boost::math::isfinite)(b))
567 return policies::raise_domain_error<T>(
569 "Argument b must be finite, but got %1%", b, pol);
574 return -float_distance(b, a, pol);
578 return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
580 return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
581 if(boost::math::sign(a) != boost::math::sign(b))
582 return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
583 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
585 // By the time we get here, both a and b must have the same sign, we want
586 // b > a and both postive for the following logic:
589 return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
591 BOOST_ASSERT(a >= 0);
592 BOOST_ASSERT(b >= a);
594 boost::intmax_t expon;
596 // Note that if a is a denorm then the usual formula fails
597 // because we actually have fewer than tools::digits<T>()
598 // significant bits in the representation:
600 expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
601 T upper = scalbn(T(1), expon);
604 // If b is greater than upper, then we *must* split the calculation
605 // as the size of the ULP changes with each order of magnitude change:
609 boost::intmax_t expon2 = 1 + ilogb(b);
610 T upper2 = scalbn(T(1), expon2 - 1);
611 result = float_distance(upper2, b);
612 result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1);
615 // Use compensated double-double addition to avoid rounding
616 // errors in the subtraction:
618 expon = std::numeric_limits<T>::digits - expon;
620 if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
623 // Special case - either one end of the range is a denormal, or else the difference is.
624 // The regular code will fail if we're using the SSE2 registers on Intel and either
625 // the FTZ or DAZ flags are set.
627 T a2 = scalbn(a, std::numeric_limits<T>::digits);
628 T b2 = scalbn(b, std::numeric_limits<T>::digits);
629 mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2);
632 y = (a2 - (x - z)) + (mb - z);
634 expon -= std::numeric_limits<T>::digits;
638 mb = -(std::min)(upper, b);
641 y = (a - (x - z)) + (mb - z);
648 result += scalbn(x, expon) + scalbn(y, expon);
650 // Result must be an integer:
652 BOOST_ASSERT(result == floor(result));
654 } // float_distance_imp
656 } // namespace detail
658 template <class T, class U, class Policy>
659 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
661 typedef typename tools::promote_args<T, U>::type result_type;
662 return detail::float_distance_imp(detail::normalize_value(static_cast<result_type>(a), typename detail::has_hidden_guard_digits<result_type>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<result_type>::type()), mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
665 template <class T, class U>
666 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
668 return boost::math::float_distance(a, b, policies::policy<>());
673 template <class T, class Policy>
674 T float_advance_imp(T val, int distance, const mpl::true_&, const Policy& pol)
680 static const char* function = "float_advance<%1%>(%1%, int)";
682 int fpclass = (boost::math::fpclassify)(val);
684 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
685 return policies::raise_domain_error<T>(
687 "Argument val must be finite, but got %1%", val, pol);
690 return -float_advance(-val, -distance, pol);
694 return float_next(val, pol);
696 return float_prior(val, pol);
698 if(fabs(val) < detail::get_min_shift_value<T>())
701 // Special case: if the value of the least significant bit is a denorm,
702 // implement in terms of float_next/float_prior.
703 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
707 do{ val = float_next(val, pol); } while(--distance);
711 do{ val = float_prior(val, pol); } while(++distance);
718 T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
719 if(val <= tools::min_value<T>())
721 limit = sign(T(distance)) * tools::min_value<T>();
723 T limit_distance = float_distance(val, limit);
724 while(fabs(limit_distance) < abs(distance))
726 distance -= itrunc(limit_distance);
738 limit_distance = float_distance(val, limit);
739 if(distance && (limit_distance == 0))
741 return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol);
744 if((0.5f == frexp(val, &expon)) && (distance < 0))
748 diff = distance * ldexp(T(1), expon - tools::digits<T>());
750 diff = distance * detail::get_smallest_value<T>();
752 } // float_advance_imp
754 // Special version for bases other than 2:
756 template <class T, class Policy>
757 T float_advance_imp(T val, int distance, const mpl::false_&, const Policy& pol)
759 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
760 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
766 static const char* function = "float_advance<%1%>(%1%, int)";
768 int fpclass = (boost::math::fpclassify)(val);
770 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
771 return policies::raise_domain_error<T>(
773 "Argument val must be finite, but got %1%", val, pol);
776 return -float_advance(-val, -distance, pol);
780 return float_next(val, pol);
782 return float_prior(val, pol);
784 if(fabs(val) < detail::get_min_shift_value<T>())
787 // Special case: if the value of the least significant bit is a denorm,
788 // implement in terms of float_next/float_prior.
789 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
793 do{ val = float_next(val, pol); } while(--distance);
797 do{ val = float_prior(val, pol); } while(++distance);
802 boost::intmax_t expon = 1 + ilogb(val);
803 T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon);
804 if(val <= tools::min_value<T>())
806 limit = sign(T(distance)) * tools::min_value<T>();
808 T limit_distance = float_distance(val, limit);
809 while(fabs(limit_distance) < abs(distance))
811 distance -= itrunc(limit_distance);
815 limit /= std::numeric_limits<T>::radix;
820 limit *= std::numeric_limits<T>::radix;
823 limit_distance = float_distance(val, limit);
824 if(distance && (limit_distance == 0))
826 return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol);
829 /*expon = 1 + ilogb(val);
830 if((1 == scalbn(val, 1 + expon)) && (distance < 0))
834 diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
836 diff = distance * detail::get_smallest_value<T>();
838 } // float_advance_imp
840 } // namespace detail
842 template <class T, class Policy>
843 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
845 typedef typename tools::promote_args<T>::type result_type;
846 return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, mpl::bool_<!std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
850 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
852 return boost::math::float_advance(val, distance, policies::policy<>());
855 }} // boost math namespaces
857 #endif // BOOST_MATH_SPECIAL_NEXT_HPP