libs/math/test/tanh_sinh_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/tools/floating_point_comparison.hpp>
  17 #include <boost/math/quadrature/tanh_sinh.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/math/special_functions/next.hpp>
  22 #include <boost/math/special_functions/gamma.hpp>
  23 #include <boost/math/special_functions/beta.hpp>
  24 #include <boost/math/special_functions/ellint_rc.hpp>
  25 #include <boost/math/special_functions/ellint_rj.hpp>
  26
  27 #ifdef BOOST_HAS_FLOAT128
  28 #include <boost/multiprecision/float128.hpp>
  29 #endif
  30
  31 #ifdef _MSC_VER
  32 #pragma warning(disable:4127)  // Conditional expression is constant
  33 #endif
  34
  35 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8)\
  36     && !defined(TEST1A) && !defined(TEST1B) && !defined(TEST2A) && !defined(TEST3A) && !defined(TEST6A) && !defined(TEST9)
  37 #  define TEST1
  38 #  define TEST2
  39 #  define TEST3
  40 #  define TEST4
  41 #  define TEST5
  42 #  define TEST6
  43 #  define TEST7
  44 #  define TEST8
  45 #  define TEST1A
  46 #  define TEST1B
  47 #  define TEST2A
  48 #  define TEST3A
  49 #  define TEST6A
  50 #  define TEST9
  51 #endif
  52
  53 using std::expm1;
  54 using std::atan;
  55 using std::tan;
  56 using std::log;
  57 using std::log1p;
  58 using std::asinh;
  59 using std::atanh;
  60 using std::sqrt;
  61 using std::isnormal;
  62 using std::abs;
  63 using std::sinh;
  64 using std::tanh;
  65 using std::cosh;
  66 using std::pow;
  67 using std::exp;
  68 using std::sin;
  69 using std::cos;
  70 using std::string;
  71 using boost::multiprecision::cpp_bin_float_50;
  72 using boost::multiprecision::cpp_bin_float_100;
  73 using boost::multiprecision::cpp_dec_float_50;
  74 using boost::multiprecision::cpp_dec_float_100;
  75 using boost::multiprecision::cpp_bin_float_quad;
  76 using boost::math::sinc_pi;
  77 using boost::math::quadrature::tanh_sinh;
  78 using boost::math::quadrature::detail::tanh_sinh_detail;
  79 using boost::math::constants::pi;
  80 using boost::math::constants::half_pi;
  81 using boost::math::constants::two_div_pi;
  82 using boost::math::constants::two_pi;
  83 using boost::math::constants::half;
  84 using boost::math::constants::third;
  85 using boost::math::constants::half;
  86 using boost::math::constants::third;
  87 using boost::math::constants::catalan;
  88 using boost::math::constants::ln_two;
  89 using boost::math::constants::root_two;
  90 using boost::math::constants::root_two_pi;
  91 using boost::math::constants::root_pi;
  92
  93 template <class T>
  94 void print_levels(const T& v, const char* suffix)
  95 {
  96    std::cout << "{\n";
  97    for (unsigned i = 0; i < v.size(); ++i)
  98    {
  99       std::cout << "      { ";
 100       for (unsigned j = 0; j < v[i].size(); ++j)
 101       {
 102          std::cout << v[i][j] << suffix << ", ";
 103       }
 104       std::cout << "},\n";
 105    }
 106    std::cout << "   };\n";
 107 }
 108
 109 template <class T>
 110 void print_complement_indexes(const T& v)
 111 {
 112    std::cout << "\n   {";
 113    for (unsigned i = 0; i < v.size(); ++i)
 114    {
 115       unsigned index = 0;
 116       while (v[i][index] >= 0)
 117          ++index;
 118       std::cout << index << ", ";
 119    }
 120    std::cout << "};\n";
 121 }
 122
 123 template <class T>
 124 void print_levels(const std::pair<T, T>& p, const char* suffix = "")
 125 {
 126    std::cout << "   static const std::vector<std::vector<Real> > abscissa = ";
 127    print_levels(p.first, suffix);
 128    std::cout << "   static const std::vector<std::vector<Real> > weights = ";
 129    print_levels(p.second, suffix);
 130    std::cout << "   static const std::vector<unsigned > indexes = ";
 131    print_complement_indexes(p.first);
 132 }
 133
 134 template <class Real>
 135 std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_index, unsigned max_rows)
 136 {
 137    using boost::math::constants::half_pi;
 138    using boost::math::constants::two_div_pi;
 139    using boost::math::constants::pi;
 140    auto g = [](Real t) { return tanh(half_pi<Real>()*sinh(t)); };
 141    auto w = [](Real t) { Real cs = cosh(half_pi<Real>() * sinh(t));  return half_pi<Real>() * cosh(t) / (cs * cs); };
 142    auto gc = [](Real t) { Real u2 = half_pi<Real>() * sinh(t);  return 1 / (exp(u2) *cosh(u2)); };
 143    auto g_inv = [](float x)->float { return asinh(two_div_pi<float>()*atanh(x)); };
 144    auto gc_inv = [](float x)
 145    {
 146       float l = log(sqrt((2 - x) / x));
 147       return log((sqrt(4 * l * l + pi<float>() * pi<float>()) + 2 * l) / pi<float>());
 148    };
 149
 150    std::vector<std::vector<Real>> abscissa, weights;
 151
 152    std::vector<Real> temp;
 153
 154    float t_crossover = gc_inv(0.5f);
 155
 156    Real h = 1;
 157    for (unsigned i = 0; i < max_index; ++i)
 158    {
 159       temp.push_back((i < t_crossover) ? g(i * h) : -gc(i * h));
 160    }
 161    abscissa.push_back(temp);
 162    temp.clear();
 163
 164    for (unsigned i = 0; i < max_index; ++i)
 165    {
 166       temp.push_back(w(i * h));
 167    }
 168    weights.push_back(temp);
 169    temp.clear();
 170
 171    for (unsigned row = 1; row < max_rows; ++row)
 172    {
 173       h /= 2;
 174       for (Real t = h; t < max_index - 1; t += 2 * h)
 175          temp.push_back((t < t_crossover) ? g(t) : -gc(t));
 176       abscissa.push_back(temp);
 177       temp.clear();
 178    }
 179    h = 1;
 180    for (unsigned row = 1; row < max_rows; ++row)
 181    {
 182       h /= 2;
 183       for (Real t = h; t < max_index - 1; t += 2 * h)
 184          temp.push_back(w(t));
 185       weights.push_back(temp);
 186       temp.clear();
 187    }
 188
 189    return std::make_pair(abscissa, weights);
 190 }
 191
 192 template <class Real>
 193 const tanh_sinh<Real>& get_integrator()
 194 {
 195    static const tanh_sinh<Real> integrator(15);
 196    return integrator;
 197 }
 198
 199 template <class Real>
 200 Real get_convergence_tolerance()
 201 {
 202    return boost::math::tools::root_epsilon<Real>();
 203 }
 204
 205
 206 template<class Real>
 207 void test_linear()
 208 {
 209     std::cout << "Testing linear functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 210     Real tol = 10*boost::math::tools::epsilon<Real>();
 211     auto integrator = get_integrator<Real>();
 212     auto f = [](const Real& x)
 213     {
 214        return 5*x + 7;
 215     };
 216     Real error;
 217     Real L1;
 218     Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 219     BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
 220     BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
 221 }
 222
 223
 224 template<class Real>
 225 void test_quadratic()
 226 {
 227     std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 228     Real tol = 10*boost::math::tools::epsilon<Real>();
 229     auto integrator = get_integrator<Real>();
 230     auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
 231     Real error;
 232     Real L1;
 233     Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 234     BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
 235     BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
 236 }
 237
 238
 239 template<class Real>
 240 void test_singular()
 241 {
 242     std::cout << "Testing integration of endpoint singularities on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 243     Real tol = 10*boost::math::tools::epsilon<Real>();
 244     Real error;
 245     Real L1;
 246     auto integrator = get_integrator<Real>();
 247     auto f = [](const Real& x) { return log(x)*log(1-x); };
 248     Real Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 249     Real Q_expected = 2 - pi<Real>()*pi<Real>()*half<Real>()*third<Real>();
 250
 251     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 252     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 253 }
 254
 255
 256 // Examples taken from
 257 //http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
 258 template<class Real>
 259 void test_ca()
 260 {
 261     std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 262     Real tol = 10 * boost::math::tools::epsilon<Real>();
 263     Real error;
 264     Real L1;
 265
 266     auto integrator = get_integrator<Real>();
 267     auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
 268     Real Q = integrator.integrate(f1, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 269     Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
 270     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 271     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 272
 273     auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
 274     Q = integrator.integrate(f2, (Real) 0 , (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 275     Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
 276     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 277     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 278
 279     auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
 280     Q = integrator.integrate(f5, (Real) 0 , (Real) 1);
 281     Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
 282     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 283
 284
 285     // Oh it suffers on this one.
 286     auto f6 = [](Real t)->Real { Real x = log(t); return x*x; };
 287     Q = integrator.integrate(f6, (Real) 0 , (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 288     Q_expected = 2;
 289
 290     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 50*tol);
 291     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 50*tol);
 292
 293
 294     // Although it doesn't get to the requested tolerance on this integral, the error bounds can be queried and are reasonable:
 295     tol = sqrt(boost::math::tools::epsilon<Real>());
 296     auto f7 = [](const Real& t) { return sqrt(tan(t)); };
 297     Q = integrator.integrate(f7, (Real) 0 , (Real) half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 298     Q_expected = pi<Real>()/root_two<Real>();
 299     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 300     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 301
 302     auto f8 = [](const Real& t) { return log(cos(t)); };
 303     Q = integrator.integrate(f8, (Real) 0 , half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 304     Q_expected = -pi<Real>()*ln_two<Real>()*half<Real>();
 305     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 306     BOOST_CHECK_CLOSE_FRACTION(L1, -Q_expected, tol);
 307 }
 308
 309
 310 template<class Real>
 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 = 10 * boost::math::tools::epsilon<Real>();
 315     Real Q;
 316     Real Q_expected;
 317
 318     auto integrator = get_integrator<Real>();
 319     // Example 1:
 320     auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
 321     Q = integrator.integrate(f1, (Real) 0 , (Real) 1);
 322     Q_expected = half<Real>()*half<Real>();
 323     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 324
 325
 326     // Example 2:
 327     auto f2 = [](const Real& t) { return t*t*atan(t); };
 328     Q = integrator.integrate(f2, (Real) 0 , (Real) 1);
 329     Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
 330     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
 331
 332     // Example 3:
 333     auto f3 = [](const Real& t) { return exp(t)*cos(t); };
 334     Q = integrator.integrate(f3, (Real) 0, half_pi<Real>());
 335     Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
 336     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 337
 338     // Example 4:
 339     auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
 340     Q = integrator.integrate(f4, (Real) 0 , (Real) 1);
 341     Q_expected = 5*pi<Real>()*pi<Real>()/96;
 342     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 343
 344     // Example 5:
 345     auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
 346     Q = integrator.integrate(f5, (Real) 0 , (Real) 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 = integrator.integrate(f6, (Real) 0 , (Real) 1);
 353     Q_expected = pi<Real>()/4;
 354     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 355 }
 356
 357
 358 template<class Real>
 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 = 10 * boost::math::tools::epsilon<Real>();
 363     Real Q;
 364     Real Q_expected;
 365     Real error;
 366     Real L1;
 367     auto integrator = get_integrator<Real>();
 368
 369     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 370     Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 371     Q_expected = pi<Real>();
 372     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 373     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 374
 375     auto f2 = [](const Real& t) { return exp(-t*t*half<Real>()); };
 376     Q = integrator.integrate(f2, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 377     Q_expected = root_two_pi<Real>();
 378     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 2);
 379     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 2);
 380
 381     // This test shows how oscillatory integrals are approximated very poorly by this method:
 382     //std::cout << "Testing sinc integral: \n";
 383     //Q = integrator.integrate(boost::math::sinc_pi<Real>, -std::numeric_limits<Real>::infinity(), std::numeric_limits<Real>::infinity(), &error, &L1);
 384     //std::cout << "Error estimate of sinc integral: " << error << std::endl;
 385     //std::cout << "L1 norm of sinc integral " << L1 << std::endl;
 386     //Q_expected = pi<Real>();
 387     //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 388
 389     auto f4 = [](const Real& t) { return 1/cosh(t);};
 390     Q = integrator.integrate(f4, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 391     Q_expected = pi<Real>();
 392     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 393     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 394
 395 }
 396
 397 template<class Real>
 398 void test_right_limit_infinite()
 399 {
 400     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";
 401     Real tol = 10 * boost::math::tools::epsilon<Real>();
 402     Real Q;
 403     Real Q_expected;
 404     Real error;
 405     Real L1;
 406     auto integrator = get_integrator<Real>();
 407
 408     // Example 11:
 409     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 410     Q = integrator.integrate(f1, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 411     Q_expected = half_pi<Real>();
 412     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 413
 414     // Example 12
 415     auto f2 = [](const Real& t) { return exp(-t)/sqrt(t); };
 416     Q = integrator.integrate(f2, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 417     Q_expected = root_pi<Real>();
 418     BOOST_CHECK_CLOSE(Q, Q_expected, 1000*tol);
 419
 420     auto f3 = [](const Real& t) { return exp(-t)*cos(t); };
 421     Q = integrator.integrate(f3, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 422     Q_expected = half<Real>();
 423     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 424
 425     auto f4 = [](const Real& t) { return 1/(1+t*t); };
 426     Q = integrator.integrate(f4, 1, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 427     Q_expected = pi<Real>()/4;
 428     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 429 }
 430
 431 template<class Real>
 432 void test_left_limit_infinite()
 433 {
 434     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";
 435     Real tol = 10 * boost::math::tools::epsilon<Real>();
 436     Real Q;
 437     Real Q_expected;
 438     auto integrator = get_integrator<Real>();
 439
 440     // Example 11:
 441     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 442     Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), Real(0));
 443     Q_expected = half_pi<Real>();
 444     BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
 445 }
 446
 447
 448 // A horrible function taken from
 449 // http://www.chebfun.org/examples/quad/GaussClenCurt.html
 450 template<class Real>
 451 void test_horrible()
 452 {
 453     std::cout << "Testing Trefenthen's horrible integral on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 454     // We only know the integral to double precision, so requesting a higher tolerance doesn't make sense.
 455     Real tol = 10 * std::numeric_limits<float>::epsilon();
 456     Real Q;
 457     Real Q_expected;
 458     Real error;
 459     Real L1;
 460     auto integrator = get_integrator<Real>();
 461
 462     auto f = [](Real x)->Real { return x*sin(2*exp(2*sin(2*exp(2*x) ) ) ); };
 463     Q = integrator.integrate(f, (Real) -1, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 464     // NIntegrate[x*Sin[2*Exp[2*Sin[2*Exp[2*x]]]], {x, -1, 1}, WorkingPrecision -> 130, MaxRecursion -> 100]
 465     Q_expected = boost::lexical_cast<Real>("0.33673283478172753598559003181355241139806404130031017259552729882281");
 466     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 467     // Over again without specifying the bounds:
 468     Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
 469     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 470 }
 471
 472 // Some examples of tough integrals from NR, section 4.5.4:
 473 template<class Real>
 474 void test_nr_examples()
 475 {
 476     std::cout << "Testing singular integrals from NR 4.5.4 on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 477     Real tol = 10 * boost::math::tools::epsilon<Real>();
 478     Real Q;
 479     Real Q_expected;
 480     Real error;
 481     Real L1;
 482     auto integrator = get_integrator<Real>();
 483
 484     auto f1 = [](Real x)->Real
 485     {
 486        return (sin(x * half<Real>()) * exp(-x) / x) / sqrt(x);
 487     };
 488     Q = integrator.integrate(f1, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 489     Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2));
 490     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 25*tol);
 491
 492     auto f2 = [](Real x)->Real { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); };
 493     Q = integrator.integrate(f2, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>());
 494     Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14);
 495     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 6);
 496
 497 }
 498
 499 // Test integrand known to fool some termination schemes:
 500 template<class Real>
 501 void test_early_termination()
 502 {
 503     std::cout << "Testing Clenshaw & Curtis's example of integrand which fools termination schemes on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 504     Real tol = 10 * boost::math::tools::epsilon<Real>();
 505     Real Q;
 506     Real Q_expected;
 507     Real error;
 508     Real L1;
 509     auto integrator = get_integrator<Real>();
 510
 511     auto f1 = [](Real x)->Real { return 23*cosh(x)/ (Real) 25 - cos(x) ; };
 512     Q = integrator.integrate(f1, (Real) -1, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 513     Q_expected = 46*sinh((Real) 1)/(Real) 25 - 2*sin((Real) 1);
 514     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 515     // Over again with no bounds:
 516     Q = integrator.integrate(f1);
 517     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 518 }
 519
 520 // Test some definite integrals from the CRC handbook, 32nd edition:
 521 template<class Real>
 522 void test_crc()
 523 {
 524     std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 525     Real tol = 10 * boost::math::tools::epsilon<Real>();
 526     Real Q;
 527     Real Q_expected;
 528     Real error;
 529     Real L1;
 530     auto integrator = get_integrator<Real>();
 531
 532     // CRC Definite integral 585
 533     auto f1 = [](Real x)->Real { Real t = log(1/x); return x*x*t*t*t; };
 534     Q = integrator.integrate(f1, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 535     Q_expected = (Real) 2/ (Real) 27;
 536     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 537
 538     // CRC 636:
 539     auto f2 = [](Real x)->Real { return sqrt(cos(x)); };
 540     Q = integrator.integrate(f2, (Real) 0, (Real) half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 541     //Q_expected = pow(two_pi<Real>(), 3*half<Real>())/pow(boost::math::tgamma((Real) 1/ (Real) 4), 2);
 542     Q_expected = boost::lexical_cast<Real>("1.198140234735592207439922492280323878227212663215651558263674952946405214143915670835885556489793389375907225");
 543     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 544
 545     // CRC Section 5.5, integral 585:
 546     for (int n = 0; n < 3; ++n) {
 547         for (int m = 0; m < 3; ++m) {
 548             auto f = [&](Real x)->Real { return pow(x, Real(m))*pow(log(1/x), Real(n)); };
 549             Q = integrator.integrate(f, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 550             // Calculation of the tgamma function is not exact, giving spurious failures.
 551             // Casting to cpp_bin_float_100 beforehand fixes most of them.
 552             cpp_bin_float_100 np1 = n + 1;
 553             cpp_bin_float_100 mp1 = m + 1;
 554             Q_expected = boost::lexical_cast<Real>((tgamma(np1)/pow(mp1, np1)).str());
 555             BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 556         }
 557     }
 558
 559     // CRC Section 5.5, integral 591
 560     // The parameter p allows us to control the strength of the singularity.
 561     // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
 562     // This converges only when our test type has an extended exponent range as all the area of the integral
 563     // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
 564     // "There's a lot of room at the bottom".
 565     // We also use a 2 argument functor so that 1-x is evaluated accurately:
 566     if (std::numeric_limits<Real>::max_exponent > std::numeric_limits<double>::max_exponent)
 567     {
 568        for (Real p = Real (-0.99); p < 1; p += Real(0.1)) {
 569           auto f = [&](Real x, Real cx)->Real
 570           {
 571              //return pow(x, p) / pow(1 - x, p);
 572              return cx < 0 ? exp(p * (log(x) - boost::math::log1p(-x))) : pow(x, p) / pow(cx, p);
 573           };
 574           Q = integrator.integrate(f, (Real)0, (Real)1, get_convergence_tolerance<Real>(), &error, &L1);
 575           Q_expected = 1 / sinc_pi(p*pi<Real>());
 576           BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
 577        }
 578     }
 579     // There is an alternative way to evaluate the above integral: by noticing that all the area of the integral
 580     // is near zero for p < 0 and near 1 for p > 0 we can substitute exp(-x) for x and remap the integral to the
 581     // domain (0, INF).  Internally we need to expand out the terms and evaluate using logs to avoid spurous overflow,
 582     // this gives us
 583     // for p > 0:
 584     for (Real p = Real(0.99); p > 0; p -= Real(0.1)) {
 585        auto f = [&](Real x)->Real
 586        {
 587           return exp(-x * (1 - p) + p * log(-boost::math::expm1(-x)));
 588        };
 589        Q = integrator.integrate(f, 0, boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 590        Q_expected = 1 / sinc_pi(p*pi<Real>());
 591        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
 592     }
 593     // and for p < 1:
 594     for (Real p = Real (-0.99); p < 0; p += Real(0.1)) {
 595        auto f = [&](Real x)->Real
 596        {
 597           return exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x);
 598        };
 599        Q = integrator.integrate(f, 0, boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 600        Q_expected = 1 / sinc_pi(p*pi<Real>());
 601        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
 602     }
 603
 604     // CRC Section 5.5, integral 635
 605     for (int m = 0; m < 10; ++m) {
 606         auto f = [&](Real x)->Real { return 1/(1 + pow(tan(x), m)); };
 607         Q = integrator.integrate(f, (Real) 0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 608         Q_expected = half_pi<Real>()/2;
 609         BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 610     }
 611
 612     // CRC Section 5.5, integral 637:
 613     //
 614     // When h gets very close to 1, the strength of the singularity gradually increases until we
 615     // no longer have enough exponent range to evaluate the integral....
 616     // We also have to use the 2-argument functor version of the integrator to avoid
 617     // cancellation error, since the singularity is near PI/2.
 618     //
 619     Real limit = std::numeric_limits<Real>::max_exponent > std::numeric_limits<double>::max_exponent
 620        ? .95f : .85f;
 621     for (Real h = Real(0.01); h < limit; h += Real(0.1)) {
 622         auto f = [&](Real x, Real xc)->Real { return xc > 0 ? pow(1/tan(xc), h) : pow(tan(x), h); };
 623         Q = integrator.integrate(f, (Real) 0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 624         Q_expected = half_pi<Real>()/cos(h*half_pi<Real>());
 625         BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 626     }
 627     // CRC Section 5.5, integral 637:
 628     //
 629     // Over again, but with argument transformation, we want:
 630     //
 631     // Integral of tan(x)^h over (0, PI/2)
 632     //
 633     // Note that the bulk of the area is next to the singularity at PI/2,
 634     // so we'll start by replacing x by PI/2 - x, and that tan(PI/2 - x) == 1/tan(x)
 635     // so we now have:
 636     //
 637     // Integral of 1/tan(x)^h over (0, PI/2)
 638     //
 639     // Which is almost the ideal form, except that when h is very close to 1
 640     // we run out of exponent range in evaluating the integral arbitrarily close to 0.
 641     // So now we substitute exp(-x) for x: this stretches out the range of the
 642     // integral to (-log(PI/2), +INF) with the singularity at +INF giving:
 643     //
 644     // Integral of exp(-x)/tan(exp(-x))^h over (-log(PI/2), +INF)
 645     //
 646     // We just need a way to evaluate the function without spurious under/overflow
 647     // in the exp terms.  Note that for small x: tan(x) ~= x, so making this
 648     // substitution and evaluating by logs we have:
 649     //
 650     // exp(-x)/tan(exp(-x))^h ~= exp((h - 1) * x)  for x > -log(epsion);
 651     //
 652     // Here's how that looks in code:
 653     //
 654     for (Real i = 80; i < 100; ++i) {
 655        Real h = i / 100;
 656        auto f = [&](Real x)->Real
 657        {
 658           if (x > -log(boost::math::tools::epsilon<Real>()))
 659              return exp((h - 1) * x);
 660           else
 661           {
 662              Real et = exp(-x);
 663              // Need to deal with numeric instability causing et to be greater than PI/2:
 664              return et >= boost::math::constants::half_pi<Real>() ? 0 : et * pow(1 / tan(et), h);
 665           }
 666        };
 667        Q = integrator.integrate(f, -log(half_pi<Real>()), boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 668        Q_expected = half_pi<Real>() / cos(h*half_pi<Real>());
 669        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 5 * tol);
 670     }
 671
 672     // CRC Section 5.5, integral 670:
 673
 674     auto f3 = [](Real x)->Real { return sqrt(log(1/x)); };
 675     Q = integrator.integrate(f3, (Real) 0, (Real) 1, get_convergence_tolerance<Real>(), &error, &L1);
 676     Q_expected = root_pi<Real>()/2;
 677     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 678
 679 }
 680
 681 template <class Real>
 682 void test_sf()
 683 {
 684    using std::sqrt;
 685    // Test some special functions that we already know how to evaluate:
 686    std::cout << "Testing special functions on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 687    Real tol = 10 * boost::math::tools::epsilon<Real>();
 688    auto integrator = get_integrator<Real>();
 689
 690    // incomplete beta:
 691    if (std::numeric_limits<Real>::digits10 < 37) // Otherwise too slow
 692    {
 693       Real a(100), b(15);
 694       auto f = [&](Real x)->Real { return boost::math::ibeta_derivative(a, b, x); };
 695       BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, Real(0.25)), boost::math::ibeta(100, 15, Real(0.25)), tol * 10);
 696       // Check some really extreme versions:
 697       a = 1000;
 698       b = 500;
 699       BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, 1), Real(1), tol * 15);
 700       //
 701       // This is as extreme as we can get in this domain: otherwise the function has all it's
 702       // area so close to zero we never get in there no matter how many levels we go down:
 703       //
 704       a = Real(1) / 15;
 705       b = 30;
 706       BOOST_CHECK_CLOSE_FRACTION(integrator.integrate(f, 0, 1), Real(1), tol * 25);
 707    }
 708
 709    Real x = 2, y = 3, z = 0.5, p = 0.25;
 710    // Elliptic integral RC:
 711    Real Q = integrator.integrate([&](const Real& t)->Real { return 1 / (sqrt(t + x) * (t + y)); }, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>()) / 2;
 712    BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::ellint_rc(x, y), tol);
 713    // Elliptic Integral RJ:
 714    BOOST_CHECK_CLOSE_FRACTION(Real((Real(3) / 2) * integrator.integrate([&](Real t)->Real { return 1 / (sqrt((t + x) * (t + y) * (t + z)) * (t + p)); }, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>())), boost::math::ellint_rj(x, y, z, p), tol);
 715
 716    z = 5.5;
 717    if (std::numeric_limits<Real>::digits10 > 40)
 718       tol *= 200;
 719    else if (!std::numeric_limits<Real>::is_specialized)
 720       tol *= 3;
 721    // tgamma expressed as an integral:
 722    BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t)->Real { using std::pow; using std::exp; return t > 10000 ? Real(0) : Real(pow(t, z - 1) * exp(-t)); }, 0, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>()),
 723       boost::math::tgamma(z), tol);
 724    BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](const Real& t)->Real {  using std::exp; return exp(-t*t); }, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>()),
 725       boost::math::constants::root_pi<Real>(), tol);
 726 }
 727
 728 template <class Real>
 729 void test_2_arg()
 730 {
 731    BOOST_MATH_STD_USING
 732    std::cout << "Testing 2 argument functors on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 733    Real tol = 10 * boost::math::tools::epsilon<Real>();
 734
 735    auto integrator = get_integrator<Real>();
 736
 737    //
 738    // There are a whole family of integrals of the general form
 739    // x(1-x)^-N ; N < 1
 740    // which have all the interesting behaviour near the 2 singularities
 741    // and all converge, see: http://www.wolframalpha.com/input/?i=integrate+(x+*+(1-x))%5E-1%2FN+from+0+to+1
 742    //
 743    Real Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 744    {
 745       return tc < 0 ? 1 / sqrt(t * (1-t)) : 1 / sqrt(t * tc);
 746    }, 0, 1);
 747    BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
 748    Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 749    {
 750       return tc < 0 ? 1 / boost::math::cbrt(t * (1-t)) : 1 / boost::math::cbrt(t * tc);
 751    }, 0, 1);
 752    BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::pow<2>(boost::math::tgamma(Real(2) / 3)) / boost::math::tgamma(Real(4) / 3), tol * 3);
 753    //
 754    // We can do the same thing with ((1+x)(1-x))^-N ; N < 1
 755    //
 756    Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 757    {
 758       return t < 0 ? 1 / sqrt(-tc * (1-t)) : 1 / sqrt((t + 1) * tc);
 759    }, -1, 1);
 760    BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
 761    Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 762    {
 763       return t < 0 ? 1 / sqrt(-tc * (1-t)) : 1 / sqrt((t + 1) * tc);
 764    });
 765    BOOST_CHECK_CLOSE_FRACTION(Q, boost::math::constants::pi<Real>(), tol);
 766    Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 767    {
 768       return t < 0 ? 1 / boost::math::cbrt(-tc * (1-t)) : 1 / boost::math::cbrt((t + 1) * tc);
 769    }, -1, 1);
 770    BOOST_CHECK_CLOSE_FRACTION(Q, sqrt(boost::math::constants::pi<Real>()) * boost::math::tgamma(Real(2) / 3) / boost::math::tgamma(Real(7) / 6), tol * 10);
 771    Q = integrator.integrate([&](const Real& t, const Real & tc)->Real
 772    {
 773       return t < 0 ? 1 / boost::math::cbrt(-tc * (1-t)) : 1 / boost::math::cbrt((t + 1) * tc);
 774    });
 775    BOOST_CHECK_CLOSE_FRACTION(Q, sqrt(boost::math::constants::pi<Real>()) * boost::math::tgamma(Real(2) / 3) / boost::math::tgamma(Real(7) / 6), tol * 10);
 776    //
 777    // These are taken from above, and do not get to full precision via the single arg functor:
 778    //
 779    auto f7 = [](const Real& t, const Real& tc) { return t < 1 ? sqrt(tan(t)) : sqrt(1/tan(tc)); };
 780    Real error, L1;
 781    Q = integrator.integrate(f7, (Real)0, (Real)half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 782    Real Q_expected = pi<Real>() / root_two<Real>();
 783    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 784    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 785
 786    auto f8 = [](const Real& t, const Real& tc) { return t < 1 ? log(cos(t)) : log(sin(tc)); };
 787    Q = integrator.integrate(f8, (Real)0, half_pi<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 788    Q_expected = -pi<Real>()*ln_two<Real>()*half<Real>();
 789    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 790    BOOST_CHECK_CLOSE_FRACTION(L1, -Q_expected, tol);
 791 }
 792
 793 template <class Complex>
 794 void test_complex()
 795 {
 796    typedef typename Complex::value_type value_type;
 797    //
 798    // Integral version of the confluent hypergeometric function:
 799    // https://dlmf.nist.gov/13.4#E1
 800    //
 801    value_type tol = 10 * boost::math::tools::epsilon<value_type>();
 802    Complex a(2, 3), b(3, 4), z(0.5, -2);
 803
 804    auto f = [&](value_type t)
 805    {
 806       return exp(z * t) * pow(t, a - value_type(1)) * pow(value_type(1) - t, b - a - value_type(1));
 807    };
 808
 809    auto integrator = get_integrator<value_type>();
 810    auto Q = integrator.integrate(f, value_type(0), value_type(1), get_convergence_tolerance<value_type>());
 811    //
 812    // Expected result computed from http://www.wolframalpha.com/input/?i=1F1%5B(2%2B3i),+(3%2B4i);+(0.5-2i)%5D+*+gamma(2%2B3i)+*+gamma(1%2Bi)+%2F+gamma(3%2B4i)
 813    //
 814    Complex expected(boost::lexical_cast<value_type>("-0.2911081612888249710582867318081776512805281815037891183828405999609246645054069649838607112484426042883371996"),
 815       boost::lexical_cast<value_type>("0.4507983563969959578849120188097153649211346293694903758252662015991543519595834937475296809912196906074655385"));
 816
 817    value_type error = abs(expected - Q);
 818    BOOST_CHECK_LE(error, tol);
 819
 820    //
 821    // Sin Integral https://dlmf.nist.gov/6.2#E9
 822    //
 823    auto f2 = [z](value_type t)
 824    {
 825       return -exp(-z * cos(t)) * cos(z * sin(t));
 826    };
 827    Q = integrator.integrate(f2, value_type(0), boost::math::constants::half_pi<value_type>(), get_convergence_tolerance<value_type>());
 828
 829    expected = Complex(boost::lexical_cast<value_type>("0.8893822921008980697856313681734926564752476188106405688951257340480164694708337246829840859633322683740376134733"),
 830       -boost::lexical_cast<value_type>("2.381380802906111364088958767973164614925936185337231718483495612539455538280372745733208000514737758457795502168"));
 831    expected -= boost::math::constants::half_pi<value_type>();
 832
 833    error = abs(expected - Q);
 834    BOOST_CHECK_LE(error, tol);
 835 }
 836
 837
 838 BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
 839 {
 840 #ifdef GENERATE_CONSTANTS
 841    //
 842    // Generate pre-computed coefficients:
 843    std::cout << std::setprecision(35);
 844    print_levels(generate_constants<cpp_bin_float_100>(10, 8), "L");
 845
 846 #else
 847
 848 #ifdef TEST1
 849
 850     test_right_limit_infinite<float>();
 851     test_left_limit_infinite<float>();
 852     test_linear<float>();
 853     test_quadratic<float>();
 854     test_singular<float>();
 855     test_ca<float>();
 856     test_three_quadrature_schemes_examples<float>();
 857     test_horrible<float>();
 858     test_integration_over_real_line<float>();
 859     test_nr_examples<float>();
 860 #endif
 861 #ifdef TEST1A
 862     test_early_termination<float>();
 863     test_2_arg<float>();
 864 #endif
 865 #ifdef TEST1B
 866     test_crc<float>();
 867 #endif
 868 #ifdef TEST2
 869     test_right_limit_infinite<double>();
 870     test_left_limit_infinite<double>();
 871     test_linear<double>();
 872     test_quadratic<double>();
 873     test_singular<double>();
 874     test_ca<double>();
 875     test_three_quadrature_schemes_examples<double>();
 876     test_horrible<double>();
 877     test_integration_over_real_line<double>();
 878     test_nr_examples<double>();
 879     test_early_termination<double>();
 880     test_sf<double>();
 881     test_2_arg<double>();
 882 #endif
 883 #ifdef TEST2A
 884     test_crc<double>();
 885 #endif
 886
 887 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 888
 889 #ifdef TEST3
 890     test_right_limit_infinite<long double>();
 891     test_left_limit_infinite<long double>();
 892     test_linear<long double>();
 893     test_quadratic<long double>();
 894     test_singular<long double>();
 895     test_ca<long double>();
 896     test_three_quadrature_schemes_examples<long double>();
 897     test_horrible<long double>();
 898     test_integration_over_real_line<long double>();
 899     test_nr_examples<long double>();
 900     test_early_termination<long double>();
 901     test_sf<long double>();
 902     test_2_arg<long double>();
 903 #endif
 904 #ifdef TEST3A
 905     test_crc<long double>();
 906
 907 #endif
 908 #endif
 909
 910 #ifdef TEST4
 911     test_right_limit_infinite<cpp_bin_float_quad>();
 912     test_left_limit_infinite<cpp_bin_float_quad>();
 913     test_linear<cpp_bin_float_quad>();
 914     test_quadratic<cpp_bin_float_quad>();
 915     test_singular<cpp_bin_float_quad>();
 916     test_ca<cpp_bin_float_quad>();
 917     test_three_quadrature_schemes_examples<cpp_bin_float_quad>();
 918     test_horrible<cpp_bin_float_quad>();
 919     test_nr_examples<cpp_bin_float_quad>();
 920     test_early_termination<cpp_bin_float_quad>();
 921     test_crc<cpp_bin_float_quad>();
 922     test_sf<cpp_bin_float_quad>();
 923     test_2_arg<cpp_bin_float_quad>();
 924
 925 #endif
 926 #ifdef TEST5
 927
 928     test_sf<cpp_bin_float_50>();
 929     test_sf<cpp_bin_float_100>();
 930     test_sf<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<150> > >();
 931
 932 #endif
 933 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 934 #ifdef TEST6
 935
 936     test_right_limit_infinite<boost::math::concepts::real_concept>();
 937     test_left_limit_infinite<boost::math::concepts::real_concept>();
 938     test_linear<boost::math::concepts::real_concept>();
 939     test_quadratic<boost::math::concepts::real_concept>();
 940     test_singular<boost::math::concepts::real_concept>();
 941     test_ca<boost::math::concepts::real_concept>();
 942     test_three_quadrature_schemes_examples<boost::math::concepts::real_concept>();
 943     test_horrible<boost::math::concepts::real_concept>();
 944     test_integration_over_real_line<boost::math::concepts::real_concept>();
 945     test_nr_examples<boost::math::concepts::real_concept>();
 946     test_early_termination<boost::math::concepts::real_concept>();
 947     test_sf<boost::math::concepts::real_concept>();
 948     test_2_arg<boost::math::concepts::real_concept>();
 949 #endif
 950 #ifdef TEST6A
 951     test_crc<boost::math::concepts::real_concept>();
 952
 953 #endif
 954 #endif
 955 #ifdef TEST7
 956     test_sf<cpp_dec_float_50>();
 957 #endif
 958 #if defined(TEST8) && defined(BOOST_HAS_FLOAT128)
 959
 960     test_right_limit_infinite<boost::multiprecision::float128>();
 961     test_left_limit_infinite<boost::multiprecision::float128>();
 962     test_linear<boost::multiprecision::float128>();
 963     test_quadratic<boost::multiprecision::float128>();
 964     test_singular<boost::multiprecision::float128>();
 965     test_ca<boost::multiprecision::float128>();
 966     test_three_quadrature_schemes_examples<boost::multiprecision::float128>();
 967     test_horrible<boost::multiprecision::float128>();
 968     test_integration_over_real_line<boost::multiprecision::float128>();
 969     test_nr_examples<boost::multiprecision::float128>();
 970     test_early_termination<boost::multiprecision::float128>();
 971     test_crc<boost::multiprecision::float128>();
 972     test_sf<boost::multiprecision::float128>();
 973     test_2_arg<boost::multiprecision::float128>();
 974
 975 #endif
 976 #ifdef TEST9
 977     test_complex<std::complex<double> >();
 978     test_complex<std::complex<float> >();
 979 #endif
 980
 981
 982 #endif
 983 }
 984
 985 #else
 986
 987 int main() { return 0; }
 988
 989 #endif