• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 //!file
2 //! \brief floating-point comparison from Boost.Test
3 // Copyright Paul A. Bristow 2015.
4 // Copyright John Maddock 2015.
5 
6 // Use, modification and distribution are subject to the
7 // Boost Software License, Version 1.0.
8 // (See accompanying file LICENSE_1_0.txt
9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
10 
11 // Note that this file contains Quickbook mark-up as well as code
12 // and comments, don't change any of the special comment mark-ups!
13 
14 #include <boost/math/special_functions/relative_difference.hpp>
15 #include <boost/math/special_functions/next.hpp>
16 
17 #include <iostream>
18 #include <limits> // for std::numeric_limits<T>::epsilon().
19 
main()20 int main()
21 {
22   std::cout << "Compare floats using Boost.Math functions/classes" << std::endl;
23 
24 
25 //[compare_floats_using
26 /*`Some using statements will ensure that the functions we need are accessible.
27 */
28 
29   using namespace boost::math;
30 
31 //`or
32 
33   using boost::math::relative_difference;
34   using boost::math::epsilon_difference;
35   using boost::math::float_next;
36   using boost::math::float_prior;
37 
38 //] [/compare_floats_using]
39 
40 
41 //[compare_floats_example_1
42 /*`The following examples display values with all possibly significant digits.
43 Newer compilers should provide `std::numeric_limits<FPT>::max_digits10`
44 for this purpose, and here we use `float` precision where `max_digits10` = 9
45 to avoid displaying a distracting number of decimal digits.
46 
47 [note Older compilers can use this formula to calculate `max_digits10`
48 from `std::numeric_limits<FPT>::digits10`:
49 __spaces `int max_digits10 = 2 + std::numeric_limits<FPT>::digits10 * 3010/10000;`
50 ] [/note]
51 
52 One can set the display including all trailing zeros
53 (helpful for this example to show all potentially significant digits),
54 and also to display `bool` values as words rather than integers:
55 */
56   std::cout.precision(std::numeric_limits<float>::max_digits10);
57   std::cout << std::boolalpha << std::showpoint << std::endl;
58 
59 //] [/compare_floats_example_1]
60 
61 //[compare_floats_example_2]
62 /*`
63 When comparing values that are ['quite close] or ['approximately equal],
64 we could use either `float_distance` or `relative_difference`/`epsilon_difference`, for example
65 with type `float`, these two values are adjacent to each other:
66 */
67 
68   float a = 1;
69   float b = 1 + std::numeric_limits<float>::epsilon();
70   std::cout << "a = " << a << std::endl;
71   std::cout << "b = " << b << std::endl;
72   std::cout << "float_distance = " << float_distance(a, b) << std::endl;
73   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
74   std::cout << "epsilon_difference = " << epsilon_difference(a, b) << std::endl;
75 
76 /*`
77 Which produces the output:
78 
79 [pre
80 a = 1.00000000
81 b = 1.00000012
82 float_distance = 1.00000000
83 relative_difference = 1.19209290e-007
84 epsilon_difference = 1.00000000
85 ]
86 */
87   //] [/compare_floats_example_2]
88 
89 //[compare_floats_example_3]
90 /*`
91 In the example above, it just so happens that the edit distance as measured by `float_distance`, and the
92 difference measured in units of epsilon were equal.  However, due to the way floating point
93 values are represented, that is not always the case:*/
94 
95   a = 2.0f / 3.0f;   // 2/3 inexactly represented as a float
96   b = float_next(float_next(float_next(a))); // 3 floating point values above a
97   std::cout << "a = " << a << std::endl;
98   std::cout << "b = " << b << std::endl;
99   std::cout << "float_distance = " << float_distance(a, b) << std::endl;
100   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
101   std::cout << "epsilon_difference = " << epsilon_difference(a, b) << std::endl;
102 
103 /*`
104 Which produces the output:
105 
106 [pre
107 a = 0.666666687
108 b = 0.666666865
109 float_distance = 3.00000000
110 relative_difference = 2.68220901e-007
111 epsilon_difference = 2.25000000
112 ]
113 
114 There is another important difference between `float_distance` and the
115 `relative_difference/epsilon_difference` functions in that `float_distance`
116 returns a signed result that reflects which argument is larger in magnitude,
117 where as `relative_difference/epsilon_difference` simply return an unsigned
118 value that represents how far apart the values are.  For example if we swap
119 the order of the arguments:
120 */
121 
122   std::cout << "float_distance = " << float_distance(b, a) << std::endl;
123   std::cout << "relative_difference = " << relative_difference(b, a) << std::endl;
124   std::cout << "epsilon_difference = " << epsilon_difference(b, a) << std::endl;
125 
126   /*`
127   The output is now:
128 
129   [pre
130   float_distance = -3.00000000
131   relative_difference = 2.68220901e-007
132   epsilon_difference = 2.25000000
133   ]
134 */
135   //] [/compare_floats_example_3]
136 
137 //[compare_floats_example_4]
138 /*`
139 Zeros are always treated as equal, as are infinities as long as they have the same sign:*/
140 
141   a = 0;
142   b = -0;  // signed zero
143   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
144   a = b = std::numeric_limits<float>::infinity();
145   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
146   std::cout << "relative_difference = " << relative_difference(a, -b) << std::endl;
147 
148 /*`
149 Which produces the output:
150 
151 [pre
152 relative_difference = 0.000000000
153 relative_difference = 0.000000000
154 relative_difference = 3.40282347e+038
155 ]
156 */
157 //] [/compare_floats_example_4]
158 
159 //[compare_floats_example_5]
160 /*`
161 Note that finite values are always infinitely far away from infinities even if those finite values are very large:*/
162 
163   a = (std::numeric_limits<float>::max)();
164   b = std::numeric_limits<float>::infinity();
165   std::cout << "a = " << a << std::endl;
166   std::cout << "b = " << b << std::endl;
167   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
168   std::cout << "epsilon_difference = " << epsilon_difference(a, b) << std::endl;
169 
170 /*`
171 Which produces the output:
172 
173 [pre
174 a = 3.40282347e+038
175 b = 1.#INF0000
176 relative_difference = 3.40282347e+038
177 epsilon_difference = 3.40282347e+038
178 ]
179 */
180 //] [/compare_floats_example_5]
181 
182 //[compare_floats_example_6]
183 /*`
184 Finally, all denormalized values and zeros are treated as being effectively equal:*/
185 
186   a = std::numeric_limits<float>::denorm_min();
187   b = a * 2;
188   std::cout << "a = " << a << std::endl;
189   std::cout << "b = " << b << std::endl;
190   std::cout << "float_distance = " << float_distance(a, b) << std::endl;
191   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
192   std::cout << "epsilon_difference = " << epsilon_difference(a, b) << std::endl;
193   a = 0;
194   std::cout << "a = " << a << std::endl;
195   std::cout << "b = " << b << std::endl;
196   std::cout << "float_distance = " << float_distance(a, b) << std::endl;
197   std::cout << "relative_difference = " << relative_difference(a, b) << std::endl;
198   std::cout << "epsilon_difference = " << epsilon_difference(a, b) << std::endl;
199 
200 /*`
201 Which produces the output:
202 
203 [pre
204 a = 1.40129846e-045
205 b = 2.80259693e-045
206 float_distance = 1.00000000
207 relative_difference = 0.000000000
208 epsilon_difference = 0.000000000
209 a = 0.000000000
210 b = 2.80259693e-045
211 float_distance = 2.00000000
212 relative_difference = 0.000000000
213 epsilon_difference = 0.000000000]
214 
215 Notice how, in the above example, two denormalized values that are a factor of 2 apart are
216 none the less only one representation apart!
217 
218 */
219 //] [/compare_floats_example_6]
220 
221 
222 #if 0
223 //[old_compare_floats_example_3
224 //`The simplest use is to compare two values with a tolerance thus:
225 
226   bool is_close = is_close_to(1.F, 1.F + epsilon, epsilon); // One epsilon apart is close enough.
227   std::cout << "is_close_to(1.F, 1.F + epsilon, epsilon); is " << is_close << std::endl; // true
228 
229   is_close = is_close_to(1.F, 1.F + 2 * epsilon, epsilon); // Two epsilon apart isn't close enough.
230   std::cout << "is_close_to(1.F, 1.F + epsilon, epsilon); is " << is_close << std::endl; // false
231 
232 /*`
233 [note The type FPT of the tolerance and the type of the values [*must match].
234 
235 So `is_close(0.1F, 1., 1.)` will fail to compile because "template parameter 'FPT' is ambiguous".
236 Always provide the same type, using `static_cast<FPT>` if necessary.]
237 */
238 
239 
240 /*`An instance of class `close_at_tolerance` is more convenient
241 when multiple tests with the same conditions are planned.
242 A class that stores a tolerance of three epsilon (and the default ['strong] test) is:
243 */
244 
245   close_at_tolerance<float> three_rounds(3 * epsilon); // 'strong' by default.
246 
247 //`and we can confirm these settings:
248 
249   std::cout << "fraction_tolerance = "
250     << three_rounds.fraction_tolerance()
251     << std::endl; // +3.57627869e-007
252   std::cout << "strength = "
253     << (three_rounds.strength() == FPC_STRONG ? "strong" : "weak")
254     << std::endl; // strong
255 
256 //`To start, let us use two values that are truly equal (having identical bit patterns)
257 
258   float a = 1.23456789F;
259   float b = 1.23456789F;
260 
261 //`and make a comparison using our 3*epsilon `three_rounds` functor:
262 
263   bool close = three_rounds(a, b);
264   std::cout << "three_rounds(a, b) = " << close << std::endl; // true
265 
266 //`Unsurprisingly, the result is true, and the failed fraction is zero.
267 
268   std::cout << "failed_fraction = " << three_rounds.failed_fraction() << std::endl;
269 
270 /*`To get some nearby values, it is convenient to use the Boost.Math __next_float functions,
271 for which we need an include
272 
273   #include <boost/math/special_functions/next.hpp>
274 
275 and some using declarations:
276 */
277 
278   using boost::math::float_next;
279   using boost::math::float_prior;
280   using boost::math::nextafter;
281   using boost::math::float_distance;
282 
283 //`To add a few __ulp to one value:
284   b = float_next(a); // Add just one ULP to a.
285   b = float_next(b); // Add another one ULP.
286   b = float_next(b); // Add another one ULP.
287   // 3 epsilon would pass.
288   b = float_next(b); // Add another one ULP.
289 
290 //`and repeat our comparison:
291 
292   close = three_rounds(a, b);
293   std::cout << "three_rounds(a, b) = " << close << std::endl; // false
294   std::cout << "failed_fraction = " << three_rounds.failed_fraction()
295     << std::endl;  // abs(u-v) / abs(v) = 3.86237957e-007
296 
297 //`We can also 'measure' the number of bits different using the `float_distance` function:
298 
299   std::cout << "float_distance = " << float_distance(a, b) << std::endl; // 4
300 
301 /*`Now consider two values that are much further apart
302 than one might expect from ['computational noise],
303 perhaps the result of two measurements of some physical property like length
304 where an uncertainty of a percent or so might be expected.
305 */
306   float fp1 = 0.01000F;
307   float fp2 = 0.01001F; // Slightly different.
308 
309   float tolerance = 0.0001F;
310 
311   close_at_tolerance<float> strong(epsilon); // Default is strong.
312   bool rs = strong(fp1, fp2);
313   std::cout << "strong(fp1, fp2) is " << rs << std::endl;
314 
315 //`Or we could contrast using the ['weak] criterion:
316   close_at_tolerance<float> weak(epsilon, FPC_WEAK); // Explicitly weak.
317   bool rw = weak(fp1, fp2); //
318   std::cout << "weak(fp1, fp2) is " << rw << std::endl;
319 
320 //`We can also construct, setting tolerance and strength, and compare in one statement:
321 
322   std::cout << a << " #= " << b << " is "
323     << close_at_tolerance<float>(epsilon, FPC_STRONG)(a, b) << std::endl;
324   std::cout << a << " ~= " << b << " is "
325     << close_at_tolerance<float>(epsilon, FPC_WEAK)(a, b) << std::endl;
326 
327 //`but this has little advantage over using function `is_close_to` directly.
328 
329 //] [/old_compare_floats_example_3]
330 
331 
332 /*When the floating-point values become very small and near zero, using
333 //a relative test becomes unhelpful because one is dividing by zero or tiny,
334 
335 //Instead, an absolute test is needed, comparing one (or usually both) values with zero,
336 //using a tolerance.
337 //This is provided by the `small_with_tolerance` class and `is_small` function.
338 
339   namespace boost {
340   namespace math {
341   namespace fpc {
342 
343 
344   template<typename FPT>
345   class small_with_tolerance
346   {
347   public:
348   // Public typedefs.
349   typedef bool result_type;
350 
351   // Constructor.
352   explicit small_with_tolerance(FPT tolerance); // tolerance >= 0
353 
354   // Functor
355   bool operator()(FPT value) const; // return true if <= absolute tolerance (near zero).
356   };
357 
358   template<typename FPT>
359   bool
360   is_small(FPT value, FPT tolerance); // return true if value <= absolute tolerance (near zero).
361 
362   }}} // namespaces.
363 
364 /*`
365 [note The type FPT of the tolerance and the type of the value [*must match].
366 
367 So `is_small(0.1F, 0.000001)` will fail to compile because "template parameter 'FPT' is ambiguous".
368 Always provide the same type, using `static_cast<FPT>` if necessary.]
369 
370 A few values near zero are tested with varying tolerance below.
371 */
372 //[compare_floats_small_1
373 
374   float c = 0;
375   std::cout << "0 is_small " << is_small(c, epsilon) << std::endl; // true
376 
377   c = std::numeric_limits<float>::denorm_min(); // 1.40129846e-045
378   std::cout << "denorm_ min =" << c << ", is_small is " << is_small(c, epsilon) << std::endl; // true
379 
380   c = (std::numeric_limits<float>::min)(); // 1.17549435e-038
381   std::cout << "min = " << c << ", is_small is " << is_small(c, epsilon) << std::endl; // true
382 
383   c = 1 * epsilon; // 1.19209290e-007
384   std::cout << "epsilon = " << c << ", is_small is " << is_small(c, epsilon) << std::endl; // false
385 
386   c = 1 * epsilon; // 1.19209290e-007
387   std::cout << "2 epsilon = " << c << ", is_small is " << is_small(c, 2 * epsilon) << std::endl; // true
388 
389   c = 2 * epsilon; //2.38418579e-007
390   std::cout << "4 epsilon = " << c << ", is_small is " << is_small(c, 2 * epsilon) << std::endl; // false
391 
392   c = 0.00001F;
393   std::cout << "0.00001 = " << c << ", is_small is " << is_small(c, 0.0001F) << std::endl; // true
394 
395   c = -0.00001F;
396   std::cout << "0.00001 = " << c << ", is_small is " << is_small(c, 0.0001F) << std::endl; // true
397 
398 /*`Using the class `small_with_tolerance` allows storage of the tolerance,
399 convenient if you make repeated tests with the same tolerance.
400 */
401 
402   small_with_tolerance<float>my_test(0.01F);
403 
404   std::cout << "my_test(0.001F) is " << my_test(0.001F) << std::endl; // true
405   std::cout << "my_test(0.001F) is " << my_test(0.01F) << std::endl; // false
406 
407   //] [/compare_floats_small_1]
408 #endif
409   return 0;
410 }  // int main()
411 
412 /*
413 
414 Example output is:
415 
416 //[compare_floats_output
417 Compare floats using Boost.Test functions/classes
418 
419 float epsilon = 1.19209290e-007
420 is_close_to(1.F, 1.F + epsilon, epsilon); is true
421 is_close_to(1.F, 1.F + epsilon, epsilon); is false
422 fraction_tolerance = 3.57627869e-007
423 strength = strong
424 three_rounds(a, b) = true
425 failed_fraction = 0.000000000
426 three_rounds(a, b) = false
427 failed_fraction = 3.86237957e-007
428 float_distance = 4.00000000
429 strong(fp1, fp2) is false
430 weak(fp1, fp2) is false
431 1.23456788 #= 1.23456836 is false
432 1.23456788 ~= 1.23456836 is false
433 0 is_small true
434 denorm_ min =1.40129846e-045, is_small is true
435 min = 1.17549435e-038, is_small is true
436 epsilon = 1.19209290e-007, is_small is false
437 2 epsilon = 1.19209290e-007, is_small is true
438 4 epsilon = 2.38418579e-007, is_small is false
439 0.00001 = 9.99999975e-006, is_small is true
440 0.00001 = -9.99999975e-006, is_small is true
441 my_test(0.001F) is true
442 
443 my_test(0.001F) is false//] [/compare_floats_output]
444 */
445