ceph/src/boost/libs/math/test/exp_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 exp_sinh_quadrature_test
   8
   9 #include <boost/math/concepts/real_concept.hpp>
  10 #include <boost/test/included/unit_test.hpp>
  11 #include <boost/test/floating_point_comparison.hpp>
  12 #include <boost/math/quadrature/exp_sinh.hpp>
  13 #include <boost/math/special_functions/sinc.hpp>
  14 #include <boost/math/special_functions/bessel.hpp>
  15 #include <boost/multiprecision/cpp_bin_float.hpp>
  16 #include <boost/multiprecision/cpp_dec_float.hpp>
  17 #include <boost/math/special_functions/next.hpp>
  18 #include <boost/math/special_functions/gamma.hpp>
  19 #include <boost/math/special_functions/sinc.hpp>
  20 #include <boost/type_traits/is_class.hpp>
  21
  22 using std::exp;
  23 using std::cos;
  24 using std::tan;
  25 using std::log;
  26 using std::sqrt;
  27 using std::abs;
  28 using std::sinh;
  29 using std::cosh;
  30 using std::pow;
  31 using std::atan;
  32 using boost::multiprecision::cpp_bin_float_50;
  33 using boost::multiprecision::cpp_bin_float_100;
  34 using boost::multiprecision::cpp_bin_float_quad;
  35 using boost::math::constants::pi;
  36 using boost::math::constants::half_pi;
  37 using boost::math::constants::two_div_pi;
  38 using boost::math::constants::half;
  39 using boost::math::constants::third;
  40 using boost::math::constants::half;
  41 using boost::math::constants::third;
  42 using boost::math::constants::catalan;
  43 using boost::math::constants::ln_two;
  44 using boost::math::constants::root_two;
  45 using boost::math::constants::root_two_pi;
  46 using boost::math::constants::root_pi;
  47 using boost::math::quadrature::exp_sinh;
  48
  49 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7)
  50 #  define TEST1
  51 #  define TEST2
  52 #  define TEST3
  53 #  define TEST4
  54 #  define TEST5
  55 #  define TEST6
  56 #  define TEST7
  57 #endif
  58
  59 #ifdef BOOST_MSVC
  60 #pragma warning (disable:4127)
  61 #endif
  62
  63 //
  64 // Coefficient generation code:
  65 //
  66 template <class T>
  67 void print_levels(const T& v, const char* suffix)
  68 {
  69    std::cout << "{\n";
  70    for (unsigned i = 0; i < v.size(); ++i)
  71    {
  72       std::cout << "      { ";
  73       for (unsigned j = 0; j < v[i].size(); ++j)
  74       {
  75          std::cout << v[i][j] << suffix << ", ";
  76       }
  77       std::cout << "},\n";
  78    }
  79    std::cout << "   };\n";
  80 }
  81
  82 template <class T>
  83 void print_levels(const std::pair<T, T>& p, const char* suffix = "")
  84 {
  85    std::cout << "   static const std::vector<std::vector<Real> > abscissa = ";
  86    print_levels(p.first, suffix);
  87    std::cout << "   static const std::vector<std::vector<Real> > weights = ";
  88    print_levels(p.second, suffix);
  89 }
  90
  91 template <class Real, class TargetType>
  92 std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
  93 {
  94    using boost::math::constants::half_pi;
  95    using boost::math::constants::two_div_pi;
  96    using boost::math::constants::pi;
  97    auto g = [](Real t) { return exp(half_pi<Real>()*sinh(t)); };
  98    auto w = [](Real t) { return cosh(t)*half_pi<Real>()*exp(half_pi<Real>()*sinh(t)); };
  99
 100    std::vector<std::vector<Real>> abscissa, weights;
 101
 102    std::vector<Real> temp;
 103
 104    Real tmp = (Real(boost::math::tools::log_min_value<TargetType>()) + log(Real(boost::math::tools::epsilon<TargetType>())))*0.5f;
 105    Real t_min = asinh(two_div_pi<Real>()*tmp);
 106    // truncate t_min to an exact binary value:
 107    t_min = floor(t_min * 128) / 128;
 108
 109    std::cout << "m_t_min = " << t_min << ";\n";
 110
 111    // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g).
 112    // This will allow some flexibility on the users part; they can at least square a number function without overflow.
 113    // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands.
 114    const Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(Real(boost::math::tools::max_value<TargetType>()))));
 115
 116    Real h = 1;
 117    for (Real t = t_min; t < t_max; t += h)
 118    {
 119       temp.push_back(g(t));
 120    }
 121    abscissa.push_back(temp);
 122    temp.clear();
 123
 124    for (Real t = t_min; t < t_max; t += h)
 125    {
 126       temp.push_back(w(t * h));
 127    }
 128    weights.push_back(temp);
 129    temp.clear();
 130
 131    for (unsigned row = 1; row < max_rows; ++row)
 132    {
 133       h /= 2;
 134       for (Real t = t_min + h; t < t_max; t += 2 * h)
 135          temp.push_back(g(t));
 136       abscissa.push_back(temp);
 137       temp.clear();
 138    }
 139    h = 1;
 140    for (unsigned row = 1; row < max_rows; ++row)
 141    {
 142       h /= 2;
 143       for (Real t = t_min + h; t < t_max; t += 2 * h)
 144          temp.push_back(w(t));
 145       weights.push_back(temp);
 146       temp.clear();
 147    }
 148
 149    return std::make_pair(abscissa, weights);
 150 }
 151
 152
 153 template <class Real>
 154 const exp_sinh<Real>& get_integrator()
 155 {
 156    static const exp_sinh<Real> integrator(14);
 157    return integrator;
 158 }
 159
 160 template <class Real>
 161 Real get_convergence_tolerance()
 162 {
 163    return boost::math::tools::root_epsilon<Real>();
 164 }
 165
 166 template<class Real>
 167 void test_right_limit_infinite()
 168 {
 169     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";
 170     Real tol = 10 * boost::math::tools::epsilon<Real>();
 171     Real Q;
 172     Real Q_expected;
 173     Real error;
 174     Real L1;
 175     auto integrator = get_integrator<Real>();
 176
 177     // Example 12
 178     const auto f2 = [](const Real& t) { return exp(-t)/sqrt(t); };
 179     Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
 180     Q_expected = root_pi<Real>();
 181     Real tol_mult = 1;
 182     // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
 183     if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
 184        tol_mult = 12;
 185     else if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
 186        tol_mult = 5;
 187     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * tol_mult);
 188     // The integrand is strictly positive, so it coincides with the value of the integral:
 189     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * tol_mult);
 190
 191     auto f3 = [](Real t)->Real { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); };
 192     Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
 193     Q_expected = half<Real>();
 194     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 195     Q = integrator.integrate(f3, 10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 196     Q_expected = boost::lexical_cast<Real>("-6.6976341310426674140007086979326069121526743314567805278252392932e-6");
 197     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
 198     // Integrating through zero risks precision loss:
 199     Q = integrator.integrate(f3, -10, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 200     Q_expected = boost::lexical_cast<Real>("-15232.3213626280525704332288302799653087046646639974940243044623285817777006");
 201     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, std::numeric_limits<Real>::digits10 > 30 ? 1000 * tol : tol);
 202
 203     auto f4 = [](Real t)->Real { return 1/(1+t*t); };
 204     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);
 205     Q_expected = pi<Real>()/4;
 206     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 207     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 208     Q = integrator.integrate(f4, 20, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 209     Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
 210     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 211     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 212     Q = integrator.integrate(f4, 500, std::numeric_limits<Real>::has_infinity ? std::numeric_limits<Real>::infinity() : boost::math::tools::max_value<Real>(), get_convergence_tolerance<Real>(), &error, &L1);
 213     Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
 214     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 215     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 216 }
 217
 218 template<class Real>
 219 void test_left_limit_infinite()
 220 {
 221     std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 222     Real tol = 10 * boost::math::tools::epsilon<Real>();
 223     Real Q;
 224     Real Q_expected;
 225     Real error;
 226     Real L1;
 227     auto integrator = get_integrator<Real>();
 228
 229     // Example 11:
 230     auto f1 = [](const Real& t) { return 1/(1+t*t);};
 231     Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), 0, get_convergence_tolerance<Real>(), &error, &L1);
 232     Q_expected = half_pi<Real>();
 233     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 234     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 235     Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -20, get_convergence_tolerance<Real>(), &error, &L1);
 236     Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
 237     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 238     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 239     Q = integrator.integrate(f1, std::numeric_limits<Real>::has_infinity ? -std::numeric_limits<Real>::infinity() : -boost::math::tools::max_value<Real>(), -500, get_convergence_tolerance<Real>(), &error, &L1);
 240     Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
 241     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 242     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 243 }
 244
 245
 246 // Some examples of tough integrals from NR, section 4.5.4:
 247 template<class Real>
 248 void test_nr_examples()
 249 {
 250     using std::sin;
 251     using std::cos;
 252     using std::pow;
 253     using std::exp;
 254     using std::sqrt;
 255     std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 256     Real tol = 10 * boost::math::tools::epsilon<Real>();
 257     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
 258     Real Q;
 259     Real Q_expected;
 260     Real L1;
 261     Real error;
 262     auto integrator = get_integrator<Real>();
 263
 264     auto f0 = [] (Real)->Real { return (Real) 0; };
 265     Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
 266     Q_expected = 0;
 267     BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol);
 268     BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol);
 269
 270     auto f = [](const Real& x) { return 1/(1+x*x); };
 271     Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
 272     Q_expected = half_pi<Real>();
 273     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 274     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
 275
 276     auto f1 = [](Real x)->Real {
 277         Real z1 = exp(-x);
 278         if (z1 == 0)
 279         {
 280             return (Real) 0;
 281         }
 282         Real z2 = pow(x, -3*half<Real>())*z1;
 283         if (z2 == 0)
 284         {
 285             return (Real) 0;
 286         }
 287         return sin(x*half<Real>())*z2;
 288     };
 289
 290     Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
 291     Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2));
 292
 293     // The integrand is oscillatory; the accuracy is low.
 294     Real tol_mul = 1;
 295     if (std::numeric_limits<Real>::digits10 > 40)
 296        tol_mul = 500000;
 297     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
 298
 299     auto f2 = [](Real x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(pow(x, -(Real) 2/(Real) 7)*exp(-x*x)); };
 300     Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
 301     Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14);
 302     tol_mul = 1;
 303     if (std::numeric_limits<Real>::is_specialized == false)
 304        tol_mul = 3;
 305     if (std::numeric_limits<Real>::digits10 > 40)
 306        tol_mul = 100;
 307     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
 308     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);
 309
 310     auto f3 = [](Real x)->Real { return (Real) 1/ (sqrt(x)*(1+x)); };
 311     Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
 312     Q_expected = pi<Real>();
 313
 314     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon<Real>());
 315     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon<Real>());
 316
 317     auto f4 = [](const Real& t) { return  t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t*half<Real>())); };
 318     Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
 319     Q_expected = root_two_pi<Real>()/2;
 320     tol_mul = 1;
 321     if (std::numeric_limits<Real>::digits10 > 40)
 322        tol_mul = 5000;
 323     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
 324     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);
 325
 326     auto f5 = [](const Real& t) { return 1/cosh(t);};
 327     Q = integrator.integrate(f5, get_convergence_tolerance<Real>(), &error, &L1);
 328     Q_expected = half_pi<Real>();
 329     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 12);   // Fails at float precision without higher error rate
 330     BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 12);
 331 }
 332
 333 // Definite integrals found in the CRC Handbook of Mathematical Formulas
 334 template<class Real>
 335 void test_crc()
 336 {
 337     using std::sin;
 338     using std::pow;
 339     using std::exp;
 340     using std::sqrt;
 341     using std::log;
 342     using std::cos;
 343     std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 344     Real tol = 10 * boost::math::tools::epsilon<Real>();
 345     std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
 346     Real Q;
 347     Real Q_expected;
 348     Real L1;
 349     Real error;
 350     auto integrator = get_integrator<Real>();
 351
 352     auto f0 = [](const Real& x) { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(log(x)*exp(-x)); };
 353     Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
 354     Q_expected = -boost::math::constants::euler<Real>();
 355     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 356
 357     // Test the integral representation of the gamma function:
 358     auto f1 = [](Real t)->Real { Real x = exp(-t);
 359         if(x == 0)
 360         {
 361             return (Real) 0;
 362         }
 363         return pow(t, (Real) 12 - 1)*x;
 364     };
 365
 366     Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
 367     Q_expected = boost::math::tgamma(12.0f);
 368     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 369
 370     // Integral representation of the modified bessel function:
 371     // K_5(12)
 372     auto f2 = [](Real t)->Real {
 373         Real x = 12*cosh(t);
 374         if (x > boost::math::tools::log_max_value<Real>())
 375         {
 376             return (Real) 0;
 377         }
 378         return exp(-x)*cosh(5*t);
 379     };
 380     Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
 381     Q_expected = boost::math::cyl_bessel_k<int, Real>(5, 12);
 382     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 383     // Laplace transform of cos(at)
 384     Real a = 20;
 385     Real s = 1;
 386     auto f3 = [&](Real t)->Real {
 387         Real x = s * t;
 388         if (x > boost::math::tools::log_max_value<Real>())
 389         {
 390             return (Real) 0;
 391         }
 392         return cos(a * t) * exp(-x);
 393     };
 394
 395     // For high oscillation frequency, the quadrature sum is ill-conditioned.
 396     Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
 397     Q_expected = s/(a*a+s*s);
 398     // Since the integrand is oscillatory, we increase the tolerance:
 399     Real tol_mult = 10;
 400     // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
 401     if (!boost::is_class<Real>::value)
 402     {
 403        if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
 404           tol_mult = 5000; // we should really investigate this more??
 405        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult*tol);
 406     }
 407
 408     //
 409     // This one doesn't pass for real_concept..
 410     //
 411     if (std::numeric_limits<Real>::is_specialized)
 412     {
 413        // Laplace transform of J_0(t):
 414        auto f4 = [&](Real t)->Real {
 415           Real x = s * t;
 416           if (x > boost::math::tools::log_max_value<Real>())
 417           {
 418              return (Real)0;
 419           }
 420           return boost::math::cyl_bessel_j(0, t) * exp(-x);
 421        };
 422
 423        Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
 424        Q_expected = 1 / sqrt(1 + s*s);
 425        tol_mult = 3;
 426        // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
 427        if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
 428           tol_mult = 750;
 429        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult * tol);
 430     }
 431     auto f6 = [](const Real& t) { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t)*log(t));};
 432     Q = integrator.integrate(f6, get_convergence_tolerance<Real>(), &error, &L1);
 433     Q_expected = -boost::math::constants::root_pi<Real>()*(boost::math::constants::euler<Real>() + 2*ln_two<Real>())/4;
 434     BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 435
 436     // CRC Section 5.5, integral 591
 437     // The parameter p allows us to control the strength of the singularity.
 438     // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
 439     // This converges only when our test type has an extended exponent range as all the area of the integral
 440     // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
 441     // "There's a lot of room at the bottom".
 442     // This version is transformed via argument substitution (exp(-x) for x) so that the integral is spread
 443     // over (0, INF).
 444     tol *= boost::math::tools::digits<Real>() > 100 ? 100000 : 75;
 445     for (Real pn = 99; pn > 0; pn -= 10) {
 446        Real p = pn / 100;
 447        auto f = [&](Real x)->Real
 448        {
 449           return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-x * (1 - p) + p * log(-boost::math::expm1(-x))));
 450        };
 451        Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
 452        Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
 453        /*
 454        std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << p << std::endl;
 455        std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q << std::endl;
 456        std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q_expected << std::endl;
 457        std::cout << fabs((Q - Q_expected) / Q_expected) << std::endl;
 458        */
 459        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 460     }
 461     // and for p < 1:
 462     for (Real p = -0.99; p < 0; p += 0.1) {
 463        auto f = [&](Real x)->Real
 464        {
 465           return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x));
 466        };
 467        Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
 468        Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
 469        BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
 470     }
 471 }
 472
 473
 474 BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
 475 {
 476    //
 477    // Uncomment to generate the coefficients:
 478    //
 479
 480    /*
 481    std::cout << std::scientific << std::setprecision(8);
 482    print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
 483    std::cout << std::setprecision(18);
 484    print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
 485    std::cout << std::setprecision(35);
 486    print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
 487    */
 488
 489
 490 #ifdef TEST1
 491     test_left_limit_infinite<float>();
 492     test_right_limit_infinite<float>();
 493     test_nr_examples<float>();
 494     test_crc<float>();
 495 #endif
 496 #ifdef TEST2
 497     test_left_limit_infinite<double>();
 498     test_right_limit_infinite<double>();
 499     test_nr_examples<double>();
 500     test_crc<double>();
 501 #endif
 502 #ifdef TEST3
 503     test_left_limit_infinite<long double>();
 504     test_right_limit_infinite<long double>();
 505     test_nr_examples<long double>();
 506     test_crc<long double>();
 507 #endif
 508 #ifdef TEST4
 509     test_left_limit_infinite<cpp_bin_float_quad>();
 510     test_right_limit_infinite<cpp_bin_float_quad>();
 511     test_nr_examples<cpp_bin_float_quad>();
 512     test_crc<cpp_bin_float_quad>();
 513 #endif
 514
 515 #ifdef TEST5
 516     test_left_limit_infinite<boost::math::concepts::real_concept>();
 517     test_right_limit_infinite<boost::math::concepts::real_concept>();
 518     test_nr_examples<boost::math::concepts::real_concept>();
 519     test_crc<boost::math::concepts::real_concept>();
 520 #endif
 521
 522 #ifdef TEST6
 523     test_left_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
 524     test_right_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
 525     test_nr_examples<boost::multiprecision::cpp_bin_float_50>();
 526     test_crc<boost::multiprecision::cpp_bin_float_50>();
 527 #endif
 528 #ifdef TEST7
 529     test_left_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
 530     test_right_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
 531     test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
 532     test_crc<boost::multiprecision::cpp_dec_float_50>();
 533 #endif
 534 }