...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
template<class Real> class tanh_sinh { public: tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4) template<class F> auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)>decltype(std::declval<F>()(std::declval<Real>())) const; template<class F> auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)>decltype(std::declval<F>()(std::declval<Real>())) const; };
The tanhsinh
quadrature routine provided by boost
is a rapidly convergent numerical integration scheme for holomorphic integrands.
By this we mean that the integrand is the restriction to the real line of
a complexdifferentiable function which is bounded on the interior of the
unit disk z < 1, so that it lies within the socalled
Hardy space.
If your integrand obeys these conditions, it can be shown that tanhsinh
integration is optimal, in the sense that it requires the fewest function
evaluations for a given accuracy of any quadrature algorithm for a random
element from the Hardy space.
A basic example of how to use the tanhsinh
quadrature
is shown below:
tanh_sinh<double> integrator; auto f = [](double x) { return 5*x + 7; }; // Integrate over native bounds of (1,1): double Q = integrator.integrate(f); // Integrate over (0,1.1) instead: Q = integrator.integrate(f, 0.0, 1.1);
The basic idea of tanhsinh
quadrature is that a variable transformation
can cause the endpoint derivatives to decay rapidly. When the derivatives
at the endpoints decay much faster than the Bernoulli numbers grow, the EulerMaclaurin
summation formula tells us that simple trapezoidal quadrature converges faster
than any power of h. That means the number of correct
digits of the result should roughly double with each new level of integration
(halving of h), Hence the default termination condition
for integration is usually set to the square root of machine epsilon. Most
wellbehaved integrals should converge to full machine precision with this
termination condition, and in 6 or fewer levels at double precision, or 7
or fewer levels for quad precision.
One very nice property of tanhsinh quadrature is that it can handle singularities at the endpoints of the integration domain. For instance, the following integrand, singular at both endpoints, can be efficiently evaluated to 100 binary digits:
auto f = [](Real x) { return log(x)*log1p(x); }; Real Q = integrator.integrate(f, (Real) 0, (Real) 1);
Now onto the caveats: As stated before, the integrands must lie in a Hardy space to ensure rapid convergence. Attempting to integrate a function which is not bounded on the unit disk by tanhsinh can lead to very slow convergence. For example, take the Runge function:
auto f1 = [](double t) { return 1/(1+25*t*t); }; Q = integrator.integrate(f1, (double) 1, (double) 1);
This function has poles at ± ⅈ/5, and as such it is not bounded on the unit disk. However, the related function
auto f2 = [](double t) { return 1/(1+0.04*t*t); }; Q = integrator.integrate(f2, (double) 1, (double) 1);
has poles outside the unit disk (at ± 5ⅈ), and is therefore in the Hardy space. Our benchmarks show that the second integration is performed 22x faster than the first! If you do not understand the structure of your integrand in the complex plane, do performance testing before deployment.
Like the trapezoidal quadrature, the tanhsinh quadrature produces an estimate of the L_{1} norm of the integral along with the requested integral. This is to establish a scale against which to measure the tolerance, and to provide an estimate of the condition number of the summation. This can be queried as follows:
tanh_sinh<double> integrator; auto f = [](double x) { return 5*x + 7; }; double termination = std::sqrt(std::numeric_limits<double>::epsilon()); double error; double L1; size_t levels; double Q = integrator.integrate(f, 0.0, 1.0, termination, &error, &L1, &levels); double condition_number = L1/std::abs(Q);
If the condition number is large, the computed integral is worthless: typically one can assume that Q has lost one digit of precision when the condition number of O(10^Q). The returned error term is not the actual error in the result, but merely an a posteriori error estimate. It is the absolute difference between the last two approximations, and for well behaved integrals, the actual error should be very much smaller than this. The following table illustrates how the errors and conditioning vary for few sample integrals, in each case the termination condition was set to the square root of epsilon, and all tests were conducted in double precision:
Integral 
Range 
Error 
Actual measured error 
Levels 
Condition Number 
Comments 


(0,1) 
3.5e15 
0 
5 
1 
This trivial case shows just how accurate these methods can be. 

0, 1) 
0 
0 
5 
1 
This is an example of an integral that Gaussian integrators fail to handle. 

(0,+∞) 
8.0e10 
1.1e15 
5 
1 
Gaussian integrators typically fail to handle the singularities at the endpoints of this one. 

(1,1) 
7.2e16 
4.9e17 
9 
1.89 
This is a truly horrible integral that oscillates wildly and unpredictably with some very sharp "spikes" in it's graph. The higher number of levels used reflects the difficulty of sampling the more extreme features. 

(∞, ∞) 
3.0e1 
4.0e1 
15 
159 
This highly oscillatory integral isn't handled at all well by tanhsinh quadrature: there is so much cancellation in the sum that the result is essentially worthless. The argument transformation of the infinite integral behaves somewhat badly as well, in fact we do slightly better integrating over 2 symmetrical and large finite limits. 

(0,1) 
1e8 
1e8 
5 
1 
This an example of an integral that has all its area close to a nonzero endpoint, the problem here is that the function being integrated returns "garbage" values for x very close to 1. We can easily fix this issue by passing a 2 argument functor to the integrator: the second argument gives the distance to the nearest endpoint, and we can use that information to return accurate values, and thus fix the integral calculation. 

(0,1) 
0 
0 
5 
1 
This is the 2argument version of the previous integral, the second
argument xc is 
Although the tanhsinh
quadrature can compute integral over
infinite domains by variable transformations, these transformations can create
a very poorly behaved integrand. For this reason, doubleexponential variable
transformations have been provided that allow stable computation over infinite
domains; these being the expsinh and sinhsinh quadrature.
The tanh_sinh
integrator
supports integration of functions which return complex results, for example
the sineintegral Si(z)
has
the integral representation:
Which we can code up directly as:
template <class Complex> Complex Si(Complex z) { typedef typename Complex::value_type value_type; using std::sin; using std::cos; using std::exp; auto f = [&z](value_type t) { return exp(z * cos(t)) * cos(z * sin(t)); }; boost::math::quadrature::tanh_sinh<value_type> integrator; return integrator.integrate(f, 0, boost::math::constants::half_pi<value_type>()) + boost::math::constants::half_pi<value_type>(); }