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 tanh_sinh_quadrature_test
8
9 #include <complex>
10 //#include <boost/multiprecision/mpc.hpp>
11 #include <boost/config.hpp>
12 #include <boost/detail/workaround.hpp>
13
14 #if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
15
16 #include <boost/math/concepts/real_concept.hpp>
17 #include <boost/test/included/unit_test.hpp>
18 #include <boost/test/tools/floating_point_comparison.hpp>
19 #include <boost/math/quadrature/gauss.hpp>
20 #include <boost/math/special_functions/sinc.hpp>
21 #include <boost/multiprecision/cpp_bin_float.hpp>
22 #include <boost/multiprecision/cpp_complex.hpp>
23
24 #ifdef BOOST_HAS_FLOAT128
25 #include <boost/multiprecision/complex128.hpp>
26 #endif
27
28 #ifdef _MSC_VER
29 #pragma warning(disable:4127) // Conditional expression is constant
30 #endif
31
32 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3)
33 # define TEST1
34 # define TEST2
35 # define TEST3
36 #endif
37
38 using std::expm1;
39 using std::atan;
40 using std::tan;
41 using std::log;
42 using std::log1p;
43 using std::asinh;
44 using std::atanh;
45 using std::sqrt;
46 using std::isnormal;
47 using std::abs;
48 using std::sinh;
49 using std::tanh;
50 using std::cosh;
51 using std::pow;
52 using std::exp;
53 using std::sin;
54 using std::cos;
55 using std::string;
56 using boost::math::quadrature::gauss;
57 using boost::math::constants::pi;
58 using boost::math::constants::half_pi;
59 using boost::math::constants::two_div_pi;
60 using boost::math::constants::two_pi;
61 using boost::math::constants::half;
62 using boost::math::constants::third;
63 using boost::math::constants::half;
64 using boost::math::constants::third;
65 using boost::math::constants::catalan;
66 using boost::math::constants::ln_two;
67 using boost::math::constants::root_two;
68 using boost::math::constants::root_two_pi;
69 using boost::math::constants::root_pi;
70 using boost::multiprecision::cpp_bin_float_quad;
71
72 //
73 // Error rates depend only on the number of points in the approximation, not the type being tested,
74 // define all our expected errors here:
75 //
76
77 enum
78 {
79 test_ca_error_id,
80 test_ca_error_id_2,
81 test_three_quad_error_id,
82 test_three_quad_error_id_2,
83 test_integration_over_real_line_error_id,
84 test_right_limit_infinite_error_id,
85 test_left_limit_infinite_error_id
86 };
87
88 template <unsigned Points>
expected_error(unsigned)89 double expected_error(unsigned)
90 {
91 return 0; // placeholder, all tests will fail
92 }
93
94 template <>
expected_error(unsigned id)95 double expected_error<7>(unsigned id)
96 {
97 switch (id)
98 {
99 case test_ca_error_id:
100 return 1e-7;
101 case test_ca_error_id_2:
102 return 2e-5;
103 case test_three_quad_error_id:
104 return 1e-8;
105 case test_three_quad_error_id_2:
106 return 3.5e-3;
107 case test_integration_over_real_line_error_id:
108 return 6e-3;
109 case test_right_limit_infinite_error_id:
110 case test_left_limit_infinite_error_id:
111 return 1e-5;
112 }
113 return 0; // placeholder, all tests will fail
114 }
115
116 template <>
expected_error(unsigned id)117 double expected_error<9>(unsigned id)
118 {
119 switch (id)
120 {
121 case test_ca_error_id:
122 return 1e-7;
123 case test_ca_error_id_2:
124 return 2e-5;
125 case test_three_quad_error_id:
126 return 1e-8;
127 case test_three_quad_error_id_2:
128 return 3.5e-3;
129 case test_integration_over_real_line_error_id:
130 return 6e-3;
131 case test_right_limit_infinite_error_id:
132 case test_left_limit_infinite_error_id:
133 return 1e-5;
134 }
135 return 0; // placeholder, all tests will fail
136 }
137
138 template <>
expected_error(unsigned id)139 double expected_error<10>(unsigned id)
140 {
141 switch (id)
142 {
143 case test_ca_error_id:
144 return 1e-12;
145 case test_ca_error_id_2:
146 return 3e-6;
147 case test_three_quad_error_id:
148 return 2e-13;
149 case test_three_quad_error_id_2:
150 return 2e-3;
151 case test_integration_over_real_line_error_id:
152 return 6e-3; // doesn't get any better with more points!
153 case test_right_limit_infinite_error_id:
154 case test_left_limit_infinite_error_id:
155 return 5e-8;
156 }
157 return 0; // placeholder, all tests will fail
158 }
159
160 template <>
expected_error(unsigned id)161 double expected_error<15>(unsigned id)
162 {
163 switch (id)
164 {
165 case test_ca_error_id:
166 return 6e-20;
167 case test_ca_error_id_2:
168 return 3e-7;
169 case test_three_quad_error_id:
170 return 1e-19;
171 case test_three_quad_error_id_2:
172 return 6e-4;
173 case test_integration_over_real_line_error_id:
174 return 6e-3; // doesn't get any better with more points!
175 case test_right_limit_infinite_error_id:
176 case test_left_limit_infinite_error_id:
177 return 5e-11;
178 }
179 return 0; // placeholder, all tests will fail
180 }
181
182 template <>
expected_error(unsigned id)183 double expected_error<20>(unsigned id)
184 {
185 switch (id)
186 {
187 case test_ca_error_id:
188 return 1e-26;
189 case test_ca_error_id_2:
190 return 1e-7;
191 case test_three_quad_error_id:
192 return 3e-27;
193 case test_three_quad_error_id_2:
194 return 3e-4;
195 case test_integration_over_real_line_error_id:
196 return 5e-5; // doesn't get any better with more points!
197 case test_right_limit_infinite_error_id:
198 case test_left_limit_infinite_error_id:
199 return 1e-15;
200 }
201 return 0; // placeholder, all tests will fail
202 }
203
204 template <>
expected_error(unsigned id)205 double expected_error<25>(unsigned id)
206 {
207 switch (id)
208 {
209 case test_ca_error_id:
210 return 5e-33;
211 case test_ca_error_id_2:
212 return 1e-8;
213 case test_three_quad_error_id:
214 return 1e-32;
215 case test_three_quad_error_id_2:
216 return 3e-4;
217 case test_integration_over_real_line_error_id:
218 return 1e-14;
219 case test_right_limit_infinite_error_id:
220 case test_left_limit_infinite_error_id:
221 return 3e-19;
222 }
223 return 0; // placeholder, all tests will fail
224 }
225
226 template <>
expected_error(unsigned id)227 double expected_error<30>(unsigned id)
228 {
229 switch (id)
230 {
231 case test_ca_error_id:
232 return 2e-34;
233 case test_ca_error_id_2:
234 return 5e-9;
235 case test_three_quad_error_id:
236 return 4e-34;
237 case test_three_quad_error_id_2:
238 return 1e-4;
239 case test_integration_over_real_line_error_id:
240 return 1e-16;
241 case test_right_limit_infinite_error_id:
242 case test_left_limit_infinite_error_id:
243 return 3e-23;
244 }
245 return 0; // placeholder, all tests will fail
246 }
247
248
249 template<class Real, unsigned Points>
test_linear()250 void test_linear()
251 {
252 std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
253 Real tol = boost::math::tools::epsilon<Real>() * 10;
254 auto f = [](const Real& x)
255 {
256 return 5*x + 7;
257 };
258 Real L1;
259 Real Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 1, &L1);
260 BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
261 BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
262 Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 0, &L1);
263 BOOST_CHECK_CLOSE(Q, 0, tol);
264 Q = gauss<Real, Points>::integrate(f, (Real) 1, (Real) 0, &L1);
265 BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
266 }
267
268 template<class Real, unsigned Points>
test_quadratic()269 void test_quadratic()
270 {
271 std::cout << "Testing quadratic functions are integrated properly by Gaussian quadrature on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
272 Real tol = boost::math::tools::epsilon<Real>() * 10;
273
274 auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
275 Real L1;
276 Real Q = gauss<Real, Points>::integrate(f, 0, 1, &L1);
277 BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
278 BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
279 }
280
281 // Examples taken from
282 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
283 template<class Real, unsigned Points>
test_ca()284 void test_ca()
285 {
286 std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
287 Real tol = expected_error<Points>(test_ca_error_id);
288 Real L1;
289
290 auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
291 Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
292 Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
293 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
294 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
295
296 auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
297 Q = gauss<Real, Points>::integrate(f2, 0 , 1, &L1);
298 Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
299 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
300 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
301
302 tol = expected_error<Points>(test_ca_error_id_2);
303 auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
304 Q = gauss<Real, Points>::integrate(f5, 0 , 1);
305 Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
306 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
307 }
308
309 template<class Real, unsigned Points>
test_three_quadrature_schemes_examples()310 void test_three_quadrature_schemes_examples()
311 {
312 std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
313 Real tol = expected_error<Points>(test_three_quad_error_id);
314 Real Q;
315 Real Q_expected;
316
317 // Example 1:
318 auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
319 Q = gauss<Real, Points>::integrate(f1, 0 , 1);
320 Q_expected = half<Real>()*half<Real>();
321 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
322
323
324 // Example 2:
325 auto f2 = [](const Real& t) { return t*t*atan(t); };
326 Q = gauss<Real, Points>::integrate(f2, 0 , 1);
327 Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
328 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
329
330 // Example 3:
331 auto f3 = [](const Real& t) { return exp(t)*cos(t); };
332 Q = gauss<Real, Points>::integrate(f3, 0, half_pi<Real>());
333 Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
334 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
335
336 // Example 4:
337 auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
338 Q = gauss<Real, Points>::integrate(f4, 0 , 1);
339 Q_expected = 5*pi<Real>()*pi<Real>()/96;
340 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
341
342 tol = expected_error<Points>(test_three_quad_error_id_2);
343 // Example 5:
344 auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
345 Q = gauss<Real, Points>::integrate(f5, 0 , 1);
346 Q_expected = -4/ (Real) 9;
347 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
348
349 // Example 6:
350 auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
351 Q = gauss<Real, Points>::integrate(f6, 0 , 1);
352 Q_expected = pi<Real>()/4;
353 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
354 }
355
356
357 template<class Real, unsigned Points>
test_integration_over_real_line()358 void test_integration_over_real_line()
359 {
360 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";
361 Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
362 Real Q;
363 Real Q_expected;
364 Real L1;
365
366 auto f1 = [](const Real& t) { return 1/(1+t*t);};
367 Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
368 Q_expected = pi<Real>();
369 BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
370 BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
371 }
372
373 template<class Real, unsigned Points>
test_right_limit_infinite()374 void test_right_limit_infinite()
375 {
376 std::cout << "Testing right limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
377 Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
378 Real Q;
379 Real Q_expected;
380 Real L1;
381
382 // Example 11:
383 auto f1 = [](const Real& t) { return 1/(1+t*t);};
384 Q = gauss<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), &L1);
385 Q_expected = half_pi<Real>();
386 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
387
388 auto f4 = [](const Real& t) { return 1/(1+t*t); };
389 Q = gauss<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), &L1);
390 Q_expected = pi<Real>()/4;
391 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
392 }
393
394 template<class Real, unsigned Points>
test_left_limit_infinite()395 void test_left_limit_infinite()
396 {
397 std::cout << "Testing left limit infinite for Gaussian quadrature 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);
399 Real Q;
400 Real Q_expected;
401
402 // Example 11:
403 auto f1 = [](const Real& t) { return 1/(1+t*t);};
404 Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0));
405 Q_expected = half_pi<Real>();
406 BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
407 }
408
409 template<class Complex>
test_complex_lambert_w()410 void test_complex_lambert_w()
411 {
412 std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
413 typedef typename Complex::value_type Real;
414 Real tol = 10e-9;
415 using boost::math::constants::pi;
416 Complex z{2, 3};
417 auto lw = [&z](Real v)->Complex {
418 using std::cos;
419 using std::sin;
420 using std::exp;
421 Real sinv = sin(v);
422 Real cosv = cos(v);
423
424 Real cotv = cosv/sinv;
425 Real cscv = 1/sinv;
426 Real t = (1-v*cotv)*(1-v*cotv) + v*v;
427 Real x = v*cscv*exp(-v*cotv);
428 Complex den = z + x;
429 Complex num = t*(z/pi<Real>());
430 Complex res = num/den;
431 return res;
432 };
433
434 //N[ProductLog[2+3*I], 150]
435 Complex Q = gauss<Real, 30>::integrate(lw, (Real) 0, pi<Real>());
436 BOOST_CHECK_CLOSE_FRACTION(Q.real(), boost::lexical_cast<Real>("1.09007653448579084630177782678166964987102108635357778056449870727913321296238687023915522935120701763447787503167111962008709116746523970476893277703"), tol);
437 BOOST_CHECK_CLOSE_FRACTION(Q.imag(), boost::lexical_cast<Real>("0.530139720774838801426860213574121741928705631382703178297940568794784362495390544411799468140433404536019992695815009036975117285537382995180319280835"), tol);
438 }
439
BOOST_AUTO_TEST_CASE(gauss_quadrature_test)440 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
441 {
442
443 #ifdef TEST1
444 test_linear<double, 7>();
445 test_quadratic<double, 7>();
446 test_ca<double, 7>();
447 test_three_quadrature_schemes_examples<double, 7>();
448 test_integration_over_real_line<double, 7>();
449 test_right_limit_infinite<double, 7>();
450 test_left_limit_infinite<double, 7>();
451
452 test_linear<double, 9>();
453 test_quadratic<double, 9>();
454 test_ca<double, 9>();
455 test_three_quadrature_schemes_examples<double, 9>();
456 test_integration_over_real_line<double, 9>();
457 test_right_limit_infinite<double, 9>();
458 test_left_limit_infinite<double, 9>();
459
460 test_linear<cpp_bin_float_quad, 10>();
461 test_quadratic<cpp_bin_float_quad, 10>();
462 test_ca<cpp_bin_float_quad, 10>();
463 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 10>();
464 test_integration_over_real_line<cpp_bin_float_quad, 10>();
465 test_right_limit_infinite<cpp_bin_float_quad, 10>();
466 test_left_limit_infinite<cpp_bin_float_quad, 10>();
467 #endif
468 #ifdef TEST2
469 test_linear<cpp_bin_float_quad, 15>();
470 test_quadratic<cpp_bin_float_quad, 15>();
471 test_ca<cpp_bin_float_quad, 15>();
472 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 15>();
473 test_integration_over_real_line<cpp_bin_float_quad, 15>();
474 test_right_limit_infinite<cpp_bin_float_quad, 15>();
475 test_left_limit_infinite<cpp_bin_float_quad, 15>();
476
477 test_linear<cpp_bin_float_quad, 20>();
478 test_quadratic<cpp_bin_float_quad, 20>();
479 test_ca<cpp_bin_float_quad, 20>();
480 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 20>();
481 test_integration_over_real_line<cpp_bin_float_quad, 20>();
482 test_right_limit_infinite<cpp_bin_float_quad, 20>();
483 test_left_limit_infinite<cpp_bin_float_quad, 20>();
484
485 test_linear<cpp_bin_float_quad, 25>();
486 test_quadratic<cpp_bin_float_quad, 25>();
487 test_ca<cpp_bin_float_quad, 25>();
488 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 25>();
489 test_integration_over_real_line<cpp_bin_float_quad, 25>();
490 test_right_limit_infinite<cpp_bin_float_quad, 25>();
491 test_left_limit_infinite<cpp_bin_float_quad, 25>();
492
493 test_linear<cpp_bin_float_quad, 30>();
494 test_quadratic<cpp_bin_float_quad, 30>();
495 test_ca<cpp_bin_float_quad, 30>();
496 test_three_quadrature_schemes_examples<cpp_bin_float_quad, 30>();
497 test_integration_over_real_line<cpp_bin_float_quad, 30>();
498 test_right_limit_infinite<cpp_bin_float_quad, 30>();
499 test_left_limit_infinite<cpp_bin_float_quad, 30>();
500
501
502 #endif
503 #ifdef TEST3
504 test_left_limit_infinite<cpp_bin_float_quad, 30>();
505 test_complex_lambert_w<std::complex<double>>();
506 test_complex_lambert_w<std::complex<long double>>();
507 #ifdef BOOST_HAS_FLOAT128
508 test_left_limit_infinite<boost::multiprecision::float128, 30>();
509 test_complex_lambert_w<boost::multiprecision::complex128>();
510 #endif
511 test_complex_lambert_w<boost::multiprecision::cpp_complex_quad>();
512 #endif
513 }
514
515 #else
516
main()517 int main() { return 0; }
518
519 #endif
520