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
19#include <float.h>
20
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
25#endif
26#endif
27
28namespace boost{ namespace math{
29
30 namespace concepts {
31
32 class real_concept;
33 class std_real_concept;
34
35 }
36
37namespace detail{
38
39template <class T>
40struct has_hidden_guard_digits;
41template <>
42struct has_hidden_guard_digits<float> : public boost::false_type {};
43template <>
44struct has_hidden_guard_digits<double> : public boost::false_type {};
45template <>
46struct has_hidden_guard_digits<long double> : public boost::false_type {};
47#ifdef BOOST_HAS_FLOAT128
48template <>
49struct has_hidden_guard_digits<__float128> : public boost::false_type {};
50#endif
51template <>
52struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public boost::false_type {};
53template <>
54struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public boost::false_type {};
55
56template <class T, bool b>
57struct has_hidden_guard_digits_10 : public boost::false_type {};
58template <class T>
59struct has_hidden_guard_digits_10<T, true> : public boost::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {};
60
61template <class T>
62struct 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) >
66{};
67
68template <class T>
69inline const T& normalize_value(const T& val, const boost::false_type&) { return val; }
70template <class T>
71inline T normalize_value(const T& val, const boost::true_type&)
72{
73 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
74 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
75
76 boost::intmax_t shift = (boost::intmax_t)std::numeric_limits<T>::digits - (boost::intmax_t)ilogb(val) - 1;
77 T result = scalbn(val, shift);
78 result = round(result);
79 return scalbn(result, -shift);
80}
81
82template <class T>
83inline T get_smallest_value(boost::true_type const&)
84{
85 //
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.
89 //
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;;
93#else
94 return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m;
95#endif
96}
97
98template <class T>
99inline T get_smallest_value(boost::false_type const&)
100{
101 return tools::min_value<T>();
102}
103
104template <class T>
105inline T get_smallest_value()
106{
107#if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
108 return get_smallest_value<T>(boost::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
109#else
110 return get_smallest_value<T>(boost::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
111#endif
112}
113
114//
115// Returns the smallest value that won't generate denorms when
116// we calculate the value of the least-significant-bit:
117//
118template <class T>
119T get_min_shift_value();
120
121template <class T>
122struct min_shift_initializer
123{
124 struct init
125 {
126 init()
127 {
128 do_init();
129 }
130 static void do_init()
131 {
132 get_min_shift_value<T>();
133 }
134 void force_instantiate()const{}
135 };
136 static const init initializer;
137 static void force_instantiate()
138 {
139 initializer.force_instantiate();
140 }
141};
142
143template <class T>
144const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
145
146template <class T>
147inline T calc_min_shifted(const boost::true_type&)
148{
149 BOOST_MATH_STD_USING
150 return ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
151}
152template <class T>
153inline T calc_min_shifted(const boost::false_type&)
154{
155 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
156 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
157
158 return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1);
159}
160
161
162template <class T>
163inline T get_min_shift_value()
164{
165 static const T val = calc_min_shifted<T>(boost::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>());
166 min_shift_initializer<T>::force_instantiate();
167
168 return val;
169}
170
171template <class T, class Policy>
172T float_next_imp(const T& val, const boost::true_type&, const Policy& pol)
173{
174 BOOST_MATH_STD_USING
175 int expon;
176 static const char* function = "float_next<%1%>(%1%)";
177
178 int fpclass = (boost::math::fpclassify)(val);
179
180 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
181 {
182 if(val < 0)
183 return -tools::max_value<T>();
184 return policies::raise_domain_error<T>(
185 function,
186 "Argument must be finite, but got %1%", val, pol);
187 }
188
189 if(val >= tools::max_value<T>())
190 return policies::raise_overflow_error<T>(function, 0, pol);
191
192 if(val == 0)
193 return detail::get_smallest_value<T>();
194
195 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
196 {
197 //
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.
201 //
202 return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
203 }
204
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>());
208 if(diff == 0)
209 diff = detail::get_smallest_value<T>();
210 return val + diff;
211} // float_next_imp
212//
213// Special version for some base other than 2:
214//
215template <class T, class Policy>
216T float_next_imp(const T& val, const boost::false_type&, const Policy& pol)
217{
218 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
219 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
220
221 BOOST_MATH_STD_USING
222 boost::intmax_t expon;
223 static const char* function = "float_next<%1%>(%1%)";
224
225 int fpclass = (boost::math::fpclassify)(val);
226
227 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
228 {
229 if(val < 0)
230 return -tools::max_value<T>();
231 return policies::raise_domain_error<T>(
232 function,
233 "Argument must be finite, but got %1%", val, pol);
234 }
235
236 if(val >= tools::max_value<T>())
237 return policies::raise_overflow_error<T>(function, 0, pol);
238
239 if(val == 0)
240 return detail::get_smallest_value<T>();
241
242 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
243 {
244 //
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.
248 //
249 return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
250 }
251
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);
256 if(diff == 0)
257 diff = detail::get_smallest_value<T>();
258 return val + diff;
259} // float_next_imp
260
261} // namespace detail
262
263template <class T, class Policy>
264inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
265{
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()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
268}
269
270#if 0 //def BOOST_MSVC
271//
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.
275//
276template <class Policy>
277inline double float_next(const double& val, const Policy& pol)
278{
279 static const char* function = "float_next<%1%>(%1%)";
280
281 if(!(boost::math::isfinite)(val) && (val > 0))
282 return policies::raise_domain_error<double>(
283 function,
284 "Argument must be finite, but got %1%", val, pol);
285
286 if(val >= tools::max_value<double>())
287 return policies::raise_overflow_error<double>(function, 0, pol);
288
289 return ::_nextafter(val, tools::max_value<double>());
290}
291#endif
292
293template <class T>
294inline typename tools::promote_args<T>::type float_next(const T& val)
295{
296 return float_next(val, policies::policy<>());
297}
298
299namespace detail{
300
301template <class T, class Policy>
302T float_prior_imp(const T& val, const boost::true_type&, const Policy& pol)
303{
304 BOOST_MATH_STD_USING
305 int expon;
306 static const char* function = "float_prior<%1%>(%1%)";
307
308 int fpclass = (boost::math::fpclassify)(val);
309
310 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
311 {
312 if(val > 0)
313 return tools::max_value<T>();
314 return policies::raise_domain_error<T>(
315 function,
316 "Argument must be finite, but got %1%", val, pol);
317 }
318
319 if(val <= -tools::max_value<T>())
320 return -policies::raise_overflow_error<T>(function, 0, pol);
321
322 if(val == 0)
323 return -detail::get_smallest_value<T>();
324
325 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
326 {
327 //
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.
331 //
332 return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
333 }
334
335 T remain = frexp(val, &expon);
336 if(remain == 0.5f)
337 --expon; // when val is a power of two we must reduce the exponent
338 T diff = ldexp(T(1), expon - tools::digits<T>());
339 if(diff == 0)
340 diff = detail::get_smallest_value<T>();
341 return val - diff;
342} // float_prior_imp
343//
344// Special version for bases other than 2:
345//
346template <class T, class Policy>
347T float_prior_imp(const T& val, const boost::false_type&, const Policy& pol)
348{
349 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
350 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
351
352 BOOST_MATH_STD_USING
353 boost::intmax_t expon;
354 static const char* function = "float_prior<%1%>(%1%)";
355
356 int fpclass = (boost::math::fpclassify)(val);
357
358 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
359 {
360 if(val > 0)
361 return tools::max_value<T>();
362 return policies::raise_domain_error<T>(
363 function,
364 "Argument must be finite, but got %1%", val, pol);
365 }
366
367 if(val <= -tools::max_value<T>())
368 return -policies::raise_overflow_error<T>(function, 0, pol);
369
370 if(val == 0)
371 return -detail::get_smallest_value<T>();
372
373 if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
374 {
375 //
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.
379 //
380 return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits);
381 }
382
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);
388 if(diff == 0)
389 diff = detail::get_smallest_value<T>();
390 return val - diff;
391} // float_prior_imp
392
393} // namespace detail
394
395template <class T, class Policy>
396inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
397{
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()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
400}
401
402#if 0 //def BOOST_MSVC
403//
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.
407//
408template <class Policy>
409inline double float_prior(const double& val, const Policy& pol)
410{
411 static const char* function = "float_prior<%1%>(%1%)";
412
413 if(!(boost::math::isfinite)(val) && (val < 0))
414 return policies::raise_domain_error<double>(
415 function,
416 "Argument must be finite, but got %1%", val, pol);
417
418 if(val <= -tools::max_value<double>())
419 return -policies::raise_overflow_error<double>(function, 0, pol);
420
421 return ::_nextafter(val, -tools::max_value<double>());
422}
423#endif
424
425template <class T>
426inline typename tools::promote_args<T>::type float_prior(const T& val)
427{
428 return float_prior(val, policies::policy<>());
429}
430
431template <class T, class U, class Policy>
432inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
433{
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);
436}
437
438template <class T, class U>
439inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
440{
441 return nextafter(val, direction, policies::policy<>());
442}
443
444namespace detail{
445
446template <class T, class Policy>
447T float_distance_imp(const T& a, const T& b, const boost::true_type&, const Policy& pol)
448{
449 BOOST_MATH_STD_USING
450 //
451 // Error handling:
452 //
453 static const char* function = "float_distance<%1%>(%1%, %1%)";
454 if(!(boost::math::isfinite)(a))
455 return policies::raise_domain_error<T>(
456 function,
457 "Argument a must be finite, but got %1%", a, pol);
458 if(!(boost::math::isfinite)(b))
459 return policies::raise_domain_error<T>(
460 function,
461 "Argument b must be finite, but got %1%", b, pol);
462 //
463 // Special cases:
464 //
465 if(a > b)
466 return -float_distance(b, a, pol);
467 if(a == b)
468 return T(0);
469 if(a == 0)
470 return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
471 if(b == 0)
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));
476 //
477 // By the time we get here, both a and b must have the same sign, we want
478 // b > a and both positive for the following logic:
479 //
480 if(a < 0)
481 return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
482
483 BOOST_ASSERT(a >= 0);
484 BOOST_ASSERT(b >= a);
485
486 int expon;
487 //
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:
491 //
492 (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
493 T upper = ldexp(T(1), expon);
494 T result = T(0);
495 //
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:
498 //
499 if(b > upper)
500 {
501 int expon2;
502 (void)frexp(b, &expon2);
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);
506 }
507 //
508 // Use compensated double-double addition to avoid rounding
509 // errors in the subtraction:
510 //
511 expon = tools::digits<T>() - expon;
512 T mb, x, y, z;
513 if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
514 {
515 //
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.
519 //
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);
523 x = a2 + mb;
524 z = x - a2;
525 y = (a2 - (x - z)) + (mb - z);
526
527 expon -= tools::digits<T>();
528 }
529 else
530 {
531 mb = -(std::min)(upper, b);
532 x = a + mb;
533 z = x - a;
534 y = (a - (x - z)) + (mb - z);
535 }
536 if(x < 0)
537 {
538 x = -x;
539 y = -y;
540 }
541 result += ldexp(x, expon) + ldexp(y, expon);
542 //
543 // Result must be an integer:
544 //
545 BOOST_ASSERT(result == floor(result));
546 return result;
547} // float_distance_imp
548//
549// Special versions for bases other than 2:
550//
551template <class T, class Policy>
552T float_distance_imp(const T& a, const T& b, const boost::false_type&, const Policy& pol)
553{
554 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
555 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
556
557 BOOST_MATH_STD_USING
558 //
559 // Error handling:
560 //
561 static const char* function = "float_distance<%1%>(%1%, %1%)";
562 if(!(boost::math::isfinite)(a))
563 return policies::raise_domain_error<T>(
564 function,
565 "Argument a must be finite, but got %1%", a, pol);
566 if(!(boost::math::isfinite)(b))
567 return policies::raise_domain_error<T>(
568 function,
569 "Argument b must be finite, but got %1%", b, pol);
570 //
571 // Special cases:
572 //
573 if(a > b)
574 return -float_distance(b, a, pol);
575 if(a == b)
576 return T(0);
577 if(a == 0)
578 return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
579 if(b == 0)
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));
584 //
585 // By the time we get here, both a and b must have the same sign, we want
586 // b > a and both positive for the following logic:
587 //
588 if(a < 0)
589 return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
590
591 BOOST_ASSERT(a >= 0);
592 BOOST_ASSERT(b >= a);
593
594 boost::intmax_t expon;
595 //
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:
599 //
600 expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a);
601 T upper = scalbn(T(1), expon);
602 T result = T(0);
603 //
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:
606 //
607 if(b > upper)
608 {
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);
613 }
614 //
615 // Use compensated double-double addition to avoid rounding
616 // errors in the subtraction:
617 //
618 expon = std::numeric_limits<T>::digits - expon;
619 T mb, x, y, z;
620 if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
621 {
622 //
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.
626 //
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);
630 x = a2 + mb;
631 z = x - a2;
632 y = (a2 - (x - z)) + (mb - z);
633
634 expon -= std::numeric_limits<T>::digits;
635 }
636 else
637 {
638 mb = -(std::min)(upper, b);
639 x = a + mb;
640 z = x - a;
641 y = (a - (x - z)) + (mb - z);
642 }
643 if(x < 0)
644 {
645 x = -x;
646 y = -y;
647 }
648 result += scalbn(x, expon) + scalbn(y, expon);
649 //
650 // Result must be an integer:
651 //
652 BOOST_ASSERT(result == floor(result));
653 return result;
654} // float_distance_imp
655
656} // namespace detail
657
658template <class T, class U, class Policy>
659inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
660{
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()), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
663}
664
665template <class T, class U>
666typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
667{
668 return boost::math::float_distance(a, b, policies::policy<>());
669}
670
671namespace detail{
672
673template <class T, class Policy>
674T float_advance_imp(T val, int distance, const boost::true_type&, const Policy& pol)
675{
676 BOOST_MATH_STD_USING
677 //
678 // Error handling:
679 //
680 static const char* function = "float_advance<%1%>(%1%, int)";
681
682 int fpclass = (boost::math::fpclassify)(val);
683
684 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
685 return policies::raise_domain_error<T>(
686 function,
687 "Argument val must be finite, but got %1%", val, pol);
688
689 if(val < 0)
690 return -float_advance(-val, -distance, pol);
691 if(distance == 0)
692 return val;
693 if(distance == 1)
694 return float_next(val, pol);
695 if(distance == -1)
696 return float_prior(val, pol);
697
698 if(fabs(val) < detail::get_min_shift_value<T>())
699 {
700 //
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.
704 //
705 if(distance > 0)
706 {
707 do{ val = float_next(val, pol); } while(--distance);
708 }
709 else
710 {
711 do{ val = float_prior(val, pol); } while(++distance);
712 }
713 return val;
714 }
715
716 int expon;
717 (void)frexp(val, &expon);
718 T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
719 if(val <= tools::min_value<T>())
720 {
721 limit = sign(T(distance)) * tools::min_value<T>();
722 }
723 T limit_distance = float_distance(val, limit);
724 while(fabs(limit_distance) < abs(distance))
725 {
726 distance -= itrunc(limit_distance);
727 val = limit;
728 if(distance < 0)
729 {
730 limit /= 2;
731 expon--;
732 }
733 else
734 {
735 limit *= 2;
736 expon++;
737 }
738 limit_distance = float_distance(val, limit);
739 if(distance && (limit_distance == 0))
740 {
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);
742 }
743 }
744 if((0.5f == frexp(val, &expon)) && (distance < 0))
745 --expon;
746 T diff = 0;
747 if(val != 0)
748 diff = distance * ldexp(T(1), expon - tools::digits<T>());
749 if(diff == 0)
750 diff = distance * detail::get_smallest_value<T>();
751 return val += diff;
752} // float_advance_imp
753//
754// Special version for bases other than 2:
755//
756template <class T, class Policy>
757T float_advance_imp(T val, int distance, const boost::false_type&, const Policy& pol)
758{
759 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
760 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
761
762 BOOST_MATH_STD_USING
763 //
764 // Error handling:
765 //
766 static const char* function = "float_advance<%1%>(%1%, int)";
767
768 int fpclass = (boost::math::fpclassify)(val);
769
770 if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE))
771 return policies::raise_domain_error<T>(
772 function,
773 "Argument val must be finite, but got %1%", val, pol);
774
775 if(val < 0)
776 return -float_advance(-val, -distance, pol);
777 if(distance == 0)
778 return val;
779 if(distance == 1)
780 return float_next(val, pol);
781 if(distance == -1)
782 return float_prior(val, pol);
783
784 if(fabs(val) < detail::get_min_shift_value<T>())
785 {
786 //
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.
790 //
791 if(distance > 0)
792 {
793 do{ val = float_next(val, pol); } while(--distance);
794 }
795 else
796 {
797 do{ val = float_prior(val, pol); } while(++distance);
798 }
799 return val;
800 }
801
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>())
805 {
806 limit = sign(T(distance)) * tools::min_value<T>();
807 }
808 T limit_distance = float_distance(val, limit);
809 while(fabs(limit_distance) < abs(distance))
810 {
811 distance -= itrunc(limit_distance);
812 val = limit;
813 if(distance < 0)
814 {
815 limit /= std::numeric_limits<T>::radix;
816 expon--;
817 }
818 else
819 {
820 limit *= std::numeric_limits<T>::radix;
821 expon++;
822 }
823 limit_distance = float_distance(val, limit);
824 if(distance && (limit_distance == 0))
825 {
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);
827 }
828 }
829 /*expon = 1 + ilogb(val);
830 if((1 == scalbn(val, 1 + expon)) && (distance < 0))
831 --expon;*/
832 T diff = 0;
833 if(val != 0)
834 diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits);
835 if(diff == 0)
836 diff = distance * detail::get_smallest_value<T>();
837 return val += diff;
838} // float_advance_imp
839
840} // namespace detail
841
842template <class T, class Policy>
843inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
844{
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, boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
847}
848
849template <class T>
850inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
851{
852 return boost::math::float_advance(val, distance, policies::policy<>());
853}
854
855}} // boost math namespaces
856
857#endif // BOOST_MATH_SPECIAL_NEXT_HPP
858
859

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