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