diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp index bf2731dad..7e2a8afc4 100644 --- a/include/boost/math/special_functions/lambert_w.hpp +++ b/include/boost/math/special_functions/lambert_w.hpp @@ -162,6 +162,20 @@ T maybe_reduce_to_double(const T& z, const mpl::false_&) return z; } +template +inline +double must_reduce_to_double(const T& z, const mpl::true_&) +{ + return static_cast(z); // Reduce to double precision. +} + +template +inline +double must_reduce_to_double(const T& z, const mpl::false_&) +{ // try a lexical_cast and hope for the best: + return boost::lexical_cast(z); +} + //! \brief Schroeder method, fifth-order update formula, //! \details See T. Fukushima page 80-81, and //! A. Householder, The Numerical Treatment of a Single Nonlinear Equation, @@ -173,7 +187,7 @@ T maybe_reduce_to_double(const T& z, const mpl::false_&) //! \returns Refined estimate of Lambert w. // Schroeder refinement, called unless NOT required by precision policy. -template +template inline T schroeder_update(const T w, const T y) { @@ -229,7 +243,7 @@ T schroeder_update(const T w, const T y) //! are at least 21 digits precision == max_digits10 for long double. //! Longer decimal digits strings are rationals evaluated using Wolfram. -template +template T lambert_w_singularity_series(const T p) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES @@ -492,10 +506,17 @@ T lambert_w0_small_z(T x, const Policy& pol) { //std::numeric_limits::max_digits10 == 36 ? 3 : // 128-bit long double. typedef boost::mpl::int_ < +#ifndef BOOST_NO_CXX11_NUMERIC_LIMITS std::numeric_limits::max_digits10 <= 9 ? 0 : // for float 32-bit. std::numeric_limits::max_digits10 <= 17 ? 1 : // for double 64-bit. std::numeric_limits::max_digits10 <= 22 ? 2 : // for 80-bit double extended. std::numeric_limits::max_digits10 < 37 ? 4 // for both 128-bit long double (3) and 128-bit quad suffix Q type (4). +#else + ((std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= 24)) ? 0 : // for float 32-bit. + ((std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= 53)) ? 1 : // for double 64-bit. + ((std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= 64)) ? 2 : // for 80-bit double extended. + ((std::numeric_limits::radix == 2) && (std::numeric_limits::digits <= 113)) ? 4 // for both 128-bit long double (3) and 128-bit quad suffix Q type (4). +#endif : 5 // All Generic multiprecision types. > tag_type; // std::cout << "\ntag type = " << tag_type << std::endl; // error C2275: 'tag_type': illegal use of this type as an expression. @@ -953,6 +974,7 @@ T lambert_w0_approx(T z) template inline T get_near_singularity_param(T z) { + BOOST_MATH_STD_USING const T p2 = 2 * (boost::math::constants::e() * z + 1); const T p = sqrt(p2); return p; @@ -963,7 +985,7 @@ inline float get_near_singularity_param(float z) } inline double get_near_singularity_param(double z) { - return static_cast(get_near_singularity_param((long double)z)); + return static_cast(get_near_singularity_param((long double)z)); } // Forward declarations: @@ -1681,7 +1703,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) // Integral types should be promoted to double by user Lambert w functions. // If integral type provided to user function lambert_w0 or lambert_wm1, // then should already have been promoted to double. - BOOST_STATIC_ASSERT_MSG(!std::is_integral::value, + BOOST_STATIC_ASSERT_MSG(!boost::is_integral::value, "Must be floating-point or fixed type (not integer type), for example: lambert_wm1(1.), not lambert_wm1(1)!"); BOOST_MATH_STD_USING // Aid argument dependent lookup (ADL) of abs. @@ -1847,7 +1869,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) using boost::math::policies::digits2; using boost::math::policies::policy; // Compute a 50-bit precision approximate W0 in a double (no Halley refinement). - T double_approx = static_cast(lambert_wm1_imp(static_cast(z), policy >())); + T double_approx(lambert_wm1_imp(must_reduce_to_double(z, boost::is_constructible()), policy >())); #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_WM1_NOT_BUILTIN std::streamsize saved_precision = std::cout.precision(std::numeric_limits::max_digits10); std::cout << "Lambert_wm1 Argument Type " << typeid(T).name() << " approximation double = " << double_approx << std::endl; @@ -1936,13 +1958,16 @@ T lambert_wm1_imp(const T z, const Policy& pol) // and use later Halley refinement for higher precisions. using lambert_w_lookup::halves; using lambert_w_lookup::sqrtwm1s; - lookup_t w = -static_cast(n); // Equation 25, - lookup_t y = static_cast(z * wm1es[n - 1]); // Equation 26, + + typedef typename mpl::if_c::value, lookup_t, T>::type calc_type; + + calc_type w = -static_cast(n); // Equation 25, + calc_type y = static_cast(z * wm1es[n - 1]); // Equation 26, // Perform the bisections fractional bisections for necessary precision. for (int j = 0; j < bisections; ++j) { // Equation 27. - lookup_t wj = w - halves[j]; // Subtract 1/2, 1/4, 1/8 ... - lookup_t yj = y * sqrtwm1s[j]; // Multiply by sqrt(1/e), ... + calc_type wj = w - halves[j]; // Subtract 1/2, 1/4, 1/8 ... + calc_type yj = y * sqrtwm1s[j]; // Multiply by sqrt(1/e), ... if (wj < yj) { w = wj; @@ -1994,7 +2019,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) } // lambert_w0(T z, const Policy& pol) //! Lambert W0 using default policy. - template > + template inline typename tools::promote_args::type lambert_w0(T z) @@ -2005,7 +2030,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) // Work out what precision has been selected, based on the Policy and the number type. // For the default policy version, we want the *default policy* precision for T. - typedef typename policies::precision::type precision_type; + typedef typename policies::precision >::type precision_type; // and then select the correct implementation based on that (not the type T): typedef mpl::int_< (precision_type::value == 0) || (precision_type::value > 53) ? diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 709274e71..e22dad433 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -416,7 +416,17 @@ test-suite special_fun : [ run test_jacobi.cpp pch_light ../../test/build//boost_unit_test_framework ] [ run test_laguerre.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] - [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework ] + [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework ] + [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework : : : BOOST_MATH_TEST_MULTIPRECISION=1 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] : test_lambert_w_multiprecision_1 ] + [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework : : : BOOST_MATH_TEST_MULTIPRECISION=2 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] : test_lambert_w_multiprecision_2 ] + [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework : : : BOOST_MATH_TEST_MULTIPRECISION=3 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] : test_lambert_w_multiprecision_3 ] + [ run test_lambert_w.cpp ../../test/build//boost_unit_test_framework : : : BOOST_MATH_TEST_MULTIPRECISION=4 BOOST_MATH_TEST_FLOAT128 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] : test_lambert_w_multiprecision_4 ] + [ run test_lambert_w_integrals_float128.cpp ../../test/build//boost_unit_test_framework : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] ] + [ run test_lambert_w_integrals_quad.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] + [ run test_lambert_w_integrals_long_double.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] ] + [ run test_lambert_w_integrals_double.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] ] + [ run test_lambert_w_integrals_float.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax sfinae_expr ] ] + [ run test_lambert_w_derivative.cpp ../../test/build//boost_unit_test_framework : : : BOOST_MATH_TEST_MULTIPRECISION [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] [ run test_legendre.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ] ] [ run test_ldouble_simple.cpp ../../test/build//boost_unit_test_framework ] diff --git a/test/test_lambert_w.cpp b/test/test_lambert_w.cpp index 02b13e4dc..6234f6f87 100644 --- a/test/test_lambert_w.cpp +++ b/test/test_lambert_w.cpp @@ -79,7 +79,6 @@ using boost::math::lambert_w0; #include #include #include -#include #include std::string show_versions(void); @@ -298,6 +297,8 @@ void wolfram_test_near_singularity() double tolerance = boost::math::tools::epsilon() * 5; if (std::numeric_limits::digits >= std::numeric_limits::digits) tolerance *= 1e5; + else if (std::numeric_limits::digits * 2 >= std::numeric_limits::digits) + tolerance *= 5e4; double endpoint = -boost::math::constants::exp_minus_one(); for (unsigned i = 0; i < wolfram_test_near_singularity_data.size(); ++i) { @@ -358,9 +359,11 @@ void test_spots(RealType) } RealType tolerance = boost::math::tools::epsilon() * epsilons; // 2 eps as a fraction. std::cout << "Tolerance " << epsilons << " * epsilon == " << tolerance << std::endl; +#ifndef BOOST_NO_CXX11_NUMERIC_LIMITS std::cout << "Precision " << std::numeric_limits::digits10 << " decimal digits, max_digits10 = " << std::numeric_limits ::max_digits10<< std::endl; // std::cout.precision(std::numeric_limits::digits10); std::cout.precision(std::numeric_limits ::max_digits10); +#endif std::cout.setf(std::ios_base::showpoint); // show trailing significant zeros. std::cout << "-exp(-1) = " << -exp_minus_one() << std::endl; @@ -765,7 +768,7 @@ void test_spots(RealType) BOOST_CHECK_THROW(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e30)), std::domain_error); // Too negative - BOOST_CHECK_THROW(lambert_wm1(-0.5), std::domain_error); + BOOST_CHECK_THROW(lambert_wm1(RealType(-0.5)), std::domain_error); // This fails for fixed_point type used for other tests because out of range? //BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, 1.0e6)), @@ -796,7 +799,7 @@ BOOST_AUTO_TEST_CASE( test_types ) boost::unit_test_framework::unit_test_log.set_threshold_level(boost::unit_test_framework::log_messages); BOOST_TEST_MESSAGE("\nTest Lambert W function for several types."); BOOST_TEST_MESSAGE(show_versions()); // Full version of Boost, STL and compiler info. - +#ifndef BOOST_MATH_TEST_MULTIPRECISION // Fundamental built-in types: test_spots(0.0F); // float test_spots(0.0); // double @@ -805,12 +808,18 @@ BOOST_AUTO_TEST_CASE( test_types ) test_spots(0.0L); // long double } - #ifdef BOOST_MATH_TEST_MULTIPRECISION + #else // BOOST_MATH_TEST_MULTIPRECISION // Multiprecision types: +#if BOOST_MATH_TEST_MULTIPRECISION == 1 test_spots(static_cast(0)); +#endif +#if BOOST_MATH_TEST_MULTIPRECISION == 2 test_spots(static_cast(0)); +#endif +#if BOOST_MATH_TEST_MULTIPRECISION == 3 test_spots(static_cast(0)); - #endif // ifdef BOOST_MATH_TEST_MULTIPRECISION +#endif +#endif // ifdef BOOST_MATH_TEST_MULTIPRECISION #ifdef BOOST_MATH_TEST_FLOAT128 std::cout << "\nBOOST_MATH_TEST_FLOAT128 defined for float128 tests." << std::endl; @@ -857,8 +866,9 @@ BOOST_AUTO_TEST_CASE( test_range_of_double_values ) int epsilons = 1; RealType tolerance = boost::math::tools::epsilon() * epsilons; // 2 eps as a fraction. - std::cout << "Tolerance " << epsilons << " * epsilon == " << tolerance << std::endl; + std::cout << "Tolerance " << epsilons << " * epsilon == " << tolerance << std::endl; +#ifndef BOOST_MATH_TEST_MULTIPRECISION BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, 1.0e-6)), BOOST_MATH_TEST_VALUE(RealType, 9.9999900000149999733333854165586669000967020964243e-7), // Output from https://www.wolframalpha.com/input/ N[lambert_w[1e-6],50]) @@ -1020,104 +1030,108 @@ BOOST_AUTO_TEST_CASE( test_range_of_double_values ) BOOST_MATH_TEST_VALUE(RealType, -0.99999997419043196), tolerance);// diff 2.30785e-09 v 2.2204460492503131e-16 -#ifdef BOOST_MATH_TEST_MULTIPRECISION + // z increasingly close to singularity. + BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.36)), + BOOST_MATH_TEST_VALUE(RealType, -0.8060843159708177782855213616209920019974599683466713016), + 2 * tolerance); // -0.806084335 + + BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.365)), + BOOST_MATH_TEST_VALUE(RealType, -0.8798200914159538111724840007674053239388642469453350954), + 5 * tolerance); // Note 5 * tolerance + + BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.3678)), + BOOST_MATH_TEST_VALUE(RealType, -0.9793607149578284774761844434886481686055949229547379368), + 15 * tolerance); // Note 15 * tolerance when this close to singularity. + + // Just using series approximation (Fukushima switch at -0.35, but JM at 0.01 of singularity < -0.3679). + // N[productlog(-0.351), 50] = -0.72398644140937651483634596143951001600417138085814 + // N[productlog(-0.351), 55] = -0.7239864414093765148363459614395100160041713808581379727 + BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.351)), + BOOST_MATH_TEST_VALUE(RealType, -0.72398644140937651483634596143951001600417138085814), + 10 * tolerance); // Note was 2 * tolerance + + // Check value just not using near_singularity series approximation (and using rational polynomial instead). + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.3)), + BOOST_MATH_TEST_VALUE(RealType, -1.7813370234216276119741702815127452608215583564545), + // Output from https://www.wolframalpha.com/input/ + //N[productlog(-1, -0.3), 50] = -1.7813370234216276119741702815127452608215583564545 + tolerance); + + // Using table lookup and schroeder with decreasing z to zero. + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.2)), + BOOST_MATH_TEST_VALUE(RealType, -2.5426413577735264242938061566618482901614749075294), + // N[productlog[-1, -0.2],50] -2.5426413577735264242938061566618482901614749075294 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.1)), + BOOST_MATH_TEST_VALUE(RealType, -3.5771520639572972184093919635119948804017962577931), + //N[productlog(-1, -0.1), 50] = -3.5771520639572972184093919635119948804017962577931 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.001)), + BOOST_MATH_TEST_VALUE(RealType, -9.1180064704027401212583371820468142742704349737639), + // N[productlog(-1, -0.001), 50] = -9.1180064704027401212583371820468142742704349737639 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.000001)), + BOOST_MATH_TEST_VALUE(RealType, -16.626508901372473387706432163984684996461726803805), + // N[productlog(-1, -0.000001), 50] = -16.626508901372473387706432163984684996461726803805 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-6)), + BOOST_MATH_TEST_VALUE(RealType, -16.626508901372473387706432163984684996461726803805), + // N[productlog(-1, -10 ^ -6), 50] = -16.626508901372473387706432163984684996461726803805 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1.0e-26)), + BOOST_MATH_TEST_VALUE(RealType, -64.026509628385889681156090340691637712441162092868), + // Output from https://www.wolframalpha.com/input/ + // N[productlog(-1, -1 * 10^-26 ), 50] = -64.026509628385889681156090340691637712441162092868 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -2e-26)), + BOOST_MATH_TEST_VALUE(RealType, -63.322302839923597803393585145387854867226970485197), + // N[productlog[-1, -2*10^-26],50] = -63.322302839923597803393585145387854867226970485197 + tolerance * 2); + + // Smaller than lookup table, so must use approx and Halley refinements. + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-30)), + BOOST_MATH_TEST_VALUE(RealType, -73.373110313822976797067478758120874529181611813766), + // N[productlog(-1, -10 ^ -30), 50] = -73.373110313822976797067478758120874529181611813766 + tolerance); + + // std::numeric_limits::min +#ifndef BOOST_NO_CXX11_NUMERIC_LIMITS + std::cout.precision(std::numeric_limits::max_digits10); +#endif + std::cout << "(std::numeric_limits::min)() " << (std::numeric_limits::min)() << std::endl; + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -2.2250738585072014e-308)), + BOOST_MATH_TEST_VALUE(RealType, -714.96865723796647086868547560654825435542227693935), + // N[productlog[-1, -2.2250738585072014e-308],50] = -714.96865723796647086868547560654825435542227693935 + tolerance); + + // For z = 0, W = -infinity + if (std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, 0.)), + -std::numeric_limits::infinity()); + } + +#elif BOOST_MATH_TEST_MULTIPRECISION == 2 + // Comparison with Wolfram N[productlog(0,-0.36787944117144228 ), 17] // Using conversion from double to higher precision cpp_bin_float_quad using boost::multiprecision::cpp_bin_float_quad; BOOST_CHECK_CLOSE_FRACTION( // Check float_next(-exp(-1) ) lambert_w0(BOOST_MATH_TEST_VALUE(cpp_bin_float_quad, -0.36787944117144228)), - BOOST_MATH_TEST_VALUE(RealType, -0.99999998496215738), + BOOST_MATH_TEST_VALUE(cpp_bin_float_quad, -0.99999998496215738), tolerance); // OK BOOST_CHECK_CLOSE_FRACTION( // Check float_next(float_next(-exp(-1) )) lambert_w0(BOOST_MATH_TEST_VALUE(cpp_bin_float_quad, -0.36787944117144222)), - BOOST_MATH_TEST_VALUE(RealType, -0.99999997649828679), + BOOST_MATH_TEST_VALUE(cpp_bin_float_quad, -0.99999997649828679), tolerance);// OK #endif - // z increasingly close to singularity. - BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.36)), - BOOST_MATH_TEST_VALUE(RealType, -0.8060843159708177782855213616209920019974599683466713016), - 2 * tolerance); // -0.806084335 - - BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.365)), - BOOST_MATH_TEST_VALUE(RealType, -0.8798200914159538111724840007674053239388642469453350954), - 5 * tolerance); // Note 5 * tolerance - - BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.3678)), - BOOST_MATH_TEST_VALUE(RealType, -0.9793607149578284774761844434886481686055949229547379368), - 15 * tolerance); // Note 15 * tolerance when this close to singularity. - - // Just using series approximation (Fukushima switch at -0.35, but JM at 0.01 of singularity < -0.3679). - // N[productlog(-0.351), 50] = -0.72398644140937651483634596143951001600417138085814 - // N[productlog(-0.351), 55] = -0.7239864414093765148363459614395100160041713808581379727 - BOOST_CHECK_CLOSE_FRACTION(lambert_w0(BOOST_MATH_TEST_VALUE(RealType, -0.351)), - BOOST_MATH_TEST_VALUE(RealType, -0.72398644140937651483634596143951001600417138085814), - 10 * tolerance); // Note was 2 * tolerance - - // Check value just not using near_singularity series approximation (and using rational polynomial instead). - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.3)), - BOOST_MATH_TEST_VALUE(RealType, -1.7813370234216276119741702815127452608215583564545), - // Output from https://www.wolframalpha.com/input/ - //N[productlog(-1, -0.3), 50] = -1.7813370234216276119741702815127452608215583564545 - tolerance); - - // Using table lookup and schroeder with decreasing z to zero. - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.2)), - BOOST_MATH_TEST_VALUE(RealType, -2.5426413577735264242938061566618482901614749075294), - // N[productlog[-1, -0.2],50] -2.5426413577735264242938061566618482901614749075294 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.1)), - BOOST_MATH_TEST_VALUE(RealType, -3.5771520639572972184093919635119948804017962577931), - //N[productlog(-1, -0.1), 50] = -3.5771520639572972184093919635119948804017962577931 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.001)), - BOOST_MATH_TEST_VALUE(RealType, -9.1180064704027401212583371820468142742704349737639), - // N[productlog(-1, -0.001), 50] = -9.1180064704027401212583371820468142742704349737639 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -0.000001)), - BOOST_MATH_TEST_VALUE(RealType, -16.626508901372473387706432163984684996461726803805), - // N[productlog(-1, -0.000001), 50] = -16.626508901372473387706432163984684996461726803805 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-6)), - BOOST_MATH_TEST_VALUE(RealType, -16.626508901372473387706432163984684996461726803805), - // N[productlog(-1, -10 ^ -6), 50] = -16.626508901372473387706432163984684996461726803805 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1.0e-26)), - BOOST_MATH_TEST_VALUE(RealType, -64.026509628385889681156090340691637712441162092868), - // Output from https://www.wolframalpha.com/input/ - // N[productlog(-1, -1 * 10^-26 ), 50] = -64.026509628385889681156090340691637712441162092868 - tolerance); - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -2e-26)), - BOOST_MATH_TEST_VALUE(RealType, -63.322302839923597803393585145387854867226970485197), - // N[productlog[-1, -2*10^-26],50] = -63.322302839923597803393585145387854867226970485197 - tolerance *2); - - // Smaller than lookup table, so must use approx and Halley refinements. - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -1e-30)), - BOOST_MATH_TEST_VALUE(RealType, -73.373110313822976797067478758120874529181611813766), - // N[productlog(-1, -10 ^ -30), 50] = -73.373110313822976797067478758120874529181611813766 - tolerance); - - // std::numeric_limits::min - std::cout.precision(std::numeric_limits::max_digits10); - std::cout << "(std::numeric_limits::min)() " << (std::numeric_limits::min)() << std::endl; - - BOOST_CHECK_CLOSE_FRACTION(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, -2.2250738585072014e-308)), - BOOST_MATH_TEST_VALUE(RealType, -714.96865723796647086868547560654825435542227693935), - // N[productlog[-1, -2.2250738585072014e-308],50] = -714.96865723796647086868547560654825435542227693935 - tolerance); - - // For z = 0, W = -infinity - if (std::numeric_limits::has_infinity) - { - BOOST_CHECK_EQUAL(lambert_wm1(BOOST_MATH_TEST_VALUE(RealType, 0.)), - -std::numeric_limits::infinity()); - } } // BOOST_AUTO_TEST_CASE(test_range_of_double_values) diff --git a/test/test_lambert_w_derivative.cpp b/test/test_lambert_w_derivative.cpp new file mode 100644 index 000000000..aa8946efb --- /dev/null +++ b/test/test_lambert_w_derivative.cpp @@ -0,0 +1,124 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w.cpp +//! \brief Basic sanity tests for Lambert W derivative. + +#ifdef BOOST_MATH_TEST_FLOAT128 +#include // For float_64_t, float128_t. Must be first include! +#endif // #ifdef #ifdef BOOST_MATH_TEST_FLOAT128 +// Needs gnu++17 for BOOST_HAS_FLOAT128 +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include + +#ifdef BOOST_MATH_TEST_MULTIPRECISION +#include // boost::multiprecision::cpp_dec_float_50 +using boost::multiprecision::cpp_dec_float_50; + +#include +using boost::multiprecision::cpp_bin_float_quad; + +#ifdef BOOST_MATH_TEST_FLOAT128 + +#ifdef BOOST_HAS_FLOAT128 +// Including this header below without float128 triggers: +// fatal error C1189: #error: "Sorry compiler is neither GCC, not Intel, don't know how to configure this header." +#include +using boost::multiprecision::float128; +#endif // ifdef BOOST_HAS_FLOAT128 +#endif // #ifdef #ifdef BOOST_MATH_TEST_FLOAT128 + +#endif // #ifdef BOOST_MATH_TEST_MULTIPRECISION + +//#include // If available. + +#include // for real_concept tests. +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include + +std::string show_versions(void); + +BOOST_AUTO_TEST_CASE( Derivatives_of_lambert_w ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W_DERIVATIVES to test 1st derivatives is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W function 1st differentials."); + + using boost::math::constants::exp_minus_one; + using boost::math::lambert_w0_prime; + using boost::math::lambert_wm1_prime; + + // Derivatives + // https://www.wolframalpha.com/input/?i=derivative+of+productlog(0,+x) + // d/dx(W_0(x)) = W(x)/(x W(x) + x) + // https://www.wolframalpha.com/input/?i=derivative+of+productlog(-1,+x) + // d/dx(W_(-1)(x)) = (W_(-1)(x))/(x W_(-1)(x) + x) + + // 55 decimal digit values added to allow future testing using multiprecision. + + typedef double RealType; + + int epsilons = 1; + RealType tolerance = boost::math::tools::epsilon() * epsilons; // 2 eps as a fraction. + + // derivative of productlog(-1, x) at x = -0.1 == -13.8803 + // (derivative of productlog(-1, x) ) at x = N[-0.1, 55] - but the result disappears! + // (derivative of N[productlog(-1, x), 55] ) at x = N[-0.1, 55] + + // W0 branch + BOOST_CHECK_CLOSE_FRACTION(lambert_w0_prime(BOOST_MATH_TEST_VALUE(RealType, -0.2)), + // BOOST_MATH_TEST_VALUE(RealType, 1.7491967609218355), + BOOST_MATH_TEST_VALUE(RealType, 1.7491967609218358355273514903396335693828167746571404), + tolerance); // 1.7491967609218358355273514903396335693828167746571404 + + BOOST_CHECK_CLOSE_FRACTION(lambert_w0_prime(BOOST_MATH_TEST_VALUE(RealType, 10.)), + BOOST_MATH_TEST_VALUE(RealType, 0.063577133469345105142021311010780887641928338458371618), + tolerance); + +// W-1 branch + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1_prime(BOOST_MATH_TEST_VALUE(RealType, -0.1)), + BOOST_MATH_TEST_VALUE(RealType, -13.880252213229780748699361486619519025203815492277715), + tolerance); + // Lambert W_prime -13.880252213229780748699361486619519025203815492277715, double -13.880252213229781 + + BOOST_CHECK_CLOSE_FRACTION(lambert_wm1_prime(BOOST_MATH_TEST_VALUE(RealType, -0.2)), + BOOST_MATH_TEST_VALUE(RealType, -8.2411940564179044961885598641955579728547896392013239), + tolerance); + // Lambert W_prime -8.2411940564179044961885598641955579728547896392013239, double -8.2411940564179051 + + // Lambert W_prime 0.063577133469345105142021311010780887641928338458371618, double 0.063577133469345098 +}; // BOOST_AUTO_TEST_CASE("Derivatives of lambert_w") + + diff --git a/test/test_lambert_w_integrals_double.cpp b/test/test_lambert_w_integrals_double.cpp new file mode 100644 index 000000000..c286d4b94 --- /dev/null +++ b/test/test_lambert_w_integrals_double.cpp @@ -0,0 +1,266 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w_integrals.cpp +//! \brief quadrature tests that cover the whole range of the Lambert W0 function. + +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include +#include + +std::string show_versions(void); + +// Added code and test for Integral of the Lambert W function: by Nick Thompson. +// https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals + +#include // for integral tests. +#include // for integral tests. +#include // for integral tests. + + using boost::math::policies::policy; + using boost::math::policies::make_policy; + +// using statements needed for changing error handling policy. +using boost::math::policies::evaluation_error; +using boost::math::policies::domain_error; +using boost::math::policies::overflow_error; +using boost::math::policies::ignore_error; +using boost::math::policies::throw_on_error; + +typedef policy< + domain_error, + overflow_error +> no_throw_policy; + +// Assumes that function has a throw policy, for example: +// NOT lambert_w0(1 / (x * x), no_throw_policy()); +// Error in function boost::math::quadrature::exp_sinh::integrate: +// The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. +// Please ensure your function evaluates to a finite number of its entire domain. +template +T debug_integration_proc(T x) +{ + T result; // warning C4701: potentially uninitialized local variable 'result' used + // T result = 0 ; // But result may not be assigned below? + try + { + // Assign function call to result in here... + if (x <= sqrt(boost::math::tools::min_value()) ) + { + result = 0; + } + else + { + result = lambert_w0(1 / (x * x)); + } + // result = lambert_w0(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is: + // Error in function boost::math::quadrature::exp_sinh::integrate: + // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. + // Please ensure your function evaluates to a finite number of its entire domain. + + } // try + catch (const std::exception& e) + { + std::cout << "Exception " << e.what() << std::endl; + // set breakpoint here: + std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl; + if (!std::isfinite(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + if (std::isnan(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + } // catch + return result; +} // T debug_integration_proc(T x) + +template +void test_integrals() +{ + // Integral of the Lambert W function: + // https://en.wikipedia.org/wiki/Lambert_W_function + using boost::math::quadrature::tanh_sinh; + using boost::math::quadrature::exp_sinh; + // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html + using std::sqrt; + + std::cout << "Integration of type " << typeid(Real).name() << std::endl; + + Real tol = std::numeric_limits::epsilon(); + { // // Integrate for function lambert_W0(z); + tanh_sinh ts; + Real a = 0; + Real b = boost::math::constants::e(); + auto f = [](Real z)->Real + { + return lambert_w0(z); + }; + Real z = ts.integrate(f, a, b); // OK without any decltype(f) + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e() - 1, tol); + } + { + // Integrate for function lambert_W0(z/(z sqrt(z)). + exp_sinh es; + auto f = [](Real z)->Real + { + return lambert_w0(z)/(z * sqrt(z)); + }; + Real z = es.integrate(f); // OK + BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi(), tol); + } + { + // Integrate for function lambert_W0(1/z^2). + exp_sinh es; + //const Real sqrt_min = sqrt(boost::math::tools::min_value()); // 1.08420217e-19 fo 32-bit float. + // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified + auto f = [](Real z)->Real + { + if (z <= sqrt(boost::math::tools::min_value()) ) + { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter. + return static_cast(0); + } + else + { + return lambert_w0(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen. + } + }; + Real z = es.integrate(f); + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi(), tol); + } +} // template void test_integrals() + + +BOOST_AUTO_TEST_CASE( integrals ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals."); + try + { + // using statements needed to change precision policy. + using boost::math::policies::policy; + using boost::math::policies::make_policy; + using boost::math::policies::precision; + using boost::math::policies::digits2; + using boost::math::policies::digits10; + + // using statements needed for changing error handling policy. + using boost::math::policies::evaluation_error; + using boost::math::policies::domain_error; + using boost::math::policies::overflow_error; + using boost::math::policies::ignore_error; + using boost::math::policies::throw_on_error; + + typedef policy< + domain_error, + overflow_error + > no_throw_policy; + + /* + // Experiment with better diagnostics. + typedef float Real; + + Real inf = std::numeric_limits::infinity(); + Real max = (std::numeric_limits::max)(); + std::cout.precision(std::numeric_limits::max_digits10); + //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308 + std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf + std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227 + //std::cout << lambert_w0(inf) << std::endl; // inf - will throw. + std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0 + std::cout << "lambert_w0(std::numeric_limits::denorm_min()) = " << lambert_w0(std::numeric_limits::denorm_min()) << std::endl; // 4.94066e-324 + std::cout << "lambert_w0(std::numeric_limits::min()) = " << lambert_w0((std::numeric_limits::min)()) << std::endl; // 2.22507e-308 + + // Approximate the largest lambert_w you can get for type T? + float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max_f " << max_w_f << std::endl; // 84.2879 + Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max " << max_w << std::endl; // 703.227 + + std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; // + std::cout << "test integral 1/z^2" << std::endl; + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1e-10, policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "epsilon = " << std::numeric_limits::epsilon() << std::endl; // + std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19 + std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19 + + + +// Demo debug version. +Real tol = std::numeric_limits::epsilon(); +Real x; +{ + using boost::math::quadrature::exp_sinh; + exp_sinh es; + // Function to be integrated, lambert_w0(1/z^2). + + //auto f = [](Real z)->Real + //{ // Naive - no protection against underflow and subsequent divide by zero. + // return lambert_w0(1 / (z * z)); + //}; + // Diagnostic is: + // Error in function boost::math::lambert_w0: Expected a finite value but got inf + + auto f = [](Real z)->Real + { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things. + return debug_integration_proc(z); + }; + // Exception Error in function boost::math::lambert_w0: Expected a finite value but got inf. + + // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163. + // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23. + x = es.integrate(f); + std::cout << "es.integrate(f) = " << x << std::endl; + BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi(), tol); + // root_two_pi(); + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } +} + diff --git a/test/test_lambert_w_integrals_float.cpp b/test/test_lambert_w_integrals_float.cpp new file mode 100644 index 000000000..c7a80a1f7 --- /dev/null +++ b/test/test_lambert_w_integrals_float.cpp @@ -0,0 +1,266 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w_integrals.cpp +//! \brief quadrature tests that cover the whole range of the Lambert W0 function. + +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include +#include + +std::string show_versions(void); + +// Added code and test for Integral of the Lambert W function: by Nick Thompson. +// https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals + +#include // for integral tests. +#include // for integral tests. +#include // for integral tests. + + using boost::math::policies::policy; + using boost::math::policies::make_policy; + +// using statements needed for changing error handling policy. +using boost::math::policies::evaluation_error; +using boost::math::policies::domain_error; +using boost::math::policies::overflow_error; +using boost::math::policies::ignore_error; +using boost::math::policies::throw_on_error; + +typedef policy< + domain_error, + overflow_error +> no_throw_policy; + +// Assumes that function has a throw policy, for example: +// NOT lambert_w0(1 / (x * x), no_throw_policy()); +// Error in function boost::math::quadrature::exp_sinh::integrate: +// The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. +// Please ensure your function evaluates to a finite number of its entire domain. +template +T debug_integration_proc(T x) +{ + T result; // warning C4701: potentially uninitialized local variable 'result' used + // T result = 0 ; // But result may not be assigned below? + try + { + // Assign function call to result in here... + if (x <= sqrt(boost::math::tools::min_value()) ) + { + result = 0; + } + else + { + result = lambert_w0(1 / (x * x)); + } + // result = lambert_w0(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is: + // Error in function boost::math::quadrature::exp_sinh::integrate: + // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. + // Please ensure your function evaluates to a finite number of its entire domain. + + } // try + catch (const std::exception& e) + { + std::cout << "Exception " << e.what() << std::endl; + // set breakpoint here: + std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl; + if (!std::isfinite(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + if (std::isnan(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + } // catch + return result; +} // T debug_integration_proc(T x) + +template +void test_integrals() +{ + // Integral of the Lambert W function: + // https://en.wikipedia.org/wiki/Lambert_W_function + using boost::math::quadrature::tanh_sinh; + using boost::math::quadrature::exp_sinh; + // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html + using std::sqrt; + + std::cout << "Integration of type " << typeid(Real).name() << std::endl; + + Real tol = std::numeric_limits::epsilon(); + { // // Integrate for function lambert_W0(z); + tanh_sinh ts; + Real a = 0; + Real b = boost::math::constants::e(); + auto f = [](Real z)->Real + { + return lambert_w0(z); + }; + Real z = ts.integrate(f, a, b); // OK without any decltype(f) + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e() - 1, tol); + } + { + // Integrate for function lambert_W0(z/(z sqrt(z)). + exp_sinh es; + auto f = [](Real z)->Real + { + return lambert_w0(z)/(z * sqrt(z)); + }; + Real z = es.integrate(f); // OK + BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi(), tol); + } + { + // Integrate for function lambert_W0(1/z^2). + exp_sinh es; + //const Real sqrt_min = sqrt(boost::math::tools::min_value()); // 1.08420217e-19 fo 32-bit float. + // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified + auto f = [](Real z)->Real + { + if (z <= sqrt(boost::math::tools::min_value()) ) + { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter. + return static_cast(0); + } + else + { + return lambert_w0(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen. + } + }; + Real z = es.integrate(f); + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi(), tol); + } +} // template void test_integrals() + + +BOOST_AUTO_TEST_CASE( integrals ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals."); + try + { + // using statements needed to change precision policy. + using boost::math::policies::policy; + using boost::math::policies::make_policy; + using boost::math::policies::precision; + using boost::math::policies::digits2; + using boost::math::policies::digits10; + + // using statements needed for changing error handling policy. + using boost::math::policies::evaluation_error; + using boost::math::policies::domain_error; + using boost::math::policies::overflow_error; + using boost::math::policies::ignore_error; + using boost::math::policies::throw_on_error; + + typedef policy< + domain_error, + overflow_error + > no_throw_policy; + + /* + // Experiment with better diagnostics. + typedef float Real; + + Real inf = std::numeric_limits::infinity(); + Real max = (std::numeric_limits::max)(); + std::cout.precision(std::numeric_limits::max_digits10); + //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308 + std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf + std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227 + //std::cout << lambert_w0(inf) << std::endl; // inf - will throw. + std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0 + std::cout << "lambert_w0(std::numeric_limits::denorm_min()) = " << lambert_w0(std::numeric_limits::denorm_min()) << std::endl; // 4.94066e-324 + std::cout << "lambert_w0(std::numeric_limits::min()) = " << lambert_w0((std::numeric_limits::min)()) << std::endl; // 2.22507e-308 + + // Approximate the largest lambert_w you can get for type T? + float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max_f " << max_w_f << std::endl; // 84.2879 + Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max " << max_w << std::endl; // 703.227 + + std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; // + std::cout << "test integral 1/z^2" << std::endl; + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1e-10, policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "epsilon = " << std::numeric_limits::epsilon() << std::endl; // + std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19 + std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19 + + + +// Demo debug version. +Real tol = std::numeric_limits::epsilon(); +Real x; +{ + using boost::math::quadrature::exp_sinh; + exp_sinh es; + // Function to be integrated, lambert_w0(1/z^2). + + //auto f = [](Real z)->Real + //{ // Naive - no protection against underflow and subsequent divide by zero. + // return lambert_w0(1 / (z * z)); + //}; + // Diagnostic is: + // Error in function boost::math::lambert_w0: Expected a finite value but got inf + + auto f = [](Real z)->Real + { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things. + return debug_integration_proc(z); + }; + // Exception Error in function boost::math::lambert_w0: Expected a finite value but got inf. + + // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163. + // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23. + x = es.integrate(f); + std::cout << "es.integrate(f) = " << x << std::endl; + BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi(), tol); + // root_two_pi(); + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } +} + diff --git a/test/test_lambert_w_integrals_float128.cpp b/test/test_lambert_w_integrals_float128.cpp new file mode 100644 index 000000000..9db4fc199 --- /dev/null +++ b/test/test_lambert_w_integrals_float128.cpp @@ -0,0 +1,269 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w_integrals.cpp +//! \brief quadrature tests that cover the whole range of the Lambert W0 function. + +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include + +#include + +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include +#include + +std::string show_versions(void); + +// Added code and test for Integral of the Lambert W function: by Nick Thompson. +// https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals + +#include // for integral tests. +#include // for integral tests. +#include // for integral tests. + + using boost::math::policies::policy; + using boost::math::policies::make_policy; + +// using statements needed for changing error handling policy. +using boost::math::policies::evaluation_error; +using boost::math::policies::domain_error; +using boost::math::policies::overflow_error; +using boost::math::policies::ignore_error; +using boost::math::policies::throw_on_error; + +typedef policy< + domain_error, + overflow_error +> no_throw_policy; + +// Assumes that function has a throw policy, for example: +// NOT lambert_w0(1 / (x * x), no_throw_policy()); +// Error in function boost::math::quadrature::exp_sinh::integrate: +// The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. +// Please ensure your function evaluates to a finite number of its entire domain. +template +T debug_integration_proc(T x) +{ + T result; // warning C4701: potentially uninitialized local variable 'result' used + // T result = 0 ; // But result may not be assigned below? + try + { + // Assign function call to result in here... + if (x <= sqrt(boost::math::tools::min_value()) ) + { + result = 0; + } + else + { + result = lambert_w0(1 / (x * x)); + } + // result = lambert_w0(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is: + // Error in function boost::math::quadrature::exp_sinh::integrate: + // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. + // Please ensure your function evaluates to a finite number of its entire domain. + + } // try + catch (const std::exception& e) + { + std::cout << "Exception " << e.what() << std::endl; + // set breakpoint here: + std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl; + if (!std::isfinite(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + if (std::isnan(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + } // catch + return result; +} // T debug_integration_proc(T x) + +template +void test_integrals() +{ + // Integral of the Lambert W function: + // https://en.wikipedia.org/wiki/Lambert_W_function + using boost::math::quadrature::tanh_sinh; + using boost::math::quadrature::exp_sinh; + // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html + using std::sqrt; + + std::cout << "Integration of type " << typeid(Real).name() << std::endl; + + Real tol = std::numeric_limits::epsilon(); + { // // Integrate for function lambert_W0(z); + tanh_sinh ts; + Real a = 0; + Real b = boost::math::constants::e(); + auto f = [](Real z)->Real + { + return lambert_w0(z); + }; + Real z = ts.integrate(f, a, b); // OK without any decltype(f) + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e() - 1, tol); + } + { + // Integrate for function lambert_W0(z/(z sqrt(z)). + exp_sinh es; + auto f = [](Real z)->Real + { + return lambert_w0(z)/(z * sqrt(z)); + }; + Real z = es.integrate(f); // OK + BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi(), tol); + } + { + // Integrate for function lambert_W0(1/z^2). + exp_sinh es; + //const Real sqrt_min = sqrt(boost::math::tools::min_value()); // 1.08420217e-19 fo 32-bit float. + // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified + auto f = [](Real z)->Real + { + if (z <= sqrt(boost::math::tools::min_value()) ) + { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter. + return static_cast(0); + } + else + { + return lambert_w0(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen. + } + }; + Real z = es.integrate(f); + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi(), tol); + } +} // template void test_integrals() + + +BOOST_AUTO_TEST_CASE( integrals ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals."); + try + { + // using statements needed to change precision policy. + using boost::math::policies::policy; + using boost::math::policies::make_policy; + using boost::math::policies::precision; + using boost::math::policies::digits2; + using boost::math::policies::digits10; + + // using statements needed for changing error handling policy. + using boost::math::policies::evaluation_error; + using boost::math::policies::domain_error; + using boost::math::policies::overflow_error; + using boost::math::policies::ignore_error; + using boost::math::policies::throw_on_error; + + typedef policy< + domain_error, + overflow_error + > no_throw_policy; + + /* + // Experiment with better diagnostics. + typedef float Real; + + Real inf = std::numeric_limits::infinity(); + Real max = (std::numeric_limits::max)(); + std::cout.precision(std::numeric_limits::max_digits10); + //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308 + std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf + std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227 + //std::cout << lambert_w0(inf) << std::endl; // inf - will throw. + std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0 + std::cout << "lambert_w0(std::numeric_limits::denorm_min()) = " << lambert_w0(std::numeric_limits::denorm_min()) << std::endl; // 4.94066e-324 + std::cout << "lambert_w0(std::numeric_limits::min()) = " << lambert_w0((std::numeric_limits::min)()) << std::endl; // 2.22507e-308 + + // Approximate the largest lambert_w you can get for type T? + float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max_f " << max_w_f << std::endl; // 84.2879 + Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max " << max_w << std::endl; // 703.227 + + std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; // + std::cout << "test integral 1/z^2" << std::endl; + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1e-10, policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "epsilon = " << std::numeric_limits::epsilon() << std::endl; // + std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19 + std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19 + + + +// Demo debug version. +Real tol = std::numeric_limits::epsilon(); +Real x; +{ + using boost::math::quadrature::exp_sinh; + exp_sinh es; + // Function to be integrated, lambert_w0(1/z^2). + + //auto f = [](Real z)->Real + //{ // Naive - no protection against underflow and subsequent divide by zero. + // return lambert_w0(1 / (z * z)); + //}; + // Diagnostic is: + // Error in function boost::math::lambert_w0: Expected a finite value but got inf + + auto f = [](Real z)->Real + { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things. + return debug_integration_proc(z); + }; + // Exception Error in function boost::math::lambert_w0: Expected a finite value but got inf. + + // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163. + // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23. + x = es.integrate(f); + std::cout << "es.integrate(f) = " << x << std::endl; + BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi(), tol); + // root_two_pi(); + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } +} + diff --git a/test/test_lambert_w_integrals_long_double.cpp b/test/test_lambert_w_integrals_long_double.cpp new file mode 100644 index 000000000..180fb330f --- /dev/null +++ b/test/test_lambert_w_integrals_long_double.cpp @@ -0,0 +1,266 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w_integrals.cpp +//! \brief quadrature tests that cover the whole range of the Lambert W0 function. + +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include +#include + +std::string show_versions(void); + +// Added code and test for Integral of the Lambert W function: by Nick Thompson. +// https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals + +#include // for integral tests. +#include // for integral tests. +#include // for integral tests. + + using boost::math::policies::policy; + using boost::math::policies::make_policy; + +// using statements needed for changing error handling policy. +using boost::math::policies::evaluation_error; +using boost::math::policies::domain_error; +using boost::math::policies::overflow_error; +using boost::math::policies::ignore_error; +using boost::math::policies::throw_on_error; + +typedef policy< + domain_error, + overflow_error +> no_throw_policy; + +// Assumes that function has a throw policy, for example: +// NOT lambert_w0(1 / (x * x), no_throw_policy()); +// Error in function boost::math::quadrature::exp_sinh::integrate: +// The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. +// Please ensure your function evaluates to a finite number of its entire domain. +template +T debug_integration_proc(T x) +{ + T result; // warning C4701: potentially uninitialized local variable 'result' used + // T result = 0 ; // But result may not be assigned below? + try + { + // Assign function call to result in here... + if (x <= sqrt(boost::math::tools::min_value()) ) + { + result = 0; + } + else + { + result = lambert_w0(1 / (x * x)); + } + // result = lambert_w0(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is: + // Error in function boost::math::quadrature::exp_sinh::integrate: + // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. + // Please ensure your function evaluates to a finite number of its entire domain. + + } // try + catch (const std::exception& e) + { + std::cout << "Exception " << e.what() << std::endl; + // set breakpoint here: + std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl; + if (!std::isfinite(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + if (std::isnan(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + } // catch + return result; +} // T debug_integration_proc(T x) + +template +void test_integrals() +{ + // Integral of the Lambert W function: + // https://en.wikipedia.org/wiki/Lambert_W_function + using boost::math::quadrature::tanh_sinh; + using boost::math::quadrature::exp_sinh; + // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html + using std::sqrt; + + std::cout << "Integration of type " << typeid(Real).name() << std::endl; + + Real tol = std::numeric_limits::epsilon(); + { // // Integrate for function lambert_W0(z); + tanh_sinh ts; + Real a = 0; + Real b = boost::math::constants::e(); + auto f = [](Real z)->Real + { + return lambert_w0(z); + }; + Real z = ts.integrate(f, a, b); // OK without any decltype(f) + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e() - 1, tol); + } + { + // Integrate for function lambert_W0(z/(z sqrt(z)). + exp_sinh es; + auto f = [](Real z)->Real + { + return lambert_w0(z)/(z * sqrt(z)); + }; + Real z = es.integrate(f); // OK + BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi(), tol); + } + { + // Integrate for function lambert_W0(1/z^2). + exp_sinh es; + //const Real sqrt_min = sqrt(boost::math::tools::min_value()); // 1.08420217e-19 fo 32-bit float. + // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified + auto f = [](Real z)->Real + { + if (z <= sqrt(boost::math::tools::min_value()) ) + { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter. + return static_cast(0); + } + else + { + return lambert_w0(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen. + } + }; + Real z = es.integrate(f); + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi(), tol); + } +} // template void test_integrals() + + +BOOST_AUTO_TEST_CASE( integrals ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals."); + try + { + // using statements needed to change precision policy. + using boost::math::policies::policy; + using boost::math::policies::make_policy; + using boost::math::policies::precision; + using boost::math::policies::digits2; + using boost::math::policies::digits10; + + // using statements needed for changing error handling policy. + using boost::math::policies::evaluation_error; + using boost::math::policies::domain_error; + using boost::math::policies::overflow_error; + using boost::math::policies::ignore_error; + using boost::math::policies::throw_on_error; + + typedef policy< + domain_error, + overflow_error + > no_throw_policy; + + /* + // Experiment with better diagnostics. + typedef float Real; + + Real inf = std::numeric_limits::infinity(); + Real max = (std::numeric_limits::max)(); + std::cout.precision(std::numeric_limits::max_digits10); + //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308 + std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf + std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227 + //std::cout << lambert_w0(inf) << std::endl; // inf - will throw. + std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0 + std::cout << "lambert_w0(std::numeric_limits::denorm_min()) = " << lambert_w0(std::numeric_limits::denorm_min()) << std::endl; // 4.94066e-324 + std::cout << "lambert_w0(std::numeric_limits::min()) = " << lambert_w0((std::numeric_limits::min)()) << std::endl; // 2.22507e-308 + + // Approximate the largest lambert_w you can get for type T? + float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max_f " << max_w_f << std::endl; // 84.2879 + Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max " << max_w << std::endl; // 703.227 + + std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; // + std::cout << "test integral 1/z^2" << std::endl; + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1e-10, policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "epsilon = " << std::numeric_limits::epsilon() << std::endl; // + std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19 + std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19 + + + +// Demo debug version. +Real tol = std::numeric_limits::epsilon(); +Real x; +{ + using boost::math::quadrature::exp_sinh; + exp_sinh es; + // Function to be integrated, lambert_w0(1/z^2). + + //auto f = [](Real z)->Real + //{ // Naive - no protection against underflow and subsequent divide by zero. + // return lambert_w0(1 / (z * z)); + //}; + // Diagnostic is: + // Error in function boost::math::lambert_w0: Expected a finite value but got inf + + auto f = [](Real z)->Real + { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things. + return debug_integration_proc(z); + }; + // Exception Error in function boost::math::lambert_w0: Expected a finite value but got inf. + + // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163. + // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23. + x = es.integrate(f); + std::cout << "es.integrate(f) = " << x << std::endl; + BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi(), tol); + // root_two_pi(); + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } +} + diff --git a/test/test_lambert_w_integrals_quad.cpp b/test/test_lambert_w_integrals_quad.cpp new file mode 100644 index 000000000..907b43cff --- /dev/null +++ b/test/test_lambert_w_integrals_quad.cpp @@ -0,0 +1,269 @@ +// Copyright Paul A. Bristow 2016, 2017, 2018. +// Copyright John Maddock 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +// test_lambert_w_integrals.cpp +//! \brief quadrature tests that cover the whole range of the Lambert W0 function. + +#include // for BOOST_MSVC definition etc. +#include // for BOOST_MSVC versions. + +// Boost macros +#define BOOST_TEST_MAIN +#define BOOST_LIB_DIAGNOSTIC "on" // Report library file details. +#include // Boost.Test +// #include // Boost.Test +#include + +#include +#include +#include +#include +using boost::multiprecision::cpp_bin_float_quad; + +#include // isnan, ifinite. +#include // float_next, float_prior +using boost::math::float_next; +using boost::math::float_prior; +#include // ulp + +#include // for create_test_value and macro BOOST_MATH_TEST_VALUE. +#include +using boost::math::policies::digits2; +using boost::math::policies::digits10; +#include // For Lambert W lambert_w function. +using boost::math::lambert_wm1; +using boost::math::lambert_w0; + +#include +#include +#include +#include +#include +#include + +std::string show_versions(void); + +// Added code and test for Integral of the Lambert W function: by Nick Thompson. +// https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals + +#include // for integral tests. +#include // for integral tests. +#include // for integral tests. + + using boost::math::policies::policy; + using boost::math::policies::make_policy; + +// using statements needed for changing error handling policy. +using boost::math::policies::evaluation_error; +using boost::math::policies::domain_error; +using boost::math::policies::overflow_error; +using boost::math::policies::ignore_error; +using boost::math::policies::throw_on_error; + +typedef policy< + domain_error, + overflow_error +> no_throw_policy; + +// Assumes that function has a throw policy, for example: +// NOT lambert_w0(1 / (x * x), no_throw_policy()); +// Error in function boost::math::quadrature::exp_sinh::integrate: +// The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. +// Please ensure your function evaluates to a finite number of its entire domain. +template +T debug_integration_proc(T x) +{ + T result; // warning C4701: potentially uninitialized local variable 'result' used + // T result = 0 ; // But result may not be assigned below? + try + { + // Assign function call to result in here... + if (x <= sqrt(boost::math::tools::min_value()) ) + { + result = 0; + } + else + { + result = lambert_w0(1 / (x * x)); + } + // result = lambert_w0(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is: + // Error in function boost::math::quadrature::exp_sinh::integrate: + // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf. + // Please ensure your function evaluates to a finite number of its entire domain. + + } // try + catch (const std::exception& e) + { + std::cout << "Exception " << e.what() << std::endl; + // set breakpoint here: + std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl; + if (!std::isfinite(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + if (std::isnan(result)) + { + // set breakpoint here: + std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl; + } + } // catch + return result; +} // T debug_integration_proc(T x) + +template +void test_integrals() +{ + // Integral of the Lambert W function: + // https://en.wikipedia.org/wiki/Lambert_W_function + using boost::math::quadrature::tanh_sinh; + using boost::math::quadrature::exp_sinh; + // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html + using std::sqrt; + + std::cout << "Integration of type " << typeid(Real).name() << std::endl; + + Real tol = std::numeric_limits::epsilon(); + { // // Integrate for function lambert_W0(z); + tanh_sinh ts; + Real a = 0; + Real b = boost::math::constants::e(); + auto f = [](Real z)->Real + { + return lambert_w0(z); + }; + Real z = ts.integrate(f, a, b); // OK without any decltype(f) + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e() - 1, tol); + } + { + // Integrate for function lambert_W0(z/(z sqrt(z)). + exp_sinh es; + auto f = [](Real z)->Real + { + return lambert_w0(z)/(z * sqrt(z)); + }; + Real z = es.integrate(f); // OK + BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi(), tol); + } + { + // Integrate for function lambert_W0(1/z^2). + exp_sinh es; + //const Real sqrt_min = sqrt(boost::math::tools::min_value()); // 1.08420217e-19 fo 32-bit float. + // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified + auto f = [](Real z)->Real + { + if (z <= sqrt(boost::math::tools::min_value()) ) + { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter. + return static_cast(0); + } + else + { + return lambert_w0(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen. + } + }; + Real z = es.integrate(f); + BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi(), tol); + } +} // template void test_integrals() + + +BOOST_AUTO_TEST_CASE( integrals ) +{ + std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl; + BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals."); + try + { + // using statements needed to change precision policy. + using boost::math::policies::policy; + using boost::math::policies::make_policy; + using boost::math::policies::precision; + using boost::math::policies::digits2; + using boost::math::policies::digits10; + + // using statements needed for changing error handling policy. + using boost::math::policies::evaluation_error; + using boost::math::policies::domain_error; + using boost::math::policies::overflow_error; + using boost::math::policies::ignore_error; + using boost::math::policies::throw_on_error; + + typedef policy< + domain_error, + overflow_error + > no_throw_policy; + + /* + // Experiment with better diagnostics. + typedef float Real; + + Real inf = std::numeric_limits::infinity(); + Real max = (std::numeric_limits::max)(); + std::cout.precision(std::numeric_limits::max_digits10); + //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308 + std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf + std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227 + //std::cout << lambert_w0(inf) << std::endl; // inf - will throw. + std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0 + std::cout << "lambert_w0(std::numeric_limits::denorm_min()) = " << lambert_w0(std::numeric_limits::denorm_min()) << std::endl; // 4.94066e-324 + std::cout << "lambert_w0(std::numeric_limits::min()) = " << lambert_w0((std::numeric_limits::min)()) << std::endl; // 2.22507e-308 + + // Approximate the largest lambert_w you can get for type T? + float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max_f " << max_w_f << std::endl; // 84.2879 + Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. + std::cout << "w max " << max_w << std::endl; // 703.227 + + std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; // + std::cout << "test integral 1/z^2" << std::endl; + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1e-10, policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "ULP = " << boost::math::ulp(1., policy >()) << std::endl; // ULP = 2.2204460492503131e-16 + std::cout << "epsilon = " << std::numeric_limits::epsilon() << std::endl; // + std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19 + std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19 + + + +// Demo debug version. +Real tol = std::numeric_limits::epsilon(); +Real x; +{ + using boost::math::quadrature::exp_sinh; + exp_sinh es; + // Function to be integrated, lambert_w0(1/z^2). + + //auto f = [](Real z)->Real + //{ // Naive - no protection against underflow and subsequent divide by zero. + // return lambert_w0(1 / (z * z)); + //}; + // Diagnostic is: + // Error in function boost::math::lambert_w0: Expected a finite value but got inf + + auto f = [](Real z)->Real + { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things. + return debug_integration_proc(z); + }; + // Exception Error in function boost::math::lambert_w0: Expected a finite value but got inf. + + // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163. + // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23. + x = es.integrate(f); + std::cout << "es.integrate(f) = " << x << std::endl; + BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi(), tol); + // root_two_pi(); + } + catch (std::exception& ex) + { + std::cout << ex.what() << std::endl; + } +} +