1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
|
// 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]
*/
|
