ceph/src/boost/libs/math/test/test_lambert_w_integrals_float128.cpp

   1 // Copyright Paul A. Bristow 2016, 2017, 2018.
   2 // Copyright John Maddock 2016.
   3
   4 // Use, modification and distribution are subject to the
   5 // Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt
   7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 // test_lambert_w_integrals.cpp
  10 //! \brief quadrature tests that cover the whole range of the Lambert W0 function.
  11
  12 #include <boost/config.hpp>   // for BOOST_MSVC definition etc.
  13 #include <boost/version.hpp>   // for BOOST_MSVC versions.
  14 #include <climits>
  15
  16 #if defined(BOOST_HAS_FLOAT128) && (LDBL_MANT_DIG > 100)
  17 //
  18 // Mixing __float128 and long double results in:
  19 // error: __float128 and long double cannot be used in the same expression
  20 // whenever long double is a [possibly quasi-] quad precision type.
  21 //
  22 #undef BOOST_HAS_FLOAT128
  23 #endif
  24
  25
  26 #ifdef BOOST_HAS_FLOAT128
  27
  28 // Boost macros
  29 #define BOOST_TEST_MAIN
  30 #define BOOST_LIB_DIAGNOSTIC "on" // Report library file details.
  31 #include <boost/test/included/unit_test.hpp> // Boost.Test
  32 #include <boost/test/tools/floating_point_comparison.hpp>
  33
  34 #include <boost/array.hpp>
  35 #include <boost/type_traits/is_constructible.hpp>
  36
  37 #include <boost/multiprecision/float128.hpp>
  38
  39 #include <boost/math/special_functions/fpclassify.hpp> // isnan, isfinite.
  40 #include <boost/math/special_functions/next.hpp> // float_next, float_prior
  41 using boost::math::float_next;
  42 using boost::math::float_prior;
  43 #include <boost/math/special_functions/ulp.hpp>  // ulp
  44
  45 #include <boost/math/tools/test_value.hpp>  // for create_test_value and macro BOOST_MATH_TEST_VALUE.
  46 #include <boost/math/policies/policy.hpp>
  47 using boost::math::policies::digits2;
  48 using boost::math::policies::digits10;
  49 #include <boost/math/special_functions/lambert_w.hpp> // For Lambert W lambert_w function.
  50 using boost::math::lambert_wm1;
  51 using boost::math::lambert_w0;
  52
  53 #include <limits>
  54 #include <cmath>
  55 #include <typeinfo>
  56 #include <iostream>
  57 #include <type_traits>
  58 #include <exception>
  59
  60 std::string show_versions(void);
  61
  62 // Added code and test for Integral of the Lambert W function: by Nick Thompson.
  63 // https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals
  64
  65 #include <boost/math/constants/constants.hpp> // for integral tests.
  66 #include <boost/math/quadrature/tanh_sinh.hpp> // for integral tests.
  67 #include <boost/math/quadrature/exp_sinh.hpp> // for integral tests.
  68
  69   using boost::math::policies::policy;
  70   using boost::math::policies::make_policy;
  71
  72 // using statements needed for changing error handling policy.
  73 using boost::math::policies::evaluation_error;
  74 using boost::math::policies::domain_error;
  75 using boost::math::policies::overflow_error;
  76 using boost::math::policies::ignore_error;
  77 using boost::math::policies::throw_on_error;
  78
  79 typedef policy<
  80   domain_error<throw_on_error>,
  81   overflow_error<ignore_error>
  82 > no_throw_policy;
  83
  84 // Assumes that function has a throw policy, for example:
  85 //    NOT lambert_w0<T>(1 / (x * x), no_throw_policy());
  86 // Error in function boost::math::quadrature::exp_sinh<double>::integrate:
  87 // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf.
  88 // Please ensure your function evaluates to a finite number of its entire domain.
  89 template <typename T>
  90 T debug_integration_proc(T x)
  91 {
  92    T result; // warning C4701: potentially uninitialized local variable 'result' used
  93   // T result = 0 ; // But result may not be assigned below?
  94   try
  95   {
  96    // Assign function call to result in here...
  97     if (x <= sqrt(boost::math::tools::min_value<T>()) )
  98     {
  99       result = 0;
 100     }
 101     else
 102     {
 103       result = lambert_w0<T>(1 / (x * x));
 104     }
 105    // result = lambert_w0<T>(1 / (x * x), no_throw_policy());  // Bad idea, less helpful diagnostic message is:
 106     // Error in function boost::math::quadrature::exp_sinh<double>::integrate:
 107     // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf.
 108     // Please ensure your function evaluates to a finite number of its entire domain.
 109
 110   } // try
 111   catch (const std::exception& e)
 112   {
 113     std::cout << "Exception " << e.what() << std::endl;
 114     // set breakpoint here:
 115     std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl;
 116     if (!std::isfinite(result))
 117     {
 118       // set breakpoint here:
 119       std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl;
 120     }
 121     if (std::isnan(result))
 122     {
 123       // set breakpoint here:
 124       std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl;
 125     }
 126   } // catch
 127   return result;
 128 } // T debug_integration_proc(T x)
 129
 130 template<class Real>
 131 void test_integrals()
 132 {
 133   // Integral of the Lambert W function:
 134   // https://en.wikipedia.org/wiki/Lambert_W_function
 135   using boost::math::quadrature::tanh_sinh;
 136   using boost::math::quadrature::exp_sinh;
 137   // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html
 138   using std::sqrt;
 139
 140   std::cout << "Integration of type " << typeid(Real).name()  << std::endl;
 141
 142   Real tol = std::numeric_limits<Real>::epsilon();
 143   { //  // Integrate for function lambert_W0(z);
 144     tanh_sinh<Real> ts;
 145     Real a = 0;
 146     Real b = boost::math::constants::e<Real>();
 147     auto f = [](Real z)->Real
 148     {
 149       return lambert_w0<Real>(z);
 150     };
 151     Real z = ts.integrate(f, a, b); // OK without any decltype(f)
 152     BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e<Real>() - 1, tol);
 153   }
 154   {
 155     // Integrate for function lambert_W0(z/(z sqrt(z)).
 156     exp_sinh<Real> es;
 157     auto f = [](Real z)->Real
 158     {
 159       return lambert_w0<Real>(z)/(z * sqrt(z));
 160     };
 161     Real z = es.integrate(f); // OK
 162     BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi<Real>(), tol);
 163   }
 164   {
 165     // Integrate for function lambert_W0(1/z^2).
 166     exp_sinh<Real> es;
 167     //const Real sqrt_min = sqrt(boost::math::tools::min_value<Real>()); // 1.08420217e-19 fo 32-bit float.
 168     // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified
 169     auto f = [](Real z)->Real
 170     {
 171       if (z <= sqrt(boost::math::tools::min_value<Real>()) )
 172       { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter.
 173         return static_cast<Real>(0);
 174       }
 175       else
 176       {
 177         return lambert_w0<Real>(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen.
 178       }
 179     };
 180     Real z = es.integrate(f);
 181     BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi<Real>(), tol);
 182   }
 183 } // template<class Real> void test_integrals()
 184
 185
 186 BOOST_AUTO_TEST_CASE( integrals )
 187 {
 188   std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl;
 189   BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals.");
 190   try
 191   {
 192   // using statements needed to change precision policy.
 193   using boost::math::policies::policy;
 194   using boost::math::policies::make_policy;
 195   using boost::math::policies::precision;
 196   using boost::math::policies::digits2;
 197   using boost::math::policies::digits10;
 198
 199   // using statements needed for changing error handling policy.
 200   using boost::math::policies::evaluation_error;
 201   using boost::math::policies::domain_error;
 202   using boost::math::policies::overflow_error;
 203   using boost::math::policies::ignore_error;
 204   using boost::math::policies::throw_on_error;
 205
 206   /*
 207   typedef policy<
 208     domain_error<throw_on_error>,
 209     overflow_error<ignore_error>
 210   > no_throw_policy;
 211
 212
 213   // Experiment with better diagnostics.
 214   typedef float Real;
 215
 216   Real inf = std::numeric_limits<Real>::infinity();
 217   Real max = (std::numeric_limits<Real>::max)();
 218   std::cout.precision(std::numeric_limits<Real>::max_digits10);
 219   //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308
 220   std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf
 221   std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227
 222   //std::cout << lambert_w0(inf) << std::endl; // inf - will throw.
 223   std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0
 224   std::cout << "lambert_w0(std::numeric_limits<Real>::denorm_min()) = " << lambert_w0(std::numeric_limits<Real>::denorm_min()) << std::endl; // 4.94066e-324
 225   std::cout << "lambert_w0(std::numeric_limits<Real>::min()) = " << lambert_w0((std::numeric_limits<Real>::min)()) << std::endl; // 2.22507e-308
 226
 227   // Approximate the largest lambert_w you can get for type T?
 228   float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits<float>::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162.
 229   std::cout << "w max_f " << max_w_f << std::endl; // 84.2879
 230   Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits<Real>::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162.
 231   std::cout << "w max " << max_w << std::endl; // 703.227
 232
 233   std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; //
 234   std::cout << "test integral 1/z^2" << std::endl;
 235   std::cout << "ULP = " << boost::math::ulp(1., policy<digits2<> >()) << std::endl; // ULP = 2.2204460492503131e-16
 236   std::cout << "ULP = " << boost::math::ulp(1e-10, policy<digits2<> >()) << std::endl; // ULP = 2.2204460492503131e-16
 237   std::cout << "ULP = " << boost::math::ulp(1., policy<digits2<11> >()) << std::endl; // ULP = 2.2204460492503131e-16
 238   std::cout << "epsilon =  " << std::numeric_limits<Real>::epsilon() << std::endl; //
 239   std::cout << "sqrt(max) =  " << sqrt(boost::math::tools::max_value<float>() ) << std::endl; // sqrt(max) =  1.8446742974197924e+19
 240   std::cout << "sqrt(min) =  " << sqrt(boost::math::tools::min_value<float>() ) << std::endl; // sqrt(min) =  1.0842021724855044e-19
 241
 242
 243
 244 // Demo debug version.
 245 Real tol = std::numeric_limits<Real>::epsilon();
 246 Real x;
 247 {
 248   using boost::math::quadrature::exp_sinh;
 249   exp_sinh<Real> es;
 250   // Function to be integrated, lambert_w0(1/z^2).
 251
 252     //auto f = [](Real z)->Real
 253     //{ // Naive - no protection against underflow and subsequent divide by zero.
 254     //  return lambert_w0<Real>(1 / (z * z));
 255     //};
 256     // Diagnostic is:
 257     // Error in function boost::math::lambert_w0<Real>: Expected a finite value but got inf
 258
 259     auto f = [](Real z)->Real
 260     { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things.
 261       return debug_integration_proc(z);
 262     };
 263     // Exception Error in function boost::math::lambert_w0<double>: Expected a finite value but got inf.
 264
 265     // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163.
 266     // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23.
 267     x = es.integrate(f);
 268     std::cout << "es.integrate(f) = " << x << std::endl;
 269     BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi<Real>(), tol);
 270     // root_two_pi<double = 2.506628274631000502
 271   }
 272     */
 273
 274   test_integrals<boost::multiprecision::float128>();
 275   }
 276   catch (std::exception& ex)
 277   {
 278     std::cout << ex.what() << std::endl;
 279   }
 280 }
 281
 282 #else
 283
 284 int main() { return 0; }
 285
 286 #endif