123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251 |
- // Copyright Paul A. Bristow 2016, 2017.
- // Distributed under 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).
- // Build and run a simple examples of Lambert W function.
- // Some macros that will show some(or much) diagnostic values if #defined.
- //#define-able macros
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W0 branch diagnostics.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_Wm1 // W1 branch diagnostics.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY // Halley refinement diagnostics.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER // Schroeder refinement diagnostics.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS // Number of terms used for near-singularity series.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_NOT_BUILTIN // higher than built-in precision types approximation and refinement.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES // Show evaluation of series near branch singularity.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of series for small z.
- //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_LOOKUP // Show results from lookup table.
- #include <boost/config.hpp> // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
- #include <boost/version.hpp> // for BOOST_MSVC versions.
- #include <boost/cstdint.hpp>
- #include <boost/exception/exception.hpp> // boost::exception
- #include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
- #include <boost/math/policies/policy.hpp>
- #include <boost/multiprecision/cpp_dec_float.hpp> // boost::multiprecision::cpp_dec_float_50
- using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type.
- using boost::multiprecision::cpp_dec_float_100; // 100 decimal digits type.
- using boost::multiprecision::backends::cpp_dec_float;
- using boost::multiprecision::number;
- typedef number<cpp_dec_float<1000> > cpp_dec_float_1000; // 1000 decimal digit types
- #include <boost/multiprecision/cpp_bin_float.hpp>
- using boost::multiprecision::cpp_bin_float_double; // == double
- using boost::multiprecision::cpp_bin_float_double_extended; // 80-bit long double emulation.
- using boost::multiprecision::cpp_bin_float_quad; // 128-bit quad precision.
- //[lambert_w_simple_examples_includes
- #include <boost/math/special_functions/lambert_w.hpp> // For lambert_w function.
- using boost::math::lambert_w0;
- using boost::math::lambert_wm1;
- //] //[/lambert_w_simple_examples_includes]
- #include <iostream>
- // using std::cout;
- // using std::endl;
- #include <exception>
- #include <stdexcept>
- #include <string>
- #include <limits> // For std::numeric_limits.
- //! Show value of z to the full possibly-significant max_digits10 precision of type T.
- template<typename T>
- void show_value(T z)
- {
- std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10); // Save.
- std::cout.precision(std::numeric_limits<T>::max_digits10); // Show all posssibly significant digits.
- std::ios::fmtflags flags(std::cout.flags());
- std::cout.setf(std::ios_base::showpoint); // Include any trailing zeros.
- std::cout << z;
- // Restore:
- std::cout.precision(precision);
- std::cout.flags(flags);
- } // template<typename T> void show_value(T z)
- int main()
- {
- try
- {
- std::cout << "Lambert W simple examples." << std::endl;
- using boost::math::constants::exp_minus_one; //-1/e, the branch point, a singularity ~= -0.367879.
- // using statements needed for changing error handling policy.
- using boost::math::policies::policy;
- using boost::math::policies::make_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;
- {
- //[lambert_w_simple_examples_0
- std::cout.precision(std::numeric_limits<double>::max_digits10);
- // Show all potentially significant decimal digits,
- std::cout << std::showpoint << std::endl;
- // and show significant trailing zeros too.
- double z = 10.;
- double r = lambert_w0(z); // Default policy for double.
- std::cout << "lambert_w0(z) = " << r << std::endl;
- // lambert_w0(z) = 1.7455280027406994
- //] [/lambert_w_simple_examples_0]
- }
- {
- // Other floating-point types can be used too, here float.
- // It is convenient to use a function like `show_value`
- // to display all potentially significant decimal digits
- // for the type, including any significant trailing zeros.
- //[lambert_w_simple_examples_1
- float z = 10.F;
- float r;
- r = lambert_w0(z); // Default policy digits10 = 7, digits2 = 24
- std::cout << "lambert_w0(";
- show_value(z);
- std::cout << ") = ";
- show_value(r);
- std::cout << std::endl; // lambert_w0(10.0000000) = 1.74552798
- //] //[/lambert_w_simple_examples_1]
- }
- {
- // Example of an integer argument to lambert_w,
- // showing that an integer is correctly promoted to a double.
- //[lambert_w_simple_examples_2
- std::cout.precision(std::numeric_limits<double>::max_digits10);
- double r = lambert_w0(10); // Pass an int argument "10" that should be promoted to double argument.
- std::cout << "lambert_w0(10) = " << r << std::endl; // lambert_w0(10) = 1.7455280027406994
- double rp = lambert_w0(10);
- std::cout << "lambert_w0(10) = " << rp << std::endl;
- // lambert_w0(10) = 1.7455280027406994
- auto rr = lambert_w0(10); // C++11 needed.
- std::cout << "lambert_w0(10) = " << rr << std::endl;
- // lambert_w0(10) = 1.7455280027406994 too, showing that rr has been promoted to double.
- //] //[/lambert_w_simple_examples_2]
- }
- {
- // Using multiprecision types to get much higher precision is painless.
- //[lambert_w_simple_examples_3
- cpp_dec_float_50 z("10");
- // Note construction using a decimal digit string "10",
- // NOT a floating-point double literal 10.
- cpp_dec_float_50 r;
- r = lambert_w0(z);
- std::cout << "lambert_w0("; show_value(z); std::cout << ") = ";
- show_value(r);
- std::cout << std::endl;
- // lambert_w0(10.000000000000000000000000000000000000000000000000000000000000000000000000000000) =
- // 1.7455280027406993830743012648753899115352881290809413313533156980404446940000000
- //] //[/lambert_w_simple_examples_3]
- }
- // Using multiprecision types to get multiprecision precision wrong!
- {
- //[lambert_w_simple_examples_4
- cpp_dec_float_50 z(0.7777777777777777777777777777777777777777777777777777777777777777777777777);
- // Compiler evaluates the nearest double-precision binary representation,
- // from the max_digits10 of the floating_point literal double 0.7777777777777777777777777777...,
- // so any extra digits in the multiprecision type
- // beyond max_digits10 (usually 17) are random and meaningless.
- cpp_dec_float_50 r;
- r = lambert_w0(z);
- std::cout << "lambert_w0(";
- show_value(z);
- std::cout << ") = "; show_value(r);
- std::cout << std::endl;
- // lambert_w0(0.77777777777777779011358916250173933804035186767578125000000000000000000000000000)
- // = 0.48086152073210493501934682309060873341910109230469724725005039758139532631901386
- //] //[/lambert_w_simple_examples_4]
- }
- {
- //[lambert_w_simple_examples_4a
- cpp_dec_float_50 z(0.9); // Construct from floating_point literal double 0.9.
- cpp_dec_float_50 r;
- r = lambert_w0(0.9);
- std::cout << "lambert_w0(";
- show_value(z);
- std::cout << ") = "; show_value(r);
- std::cout << std::endl;
- // lambert_w0(0.90000000000000002220446049250313080847263336181640625000000000000000000000000000)
- // = 0.52983296563343440510607251781038939952850341796875000000000000000000000000000000
- std::cout << "lambert_w0(0.9) = " << lambert_w0(static_cast<double>(0.9))
- // lambert_w0(0.9)
- // = 0.52983296563343441
- << std::endl;
- //] //[/lambert_w_simple_examples_4a]
- }
- {
- // Using multiprecision types to get multiprecision precision right!
- //[lambert_w_simple_examples_4b
- cpp_dec_float_50 z("0.9"); // Construct from decimal digit string.
- cpp_dec_float_50 r;
- r = lambert_w0(z);
- std::cout << "lambert_w0(";
- show_value(z);
- std::cout << ") = "; show_value(r);
- std::cout << std::endl;
- // 0.90000000000000000000000000000000000000000000000000000000000000000000000000000000)
- // = 0.52983296563343441213336643954546304857788132269804249284012528304239956413801252
- //] //[/lambert_w_simple_examples_4b]
- }
- // Getting extreme precision (1000 decimal digits) Lambert W values.
- {
- std::cout.precision(std::numeric_limits<cpp_dec_float_1000>::digits10);
- cpp_dec_float_1000 z("2.0");
- cpp_dec_float_1000 r;
- r = lambert_w0(z);
- std::cout << "lambert_w0(z) = " << r << std::endl;
- // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
- // Wolfram alpha command N[productlog[0, 2.0],1000] gives the identical result:
- // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
- }
- {
- //[lambert_w_simple_examples_error_policies
- // Define an error handling policy:
- typedef policy<
- domain_error<throw_on_error>,
- overflow_error<ignore_error> // possibly unwise?
- > my_throw_policy;
- std::cout.precision(std::numeric_limits<double>::max_digits10);
- // Show all potentially significant decimal digits,
- std::cout << std::showpoint << std::endl;
- // and show significant trailing zeros too.
- double z = +1;
- std::cout << "Lambert W (" << z << ") = " << lambert_w0(z) << std::endl;
- // Lambert W (1.0000000000000000) = 0.56714329040978384
- std::cout << "\nLambert W (" << z << ", my_throw_policy()) = "
- << lambert_w0(z, my_throw_policy()) << std::endl;
- // Lambert W (1.0000000000000000, my_throw_policy()) = 0.56714329040978384
- //] //[/lambert_w_simple_example_error_policies]
- }
- {
- // Show error reporting from passing a value to lambert_wm1 that is out of range,
- // (and probably was meant to be passed to lambert_0 instead).
- //[lambert_w_simple_examples_out_of_range
- double z = +1.;
- double r = lambert_wm1(z);
- std::cout << "lambert_wm1(+1.) = " << r << std::endl;
- //] [/lambert_w_simple_examples_out_of_range]
- // Error in function boost::math::lambert_wm1<RealType>(<RealType>):
- // Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
- }
- }
- catch (std::exception& ex)
- {
- std::cout << ex.what() << std::endl;
- }
- } // int main()
- /*
- Output:
- //[lambert_w_simple_examples_error_message_1
- Error in function boost::math::lambert_wm1<RealType>(<RealType>):
- Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
- //] [/lambert_w_simple_examples_error_message_1]
- */
|