1 // Copyright Paul A. Bristow 2016, 2017, 2018.
2 // Copyright John Maddock 2016.
3
4 // Use, modification and distribution are subject to the
5 // Boost Software License, Version 1.0.
6 // (See accompanying file LICENSE_1_0.txt
7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
8
9 // test_lambert_w_integrals.cpp
10 //! \brief quadrature tests that cover the whole range of the Lambert W0 function.
11
12 #include <boost/config.hpp> // for BOOST_MSVC definition etc.
13 #include <boost/version.hpp> // for BOOST_MSVC versions.
14
15 #ifdef BOOST_HAS_FLOAT128
16
17 // Boost macros
18 #define BOOST_TEST_MAIN
19 #define BOOST_LIB_DIAGNOSTIC "on" // Report library file details.
20 #include <boost/test/included/unit_test.hpp> // Boost.Test
21 #include <boost/test/tools/floating_point_comparison.hpp>
22
23 #include <boost/array.hpp>
24 #include <boost/lexical_cast.hpp>
25 #include <boost/type_traits/is_constructible.hpp>
26
27 #include <boost/multiprecision/float128.hpp>
28
29 #include <boost/math/special_functions/fpclassify.hpp> // isnan, isfinite.
30 #include <boost/math/special_functions/next.hpp> // float_next, float_prior
31 using boost::math::float_next;
32 using boost::math::float_prior;
33 #include <boost/math/special_functions/ulp.hpp> // ulp
34
35 #include <boost/math/tools/test_value.hpp> // for create_test_value and macro BOOST_MATH_TEST_VALUE.
36 #include <boost/math/policies/policy.hpp>
37 using boost::math::policies::digits2;
38 using boost::math::policies::digits10;
39 #include <boost/math/special_functions/lambert_w.hpp> // For Lambert W lambert_w function.
40 using boost::math::lambert_wm1;
41 using boost::math::lambert_w0;
42
43 #include <limits>
44 #include <cmath>
45 #include <typeinfo>
46 #include <iostream>
47 #include <type_traits>
48 #include <exception>
49
50 std::string show_versions(void);
51
52 // Added code and test for Integral of the Lambert W function: by Nick Thompson.
53 // https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals
54
55 #include <boost/math/constants/constants.hpp> // for integral tests.
56 #include <boost/math/quadrature/tanh_sinh.hpp> // for integral tests.
57 #include <boost/math/quadrature/exp_sinh.hpp> // for integral tests.
58
59 using boost::math::policies::policy;
60 using boost::math::policies::make_policy;
61
62 // using statements needed for changing error handling policy.
63 using boost::math::policies::evaluation_error;
64 using boost::math::policies::domain_error;
65 using boost::math::policies::overflow_error;
66 using boost::math::policies::ignore_error;
67 using boost::math::policies::throw_on_error;
68
69 typedef policy<
70 domain_error<throw_on_error>,
71 overflow_error<ignore_error>
72 > no_throw_policy;
73
74 // Assumes that function has a throw policy, for example:
75 // NOT lambert_w0<T>(1 / (x * x), no_throw_policy());
76 // Error in function boost::math::quadrature::exp_sinh<double>::integrate:
77 // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf.
78 // Please ensure your function evaluates to a finite number of its entire domain.
79 template <typename T>
debug_integration_proc(T x)80 T debug_integration_proc(T x)
81 {
82 T result; // warning C4701: potentially uninitialized local variable 'result' used
83 // T result = 0 ; // But result may not be assigned below?
84 try
85 {
86 // Assign function call to result in here...
87 if (x <= sqrt(boost::math::tools::min_value<T>()) )
88 {
89 result = 0;
90 }
91 else
92 {
93 result = lambert_w0<T>(1 / (x * x));
94 }
95 // result = lambert_w0<T>(1 / (x * x), no_throw_policy()); // Bad idea, less helpful diagnostic message is:
96 // Error in function boost::math::quadrature::exp_sinh<double>::integrate:
97 // The exp_sinh quadrature evaluated your function at a singular point and resulted in inf.
98 // Please ensure your function evaluates to a finite number of its entire domain.
99
100 } // try
101 catch (const std::exception& e)
102 {
103 std::cout << "Exception " << e.what() << std::endl;
104 // set breakpoint here:
105 std::cout << "Unexpected exception thrown in integration code at abscissa (x): " << x << "." << std::endl;
106 if (!std::isfinite(result))
107 {
108 // set breakpoint here:
109 std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl;
110 }
111 if (std::isnan(result))
112 {
113 // set breakpoint here:
114 std::cout << "Unexpected non-finite result in integration code at abscissa (x): " << x << "." << std::endl;
115 }
116 } // catch
117 return result;
118 } // T debug_integration_proc(T x)
119
120 template<class Real>
test_integrals()121 void test_integrals()
122 {
123 // Integral of the Lambert W function:
124 // https://en.wikipedia.org/wiki/Lambert_W_function
125 using boost::math::quadrature::tanh_sinh;
126 using boost::math::quadrature::exp_sinh;
127 // file:///I:/modular-boost/libs/math/doc/html/math_toolkit/quadrature/double_exponential/de_tanh_sinh.html
128 using std::sqrt;
129
130 std::cout << "Integration of type " << typeid(Real).name() << std::endl;
131
132 Real tol = std::numeric_limits<Real>::epsilon();
133 { // // Integrate for function lambert_W0(z);
134 tanh_sinh<Real> ts;
135 Real a = 0;
136 Real b = boost::math::constants::e<Real>();
137 auto f = [](Real z)->Real
138 {
139 return lambert_w0<Real>(z);
140 };
141 Real z = ts.integrate(f, a, b); // OK without any decltype(f)
142 BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::e<Real>() - 1, tol);
143 }
144 {
145 // Integrate for function lambert_W0(z/(z sqrt(z)).
146 exp_sinh<Real> es;
147 auto f = [](Real z)->Real
148 {
149 return lambert_w0<Real>(z)/(z * sqrt(z));
150 };
151 Real z = es.integrate(f); // OK
152 BOOST_CHECK_CLOSE_FRACTION(z, 2 * boost::math::constants::root_two_pi<Real>(), tol);
153 }
154 {
155 // Integrate for function lambert_W0(1/z^2).
156 exp_sinh<Real> es;
157 //const Real sqrt_min = sqrt(boost::math::tools::min_value<Real>()); // 1.08420217e-19 fo 32-bit float.
158 // error C3493: 'sqrt_min' cannot be implicitly captured because no default capture mode has been specified
159 auto f = [](Real z)->Real
160 {
161 if (z <= sqrt(boost::math::tools::min_value<Real>()) )
162 { // Too small would underflow z * z and divide by zero to overflow 1/z^2 for lambert_w0 z parameter.
163 return static_cast<Real>(0);
164 }
165 else
166 {
167 return lambert_w0<Real>(1 / (z * z)); // warning C4756: overflow in constant arithmetic, even though cannot happen.
168 }
169 };
170 Real z = es.integrate(f);
171 BOOST_CHECK_CLOSE_FRACTION(z, boost::math::constants::root_two_pi<Real>(), tol);
172 }
173 } // template<class Real> void test_integrals()
174
175
BOOST_AUTO_TEST_CASE(integrals)176 BOOST_AUTO_TEST_CASE( integrals )
177 {
178 std::cout << "Macro BOOST_MATH_LAMBERT_W0_INTEGRALS is defined." << std::endl;
179 BOOST_TEST_MESSAGE("\nTest Lambert W0 integrals.");
180 try
181 {
182 // using statements needed to change precision policy.
183 using boost::math::policies::policy;
184 using boost::math::policies::make_policy;
185 using boost::math::policies::precision;
186 using boost::math::policies::digits2;
187 using boost::math::policies::digits10;
188
189 // using statements needed for changing error handling policy.
190 using boost::math::policies::evaluation_error;
191 using boost::math::policies::domain_error;
192 using boost::math::policies::overflow_error;
193 using boost::math::policies::ignore_error;
194 using boost::math::policies::throw_on_error;
195
196 /*
197 typedef policy<
198 domain_error<throw_on_error>,
199 overflow_error<ignore_error>
200 > no_throw_policy;
201
202
203 // Experiment with better diagnostics.
204 typedef float Real;
205
206 Real inf = std::numeric_limits<Real>::infinity();
207 Real max = (std::numeric_limits<Real>::max)();
208 std::cout.precision(std::numeric_limits<Real>::max_digits10);
209 //std::cout << "lambert_w0(inf) = " << lambert_w0(inf) << std::endl; // lambert_w0(inf) = 1.79769e+308
210 std::cout << "lambert_w0(inf, throw_policy()) = " << lambert_w0(inf, no_throw_policy()) << std::endl; // inf
211 std::cout << "lambert_w0(max) = " << lambert_w0(max) << std::endl; // lambert_w0(max) = 703.227
212 //std::cout << lambert_w0(inf) << std::endl; // inf - will throw.
213 std::cout << "lambert_w0(0) = " << lambert_w0(0.) << std::endl; // 0
214 std::cout << "lambert_w0(std::numeric_limits<Real>::denorm_min()) = " << lambert_w0(std::numeric_limits<Real>::denorm_min()) << std::endl; // 4.94066e-324
215 std::cout << "lambert_w0(std::numeric_limits<Real>::min()) = " << lambert_w0((std::numeric_limits<Real>::min)()) << std::endl; // 2.22507e-308
216
217 // Approximate the largest lambert_w you can get for type T?
218 float max_w_f = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits<float>::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162.
219 std::cout << "w max_f " << max_w_f << std::endl; // 84.2879
220 Real max_w = boost::math::lambert_w_detail::lambert_w0_approx((std::numeric_limits<Real>::max)()); // Corless equation 4.19, page 349, and Chapeau-Blondeau equation 20, page 2162.
221 std::cout << "w max " << max_w << std::endl; // 703.227
222
223 std::cout << "lambert_w0(7.2416706213544837e-163) = " << lambert_w0(7.2416706213544837e-163) << std::endl; //
224 std::cout << "test integral 1/z^2" << std::endl;
225 std::cout << "ULP = " << boost::math::ulp(1., policy<digits2<> >()) << std::endl; // ULP = 2.2204460492503131e-16
226 std::cout << "ULP = " << boost::math::ulp(1e-10, policy<digits2<> >()) << std::endl; // ULP = 2.2204460492503131e-16
227 std::cout << "ULP = " << boost::math::ulp(1., policy<digits2<11> >()) << std::endl; // ULP = 2.2204460492503131e-16
228 std::cout << "epsilon = " << std::numeric_limits<Real>::epsilon() << std::endl; //
229 std::cout << "sqrt(max) = " << sqrt(boost::math::tools::max_value<float>() ) << std::endl; // sqrt(max) = 1.8446742974197924e+19
230 std::cout << "sqrt(min) = " << sqrt(boost::math::tools::min_value<float>() ) << std::endl; // sqrt(min) = 1.0842021724855044e-19
231
232
233
234 // Demo debug version.
235 Real tol = std::numeric_limits<Real>::epsilon();
236 Real x;
237 {
238 using boost::math::quadrature::exp_sinh;
239 exp_sinh<Real> es;
240 // Function to be integrated, lambert_w0(1/z^2).
241
242 //auto f = [](Real z)->Real
243 //{ // Naive - no protection against underflow and subsequent divide by zero.
244 // return lambert_w0<Real>(1 / (z * z));
245 //};
246 // Diagnostic is:
247 // Error in function boost::math::lambert_w0<Real>: Expected a finite value but got inf
248
249 auto f = [](Real z)->Real
250 { // Debug with diagnostics for underflow and subsequent divide by zero and other bad things.
251 return debug_integration_proc(z);
252 };
253 // Exception Error in function boost::math::lambert_w0<double>: Expected a finite value but got inf.
254
255 // Unexpected exception thrown in integration code at abscissa: 7.2416706213544837e-163.
256 // Unexpected exception thrown in integration code at abscissa (x): 3.478765835953569e-23.
257 x = es.integrate(f);
258 std::cout << "es.integrate(f) = " << x << std::endl;
259 BOOST_CHECK_CLOSE_FRACTION(x, boost::math::constants::root_two_pi<Real>(), tol);
260 // root_two_pi<double = 2.506628274631000502
261 }
262 */
263
264 test_integrals<boost::multiprecision::float128>();
265 }
266 catch (std::exception& ex)
267 {
268 std::cout << ex.what() << std::endl;
269 }
270 }
271
272 #else
273
main()274 int main() { return 0; }
275
276 #endif
277