Boost C++ Libraries

template<class Real>
class exp_sinh
{
public:
    exp_sinh(size_t max_refinements = 9);

    template<class F>
    auto integrate(const F f, Real a, Real b,
                   Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
                   Real* error = nullptr,
                   Real* L1 = nullptr,
                   size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
    template<class F>
    auto integrate(const F f,
                   Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
                   Real* error = nullptr,
                   Real* L1 = nullptr,
                   size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
};

For half-infinite intervals, the exp-sinh quadrature is provided:

exp_sinh<double> integrator;
auto f = [](double x) { return exp(-3*x); };
double termination = sqrt(std::numeric_limits<double>::epsilon());
double error;
double L1;
double Q = integrator.integrate(f, termination, &error, &L1);

The native integration range of this integrator is (0, ∞), but we also support /(a, ∞), (-∞, 0)/ and /(-∞, b)/ via argument transformations.

Endpoint singularities and complex-valued integrands are supported by exp-sinh.

For example, the modified Bessel function K can be represented via:

Which we can code up as:

template <class Complex>
Complex bessel_K(Complex alpha, Complex z)
{
   typedef typename Complex::value_type value_type;
   using std::cosh;  using std::exp;
   auto f = [&, alpha, z](value_type t)
   {
      value_type ct = cosh(t);
      if (ct > log(std::numeric_limits<value_type>::max()))
         return Complex(0);
      return exp(-z * ct) * cosh(alpha * t);
   };
   boost::math::quadrature::exp_sinh<value_type> integrator;
   return integrator.integrate(f);
}

The only wrinkle in the above code is the need to check for large cosh(t) in which case we assume that exp(-x cosh(t)) tends to zero faster than cosh(alpha x) tends to infinity and return 0. Without that check we end up with 0 * Infinity as the result (a NaN).