1 // (C) Copyright John Maddock 2015.
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_ULP_HPP
7 #define BOOST_MATH_SPECIAL_ULP_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/next.hpp>
17
18 namespace boost{ namespace math{ namespace detail{
19
20 template <class T, class Policy>
21 T ulp_imp(const T& val, const boost::true_type&, const Policy& pol)
22 {
23 BOOST_MATH_STD_USING
24 int expon;
25 static const char* function = "ulp<%1%>(%1%)";
26
27 int fpclass = (boost::math::fpclassify)(val);
28
29 if(fpclass == (int)FP_NAN)
30 {
31 return policies::raise_domain_error<T>(
32 function,
33 "Argument must be finite, but got %1%", val, pol);
34 }
35 else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
36 {
37 return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
38 }
39 else if(fpclass == FP_ZERO)
40 return detail::get_smallest_value<T>();
41 //
42 // This code is almost the same as that for float_next, except for negative integers,
43 // where we preserve the relation ulp(x) == ulp(-x) as does Java:
44 //
45 frexp(fabs(val), &expon);
46 T diff = ldexp(T(1), expon - tools::digits<T>());
47 if(diff == 0)
48 diff = detail::get_smallest_value<T>();
49 return diff;
50 }
51 // non-binary version:
52 template <class T, class Policy>
53 T ulp_imp(const T& val, const boost::false_type&, const Policy& pol)
54 {
55 BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
56 BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
57 BOOST_MATH_STD_USING
58 int expon;
59 static const char* function = "ulp<%1%>(%1%)";
60
61 int fpclass = (boost::math::fpclassify)(val);
62
63 if(fpclass == (int)FP_NAN)
64 {
65 return policies::raise_domain_error<T>(
66 function,
67 "Argument must be finite, but got %1%", val, pol);
68 }
69 else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
70 {
71 return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
72 }
73 else if(fpclass == FP_ZERO)
74 return detail::get_smallest_value<T>();
75 //
76 // This code is almost the same as that for float_next, except for negative integers,
77 // where we preserve the relation ulp(x) == ulp(-x) as does Java:
78 //
79 expon = 1 + ilogb(fabs(val));
80 T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
81 if(diff == 0)
82 diff = detail::get_smallest_value<T>();
83 return diff;
84 }
85
86 }
87
88 template <class T, class Policy>
ulp(const T & val,const Policy & pol)89 inline typename tools::promote_args<T>::type ulp(const T& val, const Policy& pol)
90 {
91 typedef typename tools::promote_args<T>::type result_type;
92 return detail::ulp_imp(static_cast<result_type>(val), boost::integral_constant<bool, !std::numeric_limits<result_type>::is_specialized || (std::numeric_limits<result_type>::radix == 2)>(), pol);
93 }
94
95 template <class T>
ulp(const T & val)96 inline typename tools::promote_args<T>::type ulp(const T& val)
97 {
98 return ulp(val, policies::policy<>());
99 }
100
101
102 }} // namespaces
103
104 #endif // BOOST_MATH_SPECIAL_ULP_HPP
105
106