252 lines
13 KiB
C++
252 lines
13 KiB
C++
// 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]
|
|
|
|
*/
|