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