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/math/tools/test_value.hpp>
17 #include <boost/test/included/unit_test.hpp>
18 #include <boost/test/tools/floating_point_comparison.hpp>
19 #include <boost/math/quadrature/gauss_kronrod.hpp>
20 #include <boost/math/special_functions/sinc.hpp>
21 #include <boost/multiprecision/cpp_bin_float.hpp>
22 #include <boost/multiprecision/cpp_dec_float.hpp>
23 #include <boost/multiprecision/debug_adaptor.hpp>
25 #ifdef BOOST_HAS_FLOAT128
26 #include <boost/multiprecision/complex128.hpp>
29 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
37 #pragma warning(disable:4127) // Conditional expression is constant
58 using boost::math::quadrature::gauss_kronrod
;
59 using boost::math::constants::pi
;
60 using boost::math::constants::half_pi
;
61 using boost::math::constants::two_div_pi
;
62 using boost::math::constants::two_pi
;
63 using boost::math::constants::half
;
64 using boost::math::constants::third
;
65 using boost::math::constants::half
;
66 using boost::math::constants::third
;
67 using boost::math::constants::catalan
;
68 using boost::math::constants::ln_two
;
69 using boost::math::constants::root_two
;
70 using boost::math::constants::root_two_pi
;
71 using boost::math::constants::root_pi
;
72 using boost::multiprecision::cpp_bin_float_quad
;
73 using boost::multiprecision::cpp_dec_float_50
;
74 using boost::multiprecision::debug_adaptor
;
75 using boost::multiprecision::number
;
78 // Error rates depend only on the number of points in the approximation, not the type being tested,
79 // define all our expected errors here:
86 test_three_quad_error_id
,
87 test_three_quad_error_id_2
,
88 test_integration_over_real_line_error_id
,
89 test_right_limit_infinite_error_id
,
90 test_left_limit_infinite_error_id
93 template <unsigned Points
>
94 double expected_error(unsigned)
96 return 0; // placeholder, all tests will fail
100 double expected_error
<15>(unsigned id
)
104 case test_ca_error_id
:
106 case test_ca_error_id_2
:
108 case test_three_quad_error_id
:
110 case test_three_quad_error_id_2
:
112 case test_integration_over_real_line_error_id
:
114 case test_right_limit_infinite_error_id
:
115 case test_left_limit_infinite_error_id
:
118 return 0; // placeholder, all tests will fail
122 double expected_error
<17>(unsigned id
)
126 case test_ca_error_id
:
128 case test_ca_error_id_2
:
130 case test_three_quad_error_id
:
132 case test_three_quad_error_id_2
:
134 case test_integration_over_real_line_error_id
:
136 case test_right_limit_infinite_error_id
:
137 case test_left_limit_infinite_error_id
:
140 return 0; // placeholder, all tests will fail
144 double expected_error
<21>(unsigned id
)
148 case test_ca_error_id
:
150 case test_ca_error_id_2
:
152 case test_three_quad_error_id
:
154 case test_three_quad_error_id_2
:
156 case test_integration_over_real_line_error_id
:
157 return 6e-3; // doesn't get any better with more points!
158 case test_right_limit_infinite_error_id
:
159 case test_left_limit_infinite_error_id
:
162 return 0; // placeholder, all tests will fail
166 double expected_error
<31>(unsigned id
)
170 case test_ca_error_id
:
172 case test_ca_error_id_2
:
174 case test_three_quad_error_id
:
176 case test_three_quad_error_id_2
:
178 case test_integration_over_real_line_error_id
:
179 return 6e-3; // doesn't get any better with more points!
180 case test_right_limit_infinite_error_id
:
181 case test_left_limit_infinite_error_id
:
184 return 0; // placeholder, all tests will fail
188 double expected_error
<41>(unsigned id
)
192 case test_ca_error_id
:
194 case test_ca_error_id_2
:
196 case test_three_quad_error_id
:
198 case test_three_quad_error_id_2
:
200 case test_integration_over_real_line_error_id
:
201 return 5e-5; // doesn't get any better with more points!
202 case test_right_limit_infinite_error_id
:
203 case test_left_limit_infinite_error_id
:
206 return 0; // placeholder, all tests will fail
210 double expected_error
<51>(unsigned id
)
214 case test_ca_error_id
:
216 case test_ca_error_id_2
:
218 case test_three_quad_error_id
:
220 case test_three_quad_error_id_2
:
222 case test_integration_over_real_line_error_id
:
224 case test_right_limit_infinite_error_id
:
225 case test_left_limit_infinite_error_id
:
228 return 0; // placeholder, all tests will fail
232 double expected_error
<61>(unsigned id
)
236 case test_ca_error_id
:
238 case test_ca_error_id_2
:
240 case test_three_quad_error_id
:
242 case test_three_quad_error_id_2
:
244 case test_integration_over_real_line_error_id
:
246 case test_right_limit_infinite_error_id
:
247 case test_left_limit_infinite_error_id
:
250 return 0; // placeholder, all tests will fail
254 template<class Real
, unsigned Points
>
257 std::cout
<< "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
258 Real tol
= boost::math::tools::epsilon
<Real
>() * 10;
260 auto f
= [](const Real
& x
)->Real
265 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, (Real
) 0, (Real
) 1, 0, 0, &error
, &L1
);
266 BOOST_CHECK_CLOSE_FRACTION(Q
, 9.5, tol
);
267 BOOST_CHECK_CLOSE_FRACTION(L1
, 9.5, tol
);
269 Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, (Real
) 1, (Real
) 0, 0, 0, &error
, &L1
);
270 BOOST_CHECK_CLOSE_FRACTION(Q
, -9.5, tol
);
271 BOOST_CHECK_CLOSE_FRACTION(L1
, 9.5, tol
);
273 Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, (Real
) 0, (Real
) 0, 0, 0, &error
, &L1
);
274 BOOST_CHECK_CLOSE(Q
, Real(0), tol
);
277 template<class Real
, unsigned Points
>
278 void test_quadratic()
280 std::cout
<< "Testing quadratic functions are integrated properly by Gauss Kronrod on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
281 Real tol
= boost::math::tools::epsilon
<Real
>() * 10;
284 auto f
= [](const Real
& x
)->Real
{ return 5*x
*x
+ 7*x
+ 12; };
286 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f
, 0, 1, 0, 0, &error
, &L1
);
287 BOOST_CHECK_CLOSE_FRACTION(Q
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
288 BOOST_CHECK_CLOSE_FRACTION(L1
, (Real
) 17 + half
<Real
>()*third
<Real
>(), tol
);
291 // Examples taken from
292 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
293 template<class Real
, unsigned Points
>
296 std::cout
<< "Testing integration of C(a) on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
297 Real tol
= expected_error
<Points
>(test_ca_error_id
);
301 auto f1
= [](const Real
& x
)->Real
{ return atan(x
)/(x
*(x
*x
+ 1)) ; };
302 Real Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0, 1, 0, 0, &error
, &L1
);
303 Real Q_expected
= pi
<Real
>()*ln_two
<Real
>()/8 + catalan
<Real
>()*half
<Real
>();
304 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
305 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
307 auto f2
= [](Real x
)->Real
{ Real t0
= x
*x
+ 1; Real t1
= sqrt(t0
); return atan(t1
)/(t0
*t1
); };
308 Q
= gauss_kronrod
<Real
, Points
>::integrate(f2
, 0 , 1, 0, 0, &error
, &L1
);
309 Q_expected
= pi
<Real
>()/4 - pi
<Real
>()/root_two
<Real
>() + 3*atan(root_two
<Real
>())/root_two
<Real
>();
310 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
311 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
313 tol
= expected_error
<Points
>(test_ca_error_id_2
);
314 auto f5
= [](Real t
)->Real
{ return t
*t
*log(t
)/((t
*t
- 1)*(t
*t
*t
*t
+ 1)); };
315 Q
= gauss_kronrod
<Real
, Points
>::integrate(f5
, 0, 1, 0);
316 Q_expected
= pi
<Real
>()*pi
<Real
>()*(2 - root_two
<Real
>())/32;
317 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
320 template<class Real
, unsigned Points
>
321 void test_three_quadrature_schemes_examples()
323 std::cout
<< "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
324 Real tol
= expected_error
<Points
>(test_three_quad_error_id
);
329 auto f1
= [](const Real
& t
)->Real
{ return t
*boost::math::log1p(t
); };
330 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0 , 1, 0);
331 Q_expected
= half
<Real
>()*half
<Real
>();
332 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
336 auto f2
= [](const Real
& t
)->Real
{ return t
*t
*atan(t
); };
337 Q
= gauss_kronrod
<Real
, Points
>::integrate(f2
, 0 , 1, 0);
338 Q_expected
= (pi
<Real
>() -2 + 2*ln_two
<Real
>())/12;
339 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, 2 * tol
);
342 auto f3
= [](const Real
& t
)->Real
{ return exp(t
)*cos(t
); };
343 Q
= gauss_kronrod
<Real
, Points
>::integrate(f3
, 0, half_pi
<Real
>(), 0);
344 Q_expected
= boost::math::expm1(half_pi
<Real
>())*half
<Real
>();
345 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
348 auto f4
= [](Real x
)->Real
{ Real t0
= sqrt(x
*x
+ 2); return atan(t0
)/(t0
*(x
*x
+1)); };
349 Q
= gauss_kronrod
<Real
, Points
>::integrate(f4
, 0 , 1, 0);
350 Q_expected
= 5*pi
<Real
>()*pi
<Real
>()/96;
351 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
353 tol
= expected_error
<Points
>(test_three_quad_error_id_2
);
355 auto f5
= [](const Real
& t
)->Real
{ return sqrt(t
)*log(t
); };
356 Q
= gauss_kronrod
<Real
, Points
>::integrate(f5
, 0 , 1, 0);
357 Q_expected
= -4/ (Real
) 9;
358 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
361 auto f6
= [](const Real
& t
)->Real
{ return sqrt(1 - t
*t
); };
362 Q
= gauss_kronrod
<Real
, Points
>::integrate(f6
, 0 , 1, 0);
363 Q_expected
= pi
<Real
>()/4;
364 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
368 template<class Real
, unsigned Points
>
369 void test_integration_over_real_line()
371 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";
372 Real tol
= expected_error
<Points
>(test_integration_over_real_line_error_id
);
378 auto f1
= [](const Real
& t
)->Real
{ return 1/(1+t
*t
);};
379 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, -boost::math::tools::max_value
<Real
>(), boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
380 Q_expected
= pi
<Real
>();
381 BOOST_CHECK_CLOSE_FRACTION(Q
, Q_expected
, tol
);
382 BOOST_CHECK_CLOSE_FRACTION(L1
, Q_expected
, tol
);
385 template<class Real
, unsigned Points
>
386 void test_right_limit_infinite()
388 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";
389 Real tol
= expected_error
<Points
>(test_right_limit_infinite_error_id
);
396 auto f1
= [](const Real
& t
)->Real
{ return 1/(1+t
*t
);};
397 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, 0, boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
398 Q_expected
= half_pi
<Real
>();
399 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
401 auto f4
= [](const Real
& t
)->Real
{ return 1/(1+t
*t
); };
402 Q
= gauss_kronrod
<Real
, Points
>::integrate(f4
, 1, boost::math::tools::max_value
<Real
>(), 0, 0, &error
, &L1
);
403 Q_expected
= pi
<Real
>()/4;
404 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
407 template<class Real
, unsigned Points
>
408 void test_left_limit_infinite()
410 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";
411 Real tol
= expected_error
<Points
>(test_left_limit_infinite_error_id
);
416 auto f1
= [](const Real
& t
)->Real
{ return 1/(1+t
*t
);};
417 Q
= gauss_kronrod
<Real
, Points
>::integrate(f1
, -boost::math::tools::max_value
<Real
>(), Real(0), 0);
418 Q_expected
= half_pi
<Real
>();
419 BOOST_CHECK_CLOSE(Q
, Q_expected
, 100*tol
);
422 template<class Complex
>
423 void test_complex_lambert_w()
425 std::cout
<< "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id
<Complex
>().pretty_name() << "\n";
426 typedef typename
Complex::value_type Real
;
428 using boost::math::constants::pi
;
430 auto lw
= [&z
](Real v
)->Complex
{
437 Real cotv
= cosv
/sinv
;
439 Real t
= (1-v
*cotv
)*(1-v
*cotv
) + v
*v
;
440 Real x
= v
*cscv
*exp(-v
*cotv
);
442 Complex num
= t
*(z
/pi
<Real
>());
443 Complex res
= num
/den
;
447 //N[ProductLog[2+3*I], 150]
448 boost::math::quadrature::gauss_kronrod
<Real
, 61> integrator
;
449 Complex Q
= integrator
.integrate(lw
, (Real
) 0, pi
<Real
>());
450 BOOST_CHECK_CLOSE_FRACTION(Q
.real(), BOOST_MATH_TEST_VALUE(Real
, 1.0900765344857908463017778267816696498710210863535777805644), tol
);
451 BOOST_CHECK_CLOSE_FRACTION(Q
.imag(), BOOST_MATH_TEST_VALUE(Real
, 0.5301397207748388014268602135741217419287056313827031782979), tol
);
454 BOOST_AUTO_TEST_CASE(gauss_quadrature_test
)
457 std::cout
<< "Testing 15 point approximation:\n";
458 test_linear
<double, 15>();
459 test_quadratic
<double, 15>();
460 test_ca
<double, 15>();
461 test_three_quadrature_schemes_examples
<double, 15>();
462 test_integration_over_real_line
<double, 15>();
463 test_right_limit_infinite
<double, 15>();
464 test_left_limit_infinite
<double, 15>();
466 // test one case where we do not have pre-computed constants:
467 std::cout
<< "Testing 17 point approximation:\n";
468 test_linear
<double, 17>();
469 test_quadratic
<double, 17>();
470 test_ca
<double, 17>();
471 test_three_quadrature_schemes_examples
<double, 17>();
472 test_integration_over_real_line
<double, 17>();
473 test_right_limit_infinite
<double, 17>();
474 test_left_limit_infinite
<double, 17>();
475 test_complex_lambert_w
<std::complex<double>>();
476 test_complex_lambert_w
<std::complex<long double>>();
479 std::cout
<< "Testing 21 point approximation:\n";
480 test_linear
<cpp_bin_float_quad
, 21>();
481 test_quadratic
<cpp_bin_float_quad
, 21>();
482 test_ca
<cpp_bin_float_quad
, 21>();
483 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 21>();
484 test_integration_over_real_line
<cpp_bin_float_quad
, 21>();
485 test_right_limit_infinite
<cpp_bin_float_quad
, 21>();
486 test_left_limit_infinite
<cpp_bin_float_quad
, 21>();
488 std::cout
<< "Testing 31 point approximation:\n";
489 test_linear
<cpp_bin_float_quad
, 31>();
490 test_quadratic
<cpp_bin_float_quad
, 31>();
491 test_ca
<cpp_bin_float_quad
, 31>();
492 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 31>();
493 test_integration_over_real_line
<cpp_bin_float_quad
, 31>();
494 test_right_limit_infinite
<cpp_bin_float_quad
, 31>();
495 test_left_limit_infinite
<cpp_bin_float_quad
, 31>();
498 std::cout
<< "Testing 41 point approximation:\n";
499 test_linear
<cpp_bin_float_quad
, 41>();
500 test_quadratic
<cpp_bin_float_quad
, 41>();
501 test_ca
<cpp_bin_float_quad
, 41>();
502 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 41>();
503 test_integration_over_real_line
<cpp_bin_float_quad
, 41>();
504 test_right_limit_infinite
<cpp_bin_float_quad
, 41>();
505 test_left_limit_infinite
<cpp_bin_float_quad
, 41>();
507 std::cout
<< "Testing 51 point approximation:\n";
508 test_linear
<cpp_bin_float_quad
, 51>();
509 test_quadratic
<cpp_bin_float_quad
, 51>();
510 test_ca
<cpp_bin_float_quad
, 51>();
511 test_three_quadrature_schemes_examples
<cpp_bin_float_quad
, 51>();
512 test_integration_over_real_line
<cpp_bin_float_quad
, 51>();
513 test_right_limit_infinite
<cpp_bin_float_quad
, 51>();
514 test_left_limit_infinite
<cpp_bin_float_quad
, 51>();
517 // Need at least one set of tests with expression templates turned on:
518 std::cout
<< "Testing 61 point approximation:\n";
519 test_linear
<cpp_dec_float_50
, 61>();
520 test_quadratic
<cpp_dec_float_50
, 61>();
521 test_ca
<cpp_dec_float_50
, 61>();
522 test_three_quadrature_schemes_examples
<cpp_dec_float_50
, 61>();
523 test_integration_over_real_line
<cpp_dec_float_50
, 61>();
524 test_right_limit_infinite
<cpp_dec_float_50
, 61>();
525 test_left_limit_infinite
<cpp_dec_float_50
, 61>();
526 #ifdef BOOST_HAS_FLOAT128
527 test_complex_lambert_w
<boost::multiprecision::complex128
>();
534 int main() { return 0; }