]>
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 | |
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 | ||
28 | namespace boost{ namespace math{ | |
29 | ||
30 | namespace detail{ | |
31 | ||
32 | template <class T> | |
33 | inline T get_smallest_value(mpl::true_ const&) | |
34 | { | |
35 | // | |
36 | // numeric_limits lies about denorms being present - particularly | |
37 | // when this can be turned on or off at runtime, as is the case | |
38 | // when using the SSE2 registers in DAZ or FTZ mode. | |
39 | // | |
40 | static const T m = std::numeric_limits<T>::denorm_min(); | |
41 | #ifdef BOOST_MATH_CHECK_SSE2 | |
42 | return (_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) ? tools::min_value<T>() : m;; | |
43 | #else | |
44 | return ((tools::min_value<T>() / 2) == 0) ? tools::min_value<T>() : m; | |
45 | #endif | |
46 | } | |
47 | ||
48 | template <class T> | |
49 | inline T get_smallest_value(mpl::false_ const&) | |
50 | { | |
51 | return tools::min_value<T>(); | |
52 | } | |
53 | ||
54 | template <class T> | |
55 | inline T get_smallest_value() | |
56 | { | |
57 | #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310) | |
58 | return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>()); | |
59 | #else | |
60 | return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>()); | |
61 | #endif | |
62 | } | |
63 | ||
64 | // | |
65 | // Returns the smallest value that won't generate denorms when | |
66 | // we calculate the value of the least-significant-bit: | |
67 | // | |
68 | template <class T> | |
69 | T get_min_shift_value(); | |
70 | ||
71 | template <class T> | |
72 | struct min_shift_initializer | |
73 | { | |
74 | struct init | |
75 | { | |
76 | init() | |
77 | { | |
78 | do_init(); | |
79 | } | |
80 | static void do_init() | |
81 | { | |
82 | get_min_shift_value<T>(); | |
83 | } | |
84 | void force_instantiate()const{} | |
85 | }; | |
86 | static const init initializer; | |
87 | static void force_instantiate() | |
88 | { | |
89 | initializer.force_instantiate(); | |
90 | } | |
91 | }; | |
92 | ||
93 | template <class T> | |
94 | const typename min_shift_initializer<T>::init min_shift_initializer<T>::initializer; | |
95 | ||
96 | ||
97 | template <class T> | |
98 | inline T get_min_shift_value() | |
99 | { | |
100 | BOOST_MATH_STD_USING | |
101 | static const T val = ldexp(tools::min_value<T>(), tools::digits<T>() + 1); | |
102 | min_shift_initializer<T>::force_instantiate(); | |
103 | ||
104 | return val; | |
105 | } | |
106 | ||
107 | template <class T, class Policy> | |
108 | T float_next_imp(const T& val, const Policy& pol) | |
109 | { | |
110 | BOOST_MATH_STD_USING | |
111 | int expon; | |
112 | static const char* function = "float_next<%1%>(%1%)"; | |
113 | ||
114 | int fpclass = (boost::math::fpclassify)(val); | |
115 | ||
116 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
117 | { | |
118 | if(val < 0) | |
119 | return -tools::max_value<T>(); | |
120 | return policies::raise_domain_error<T>( | |
121 | function, | |
122 | "Argument must be finite, but got %1%", val, pol); | |
123 | } | |
124 | ||
125 | if(val >= tools::max_value<T>()) | |
126 | return policies::raise_overflow_error<T>(function, 0, pol); | |
127 | ||
128 | if(val == 0) | |
129 | return detail::get_smallest_value<T>(); | |
130 | ||
131 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != -tools::min_value<T>())) | |
132 | { | |
133 | // | |
134 | // Special case: if the value of the least significant bit is a denorm, and the result | |
135 | // would not be a denorm, then shift the input, increment, and shift back. | |
136 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
137 | // | |
138 | return ldexp(float_next(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); | |
139 | } | |
140 | ||
141 | if(-0.5f == frexp(val, &expon)) | |
142 | --expon; // reduce exponent when val is a power of two, and negative. | |
143 | T diff = ldexp(T(1), expon - tools::digits<T>()); | |
144 | if(diff == 0) | |
145 | diff = detail::get_smallest_value<T>(); | |
146 | return val + diff; | |
147 | } | |
148 | ||
149 | } | |
150 | ||
151 | template <class T, class Policy> | |
152 | inline typename tools::promote_args<T>::type float_next(const T& val, const Policy& pol) | |
153 | { | |
154 | typedef typename tools::promote_args<T>::type result_type; | |
155 | return detail::float_next_imp(static_cast<result_type>(val), pol); | |
156 | } | |
157 | ||
158 | #if 0 //def BOOST_MSVC | |
159 | // | |
160 | // We used to use ::_nextafter here, but doing so fails when using | |
161 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own | |
162 | // - albeit slower - code instead as at least that gives the correct answer. | |
163 | // | |
164 | template <class Policy> | |
165 | inline double float_next(const double& val, const Policy& pol) | |
166 | { | |
167 | static const char* function = "float_next<%1%>(%1%)"; | |
168 | ||
169 | if(!(boost::math::isfinite)(val) && (val > 0)) | |
170 | return policies::raise_domain_error<double>( | |
171 | function, | |
172 | "Argument must be finite, but got %1%", val, pol); | |
173 | ||
174 | if(val >= tools::max_value<double>()) | |
175 | return policies::raise_overflow_error<double>(function, 0, pol); | |
176 | ||
177 | return ::_nextafter(val, tools::max_value<double>()); | |
178 | } | |
179 | #endif | |
180 | ||
181 | template <class T> | |
182 | inline typename tools::promote_args<T>::type float_next(const T& val) | |
183 | { | |
184 | return float_next(val, policies::policy<>()); | |
185 | } | |
186 | ||
187 | namespace detail{ | |
188 | ||
189 | template <class T, class Policy> | |
190 | T float_prior_imp(const T& val, const Policy& pol) | |
191 | { | |
192 | BOOST_MATH_STD_USING | |
193 | int expon; | |
194 | static const char* function = "float_prior<%1%>(%1%)"; | |
195 | ||
196 | int fpclass = (boost::math::fpclassify)(val); | |
197 | ||
198 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
199 | { | |
200 | if(val > 0) | |
201 | return tools::max_value<T>(); | |
202 | return policies::raise_domain_error<T>( | |
203 | function, | |
204 | "Argument must be finite, but got %1%", val, pol); | |
205 | } | |
206 | ||
207 | if(val <= -tools::max_value<T>()) | |
208 | return -policies::raise_overflow_error<T>(function, 0, pol); | |
209 | ||
210 | if(val == 0) | |
211 | return -detail::get_smallest_value<T>(); | |
212 | ||
213 | if((fpclass != (int)FP_SUBNORMAL) && (fpclass != (int)FP_ZERO) && (fabs(val) < detail::get_min_shift_value<T>()) && (val != tools::min_value<T>())) | |
214 | { | |
215 | // | |
216 | // Special case: if the value of the least significant bit is a denorm, and the result | |
217 | // would not be a denorm, then shift the input, increment, and shift back. | |
218 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
219 | // | |
220 | return ldexp(float_prior(T(ldexp(val, 2 * tools::digits<T>())), pol), -2 * tools::digits<T>()); | |
221 | } | |
222 | ||
223 | T remain = frexp(val, &expon); | |
224 | if(remain == 0.5) | |
225 | --expon; // when val is a power of two we must reduce the exponent | |
226 | T diff = ldexp(T(1), expon - tools::digits<T>()); | |
227 | if(diff == 0) | |
228 | diff = detail::get_smallest_value<T>(); | |
229 | return val - diff; | |
230 | } | |
231 | ||
232 | } | |
233 | ||
234 | template <class T, class Policy> | |
235 | inline typename tools::promote_args<T>::type float_prior(const T& val, const Policy& pol) | |
236 | { | |
237 | typedef typename tools::promote_args<T>::type result_type; | |
238 | return detail::float_prior_imp(static_cast<result_type>(val), pol); | |
239 | } | |
240 | ||
241 | #if 0 //def BOOST_MSVC | |
242 | // | |
243 | // We used to use ::_nextafter here, but doing so fails when using | |
244 | // the SSE2 registers if the FTZ or DAZ flags are set, so use our own | |
245 | // - albeit slower - code instead as at least that gives the correct answer. | |
246 | // | |
247 | template <class Policy> | |
248 | inline double float_prior(const double& val, const Policy& pol) | |
249 | { | |
250 | static const char* function = "float_prior<%1%>(%1%)"; | |
251 | ||
252 | if(!(boost::math::isfinite)(val) && (val < 0)) | |
253 | return policies::raise_domain_error<double>( | |
254 | function, | |
255 | "Argument must be finite, but got %1%", val, pol); | |
256 | ||
257 | if(val <= -tools::max_value<double>()) | |
258 | return -policies::raise_overflow_error<double>(function, 0, pol); | |
259 | ||
260 | return ::_nextafter(val, -tools::max_value<double>()); | |
261 | } | |
262 | #endif | |
263 | ||
264 | template <class T> | |
265 | inline typename tools::promote_args<T>::type float_prior(const T& val) | |
266 | { | |
267 | return float_prior(val, policies::policy<>()); | |
268 | } | |
269 | ||
270 | template <class T, class U, class Policy> | |
271 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction, const Policy& pol) | |
272 | { | |
273 | typedef typename tools::promote_args<T, U>::type result_type; | |
274 | return val < direction ? boost::math::float_next<result_type>(val, pol) : val == direction ? val : boost::math::float_prior<result_type>(val, pol); | |
275 | } | |
276 | ||
277 | template <class T, class U> | |
278 | inline typename tools::promote_args<T, U>::type nextafter(const T& val, const U& direction) | |
279 | { | |
280 | return nextafter(val, direction, policies::policy<>()); | |
281 | } | |
282 | ||
283 | namespace detail{ | |
284 | ||
285 | template <class T, class Policy> | |
286 | T float_distance_imp(const T& a, const T& b, const Policy& pol) | |
287 | { | |
288 | BOOST_MATH_STD_USING | |
289 | // | |
290 | // Error handling: | |
291 | // | |
292 | static const char* function = "float_distance<%1%>(%1%, %1%)"; | |
293 | if(!(boost::math::isfinite)(a)) | |
294 | return policies::raise_domain_error<T>( | |
295 | function, | |
296 | "Argument a must be finite, but got %1%", a, pol); | |
297 | if(!(boost::math::isfinite)(b)) | |
298 | return policies::raise_domain_error<T>( | |
299 | function, | |
300 | "Argument b must be finite, but got %1%", b, pol); | |
301 | // | |
302 | // Special cases: | |
303 | // | |
304 | if(a > b) | |
305 | return -float_distance(b, a, pol); | |
306 | if(a == b) | |
307 | return 0; | |
308 | if(a == 0) | |
309 | return 1 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)); | |
310 | if(b == 0) | |
311 | return 1 + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
312 | if(boost::math::sign(a) != boost::math::sign(b)) | |
313 | return 2 + fabs(float_distance(static_cast<T>((b < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), b, pol)) | |
314 | + fabs(float_distance(static_cast<T>((a < 0) ? T(-detail::get_smallest_value<T>()) : detail::get_smallest_value<T>()), a, pol)); | |
315 | // | |
316 | // By the time we get here, both a and b must have the same sign, we want | |
317 | // b > a and both postive for the following logic: | |
318 | // | |
319 | if(a < 0) | |
320 | return float_distance(static_cast<T>(-b), static_cast<T>(-a), pol); | |
321 | ||
322 | BOOST_ASSERT(a >= 0); | |
323 | BOOST_ASSERT(b >= a); | |
324 | ||
325 | int expon; | |
326 | // | |
327 | // Note that if a is a denorm then the usual formula fails | |
328 | // because we actually have fewer than tools::digits<T>() | |
329 | // significant bits in the representation: | |
330 | // | |
331 | frexp(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon); | |
332 | T upper = ldexp(T(1), expon); | |
333 | T result = 0; | |
334 | expon = tools::digits<T>() - expon; | |
335 | // | |
336 | // If b is greater than upper, then we *must* split the calculation | |
337 | // as the size of the ULP changes with each order of magnitude change: | |
338 | // | |
339 | if(b > upper) | |
340 | { | |
341 | result = float_distance(upper, b); | |
342 | } | |
343 | // | |
344 | // Use compensated double-double addition to avoid rounding | |
345 | // errors in the subtraction: | |
346 | // | |
347 | T mb, x, y, z; | |
348 | if(((boost::math::fpclassify)(a) == (int)FP_SUBNORMAL) || (b - a < tools::min_value<T>())) | |
349 | { | |
350 | // | |
351 | // Special case - either one end of the range is a denormal, or else the difference is. | |
352 | // The regular code will fail if we're using the SSE2 registers on Intel and either | |
353 | // the FTZ or DAZ flags are set. | |
354 | // | |
355 | T a2 = ldexp(a, tools::digits<T>()); | |
356 | T b2 = ldexp(b, tools::digits<T>()); | |
357 | mb = -(std::min)(T(ldexp(upper, tools::digits<T>())), b2); | |
358 | x = a2 + mb; | |
359 | z = x - a2; | |
360 | y = (a2 - (x - z)) + (mb - z); | |
361 | ||
362 | expon -= tools::digits<T>(); | |
363 | } | |
364 | else | |
365 | { | |
366 | mb = -(std::min)(upper, b); | |
367 | x = a + mb; | |
368 | z = x - a; | |
369 | y = (a - (x - z)) + (mb - z); | |
370 | } | |
371 | if(x < 0) | |
372 | { | |
373 | x = -x; | |
374 | y = -y; | |
375 | } | |
376 | result += ldexp(x, expon) + ldexp(y, expon); | |
377 | // | |
378 | // Result must be an integer: | |
379 | // | |
380 | BOOST_ASSERT(result == floor(result)); | |
381 | return result; | |
382 | } | |
383 | ||
384 | } | |
385 | ||
386 | template <class T, class U, class Policy> | |
387 | inline typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b, const Policy& pol) | |
388 | { | |
389 | typedef typename tools::promote_args<T, U>::type result_type; | |
390 | return detail::float_distance_imp(static_cast<result_type>(a), static_cast<result_type>(b), pol); | |
391 | } | |
392 | ||
393 | template <class T, class U> | |
394 | typename tools::promote_args<T, U>::type float_distance(const T& a, const U& b) | |
395 | { | |
396 | return boost::math::float_distance(a, b, policies::policy<>()); | |
397 | } | |
398 | ||
399 | namespace detail{ | |
400 | ||
401 | template <class T, class Policy> | |
402 | T float_advance_imp(T val, int distance, const Policy& pol) | |
403 | { | |
404 | BOOST_MATH_STD_USING | |
405 | // | |
406 | // Error handling: | |
407 | // | |
408 | static const char* function = "float_advance<%1%>(%1%, int)"; | |
409 | ||
410 | int fpclass = (boost::math::fpclassify)(val); | |
411 | ||
412 | if((fpclass == (int)FP_NAN) || (fpclass == (int)FP_INFINITE)) | |
413 | return policies::raise_domain_error<T>( | |
414 | function, | |
415 | "Argument val must be finite, but got %1%", val, pol); | |
416 | ||
417 | if(val < 0) | |
418 | return -float_advance(-val, -distance, pol); | |
419 | if(distance == 0) | |
420 | return val; | |
421 | if(distance == 1) | |
422 | return float_next(val, pol); | |
423 | if(distance == -1) | |
424 | return float_prior(val, pol); | |
425 | ||
426 | if(fabs(val) < detail::get_min_shift_value<T>()) | |
427 | { | |
428 | // | |
429 | // Special case: if the value of the least significant bit is a denorm, | |
430 | // implement in terms of float_next/float_prior. | |
431 | // This avoids issues with the Intel SSE2 registers when the FTZ or DAZ flags are set. | |
432 | // | |
433 | if(distance > 0) | |
434 | { | |
435 | do{ val = float_next(val, pol); } while(--distance); | |
436 | } | |
437 | else | |
438 | { | |
439 | do{ val = float_prior(val, pol); } while(++distance); | |
440 | } | |
441 | return val; | |
442 | } | |
443 | ||
444 | int expon; | |
445 | frexp(val, &expon); | |
446 | T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon); | |
447 | if(val <= tools::min_value<T>()) | |
448 | { | |
449 | limit = sign(T(distance)) * tools::min_value<T>(); | |
450 | } | |
451 | T limit_distance = float_distance(val, limit); | |
452 | while(fabs(limit_distance) < abs(distance)) | |
453 | { | |
454 | distance -= itrunc(limit_distance); | |
455 | val = limit; | |
456 | if(distance < 0) | |
457 | { | |
458 | limit /= 2; | |
459 | expon--; | |
460 | } | |
461 | else | |
462 | { | |
463 | limit *= 2; | |
464 | expon++; | |
465 | } | |
466 | limit_distance = float_distance(val, limit); | |
467 | if(distance && (limit_distance == 0)) | |
468 | { | |
469 | 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); | |
470 | } | |
471 | } | |
472 | if((0.5f == frexp(val, &expon)) && (distance < 0)) | |
473 | --expon; | |
474 | T diff = 0; | |
475 | if(val != 0) | |
476 | diff = distance * ldexp(T(1), expon - tools::digits<T>()); | |
477 | if(diff == 0) | |
478 | diff = distance * detail::get_smallest_value<T>(); | |
479 | return val += diff; | |
480 | } | |
481 | ||
482 | } | |
483 | ||
484 | template <class T, class Policy> | |
485 | inline typename tools::promote_args<T>::type float_advance(T val, int distance, const Policy& pol) | |
486 | { | |
487 | typedef typename tools::promote_args<T>::type result_type; | |
488 | return detail::float_advance_imp(static_cast<result_type>(val), distance, pol); | |
489 | } | |
490 | ||
491 | template <class T> | |
492 | inline typename tools::promote_args<T>::type float_advance(const T& val, int distance) | |
493 | { | |
494 | return boost::math::float_advance(val, distance, policies::policy<>()); | |
495 | } | |
496 | ||
497 | }} // namespaces | |
498 | ||
499 | #endif // BOOST_MATH_SPECIAL_NEXT_HPP | |
500 |