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 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
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 return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized>());
110}
111
112template <class T>
113inline 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//
121template <class T>
122T get_min_shift_value();
123
124template <class T>
125struct 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
146template <class T>
147const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
148
149template <class T>
150inline 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}
155template <class T>
156inline 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
165template <class T>
166inline 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
174template <class T, bool b = boost::math::tools::detail::has_backend_type<T>::value>
175struct exponent_type
176{
177 typedef int type;
178};
179
180template <class T>
181struct exponent_type<T, true>
182{
183 typedef typename T::backend_type::exponent_type type;
184};
185
186template <class T, class Policy>
187T 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//
232template <class T, class Policy>
233T 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
282template <class T, class Policy>
283inline 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//
295template <class Policy>
296inline 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
312template <class T>
313inline typename tools::promote_args<T>::type float_next(const T& val)
314{
315 return float_next(val, policies::policy<>());
316}
317
318namespace detail{
319
320template <class T, class Policy>
321T 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//
367template <class T, class Policy>
368T 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
418template <class T, class Policy>
419inline 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//
431template <class Policy>
432inline 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
448template <class T>
449inline typename tools::promote_args<T>::type float_prior(const T& val)
450{
451 return float_prior(val, policies::policy<>());
452}
453
454template <class T, class U, class Policy>
455inline 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
461template <class T, class U>
462inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
463{
464 return nextafter(val, direction, policies::policy<>());
465}
466
467namespace detail{
468
469template <class T, class Policy>
470T 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//
574template <class T, class Policy>
575T 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
681template <class T, class U, class Policy>
682inline 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
715template <class T, class U>
716typename 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
721namespace detail{
722
723template <class T, class Policy>
724T 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//
806template <class T, class Policy>
807T 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
892template <class T, class Policy>
893inline 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
899template <class T>
900inline 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

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