123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239 |
- // Copyright Paul A. Bristow 2016.
- // 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).
- // Test that can build and run a simple example of Lambert W function,
- // using algorithm of Thomas Luu.
- // https://svn.boost.org/trac/boost/ticket/11027
- #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.
- #define BOOST_MATH_INSTRUMENT_LAMBERT_W // #define only for diagnostic output.
- // For lambert_w function.
- #include <boost/math/special_functions/lambert_w.hpp>
- #include <iostream>
- // using std::cout;
- // using std::endl;
- #include <exception>
- #include <stdexcept>
- #include <string>
- #include <limits> // For std::numeric_limits.
- //! Show information about build, architecture, address model, platform, ...
- std::string show_versions()
- {
- std::ostringstream message;
- message << "Program: " << __FILE__ << "\n";
- #ifdef __TIMESTAMP__
- message << __TIMESTAMP__;
- #endif
- message << "\nBuildInfo:\n" " Platform " << BOOST_PLATFORM;
- // http://stackoverflow.com/questions/1505582/determining-32-vs-64-bit-in-c
- #if defined(__LP64__) || defined(_WIN64) || (defined(__x86_64__) && !defined(__ILP32__) ) || defined(_M_X64) || defined(__ia64) || defined (_M_IA64) || defined(__aarch64__) || defined(__powerpc64__)
- #define IS64BIT 1
- message << ", 64-bit.";
- #else
- #define IS32BIT 1
- message << ", 32-bit.";
- #endif
- message << "\n Compiler " BOOST_COMPILER;
- #ifdef BOOST_MSC_VER
- #ifdef _MSC_FULL_VER
- message << "\n MSVC version " << BOOST_STRINGIZE(_MSC_FULL_VER) << ".";
- #endif
- #ifdef __WIN64
- mess age << "\n WIN64" << std::endl;
- #endif // __WIN64
- #ifdef _WIN32
- message << "\n WIN32" << std::endl;
- #endif // __WIN32
- #endif
- #ifdef __GNUC__
- //PRINT_MACRO(__GNUC__);
- //PRINT_MACRO(__GNUC_MINOR__);
- //PRINT_MACRO(__GNUC_PATCH__);
- std::cout << "GCC " << __VERSION__ << std::endl;
- //PRINT_MACRO(LONG_MAX);
- #endif // __GNUC__
- message << "\n STL " << BOOST_STDLIB;
- message << "\n Boost version " << BOOST_VERSION / 100000 << "." << BOOST_VERSION / 100 % 1000 << "." << BOOST_VERSION % 100;
- #ifdef BOOST_HAS_FLOAT128
- message << ", BOOST_HAS_FLOAT128" << std::endl;
- #endif
- message << std::endl;
- return message.str();
- } // std::string versions()
- int main()
- {
- try
- {
- //std::cout << "Lambert W example basic!" << std::endl;
- //std::cout << show_versions() << std::endl;
- //std::cout << exp(1) << std::endl; // 2.71828
- //std::cout << exp(-1) << std::endl; // 0.367879
- //std::cout << std::numeric_limits<double>::epsilon() / 2 << std::endl; // 1.11022e-16
- using namespace boost::math;
- using boost::math::constants::exp_minus_one;
- double x = 1.;
- double W1 = lambert_w(1.);
- // Note, NOT integer X, for example: lambert_w(1); or will get message like
- // error C2338: Must be floating-point, not integer type, for example W(1.), not W(1)!
- //
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.567143
- // This 'golden ratio' for exponentials is http://mathworld.wolfram.com/OmegaConstant.html
- // since exp[-W(1)] = W(1)
- // A030178 Decimal expansion of LambertW(1): the solution to x*exp(x)
- // = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866
- // http://oeis.org/A030178
- double expplogone = exp(-lambert_w(1.));
- if (expplogone != W1)
- {
- std::cout << expplogone << " " << W1 << std::endl; //
- }
- //[lambert_w_example_1
- x = 0.01;
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.00990147
- //] [/lambert_w_example_1]
- x = -0.01;
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // -0.0101015
- x = -0.1;
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
- /**/
- for (double xd = 1.; xd < 1e20; xd *= 10)
- {
- // 1. 0.56714329040978387
- // 0.56714329040978384
- // 10 1.7455280027406994
- // 1.7455280027406994
- // 100 3.3856301402900502
- // 3.3856301402900502
- // 1000 5.2496028524015959
- // 5.249602852401596227126056319697306282521472386059592844451465483991362228320942832739693150854347718
- // 1e19 40.058769161984308
- // 40.05876916198431163898797971203180915622644925765346546858291325452428038208071849105889199253335063
- std::cout << "Lambert W (" << xd << ") = " << lambert_w(xd) << std::endl; //
- }
- //
- // Test near singularity.
- // http://www.wolframalpha.com/input/?i=N%5Blambert_w%5B-0.367879%5D,17%5D test value N[lambert_w[-0.367879],17]
- // -0.367879441171442321595523770161460867445811131031767834
- x = -0.367879; // < -exp(1) = -0.367879
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // Lambert W (-0.36787900000000001) = -0.99845210378080340
- // -0.99845210378080340
- // -0.99845210378072726 N[lambert_w[-0.367879],17] wolfram so very close.
- x = -0.3678794; // expect -0.99952696660756813
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- x = -0.36787944; // expect -0.99992019848408340
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- x = -0.367879441; // -0.99996947070054883
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- x = -0.36787944117; // -0.99999719977527159
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- x = -0.367879441171; // -0.99999844928821992
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- x = -exp_minus_one<double>() + std::numeric_limits<double>::epsilon();
- // Lambert W (-0.36787944117144211) = -0.99999996349975895
- // N[lambert_w[-0.36787944117144211],17] == -0.99999996608315303
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- std::cout << " 1 - sqrt(eps) = " << static_cast<double>(1) - sqrt(std::numeric_limits<double>::epsilon()) << std::endl;
- x = -exp_minus_one<double>();
- // N[lambert_w[-0.36787944117144233],17] == -1.000000000000000 + 6.7595465843924897*10^-9i
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
- // At Singularity - 0.36787944117144233 == -0.36787944117144233 returned - 1.0000000000000000
- // Lambert W(-0.36787944117144233) = -1.0000000000000000
- x = (std::numeric_limits<double>::max)()/4;
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // OK 702.023799146706
- x = (std::numeric_limits<double>::max)()/2;
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
- x = (std::numeric_limits<double>::max)();
- std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
- // Error in function boost::math::log1p<double>(double): numeric overflow
- /* */
- }
- catch (std::exception& ex)
- {
- std::cout << ex.what() << std::endl;
- }
- } // int main()
- /*
- //[lambert_w_output_1
- Output:
- 1> example_basic.cpp
- 1> Generating code
- 1> All 237 functions were compiled because no usable IPDB/IOBJ from previous compilation was found.
- 1> Finished generating code
- 1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.exe
- 1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.pdb (Full PDB)
- 1> Lambert W example basic!
- 1> Platform: Win32
- 1> Compiler: Microsoft Visual C++ version 14.0
- 1> STL : Dinkumware standard library version 650
- 1> Boost : 1.63.0
- 1> _MSC_FULL_VER = 190024123
- 1> Win32
- 1> x64
- 1> (x64)
- 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
- 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
- 1> Final 0.567143290409784 after 2 iterations, difference = 0
- 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
- 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
- 1> Final 0.567143290409784 after 2 iterations, difference = 0
- 1> Lambert W (1) = 0.567143290409784
- 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
- 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
- 1> Final 0.567143290409784 after 2 iterations, difference = 0
- 1> Iteration #0, w0 0.0099072820916067, w1 = 0.00990147384359511, difference = 5.92416060777624e-06, relative 0.000586604388734591
- 1> Final 0.00990147384359511 after 1 iterations, difference = 0
- 1> Lambert W (0.01) = 0.00990147384359511
- 1> Iteration #0, w0 -0.0101016472705154, w1 = -0.0101015271985388, difference = -1.17664437923951e-07, relative 1.18865171889748e-05
- 1> Final -0.0101015271985388 after 1 iterations, difference = 0
- 1> Lambert W (-0.01) = -0.0101015271985388
- 1> Iteration #0, w0 -0.111843322610692, w1 = -0.111832559158964, difference = -8.54817065376601e-06, relative 9.62461362694622e-05
- 1> Iteration #1, w0 -0.111832559158964, w1 = -0.111832559158963, difference = -5.68989300120393e-16, relative 6.43929354282591e-15
- 1> Final -0.111832559158963 after 2 iterations, difference = 0
- 1> Lambert W (-0.1) = -0.111832559158963
- 1> Iteration #0, w0 -0.998452103785573, w1 = -0.998452103780803, difference = -2.72004641033163e-15, relative 4.77662354114727e-12
- 1> Final -0.998452103780803 after 1 iterations, difference = 0
- 1> Lambert W (-0.367879) = -0.998452103780803
- //] [/lambert_w_output_1]
- */
|