1 // Copyright Paul A. Bristow 2016.
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 // Test that can build and run a simple example of Lambert W function,
8 // using algorithm of Thomas Luu.
9 // https://svn.boost.org/trac/boost/ticket/11027
10
11 #include <boost/config.hpp> // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
12 #include <boost/version.hpp> // for BOOST_MSVC versions.
13 #include <boost/cstdint.hpp>
14 #include <boost/exception/exception.hpp> // boost::exception
15 #include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
16
17 #define BOOST_MATH_INSTRUMENT_LAMBERT_W // #define only for diagnostic output.
18
19 // For lambert_w function.
20 #include <boost/math/special_functions/lambert_w.hpp>
21
22 #include <iostream>
23 // using std::cout;
24 // using std::endl;
25 #include <exception>
26 #include <stdexcept>
27 #include <string>
28 #include <limits> // For std::numeric_limits.
29
30 //! Show information about build, architecture, address model, platform, ...
show_versions()31 std::string show_versions()
32 {
33 std::ostringstream message;
34
35 message << "Program: " << __FILE__ << "\n";
36 #ifdef __TIMESTAMP__
37 message << __TIMESTAMP__;
38 #endif
39 message << "\nBuildInfo:\n" " Platform " << BOOST_PLATFORM;
40 // http://stackoverflow.com/questions/1505582/determining-32-vs-64-bit-in-c
41 #if defined(__LP64__) || defined(_WIN64) || (defined(__x86_64__) && !defined(__ILP32__) ) || defined(_M_X64) || defined(__ia64) || defined (_M_IA64) || defined(__aarch64__) || defined(__powerpc64__)
42 #define IS64BIT 1
43 message << ", 64-bit.";
44 #else
45 #define IS32BIT 1
46 message << ", 32-bit.";
47 #endif
48
49 message << "\n Compiler " BOOST_COMPILER;
50 #ifdef BOOST_MSC_VER
51 #ifdef _MSC_FULL_VER
52 message << "\n MSVC version " << BOOST_STRINGIZE(_MSC_FULL_VER) << ".";
53 #endif
54 #ifdef __WIN64
55 mess age << "\n WIN64" << std::endl;
56 #endif // __WIN64
57 #ifdef _WIN32
58 message << "\n WIN32" << std::endl;
59 #endif // __WIN32
60 #endif
61 #ifdef __GNUC__
62 //PRINT_MACRO(__GNUC__);
63 //PRINT_MACRO(__GNUC_MINOR__);
64 //PRINT_MACRO(__GNUC_PATCH__);
65 std::cout << "GCC " << __VERSION__ << std::endl;
66 //PRINT_MACRO(LONG_MAX);
67 #endif // __GNUC__
68
69 message << "\n STL " << BOOST_STDLIB;
70
71 message << "\n Boost version " << BOOST_VERSION / 100000 << "." << BOOST_VERSION / 100 % 1000 << "." << BOOST_VERSION % 100;
72
73 #ifdef BOOST_HAS_FLOAT128
74 message << ", BOOST_HAS_FLOAT128" << std::endl;
75 #endif
76 message << std::endl;
77 return message.str();
78 } // std::string versions()
79
main()80 int main()
81 {
82 try
83 {
84 //std::cout << "Lambert W example basic!" << std::endl;
85 //std::cout << show_versions() << std::endl;
86
87 //std::cout << exp(1) << std::endl; // 2.71828
88 //std::cout << exp(-1) << std::endl; // 0.367879
89 //std::cout << std::numeric_limits<double>::epsilon() / 2 << std::endl; // 1.11022e-16
90
91 using namespace boost::math;
92 using boost::math::constants::exp_minus_one;
93 double x = 1.;
94
95 double W1 = lambert_w(1.);
96 // Note, NOT integer X, for example: lambert_w(1); or will get message like
97 // error C2338: Must be floating-point, not integer type, for example W(1.), not W(1)!
98 //
99
100 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.567143
101 // This 'golden ratio' for exponentials is http://mathworld.wolfram.com/OmegaConstant.html
102 // since exp[-W(1)] = W(1)
103 // A030178 Decimal expansion of LambertW(1): the solution to x*exp(x)
104 // = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866
105 // http://oeis.org/A030178
106
107 double expplogone = exp(-lambert_w(1.));
108 if (expplogone != W1)
109 {
110 std::cout << expplogone << " " << W1 << std::endl; //
111 }
112
113
114 //[lambert_w_example_1
115
116 x = 0.01;
117 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.00990147
118 //] [/lambert_w_example_1]
119 x = -0.01;
120 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // -0.0101015
121 x = -0.1;
122 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
123 /**/
124
125 for (double xd = 1.; xd < 1e20; xd *= 10)
126 {
127
128 // 1. 0.56714329040978387
129 // 0.56714329040978384
130
131 // 10 1.7455280027406994
132 // 1.7455280027406994
133
134 // 100 3.3856301402900502
135 // 3.3856301402900502
136 // 1000 5.2496028524015959
137 // 5.249602852401596227126056319697306282521472386059592844451465483991362228320942832739693150854347718
138
139 // 1e19 40.058769161984308
140 // 40.05876916198431163898797971203180915622644925765346546858291325452428038208071849105889199253335063
141 std::cout << "Lambert W (" << xd << ") = " << lambert_w(xd) << std::endl; //
142 }
143 //
144 // Test near singularity.
145
146 // http://www.wolframalpha.com/input/?i=N%5Blambert_w%5B-0.367879%5D,17%5D test value N[lambert_w[-0.367879],17]
147 // -0.367879441171442321595523770161460867445811131031767834
148 x = -0.367879; // < -exp(1) = -0.367879
149 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // Lambert W (-0.36787900000000001) = -0.99845210378080340
150 // -0.99845210378080340
151 // -0.99845210378072726 N[lambert_w[-0.367879],17] wolfram so very close.
152
153 x = -0.3678794; // expect -0.99952696660756813
154 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
155 x = -0.36787944; // expect -0.99992019848408340
156 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
157 x = -0.367879441; // -0.99996947070054883
158 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
159 x = -0.36787944117; // -0.99999719977527159
160 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
161 x = -0.367879441171; // -0.99999844928821992
162 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
163
164 x = -exp_minus_one<double>() + std::numeric_limits<double>::epsilon();
165 // Lambert W (-0.36787944117144211) = -0.99999996349975895
166 // N[lambert_w[-0.36787944117144211],17] == -0.99999996608315303
167 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
168 std::cout << " 1 - sqrt(eps) = " << static_cast<double>(1) - sqrt(std::numeric_limits<double>::epsilon()) << std::endl;
169 x = -exp_minus_one<double>();
170 // N[lambert_w[-0.36787944117144233],17] == -1.000000000000000 + 6.7595465843924897*10^-9i
171 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // 0.0
172 // At Singularity - 0.36787944117144233 == -0.36787944117144233 returned - 1.0000000000000000
173 // Lambert W(-0.36787944117144233) = -1.0000000000000000
174
175
176 x = (std::numeric_limits<double>::max)()/4;
177 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; // OK 702.023799146706
178 x = (std::numeric_limits<double>::max)()/2;
179 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
180 x = (std::numeric_limits<double>::max)();
181 std::cout << "Lambert W (" << x << ") = " << lambert_w(x) << std::endl; //
182 // Error in function boost::math::log1p<double>(double): numeric overflow
183 /* */
184
185 }
186 catch (std::exception& ex)
187 {
188 std::cout << ex.what() << std::endl;
189 }
190
191
192 } // int main()
193
194 /*
195
196 //[lambert_w_output_1
197 Output:
198
199 1> example_basic.cpp
200 1> Generating code
201 1> All 237 functions were compiled because no usable IPDB/IOBJ from previous compilation was found.
202 1> Finished generating code
203 1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.exe
204 1> LambertW.vcxproj -> J:\Cpp\Misc\x64\Release\LambertW.pdb (Full PDB)
205 1> Lambert W example basic!
206 1> Platform: Win32
207 1> Compiler: Microsoft Visual C++ version 14.0
208 1> STL : Dinkumware standard library version 650
209 1> Boost : 1.63.0
210 1> _MSC_FULL_VER = 190024123
211 1> Win32
212 1> x64
213 1> (x64)
214 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
215 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
216 1> Final 0.567143290409784 after 2 iterations, difference = 0
217 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
218 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
219 1> Final 0.567143290409784 after 2 iterations, difference = 0
220 1> Lambert W (1) = 0.567143290409784
221 1> Iteration #0, w0 0.577547206058041, w1 = 0.567143616915443, difference = 0.0289944962755619, relative 0.018343835374856
222 1> Iteration #1, w0 0.567143616915443, w1 = 0.567143290409784, difference = 9.02208135089566e-07, relative 5.75702234328901e-07
223 1> Final 0.567143290409784 after 2 iterations, difference = 0
224 1> Iteration #0, w0 0.0099072820916067, w1 = 0.00990147384359511, difference = 5.92416060777624e-06, relative 0.000586604388734591
225 1> Final 0.00990147384359511 after 1 iterations, difference = 0
226 1> Lambert W (0.01) = 0.00990147384359511
227 1> Iteration #0, w0 -0.0101016472705154, w1 = -0.0101015271985388, difference = -1.17664437923951e-07, relative 1.18865171889748e-05
228 1> Final -0.0101015271985388 after 1 iterations, difference = 0
229 1> Lambert W (-0.01) = -0.0101015271985388
230 1> Iteration #0, w0 -0.111843322610692, w1 = -0.111832559158964, difference = -8.54817065376601e-06, relative 9.62461362694622e-05
231 1> Iteration #1, w0 -0.111832559158964, w1 = -0.111832559158963, difference = -5.68989300120393e-16, relative 6.43929354282591e-15
232 1> Final -0.111832559158963 after 2 iterations, difference = 0
233 1> Lambert W (-0.1) = -0.111832559158963
234 1> Iteration #0, w0 -0.998452103785573, w1 = -0.998452103780803, difference = -2.72004641033163e-15, relative 4.77662354114727e-12
235 1> Final -0.998452103780803 after 1 iterations, difference = 0
236 1> Lambert W (-0.367879) = -0.998452103780803
237
238 //] [/lambert_w_output_1]
239 */
240