1 // Copyright Paul A. Bristow 2016, 2017.
2
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or
5 // copy at http ://www.boost.org/LICENSE_1_0.txt).
6
7 // Build and run a simple examples of Lambert W function.
8
9 // Some macros that will show some(or much) diagnostic values if #defined.
10 //#define-able macros
11 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0 // W0 branch diagnostics.
12 //#define BOOST_MATH_INSTRUMENT_LAMBERT_Wm1 // W1 branch diagnostics.
13 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_HALLEY // Halley refinement diagnostics.
14 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SCHROEDER // Schroeder refinement diagnostics.
15 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_TERMS // Number of terms used for near-singularity series.
16 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_NOT_BUILTIN // higher than built-in precision types approximation and refinement.
17 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SINGULARITY_SERIES // Show evaluation of series near branch singularity.
18 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of series for small z.
19 //#define BOOST_MATH_INSTRUMENT_LAMBERT_W0_LOOKUP // Show results from lookup table.
20
21 #include <boost/config.hpp> // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
22 #include <boost/version.hpp> // for BOOST_MSVC versions.
23 #include <boost/cstdint.hpp>
24 #include <boost/exception/exception.hpp> // boost::exception
25 #include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
26 #include <boost/math/policies/policy.hpp>
27
28 #include <boost/multiprecision/cpp_dec_float.hpp> // boost::multiprecision::cpp_dec_float_50
29 using boost::multiprecision::cpp_dec_float_50; // 50 decimal digits type.
30 using boost::multiprecision::cpp_dec_float_100; // 100 decimal digits type.
31 using boost::multiprecision::backends::cpp_dec_float;
32 using boost::multiprecision::number;
33 typedef number<cpp_dec_float<1000> > cpp_dec_float_1000; // 1000 decimal digit types
34
35 #include <boost/multiprecision/cpp_bin_float.hpp>
36 using boost::multiprecision::cpp_bin_float_double; // == double
37 using boost::multiprecision::cpp_bin_float_double_extended; // 80-bit long double emulation.
38 using boost::multiprecision::cpp_bin_float_quad; // 128-bit quad precision.
39
40 //[lambert_w_simple_examples_includes
41 #include <boost/math/special_functions/lambert_w.hpp> // For lambert_w function.
42
43 using boost::math::lambert_w0;
44 using boost::math::lambert_wm1;
45 //] //[/lambert_w_simple_examples_includes]
46
47 #include <iostream>
48 // using std::cout;
49 // using std::endl;
50 #include <exception>
51 #include <stdexcept>
52 #include <string>
53 #include <limits> // For std::numeric_limits.
54
55 //! Show value of z to the full possibly-significant max_digits10 precision of type T.
56 template<typename T>
show_value(T z)57 void show_value(T z)
58 {
59 std::streamsize precision = std::cout.precision(std::numeric_limits<T>::max_digits10); // Save.
60 std::cout.precision(std::numeric_limits<T>::max_digits10); // Show all possibly significant digits.
61 std::ios::fmtflags flags(std::cout.flags());
62 std::cout.setf(std::ios_base::showpoint); // Include any trailing zeros.
63 std::cout << z;
64 // Restore:
65 std::cout.precision(precision);
66 std::cout.flags(flags);
67 } // template<typename T> void show_value(T z)
68
main()69 int main()
70 {
71 try
72 {
73 std::cout << "Lambert W simple examples." << std::endl;
74
75 using boost::math::constants::exp_minus_one; //-1/e, the branch point, a singularity ~= -0.367879.
76
77 // using statements needed for changing error handling policy.
78 using boost::math::policies::policy;
79 using boost::math::policies::make_policy;
80 using boost::math::policies::evaluation_error;
81 using boost::math::policies::domain_error;
82 using boost::math::policies::overflow_error;
83 using boost::math::policies::ignore_error;
84 using boost::math::policies::throw_on_error;
85
86 {
87 //[lambert_w_simple_examples_0
88 std::cout.precision(std::numeric_limits<double>::max_digits10);
89 // Show all potentially significant decimal digits,
90 std::cout << std::showpoint << std::endl;
91 // and show significant trailing zeros too.
92
93 double z = 10.;
94 double r = lambert_w0(z); // Default policy for double.
95 std::cout << "lambert_w0(z) = " << r << std::endl;
96 // lambert_w0(z) = 1.7455280027406994
97 //] [/lambert_w_simple_examples_0]
98 }
99 {
100 // Other floating-point types can be used too, here float.
101 // It is convenient to use a function like `show_value`
102 // to display all potentially significant decimal digits
103 // for the type, including any significant trailing zeros.
104 //[lambert_w_simple_examples_1
105 float z = 10.F;
106 float r;
107 r = lambert_w0(z); // Default policy digits10 = 7, digits2 = 24
108 std::cout << "lambert_w0(";
109 show_value(z);
110 std::cout << ") = ";
111 show_value(r);
112 std::cout << std::endl; // lambert_w0(10.0000000) = 1.74552798
113 //] //[/lambert_w_simple_examples_1]
114 }
115 {
116 // Example of an integer argument to lambert_w,
117 // showing that an integer is correctly promoted to a double.
118 //[lambert_w_simple_examples_2
119 std::cout.precision(std::numeric_limits<double>::max_digits10);
120 double r = lambert_w0(10); // Pass an int argument "10" that should be promoted to double argument.
121 std::cout << "lambert_w0(10) = " << r << std::endl; // lambert_w0(10) = 1.7455280027406994
122 double rp = lambert_w0(10);
123 std::cout << "lambert_w0(10) = " << rp << std::endl;
124 // lambert_w0(10) = 1.7455280027406994
125 auto rr = lambert_w0(10); // C++11 needed.
126 std::cout << "lambert_w0(10) = " << rr << std::endl;
127 // lambert_w0(10) = 1.7455280027406994 too, showing that rr has been promoted to double.
128 //] //[/lambert_w_simple_examples_2]
129 }
130 {
131 // Using multiprecision types to get much higher precision is painless.
132 //[lambert_w_simple_examples_3
133 cpp_dec_float_50 z("10");
134 // Note construction using a decimal digit string "10",
135 // NOT a floating-point double literal 10.
136 cpp_dec_float_50 r;
137 r = lambert_w0(z);
138 std::cout << "lambert_w0("; show_value(z); std::cout << ") = ";
139 show_value(r);
140 std::cout << std::endl;
141 // lambert_w0(10.000000000000000000000000000000000000000000000000000000000000000000000000000000) =
142 // 1.7455280027406993830743012648753899115352881290809413313533156980404446940000000
143 //] //[/lambert_w_simple_examples_3]
144 }
145 // Using multiprecision types to get multiprecision precision wrong!
146 {
147 //[lambert_w_simple_examples_4
148 cpp_dec_float_50 z(0.7777777777777777777777777777777777777777777777777777777777777777777777777);
149 // Compiler evaluates the nearest double-precision binary representation,
150 // from the max_digits10 of the floating_point literal double 0.7777777777777777777777777777...,
151 // so any extra digits in the multiprecision type
152 // beyond max_digits10 (usually 17) are random and meaningless.
153 cpp_dec_float_50 r;
154 r = lambert_w0(z);
155 std::cout << "lambert_w0(";
156 show_value(z);
157 std::cout << ") = "; show_value(r);
158 std::cout << std::endl;
159 // lambert_w0(0.77777777777777779011358916250173933804035186767578125000000000000000000000000000)
160 // = 0.48086152073210493501934682309060873341910109230469724725005039758139532631901386
161 //] //[/lambert_w_simple_examples_4]
162 }
163 {
164 //[lambert_w_simple_examples_4a
165 cpp_dec_float_50 z(0.9); // Construct from floating_point literal double 0.9.
166 cpp_dec_float_50 r;
167 r = lambert_w0(0.9);
168 std::cout << "lambert_w0(";
169 show_value(z);
170 std::cout << ") = "; show_value(r);
171 std::cout << std::endl;
172 // lambert_w0(0.90000000000000002220446049250313080847263336181640625000000000000000000000000000)
173 // = 0.52983296563343440510607251781038939952850341796875000000000000000000000000000000
174 std::cout << "lambert_w0(0.9) = " << lambert_w0(static_cast<double>(0.9))
175 // lambert_w0(0.9)
176 // = 0.52983296563343441
177 << std::endl;
178 //] //[/lambert_w_simple_examples_4a]
179 }
180 {
181 // Using multiprecision types to get multiprecision precision right!
182 //[lambert_w_simple_examples_4b
183 cpp_dec_float_50 z("0.9"); // Construct from decimal digit string.
184 cpp_dec_float_50 r;
185 r = lambert_w0(z);
186 std::cout << "lambert_w0(";
187 show_value(z);
188 std::cout << ") = "; show_value(r);
189 std::cout << std::endl;
190 // 0.90000000000000000000000000000000000000000000000000000000000000000000000000000000)
191 // = 0.52983296563343441213336643954546304857788132269804249284012528304239956413801252
192 //] //[/lambert_w_simple_examples_4b]
193 }
194 // Getting extreme precision (1000 decimal digits) Lambert W values.
195 {
196 std::cout.precision(std::numeric_limits<cpp_dec_float_1000>::digits10);
197 cpp_dec_float_1000 z("2.0");
198 cpp_dec_float_1000 r;
199 r = lambert_w0(z);
200 std::cout << "lambert_w0(z) = " << r << std::endl;
201 // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
202 // Wolfram alpha command N[productlog[0, 2.0],1000] gives the identical result:
203 // 0.8526055020137254913464724146953174668984533001514035087721073946525150656742630448965773783502494847334503972691804119834761668851953598826198984364998343940330324849743119327028383008883133161249045727544669202220292076639777316648311871183719040610274221013237163543451621208284315007250267190731048119566857455987975973474411544571619699938899354169616378479326962044241495398851839432070255805880208619490399218130868317114428351234208216131218024303904457925834743326836272959669122797896855064630871955955318383064292191644322931561534814178034773896739684452724587331245831001449498844495771266728242975586931792421997636537572767708722190588748148949667744956650966402600446780664924889043543203483210769017254907808218556111831854276511280553252641907484685164978750601216344998778097446525021666473925144772131644151718261199915247932015387685261438125313159125475113124470774926288823525823567568542843625471594347837868505309329628014463491611881381186810879712667681285740515197493390563
204 }
205 {
206 //[lambert_w_simple_examples_error_policies
207 // Define an error handling policy:
208 typedef policy<
209 domain_error<throw_on_error>,
210 overflow_error<ignore_error> // possibly unwise?
211 > my_throw_policy;
212
213 std::cout.precision(std::numeric_limits<double>::max_digits10);
214 // Show all potentially significant decimal digits,
215 std::cout << std::showpoint << std::endl;
216 // and show significant trailing zeros too.
217 double z = +1;
218 std::cout << "Lambert W (" << z << ") = " << lambert_w0(z) << std::endl;
219 // Lambert W (1.0000000000000000) = 0.56714329040978384
220 std::cout << "\nLambert W (" << z << ", my_throw_policy()) = "
221 << lambert_w0(z, my_throw_policy()) << std::endl;
222 // Lambert W (1.0000000000000000, my_throw_policy()) = 0.56714329040978384
223 //] //[/lambert_w_simple_example_error_policies]
224 }
225 {
226 // Show error reporting from passing a value to lambert_wm1 that is out of range,
227 // (and probably was meant to be passed to lambert_0 instead).
228 //[lambert_w_simple_examples_out_of_range
229 double z = +1.;
230 double r = lambert_wm1(z);
231 std::cout << "lambert_wm1(+1.) = " << r << std::endl;
232 //] [/lambert_w_simple_examples_out_of_range]
233 // Error in function boost::math::lambert_wm1<RealType>(<RealType>):
234 // Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
235 }
236 }
237 catch (std::exception& ex)
238 {
239 std::cout << ex.what() << std::endl;
240 }
241 } // int main()
242
243 /*
244
245 Output:
246 //[lambert_w_simple_examples_error_message_1
247 Error in function boost::math::lambert_wm1<RealType>(<RealType>):
248 Argument z = 1 is out of range (z <= 0) for Lambert W-1 branch! (Try Lambert W0 branch?)
249 //] [/lambert_w_simple_examples_error_message_1]
250
251 */
252