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/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>
17
18 #ifdef BOOST_MSVC
19 #include <float.h>
20 #endif
21
22 namespace boost{ namespace math{
23
24 namespace detail{
25
26 template <class T>
get_smallest_value(mpl::true_ const &)27 inline T get_smallest_value(mpl::true_ const&)
28 {
29 return std::numeric_limits<T>::denorm_min();
30 }
31
32 template <class T>
get_smallest_value(mpl::false_ const &)33 inline T get_smallest_value(mpl::false_ const&)
34 {
35 return tools::min_value<T>();
36 }
37
38 template <class T>
get_smallest_value()39 inline T get_smallest_value()
40 {
41 #if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
42 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
43 #else
44 return get_smallest_value<T>(mpl::bool_<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
45 #endif
46 }
47
48 }
49
50 template <class T, class Policy>
51 T float_next(const T& val, const Policy& pol)
52 {
53 BOOST_MATH_STD_USING
54 int expon;
55 static const char* function = "float_next<%1%>(%1%)";
56
57 if(!(boost::math::isfinite)(val))
58 return policies::raise_domain_error<T>(
59 function,
60 "Argument must be finite, but got %1%", val, pol);
61
62 if(val >= tools::max_value<T>())
63 return policies::raise_overflow_error<T>(function, 0, pol);
64
65 if(val == 0)
66 return detail::get_smallest_value<T>();
67
68 if(-0.5f == frexp(val, &expon))
69 --expon; // reduce exponent when val is a power of two, and negative.
70 T diff = ldexp(T(1), expon - tools::digits<T>());
71 if(diff == 0)
72 diff = detail::get_smallest_value<T>();
73 return val + diff;
74 }
75
76 #ifdef BOOST_MSVC
77 template <class Policy>
float_next(const double & val,const Policy & pol)78 inline double float_next(const double& val, const Policy& pol)
79 {
80 static const char* function = "float_next<%1%>(%1%)";
81
82 if(!(boost::math::isfinite)(val))
83 return policies::raise_domain_error<double>(
84 function,
85 "Argument must be finite, but got %1%", val, pol);
86
87 if(val >= tools::max_value<double>())
88 return policies::raise_overflow_error<double>(function, 0, pol);
89
90 return ::_nextafter(val, tools::max_value<double>());
91 }
92 #endif
93
94 template <class T>
float_next(const T & val)95 inline T float_next(const T& val)
96 {
97 return float_next(val, policies::policy<>());
98 }
99
100 template <class T, class Policy>
101 T float_prior(const T& val, const Policy& pol)
102 {
103 BOOST_MATH_STD_USING
104 int expon;
105 static const char* function = "float_prior<%1%>(%1%)";
106
107 if(!(boost::math::isfinite)(val))
108 return policies::raise_domain_error<T>(
109 function,
110 "Argument must be finite, but got %1%", val, pol);
111
112 if(val <= -tools::max_value<T>())
113 return -policies::raise_overflow_error<T>(function, 0, pol);
114
115 if(val == 0)
116 return -detail::get_smallest_value<T>();
117
118 T remain = frexp(val, &expon);
119 if(remain == 0.5)
120 --expon; // when val is a power of two we must reduce the exponent
121 T diff = ldexp(T(1), expon - tools::digits<T>());
122 if(diff == 0)
123 diff = detail::get_smallest_value<T>();
124 return val - diff;
125 }
126
127 #ifdef BOOST_MSVC
128 template <class Policy>
float_prior(const double & val,const Policy & pol)129 inline double float_prior(const double& val, const Policy& pol)
130 {
131 static const char* function = "float_prior<%1%>(%1%)";
132
133 if(!(boost::math::isfinite)(val))
134 return policies::raise_domain_error<double>(
135 function,
136 "Argument must be finite, but got %1%", val, pol);
137
138 if(val <= -tools::max_value<double>())
139 return -policies::raise_overflow_error<double>(function, 0, pol);
140
141 return ::_nextafter(val, -tools::max_value<double>());
142 }
143 #endif
144
145 template <class T>
float_prior(const T & val)146 inline T float_prior(const T& val)
147 {
148 return float_prior(val, policies::policy<>());
149 }
150
151 template <class T, class Policy>
nextafter(const T & val,const T & direction,const Policy & pol)152 inline T nextafter(const T& val, const T& direction, const Policy& pol)
153 {
154 return val < direction ? boost::math::float_next(val, pol) : val == direction ? val : boost::math::float_prior(val, pol);
155 }
156
157 template <class T>
nextafter(const T & val,const T & direction)158 inline T nextafter(const T& val, const T& direction)
159 {
160 return nextafter(val, direction, policies::policy<>());
161 }
162
163 template <class T, class Policy>
164 T float_distance(const T& a, const T& b, const Policy& pol)
165 {
166 BOOST_MATH_STD_USING
167 //
168 // Error handling:
169 //
170 static const char* function = "float_distance<%1%>(%1%, %1%)";
171 if(!(boost::math::isfinite)(a))
172 return policies::raise_domain_error<T>(
173 function,
174 "Argument a must be finite, but got %1%", a, pol);
175 if(!(boost::math::isfinite)(b))
176 return policies::raise_domain_error<T>(
177 function,
178 "Argument b must be finite, but got %1%", b, pol);
179 //
180 // Special cases:
181 //
182 if(a > b)
183 return -float_distance(b, a);
184 if(a == b)
185 return 0;
186 if(a == 0)
187 return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol));
188 if(b == 0)
189 return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(a) * detail::get_smallest_value<T>()), a, pol));
190 if(boost::math::sign(a) != boost::math::sign(b))
191 return 2 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol))
192 + fabs(float_distance(static_cast<T>(boost::math::sign(a) * detail::get_smallest_value<T>()), a, pol));
193 //
194 // By the time we get here, both a and b must have the same sign, we want
195 // b > a and both postive for the following logic:
196 //
197 if(a < 0)
198 return float_distance(static_cast<T>(-b), static_cast<T>(-a));
199
200 BOOST_ASSERT(a >= 0);
201 BOOST_ASSERT(b >= a);
202
203 int expon;
204 //
205 // Note that if a is a denorm then the usual formula fails
206 // because we actually have fewer than tools::digits<T>()
207 // significant bits in the representation:
208 //
209 frexp(((boost::math::fpclassify)(a) == FP_SUBNORMAL) ? tools::min_value<T>() : a, &expon);
210 T upper = ldexp(T(1), expon);
211 T result = 0;
212 expon = tools::digits<T>() - expon;
213 //
214 // If b is greater than upper, then we *must* split the calculation
215 // as the size of the ULP changes with each order of magnitude change:
216 //
217 if(b > upper)
218 {
219 result = float_distance(upper, b);
220 }
221 //
222 // Use compensated double-double addition to avoid rounding
223 // errors in the subtraction:
224 //
225 T mb = -(std::min)(upper, b);
226 T x = a + mb;
227 T z = x - a;
228 T y = (a - (x - z)) + (mb - z);
229 if(x < 0)
230 {
231 x = -x;
232 y = -y;
233 }
234 result += ldexp(x, expon) + ldexp(y, expon);
235 //
236 // Result must be an integer:
237 //
238 BOOST_ASSERT(result == floor(result));
239 return result;
240 }
241
242 template <class T>
float_distance(const T & a,const T & b)243 T float_distance(const T& a, const T& b)
244 {
245 return boost::math::float_distance(a, b, policies::policy<>());
246 }
247
248 template <class T, class Policy>
249 T float_advance(T val, int distance, const Policy& pol)
250 {
251 //
252 // Error handling:
253 //
254 static const char* function = "float_advance<%1%>(%1%, int)";
255 if(!(boost::math::isfinite)(val))
256 return policies::raise_domain_error<T>(
257 function,
258 "Argument val must be finite, but got %1%", val, pol);
259
260 if(val < 0)
261 return -float_advance(-val, -distance, pol);
262 if(distance == 0)
263 return val;
264 if(distance == 1)
265 return float_next(val, pol);
266 if(distance == -1)
267 return float_prior(val, pol);
268 BOOST_MATH_STD_USING
269 int expon;
270 frexp(val, &expon);
271 T limit = ldexp((distance < 0 ? T(0.5f) : T(1)), expon);
272 if(val <= tools::min_value<T>())
273 {
274 limit = sign(T(distance)) * tools::min_value<T>();
275 }
276 T limit_distance = float_distance(val, limit);
277 while(fabs(limit_distance) < abs(distance))
278 {
279 distance -= itrunc(limit_distance);
280 val = limit;
281 if(distance < 0)
282 {
283 limit /= 2;
284 expon--;
285 }
286 else
287 {
288 limit *= 2;
289 expon++;
290 }
291 limit_distance = float_distance(val, limit);
292 }
293 if((0.5f == frexp(val, &expon)) && (distance < 0))
294 --expon;
295 T diff = 0;
296 if(val != 0)
297 diff = distance * ldexp(T(1), expon - tools::digits<T>());
298 if(diff == 0)
299 diff = distance * detail::get_smallest_value<T>();
300 return val += diff;
301 }
302
303 template <class T>
float_advance(const T & val,int distance)304 inline T float_advance(const T& val, int distance)
305 {
306 return boost::math::float_advance(val, distance, policies::policy<>());
307 }
308
309 }} // namespaces
310
311 #endif // BOOST_MATH_SPECIAL_NEXT_HPP
312
313