From 75706cb4c28da9ed577317930e064e3f5bb9d27c Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 9 Feb 2021 20:59:43 +0300 Subject: [PATCH 1/2] Remove MPL from lambert_w --- .../math/special_functions/lambert_w.hpp | 286 +++++++++--------- 1 file changed, 139 insertions(+), 147 deletions(-) diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp index bd40bcd88..f161fc1b5 100644 --- a/include/boost/math/special_functions/lambert_w.hpp +++ b/include/boost/math/special_functions/lambert_w.hpp @@ -58,7 +58,6 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include // series functor. //#include // polynomial. #include // evaluate_polynomial. -#include #include #include // boost::math::tools::max_value(). #include @@ -68,13 +67,16 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include #include #include +#include +#include +#include // Needed for testing and diagnostics only. #include #include #include // For float_distance. -typedef double lookup_t; // Type for lookup table (double or float, or even long double?) +using lookup_t = double; // Type for lookup table (double or float, or even long double?) //#include "J:\Cpp\Misc\lambert_w_lookup_table_generator\lambert_w_lookup_table.ipp" // #include "lambert_w_lookup_table.ipp" // Boost.Math version. @@ -104,14 +106,14 @@ namespace lambert_w_detail { //! \param z Argument z for Lambert_w function. //! \returns New estimate of Lambert W, hopefully improved. //! -template +template inline T lambert_w_halley_step(T w_est, const T z) { BOOST_MATH_STD_USING T e = exp(w_est); w_est -= 2 * (w_est + 1) * (e * w_est - z) / (z * (w_est + 2) + e * (w_est * (w_est + 2) + 2)); return w_est; -} // template lambert_w_halley_step(T w_est, T z) +} // template lambert_w_halley_step(T w_est, T z) //! \brief Halley iterate to refine Lambert_w estimate, //! taking at least one Halley_step. @@ -121,9 +123,8 @@ inline T lambert_w_halley_step(T w_est, const T z) //! \tparam T floating-point (or fixed-point) type. //! \param z Argument z for Lambert_w function. //! \param w_est Lambert w estimate. -template -inline - T lambert_w_halley_iterate(T w_est, const T z) +template +inline T lambert_w_halley_iterate(T w_est, const T z) { BOOST_MATH_STD_USING static const T max_diff = boost::math::tools::root_epsilon() * fabs(w_est); @@ -137,21 +138,19 @@ inline diff = fabs(w_est - w_new); } return w_new; -} // template lambert_w_halley_iterate(T w_est, T z) +} // template lambert_w_halley_iterate(T w_est, T z) // Two Halley function versions that either -// single step (if boost::false_type) or iterate (if boost::true_type). +// single step (if std::false_type) or iterate (if std::true_type). // Selected at compile-time using parameter 3. -template -inline -T lambert_w_maybe_halley_iterate(T z, T w, boost::false_type const&) +template +inline T lambert_w_maybe_halley_iterate(T z, T w, std::false_type const&) { return lambert_w_halley_step(z, w); // Single step. } -template -inline -T lambert_w_maybe_halley_iterate(T z, T w, boost::true_type const&) +template +inline T lambert_w_maybe_halley_iterate(T z, T w, std::true_type const&) { return lambert_w_halley_iterate(z, w); // Iterate steps. } @@ -161,30 +160,26 @@ T lambert_w_maybe_halley_iterate(T z, T w, boost::true_type const&) //! reduce argument z to double precision (if true_type). //! Version is selected at compile-time using parameter 2. -template -inline -double maybe_reduce_to_double(const T& z, const boost::true_type&) +template +inline double maybe_reduce_to_double(const T& z, const std::true_type&) { return static_cast(z); // Reduce to double precision. } -template -inline -T maybe_reduce_to_double(const T& z, const boost::false_type&) +template +inline T maybe_reduce_to_double(const T& z, const std::false_type&) { // Don't reduce to double. return z; } -template -inline -double must_reduce_to_double(const T& z, const boost::true_type&) +template +inline double must_reduce_to_double(const T& z, const std::true_type&) { return static_cast(z); // Reduce to double precision. } -template -inline -double must_reduce_to_double(const T& z, const boost::false_type&) +template +inline double must_reduce_to_double(const T& z, const std::false_type&) { // try a lexical_cast and hope for the best: return boost::lexical_cast(z); } @@ -201,8 +196,7 @@ double must_reduce_to_double(const T& z, const boost::false_type&) // Schroeder refinement, called unless NOT required by precision policy. template -inline -T schroeder_update(const T w, const T y) +inline T schroeder_update(const T w, const T y) { // Compute derivatives using 5th order Schroeder refinement. // Since this is the final step, it will always use the highest precision type T. @@ -496,28 +490,28 @@ T lambert_w_singularity_series(const T p) //! but for the tag_type selection to work, they all must include Policy in their signature. // Forward declaration of variants of lambert_w0_small_z. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // for float (32-bit) type. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // for float (32-bit) type. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // for double (64-bit) type. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // for double (64-bit) type. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // for long double (double extended 80-bit) type. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // for long double (double extended 80-bit) type. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // for long double (128-bit) type. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // for long double (128-bit) type. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // for float128 quadmath Q type. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // for float128 quadmath Q type. -template -T lambert_w0_small_z(T x, const Policy&, boost::integral_constant const&); // Generic multiprecision T. +template +T lambert_w0_small_z(T x, const Policy&, std::integral_constant const&); // Generic multiprecision T. // Set tag_type depending on max_digits10. -template +template T lambert_w0_small_z(T x, const Policy& pol) { //std::numeric_limits::max_digits10 == 36 ? 3 : // 128-bit long double. - typedef boost::integral_constant::is_specialized == 0 ? 5 : #ifndef BOOST_NO_CXX11_NUMERIC_LIMITS std::numeric_limits::max_digits10 <= 9 ? 0 : // for float 32-bit. @@ -531,11 +525,10 @@ T lambert_w0_small_z(T x, const Policy& pol) std::numeric_limits::digits <= 64 ? 2 : // for 80-bit double extended. 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; + : 5>; // All Generic multiprecision types. // std::cout << "\ntag type = " << tag_type << std::endl; // error C2275: 'tag_type': illegal use of this type as an expression. return lambert_w0_small_z(x, pol, tag_type()); -} // template T lambert_w0_small_z(T x) +} // template T lambert_w0_small_z(T x) //! Specialization of float (32-bit) series expansion used for small z (abs(z) < 0.05). // Only 9 Coefficients are computed to 21 decimal digits precision, ample for 32-bit float used by most platforms. @@ -543,8 +536,8 @@ T lambert_w0_small_z(T x, const Policy& pol) // N[InverseSeries[Series[z Exp[z],{z,0,34}]],50], // as proposed by Tosio Fukushima and implemented by Darko Veberic. -template -T lambert_w0_small_z(T z, const Policy&, boost::integral_constant const&) +template +T lambert_w0_small_z(T z, const Policy&, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize prec = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -568,15 +561,15 @@ T lambert_w0_small_z(T z, const Policy&, boost::integral_constant const& #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES return result; -} // template T lambert_w0_small_z(T x, boost::integral_constant const&) +} // template T lambert_w0_small_z(T x, std::integral_constant const&) //! Specialization of double (64-bit double) series expansion used for small z (abs(z) < 0.05). // 17 Coefficients are computed to 21 decimal digits precision suitable for 64-bit double used by most platforms. // Taylor series coefficients used are computed by Wolfram to 50 decimal digits using instruction // N[InverseSeries[Series[z Exp[z],{z,0,34}]],50], as proposed by Tosio Fukushima and implemented by Veberic. -template -T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant const&) +template +T lambert_w0_small_z(const T z, const Policy&, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize prec = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -608,15 +601,15 @@ T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES return result; -} // T lambert_w0_small_z(const T z, boost::integral_constant const&) +} // T lambert_w0_small_z(const T z, std::integral_constant const&) //! Specialization of long double (80-bit double extended) series expansion used for small z (abs(z) < 0.05). // 21 Coefficients are computed to 21 decimal digits precision suitable for 80-bit long double used by some // platforms including GCC and Clang when generating for Intel X86 floating-point processors with 80-bit operations enabled (the default). // (This is NOT used by Microsoft Visual Studio where double and long always both use only 64-bit type. // Nor used for 128-bit float128.) -template -T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant const&) +template +T lambert_w0_small_z(const T z, const Policy&, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize precision = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -676,7 +669,7 @@ z * (2.154990206091088289321708745358647e6L // z^20 distance -5 without term 20 std::cout.precision(precision); // Restore. #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES return result; -} // long double lambert_w0_small_z(const T z, boost::integral_constant const&) +} // long double lambert_w0_small_z(const T z, std::integral_constant const&) //! Specialization of 128-bit long double series expansion used for small z (abs(z) < 0.05). // 34 Taylor series coefficients used are computed by Wolfram to 50 decimal digits using instruction @@ -686,8 +679,8 @@ z * (2.154990206091088289321708745358647e6L // z^20 distance -5 without term 20 // nor multiprecision type cpp_bin_float_quad that can only be initialised at full precision of the type // constructed with a decimal digit string like "2.6666666666666666666666666666666666666666666666667".) -template -T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant const&) +template +T lambert_w0_small_z(const T z, const Policy&, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize precision = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -719,7 +712,7 @@ T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant std::cout.precision(precision); // Restore. #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES return result; -} // T lambert_w0_small_z(const T z, boost::integral_constant const&) +} // T lambert_w0_small_z(const T z, std::integral_constant const&) //! Specialization of 128-bit quad series expansion used for small z (abs(z) < 0.05). // 34 Taylor series coefficients used were computed by Wolfram to 50 decimal digits using instruction @@ -732,8 +725,8 @@ T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant // it is slightly slower than getting a double approximation followed by a single Halley step. #ifdef BOOST_HAS_FLOAT128 -template -T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant const&) +template +T lambert_w0_small_z(const T z, const Policy&, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize precision = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -784,24 +777,24 @@ T lambert_w0_small_z(const T z, const Policy&, boost::integral_constant #endif // BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES return result; -} // T lambert_w0_small_z(const T z, boost::integral_constant const&) float128 +} // T lambert_w0_small_z(const T z, std::integral_constant const&) float128 #else -template -inline T lambert_w0_small_z(const T z, const Policy& pol, boost::integral_constant const&) +template +inline T lambert_w0_small_z(const T z, const Policy& pol, std::integral_constant const&) { - return lambert_w0_small_z(z, pol, boost::integral_constant()); + return lambert_w0_small_z(z, pol, std::integral_constant()); } #endif // BOOST_HAS_FLOAT128 //! Series functor to compute series term using pow and factorial. //! \details Functor is called after evaluating polynomial with the coefficients as rationals below. -template +template struct lambert_w0_small_z_series_term { - typedef T result_type; + using result_type = T; //! \param _z Lambert W argument z. //! \param -term -pow<18>(z) / 6402373705728000uLL //! \param _k number of terms == initially 18 @@ -825,11 +818,11 @@ private: int k; T z; T term; -}; // template struct lambert_w0_small_z_series_term +}; // template struct lambert_w0_small_z_series_term //! Generic variant for T a User-defined types like Boost.Multiprecision. -template -inline T lambert_w0_small_z(T z, const Policy& pol, boost::integral_constant const&) +template +inline T lambert_w0_small_z(T z, const Policy& pol, std::integral_constant const&) { #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::streamsize precision = std::cout.precision(std::numeric_limits::max_digits10); // Save. @@ -913,7 +906,7 @@ inline T lambert_w0_small_z(T z, const Policy& pol, boost::integral_constant::max_digits10); T result = evaluate_polynomial(coeff, z); - // template + // template // V evaluate_polynomial(const T(&poly)[N], const V& val); // Size of coeff found from N //std::cout << "evaluate_polynomial(coeff, z); == " << result << std::endl; @@ -942,7 +935,7 @@ inline T lambert_w0_small_z(T z, const Policy& pol, boost::integral_constant(); // Max iterations from policy. + std::uintmax_t max_iter = policies::get_max_series_iterations(); // Max iterations from policy. #ifdef BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES std::cout << "max iter from policy = " << max_iter << std::endl; // // max iter from policy = 1000000 is default. @@ -950,7 +943,7 @@ inline T lambert_w0_small_z(T z, const Policy& pol, boost::integral_constant(), max_iter, result); // result == evaluate_polynomial. - //sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, const U& init_value) + //sum_series(Functor& func, int bits, std::uintmax_t& max_terms, const U& init_value) // std::cout << "sum_series(s, get_epsilon(), max_iter, result); = " << result << std::endl; //T epsilon = get_epsilon(); @@ -965,13 +958,12 @@ inline T lambert_w0_small_z(T z, const Policy& pol, boost::integral_constant inline T lambert_w0_small_z_series(T z, const Policy& pol) +} // template inline T lambert_w0_small_z_series(T z, const Policy& pol) // Approximate lambert_w0 (used for z values that are outside range of lookup table or rational functions) // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162. template -inline -T lambert_w0_approx(T z) +inline T lambert_w0_approx(T z) { BOOST_MATH_STD_USING T lz = log(z); @@ -995,7 +987,7 @@ T lambert_w0_approx(T z) //! For higher precisions, a 64-bit double approximation is computed first, //! and then refined using Halley iterations. -template +template inline T do_get_near_singularity_param(T z) { BOOST_MATH_STD_USING @@ -1003,24 +995,24 @@ inline T do_get_near_singularity_param(T z) const T p = sqrt(p2); return p; } -template +template inline T get_near_singularity_param(T z, const Policy) { - typedef typename policies::evaluation::type value_type; + using value_type = typename policies::evaluation::type; return static_cast(do_get_near_singularity_param(static_cast(z))); } // Forward declarations: -//template T lambert_w0_small_z(T z, const Policy& pol); -//template -//T lambert_w0_imp(T w, const Policy& pol, const boost::integral_constant&); // 32 bit usually float. -//template -//T lambert_w0_imp(T w, const Policy& pol, const boost::integral_constant&); // 64 bit usually double. -//template -//T lambert_w0_imp(T w, const Policy& pol, const boost::integral_constant&); // 80-bit long double. +//template T lambert_w0_small_z(T z, const Policy& pol); +//template +//T lambert_w0_imp(T w, const Policy& pol, const std::integral_constant&); // 32 bit usually float. +//template +//T lambert_w0_imp(T w, const Policy& pol, const std::integral_constant&); // 64 bit usually double. +//template +//T lambert_w0_imp(T w, const Policy& pol, const std::integral_constant&); // 80-bit long double. -template +template T lambert_w_positive_rational_float(T z) { BOOST_MATH_STD_USING @@ -1164,7 +1156,7 @@ T lambert_w_positive_rational_float(T z) } } -template +template T lambert_w_negative_rational_float(T z, const Policy& pol) { BOOST_MATH_STD_USING @@ -1223,19 +1215,19 @@ T lambert_w_negative_rational_float(T z, const Policy& pol) } //! Lambert_w0 @b 'float' implementation, selected when T is 32-bit precision. -template -inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) +template +inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static const char* function = "boost::math::lambert_w0<%1%>"; // For error messages. + static std::string function = "boost::math::lambert_w0<%1%>"; // For error messages. BOOST_MATH_STD_USING // Aid ADL of std functions. if ((boost::math::isnan)(z)) { - return boost::math::policies::raise_domain_error(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); } if ((boost::math::isinf)(z)) { - return boost::math::policies::raise_overflow_error(function, "Expected a finite value but got %1%.", z, pol); + return boost::math::policies::raise_overflow_error(function.c_str(), "Expected a finite value but got %1%.", z, pol); } if (z >= 0.05) // Fukushima switch point. @@ -1246,16 +1238,16 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); return -1; } else // z < 0.05 { return lambert_w_negative_rational_float(z, pol); } -} // T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) for 32-bit usually float. +} // T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) for 32-bit usually float. -template +template T lambert_w_positive_rational_double(T z) { BOOST_MATH_STD_USING @@ -1491,7 +1483,7 @@ T lambert_w_positive_rational_double(T z) } } -template +template T lambert_w_negative_rational_double(T z, const Policy& pol) { BOOST_MATH_STD_USING @@ -1622,10 +1614,10 @@ T lambert_w_negative_rational_double(T z, const Policy& pol) } //! Lambert_w0 @b 'double' implementation, selected when T is 64-bit precision. -template -inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) +template +inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static const char* function = "boost::math::lambert_w0<%1%>"; + static std::string function = "boost::math::lambert_w0<%1%>"; BOOST_MATH_STD_USING // Aid ADL of std functions. // Detect unusual case of 32-bit double with a wider/64-bit long double @@ -1637,11 +1629,11 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); } if ((boost::math::isinf)(z)) { - return boost::math::policies::raise_overflow_error(function, "Expected a finite value but got %1%.", z, pol); + return boost::math::policies::raise_overflow_error(function.c_str(), "Expected a finite value but got %1%.", z, pol); } if (z >= 0.05) @@ -1652,7 +1644,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } return -1; } @@ -1660,22 +1652,22 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) 64-bit precision, usually double. +} // T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) 64-bit precision, usually double. //! lambert_W0 implementation for extended precision types including //! long double (80-bit and 128-bit), ??? //! quad float128, Boost.Multiprecision types like cpp_bin_float_quad, cpp_bin_float_50... -template -inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) +template +inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static const char* function = "boost::math::lambert_w0<%1%>"; + static std::string function = "boost::math::lambert_w0<%1%>"; BOOST_MATH_STD_USING // Aid ADL of std functions. // Filter out special cases first: if ((boost::math::isnan)(z)) { - return boost::math::policies::raise_domain_error(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } if (fabs(z) <= 0.05f) { @@ -1686,7 +1678,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant(function, 0, pol); + return policies::raise_overflow_error(function.c_str(), 0, pol); // Or might return infinity if available else max_value, // but other Boost.Math special functions raise overflow. } @@ -1707,7 +1699,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } // z is very close (within 0.01) of the branch singularity at -e^-1 // so use a series approximation proposed by Corless et al. @@ -1723,27 +1715,27 @@ inline T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant::type precision_type; - typedef boost::integral_constant::type; + using tag_type = std::integral_constant 113) ? true // Unknown at compile-time, variable/arbitrary, or more than float128 or cpp_bin_quad 128-bit precision. : false // float, double, float128, cpp_bin_quad 128-bit, so single Halley step. - > tag_type; + >; // For speed, we also cast z to type double when that is possible - // if (boost::is_constructible() == true). - T w = lambert_w0_imp(maybe_reduce_to_double(z, boost::is_constructible()), pol, boost::integral_constant()); + // if (std::is_constructible() == true). + T w = lambert_w0_imp(maybe_reduce_to_double(z, std::is_constructible()), pol, std::integral_constant()); return lambert_w_maybe_halley_iterate(w, z, tag_type()); -} // T lambert_w0_imp(T z, const Policy& pol, const boost::integral_constant&) all extended precision types. +} // T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) all extended precision types. // Lambert w-1 implementation // ============================================================================================== //! Lambert W for W-1 branch, -max(z) < z <= -1/e. // TODO is -max(z) allowed? -template +template T lambert_wm1_imp(const T z, const Policy& pol) { // Catch providing an integer value as parameter x to lambert_w, for example, lambert_w(1). @@ -1758,19 +1750,19 @@ T lambert_wm1_imp(const T z, const Policy& pol) BOOST_MATH_STD_USING // Aid argument dependent lookup (ADL) of abs. - const char* function = "boost::math::lambert_wm1()"; // Used for error messages. + std::string function = "boost::math::lambert_wm1()"; // Used for error messages. // Check for edge and corner cases first: if ((boost::math::isnan)(z)) { - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z is NaN!", z, pol); } // isnan if ((boost::math::isinf)(z)) { - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z is infinite!", z, pol); } // isinf @@ -1790,7 +1782,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) { // All real types except arbitrary precision. if (!(boost::math::isnormal)(z)) { // Almost zero - might also just return infinity like z == 0 or max_value? - return policies::raise_overflow_error(function, + return policies::raise_overflow_error(function.c_str(), "Argument z = %1% is denormalized! (must be z > (std::numeric_limits::min)() or z == 0)", z, pol); } @@ -1798,7 +1790,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) if (z > static_cast(0)) { // - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z = %1% is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)", z, pol); } @@ -1806,7 +1798,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) { // z is denormalized, so cannot be computed. // -std::numeric_limits::min() is smallest for type T, // for example, for double: lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634 - return policies::raise_overflow_error(function, + return policies::raise_overflow_error(function.c_str(), "Argument z = %1% is too small (z < -std::numeric_limits::min so denormalized) for Lambert W-1 branch!", z, pol); } @@ -1817,7 +1809,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) // z is too negative for the W-1 (or W0) branch. if (z < -boost::math::constants::exp_minus_one()) // > singularity/branch point z = -exp(-1) = -3.6787944. { - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z = %1% is out of range (z < -exp(-1) = -3.6787944... <= 0) for Lambert W-1 (or W0) branch!", z, pol); } @@ -1843,7 +1835,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) return w_series; } // Should not get here. - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z = %1% is out of range for Lambert W-1 branch. (Should not get here - please report!)", z, pol); } // if (z < -0.35) @@ -1919,7 +1911,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(must_reduce_to_double(z, boost::is_constructible()), policy >()))); + T double_approx(static_cast(lambert_wm1_imp(must_reduce_to_double(z, std::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; @@ -1955,7 +1947,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) } // else z < g[63] == -1.0264389699511303e-26, so Lambert W-1 integer part > 64. // This should not now occur (should be caught by test and code above) so should be a logic_error? - return policies::raise_domain_error(function, + return policies::raise_domain_error(function.c_str(), "Argument z = %1% is too small (< -1.026439e-26) (logic error - please report!)", z, pol); overshot: @@ -2009,7 +2001,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) using lambert_w_lookup::halves; using lambert_w_lookup::sqrtwm1s; - typedef typename mpl::if_c::value, lookup_t, T>::type calc_type; + using calc_type = typename std::conditional::value, lookup_t, T>::type; calc_type w = -static_cast(n); // Equation 25, calc_type y = static_cast(z * wm1es[n - 1]); // Equation 26, @@ -2045,82 +2037,82 @@ T lambert_wm1_imp(const T z, const Policy& pol) ///////////////////////////// User Lambert w functions. ////////////////////////////// //! Lambert W0 using User-defined policy. - template + template inline typename boost::math::tools::promote_args::type lambert_w0(T z, const Policy& pol) { // Promote integer or expression template arguments to double, // without doing any other internal promotion like float to double. - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; // Work out what precision has been selected, // based on the Policy and the number type. - typedef typename policies::precision::type precision_type; + using precision_type = typename policies::precision::type; // and then select the correct implementation based on that precision (not the type T): - typedef boost::integral_constant 53) ? 0 // either variable precision (0), or greater than 64-bit precision. : (precision_type::value <= 24) ? 1 // 32-bit (probably float) precision. : 2 // 64-bit (probably double) precision. - > tag_type; + >; return lambert_w_detail::lambert_w0_imp(result_type(z), pol, tag_type()); // } // lambert_w0(T z, const Policy& pol) //! Lambert W0 using default policy. - template + template inline typename tools::promote_args::type lambert_w0(T z) { // Promote integer or expression template arguments to double, // without doing any other internal promotion like float to double. - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; // 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; + using precision_type = typename policies::precision>::type; // and then select the correct implementation based on that (not the type T): - typedef boost::integral_constant 53) ? 0 // either variable precision (0), or greater than 64-bit precision. : (precision_type::value <= 24) ? 1 // 32-bit (probably float) precision. : 2 // 64-bit (probably double) precision. - > tag_type; + >; return lambert_w_detail::lambert_w0_imp(result_type(z), policies::policy<>(), tag_type()); } // lambert_w0(T z) using default policy. //! W-1 branch (-max(z) < z <= -1/e). //! Lambert W-1 using User-defined policy. - template + template inline typename tools::promote_args::type lambert_wm1(T z, const Policy& pol) { // Promote integer or expression template arguments to double, // without doing any other internal promotion like float to double. - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; return lambert_w_detail::lambert_wm1_imp(result_type(z), pol); // } //! Lambert W-1 using default policy. - template + template inline typename tools::promote_args::type lambert_wm1(T z) { - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; return lambert_w_detail::lambert_wm1_imp(result_type(z), policies::policy<>()); } // lambert_wm1(T z) // First derivative of Lambert W0 and W-1. - template + template inline typename tools::promote_args::type lambert_w0_prime(T z, const Policy& pol) { - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; using std::numeric_limits; if (z == 0) { @@ -2144,19 +2136,19 @@ T lambert_wm1_imp(const T z, const Policy& pol) return w / (z * (1 + w)); } // lambert_w0_prime(T z) - template + template inline typename tools::promote_args::type lambert_w0_prime(T z) { return lambert_w0_prime(z, policies::policy<>()); } - template + template inline typename tools::promote_args::type lambert_wm1_prime(T z, const Policy& pol) { using std::numeric_limits; - typedef typename tools::promote_args::type result_type; + using result_type = typename tools::promote_args::type; //if (z == 0) //{ // return static_cast(1); @@ -2171,7 +2163,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) return w/(z*(1+w)); } // lambert_wm1_prime(T z) - template + template inline typename tools::promote_args::type lambert_wm1_prime(T z) { From b5031bc9c436438fdfe3496b384490c6414dcbb6 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 10 Feb 2021 19:35:04 +0300 Subject: [PATCH 2/2] revert conversion to string [CI SKIP] --- .../math/special_functions/lambert_w.hpp | 43 +++++++++---------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/include/boost/math/special_functions/lambert_w.hpp b/include/boost/math/special_functions/lambert_w.hpp index f161fc1b5..09ebbcc90 100644 --- a/include/boost/math/special_functions/lambert_w.hpp +++ b/include/boost/math/special_functions/lambert_w.hpp @@ -68,7 +68,6 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include #include #include -#include #include // Needed for testing and diagnostics only. @@ -1218,16 +1217,16 @@ T lambert_w_negative_rational_float(T z, const Policy& pol) template inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static std::string function = "boost::math::lambert_w0<%1%>"; // For error messages. + static const char* function = "boost::math::lambert_w0<%1%>"; // For error messages. BOOST_MATH_STD_USING // Aid ADL of std functions. if ((boost::math::isnan)(z)) { - return boost::math::policies::raise_domain_error(function.c_str(), "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); } if ((boost::math::isinf)(z)) { - return boost::math::policies::raise_overflow_error(function.c_str(), "Expected a finite value but got %1%.", z, pol); + return boost::math::policies::raise_overflow_error(function, "Expected a finite value but got %1%.", z, pol); } if (z >= 0.05) // Fukushima switch point. @@ -1238,7 +1237,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); return -1; } else // z < 0.05 @@ -1617,7 +1616,7 @@ T lambert_w_negative_rational_double(T z, const Policy& pol) template inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static std::string function = "boost::math::lambert_w0<%1%>"; + static const char* function = "boost::math::lambert_w0<%1%>"; BOOST_MATH_STD_USING // Aid ADL of std functions. // Detect unusual case of 32-bit double with a wider/64-bit long double @@ -1629,11 +1628,11 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant(function.c_str(), "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected a value > -e^-1 (-0.367879...) but got %1%.", z, pol); } if ((boost::math::isinf)(z)) { - return boost::math::policies::raise_overflow_error(function.c_str(), "Expected a finite value but got %1%.", z, pol); + return boost::math::policies::raise_overflow_error(function, "Expected a finite value but got %1%.", z, pol); } if (z >= 0.05) @@ -1644,7 +1643,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } return -1; } @@ -1661,13 +1660,13 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant&) { - static std::string function = "boost::math::lambert_w0<%1%>"; + static const char* function = "boost::math::lambert_w0<%1%>"; BOOST_MATH_STD_USING // Aid ADL of std functions. // Filter out special cases first: if ((boost::math::isnan)(z)) { - return boost::math::policies::raise_domain_error(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } if (fabs(z) <= 0.05f) { @@ -1678,7 +1677,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant(function.c_str(), 0, pol); + return policies::raise_overflow_error(function, 0, pol); // Or might return infinity if available else max_value, // but other Boost.Math special functions raise overflow. } @@ -1699,7 +1698,7 @@ inline T lambert_w0_imp(T z, const Policy& pol, const std::integral_constant(function.c_str(), "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); + return boost::math::policies::raise_domain_error(function, "Expected z >= -e^-1 (-0.367879...) but got %1%.", z, pol); } // z is very close (within 0.01) of the branch singularity at -e^-1 // so use a series approximation proposed by Corless et al. @@ -1750,19 +1749,19 @@ T lambert_wm1_imp(const T z, const Policy& pol) BOOST_MATH_STD_USING // Aid argument dependent lookup (ADL) of abs. - std::string function = "boost::math::lambert_wm1()"; // Used for error messages. + const char* function = "boost::math::lambert_wm1()"; // Used for error messages. // Check for edge and corner cases first: if ((boost::math::isnan)(z)) { - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z is NaN!", z, pol); } // isnan if ((boost::math::isinf)(z)) { - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z is infinite!", z, pol); } // isinf @@ -1782,7 +1781,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) { // All real types except arbitrary precision. if (!(boost::math::isnormal)(z)) { // Almost zero - might also just return infinity like z == 0 or max_value? - return policies::raise_overflow_error(function.c_str(), + return policies::raise_overflow_error(function, "Argument z = %1% is denormalized! (must be z > (std::numeric_limits::min)() or z == 0)", z, pol); } @@ -1790,7 +1789,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) if (z > static_cast(0)) { // - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z = %1% is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)", z, pol); } @@ -1798,7 +1797,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) { // z is denormalized, so cannot be computed. // -std::numeric_limits::min() is smallest for type T, // for example, for double: lambert_wm1(-2.2250738585072014e-308) = -714.96865723796634 - return policies::raise_overflow_error(function.c_str(), + return policies::raise_overflow_error(function, "Argument z = %1% is too small (z < -std::numeric_limits::min so denormalized) for Lambert W-1 branch!", z, pol); } @@ -1809,7 +1808,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) // z is too negative for the W-1 (or W0) branch. if (z < -boost::math::constants::exp_minus_one()) // > singularity/branch point z = -exp(-1) = -3.6787944. { - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z = %1% is out of range (z < -exp(-1) = -3.6787944... <= 0) for Lambert W-1 (or W0) branch!", z, pol); } @@ -1835,7 +1834,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) return w_series; } // Should not get here. - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z = %1% is out of range for Lambert W-1 branch. (Should not get here - please report!)", z, pol); } // if (z < -0.35) @@ -1947,7 +1946,7 @@ T lambert_wm1_imp(const T z, const Policy& pol) } // else z < g[63] == -1.0264389699511303e-26, so Lambert W-1 integer part > 64. // This should not now occur (should be caught by test and code above) so should be a logic_error? - return policies::raise_domain_error(function.c_str(), + return policies::raise_domain_error(function, "Argument z = %1% is too small (< -1.026439e-26) (logic error - please report!)", z, pol); overshot: