ceph/src/boost/libs/math/example/lambert_w_precision_example.cpp

   1 // Copyright Paul A. Bristow 2016, 2018.
   2
   3 // Distributed under the Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt or
   5 //  copy at http ://www.boost.org/LICENSE_1_0.txt).
   6
   7 //! Lambert W examples of controlling precision
   8
   9 // #define BOOST_MATH_INSTRUMENT_LAMBERT_W  // #define only for (much) diagnostic output.
  10
  11 #include <boost/config.hpp> // for BOOST_PLATFORM, BOOST_COMPILER,  BOOST_STDLIB ...
  12 #include <boost/version.hpp>   // for BOOST_MSVC versions.
  13 #include <boost/cstdint.hpp>
  14 #include <boost/exception/exception.hpp>  // boost::exception
  15 #include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
  16 #include <boost/math/policies/policy.hpp>
  17 #include <boost/math/special_functions/next.hpp>  // for float_distance.
  18 #include <boost/math/special_functions/relative_difference.hpp> // for relative and epsilon difference.
  19
  20 // Built-in/fundamental GCC float128 or Intel Quad 128-bit type, if available.
  21 #ifdef BOOST_HAS_FLOAT128
  22 #include <boost/multiprecision/float128.hpp> // Not available for MSVC.
  23 // sets BOOST_MP_USE_FLOAT128 for GCC
  24 using boost::multiprecision::float128;
  25 #endif //# NOT _MSC_VER
  26
  27 #include <boost/multiprecision/cpp_dec_float.hpp> // boost::multiprecision::cpp_dec_float_50
  28 using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type.
  29 using boost::multiprecision::cpp_dec_float_100; // 100 decimal digits type.
  30
  31 #include <boost/multiprecision/cpp_bin_float.hpp>
  32 using boost::multiprecision::cpp_bin_float_double_extended;
  33 using boost::multiprecision::cpp_bin_float_double;
  34 using boost::multiprecision::cpp_bin_float_quad;
  35 // For lambert_w function.
  36 #include <boost/math/special_functions/lambert_w.hpp>
  37 // using boost::math::lambert_w0;
  38 // using boost::math::lambert_wm1;
  39
  40 #include <iostream>
  41 #include <exception>
  42 #include <stdexcept>
  43 #include <string>
  44 #include <limits>  // For std::numeric_limits.
  45
  46 int main()
  47 {
  48   try
  49   {
  50     std::cout << "Lambert W examples of precision control." << std::endl;
  51     std::cout.precision(std::numeric_limits<double>::max_digits10);
  52     std::cout << std::showpoint << std::endl; // Show any trailing zeros.
  53
  54     using boost::math::constants::exp_minus_one;
  55
  56     using boost::math::lambert_w0;
  57     using boost::math::lambert_wm1;
  58
  59     // Error handling policy examples.
  60     using namespace boost::math::policies;
  61     using boost::math::policies::make_policy;
  62     using boost::math::policies::policy;
  63     using boost::math::policies::evaluation_error;
  64     using boost::math::policies::domain_error;
  65     using boost::math::policies::overflow_error;
  66     using boost::math::policies::throw_on_error;
  67
  68 //[lambert_w_precision_reference_w
  69
  70     using boost::multiprecision::cpp_bin_float_50;
  71     using boost::math::float_distance;
  72
  73     cpp_bin_float_50 z("10."); // Note use a decimal digit string, not a double 10.
  74     cpp_bin_float_50 r;
  75     std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10);
  76
  77     r = lambert_w0(z); // Default policy.
  78     std::cout << "lambert_w0(z) cpp_bin_float_50  = " << r << std::endl;
  79     //lambert_w0(z) cpp_bin_float_50  = 1.7455280027406993830743012648753899115352881290809
  80     //       [N[productlog[10], 50]] == 1.7455280027406993830743012648753899115352881290809
  81     std::cout.precision(std::numeric_limits<double>::max_digits10);
  82     std::cout << "lambert_w0(z) static_cast from cpp_bin_float_50  = "
  83       << static_cast<double>(r) << std::endl;
  84     // double lambert_w0(z) static_cast from cpp_bin_float_50  = 1.7455280027406994
  85     // [N[productlog[10], 17]]                                == 1.7455280027406994
  86    std::cout << "bits different from Wolfram = "
  87      << static_cast<int>(float_distance(static_cast<double>(r), 1.7455280027406994))
  88      << std::endl; // 0
  89
  90
  91 //] [/lambert_w_precision_reference_w]
  92
  93 //[lambert_w_precision_0
  94     std::cout.precision(std::numeric_limits<float>::max_digits10); // Show all potentially significant decimal digits,
  95     std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
  96
  97     float x = 10.;
  98     std::cout << "Lambert W (" << x << ") = " << lambert_w0(x) << std::endl;
  99 //] [/lambert_w_precision_0]
 100
 101 /*
 102 //[lambert_w_precision_output_0
 103 Lambert W (10.0000000) = 1.74552800
 104 //] [/lambert_w_precision_output_0]
 105 */
 106   { // Lambert W0 Halley step example
 107 //[lambert_w_precision_1
 108     using boost::math::lambert_w_detail::lambert_w_halley_step;
 109     using boost::math::epsilon_difference;
 110     using boost::math::relative_difference;
 111
 112     std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
 113     std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
 114
 115     cpp_bin_float_50 z50("1.23"); // Note: use a decimal digit string, not a double 1.23!
 116     double z = static_cast<double>(z50);
 117     cpp_bin_float_50 w50;
 118     w50 = lambert_w0(z50);
 119     std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
 120     std::cout << "Reference Lambert W (" << z << ") =\n                                              "
 121       << w50 << std::endl;
 122     std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
 123     double wr = static_cast<double>(w50);
 124     std::cout << "Reference Lambert W (" << z << ") =    " << wr << std::endl;
 125
 126     double w = lambert_w0(z);
 127     std::cout << "Rat/poly Lambert W  (" << z << ")  =   " << lambert_w0(z) << std::endl;
 128     // Add a Halley step to the value obtained from rational polynomial approximation.
 129     double ww = lambert_w_halley_step(lambert_w0(z), z);
 130     std::cout << "Halley Step Lambert W (" << z << ") =  " << lambert_w_halley_step(lambert_w0(z), z) << std::endl;
 131
 132     std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
 133     std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
 134     std::cout << "epsilon difference from Halley step  = " << epsilon_difference(w, ww) << std::endl;
 135     std::cout << "epsilon for float =                    " << std::numeric_limits<double>::epsilon() << std::endl;
 136     std::cout << "bits different from Halley step  =     " << static_cast<int>(float_distance(w, ww)) << std::endl;
 137 //] [/lambert_w_precision_1]
 138
 139
 140 /*
 141 //[lambert_w_precision_output_1
 142   Reference Lambert W (1.2299999999999999822364316059974953532218933105468750) =
 143   0.64520356959320237759035605255334853830173300262666480
 144   Reference Lambert W (1.2300000000000000) =    0.64520356959320235
 145   Rat/poly Lambert W  (1.2300000000000000)  =   0.64520356959320224
 146   Halley Step Lambert W (1.2300000000000000) =  0.64520356959320235
 147   absolute difference from Halley step = -1.1102230246251565e-16
 148   relative difference from Halley step = 1.7207329236029286e-16
 149   epsilon difference from Halley step  = 0.77494921535422934
 150   epsilon for float =                    2.2204460492503131e-16
 151   bits different from Halley step  =     1
 152 //] [/lambert_w_precision_output_1]
 153 */
 154
 155   } // Lambert W0 Halley step example
 156
 157   { // Lambert W-1 Halley step example
 158     //[lambert_w_precision_2
 159     using boost::math::lambert_w_detail::lambert_w_halley_step;
 160     using boost::math::epsilon_difference;
 161     using boost::math::relative_difference;
 162
 163     std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
 164     std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
 165
 166     cpp_bin_float_50 z50("-0.123"); // Note: use a decimal digit string, not a double -1.234!
 167     double z = static_cast<double>(z50);
 168     cpp_bin_float_50 wm1_50;
 169     wm1_50 = lambert_wm1(z50);
 170     std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
 171     std::cout << "Reference Lambert W-1 (" << z << ") =\n                                                  "
 172       << wm1_50 << std::endl;
 173     std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
 174     double wr = static_cast<double>(wm1_50);
 175     std::cout << "Reference Lambert W-1 (" << z << ") =    " << wr << std::endl;
 176
 177     double w = lambert_wm1(z);
 178     std::cout << "Rat/poly Lambert W-1 (" << z << ")  =    " << lambert_wm1(z) << std::endl;
 179     // Add a Halley step to the value obtained from rational polynomial approximation.
 180     double ww = lambert_w_halley_step(lambert_wm1(z), z);
 181     std::cout << "Halley Step Lambert W (" << z << ") =    " << lambert_w_halley_step(lambert_wm1(z), z) << std::endl;
 182
 183     std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
 184     std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
 185     std::cout << "epsilon difference from Halley step  = " << epsilon_difference(w, ww) << std::endl;
 186     std::cout << "epsilon for float =                    " << std::numeric_limits<double>::epsilon() << std::endl;
 187     std::cout << "bits different from Halley step  =     " << static_cast<int>(float_distance(w, ww)) << std::endl;
 188     //] [/lambert_w_precision_2]
 189   }
 190   /*
 191   //[lambert_w_precision_output_2
 192     Reference Lambert W-1 (-0.12299999999999999822364316059974953532218933105468750) =
 193     -3.2849102557740360179084675531714935199110302996513384
 194     Reference Lambert W-1 (-0.12300000000000000) =    -3.2849102557740362
 195     Rat/poly Lambert W-1 (-0.12300000000000000)  =    -3.2849102557740357
 196     Halley Step Lambert W (-0.12300000000000000) =    -3.2849102557740362
 197     absolute difference from Halley step = 4.4408920985006262e-16
 198     relative difference from Halley step = 1.3519066740696092e-16
 199     epsilon difference from Halley step  = 0.60884463935795785
 200     epsilon for float =                    2.2204460492503131e-16
 201     bits different from Halley step  =     -1
 202   //] [/lambert_w_precision_output_2]
 203   */
 204
 205
 206
 207   // Similar example using cpp_bin_float_quad (128-bit floating-point types).
 208
 209   cpp_bin_float_quad zq = 10.;
 210   std::cout << "\nTest evaluation of cpp_bin_float_quad Lambert W(" << zq << ")"
 211     << std::endl;
 212   std::cout << std::setprecision(3) << "std::numeric_limits<cpp_bin_float_quad>::digits = " << std::numeric_limits<cpp_bin_float_quad>::digits << std::endl;
 213   std::cout << std::setprecision(3) << "std::numeric_limits<cpp_bin_float_quad>::epsilon() = " << std::numeric_limits<cpp_bin_float_quad>::epsilon() << std::endl;
 214   std::cout << std::setprecision(3) << "std::numeric_limits<cpp_bin_float_quad>::max_digits10 = " << std::numeric_limits<cpp_bin_float_quad>::max_digits10 << std::endl;
 215   std::cout << std::setprecision(3) << "std::numeric_limits<cpp_bin_float_quad>::digits10 = " << std::numeric_limits<cpp_bin_float_quad>::digits10 << std::endl;
 216   std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::max_digits10);
 217   // All are same precision because double precision first approximation used before Halley.
 218
 219   /*
 220
 221       */
 222
 223   { // Reference value for lambert_w0(10)
 224     cpp_dec_float_50 z("10");
 225     cpp_dec_float_50 r;
 226     std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
 227
 228     r = lambert_w0(z); // Default policy.
 229     std::cout << "lambert_w0(z) cpp_dec_float_50                                   = " << r << std::endl; //  0.56714329040978387299996866221035554975381578718651
 230     std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::max_digits10);
 231
 232     std::cout << "lambert_w0(z) cpp_dec_float_50 cast to quad (max_digits10(" << std::numeric_limits<cpp_bin_float_quad>::max_digits10 <<
 233       " )   = " << static_cast<cpp_bin_float_quad>(r) << std::endl;         // 1.7455280027406993830743012648753899115352881290809
 234     std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::digits10); // 1.745528002740699383074301264875389837
 235     std::cout << "lambert_w0(z) cpp_dec_float_50 cast to quad (digits10(" << std::numeric_limits<cpp_bin_float_quad>::digits10 <<
 236       " )       = " << static_cast<cpp_bin_float_quad>(r) << std::endl;    // 1.74552800274069938307430126487539
 237     std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::digits10 + 1); //
 238
 239     std::cout << "lambert_w0(z) cpp_dec_float_50 cast to quad (digits10(" << std::numeric_limits<cpp_bin_float_quad>::digits10 <<
 240       " )       = " << static_cast<cpp_bin_float_quad>(r) << std::endl;    // 1.74552800274069938307430126487539
 241
 242     // [N[productlog[10], 50]] == 1.7455280027406993830743012648753899115352881290809
 243
 244     // [N[productlog[10], 37]] == 1.745528002740699383074301264875389912
 245     // [N[productlog[10], 34]] == 1.745528002740699383074301264875390
 246     // [N[productlog[10], 33]] == 1.74552800274069938307430126487539
 247
 248     // lambert_w0(z) cpp_dec_float_50 cast to quad = 1.745528002740699383074301264875389837
 249
 250     // lambert_w0(z) cpp_dec_float_50 = 1.7455280027406993830743012648753899115352881290809
 251     // lambert_w0(z) cpp_dec_float_50 cast to quad = 1.745528002740699383074301264875389837
 252     // lambert_w0(z) cpp_dec_float_50 cast to quad = 1.74552800274069938307430126487539
 253     }
 254   }
 255   catch (std::exception& ex)
 256   {
 257     std::cout << ex.what() << std::endl;
 258   }
 259 }  // int main()
 260
 261
 262
 263