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
9 #include <boost/config.hpp>
10 #include <boost/detail/workaround.hpp>
12 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
14 #include <boost/math/concepts/real_concept.hpp>
15 #include <boost/test/included/unit_test.hpp>
16 #include <boost/test/floating_point_comparison.hpp>
17 #include <boost/math/quadrature/gauss_kronrod.hpp>
18 #include <boost/math/special_functions/sinc.hpp>
19 #include <boost/multiprecision/cpp_bin_float.hpp>
20 #include <boost/multiprecision/cpp_dec_float.hpp>
21 #include <boost/multiprecision/debug_adaptor.hpp>
23 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
31 #pragma warning(disable:4127) // Conditional expression is constant
52 using boost::math::quadrature::gauss_kronrod
;
53 using boost::math::constants::pi
;
54 using boost::math::constants::half_pi
;
55 using boost::math::constants::two_div_pi
;
56 using boost::math::constants::two_pi
;
57 using boost::math::constants::half
;
58 using boost::math::constants::third
;
59 using boost::math::constants::half
;
60 using boost::math::constants::third
;
61 using boost::math::constants::catalan
;
62 using boost::math::constants::ln_two
;
63 using boost::math::constants::root_two
;
64 using boost::math::constants::root_two_pi
;
65 using boost::math::constants::root_pi
;
66 using boost::multiprecision::cpp_bin_float_quad
;
67 using boost::multiprecision::cpp_dec_float_50
;
68 using boost::multiprecision::debug_adaptor
;
69 using boost::multiprecision::number
;
72 // Error rates depend only on the number of points in the approximation, not the type being tested,
73 // define all our expected errors here:
80 test_three_quad_error_id
,
81 test_three_quad_error_id_2
,
82 test_integration_over_real_line_error_id
,
83 test_right_limit_infinite_error_id
,
84 test_left_limit_infinite_error_id
87 template <unsigned Points
>
88 double expected_error(unsigned)
90 return 0; // placeholder, all tests will fail
94 double expected_error
<15>(unsigned id
)
98 case test_ca_error_id
:
100 case test_ca_error_id_2
:
102 case test_three_quad_error_id
:
104 case test_three_quad_error_id_2
:
106 case test_integration_over_real_line_error_id
:
108 case test_right_limit_infinite_error_id
:
109 case test_left_limit_infinite_error_id
:
112 return 0; // placeholder, all tests will fail
116 double expected_error
<17>(unsigned id
)
120 case test_ca_error_id
:
122 case test_ca_error_id_2
:
124 case test_three_quad_error_id
:
126 case test_three_quad_error_id_2
:
128 case test_integration_over_real_line_error_id
:
130 case test_right_limit_infinite_error_id
:
131 case test_left_limit_infinite_error_id
:
134 return 0; // placeholder, all tests will fail
138 double expected_error
<21>(unsigned id
)
142 case test_ca_error_id
:
144 case test_ca_error_id_2
:
146 case test_three_quad_error_id
:
148 case test_three_quad_error_id_2
:
150 case test_integration_over_real_line_error_id
:
151 return 6e-3; // doesn't get any better with more points!
152 case test_right_limit_infinite_error_id
:
153 case test_left_limit_infinite_error_id
:
156 return 0; // placeholder, all tests will fail
160 double expected_error
<31>(unsigned id
)
164 case test_ca_error_id
:
166 case test_ca_error_id_2
:
168 case test_three_quad_error_id
:
170 case test_three_quad_error_id_2
:
172 case test_integration_over_real_line_error_id
:
173 return 6e-3; // doesn't get any better with more points!
174 case test_right_limit_infinite_error_id
:
175 case test_left_limit_infinite_error_id
:
178 return 0; // placeholder, all tests will fail
182 double expected_error
<41>(unsigned id
)
186 case test_ca_error_id
:
188 case test_ca_error_id_2
:
190 case test_three_quad_error_id
:
192 case test_three_quad_error_id_2
:
194 case test_integration_over_real_line_error_id
:
195 return 5e-5; // doesn't get any better with more points!
196 case test_right_limit_infinite_error_id
:
197 case test_left_limit_infinite_error_id
:
200 return 0; // placeholder, all tests will fail
204 double expected_error
<51>(unsigned id
)
208 case test_ca_error_id
:
210 case test_ca_error_id_2
:
212 case test_three_quad_error_id
:
214 case test_three_quad_error_id_2
:
216 case test_integration_over_real_line_error_id
:
218 case test_right_limit_infinite_error_id
:
219 case test_left_limit_infinite_error_id
:
222 return 0; // placeholder, all tests will fail
226 double expected_error
<61>(unsigned id
)
230 case test_ca_error_id
:
232 case test_ca_error_id_2
:
234 case test_three_quad_error_id
:
236 case test_three_quad_error_id_2
:
238 case test_integration_over_real_line_error_id
:
240 case test_right_limit_infinite_error_id
:
241 case test_left_limit_infinite_error_id
:
244 return 0; // placeholder, all tests will fail
248 template<class Real
, unsigned Points
>
251 std::cout
<< "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
252 Real tol
= boost::math::tools::epsilon
<Real
>() * 10;
254 auto f
= [](const Real
& x
)
259 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, (Real
) 0, (Real
) 1, 0, 0, &error
, &L1
);
260 BOOST_CHECK_CLOSE_FRACTION(Q
, 9.5, tol
);
261 BOOST_CHECK_CLOSE_FRACTION(L1
, 9.5, tol
);
264 template<class Real
, unsigned Points
>
265 void test_quadratic()
267 std::cout
<< "Testing quadratic functions are integrated properly by Gauss Kronrod on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
268 Real tol
= boost::math::tools::epsilon
<Real
>() * 10;
271 auto f
= [](const Real
& x
) { return 5*x
*x
+ 7*x
+ 12; };
273 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, 0, 1, 0, 0, &error
, &L1
);
274 BOOST_CHECK_CLOSE_FRACTION(Q
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
275 BOOST_CHECK_CLOSE_FRACTION(L1
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
278 // Examples taken from
279 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
280 template<class Real
, unsigned Points
>
283 std::cout
<< "Testing integration of C(a) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
284 Real tol
= expected_error
<Points
>(test_ca_error_id
);
288 auto f1
= [](const Real
& x
) { return atan(x
)/(x
*(x
*x
+ 1)) ; };
289 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0, 1, 0, 0, &error
, &L1
);
290 Real Q_expected
= pi
<Real
>()*ln_two
<Real
>()/8 + catalan
<Real
>()*half
<Real
>();
291 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
292 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
294 auto f2
= [](Real x
)->Real
{ Real t0
= x
*x
+ 1; Real t1
= sqrt(t0
); return atan(t1
)/(t0
*t1
); };
295 Q
= gauss_kronrod
<Real
, Points
>::integrate(f2
, 0 , 1, 0, 0, &error
, &L1
);
296 Q_expected
= pi
<Real
>()/4 - pi
<Real
>()/root_two
<Real
>() + 3*atan(root_two
<Real
>())/root_two
<Real
>();
297 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
298 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
300 tol
= expected_error
<Points
>(test_ca_error_id_2
);
301 auto f5
= [](Real t
)->Real
{ return t
*t
*log(t
)/((t
*t
- 1)*(t
*t
*t
*t
+ 1)); };
302 Q
= gauss_kronrod
<Real
, Points
>::integrate(f5
, 0, 1, 0);
303 Q_expected
= pi
<Real
>()*pi
<Real
>()*(2 - root_two
<Real
>())/32;
304 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
307 template<class Real
, unsigned Points
>
308 void test_three_quadrature_schemes_examples()
310 std::cout
<< "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
311 Real tol
= expected_error
<Points
>(test_three_quad_error_id
);
316 auto f1
= [](const Real
& t
)->Real
{ return t
*boost::math::log1p(t
); };
317 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0 , 1, 0);
318 Q_expected
= half
<Real
>()*half
<Real
>();
319 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
323 auto f2
= [](const Real
& t
)->Real
{ return t
*t
*atan(t
); };
324 Q
= gauss_kronrod
<Real
, Points
>::integrate(f2
, 0 , 1, 0);
325 Q_expected
= (pi
<Real
>() -2 + 2*ln_two
<Real
>())/12;
326 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 2 * tol
);
329 auto f3
= [](const Real
& t
)->Real
{ return exp(t
)*cos(t
); };
330 Q
= gauss_kronrod
<Real
, Points
>::integrate(f3
, 0, half_pi
<Real
>(), 0);
331 Q_expected
= boost::math::expm1(half_pi
<Real
>())*half
<Real
>();
332 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
335 auto f4
= [](Real x
)->Real
{ Real t0
= sqrt(x
*x
+ 2); return atan(t0
)/(t0
*(x
*x
+1)); };
336 Q
= gauss_kronrod
<Real
, Points
>::integrate(f4
, 0 , 1, 0);
337 Q_expected
= 5*pi
<Real
>()*pi
<Real
>()/96;
338 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
340 tol
= expected_error
<Points
>(test_three_quad_error_id_2
);
342 auto f5
= [](const Real
& t
)->Real
{ return sqrt(t
)*log(t
); };
343 Q
= gauss_kronrod
<Real
, Points
>::integrate(f5
, 0 , 1, 0);
344 Q_expected
= -4/ (Real
) 9;
345 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
348 auto f6
= [](const Real
& t
)->Real
{ return sqrt(1 - t
*t
); };
349 Q
= gauss_kronrod
<Real
, Points
>::integrate(f6
, 0 , 1, 0);
350 Q_expected
= pi
<Real
>()/4;
351 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
355 template<class Real
, unsigned Points
>
356 void test_integration_over_real_line()
358 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";
359 Real tol
= expected_error
<Points
>(test_integration_over_real_line_error_id
);
365 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
366 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, -boost::math::tools::max_value
<Real
>(), boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
367 Q_expected
= pi
<Real
>();
368 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
369 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
372 template<class Real
, unsigned Points
>
373 void test_right_limit_infinite()
375 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";
376 Real tol
= expected_error
<Points
>(test_right_limit_infinite_error_id
);
383 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
384 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0, boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
385 Q_expected
= half_pi
<Real
>();
386 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
388 auto f4
= [](const Real
& t
) { return 1/(1+t
*t
); };
389 Q
= gauss_kronrod
<Real
, Points
>::integrate(f4
, 1, boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
390 Q_expected
= pi
<Real
>()/4;
391 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
394 template<class Real
, unsigned Points
>
395 void test_left_limit_infinite()
397 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";
398 Real tol
= expected_error
<Points
>(test_left_limit_infinite_error_id
);
403 auto f1
= [](const Real
& t
) { return 1/(1+t
*t
);};
404 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, -boost::math::tools::max_value
<Real
>(), Real(0), 0);
405 Q_expected
= half_pi
<Real
>();
406 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
409 BOOST_AUTO_TEST_CASE(gauss_quadrature_test
)
412 std::cout
<< "Testing 15 point approximation:\n";
413 test_linear
<double, 15>();
414 test_quadratic
<double, 15>();
415 test_ca
<double, 15>();
416 test_three_quadrature_schemes_examples
<double, 15>();
417 test_integration_over_real_line
<double, 15>();
418 test_right_limit_infinite
<double, 15>();
419 test_left_limit_infinite
<double, 15>();
421 // test one case where we do not have pre-computed constants:
422 std::cout
<< "Testing 17 point approximation:\n";
423 test_linear
<double, 17>();
424 test_quadratic
<double, 17>();
425 test_ca
<double, 17>();
426 test_three_quadrature_schemes_examples
<double, 17>();
427 test_integration_over_real_line
<double, 17>();
428 test_right_limit_infinite
<double, 17>();
429 test_left_limit_infinite
<double, 17>();
432 std::cout
<< "Testing 21 point approximation:\n";
433 test_linear
<cpp_bin_float_quad
, 21>();
434 test_quadratic
<cpp_bin_float_quad
, 21>();
435 test_ca
<cpp_bin_float_quad
, 21>();
436 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 21>();
437 test_integration_over_real_line
<cpp_bin_float_quad
, 21>();
438 test_right_limit_infinite
<cpp_bin_float_quad
, 21>();
439 test_left_limit_infinite
<cpp_bin_float_quad
, 21>();
441 std::cout
<< "Testing 31 point approximation:\n";
442 test_linear
<cpp_bin_float_quad
, 31>();
443 test_quadratic
<cpp_bin_float_quad
, 31>();
444 test_ca
<cpp_bin_float_quad
, 31>();
445 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 31>();
446 test_integration_over_real_line
<cpp_bin_float_quad
, 31>();
447 test_right_limit_infinite
<cpp_bin_float_quad
, 31>();
448 test_left_limit_infinite
<cpp_bin_float_quad
, 31>();
451 std::cout
<< "Testing 41 point approximation:\n";
452 test_linear
<cpp_bin_float_quad
, 41>();
453 test_quadratic
<cpp_bin_float_quad
, 41>();
454 test_ca
<cpp_bin_float_quad
, 41>();
455 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 41>();
456 test_integration_over_real_line
<cpp_bin_float_quad
, 41>();
457 test_right_limit_infinite
<cpp_bin_float_quad
, 41>();
458 test_left_limit_infinite
<cpp_bin_float_quad
, 41>();
460 std::cout
<< "Testing 51 point approximation:\n";
461 test_linear
<cpp_bin_float_quad
, 51>();
462 test_quadratic
<cpp_bin_float_quad
, 51>();
463 test_ca
<cpp_bin_float_quad
, 51>();
464 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 51>();
465 test_integration_over_real_line
<cpp_bin_float_quad
, 51>();
466 test_right_limit_infinite
<cpp_bin_float_quad
, 51>();
467 test_left_limit_infinite
<cpp_bin_float_quad
, 51>();
470 // Need at least one set of tests with expression templates turned on:
471 std::cout
<< "Testing 61 point approximation:\n";
472 test_linear
<cpp_dec_float_50
, 61>();
473 test_quadratic
<cpp_dec_float_50
, 61>();
474 test_ca
<cpp_dec_float_50
, 61>();
475 test_three_quadrature_schemes_examples
<cpp_dec_float_50
, 61>();
476 test_integration_over_real_line
<cpp_dec_float_50
, 61>();
477 test_right_limit_infinite
<cpp_dec_float_50
, 61>();
478 test_left_limit_infinite
<cpp_dec_float_50
, 61>();
484 int main() { return 0; }