1[section:next_float Floating-Point Representation Distance (ULP), 2 and Finding Adjacent Floating-Point Values] 3 4[@http://en.wikipedia.org/wiki/Unit_in_the_last_place Unit of Least Precision or Unit in the Last Place] 5is the gap between two different, but as close as possible, floating-point numbers. 6 7Most decimal values, for example 0.1, cannot be exactly represented as floating-point values, 8but will be stored as the 9[@http://en.wikipedia.org/wiki/Floating_point#Representable_numbers.2C_conversion_and_rounding 10closest representable floating-point]. 11 12Functions are provided for finding adjacent greater and lesser floating-point values, 13and estimating the number of gaps between any two floating-point values. 14 15The floating-point type (FPT) must have has a fixed number of bits in the representation. 16The number of bits may set at runtime, but must be the same for all numbers. 17For example, __NTL_quad_float type (fixed 128-bit representation), 18__NTL_RR type (arbitrary but fixed decimal digits, default 150) or 19__multiprecision __cpp_dec_float and__cpp_bin_float are fixed at runtime, 20but [*not] a type that extends the representation to provide an exact representation 21for any number, for example [@http://keithbriggs.info/xrc.html XRC eXact Real in C]. 22 23The accuracy of mathematical functions can be assessed and displayed in terms of __ULP, 24often as a ulps plot or by binning the differences as a histogram. 25Samples are evaluated using the implementation under test and compared with 'known good' 26representation obtained using a more accurate method. Other implementations, often using 27arbitrary precision arithmetic, for example __WolframAlpha are one source of references 28values. The other method, used widely in Boost.Math special functions, it to carry out 29the same algorithm, but using a higher precision type, typically using Boost.Multiprecision 30types like `cpp_bin_float_quad` for 128-bit (about 35 decimal digit precision), or 31`cpp_bin_float_50` (for 50 decimal digit precision). 32 33When converted to a particular machine representation, say `double`, say using a `static_cast`, 34the value is the nearest representation possible for the `double` type. This value 35cannot be 'wrong' by more than half a __ulp, and can be obtained using the Boost.Math function `ulp`. 36(Unless the algorithm is fundamentally flawed, something that should be revealed by 'sanity' 37checks using some independent sources). 38 39See some discussion and example plots by Cleve Moler of Mathworks 40[@https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/ 41ulps plots reveal math-function accuracy]. 42 43[section:nextafter Finding the Next Representable Value in a Specific Direction (nextafter)] 44 45[h4 Synopsis] 46 47`` 48#include <boost/math/special_functions/next.hpp> 49`` 50 51 namespace boost{ namespace math{ 52 53 template <class FPT> 54 FPT nextafter(FPT val, FPT direction); 55 56 }} // namespaces 57 58[h4 Description - nextafter] 59 60This is an implementation of the `nextafter` function included in the C99 standard. 61(It is also effectively an implementation of the C99 `nexttoward` legacy function 62which differs only having a `long double` direction, 63and can generally serve in its place if required). 64 65[note The C99 functions must use suffixes f and l to distinguish `float` and `long double` versions. 66C++ uses the template mechanism instead.] 67 68Returns the next representable value after /x/ in the direction of /y/. If 69`x == y` then returns /x/. If /x/ is non-finite then returns the result of 70a __domain_error. If there is no such value in the direction of /y/ then 71returns an __overflow_error. 72 73[warning The template parameter FTP must be a floating-point type. 74An integer type, for example, will produce an unhelpful error message.] 75 76[tip Nearly always, you just want the next or prior representable value, 77so instead use `float_next` or `float_prior` below.] 78 79[h4 Examples - nextafter] 80 81The two representations using a 32-bit float either side of unity are: 82`` 83The nearest (exact) representation of 1.F is 1.00000000 84nextafter(1.F, 999) is 1.00000012 85nextafter(1/f, -999) is 0.99999994 86 87The nearest (not exact) representation of 0.1F is 0.100000001 88nextafter(0.1F, 10) is 0.100000009 89nextafter(0.1F, 10) is 0.099999994 90`` 91 92[endsect] [/section:nextafter Finding the Next Representable Value in a Specific Direction (nextafter)] 93 94[section:float_next Finding the Next Greater Representable Value (float_next)] 95 96[h4 Synopsis] 97 98`` 99#include <boost/math/special_functions/next.hpp> 100`` 101 102 namespace boost{ namespace math{ 103 104 template <class FPT> 105 FPT float_next(FPT val); 106 107 }} // namespaces 108 109[h4 Description - float_next] 110 111Returns the next representable value which is greater than /x/. 112If /x/ is non-finite then returns the result of 113a __domain_error. If there is no such value greater than /x/ then 114returns an __overflow_error. 115 116Has the same effect as 117 118 nextafter(val, (std::numeric_limits<FPT>::max)()); 119 120[endsect] [/section:float_next Finding the Next Greater Representable Value (float_prior)] 121 122[section:float_prior Finding the Next Smaller Representable Value (float_prior)] 123 124[h4 Synopsis] 125 126`` 127#include <boost/math/special_functions/next.hpp> 128`` 129 130 namespace boost{ namespace math{ 131 132 template <class FPT> 133 FPT float_prior(FPT val); 134 135 }} // namespaces 136 137 138[h4 Description - float_prior] 139 140Returns the next representable value which is less than /x/. 141If /x/ is non-finite then returns the result of 142a __domain_error. If there is no such value less than /x/ then 143returns an __overflow_error. 144 145Has the same effect as 146 147 nextafter(val, -(std::numeric_limits<FPT>::max)()); // Note most negative value -max. 148 149[endsect] [/section:float_prior Finding the Next Smaller Representable Value (float_prior)] 150 151[section:float_distance Calculating the Representation Distance 152 Between Two floating-point Values (ULP) float_distance] 153 154Function float_distance finds the number of gaps/bits/ULP between any two floating-point values. 155If the significands of floating-point numbers are viewed as integers, 156then their difference is the number of ULP/gaps/bits different. 157 158[h4 Synopsis] 159 160`` 161#include <boost/math/special_functions/next.hpp> 162`` 163 164 namespace boost{ namespace math{ 165 166 template <class FPT> 167 FPT float_distance(FPT a, FPT b); 168 169 }} // namespaces 170 171[h4 Description - float_distance] 172 173Returns the distance between /a/ and /b/: the result is always 174a signed integer value (stored in floating-point type FPT) 175representing the number of distinct representations between /a/ and /b/. 176 177Note that 178 179* `float_distance(a, a)` always returns 0. 180* `float_distance(float_next(a), a)` always returns -1. 181* `float_distance(float_prior(a), a)` always returns 1. 182 183The function `float_distance` is equivalent to calculating the number 184of ULP (Units in the Last Place) between /a/ and /b/ except that it 185returns a signed value indicating whether `a > b` or not. 186 187If the distance is too great then it may not be able 188to be represented as an exact integer by type FPT, 189but in practice this is unlikely to be a issue. 190 191[endsect] [/section:float_distance Calculating the Representation Distance 192 Between Two floating-point Values (ULP) float_distance] 193 194[section:float_advance Advancing a floating-point Value by a Specific 195Representation Distance (ULP) float_advance] 196 197Function `float_advance` advances a floating-point number by a specified number 198of ULP. 199 200[h4 Synopsis] 201 202`` 203#include <boost/math/special_functions/next.hpp> 204`` 205 206 namespace boost{ namespace math{ 207 208 template <class FPT> 209 FPT float_advance(FPT val, int distance); 210 211 }} // namespaces 212 213[h4 Description - float_advance] 214 215Returns a floating-point number /r/ such that `float_distance(val, r) == distance`. 216 217[endsect] [/section:float_advance] 218 219[section:ulp Obtaining the Size of a Unit In the Last Place - ULP] 220 221Function `ulp` gives the size of a unit-in-the-last-place for a specified floating-point value. 222 223[h4 Synopsis] 224 225`` 226#include <boost/math/special_functions/ulp.hpp> 227`` 228 229 namespace boost{ namespace math{ 230 231 template <class FPT> 232 FPT ulp(const FPT& x); 233 234 template <class FPT, class Policy> 235 FPT ulp(const FPT& x, const Policy&); 236 237 }} // namespaces 238 239[h4 Description - ulp] 240 241Returns one [@http://en.wikipedia.org/wiki/Unit_in_the_last_place unit in the last place] of ['x]. 242 243Corner cases are handled as follows: 244 245* If the argument is a NaN, then raises a __domain_error. 246* If the argument is an infinity, then raises an __overflow_error. 247* If the argument is zero then returns the smallest representable value: for example for type 248`double` this would be either `std::numeric_limits<double>::min()` or `std::numeric_limits<double>::denorm_min()` 249depending whether denormals are supported (which have the values `2.2250738585072014e-308` and `4.9406564584124654e-324` respectively). 250* If the result is too small to represent, then returns the smallest representable value. 251* Always returns a positive value such that `ulp(x) == ulp(-x)`. 252 253[*Important:] The behavior of this function is aligned to that of [@http://docs.oracle.com/javase/7/docs/api/java/lang/Math.html#ulp%28double%29 254Java's ulp function], please note 255however that this function should only ever be used for rough and ready calculations as there are enough 256corner cases to trap even careful programmers. In particular: 257 258* The function is asymmetrical, which is to say, given `u = ulp(x)` if `x > 0` then `x + u` is the 259next floating-point value, but `x - u` is not necessarily the previous value. Similarly, if 260`x < 0` then `x - u` is the previous floating-point value, but `x + u` is not necessarily the next 261value. The corner cases occur at power of 2 boundaries. 262* When the argument becomes very small, it may be that there is no floating-point value that 263represents one ULP. Whether this is the case or not depends not only on whether the hardware 264may ['sometimes] support denormals (as signalled by `std::numeric_limits<FPT>::has_denorm`), but also whether these are 265currently enabled at runtime (for example on SSE hardware, the DAZ or FTZ flags will disable denormal support). 266In this situation, the `ulp` function may return a value that is many orders of magnitude too large. 267 268In light of the issues above, we recommend that: 269 270* To move between adjacent floating-point values always use __float_next, __float_prior or __nextafter (`std::nextafter` 271is another candidate, but our experience is that this also often breaks depending which optimizations and 272hardware flags are in effect). 273* To move several floating-point values away use __float_advance. 274* To calculate the edit distance between two floats use __float_distance. 275 276There is none the less, one important use case for this function: 277 278If it is known that the true result of some function is x[sub t] and the calculated result 279is x[sub c], then the error measured in ulp is simply [^fabs(x[sub t] - x[sub c]) / ulp(x[sub t])]. 280 281[endsect] [/section ulp] 282 283[endsect] [/ section:next_float Floating-Point Representation Distance (ULP), 284 and Finding Adjacent Floating-Point Values] 285 286[/ 287 Copyright 2008 John Maddock and Paul A. Bristow. 288 Distributed under the Boost Software License, Version 1.0. 289 (See accompanying file LICENSE_1_0.txt or copy at 290 http://www.boost.org/LICENSE_1_0.txt). 291] 292 293