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