]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Copyright Nick Thompson, 2017 |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. | |
4 | // (See accompanying file LICENSE_1_0.txt | |
5 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
7 | #define BOOST_TEST_MODULE exp_sinh_quadrature_test | |
8 | ||
92f5a8d4 | 9 | #include <complex> |
1e59de90 TL |
10 | #include <type_traits> |
11 | #include <boost/math/tools/config.hpp> | |
12 | #include <boost/math/tools/test_value.hpp> | |
92f5a8d4 | 13 | #include <boost/multiprecision/cpp_complex.hpp> |
b32b8144 FG |
14 | #include <boost/math/concepts/real_concept.hpp> |
15 | #include <boost/test/included/unit_test.hpp> | |
92f5a8d4 | 16 | #include <boost/test/tools/floating_point_comparison.hpp> |
b32b8144 FG |
17 | #include <boost/math/quadrature/exp_sinh.hpp> |
18 | #include <boost/math/special_functions/sinc.hpp> | |
19 | #include <boost/math/special_functions/bessel.hpp> | |
20 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
21 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
22 | #include <boost/math/special_functions/next.hpp> | |
23 | #include <boost/math/special_functions/gamma.hpp> | |
24 | #include <boost/math/special_functions/sinc.hpp> | |
25 | #include <boost/type_traits/is_class.hpp> | |
26 | ||
92f5a8d4 TL |
27 | #ifdef BOOST_HAS_FLOAT128 |
28 | #include <boost/multiprecision/complex128.hpp> | |
29 | #endif | |
30 | ||
b32b8144 FG |
31 | using std::exp; |
32 | using std::cos; | |
33 | using std::tan; | |
34 | using std::log; | |
35 | using std::sqrt; | |
36 | using std::abs; | |
37 | using std::sinh; | |
38 | using std::cosh; | |
39 | using std::pow; | |
40 | using std::atan; | |
41 | using boost::multiprecision::cpp_bin_float_50; | |
42 | using boost::multiprecision::cpp_bin_float_100; | |
43 | using boost::multiprecision::cpp_bin_float_quad; | |
44 | using boost::math::constants::pi; | |
45 | using boost::math::constants::half_pi; | |
46 | using boost::math::constants::two_div_pi; | |
47 | using boost::math::constants::half; | |
48 | using boost::math::constants::third; | |
49 | using boost::math::constants::half; | |
50 | using boost::math::constants::third; | |
51 | using boost::math::constants::catalan; | |
52 | using boost::math::constants::ln_two; | |
53 | using boost::math::constants::root_two; | |
54 | using boost::math::constants::root_two_pi; | |
55 | using boost::math::constants::root_pi; | |
56 | using boost::math::quadrature::exp_sinh; | |
57 | ||
1e59de90 | 58 | #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8) && !defined(TEST9) && !defined(TEST10) |
b32b8144 FG |
59 | # define TEST1 |
60 | # define TEST2 | |
61 | # define TEST3 | |
62 | # define TEST4 | |
63 | # define TEST5 | |
64 | # define TEST6 | |
65 | # define TEST7 | |
92f5a8d4 | 66 | # define TEST8 |
f67539c2 | 67 | # define TEST9 |
1e59de90 | 68 | # define TEST10 |
b32b8144 FG |
69 | #endif |
70 | ||
1e59de90 | 71 | #ifdef _MSC_VER |
b32b8144 FG |
72 | #pragma warning (disable:4127) |
73 | #endif | |
74 | ||
75 | // | |
76 | // Coefficient generation code: | |
77 | // | |
78 | template <class T> | |
79 | void print_levels(const T& v, const char* suffix) | |
80 | { | |
81 | std::cout << "{\n"; | |
82 | for (unsigned i = 0; i < v.size(); ++i) | |
83 | { | |
84 | std::cout << " { "; | |
85 | for (unsigned j = 0; j < v[i].size(); ++j) | |
86 | { | |
87 | std::cout << v[i][j] << suffix << ", "; | |
88 | } | |
89 | std::cout << "},\n"; | |
90 | } | |
91 | std::cout << " };\n"; | |
92 | } | |
93 | ||
94 | template <class T> | |
95 | void print_levels(const std::pair<T, T>& p, const char* suffix = "") | |
96 | { | |
97 | std::cout << " static const std::vector<std::vector<Real> > abscissa = "; | |
98 | print_levels(p.first, suffix); | |
99 | std::cout << " static const std::vector<std::vector<Real> > weights = "; | |
100 | print_levels(p.second, suffix); | |
101 | } | |
102 | ||
103 | template <class Real, class TargetType> | |
104 | std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows) | |
105 | { | |
106 | using boost::math::constants::half_pi; | |
107 | using boost::math::constants::two_div_pi; | |
108 | using boost::math::constants::pi; | |
92f5a8d4 TL |
109 | auto g = [](Real t)->Real { return exp(half_pi<Real>()*sinh(t)); }; |
110 | auto w = [](Real t)->Real { return cosh(t)*half_pi<Real>()*exp(half_pi<Real>()*sinh(t)); }; | |
b32b8144 FG |
111 | |
112 | std::vector<std::vector<Real>> abscissa, weights; | |
113 | ||
114 | std::vector<Real> temp; | |
115 | ||
116 | Real tmp = (Real(boost::math::tools::log_min_value<TargetType>()) + log(Real(boost::math::tools::epsilon<TargetType>())))*0.5f; | |
117 | Real t_min = asinh(two_div_pi<Real>()*tmp); | |
118 | // truncate t_min to an exact binary value: | |
119 | t_min = floor(t_min * 128) / 128; | |
120 | ||
121 | std::cout << "m_t_min = " << t_min << ";\n"; | |
122 | ||
123 | // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g). | |
124 | // This will allow some flexibility on the users part; they can at least square a number function without overflow. | |
125 | // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands. | |
126 | const Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(Real(boost::math::tools::max_value<TargetType>())))); | |
127 | ||
128 | Real h = 1; | |
129 | for (Real t = t_min; t < t_max; t += h) | |
130 | { | |
131 | temp.push_back(g(t)); | |
132 | } | |
133 | abscissa.push_back(temp); | |
134 | temp.clear(); | |
135 | ||
136 | for (Real t = t_min; t < t_max; t += h) | |
137 | { | |
138 | temp.push_back(w(t * h)); | |
139 | } | |
140 | weights.push_back(temp); | |
141 | temp.clear(); | |
142 | ||
143 | for (unsigned row = 1; row < max_rows; ++row) | |
144 | { | |
145 | h /= 2; | |
146 | for (Real t = t_min + h; t < t_max; t += 2 * h) | |
147 | temp.push_back(g(t)); | |
148 | abscissa.push_back(temp); | |
149 | temp.clear(); | |
150 | } | |
151 | h = 1; | |
152 | for (unsigned row = 1; row < max_rows; ++row) | |
153 | { | |
154 | h /= 2; | |
155 | for (Real t = t_min + h; t < t_max; t += 2 * h) | |
156 | temp.push_back(w(t)); | |
157 | weights.push_back(temp); | |
158 | temp.clear(); | |
159 | } | |
160 | ||
161 | return std::make_pair(abscissa, weights); | |
162 | } | |
163 | ||
164 | ||
165 | template <class Real> | |
166 | const exp_sinh<Real>& get_integrator() | |
167 | { | |
168 | static const exp_sinh<Real> integrator(14); | |
169 | return integrator; | |
170 | } | |
171 | ||
172 | template <class Real> | |
173 | Real get_convergence_tolerance() | |
174 | { | |
175 | return boost::math::tools::root_epsilon<Real>(); | |
176 | } | |
177 | ||
178 | template<class Real> | |
179 | void test_right_limit_infinite() | |
180 | { | |
181 | std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
182 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
183 | Real Q; | |
184 | Real Q_expected; | |
185 | Real error; | |
186 | Real L1; | |
187 | auto integrator = get_integrator<Real>(); | |
188 | ||
189 | // Example 12 | |
92f5a8d4 | 190 | const auto f2 = [](const Real& t)->Real { return exp(-t)/sqrt(t); }; |
b32b8144 FG |
191 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); |
192 | Q_expected = root_pi<Real>(); | |
193 | Real tol_mult = 1; | |
194 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
195 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10) | |
196 | tol_mult = 12; | |
197 | else if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10) | |
198 | tol_mult = 5; | |
199 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * tol_mult); | |
200 | // The integrand is strictly positive, so it coincides with the value of the integral: | |
201 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * tol_mult); | |
202 | ||
1e59de90 TL |
203 | #ifdef BOOST_MATH_STANDALONE |
204 | BOOST_IF_CONSTEXPR (std::is_fundamental<Real>::value) | |
205 | #endif | |
206 | { | |
207 | auto f3 = [](Real t)->Real { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); }; | |
208 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
209 | Q_expected = half<Real>(); | |
210 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
211 | Q = integrator.integrate(f3, 10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
212 | Q_expected = BOOST_MATH_TEST_VALUE(Real, -6.6976341310426674140007086979326069121526743314567805278252392932e-6); | |
213 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol); | |
214 | // Integrating through zero risks precision loss: | |
215 | Q = integrator.integrate(f3, -10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
216 | Q_expected = BOOST_MATH_TEST_VALUE(Real, -15232.3213626280525704332288302799653087046646639974940243044623285817777006); | |
217 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, std::numeric_limits<Real>::digits10 > 30 ? 1000 * tol : tol); | |
218 | ||
219 | auto f4 = [](Real t)->Real { return 1/(1+t*t); }; | |
220 | Q = integrator.integrate(f4, 1, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
221 | Q_expected = pi<Real>()/4; | |
222 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
223 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
224 | Q = integrator.integrate(f4, 20, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
225 | Q_expected = BOOST_MATH_TEST_VALUE(Real, 0.0499583957219427614100062870348448814912770804235071744108534548299835954767); | |
226 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
227 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
228 | Q = integrator.integrate(f4, 500, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1); | |
229 | Q_expected = BOOST_MATH_TEST_VALUE(Real, 0.0019999973333397333150476759363217553199063513829126652556286269630); | |
230 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
231 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
232 | } | |
b32b8144 FG |
233 | } |
234 | ||
235 | template<class Real> | |
236 | void test_left_limit_infinite() | |
237 | { | |
238 | std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
239 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
240 | Real Q; | |
241 | Real Q_expected; | |
242 | Real error; | |
243 | Real L1; | |
244 | auto integrator = get_integrator<Real>(); | |
245 | ||
246 | // Example 11: | |
1e59de90 TL |
247 | #ifdef BOOST_MATH_STANDALONE |
248 | BOOST_IF_CONSTEXPR (std::is_fundamental<Real>::value) | |
249 | #endif | |
250 | { | |
251 | auto f1 = [](const Real& t)->Real { return 1/(1+t*t);}; | |
252 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), 0, get_convergence_tolerance<Real>(), &error, &L1); | |
253 | Q_expected = half_pi<Real>(); | |
254 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
255 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
256 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -20, get_convergence_tolerance<Real>(), &error, &L1); | |
257 | Q_expected = BOOST_MATH_TEST_VALUE(Real, 0.0499583957219427614100062870348448814912770804235071744108534548299835954767); | |
258 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
259 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
260 | Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -500, get_convergence_tolerance<Real>(), &error, &L1); | |
261 | Q_expected = BOOST_MATH_TEST_VALUE(Real, 0.0019999973333397333150476759363217553199063513829126652556286269630); | |
262 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
263 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
264 | } | |
b32b8144 FG |
265 | } |
266 | ||
267 | ||
268 | // Some examples of tough integrals from NR, section 4.5.4: | |
269 | template<class Real> | |
270 | void test_nr_examples() | |
271 | { | |
272 | using std::sin; | |
273 | using std::cos; | |
274 | using std::pow; | |
275 | using std::exp; | |
276 | using std::sqrt; | |
277 | std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
278 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
279 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
280 | Real Q; | |
281 | Real Q_expected; | |
282 | Real L1; | |
283 | Real error; | |
284 | auto integrator = get_integrator<Real>(); | |
285 | ||
286 | auto f0 = [] (Real)->Real { return (Real) 0; }; | |
287 | Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1); | |
288 | Q_expected = 0; | |
289 | BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol); | |
290 | BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol); | |
291 | ||
92f5a8d4 | 292 | auto f = [](const Real& x)->Real { return 1/(1+x*x); }; |
b32b8144 FG |
293 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); |
294 | Q_expected = half_pi<Real>(); | |
295 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
296 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); | |
297 | ||
298 | auto f1 = [](Real x)->Real { | |
299 | Real z1 = exp(-x); | |
300 | if (z1 == 0) | |
301 | { | |
302 | return (Real) 0; | |
303 | } | |
304 | Real z2 = pow(x, -3*half<Real>())*z1; | |
305 | if (z2 == 0) | |
306 | { | |
307 | return (Real) 0; | |
308 | } | |
309 | return sin(x*half<Real>())*z2; | |
310 | }; | |
311 | ||
312 | Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1); | |
313 | Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2)); | |
314 | ||
315 | // The integrand is oscillatory; the accuracy is low. | |
316 | Real tol_mul = 1; | |
317 | if (std::numeric_limits<Real>::digits10 > 40) | |
318 | tol_mul = 500000; | |
319 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); | |
320 | ||
321 | auto f2 = [](Real x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(pow(x, -(Real) 2/(Real) 7)*exp(-x*x)); }; | |
322 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); | |
323 | Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14); | |
324 | tol_mul = 1; | |
325 | if (std::numeric_limits<Real>::is_specialized == false) | |
92f5a8d4 TL |
326 | tol_mul = 6; |
327 | else if (std::numeric_limits<Real>::digits10 > 40) | |
b32b8144 | 328 | tol_mul = 100; |
92f5a8d4 TL |
329 | else |
330 | tol_mul = 3; | |
b32b8144 FG |
331 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); |
332 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol); | |
333 | ||
334 | auto f3 = [](Real x)->Real { return (Real) 1/ (sqrt(x)*(1+x)); }; | |
335 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
336 | Q_expected = pi<Real>(); | |
337 | ||
338 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon<Real>()); | |
339 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon<Real>()); | |
340 | ||
92f5a8d4 | 341 | auto f4 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t*half<Real>())); }; |
b32b8144 FG |
342 | Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1); |
343 | Q_expected = root_two_pi<Real>()/2; | |
344 | tol_mul = 1; | |
345 | if (std::numeric_limits<Real>::digits10 > 40) | |
346 | tol_mul = 5000; | |
347 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol); | |
348 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol); | |
349 | ||
92f5a8d4 | 350 | auto f5 = [](const Real& t)->Real { return 1/cosh(t);}; |
b32b8144 FG |
351 | Q = integrator.integrate(f5, get_convergence_tolerance<Real>(), &error, &L1); |
352 | Q_expected = half_pi<Real>(); | |
353 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 12); // Fails at float precision without higher error rate | |
354 | BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 12); | |
355 | } | |
356 | ||
357 | // Definite integrals found in the CRC Handbook of Mathematical Formulas | |
358 | template<class Real> | |
359 | void test_crc() | |
360 | { | |
361 | using std::sin; | |
362 | using std::pow; | |
363 | using std::exp; | |
364 | using std::sqrt; | |
365 | using std::log; | |
366 | using std::cos; | |
367 | std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
368 | Real tol = 10 * boost::math::tools::epsilon<Real>(); | |
369 | std::cout << std::setprecision(std::numeric_limits<Real>::digits10); | |
370 | Real Q; | |
371 | Real Q_expected; | |
372 | Real L1; | |
373 | Real error; | |
374 | auto integrator = get_integrator<Real>(); | |
375 | ||
92f5a8d4 | 376 | auto f0 = [](const Real& x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(log(x)*exp(-x)); }; |
b32b8144 FG |
377 | Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1); |
378 | Q_expected = -boost::math::constants::euler<Real>(); | |
379 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
380 | ||
381 | // Test the integral representation of the gamma function: | |
382 | auto f1 = [](Real t)->Real { Real x = exp(-t); | |
383 | if(x == 0) | |
384 | { | |
385 | return (Real) 0; | |
386 | } | |
387 | return pow(t, (Real) 12 - 1)*x; | |
388 | }; | |
389 | ||
390 | Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1); | |
391 | Q_expected = boost::math::tgamma(12.0f); | |
392 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
393 | ||
394 | // Integral representation of the modified bessel function: | |
395 | // K_5(12) | |
396 | auto f2 = [](Real t)->Real { | |
397 | Real x = 12*cosh(t); | |
398 | if (x > boost::math::tools::log_max_value<Real>()) | |
399 | { | |
400 | return (Real) 0; | |
401 | } | |
402 | return exp(-x)*cosh(5*t); | |
403 | }; | |
404 | Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1); | |
405 | Q_expected = boost::math::cyl_bessel_k<int, Real>(5, 12); | |
406 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
407 | // Laplace transform of cos(at) | |
408 | Real a = 20; | |
409 | Real s = 1; | |
410 | auto f3 = [&](Real t)->Real { | |
411 | Real x = s * t; | |
412 | if (x > boost::math::tools::log_max_value<Real>()) | |
413 | { | |
414 | return (Real) 0; | |
415 | } | |
416 | return cos(a * t) * exp(-x); | |
417 | }; | |
418 | ||
b32b8144 FG |
419 | // Since the integrand is oscillatory, we increase the tolerance: |
420 | Real tol_mult = 10; | |
421 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
1e59de90 | 422 | if (!std::is_class<Real>::value) |
b32b8144 | 423 | { |
92f5a8d4 TL |
424 | // For high oscillation frequency, the quadrature sum is ill-conditioned. |
425 | Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1); | |
426 | Q_expected = s/(a*a+s*s); | |
b32b8144 FG |
427 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10) |
428 | tol_mult = 5000; // we should really investigate this more?? | |
429 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult*tol); | |
430 | } | |
431 | ||
432 | // | |
433 | // This one doesn't pass for real_concept.. | |
434 | // | |
435 | if (std::numeric_limits<Real>::is_specialized) | |
436 | { | |
437 | // Laplace transform of J_0(t): | |
438 | auto f4 = [&](Real t)->Real { | |
439 | Real x = s * t; | |
440 | if (x > boost::math::tools::log_max_value<Real>()) | |
441 | { | |
442 | return (Real)0; | |
443 | } | |
444 | return boost::math::cyl_bessel_j(0, t) * exp(-x); | |
445 | }; | |
446 | ||
447 | Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1); | |
448 | Q_expected = 1 / sqrt(1 + s*s); | |
449 | tol_mult = 3; | |
450 | // Multiprecision type have higher error rates, probably evaluation of f() is less accurate: | |
451 | if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10) | |
452 | tol_mult = 750; | |
453 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult * tol); | |
454 | } | |
92f5a8d4 | 455 | auto f6 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t)*log(t));}; |
b32b8144 FG |
456 | Q = integrator.integrate(f6, get_convergence_tolerance<Real>(), &error, &L1); |
457 | Q_expected = -boost::math::constants::root_pi<Real>()*(boost::math::constants::euler<Real>() + 2*ln_two<Real>())/4; | |
458 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
459 | ||
460 | // CRC Section 5.5, integral 591 | |
461 | // The parameter p allows us to control the strength of the singularity. | |
462 | // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk. | |
463 | // This converges only when our test type has an extended exponent range as all the area of the integral | |
464 | // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it. | |
465 | // "There's a lot of room at the bottom". | |
466 | // This version is transformed via argument substitution (exp(-x) for x) so that the integral is spread | |
467 | // over (0, INF). | |
468 | tol *= boost::math::tools::digits<Real>() > 100 ? 100000 : 75; | |
469 | for (Real pn = 99; pn > 0; pn -= 10) { | |
470 | Real p = pn / 100; | |
471 | auto f = [&](Real x)->Real | |
472 | { | |
473 | return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)))); | |
474 | }; | |
475 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
476 | Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>()); | |
477 | /* | |
478 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << p << std::endl; | |
479 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q << std::endl; | |
480 | std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q_expected << std::endl; | |
481 | std::cout << fabs((Q - Q_expected) / Q_expected) << std::endl; | |
482 | */ | |
483 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
484 | } | |
485 | // and for p < 1: | |
486 | for (Real p = -0.99; p < 0; p += 0.1) { | |
487 | auto f = [&](Real x)->Real | |
488 | { | |
489 | return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x)); | |
490 | }; | |
491 | Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
492 | Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>()); | |
493 | BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); | |
494 | } | |
495 | } | |
496 | ||
92f5a8d4 TL |
497 | template<class Complex> |
498 | void test_complex_modified_bessel() | |
499 | { | |
500 | std::cout << "Testing complex modified Bessel function on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n"; | |
501 | typedef typename Complex::value_type Real; | |
502 | Real tol = 100 * boost::math::tools::epsilon<Real>(); | |
503 | Real error; | |
504 | Real L1; | |
505 | auto integrator = get_integrator<Real>(); | |
506 | ||
507 | // Integral Representation of Modified Complex Bessel function: | |
508 | // https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions | |
509 | Complex z{2, 3}; | |
510 | const auto f = [&z](const Real& t)->Complex | |
511 | { | |
512 | using std::cosh; | |
513 | using std::exp; | |
514 | Real cosht = cosh(t); | |
515 | if (cosht > boost::math::tools::log_max_value<Real>()) | |
516 | { | |
517 | return Complex{0, 0}; | |
518 | } | |
519 | Complex arg = -z*cosht; | |
520 | Complex res = exp(arg); | |
521 | return res; | |
522 | }; | |
523 | ||
524 | Complex K0 = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1); | |
525 | ||
526 | // Mathematica code: N[BesselK[0, 2 + 3 I], 140] | |
1e59de90 TL |
527 | #ifdef BOOST_MATH_STANDALONE |
528 | BOOST_IF_CONSTEXPR (std::is_fundamental<Complex>::value) | |
529 | #endif | |
530 | { | |
531 | Real K0_x_expected = BOOST_MATH_TEST_VALUE(Real, -0.08296852656762551490517953520589186885781541203818846830385526187936132191822538822296497597191327722262903004145527496422090506197776994); | |
532 | Real K0_y_expected = BOOST_MATH_TEST_VALUE(Real, 0.027949603635183423629723306332336002340909030265538548521150904238352846705644065168365102147901993976999717171115546662967229050834575193041); | |
533 | BOOST_CHECK_CLOSE_FRACTION(K0.real(), K0_x_expected, tol); | |
534 | BOOST_CHECK_CLOSE_FRACTION(K0.imag(), K0_y_expected, tol); | |
535 | } | |
92f5a8d4 TL |
536 | } |
537 | ||
f67539c2 TL |
538 | template<typename Complex> |
539 | void test_complex_exponential_integral_E1(){ | |
540 | std::cout << "Testing complex exponential integral on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n"; | |
541 | typedef typename Complex::value_type Real; | |
542 | Real tol = 100 * boost::math::tools::epsilon<Real>(); | |
543 | Real error; | |
544 | Real L1; | |
545 | auto integrator = get_integrator<Real>(); | |
546 | ||
547 | Complex z{1.5,0.5}; | |
548 | ||
549 | // Integral representation of exponential integral E1, valid for Re z > 0 | |
550 | // https://en.wikipedia.org/wiki/Exponential_integral#Definitions | |
551 | auto f = [&z](const Real& t)->Complex | |
552 | { | |
553 | using std::exp; | |
554 | return exp(-z*t)/t; | |
555 | }; | |
556 | ||
557 | Real inf = std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(); | |
558 | ||
559 | Complex E1 = integrator.integrate(f,1,inf,get_convergence_tolerance<Real>(),&error,&L1); | |
560 | ||
561 | // Mathematica code: N[ExpIntegral[1,1.5 + 0.5 I],140] | |
1e59de90 TL |
562 | #ifdef BOOST_MATH_STANDALONE |
563 | BOOST_IF_CONSTEXPR (std::is_fundamental<Complex>::value) | |
564 | #endif | |
565 | { | |
566 | Real E1_real_expected = BOOST_MATH_TEST_VALUE(Real, 0.071702995463938694845949672113596046091766639758473558841839765788732549949008866887694451956003503764943496943262401868244277788066634858393); | |
567 | Real E1_imag_expected = BOOST_MATH_TEST_VALUE(Real, -0.065138628279238400564373880665751377423524428792583839078600260273866805818117625959446311737353882269129094759883720722150048944193926087208); | |
568 | BOOST_CHECK_CLOSE_FRACTION(E1.real(), E1_real_expected, tol); | |
569 | BOOST_CHECK_CLOSE_FRACTION(E1.imag(), E1_imag_expected, tol); | |
570 | } | |
f67539c2 TL |
571 | } |
572 | ||
b32b8144 FG |
573 | |
574 | BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test) | |
575 | { | |
576 | // | |
577 | // Uncomment to generate the coefficients: | |
578 | // | |
579 | ||
580 | /* | |
581 | std::cout << std::scientific << std::setprecision(8); | |
582 | print_levels(generate_constants<cpp_bin_float_100, float>(8), "f"); | |
583 | std::cout << std::setprecision(18); | |
584 | print_levels(generate_constants<cpp_bin_float_100, double>(8), ""); | |
585 | std::cout << std::setprecision(35); | |
586 | print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L"); | |
587 | */ | |
588 | ||
b32b8144 FG |
589 | #ifdef TEST1 |
590 | test_left_limit_infinite<float>(); | |
591 | test_right_limit_infinite<float>(); | |
592 | test_nr_examples<float>(); | |
593 | test_crc<float>(); | |
594 | #endif | |
595 | #ifdef TEST2 | |
596 | test_left_limit_infinite<double>(); | |
597 | test_right_limit_infinite<double>(); | |
598 | test_nr_examples<double>(); | |
599 | test_crc<double>(); | |
600 | #endif | |
11fdf7f2 | 601 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
b32b8144 | 602 | #ifdef TEST3 |
1e59de90 | 603 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
b32b8144 FG |
604 | test_left_limit_infinite<long double>(); |
605 | test_right_limit_infinite<long double>(); | |
606 | test_nr_examples<long double>(); | |
607 | test_crc<long double>(); | |
608 | #endif | |
11fdf7f2 | 609 | #endif |
1e59de90 TL |
610 | #endif |
611 | #if defined(TEST4) && !defined(BOOST_MATH_NO_MP_TESTS) | |
b32b8144 FG |
612 | test_left_limit_infinite<cpp_bin_float_quad>(); |
613 | test_right_limit_infinite<cpp_bin_float_quad>(); | |
614 | test_nr_examples<cpp_bin_float_quad>(); | |
615 | test_crc<cpp_bin_float_quad>(); | |
616 | #endif | |
617 | ||
1e59de90 | 618 | #if !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) |
b32b8144 FG |
619 | #ifdef TEST5 |
620 | test_left_limit_infinite<boost::math::concepts::real_concept>(); | |
621 | test_right_limit_infinite<boost::math::concepts::real_concept>(); | |
622 | test_nr_examples<boost::math::concepts::real_concept>(); | |
623 | test_crc<boost::math::concepts::real_concept>(); | |
624 | #endif | |
11fdf7f2 | 625 | #endif |
1e59de90 | 626 | #if defined(TEST6) && !defined(BOOST_MATH_NO_MP_TESTS) |
b32b8144 FG |
627 | test_left_limit_infinite<boost::multiprecision::cpp_bin_float_50>(); |
628 | test_right_limit_infinite<boost::multiprecision::cpp_bin_float_50>(); | |
629 | test_nr_examples<boost::multiprecision::cpp_bin_float_50>(); | |
630 | test_crc<boost::multiprecision::cpp_bin_float_50>(); | |
631 | #endif | |
1e59de90 | 632 | #if defined(TEST7) && !defined(BOOST_MATH_NO_MP_TESTS) |
b32b8144 FG |
633 | test_left_limit_infinite<boost::multiprecision::cpp_dec_float_50>(); |
634 | test_right_limit_infinite<boost::multiprecision::cpp_dec_float_50>(); | |
635 | test_nr_examples<boost::multiprecision::cpp_dec_float_50>(); | |
92f5a8d4 TL |
636 | // |
637 | // This one causes stack overflows on the CI machine, but not locally, | |
f67539c2 | 638 | // assume it's due to restricted resources on the server, and <shrug> for now... |
92f5a8d4 | 639 | // |
1e59de90 | 640 | #if ! BOOST_WORKAROUND(BOOST_MSVC, == 1900) && !defined(BOOST_MATH_NO_MP_TESTS) |
b32b8144 FG |
641 | test_crc<boost::multiprecision::cpp_dec_float_50>(); |
642 | #endif | |
92f5a8d4 TL |
643 | #endif |
644 | #ifdef TEST8 | |
645 | test_complex_modified_bessel<std::complex<float>>(); | |
646 | test_complex_modified_bessel<std::complex<double>>(); | |
1e59de90 | 647 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
92f5a8d4 | 648 | test_complex_modified_bessel<std::complex<long double>>(); |
1e59de90 TL |
649 | #endif |
650 | #ifndef BOOST_MATH_NO_MP_TESTS | |
92f5a8d4 TL |
651 | test_complex_modified_bessel<boost::multiprecision::cpp_complex_quad>(); |
652 | #endif | |
1e59de90 | 653 | #endif |
f67539c2 TL |
654 | #ifdef TEST9 |
655 | test_complex_exponential_integral_E1<std::complex<float>>(); | |
656 | test_complex_exponential_integral_E1<std::complex<double>>(); | |
1e59de90 | 657 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
f67539c2 | 658 | test_complex_exponential_integral_E1<std::complex<long double>>(); |
1e59de90 TL |
659 | #endif |
660 | #ifndef BOOST_MATH_NO_MP_TESTS | |
f67539c2 TL |
661 | test_complex_exponential_integral_E1<boost::multiprecision::cpp_complex_quad>(); |
662 | #endif | |
1e59de90 TL |
663 | #endif |
664 | #ifdef TEST10 | |
665 | #if defined(BOOST_HAS_FLOAT128) && !defined(BOOST_MATH_NO_MP_TESTS) | |
666 | test_complex_modified_bessel<boost::multiprecision::complex128>(); | |
667 | test_complex_exponential_integral_E1<boost::multiprecision::complex128>(); | |
668 | #endif | |
669 | #endif | |
b32b8144 | 670 | } |