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