]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | |
1e59de90 | 12 | |
7c673cae FG |
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> | |
20effc67 | 18 | #include <boost/math/tools/traits.hpp> |
1e59de90 TL |
19 | #include <type_traits> |
20 | #include <cfloat> | |
7c673cae | 21 | |
7c673cae FG |
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 | ||
30 | namespace boost{ namespace math{ | |
31 | ||
b32b8144 FG |
32 | namespace concepts { |
33 | ||
34 | class real_concept; | |
11fdf7f2 | 35 | class std_real_concept; |
b32b8144 FG |
36 | |
37 | } | |
38 | ||
7c673cae FG |
39 | namespace detail{ |
40 | ||
b32b8144 FG |
41 | template <class T> |
42 | struct has_hidden_guard_digits; | |
43 | template <> | |
1e59de90 | 44 | struct has_hidden_guard_digits<float> : public std::false_type {}; |
b32b8144 | 45 | template <> |
1e59de90 | 46 | struct has_hidden_guard_digits<double> : public std::false_type {}; |
b32b8144 | 47 | template <> |
1e59de90 | 48 | struct has_hidden_guard_digits<long double> : public std::false_type {}; |
b32b8144 FG |
49 | #ifdef BOOST_HAS_FLOAT128 |
50 | template <> | |
1e59de90 | 51 | struct has_hidden_guard_digits<__float128> : public std::false_type {}; |
b32b8144 FG |
52 | #endif |
53 | template <> | |
1e59de90 | 54 | struct has_hidden_guard_digits<boost::math::concepts::real_concept> : public std::false_type {}; |
b32b8144 | 55 | template <> |
1e59de90 | 56 | struct has_hidden_guard_digits<boost::math::concepts::std_real_concept> : public std::false_type {}; |
b32b8144 FG |
57 | |
58 | template <class T, bool b> | |
1e59de90 | 59 | struct has_hidden_guard_digits_10 : public std::false_type {}; |
b32b8144 | 60 | template <class T> |
1e59de90 | 61 | struct has_hidden_guard_digits_10<T, true> : public std::integral_constant<bool, (std::numeric_limits<T>::digits10 != std::numeric_limits<T>::max_digits10)> {}; |
b32b8144 FG |
62 | |
63 | template <class T> | |
64 | struct 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 | ||
70 | template <class T> | |
1e59de90 | 71 | inline const T& normalize_value(const T& val, const std::false_type&) { return val; } |
b32b8144 | 72 | template <class T> |
1e59de90 | 73 | inline T normalize_value(const T& val, const std::true_type&) |
b32b8144 | 74 | { |
1e59de90 TL |
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."); | |
b32b8144 | 77 | |
1e59de90 | 78 | std::intmax_t shift = (std::intmax_t)std::numeric_limits<T>::digits - (std::intmax_t)ilogb(val) - 1; |
b32b8144 FG |
79 | T result = scalbn(val, shift); |
80 | result = round(result); | |
81 | return scalbn(result, -shift); | |
82 | } | |
83 | ||
7c673cae | 84 | template <class T> |
1e59de90 | 85 | inline T get_smallest_value(std::true_type const&) |
7c673cae FG |
86 | { |
87 | // | |
88 | // numeric_limits lies about denorms being present - particularly | |
89 | // when this can be turned on or off at runtime, as is the case | |
90 | // when using the SSE2 registers in DAZ or FTZ mode. | |
91 | // | |
92 | static const T m = std::numeric_limits<T>::denorm_min(); | |
93 | #ifdef BOOST_MATH_CHECK_SSE2 | |
20effc67 | 94 | return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m; |
7c673cae FG |
95 | #else |
96 | return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m; | |
97 | #endif | |
98 | } | |
99 | ||
100 | template <class T> | |
1e59de90 | 101 | inline T get_smallest_value(std::false_type const&) |
7c673cae FG |
102 | { |
103 | return tools::min_value<T>(); | |
104 | } | |
105 | ||
106 | template <class T> | |
107 | inline T get_smallest_value() | |
108 | { | |
109 | #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310) | |
1e59de90 | 110 | return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>()); |
7c673cae | 111 | #else |
1e59de90 | 112 | return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>()); |
7c673cae FG |
113 | #endif |
114 | } | |
115 | ||
116 | // | |
117 | // Returns the smallest value that won't generate denorms when | |
118 | // we calculate the value of the least-significant-bit: | |
119 | // | |
120 | template <class T> | |
121 | T get_min_shift_value(); | |
122 | ||
123 | template <class T> | |
124 | struct min_shift_initializer | |
125 | { | |
126 | struct init | |
127 | { | |
128 | init() | |
129 | { | |
130 | do_init(); | |
131 | } | |
132 | static void do_init() | |
133 | { | |
134 | get_min_shift_value<T>(); | |
135 | } | |
136 | void force_instantiate()const{} | |
137 | }; | |
138 | static const init initializer; | |
139 | static void force_instantiate() | |
140 | { | |
141 | initializer.force_instantiate(); | |
142 | } | |
143 | }; | |
144 | ||
145 | template <class T> | |
146 | const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer; | |
147 | ||
b32b8144 | 148 | template <class T> |
1e59de90 | 149 | inline T calc_min_shifted(const std::true_type&) |
b32b8144 FG |
150 | { |
151 | BOOST_MATH_STD_USING | |
152 | return ldexp(tools::min_value<T>(), tools::digits<T>() + 1); | |
153 | } | |
154 | template <class T> | |
1e59de90 | 155 | inline T calc_min_shifted(const std::false_type&) |
b32b8144 | 156 | { |
1e59de90 TL |
157 | static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized."); |
158 | static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized."); | |
b32b8144 FG |
159 | |
160 | return scalbn(tools::min_value<T>(), std::numeric_limits<T>::digits + 1); | |
161 | } | |
162 | ||
7c673cae FG |
163 | |
164 | template <class T> | |
165 | inline T get_min_shift_value() | |
166 | { | |
1e59de90 | 167 | static const T val = calc_min_shifted<T>(std::integral_constant<bool, !std::numeric_limits<T>::is_specialized || std::numeric_limits<T>::radix == 2>()); |
7c673cae FG |
168 | min_shift_initializer<T>::force_instantiate(); |
169 | ||
170 | return val; | |
171 | } | |
172 | ||
20effc67 TL |
173 | template <class T, bool b = boost::math::tools::detail::has_backend_type<T>::value> |
174 | struct exponent_type | |
175 | { | |
176 | typedef int type; | |
177 | }; | |
178 | ||
179 | template <class T> | |
180 | struct exponent_type<T, true> | |
181 | { | |
182 | typedef typename T::backend_type::exponent_type type; | |
183 | }; | |
184 | ||
7c673cae | 185 | template <class T, class Policy> |
1e59de90 | 186 | T float_next_imp(const T& val, const std::true_type&, const Policy& pol) |
7c673cae | 187 | { |
20effc67 TL |
188 | typedef typename exponent_type<T>::type exponent_type; |
189 | ||
7c673cae | 190 | BOOST_MATH_STD_USING |
20effc67 | 191 | exponent_type expon; |
7c673cae FG |
192 | static const char* function = "float_next<%1%>(%1%)"; |
193 | ||
194 | int fpclass = (boost::math::fpclassify)(val); | |
195 | ||
196 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
197 | { | |
198 | if(val < 0) | |
199 | return -tools::max_value<T>(); | |
200 | return policies::raise_domain_error<T>( | |
201 | function, | |
202 | "Argument must be finite, but got %1%", val, pol); | |
203 | } | |
204 | ||
205 | if(val >= tools::max_value<T>()) | |
206 | return policies::raise_overflow_error<T>(function, 0, pol); | |
207 | ||
208 | if(val == 0) | |
209 | return detail::get_smallest_value<T>(); | |
210 | ||
211 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>())) | |
212 | { | |
213 | // | |
214 | // Special case: if the value of the least significant bit is a denorm, and the result | |
215 | // would not be a denorm, then shift the input, increment, and shift back. | |
216 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
217 | // | |
218 | return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); | |
219 | } | |
220 | ||
221 | if(-0.5f == frexp(val, &expon)) | |
222 | --expon; // reduce exponent when val is a power of two, and negative. | |
223 | T diff = ldexp(T(1), expon - tools::digits<T>()); | |
224 | if(diff == 0) | |
225 | diff = detail::get_smallest_value<T>(); | |
226 | return val + diff; | |
b32b8144 FG |
227 | } // float_next_imp |
228 | // | |
229 | // Special version for some base other than 2: | |
230 | // | |
231 | template <class T, class Policy> | |
1e59de90 | 232 | T float_next_imp(const T& val, const std::false_type&, const Policy& pol) |
b32b8144 | 233 | { |
20effc67 TL |
234 | typedef typename exponent_type<T>::type exponent_type; |
235 | ||
1e59de90 TL |
236 | static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized."); |
237 | static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized."); | |
7c673cae | 238 | |
b32b8144 | 239 | BOOST_MATH_STD_USING |
20effc67 | 240 | exponent_type expon; |
b32b8144 FG |
241 | static const char* function = "float_next<%1%>(%1%)"; |
242 | ||
243 | int fpclass = (boost::math::fpclassify)(val); | |
244 | ||
245 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
246 | { | |
247 | if(val < 0) | |
248 | return -tools::max_value<T>(); | |
249 | return policies::raise_domain_error<T>( | |
250 | function, | |
251 | "Argument must be finite, but got %1%", val, pol); | |
252 | } | |
253 | ||
254 | if(val >= tools::max_value<T>()) | |
255 | return policies::raise_overflow_error<T>(function, 0, pol); | |
256 | ||
257 | if(val == 0) | |
258 | return detail::get_smallest_value<T>(); | |
259 | ||
260 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>())) | |
261 | { | |
262 | // | |
263 | // Special case: if the value of the least significant bit is a denorm, and the result | |
264 | // would not be a denorm, then shift the input, increment, and shift back. | |
265 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
266 | // | |
267 | return scalbn(float_next(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits); | |
268 | } | |
269 | ||
270 | expon = 1 + ilogb(val); | |
271 | if(-1 == scalbn(val, -expon) * std::numeric_limits<T>::radix) | |
272 | --expon; // reduce exponent when val is a power of base, and negative. | |
273 | T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits); | |
274 | if(diff == 0) | |
275 | diff = detail::get_smallest_value<T>(); | |
276 | return val + diff; | |
277 | } // float_next_imp | |
278 | ||
279 | } // namespace detail | |
7c673cae FG |
280 | |
281 | template <class T, class Policy> | |
282 | inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol) | |
283 | { | |
284 | typedef typename tools::promote_args<T>::type result_type; | |
1e59de90 | 285 | return detail::float_next_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
7c673cae FG |
286 | } |
287 | ||
288 | #if 0 //def BOOST_MSVC | |
289 | // | |
290 | // We used to use ::_nextafter here, but doing so fails when using | |
291 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own | |
292 | // - albeit slower - code instead as at least that gives the correct answer. | |
293 | // | |
294 | template <class Policy> | |
295 | inline double float_next(const double& val, const Policy& pol) | |
296 | { | |
297 | static const char* function = "float_next<%1%>(%1%)"; | |
298 | ||
299 | if(!(boost::math::isfinite)(val) && (val > 0)) | |
300 | return policies::raise_domain_error<double>( | |
301 | function, | |
302 | "Argument must be finite, but got %1%", val, pol); | |
303 | ||
304 | if(val >= tools::max_value<double>()) | |
305 | return policies::raise_overflow_error<double>(function, 0, pol); | |
306 | ||
307 | return ::_nextafter(val, tools::max_value<double>()); | |
308 | } | |
309 | #endif | |
310 | ||
311 | template <class T> | |
312 | inline typename tools::promote_args<T>::type float_next(const T& val) | |
313 | { | |
314 | return float_next(val, policies::policy<>()); | |
315 | } | |
316 | ||
317 | namespace detail{ | |
318 | ||
319 | template <class T, class Policy> | |
1e59de90 | 320 | T float_prior_imp(const T& val, const std::true_type&, const Policy& pol) |
7c673cae | 321 | { |
20effc67 TL |
322 | typedef typename exponent_type<T>::type exponent_type; |
323 | ||
7c673cae | 324 | BOOST_MATH_STD_USING |
20effc67 | 325 | exponent_type expon; |
7c673cae FG |
326 | static const char* function = "float_prior<%1%>(%1%)"; |
327 | ||
328 | int fpclass = (boost::math::fpclassify)(val); | |
329 | ||
330 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
331 | { | |
332 | if(val > 0) | |
333 | return tools::max_value<T>(); | |
334 | return policies::raise_domain_error<T>( | |
335 | function, | |
336 | "Argument must be finite, but got %1%", val, pol); | |
337 | } | |
338 | ||
339 | if(val <= -tools::max_value<T>()) | |
340 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
341 | ||
342 | if(val == 0) | |
343 | return -detail::get_smallest_value<T>(); | |
344 | ||
345 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>())) | |
346 | { | |
347 | // | |
348 | // Special case: if the value of the least significant bit is a denorm, and the result | |
349 | // would not be a denorm, then shift the input, increment, and shift back. | |
350 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
351 | // | |
352 | return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); | |
353 | } | |
354 | ||
355 | T remain = frexp(val, &expon); | |
b32b8144 | 356 | if(remain == 0.5f) |
7c673cae FG |
357 | --expon; // when val is a power of two we must reduce the exponent |
358 | T diff = ldexp(T(1), expon - tools::digits<T>()); | |
359 | if(diff == 0) | |
360 | diff = detail::get_smallest_value<T>(); | |
361 | return val - diff; | |
b32b8144 FG |
362 | } // float_prior_imp |
363 | // | |
364 | // Special version for bases other than 2: | |
365 | // | |
366 | template <class T, class Policy> | |
1e59de90 | 367 | T float_prior_imp(const T& val, const std::false_type&, const Policy& pol) |
b32b8144 | 368 | { |
20effc67 TL |
369 | typedef typename exponent_type<T>::type exponent_type; |
370 | ||
1e59de90 TL |
371 | static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized."); |
372 | static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized."); | |
7c673cae | 373 | |
b32b8144 | 374 | BOOST_MATH_STD_USING |
20effc67 | 375 | exponent_type expon; |
b32b8144 FG |
376 | static const char* function = "float_prior<%1%>(%1%)"; |
377 | ||
378 | int fpclass = (boost::math::fpclassify)(val); | |
379 | ||
380 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
381 | { | |
382 | if(val > 0) | |
383 | return tools::max_value<T>(); | |
384 | return policies::raise_domain_error<T>( | |
385 | function, | |
386 | "Argument must be finite, but got %1%", val, pol); | |
387 | } | |
388 | ||
389 | if(val <= -tools::max_value<T>()) | |
390 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
391 | ||
392 | if(val == 0) | |
393 | return -detail::get_smallest_value<T>(); | |
394 | ||
395 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>())) | |
396 | { | |
397 | // | |
398 | // Special case: if the value of the least significant bit is a denorm, and the result | |
399 | // would not be a denorm, then shift the input, increment, and shift back. | |
400 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
401 | // | |
402 | return scalbn(float_prior(T(scalbn(val, 2 * std::numeric_limits<T>::digits)), pol), -2 * std::numeric_limits<T>::digits); | |
403 | } | |
404 | ||
405 | expon = 1 + ilogb(val); | |
406 | T remain = scalbn(val, -expon); | |
407 | if(remain * std::numeric_limits<T>::radix == 1) | |
408 | --expon; // when val is a power of two we must reduce the exponent | |
409 | T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits); | |
410 | if(diff == 0) | |
411 | diff = detail::get_smallest_value<T>(); | |
412 | return val - diff; | |
413 | } // float_prior_imp | |
414 | ||
415 | } // namespace detail | |
7c673cae FG |
416 | |
417 | template <class T, class Policy> | |
418 | inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol) | |
419 | { | |
420 | typedef typename tools::promote_args<T>::type result_type; | |
1e59de90 | 421 | return detail::float_prior_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
7c673cae FG |
422 | } |
423 | ||
424 | #if 0 //def BOOST_MSVC | |
425 | // | |
426 | // We used to use ::_nextafter here, but doing so fails when using | |
427 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own | |
428 | // - albeit slower - code instead as at least that gives the correct answer. | |
429 | // | |
430 | template <class Policy> | |
431 | inline double float_prior(const double& val, const Policy& pol) | |
432 | { | |
433 | static const char* function = "float_prior<%1%>(%1%)"; | |
434 | ||
435 | if(!(boost::math::isfinite)(val) && (val < 0)) | |
436 | return policies::raise_domain_error<double>( | |
437 | function, | |
438 | "Argument must be finite, but got %1%", val, pol); | |
439 | ||
440 | if(val <= -tools::max_value<double>()) | |
441 | return -policies::raise_overflow_error<double>(function, 0, pol); | |
442 | ||
443 | return ::_nextafter(val, -tools::max_value<double>()); | |
444 | } | |
445 | #endif | |
446 | ||
447 | template <class T> | |
448 | inline typename tools::promote_args<T>::type float_prior(const T& val) | |
449 | { | |
450 | return float_prior(val, policies::policy<>()); | |
451 | } | |
452 | ||
453 | template <class T, class U, class Policy> | |
454 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol) | |
455 | { | |
456 | typedef typename tools::promote_args<T, U>::type result_type; | |
457 | return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol); | |
458 | } | |
459 | ||
460 | template <class T, class U> | |
461 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction) | |
462 | { | |
463 | return nextafter(val, direction, policies::policy<>()); | |
464 | } | |
465 | ||
466 | namespace detail{ | |
467 | ||
468 | template <class T, class Policy> | |
1e59de90 | 469 | T float_distance_imp(const T& a, const T& b, const std::true_type&, const Policy& pol) |
7c673cae FG |
470 | { |
471 | BOOST_MATH_STD_USING | |
472 | // | |
473 | // Error handling: | |
474 | // | |
475 | static const char* function = "float_distance<%1%>(%1%, %1%)"; | |
476 | if(!(boost::math::isfinite)(a)) | |
477 | return policies::raise_domain_error<T>( | |
478 | function, | |
479 | "Argument a must be finite, but got %1%", a, pol); | |
480 | if(!(boost::math::isfinite)(b)) | |
481 | return policies::raise_domain_error<T>( | |
482 | function, | |
483 | "Argument b must be finite, but got %1%", b, pol); | |
484 | // | |
485 | // Special cases: | |
486 | // | |
487 | if(a > b) | |
488 | return -float_distance(b, a, pol); | |
489 | if(a == b) | |
b32b8144 | 490 | return T(0); |
7c673cae FG |
491 | if(a == 0) |
492 | return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)); | |
493 | if(b == 0) | |
494 | return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
495 | if(boost::math::sign(a) != boost::math::sign(b)) | |
496 | return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)) | |
497 | + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
498 | // | |
499 | // By the time we get here, both a and b must have the same sign, we want | |
f67539c2 | 500 | // b > a and both positive for the following logic: |
7c673cae FG |
501 | // |
502 | if(a < 0) | |
503 | return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol); | |
504 | ||
1e59de90 TL |
505 | BOOST_MATH_ASSERT(a >= 0); |
506 | BOOST_MATH_ASSERT(b >= a); | |
7c673cae FG |
507 | |
508 | int expon; | |
509 | // | |
510 | // Note that if a is a denorm then the usual formula fails | |
511 | // because we actually have fewer than tools::digits<T>() | |
512 | // significant bits in the representation: | |
513 | // | |
92f5a8d4 | 514 | (void)frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon); |
7c673cae | 515 | T upper = ldexp(T(1), expon); |
b32b8144 | 516 | T result = T(0); |
7c673cae FG |
517 | // |
518 | // If b is greater than upper, then we *must* split the calculation | |
519 | // as the size of the ULP changes with each order of magnitude change: | |
520 | // | |
521 | if(b > upper) | |
522 | { | |
b32b8144 | 523 | int expon2; |
92f5a8d4 | 524 | (void)frexp(b, &expon2); |
b32b8144 FG |
525 | T upper2 = ldexp(T(0.5), expon2); |
526 | result = float_distance(upper2, b); | |
527 | result += (expon2 - expon - 1) * ldexp(T(1), tools::digits<T>() - 1); | |
7c673cae FG |
528 | } |
529 | // | |
530 | // Use compensated double-double addition to avoid rounding | |
531 | // errors in the subtraction: | |
532 | // | |
b32b8144 | 533 | expon = tools::digits<T>() - expon; |
7c673cae FG |
534 | T mb, x, y, z; |
535 | if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) | |
536 | { | |
537 | // | |
538 | // Special case - either one end of the range is a denormal, or else the difference is. | |
539 | // The regular code will fail if we're using the SSE2 registers on Intel and either | |
540 | // the FTZ or DAZ flags are set. | |
541 | // | |
542 | T a2 = ldexp(a, tools::digits<T>()); | |
543 | T b2 = ldexp(b, tools::digits<T>()); | |
544 | mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2); | |
545 | x = a2 + mb; | |
546 | z = x - a2; | |
547 | y = (a2 - (x - z)) + (mb - z); | |
548 | ||
549 | expon -= tools::digits<T>(); | |
550 | } | |
551 | else | |
552 | { | |
553 | mb = -(std::min)(upper, b); | |
554 | x = a + mb; | |
555 | z = x - a; | |
556 | y = (a - (x - z)) + (mb - z); | |
557 | } | |
558 | if(x < 0) | |
559 | { | |
560 | x = -x; | |
561 | y = -y; | |
562 | } | |
563 | result += ldexp(x, expon) + ldexp(y, expon); | |
564 | // | |
565 | // Result must be an integer: | |
566 | // | |
1e59de90 | 567 | BOOST_MATH_ASSERT(result == floor(result)); |
7c673cae | 568 | return result; |
b32b8144 FG |
569 | } // float_distance_imp |
570 | // | |
571 | // Special versions for bases other than 2: | |
572 | // | |
573 | template <class T, class Policy> | |
1e59de90 | 574 | T float_distance_imp(const T& a, const T& b, const std::false_type&, const Policy& pol) |
b32b8144 | 575 | { |
1e59de90 TL |
576 | static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized."); |
577 | static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized."); | |
7c673cae | 578 | |
b32b8144 FG |
579 | BOOST_MATH_STD_USING |
580 | // | |
581 | // Error handling: | |
582 | // | |
583 | static const char* function = "float_distance<%1%>(%1%, %1%)"; | |
584 | if(!(boost::math::isfinite)(a)) | |
585 | return policies::raise_domain_error<T>( | |
586 | function, | |
587 | "Argument a must be finite, but got %1%", a, pol); | |
588 | if(!(boost::math::isfinite)(b)) | |
589 | return policies::raise_domain_error<T>( | |
590 | function, | |
591 | "Argument b must be finite, but got %1%", b, pol); | |
592 | // | |
593 | // Special cases: | |
594 | // | |
595 | if(a > b) | |
596 | return -float_distance(b, a, pol); | |
597 | if(a == b) | |
598 | return T(0); | |
599 | if(a == 0) | |
600 | return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)); | |
601 | if(b == 0) | |
602 | return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
603 | if(boost::math::sign(a) != boost::math::sign(b)) | |
604 | return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)) | |
605 | + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
606 | // | |
607 | // By the time we get here, both a and b must have the same sign, we want | |
f67539c2 | 608 | // b > a and both positive for the following logic: |
b32b8144 FG |
609 | // |
610 | if(a < 0) | |
611 | return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol); | |
612 | ||
1e59de90 TL |
613 | BOOST_MATH_ASSERT(a >= 0); |
614 | BOOST_MATH_ASSERT(b >= a); | |
b32b8144 | 615 | |
1e59de90 | 616 | std::intmax_t expon; |
b32b8144 FG |
617 | // |
618 | // Note that if a is a denorm then the usual formula fails | |
619 | // because we actually have fewer than tools::digits<T>() | |
620 | // significant bits in the representation: | |
621 | // | |
622 | expon = 1 + ilogb(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a); | |
623 | T upper = scalbn(T(1), expon); | |
624 | T result = T(0); | |
625 | // | |
626 | // If b is greater than upper, then we *must* split the calculation | |
627 | // as the size of the ULP changes with each order of magnitude change: | |
628 | // | |
629 | if(b > upper) | |
630 | { | |
1e59de90 | 631 | std::intmax_t expon2 = 1 + ilogb(b); |
b32b8144 FG |
632 | T upper2 = scalbn(T(1), expon2 - 1); |
633 | result = float_distance(upper2, b); | |
634 | result += (expon2 - expon - 1) * scalbn(T(1), std::numeric_limits<T>::digits - 1); | |
635 | } | |
636 | // | |
637 | // Use compensated double-double addition to avoid rounding | |
638 | // errors in the subtraction: | |
639 | // | |
640 | expon = std::numeric_limits<T>::digits - expon; | |
641 | T mb, x, y, z; | |
642 | if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) | |
643 | { | |
644 | // | |
645 | // Special case - either one end of the range is a denormal, or else the difference is. | |
646 | // The regular code will fail if we're using the SSE2 registers on Intel and either | |
647 | // the FTZ or DAZ flags are set. | |
648 | // | |
649 | T a2 = scalbn(a, std::numeric_limits<T>::digits); | |
650 | T b2 = scalbn(b, std::numeric_limits<T>::digits); | |
651 | mb = -(std::min)(T(scalbn(upper, std::numeric_limits<T>::digits)), b2); | |
652 | x = a2 + mb; | |
653 | z = x - a2; | |
654 | y = (a2 - (x - z)) + (mb - z); | |
655 | ||
656 | expon -= std::numeric_limits<T>::digits; | |
657 | } | |
658 | else | |
659 | { | |
660 | mb = -(std::min)(upper, b); | |
661 | x = a + mb; | |
662 | z = x - a; | |
663 | y = (a - (x - z)) + (mb - z); | |
664 | } | |
665 | if(x < 0) | |
666 | { | |
667 | x = -x; | |
668 | y = -y; | |
669 | } | |
670 | result += scalbn(x, expon) + scalbn(y, expon); | |
671 | // | |
672 | // Result must be an integer: | |
673 | // | |
1e59de90 | 674 | BOOST_MATH_ASSERT(result == floor(result)); |
b32b8144 FG |
675 | return result; |
676 | } // float_distance_imp | |
677 | ||
678 | } // namespace detail | |
7c673cae FG |
679 | |
680 | template <class T, class U, class Policy> | |
681 | inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol) | |
682 | { | |
20effc67 TL |
683 | // |
684 | // We allow ONE of a and b to be an integer type, otherwise both must be the SAME type. | |
685 | // | |
1e59de90 TL |
686 | static_assert( |
687 | (std::is_same<T, U>::value | |
688 | || (std::is_integral<T>::value && !std::is_integral<U>::value) | |
689 | || (!std::is_integral<T>::value && std::is_integral<U>::value) | |
20effc67 TL |
690 | || (std::numeric_limits<T>::is_specialized && std::numeric_limits<U>::is_specialized |
691 | && (std::numeric_limits<T>::digits == std::numeric_limits<U>::digits) | |
692 | && (std::numeric_limits<T>::radix == std::numeric_limits<U>::radix) | |
693 | && !std::numeric_limits<T>::is_integer && !std::numeric_limits<U>::is_integer)), | |
694 | "Float distance between two different floating point types is undefined."); | |
695 | ||
1e59de90 | 696 | BOOST_IF_CONSTEXPR (!std::is_same<T, U>::value) |
20effc67 | 697 | { |
1e59de90 | 698 | BOOST_IF_CONSTEXPR(std::is_integral<T>::value) |
20effc67 TL |
699 | { |
700 | return float_distance(static_cast<U>(a), b, pol); | |
701 | } | |
702 | else | |
703 | { | |
704 | return float_distance(a, static_cast<T>(b), pol); | |
705 | } | |
706 | } | |
707 | else | |
708 | { | |
709 | typedef typename tools::promote_args<T, U>::type result_type; | |
1e59de90 | 710 | return detail::float_distance_imp(detail::normalize_value(static_cast<result_type>(a), typename detail::has_hidden_guard_digits<result_type>::type()), detail::normalize_value(static_cast<result_type>(b), typename detail::has_hidden_guard_digits<result_type>::type()), std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
20effc67 | 711 | } |
7c673cae FG |
712 | } |
713 | ||
714 | template <class T, class U> | |
715 | typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b) | |
716 | { | |
717 | return boost::math::float_distance(a, b, policies::policy<>()); | |
718 | } | |
719 | ||
720 | namespace detail{ | |
721 | ||
722 | template <class T, class Policy> | |
1e59de90 | 723 | T float_advance_imp(T val, int distance, const std::true_type&, const Policy& pol) |
7c673cae FG |
724 | { |
725 | BOOST_MATH_STD_USING | |
726 | // | |
727 | // Error handling: | |
728 | // | |
729 | static const char* function = "float_advance<%1%>(%1%, int)"; | |
730 | ||
731 | int fpclass = (boost::math::fpclassify)(val); | |
732 | ||
733 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
734 | return policies::raise_domain_error<T>( | |
735 | function, | |
736 | "Argument val must be finite, but got %1%", val, pol); | |
737 | ||
738 | if(val < 0) | |
739 | return -float_advance(-val, -distance, pol); | |
740 | if(distance == 0) | |
741 | return val; | |
742 | if(distance == 1) | |
743 | return float_next(val, pol); | |
744 | if(distance == -1) | |
745 | return float_prior(val, pol); | |
746 | ||
747 | if(fabs(val) < detail::get_min_shift_value<T>()) | |
748 | { | |
749 | // | |
750 | // Special case: if the value of the least significant bit is a denorm, | |
751 | // implement in terms of float_next/float_prior. | |
752 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
753 | // | |
754 | if(distance > 0) | |
755 | { | |
756 | do{ val = float_next(val, pol); } while(--distance); | |
757 | } | |
758 | else | |
759 | { | |
760 | do{ val = float_prior(val, pol); } while(++distance); | |
761 | } | |
762 | return val; | |
763 | } | |
764 | ||
765 | int expon; | |
92f5a8d4 | 766 | (void)frexp(val, &expon); |
7c673cae FG |
767 | T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon); |
768 | if(val <= tools::min_value<T>()) | |
769 | { | |
770 | limit = sign(T(distance)) * tools::min_value<T>(); | |
771 | } | |
772 | T limit_distance = float_distance(val, limit); | |
773 | while(fabs(limit_distance) < abs(distance)) | |
774 | { | |
775 | distance -= itrunc(limit_distance); | |
776 | val = limit; | |
777 | if(distance < 0) | |
778 | { | |
779 | limit /= 2; | |
780 | expon--; | |
781 | } | |
782 | else | |
783 | { | |
784 | limit *= 2; | |
785 | expon++; | |
786 | } | |
787 | limit_distance = float_distance(val, limit); | |
788 | if(distance && (limit_distance == 0)) | |
789 | { | |
790 | return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol); | |
791 | } | |
792 | } | |
793 | if((0.5f == frexp(val, &expon)) && (distance < 0)) | |
794 | --expon; | |
795 | T diff = 0; | |
796 | if(val != 0) | |
797 | diff = distance * ldexp(T(1), expon - tools::digits<T>()); | |
798 | if(diff == 0) | |
799 | diff = distance * detail::get_smallest_value<T>(); | |
800 | return val += diff; | |
b32b8144 FG |
801 | } // float_advance_imp |
802 | // | |
803 | // Special version for bases other than 2: | |
804 | // | |
805 | template <class T, class Policy> | |
1e59de90 | 806 | T float_advance_imp(T val, int distance, const std::false_type&, const Policy& pol) |
b32b8144 | 807 | { |
1e59de90 TL |
808 | static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized."); |
809 | static_assert(std::numeric_limits<T>::radix != 2, "Type T must be specialized."); | |
7c673cae | 810 | |
b32b8144 FG |
811 | BOOST_MATH_STD_USING |
812 | // | |
813 | // Error handling: | |
814 | // | |
815 | static const char* function = "float_advance<%1%>(%1%, int)"; | |
816 | ||
817 | int fpclass = (boost::math::fpclassify)(val); | |
818 | ||
819 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
820 | return policies::raise_domain_error<T>( | |
821 | function, | |
822 | "Argument val must be finite, but got %1%", val, pol); | |
823 | ||
824 | if(val < 0) | |
825 | return -float_advance(-val, -distance, pol); | |
826 | if(distance == 0) | |
827 | return val; | |
828 | if(distance == 1) | |
829 | return float_next(val, pol); | |
830 | if(distance == -1) | |
831 | return float_prior(val, pol); | |
832 | ||
833 | if(fabs(val) < detail::get_min_shift_value<T>()) | |
834 | { | |
835 | // | |
836 | // Special case: if the value of the least significant bit is a denorm, | |
837 | // implement in terms of float_next/float_prior. | |
838 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
839 | // | |
840 | if(distance > 0) | |
841 | { | |
842 | do{ val = float_next(val, pol); } while(--distance); | |
843 | } | |
844 | else | |
845 | { | |
846 | do{ val = float_prior(val, pol); } while(++distance); | |
847 | } | |
848 | return val; | |
849 | } | |
850 | ||
1e59de90 | 851 | std::intmax_t expon = 1 + ilogb(val); |
b32b8144 FG |
852 | T limit = scalbn(T(1), distance < 0 ? expon - 1 : expon); |
853 | if(val <= tools::min_value<T>()) | |
854 | { | |
855 | limit = sign(T(distance)) * tools::min_value<T>(); | |
856 | } | |
857 | T limit_distance = float_distance(val, limit); | |
858 | while(fabs(limit_distance) < abs(distance)) | |
859 | { | |
860 | distance -= itrunc(limit_distance); | |
861 | val = limit; | |
862 | if(distance < 0) | |
863 | { | |
864 | limit /= std::numeric_limits<T>::radix; | |
865 | expon--; | |
866 | } | |
867 | else | |
868 | { | |
869 | limit *= std::numeric_limits<T>::radix; | |
870 | expon++; | |
871 | } | |
872 | limit_distance = float_distance(val, limit); | |
873 | if(distance && (limit_distance == 0)) | |
874 | { | |
875 | return policies::raise_evaluation_error<T>(function, "Internal logic failed while trying to increment floating point value %1%: most likely your FPU is in non-IEEE conforming mode.", val, pol); | |
876 | } | |
877 | } | |
878 | /*expon = 1 + ilogb(val); | |
879 | if((1 == scalbn(val, 1 + expon)) && (distance < 0)) | |
880 | --expon;*/ | |
881 | T diff = 0; | |
882 | if(val != 0) | |
883 | diff = distance * scalbn(T(1), expon - std::numeric_limits<T>::digits); | |
884 | if(diff == 0) | |
885 | diff = distance * detail::get_smallest_value<T>(); | |
886 | return val += diff; | |
887 | } // float_advance_imp | |
888 | ||
889 | } // namespace detail | |
7c673cae FG |
890 | |
891 | template <class T, class Policy> | |
892 | inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol) | |
893 | { | |
894 | typedef typename tools::promote_args<T>::type result_type; | |
1e59de90 | 895 | return detail::float_advance_imp(detail::normalize_value(static_cast<result_type>(val), typename detail::has_hidden_guard_digits<result_type>::type()), distance, std::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol); |
7c673cae FG |
896 | } |
897 | ||
898 | template <class T> | |
899 | inline typename tools::promote_args<T>::type float_advance(const T& val, int distance) | |
900 | { | |
901 | return boost::math::float_advance(val, distance, policies::policy<>()); | |
902 | } | |
903 | ||
b32b8144 | 904 | }} // boost math namespaces |
7c673cae FG |
905 | |
906 | #endif // BOOST_MATH_SPECIAL_NEXT_HPP |