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