ceph/src/boost/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 <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_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>
  22
  23 #if !defined(TEST1) && !defined(TEST1A) && !defined(TEST2) && !defined(TEST3)
  24 #  define TEST1
  25 #  define TEST1A
  26 #  define TEST2
  27 #  define TEST3
  28 #endif
  29
  30 #ifdef _MSC_VER
  31 #pragma warning(disable:4127)  // Conditional expression is constant
  32 #endif
  33
  34 using std::expm1;
  35 using std::atan;
  36 using std::tan;
  37 using std::log;
  38 using std::log1p;
  39 using std::asinh;
  40 using std::atanh;
  41 using std::sqrt;
  42 using std::isnormal;
  43 using std::abs;
  44 using std::sinh;
  45 using std::tanh;
  46 using std::cosh;
  47 using std::pow;
  48 using std::exp;
  49 using std::sin;
  50 using std::cos;
  51 using std::string;
  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;
  70
  71 //
  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:
  74 //
  75
  76 enum
  77 {
  78    test_ca_error_id,
  79    test_ca_error_id_2,
  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
  85 };
  86
  87 template <unsigned Points>
  88 double expected_error(unsigned)
  89 {
  90    return 0; // placeholder, all tests will fail
  91 }
  92
  93 template <>
  94 double expected_error<15>(unsigned id)
  95 {
  96    switch (id)
  97    {
  98    case test_ca_error_id:
  99       return 1e-7;
 100    case test_ca_error_id_2:
 101       return 2e-5;
 102    case test_three_quad_error_id:
 103       return 1e-8;
 104    case test_three_quad_error_id_2:
 105       return 3.5e-3;
 106    case test_integration_over_real_line_error_id:
 107       return 6e-3;
 108    case test_right_limit_infinite_error_id:
 109    case test_left_limit_infinite_error_id:
 110       return 1e-5;
 111    }
 112    return 0;  // placeholder, all tests will fail
 113 }
 114
 115 template <>
 116 double expected_error<17>(unsigned id)
 117 {
 118    switch (id)
 119    {
 120    case test_ca_error_id:
 121       return 1e-7;
 122    case test_ca_error_id_2:
 123       return 2e-5;
 124    case test_three_quad_error_id:
 125       return 1e-8;
 126    case test_three_quad_error_id_2:
 127       return 3.5e-3;
 128    case test_integration_over_real_line_error_id:
 129       return 6e-3;
 130    case test_right_limit_infinite_error_id:
 131    case test_left_limit_infinite_error_id:
 132       return 1e-5;
 133    }
 134    return 0;  // placeholder, all tests will fail
 135 }
 136
 137 template <>
 138 double expected_error<21>(unsigned id)
 139 {
 140    switch (id)
 141    {
 142    case test_ca_error_id:
 143       return 1e-12;
 144    case test_ca_error_id_2:
 145       return 3e-6;
 146    case test_three_quad_error_id:
 147       return 2e-13;
 148    case test_three_quad_error_id_2:
 149       return 2e-3;
 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:
 154       return 5e-8;
 155    }
 156    return 0;  // placeholder, all tests will fail
 157 }
 158
 159 template <>
 160 double expected_error<31>(unsigned id)
 161 {
 162    switch (id)
 163    {
 164    case test_ca_error_id:
 165       return 6e-20;
 166    case test_ca_error_id_2:
 167       return 3e-7;
 168    case test_three_quad_error_id:
 169       return 1e-19;
 170    case test_three_quad_error_id_2:
 171       return 6e-4;
 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:
 176       return 5e-11;
 177    }
 178    return 0;  // placeholder, all tests will fail
 179 }
 180
 181 template <>
 182 double expected_error<41>(unsigned id)
 183 {
 184    switch (id)
 185    {
 186    case test_ca_error_id:
 187       return 1e-26;
 188    case test_ca_error_id_2:
 189       return 1e-7;
 190    case test_three_quad_error_id:
 191       return 3e-27;
 192    case test_three_quad_error_id_2:
 193       return 3e-4;
 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:
 198       return 1e-15;
 199    }
 200    return 0;  // placeholder, all tests will fail
 201 }
 202
 203 template <>
 204 double expected_error<51>(unsigned id)
 205 {
 206    switch (id)
 207    {
 208    case test_ca_error_id:
 209       return 5e-33;
 210    case test_ca_error_id_2:
 211       return 1e-8;
 212    case test_three_quad_error_id:
 213       return 1e-32;
 214    case test_three_quad_error_id_2:
 215       return 3e-4;
 216    case test_integration_over_real_line_error_id:
 217       return 1e-14;
 218    case test_right_limit_infinite_error_id:
 219    case test_left_limit_infinite_error_id:
 220       return 3e-19;
 221    }
 222    return 0;  // placeholder, all tests will fail
 223 }
 224
 225 template <>
 226 double expected_error<61>(unsigned id)
 227 {
 228    switch (id)
 229    {
 230    case test_ca_error_id:
 231       return 5e-34;
 232    case test_ca_error_id_2:
 233       return 5e-9;
 234    case test_three_quad_error_id:
 235       return 4e-34;
 236    case test_three_quad_error_id_2:
 237       return 1e-4;
 238    case test_integration_over_real_line_error_id:
 239       return 1e-16;
 240    case test_right_limit_infinite_error_id:
 241    case test_left_limit_infinite_error_id:
 242       return 3e-23;
 243    }
 244    return 0;  // placeholder, all tests will fail
 245 }
 246
 247
 248 template<class Real, unsigned Points>
 249 void test_linear()
 250 {
 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;
 253     Real error;
 254     auto f = [](const Real& x)
 255     {
 256        return 5*x + 7;
 257     };
 258     Real L1;
 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);
 262 }
 263
 264 template<class Real, unsigned Points>
 265 void test_quadratic()
 266 {
 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;
 269     Real error;
 270
 271     auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
 272     Real L1;
 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);
 276 }
 277
 278 // Examples taken from
 279 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
 280 template<class Real, unsigned Points>
 281 void test_ca()
 282 {
 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);
 285     Real L1;
 286     Real error;
 287
 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);
 293
 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);
 299
 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);
 305 }
 306
 307 template<class Real, unsigned Points>
 308 void test_three_quadrature_schemes_examples()
 309 {
 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);
 312     Real Q;
 313     Real Q_expected;
 314
 315     // Example 1:
 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);
 320
 321
 322     // Example 2:
 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);
 327
 328     // Example 3:
 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);
 333
 334     // Example 4:
 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);
 339
 340     tol = expected_error<Points>(test_three_quad_error_id_2);
 341     // Example 5:
 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);
 346
 347     // Example 6:
 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);
 352 }
 353
 354
 355 template<class Real, unsigned Points>
 356 void test_integration_over_real_line()
 357 {
 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);
 360     Real Q;
 361     Real Q_expected;
 362     Real L1;
 363     Real error;
 364
 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);
 370 }
 371
 372 template<class Real, unsigned Points>
 373 void test_right_limit_infinite()
 374 {
 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);
 377     Real Q;
 378     Real Q_expected;
 379     Real L1;
 380     Real error;
 381
 382     // Example 11:
 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);
 387
 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);
 392 }
 393
 394 template<class Real, unsigned Points>
 395 void test_left_limit_infinite()
 396 {
 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);
 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_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);
 407 }
 408
 409 BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
 410 {
 411 #ifdef TEST1
 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>();
 420
 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>();
 430 #endif
 431 #ifdef TEST1A
 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>();
 440
 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>();
 449 #endif
 450 #ifdef TEST2
 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>();
 459
 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>();
 468 #endif
 469 #ifdef TEST3
 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>();
 479 #endif
 480 }
 481
 482 #else
 483
 484 int main() { return 0; }
 485
 486 #endif