mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 16:32:10 +00:00
240 lines
9.8 KiB
C++
240 lines
9.8 KiB
C++
// 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^-9 i
|
||
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]
|
||
*/
|