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)
7 #define BOOST_TEST_MODULE Gauss Kronrod_quadrature_test
10 #include <boost/config.hpp>
11 #include <boost/detail/workaround.hpp>
13 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
15 #include <boost/math/concepts/real_concept.hpp>
16 #include <boost/test/included/unit_test.hpp>
17 #include <boost/test/tools/floating_point_comparison.hpp>
18 #include <boost/math/quadrature/gauss_kronrod.hpp>
19 #include <boost/math/special_functions/sinc.hpp>
20 #include <boost/multiprecision/cpp_bin_float.hpp>
21 #include <boost/multiprecision/cpp_dec_float.hpp>
22 #include <boost/multiprecision/debug_adaptor.hpp>
24 #ifdef BOOST_HAS_FLOAT128
25 #include <boost/multiprecision/complex128.hpp>
28 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
36 #pragma warning(disable:4127) // Conditional expression is constant
57 using boost::math::quadrature::gauss_kronrod;
58 using boost::math::constants::pi;
59 using boost::math::constants::half_pi;
60 using boost::math::constants::two_div_pi;
61 using boost::math::constants::two_pi;
62 using boost::math::constants::half;
63 using boost::math::constants::third;
64 using boost::math::constants::half;
65 using boost::math::constants::third;
66 using boost::math::constants::catalan;
67 using boost::math::constants::ln_two;
68 using boost::math::constants::root_two;
69 using boost::math::constants::root_two_pi;
70 using boost::math::constants::root_pi;
71 using boost::multiprecision::cpp_bin_float_quad;
72 using boost::multiprecision::cpp_dec_float_50;
73 using boost::multiprecision::debug_adaptor;
74 using boost::multiprecision::number;
77 // Error rates depend only on the number of points in the approximation, not the type being tested,
78 // define all our expected errors here:
85 test_three_quad_error_id,
86 test_three_quad_error_id_2,
87 test_integration_over_real_line_error_id,
88 test_right_limit_infinite_error_id,
89 test_left_limit_infinite_error_id
92 template <unsigned Points>
93 double expected_error(unsigned)
95 return 0; // placeholder, all tests will fail
99 double expected_error<15>(unsigned id)
103 case test_ca_error_id:
105 case test_ca_error_id_2:
107 case test_three_quad_error_id:
109 case test_three_quad_error_id_2:
111 case test_integration_over_real_line_error_id:
113 case test_right_limit_infinite_error_id:
114 case test_left_limit_infinite_error_id:
117 return 0; // placeholder, all tests will fail
121 double expected_error<17>(unsigned id)
125 case test_ca_error_id:
127 case test_ca_error_id_2:
129 case test_three_quad_error_id:
131 case test_three_quad_error_id_2:
133 case test_integration_over_real_line_error_id:
135 case test_right_limit_infinite_error_id:
136 case test_left_limit_infinite_error_id:
139 return 0; // placeholder, all tests will fail
143 double expected_error<21>(unsigned id)
147 case test_ca_error_id:
149 case test_ca_error_id_2:
151 case test_three_quad_error_id:
153 case test_three_quad_error_id_2:
155 case test_integration_over_real_line_error_id:
156 return 6e-3; // doesn't get any better with more points!
157 case test_right_limit_infinite_error_id:
158 case test_left_limit_infinite_error_id:
161 return 0; // placeholder, all tests will fail
165 double expected_error<31>(unsigned id)
169 case test_ca_error_id:
171 case test_ca_error_id_2:
173 case test_three_quad_error_id:
175 case test_three_quad_error_id_2:
177 case test_integration_over_real_line_error_id:
178 return 6e-3; // doesn't get any better with more points!
179 case test_right_limit_infinite_error_id:
180 case test_left_limit_infinite_error_id:
183 return 0; // placeholder, all tests will fail
187 double expected_error<41>(unsigned id)
191 case test_ca_error_id:
193 case test_ca_error_id_2:
195 case test_three_quad_error_id:
197 case test_three_quad_error_id_2:
199 case test_integration_over_real_line_error_id:
200 return 5e-5; // doesn't get any better with more points!
201 case test_right_limit_infinite_error_id:
202 case test_left_limit_infinite_error_id:
205 return 0; // placeholder, all tests will fail
209 double expected_error<51>(unsigned id)
213 case test_ca_error_id:
215 case test_ca_error_id_2:
217 case test_three_quad_error_id:
219 case test_three_quad_error_id_2:
221 case test_integration_over_real_line_error_id:
223 case test_right_limit_infinite_error_id:
224 case test_left_limit_infinite_error_id:
227 return 0; // placeholder, all tests will fail
231 double expected_error<61>(unsigned id)
235 case test_ca_error_id:
237 case test_ca_error_id_2:
239 case test_three_quad_error_id:
241 case test_three_quad_error_id_2:
243 case test_integration_over_real_line_error_id:
245 case test_right_limit_infinite_error_id:
246 case test_left_limit_infinite_error_id:
249 return 0; // placeholder, all tests will fail
253 template<class Real, unsigned Points>
256 std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
257 Real tol = boost::math::tools::epsilon<Real>() * 10;
259 auto f = [](const Real& x)->Real
264 Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 0, 0, &error, &L1);
265 BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
266 BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
269 template<class Real, unsigned Points>
270 void test_quadratic()
272 std::cout << "Testing quadratic functions are integrated properly by Gauss Kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
273 Real tol = boost::math::tools::epsilon<Real>() * 10;
276 auto f = [](const Real& x)->Real { return 5*x*x + 7*x + 12; };
278 Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 0, 0, &error, &L1);
279 BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
280 BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
283 // Examples taken from
284 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
285 template<class Real, unsigned Points>
288 std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
289 Real tol = expected_error<Points>(test_ca_error_id);
293 auto f1 = [](const Real& x)->Real { return atan(x)/(x*(x*x + 1)) ; };
294 Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 0, 0, &error, &L1);
295 Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
296 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
297 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
299 auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
300 Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 0, 0, &error, &L1);
301 Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
302 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
303 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
305 tol = expected_error<Points>(test_ca_error_id_2);
306 auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
307 Q = gauss_kronrod<Real, Points>::integrate(f5, 0, 1, 0);
308 Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
309 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
312 template<class Real, unsigned Points>
313 void test_three_quadrature_schemes_examples()
315 std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
316 Real tol = expected_error<Points>(test_three_quad_error_id);
321 auto f1 = [](const Real& t)->Real { return t*boost::math::log1p(t); };
322 Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1, 0);
323 Q_expected = half<Real>()*half<Real>();
324 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
328 auto f2 = [](const Real& t)->Real { return t*t*atan(t); };
329 Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 0);
330 Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
331 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
334 auto f3 = [](const Real& t)->Real { return exp(t)*cos(t); };
335 Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>(), 0);
336 Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
337 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
340 auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
341 Q = gauss_kronrod<Real, Points>::integrate(f4, 0 , 1, 0);
342 Q_expected = 5*pi<Real>()*pi<Real>()/96;
343 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
345 tol = expected_error<Points>(test_three_quad_error_id_2);
347 auto f5 = [](const Real& t)->Real { return sqrt(t)*log(t); };
348 Q = gauss_kronrod<Real, Points>::integrate(f5, 0 , 1, 0);
349 Q_expected = -4/ (Real) 9;
350 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
353 auto f6 = [](const Real& t)->Real { return sqrt(1 - t*t); };
354 Q = gauss_kronrod<Real, Points>::integrate(f6, 0 , 1, 0);
355 Q_expected = pi<Real>()/4;
356 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
360 template<class Real, unsigned Points>
361 void test_integration_over_real_line()
363 std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
364 Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
370 auto f1 = [](const Real& t)->Real { return 1/(1+t*t);};
371 Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
372 Q_expected = pi<Real>();
373 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
374 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
377 template<class Real, unsigned Points>
378 void test_right_limit_infinite()
380 std::cout << "Testing right limit infinite for Gauss Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
381 Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
388 auto f1 = [](const Real& t)->Real { return 1/(1+t*t);};
389 Q = gauss_kronrod<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
390 Q_expected = half_pi<Real>();
391 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
393 auto f4 = [](const Real& t)->Real { return 1/(1+t*t); };
394 Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 0, 0, &error, &L1);
395 Q_expected = pi<Real>()/4;
396 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
399 template<class Real, unsigned Points>
400 void test_left_limit_infinite()
402 std::cout << "Testing left limit infinite for Gauss Kronrod in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
403 Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
408 auto f1 = [](const Real& t)->Real { return 1/(1+t*t);};
409 Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0), 0);
410 Q_expected = half_pi<Real>();
411 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
414 template<class Complex>
415 void test_complex_lambert_w()
417 std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
418 typedef typename Complex::value_type Real;
420 using boost::math::constants::pi;
422 auto lw = [&z](Real v)->Complex {
429 Real cotv = cosv/sinv;
431 Real t = (1-v*cotv)*(1-v*cotv) + v*v;
432 Real x = v*cscv*exp(-v*cotv);
434 Complex num = t*(z/pi<Real>());
435 Complex res = num/den;
439 //N[ProductLog[2+3*I], 150]
440 boost::math::quadrature::gauss_kronrod<Real, 61> integrator;
441 Complex Q = integrator.integrate(lw, (Real) 0, pi<Real>());
442 BOOST_CHECK_CLOSE_FRACTION(Q.real(), boost::lexical_cast<Real>("1.09007653448579084630177782678166964987102108635357778056449870727913321296238687023915522935120701763447787503167111962008709116746523970476893277703"), tol);
443 BOOST_CHECK_CLOSE_FRACTION(Q.imag(), boost::lexical_cast<Real>("0.530139720774838801426860213574121741928705631382703178297940568794784362495390544411799468140433404536019992695815009036975117285537382995180319280835"), tol);
446 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
449 std::cout << "Testing 15 point approximation:\n";
450 test_linear<double, 15>();
451 test_quadratic<double, 15>();
452 test_ca<double, 15>();
453 test_three_quadrature_schemes_examples<double, 15>();
454 test_integration_over_real_line<double, 15>();
455 test_right_limit_infinite<double, 15>();
456 test_left_limit_infinite<double, 15>();
458 // test one case where we do not have pre-computed constants:
459 std::cout << "Testing 17 point approximation:\n";
460 test_linear<double, 17>();
461 test_quadratic<double, 17>();
462 test_ca<double, 17>();
463 test_three_quadrature_schemes_examples<double, 17>();
464 test_integration_over_real_line<double, 17>();
465 test_right_limit_infinite<double, 17>();
466 test_left_limit_infinite<double, 17>();
467 test_complex_lambert_w<std::complex<double>>();
468 test_complex_lambert_w<std::complex<long double>>();
471 std::cout << "Testing 21 point approximation:\n";
472 test_linear<cpp_bin_float_quad, 21>();
473 test_quadratic<cpp_bin_float_quad, 21>();
474 test_ca<cpp_bin_float_quad, 21>();
475 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
476 test_integration_over_real_line<cpp_bin_float_quad, 21>();
477 test_right_limit_infinite<cpp_bin_float_quad, 21>();
478 test_left_limit_infinite<cpp_bin_float_quad, 21>();
480 std::cout << "Testing 31 point approximation:\n";
481 test_linear<cpp_bin_float_quad, 31>();
482 test_quadratic<cpp_bin_float_quad, 31>();
483 test_ca<cpp_bin_float_quad, 31>();
484 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
485 test_integration_over_real_line<cpp_bin_float_quad, 31>();
486 test_right_limit_infinite<cpp_bin_float_quad, 31>();
487 test_left_limit_infinite<cpp_bin_float_quad, 31>();
490 std::cout << "Testing 41 point approximation:\n";
491 test_linear<cpp_bin_float_quad, 41>();
492 test_quadratic<cpp_bin_float_quad, 41>();
493 test_ca<cpp_bin_float_quad, 41>();
494 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
495 test_integration_over_real_line<cpp_bin_float_quad, 41>();
496 test_right_limit_infinite<cpp_bin_float_quad, 41>();
497 test_left_limit_infinite<cpp_bin_float_quad, 41>();
499 std::cout << "Testing 51 point approximation:\n";
500 test_linear<cpp_bin_float_quad, 51>();
501 test_quadratic<cpp_bin_float_quad, 51>();
502 test_ca<cpp_bin_float_quad, 51>();
503 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 51>();
504 test_integration_over_real_line<cpp_bin_float_quad, 51>();
505 test_right_limit_infinite<cpp_bin_float_quad, 51>();
506 test_left_limit_infinite<cpp_bin_float_quad, 51>();
509 // Need at least one set of tests with expression templates turned on:
510 std::cout << "Testing 61 point approximation:\n";
511 test_linear<cpp_dec_float_50, 61>();
512 test_quadratic<cpp_dec_float_50, 61>();
513 test_ca<cpp_dec_float_50, 61>();
514 test_three_quadrature_schemes_examples<cpp_dec_float_50, 61>();
515 test_integration_over_real_line<cpp_dec_float_50, 61>();
516 test_right_limit_infinite<cpp_dec_float_50, 61>();
517 test_left_limit_infinite<cpp_dec_float_50, 61>();
518 #ifdef BOOST_HAS_FLOAT128
519 test_complex_lambert_w<boost::multiprecision::complex128>();
526 int main() { return 0; }