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/policies/error_handling.hpp>
14 #include <boost/math/special_functions/fpclassify.hpp>
15 #include <boost/math/special_functions/sign.hpp>
16 #include <boost/math/special_functions/trunc.hpp>
22 namespace boost{ namespace math{
27 inline T get_smallest_value(mpl::true_ const&)
30 // numeric_limits lies about denorms being present - particularly
31 // when this can be turned on or off at runtime, as is the case
32 // when using the SSE2 registers in DAZ or FTZ mode.
34 static const T m = std::numeric_limits<T>::denorm_min();
35 return ((tools::min_value<T>() - m) == tools::min_value<T>()) ? tools::min_value<T>() : m;
39 inline T get_smallest_value(mpl::false_ const&)
41 return tools::min_value<T>();
45 inline T get_smallest_value()
47 #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
48 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
50 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
55 // Returns the smallest value that won't generate denorms when
56 // we calculate the value of the least-significant-bit:
59 T get_min_shift_value();
62 struct min_shift_initializer
72 get_min_shift_value<T>();
74 void force_instantiate()const{}
76 static const init initializer;
77 static void force_instantiate()
79 initializer.force_instantiate();
84 const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer;
88 inline T get_min_shift_value()
91 static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1);
92 min_shift_initializer<T>::force_instantiate();
97 template <class T, class Policy>
98 T float_next_imp(const T& val, const Policy& pol)
102 static const char* function = "float_next<%1%>(%1%)";
104 int fpclass = (boost::math::fpclassify)(val);
106 if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
109 return -tools::max_value<T>();
110 return policies::raise_domain_error<T>(
112 "Argument must be finite, but got %1%", val, pol);
115 if(val >= tools::max_value<T>())
116 return policies::raise_overflow_error<T>(function, 0, pol);
119 return detail::get_smallest_value<T>();
121 if((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>()))
124 // Special case: if the value of the least significant bit is a denorm, and the result
125 // would not be a denorm, then shift the input, increment, and shift back.
126 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
128 return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
131 if(-0.5f == frexp(val, &expon))
132 --expon; // reduce exponent when val is a power of two, and negative.
133 T diff = ldexp(T(1), expon - tools::digits<T>());
135 diff = detail::get_smallest_value<T>();
141 template <class T, class Policy>
142 inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol)
144 typedef typename tools::promote_args<T>::type result_type;
145 return detail::float_next_imp(static_cast<result_type>(val), pol);
148 #if 0 //def BOOST_MSVC
150 // We used to use ::_nextafter here, but doing so fails when using
151 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
152 // - albeit slower - code instead as at least that gives the correct answer.
154 template <class Policy>
155 inline double float_next(const double& val, const Policy& pol)
157 static const char* function = "float_next<%1%>(%1%)";
159 if(!(boost::math::isfinite)(val) && (val > 0))
160 return policies::raise_domain_error<double>(
162 "Argument must be finite, but got %1%", val, pol);
164 if(val >= tools::max_value<double>())
165 return policies::raise_overflow_error<double>(function, 0, pol);
167 return ::_nextafter(val, tools::max_value<double>());
172 inline typename tools::promote_args<T>::type float_next(const T& val)
174 return float_next(val, policies::policy<>());
179 template <class T, class Policy>
180 T float_prior_imp(const T& val, const Policy& pol)
184 static const char* function = "float_prior<%1%>(%1%)";
186 int fpclass = (boost::math::fpclassify)(val);
188 if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
191 return tools::max_value<T>();
192 return policies::raise_domain_error<T>(
194 "Argument must be finite, but got %1%", val, pol);
197 if(val <= -tools::max_value<T>())
198 return -policies::raise_overflow_error<T>(function, 0, pol);
201 return -detail::get_smallest_value<T>();
203 if((fpclass != FP_SUBNORMAL) && (fpclass != FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>()))
206 // Special case: if the value of the least significant bit is a denorm, and the result
207 // would not be a denorm, then shift the input, increment, and shift back.
208 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
210 return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>());
213 T remain = frexp(val, &expon);
215 --expon; // when val is a power of two we must reduce the exponent
216 T diff = ldexp(T(1), expon - tools::digits<T>());
218 diff = detail::get_smallest_value<T>();
224 template <class T, class Policy>
225 inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol)
227 typedef typename tools::promote_args<T>::type result_type;
228 return detail::float_prior_imp(static_cast<result_type>(val), pol);
231 #if 0 //def BOOST_MSVC
233 // We used to use ::_nextafter here, but doing so fails when using
234 // the SSE2 registers if the FTZ or DAZ flags are set, so use our own
235 // - albeit slower - code instead as at least that gives the correct answer.
237 template <class Policy>
238 inline double float_prior(const double& val, const Policy& pol)
240 static const char* function = "float_prior<%1%>(%1%)";
242 if(!(boost::math::isfinite)(val) && (val < 0))
243 return policies::raise_domain_error<double>(
245 "Argument must be finite, but got %1%", val, pol);
247 if(val <= -tools::max_value<double>())
248 return -policies::raise_overflow_error<double>(function, 0, pol);
250 return ::_nextafter(val, -tools::max_value<double>());
255 inline typename tools::promote_args<T>::type float_prior(const T& val)
257 return float_prior(val, policies::policy<>());
260 template <class T, class U, class Policy>
261 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol)
263 typedef typename tools::promote_args<T, U>::type result_type;
264 return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol);
267 template <class T, class U>
268 inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction)
270 return nextafter(val, direction, policies::policy<>());
275 template <class T, class Policy>
276 T float_distance_imp(const T& a, const T& b, const Policy& pol)
282 static const char* function = "float_distance<%1%>(%1%, %1%)";
283 if(!(boost::math::isfinite)(a))
284 return policies::raise_domain_error<T>(
286 "Argument a must be finite, but got %1%", a, pol);
287 if(!(boost::math::isfinite)(b))
288 return policies::raise_domain_error<T>(
290 "Argument b must be finite, but got %1%", b, pol);
295 return -float_distance(b, a, pol);
299 return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol));
301 return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
302 if(boost::math::sign(a) != boost::math::sign(b))
303 return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol))
304 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol));
306 // By the time we get here, both a and b must have the same sign, we want
307 // b > a and both postive for the following logic:
310 return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol);
312 BOOST_ASSERT(a >= 0);
313 BOOST_ASSERT(b >= a);
317 // Note that if a is a denorm then the usual formula fails
318 // because we actually have fewer than tools::digits<T>()
319 // significant bits in the representation:
321 frexp(((boost::math::fpclassify)(a) == FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
322 T upper = ldexp(T(1), expon);
324 expon = tools::digits<T>() - expon;
326 // If b is greater than upper, then we *must* split the calculation
327 // as the size of the ULP changes with each order of magnitude change:
331 result = float_distance(upper, b);
334 // Use compensated double-double addition to avoid rounding
335 // errors in the subtraction:
338 if(((boost::math::fpclassify)(a) == FP_SUBNORMAL) || (b - a < tools::min_value<T>()))
341 // Special case - either one end of the range is a denormal, or else the difference is.
342 // The regular code will fail if we're using the SSE2 registers on Intel and either
343 // the FTZ or DAZ flags are set.
345 T a2 = ldexp(a, tools::digits<T>());
346 T b2 = ldexp(b, tools::digits<T>());
347 mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2);
350 y = (a2 - (x - z)) + (mb - z);
352 expon -= tools::digits<T>();
356 mb = -(std::min)(upper, b);
359 y = (a - (x - z)) + (mb - z);
366 result += ldexp(x, expon) + ldexp(y, expon);
368 // Result must be an integer:
370 BOOST_ASSERT(result == floor(result));
376 template <class T, class U, class Policy>
377 inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol)
379 typedef typename tools::promote_args<T, U>::type result_type;
380 return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol);
383 template <class T, class U>
384 typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b)
386 return boost::math::float_distance(a, b, policies::policy<>());
391 template <class T, class Policy>
392 T float_advance_imp(T val, int distance, const Policy& pol)
398 static const char* function = "float_advance<%1%>(%1%, int)";
400 int fpclass = (boost::math::fpclassify)(val);
402 if((fpclass == FP_NAN) || (fpclass == FP_INFINITE))
403 return policies::raise_domain_error<T>(
405 "Argument val must be finite, but got %1%", val, pol);
408 return -float_advance(-val, -distance, pol);
412 return float_next(val, pol);
414 return float_prior(val, pol);
416 if(fabs(val) < detail::get_min_shift_value<T>())
419 // Special case: if the value of the least significant bit is a denorm,
420 // implement in terms of float_next/float_prior.
421 // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set.
425 do{ val = float_next(val, pol); } while(--distance);
429 do{ val = float_prior(val, pol); } while(++distance);
436 T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
437 if(val <= tools::min_value<T>())
439 limit = sign(T(distance)) * tools::min_value<T>();
441 T limit_distance = float_distance(val, limit);
442 while(fabs(limit_distance) < abs(distance))
444 distance -= itrunc(limit_distance);
456 limit_distance = float_distance(val, limit);
457 if(distance && (limit_distance == 0))
459 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);
462 if((0.5f == frexp(val, &expon)) && (distance < 0))
466 diff = distance * ldexp(T(1), expon - tools::digits<T>());
468 diff = distance * detail::get_smallest_value<T>();
474 template <class T, class Policy>
475 inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol)
477 typedef typename tools::promote_args<T>::type result_type;
478 return detail::float_advance_imp(static_cast<result_type>(val), distance, pol);
482 inline typename tools::promote_args<T>::type float_advance(const T& val, int distance)
484 return boost::math::float_advance(val, distance, policies::policy<>());
489 #endif // BOOST_MATH_SPECIAL_NEXT_HPP