Imported Upstream version 1.72.0
[platform/upstream/boost.git] / libs / math / test / gauss_kronrod_quadrature_test.cpp
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 Gauss Kronrod_quadrature_test
8
9 #include <complex>
10 #include <boost/config.hpp>
11 #include <boost/detail/workaround.hpp>
12
13 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
14
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>
23
24 #ifdef BOOST_HAS_FLOAT128
25 #include <boost/multiprecision/complex128.hpp>
26 #endif
27
28 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
29 #  define TEST1
30 #  define TEST1A
31 #  define TEST2
32 #  define TEST3
33 #endif
34
35 #ifdef _MSC_VER
36 #pragma warning(disable:4127)  // Conditional expression is constant
37 #endif
38
39 using std::expm1;
40 using std::atan;
41 using std::tan;
42 using std::log;
43 using std::log1p;
44 using std::asinh;
45 using std::atanh;
46 using std::sqrt;
47 using std::isnormal;
48 using std::abs;
49 using std::sinh;
50 using std::tanh;
51 using std::cosh;
52 using std::pow;
53 using std::exp;
54 using std::sin;
55 using std::cos;
56 using std::string;
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;
75
76 //
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:
79 //
80
81 enum
82 {
83    test_ca_error_id,
84    test_ca_error_id_2,
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
90 };
91
92 template <unsigned Points>
93 double expected_error(unsigned)
94 {
95    return 0; // placeholder, all tests will fail
96 }
97
98 template <>
99 double expected_error<15>(unsigned id)
100 {
101    switch (id)
102    {
103    case test_ca_error_id:
104       return 1e-7;
105    case test_ca_error_id_2:
106       return 2e-5;
107    case test_three_quad_error_id:
108       return 1e-8;
109    case test_three_quad_error_id_2:
110       return 3.5e-3;
111    case test_integration_over_real_line_error_id:
112       return 6e-3;
113    case test_right_limit_infinite_error_id:
114    case test_left_limit_infinite_error_id:
115       return 1e-5;
116    }
117    return 0;  // placeholder, all tests will fail
118 }
119
120 template <>
121 double expected_error<17>(unsigned id)
122 {
123    switch (id)
124    {
125    case test_ca_error_id:
126       return 1e-7;
127    case test_ca_error_id_2:
128       return 2e-5;
129    case test_three_quad_error_id:
130       return 1e-8;
131    case test_three_quad_error_id_2:
132       return 3.5e-3;
133    case test_integration_over_real_line_error_id:
134       return 6e-3;
135    case test_right_limit_infinite_error_id:
136    case test_left_limit_infinite_error_id:
137       return 1e-5;
138    }
139    return 0;  // placeholder, all tests will fail
140 }
141
142 template <>
143 double expected_error<21>(unsigned id)
144 {
145    switch (id)
146    {
147    case test_ca_error_id:
148       return 1e-12;
149    case test_ca_error_id_2:
150       return 3e-6;
151    case test_three_quad_error_id:
152       return 2e-13;
153    case test_three_quad_error_id_2:
154       return 2e-3;
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:
159       return 5e-8;
160    }
161    return 0;  // placeholder, all tests will fail
162 }
163
164 template <>
165 double expected_error<31>(unsigned id)
166 {
167    switch (id)
168    {
169    case test_ca_error_id:
170       return 6e-20;
171    case test_ca_error_id_2:
172       return 3e-7;
173    case test_three_quad_error_id:
174       return 1e-19;
175    case test_three_quad_error_id_2:
176       return 6e-4;
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:
181       return 5e-11;
182    }
183    return 0;  // placeholder, all tests will fail
184 }
185
186 template <>
187 double expected_error<41>(unsigned id)
188 {
189    switch (id)
190    {
191    case test_ca_error_id:
192       return 1e-26;
193    case test_ca_error_id_2:
194       return 1e-7;
195    case test_three_quad_error_id:
196       return 3e-27;
197    case test_three_quad_error_id_2:
198       return 3e-4;
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:
203       return 1e-15;
204    }
205    return 0;  // placeholder, all tests will fail
206 }
207
208 template <>
209 double expected_error<51>(unsigned id)
210 {
211    switch (id)
212    {
213    case test_ca_error_id:
214       return 5e-33;
215    case test_ca_error_id_2:
216       return 1e-8;
217    case test_three_quad_error_id:
218       return 1e-32;
219    case test_three_quad_error_id_2:
220       return 3e-4;
221    case test_integration_over_real_line_error_id:
222       return 1e-14;
223    case test_right_limit_infinite_error_id:
224    case test_left_limit_infinite_error_id:
225       return 3e-19;
226    }
227    return 0;  // placeholder, all tests will fail
228 }
229
230 template <>
231 double expected_error<61>(unsigned id)
232 {
233    switch (id)
234    {
235    case test_ca_error_id:
236       return 5e-34;
237    case test_ca_error_id_2:
238       return 5e-9;
239    case test_three_quad_error_id:
240       return 4e-34;
241    case test_three_quad_error_id_2:
242       return 1e-4;
243    case test_integration_over_real_line_error_id:
244       return 1e-16;
245    case test_right_limit_infinite_error_id:
246    case test_left_limit_infinite_error_id:
247       return 3e-23;
248    }
249    return 0;  // placeholder, all tests will fail
250 }
251
252
253 template<class Real, unsigned Points>
254 void test_linear()
255 {
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;
258     Real error;
259     auto f = [](const Real& x)->Real
260     {
261        return 5*x + 7;
262     };
263     Real L1;
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);
267 }
268
269 template<class Real, unsigned Points>
270 void test_quadratic()
271 {
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;
274     Real error;
275
276     auto f = [](const Real& x)->Real { return 5*x*x + 7*x + 12; };
277     Real L1;
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);
281 }
282
283 // Examples taken from
284 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
285 template<class Real, unsigned Points>
286 void test_ca()
287 {
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);
290     Real L1;
291     Real error;
292
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);
298
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);
304
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);
310 }
311
312 template<class Real, unsigned Points>
313 void test_three_quadrature_schemes_examples()
314 {
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);
317     Real Q;
318     Real Q_expected;
319
320     // Example 1:
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);
325
326
327     // Example 2:
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);
332
333     // Example 3:
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);
338
339     // Example 4:
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);
344
345     tol = expected_error<Points>(test_three_quad_error_id_2);
346     // Example 5:
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);
351
352     // Example 6:
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);
357 }
358
359
360 template<class Real, unsigned Points>
361 void test_integration_over_real_line()
362 {
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);
365     Real Q;
366     Real Q_expected;
367     Real L1;
368     Real error;
369
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);
375 }
376
377 template<class Real, unsigned Points>
378 void test_right_limit_infinite()
379 {
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);
382     Real Q;
383     Real Q_expected;
384     Real L1;
385     Real error;
386
387     // Example 11:
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);
392
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);
397 }
398
399 template<class Real, unsigned Points>
400 void test_left_limit_infinite()
401 {
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);
404     Real Q;
405     Real Q_expected;
406
407     // Example 11:
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);
412 }
413
414 template<class Complex>
415 void test_complex_lambert_w()
416 {
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;
419     Real tol = 10e-9;
420     using boost::math::constants::pi;
421     Complex z{2, 3};
422     auto lw = [&z](Real v)->Complex {
423       using std::cos;
424       using std::sin;
425       using std::exp;
426       Real sinv = sin(v);
427       Real cosv = cos(v);
428
429       Real cotv = cosv/sinv;
430       Real cscv = 1/sinv;
431       Real t = (1-v*cotv)*(1-v*cotv) + v*v;
432       Real x = v*cscv*exp(-v*cotv);
433       Complex den = z + x;
434       Complex num = t*(z/pi<Real>());
435       Complex res = num/den;
436       return res;
437     };
438
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);
444 }
445
446 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
447 {
448 #ifdef TEST1
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>();
457
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>>();
469 #endif
470 #ifdef TEST1A
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>();
479
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>();
488 #endif
489 #ifdef TEST2
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>();
498
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>();
507 #endif
508 #ifdef TEST3
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>();
520 #endif
521 #endif
522 }
523
524 #else
525
526 int main() { return 0; }
527
528 #endif