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