• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 // Copyright Paul A. Bristow 2013
2 // Copyright John Maddock 2013
3 // Copyright Christopher Kormanyos
4 
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt
8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
9 
10 // Examples of std::numeric_limits usage as snippets for multiprecision documentation at multiprecision.qbk.
11 
12 // Includes text as Quickbook comments.
13 
14 #include <boost/assert.hpp>
15 
16 #include <boost/math/constants/constants.hpp>
17 #include <boost/math/special_functions/nonfinite_num_facets.hpp>
18 
19 #include <boost/math/special_functions/factorials.hpp>
20 #include <boost/math/special_functions/next.hpp>
21 #include <boost/math/tools/precision.hpp>
22 #include <boost/multiprecision/cpp_dec_float.hpp> // is decimal.
23 #include <boost/multiprecision/cpp_bin_float.hpp> // is binary.
24 
25 #define BOOST_TEST_MAIN
26 #include <boost/test/unit_test.hpp> // Boost.Test
27 #include <boost/test/floating_point_comparison.hpp>
28 
29 #include <iostream>
30 #include <iomanip>
31 #include <string>
32 #include <sstream>
33 #include <limits> // numeric_limits
34 #include <iomanip>
35 #include <locale>
36 
37 // static long double const log10Two = 0.30102999566398119521373889472449L; // log10(2.)
38 // It is more portable useful to use a Boost macro
39 // See https://www.boost.org/doc/libs/release/libs/config/doc/html/boost_config/boost_macro_reference.html
40 BOOST_STATIC_CONSTEXPR long double log10Two = 0.30102999566398119521373889472449L;
41 // which expands to static constexpr on standard C++11 and up, but static const on earlier versions.
42 
43   /*`By default, output would only show the standard 6 decimal digits,
44  so set precision to show all 50 significant digits, including any trailing zeros.
45  This is generally useful to show the implicit precision of the type of the value.
46 */
47 
48 
49 
50 
51 template <typename T>
max_digits10()52 int max_digits10()
53 {
54    int significand_digits = std::numeric_limits<T>::digits;
55   // BOOST_CONSTEXPR_OR_CONST int significand_digits = std::numeric_limits<T>::digits;
56    return static_cast<int>(ceil(1 + significand_digits * log10Two));
57 } // template <typename T> int max_digits10()
58 
59 // Used to test max_digits10<>() function below.
60 //#define BOOST_NO_CXX11_NUMERIC_LIMITS
61 
BOOST_AUTO_TEST_CASE(test_numeric_limits_snips)62 BOOST_AUTO_TEST_CASE(test_numeric_limits_snips)
63 {
64 #if !(defined(CI_SUPPRESS_KNOWN_ISSUES) && defined(BOOST_MSVC) && (BOOST_MSVC == 1600))
65   try
66   {
67 
68 // Example of portable way to get `std::numeric_limits<T>::max_digits10`.
69 //[max_digits10_1
70 
71 /*`For example, to be portable (including obselete platforms) for type `T` where `T` may be:
72  `float`, `double`, `long double`, `128-bit quad type`, `cpp_bin_float_50` ...
73 */
74 
75   typedef float T;
76 
77 #if defined BOOST_NO_CXX11_NUMERIC_LIMITS
78    // No max_digits10 implemented.
79     std::cout.precision(max_digits10<T>());
80 #else
81   #if(_MSC_VER <= 1600)
82    //  The MSVC 2010 version had the wrong value for std::numeric_limits<float>::max_digits10.
83     std::cout.precision(max_digits10<T>());
84   #else // Use the C++11 max_digits10.
85      std::cout.precision(std::numeric_limits<T>::max_digits10);
86      std::cout.precision(std::numeric_limits<T>::digits10);
87      std::cout.setf(std::ios_base::showpoint); // Append any trailing zeros,
88      // or more memorably
89      std::cout << std::showpoint << std::endl; //
90   #endif
91 #endif
92 
93   std::cout << "std::cout.precision(max_digits10) = " << std::cout.precision() << std::endl; // 9
94 
95   double x = 1.2345678901234567889;
96 
97   std::cout << "x = " << x << std::endl; //
98 
99 /*`which should output:
100 
101   std::cout.precision(max_digits10) = 9
102   x = 1.23456789
103 */
104 
105 //] [/max_digits10_1]
106 
107   {
108 //[max_digits10_2
109 
110   double write = 2./3; // Any arbitrary value that cannot be represented exactly.
111   double read = 0;
112   std::stringstream s;
113   s.precision(std::numeric_limits<double>::digits10); // or `float64_t` for 64-bit IEE754 double.
114   s << write;
115   s >> read;
116   if(read != write)
117   {
118     std::cout <<  std::setprecision(std::numeric_limits<double>::digits10)
119       << read << " != " << write << std::endl;
120   }
121 
122 //] [/max_digits10_2]
123   // 0.666666666666667 != 0.666666666666667
124   }
125 
126   {
127 //[max_digits10_3
128 
129   double pi = boost::math::double_constants::pi;
130   std::cout.precision(std::numeric_limits<double>::max_digits10);
131   std::cout << pi << std::endl; // 3.1415926535897931
132 
133 //] [/max_digits10_3]
134   }
135   {
136 //[max_digits10_4
137 /*`and similarly for a much higher precision type:
138 */
139 
140   using namespace boost::multiprecision;
141 
142   typedef number<cpp_dec_float<50> > cpp_dec_float_50; // 50 decimal digits.
143 
144   using boost::multiprecision::cpp_dec_float_50;
145 
146   cpp_dec_float_50 pi = boost::math::constants::pi<cpp_dec_float_50>();
147   std::cout.precision(std::numeric_limits<cpp_dec_float_50>::max_digits10);
148   std::cout << pi << std::endl;
149   // 3.141592653589793238462643383279502884197169399375105820974944592307816406
150 //] [/max_digits10_4]
151   }
152 
153   {
154 //[max_digits10_5
155 
156   for (int i = 2; i < 15; i++)
157   {
158     std::cout << std::setw(std::numeric_limits<int>::max_digits10)
159       << boost::math::factorial<double>(i)  << std::endl;
160   }
161 
162 //] [/max_digits10_5]
163   }
164 
165   }
166   catch(const std::exception& ex)
167   {
168     std::cout << "Caught Exception " << ex.what() << std::endl;
169   }
170 
171   {
172 //[max_digits10_6
173 
174   typedef double T;
175 
176   bool denorm = std::numeric_limits<T>::denorm_min() < (std::numeric_limits<T>::min)();
177   BOOST_ASSERT(denorm);
178 
179 //] [/max_digits10_6]
180   }
181 
182   {
183     unsigned char c = 255;
184     std::cout << "char c = " << (int)c << std::endl;
185   }
186 
187   {
188 //[digits10_1
189     std::cout
190       << std::setw(std::numeric_limits<short>::digits10 +1 +1) // digits10+1, and +1 for sign.
191       << std::showpos << (std::numeric_limits<short>::max)() // +32767
192       << std::endl
193       << std::setw(std::numeric_limits<short>::digits10 +1 +1)
194       << (std::numeric_limits<short>::min)() << std::endl;   // -32767
195 //] [/digits10_1]
196   }
197 
198   {
199 //[digits10_2
200     std::cout
201       << std::setw(std::numeric_limits<unsigned short>::digits10 +1 +1) // digits10+1, and +1 for sign.
202       << std::showpos << (std::numeric_limits<unsigned short>::max)() //  65535
203       << std::endl
204       << std::setw(std::numeric_limits<unsigned short>::digits10 +1 +1) // digits10+1, and +1 for sign.
205       << (std::numeric_limits<unsigned short>::min)() << std::endl;   //      0
206 //] [/digits10_2]
207   }
208 
209   std::cout <<std::noshowpos << std::endl;
210 
211   {
212 //[digits10_3
213   std::cout.precision(std::numeric_limits<double>::max_digits10);
214   double d =  1e15;
215   double dp1 = d+1;
216   std::cout << d << "\n" << dp1 << std::endl;
217   // 1000000000000000
218   // 1000000000000001
219   std::cout <<  dp1 - d << std::endl; // 1
220 //] [/digits10_3]
221   }
222 
223   {
224 //[digits10_4
225   std::cout.precision(std::numeric_limits<double>::max_digits10);
226   double d =  1e16;
227   double dp1 = d+1;
228   std::cout << d << "\n" << dp1 << std::endl;
229   // 10000000000000000
230   // 10000000000000000
231     std::cout << dp1 - d << std::endl; // 0 !!!
232 //] [/digits10_4]
233   }
234 
235   {
236 //[epsilon_1
237   std::cout.precision(std::numeric_limits<double>::max_digits10);
238   double d = 1.;
239   double eps = std::numeric_limits<double>::epsilon();
240   double dpeps = d+eps;
241   std::cout << std::showpoint // Ensure all trailing zeros are shown.
242     << d << "\n"           // 1.0000000000000000
243     << dpeps << std::endl; // 2.2204460492503131e-016
244   std::cout << dpeps - d   // 1.0000000000000002
245     << std::endl;
246 //] [epsilon_1]
247   }
248 
249   {
250 //[epsilon_2
251   double one = 1.;
252   double nad = boost::math::float_next(one);
253   std::cout << nad << "\n"  //  1.0000000000000002
254     << nad - one // 2.2204460492503131e-016
255     << std::endl;
256 //] [epsilon_2]
257   }
258   {
259 //[epsilon_3
260   std::cout.precision(std::numeric_limits<double>::max_digits10);
261   double d = 1.;
262   double eps = std::numeric_limits<double>::epsilon();
263   double dpeps = d + eps/2;
264 
265   std::cout << std::showpoint // Ensure all trailing zeros are shown.
266     << dpeps << "\n"       // 1.0000000000000000
267     << eps/2 << std::endl; // 1.1102230246251565e-016
268   std::cout << dpeps - d   // 0.00000000000000000
269     << std::endl;
270 //] [epsilon_3]
271   }
272 
273   {
274     typedef double RealType;
275 //[epsilon_4
276 /*`A tolerance might be defined using this version of epsilon thus:
277 */
278     RealType tolerance = boost::math::tools::epsilon<RealType>() * 2;
279 //] [epsilon_4]
280     (void)tolerance; // warning suppression
281   }
282 
283   {
284     bool b =
285 //[digits10_5
286     -(std::numeric_limits<double>::max)() == std::numeric_limits<double>::lowest();
287 //] [/digits10_5]
288     (void)b;  // warning suppression
289   }
290 
291   {
292 //[denorm_min_1
293   std::cout.precision(std::numeric_limits<double>::max_digits10);
294   if (std::numeric_limits<double>::has_denorm == std::denorm_present)
295   {
296     double d = std::numeric_limits<double>::denorm_min();
297 
298       std::cout << d << std::endl; //  4.9406564584124654e-324
299 
300       int exponent;
301 
302       double significand = frexp(d, &exponent);
303       std::cout << "exponent = " << std::hex << exponent << std::endl; //  fffffbcf
304       std::cout << "significand = " << std::hex << significand << std::endl; // 0.50000000000000000
305   }
306   else
307   {
308     std::cout << "No denormalization. " << std::endl;
309   }
310 //] [denorm_min_1]
311   }
312 
313   {
314 //[round_error_1
315     double round_err = std::numeric_limits<double>::epsilon() // 2.2204460492503131e-016
316                      * std::numeric_limits<double>::round_error(); // 1/2
317     std::cout << round_err << std::endl; // 1.1102230246251565e-016
318 //] [/round_error_1]
319   }
320 
321   {
322     typedef double T;
323 //[tolerance_1
324 /*`For example, if we want a tolerance that might suit about 9 arithmetical operations,
325 say sqrt(9) = 3,  we could define:
326 */
327 
328     T tolerance =  3 * std::numeric_limits<T>::epsilon();
329 
330 /*`This is very widely used in Boost.Math testing
331 with Boost.Test's macro `BOOST_CHECK_CLOSE_FRACTION`
332 */
333 
334     T expected = 1.0;
335     T calculated = 1.0 + std::numeric_limits<T>::epsilon();
336 
337     BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
338 
339 //] [/tolerance_1]
340   }
341 
342 #if !(defined(CI_SUPPRESS_KNOWN_ISSUES) && defined(__GNUC__) && defined(_WIN32))
343   {
344 //[tolerance_2
345 
346   using boost::multiprecision::number;
347   using boost::multiprecision::cpp_dec_float;
348   using boost::multiprecision::et_off;
349 
350   typedef number<cpp_dec_float<50>, et_off > cpp_dec_float_50; // 50 decimal digits.
351 /*`[note that Boost.Test does not yet allow floating-point comparisons with expression templates on,
352 so the default expression template parameter has been replaced by `et_off`.]
353 */
354 
355   cpp_dec_float_50 tolerance =  3 * std::numeric_limits<cpp_dec_float_50>::epsilon();
356   cpp_dec_float_50 expected = boost::math::constants::two_pi<cpp_dec_float_50>();
357   cpp_dec_float_50 calculated = 2 * boost::math::constants::pi<cpp_dec_float_50>();
358 
359   BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
360 
361 //] [/tolerance_2]
362   }
363 
364   {
365 //[tolerance_3
366 
367   using boost::multiprecision::cpp_bin_float_quad;
368 
369   cpp_bin_float_quad tolerance =  3 * std::numeric_limits<cpp_bin_float_quad>::epsilon();
370   cpp_bin_float_quad expected = boost::math::constants::two_pi<cpp_bin_float_quad>();
371   cpp_bin_float_quad calculated = 2 * boost::math::constants::pi<cpp_bin_float_quad>();
372 
373   BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
374 
375 //] [/tolerance_3]
376   }
377 
378   {
379 //[tolerance_4
380 
381   using boost::multiprecision::cpp_bin_float_oct;
382 
383   cpp_bin_float_oct tolerance =  3 * std::numeric_limits<cpp_bin_float_oct>::epsilon();
384   cpp_bin_float_oct expected = boost::math::constants::two_pi<cpp_bin_float_oct>();
385   cpp_bin_float_oct calculated = 2 * boost::math::constants::pi<cpp_bin_float_oct>();
386 
387   BOOST_CHECK_CLOSE_FRACTION(expected, calculated, tolerance);
388 
389 //] [/tolerance_4]
390   }
391 
392   {
393 //[nan_1]
394 
395 /*`NaN can be used with binary multiprecision types like `cpp_bin_float_quad`:
396 */
397   using boost::multiprecision::cpp_bin_float_quad;
398 
399   if (std::numeric_limits<cpp_bin_float_quad>::has_quiet_NaN == true)
400   {
401     cpp_bin_float_quad NaN =  std::numeric_limits<cpp_bin_float_quad>::quiet_NaN();
402     std::cout << "cpp_bin_float_quad NaN is "  << NaN << std::endl; //   cpp_bin_float_quad NaN is nan
403 
404     cpp_bin_float_quad expected = NaN;
405     cpp_bin_float_quad calculated = 2 * NaN;
406     // Comparisons of NaN's always fail:
407     bool b = expected == calculated;
408     std::cout << b << std::endl;
409     BOOST_CHECK_NE(expected, expected);
410     BOOST_CHECK_NE(expected, calculated);
411   }
412   else
413   {
414     std::cout << "Type " << typeid(cpp_bin_float_quad).name() << " does not have NaNs!" << std::endl;
415   }
416 
417 //]  [/nan_1]
418   }
419 
420   {
421 //[facet_1]
422 
423 /*`
424 See [@boost:/libs/math/example/nonfinite_facet_sstream.cpp]
425 and we also need
426 
427   #include <boost/math/special_functions/nonfinite_num_facets.hpp>
428 
429 Then we can equally well use a multiprecision type cpp_bin_float_quad:
430 
431 */
432   using boost::multiprecision::cpp_bin_float_quad;
433 
434   typedef cpp_bin_float_quad T;
435 
436   using boost::math::nonfinite_num_put;
437   using boost::math::nonfinite_num_get;
438   {
439     std::locale old_locale;
440     std::locale tmp_locale(old_locale, new nonfinite_num_put<char>);
441     std::locale new_locale(tmp_locale, new nonfinite_num_get<char>);
442     std::stringstream ss;
443     ss.imbue(new_locale);
444     T inf = std::numeric_limits<T>::infinity();
445     ss << inf; // Write out.
446    BOOST_ASSERT(ss.str() == "inf");
447     T r;
448     ss >> r; // Read back in.
449     BOOST_ASSERT(inf == r); // Confirms that the floating-point values really are identical.
450     std::cout << "infinity output was " << ss.str() << std::endl;
451     std::cout << "infinity input was " << r << std::endl;
452   }
453 
454 /*`
455 ``
456   infinity output was inf
457   infinity input was inf
458 ``
459 Similarly we can do the same with NaN (except that we cannot use `assert` (because any comparisons with NaN always return false).
460 */
461   {
462     std::locale old_locale;
463     std::locale tmp_locale(old_locale, new nonfinite_num_put<char>);
464     std::locale new_locale(tmp_locale, new nonfinite_num_get<char>);
465     std::stringstream ss;
466     ss.imbue(new_locale);
467     T n;
468     T NaN = std::numeric_limits<T>::quiet_NaN();
469     ss << NaN; // Write out.
470     BOOST_ASSERT(ss.str() == "nan");
471     std::cout << "NaN output was " << ss.str() << std::endl;
472     ss >> n; // Read back in.
473     std::cout << "NaN input was " << n << std::endl;
474   }
475 /*`
476 ``
477   NaN output was nan
478   NaN input was nan
479 ``
480 */
481 //]  [/facet_1]
482   }
483 
484 #endif
485 #endif
486 } // BOOST_AUTO_TEST_CASE(test_numeric_limits_snips)
487 
488