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 <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