123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596 |
- // Copyright Nick Thompson, 2017
- // Use, modification and distribution are subject to the
- // Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- #define BOOST_TEST_MODULE exp_sinh_quadrature_test
- #include <complex>
- #include <boost/multiprecision/cpp_complex.hpp>
- #include <boost/math/concepts/real_concept.hpp>
- #include <boost/test/included/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/quadrature/exp_sinh.hpp>
- #include <boost/math/special_functions/sinc.hpp>
- #include <boost/math/special_functions/bessel.hpp>
- #include <boost/multiprecision/cpp_bin_float.hpp>
- #include <boost/multiprecision/cpp_dec_float.hpp>
- #include <boost/math/special_functions/next.hpp>
- #include <boost/math/special_functions/gamma.hpp>
- #include <boost/math/special_functions/sinc.hpp>
- #include <boost/type_traits/is_class.hpp>
- #ifdef BOOST_HAS_FLOAT128
- #include <boost/multiprecision/complex128.hpp>
- #endif
- using std::exp;
- using std::cos;
- using std::tan;
- using std::log;
- using std::sqrt;
- using std::abs;
- using std::sinh;
- using std::cosh;
- using std::pow;
- using std::atan;
- using boost::multiprecision::cpp_bin_float_50;
- using boost::multiprecision::cpp_bin_float_100;
- using boost::multiprecision::cpp_bin_float_quad;
- using boost::math::constants::pi;
- using boost::math::constants::half_pi;
- using boost::math::constants::two_div_pi;
- using boost::math::constants::half;
- using boost::math::constants::third;
- using boost::math::constants::half;
- using boost::math::constants::third;
- using boost::math::constants::catalan;
- using boost::math::constants::ln_two;
- using boost::math::constants::root_two;
- using boost::math::constants::root_two_pi;
- using boost::math::constants::root_pi;
- using boost::math::quadrature::exp_sinh;
- #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) && !defined(TEST7) && !defined(TEST8)
- # define TEST1
- # define TEST2
- # define TEST3
- # define TEST4
- # define TEST5
- # define TEST6
- # define TEST7
- # define TEST8
- #endif
- #ifdef BOOST_MSVC
- #pragma warning (disable:4127)
- #endif
- //
- // Coefficient generation code:
- //
- template <class T>
- void print_levels(const T& v, const char* suffix)
- {
- std::cout << "{\n";
- for (unsigned i = 0; i < v.size(); ++i)
- {
- std::cout << " { ";
- for (unsigned j = 0; j < v[i].size(); ++j)
- {
- std::cout << v[i][j] << suffix << ", ";
- }
- std::cout << "},\n";
- }
- std::cout << " };\n";
- }
- template <class T>
- void print_levels(const std::pair<T, T>& p, const char* suffix = "")
- {
- std::cout << " static const std::vector<std::vector<Real> > abscissa = ";
- print_levels(p.first, suffix);
- std::cout << " static const std::vector<std::vector<Real> > weights = ";
- print_levels(p.second, suffix);
- }
- template <class Real, class TargetType>
- std::pair<std::vector<std::vector<Real>>, std::vector<std::vector<Real>> > generate_constants(unsigned max_rows)
- {
- using boost::math::constants::half_pi;
- using boost::math::constants::two_div_pi;
- using boost::math::constants::pi;
- auto g = [](Real t)->Real { return exp(half_pi<Real>()*sinh(t)); };
- auto w = [](Real t)->Real { return cosh(t)*half_pi<Real>()*exp(half_pi<Real>()*sinh(t)); };
- std::vector<std::vector<Real>> abscissa, weights;
- std::vector<Real> temp;
- Real tmp = (Real(boost::math::tools::log_min_value<TargetType>()) + log(Real(boost::math::tools::epsilon<TargetType>())))*0.5f;
- Real t_min = asinh(two_div_pi<Real>()*tmp);
- // truncate t_min to an exact binary value:
- t_min = floor(t_min * 128) / 128;
- std::cout << "m_t_min = " << t_min << ";\n";
- // t_max is chosen to make g'(t_max) ~ sqrt(max) (g' grows faster than g).
- // This will allow some flexibility on the users part; they can at least square a number function without overflow.
- // But there is no unique choice; the further out we can evaluate the function, the better we can do on slowly decaying integrands.
- const Real t_max = log(2 * two_div_pi<Real>()*log(2 * two_div_pi<Real>()*sqrt(Real(boost::math::tools::max_value<TargetType>()))));
- Real h = 1;
- for (Real t = t_min; t < t_max; t += h)
- {
- temp.push_back(g(t));
- }
- abscissa.push_back(temp);
- temp.clear();
- for (Real t = t_min; t < t_max; t += h)
- {
- temp.push_back(w(t * h));
- }
- weights.push_back(temp);
- temp.clear();
- for (unsigned row = 1; row < max_rows; ++row)
- {
- h /= 2;
- for (Real t = t_min + h; t < t_max; t += 2 * h)
- temp.push_back(g(t));
- abscissa.push_back(temp);
- temp.clear();
- }
- h = 1;
- for (unsigned row = 1; row < max_rows; ++row)
- {
- h /= 2;
- for (Real t = t_min + h; t < t_max; t += 2 * h)
- temp.push_back(w(t));
- weights.push_back(temp);
- temp.clear();
- }
- return std::make_pair(abscissa, weights);
- }
- template <class Real>
- const exp_sinh<Real>& get_integrator()
- {
- static const exp_sinh<Real> integrator(14);
- return integrator;
- }
- template <class Real>
- Real get_convergence_tolerance()
- {
- return boost::math::tools::root_epsilon<Real>();
- }
- template<class Real>
- void test_right_limit_infinite()
- {
- 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";
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- Real Q;
- Real Q_expected;
- Real error;
- Real L1;
- auto integrator = get_integrator<Real>();
- // Example 12
- const auto f2 = [](const Real& t)->Real { return exp(-t)/sqrt(t); };
- Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = root_pi<Real>();
- Real tol_mult = 1;
- // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
- if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
- tol_mult = 12;
- else if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
- tol_mult = 5;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * tol_mult);
- // The integrand is strictly positive, so it coincides with the value of the integral:
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * tol_mult);
- auto f3 = [](Real t)->Real { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); };
- Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = half<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- 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);
- Q_expected = boost::lexical_cast<Real>("-6.6976341310426674140007086979326069121526743314567805278252392932e-6");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10 * tol);
- // Integrating through zero risks precision loss:
- 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);
- Q_expected = boost::lexical_cast<Real>("-15232.3213626280525704332288302799653087046646639974940243044623285817777006");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, std::numeric_limits<Real>::digits10 > 30 ? 1000 * tol : tol);
- auto f4 = [](Real t)->Real { return 1/(1+t*t); };
- 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);
- Q_expected = pi<Real>()/4;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- 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);
- Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- 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);
- Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- }
- template<class Real>
- void test_left_limit_infinite()
- {
- std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- Real Q;
- Real Q_expected;
- Real error;
- Real L1;
- auto integrator = get_integrator<Real>();
- // Example 11:
- auto f1 = [](const Real& t)->Real { return 1/(1+t*t);};
- 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);
- Q_expected = half_pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- 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);
- Q_expected = boost::lexical_cast<Real>("0.0499583957219427614100062870348448814912770804235071744108534548299835954767");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- 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);
- Q_expected = boost::lexical_cast<Real>("0.0019999973333397333150476759363217553199063513829126652556286269630");
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- }
- // Some examples of tough integrals from NR, section 4.5.4:
- template<class Real>
- void test_nr_examples()
- {
- using std::sin;
- using std::cos;
- using std::pow;
- using std::exp;
- using std::sqrt;
- std::cout << "Testing type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
- Real Q;
- Real Q_expected;
- Real L1;
- Real error;
- auto integrator = get_integrator<Real>();
- auto f0 = [] (Real)->Real { return (Real) 0; };
- Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = 0;
- BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol);
- auto f = [](const Real& x)->Real { return 1/(1+x*x); };
- Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = half_pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
- auto f1 = [](Real x)->Real {
- Real z1 = exp(-x);
- if (z1 == 0)
- {
- return (Real) 0;
- }
- Real z2 = pow(x, -3*half<Real>())*z1;
- if (z2 == 0)
- {
- return (Real) 0;
- }
- return sin(x*half<Real>())*z2;
- };
- Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = sqrt(pi<Real>()*(sqrt((Real) 5) - 2));
- // The integrand is oscillatory; the accuracy is low.
- Real tol_mul = 1;
- if (std::numeric_limits<Real>::digits10 > 40)
- tol_mul = 500000;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
- 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)); };
- Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = half<Real>()*boost::math::tgamma((Real) 5/ (Real) 14);
- tol_mul = 1;
- if (std::numeric_limits<Real>::is_specialized == false)
- tol_mul = 6;
- else if (std::numeric_limits<Real>::digits10 > 40)
- tol_mul = 100;
- else
- tol_mul = 3;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);
- auto f3 = [](Real x)->Real { return (Real) 1/ (sqrt(x)*(1+x)); };
- Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon<Real>());
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon<Real>());
- auto f4 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t*half<Real>())); };
- Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = root_two_pi<Real>()/2;
- tol_mul = 1;
- if (std::numeric_limits<Real>::digits10 > 40)
- tol_mul = 5000;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mul * tol);
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol_mul * tol);
- auto f5 = [](const Real& t)->Real { return 1/cosh(t);};
- Q = integrator.integrate(f5, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = half_pi<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol * 12); // Fails at float precision without higher error rate
- BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol * 12);
- }
- // Definite integrals found in the CRC Handbook of Mathematical Formulas
- template<class Real>
- void test_crc()
- {
- using std::sin;
- using std::pow;
- using std::exp;
- using std::sqrt;
- using std::log;
- using std::cos;
- std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
- Real tol = 10 * boost::math::tools::epsilon<Real>();
- std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
- Real Q;
- Real Q_expected;
- Real L1;
- Real error;
- auto integrator = get_integrator<Real>();
- auto f0 = [](const Real& x)->Real { return x > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(log(x)*exp(-x)); };
- Q = integrator.integrate(f0, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = -boost::math::constants::euler<Real>();
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- // Test the integral representation of the gamma function:
- auto f1 = [](Real t)->Real { Real x = exp(-t);
- if(x == 0)
- {
- return (Real) 0;
- }
- return pow(t, (Real) 12 - 1)*x;
- };
- Q = integrator.integrate(f1, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = boost::math::tgamma(12.0f);
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- // Integral representation of the modified bessel function:
- // K_5(12)
- auto f2 = [](Real t)->Real {
- Real x = 12*cosh(t);
- if (x > boost::math::tools::log_max_value<Real>())
- {
- return (Real) 0;
- }
- return exp(-x)*cosh(5*t);
- };
- Q = integrator.integrate(f2, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = boost::math::cyl_bessel_k<int, Real>(5, 12);
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- // Laplace transform of cos(at)
- Real a = 20;
- Real s = 1;
- auto f3 = [&](Real t)->Real {
- Real x = s * t;
- if (x > boost::math::tools::log_max_value<Real>())
- {
- return (Real) 0;
- }
- return cos(a * t) * exp(-x);
- };
- // Since the integrand is oscillatory, we increase the tolerance:
- Real tol_mult = 10;
- // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
- if (!boost::is_class<Real>::value)
- {
- // For high oscillation frequency, the quadrature sum is ill-conditioned.
- Q = integrator.integrate(f3, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = s/(a*a+s*s);
- if (std::numeric_limits<Real>::digits10 > std::numeric_limits<double>::digits10)
- tol_mult = 5000; // we should really investigate this more??
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult*tol);
- }
- //
- // This one doesn't pass for real_concept..
- //
- if (std::numeric_limits<Real>::is_specialized)
- {
- // Laplace transform of J_0(t):
- auto f4 = [&](Real t)->Real {
- Real x = s * t;
- if (x > boost::math::tools::log_max_value<Real>())
- {
- return (Real)0;
- }
- return boost::math::cyl_bessel_j(0, t) * exp(-x);
- };
- Q = integrator.integrate(f4, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = 1 / sqrt(1 + s*s);
- tol_mult = 3;
- // Multiprecision type have higher error rates, probably evaluation of f() is less accurate:
- if (std::numeric_limits<Real>::digits10 > std::numeric_limits<long double>::digits10)
- tol_mult = 750;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol_mult * tol);
- }
- auto f6 = [](const Real& t)->Real { return t > boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-t*t)*log(t));};
- Q = integrator.integrate(f6, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = -boost::math::constants::root_pi<Real>()*(boost::math::constants::euler<Real>() + 2*ln_two<Real>())/4;
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- // CRC Section 5.5, integral 591
- // The parameter p allows us to control the strength of the singularity.
- // Rapid convergence is not guaranteed for this function, as the branch cut makes it non-analytic on a disk.
- // This converges only when our test type has an extended exponent range as all the area of the integral
- // occurs so close to 0 (or 1) that we need abscissa values exceptionally small to find it.
- // "There's a lot of room at the bottom".
- // This version is transformed via argument substitution (exp(-x) for x) so that the integral is spread
- // over (0, INF).
- tol *= boost::math::tools::digits<Real>() > 100 ? 100000 : 75;
- for (Real pn = 99; pn > 0; pn -= 10) {
- Real p = pn / 100;
- auto f = [&](Real x)->Real
- {
- return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-x * (1 - p) + p * log(-boost::math::expm1(-x))));
- };
- Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
- /*
- std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << p << std::endl;
- std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q << std::endl;
- std::cout << std::setprecision(std::numeric_limits<Real>::max_digits10) << Q_expected << std::endl;
- std::cout << fabs((Q - Q_expected) / Q_expected) << std::endl;
- */
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- }
- // and for p < 1:
- for (Real p = -0.99; p < 0; p += 0.1) {
- auto f = [&](Real x)->Real
- {
- return x > 1000 * boost::math::tools::log_max_value<Real>() ? Real(0) : Real(exp(-p * log(-boost::math::expm1(-x)) - (1 + p) * x));
- };
- Q = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
- Q_expected = 1 / boost::math::sinc_pi(p*pi<Real>());
- BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
- }
- }
- template<class Complex>
- void test_complex_modified_bessel()
- {
- std::cout << "Testing complex modified Bessel function on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
- typedef typename Complex::value_type Real;
- Real tol = 100 * boost::math::tools::epsilon<Real>();
- Real error;
- Real L1;
- auto integrator = get_integrator<Real>();
- // Integral Representation of Modified Complex Bessel function:
- // https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions
- Complex z{2, 3};
- const auto f = [&z](const Real& t)->Complex
- {
- using std::cosh;
- using std::exp;
- Real cosht = cosh(t);
- if (cosht > boost::math::tools::log_max_value<Real>())
- {
- return Complex{0, 0};
- }
- Complex arg = -z*cosht;
- Complex res = exp(arg);
- return res;
- };
- Complex K0 = integrator.integrate(f, get_convergence_tolerance<Real>(), &error, &L1);
- // Mathematica code: N[BesselK[0, 2 + 3 I], 140]
- Real K0_x_expected = boost::lexical_cast<Real>("-0.08296852656762551490517953520589186885781541203818846830385526187936132191822538822296497597191327722262903004145527496422090506197776994");
- Real K0_y_expected = boost::lexical_cast<Real>("0.027949603635183423629723306332336002340909030265538548521150904238352846705644065168365102147901993976999717171115546662967229050834575193041");
- BOOST_CHECK_CLOSE_FRACTION(K0.real(), K0_x_expected, tol);
- BOOST_CHECK_CLOSE_FRACTION(K0.imag(), K0_y_expected, tol);
- }
- BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
- {
- //
- // Uncomment to generate the coefficients:
- //
- /*
- std::cout << std::scientific << std::setprecision(8);
- print_levels(generate_constants<cpp_bin_float_100, float>(8), "f");
- std::cout << std::setprecision(18);
- print_levels(generate_constants<cpp_bin_float_100, double>(8), "");
- std::cout << std::setprecision(35);
- print_levels(generate_constants<cpp_bin_float_100, cpp_bin_float_quad>(8), "L");
- */
- #ifdef TEST1
- test_left_limit_infinite<float>();
- test_right_limit_infinite<float>();
- test_nr_examples<float>();
- test_crc<float>();
- #endif
- #ifdef TEST2
- test_left_limit_infinite<double>();
- test_right_limit_infinite<double>();
- test_nr_examples<double>();
- test_crc<double>();
- #endif
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- #ifdef TEST3
- test_left_limit_infinite<long double>();
- test_right_limit_infinite<long double>();
- test_nr_examples<long double>();
- test_crc<long double>();
- #endif
- #endif
- #ifdef TEST4
- test_left_limit_infinite<cpp_bin_float_quad>();
- test_right_limit_infinite<cpp_bin_float_quad>();
- test_nr_examples<cpp_bin_float_quad>();
- test_crc<cpp_bin_float_quad>();
- #endif
- #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
- #ifdef TEST5
- test_left_limit_infinite<boost::math::concepts::real_concept>();
- test_right_limit_infinite<boost::math::concepts::real_concept>();
- test_nr_examples<boost::math::concepts::real_concept>();
- test_crc<boost::math::concepts::real_concept>();
- #endif
- #endif
- #ifdef TEST6
- test_left_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
- test_right_limit_infinite<boost::multiprecision::cpp_bin_float_50>();
- test_nr_examples<boost::multiprecision::cpp_bin_float_50>();
- test_crc<boost::multiprecision::cpp_bin_float_50>();
- #endif
- #ifdef TEST7
- test_left_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
- test_right_limit_infinite<boost::multiprecision::cpp_dec_float_50>();
- test_nr_examples<boost::multiprecision::cpp_dec_float_50>();
- //
- // This one causes stack overflows on the CI machine, but not locally,
- // assume it's due to resticted resources on the server, and <shrug> for now...
- //
- #if ! BOOST_WORKAROUND(BOOST_MSVC, == 1900)
- test_crc<boost::multiprecision::cpp_dec_float_50>();
- #endif
- #endif
- #ifdef TEST8
- test_complex_modified_bessel<std::complex<float>>();
- test_complex_modified_bessel<std::complex<double>>();
- test_complex_modified_bessel<std::complex<long double>>();
- #ifdef BOOST_HAS_FLOAT128
- test_complex_modified_bessel<boost::multiprecision::complex128>();
- #endif
- test_complex_modified_bessel<boost::multiprecision::cpp_complex_quad>();
- #endif
- }
|