2 Copyright (c) 2017 Nick Thompson
3 Use, modification and distribution are subject to the
4 Boost Software License, Version 1.0. (See accompanying file
5 LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
8 [section:double_exponential Double-exponential quadrature]
10 [section:de_overview Overview]
15 #include <boost/math/quadrature/tanh_sinh.hpp>
16 #include <boost/math/quadrature/exp_sinh.hpp>
17 #include <boost/math/quadrature/sinh_sinh.hpp>
19 namespace boost{ namespace math{
25 tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
28 auto integrate(const F f, Real a, Real b,
29 Real tolerance = tools::root_epsilon<Real>(),
30 Real* error = nullptr,
32 std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
35 auto integrate(const F f, Real
36 tolerance = tools::root_epsilon<Real>(),
37 Real* error = nullptr,
39 std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
47 exp_sinh(size_t max_refinements = 9);
50 auto integrate(const F f, Real a, Real b,
51 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
52 Real* error = nullptr,
54 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
56 auto integrate(const F f,
57 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
58 Real* error = nullptr,
60 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
67 sinh_sinh(size_t max_refinements = 9);
70 auto integrate(const F f,
71 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
72 Real* error = nullptr,
74 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
80 These three integration routines provide robust general purpose quadrature, each having a "native" range over which
81 quadrature is performed.
82 For example, the `sinh_sinh` quadrature integrates over the entire real line, the `tanh_sinh` over (-1, 1),
83 and the `exp_sinh` over (0, [infin]).
84 The latter integrators also have auxilliary ranges which are handled via a change of variables on the function being integrated,
85 so that the `tanh_sinh` can handle integration over /(a, b)/, and `exp_sinh` over /(a, [infin]) and(-[infin], b)/.
87 Like the other quadrature routines in Boost, these routines support both real and complex-valued integrands.
89 The `integrate` methods which do not specify a range always integrate over the native range of the method, and generally
90 are the most efficient and produce the smallest code, on the other hand the methods which do specify the bounds of integration
91 are the most general, and use argument transformations which are generally very robust. The following table summarizes
92 the ranges supported by each method:
95 [[Integrator][Native range][Other supported ranges][Comments]]
96 [[tanh_sinh] [(-1,1)] [(a,b)[br](a,[infin])[br](-[infin],b)[br](-[infin],[infin])]
97 [Special care is taken for endpoints at or near zero to ensure that abscissa values are calculated without the loss of precision
98 that would normally occur. Likewise when transforming to an infinite endpoint, the additional information which tanh_sinh has
99 internally on abscissa values is used to ensure no loss of precision during the transformation.]]
100 [[exp_sinh] [(0,[infin])] [(a,[infin])[br](-[infin],0)[br](-[infin],b)] []]
101 [[sinh_sinh] [(-[infin],[infin])] [][]]
104 [endsect] [/section:de_overview Overview]
106 [section:de_tanh_sinh tanh_sinh]
112 tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
115 auto integrate(const F f, Real a, Real b,
116 Real tolerance = tools::root_epsilon<Real>(),
117 Real* error = nullptr,
119 std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
122 auto integrate(const F f, Real
123 tolerance = tools::root_epsilon<Real>(),
124 Real* error = nullptr,
126 std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
130 The `tanh-sinh` quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands.
131 By this we mean that the integrand is the restriction to the real line of a complex-differentiable function which is bounded on the interior of the unit disk /|z| < 1/,
132 so that it lies within the so-called [@https://en.wikipedia.org/wiki/Hardy_space Hardy space].
133 If your integrand obeys these conditions, it can be shown that `tanh-sinh` integration is optimal,
134 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.
136 A basic example of how to use the `tanh-sinh` quadrature is shown below:
138 tanh_sinh<double> integrator;
139 auto f = [](double x) { return 5*x + 7; };
140 // Integrate over native bounds of (-1,1):
141 double Q = integrator.integrate(f);
142 // Integrate over (0,1.1) instead:
143 Q = integrator.integrate(f, 0.0, 1.1);
145 The basic idea of `tanh-sinh` quadrature is that a variable transformation can cause the endpoint derivatives to decay rapidly.
146 When the derivatives at the endpoints decay much faster than the Bernoulli numbers grow,
147 the Euler-Maclaurin summation formula tells us that simple trapezoidal quadrature converges faster than any power of /h/.
148 That means the number of correct digits of the result should roughly double with each new level of integration (halving of /h/),
149 Hence the default termination condition for integration is usually set to the square root of machine epsilon.
150 Most well-behaved integrals should converge to full machine precision with this termination condition,
151 and in 6 or fewer levels at double precision, or 7 or fewer levels for quad precision.
153 One very nice property of tanh-sinh quadrature is that it can handle singularities at the endpoints of the integration domain.
154 For instance, the following integrand, singular at both endpoints, can be efficiently evaluated to 100 binary digits:
156 auto f = [](Real x) { return log(x)*log1p(-x); };
157 Real Q = integrator.integrate(f, (Real) 0, (Real) 1);
159 Now onto the caveats: As stated before, the integrands must lie in a Hardy space to ensure rapid convergence.
160 Attempting to integrate a function which is not bounded on the unit disk by tanh-sinh can lead to very slow convergence.
161 For example, take the Runge function:
163 auto f1 = [](double t) { return 1/(1+25*t*t); };
164 Q = integrator.integrate(f1, (double) -1, (double) 1);
166 This function has poles at \u00B1 \u2148/5, and as such it is not bounded on the unit disk.
167 However, the related function
169 auto f2 = [](double t) { return 1/(1+0.04*t*t); };
170 Q = integrator.integrate(f2, (double) -1, (double) 1);
172 has poles outside the unit disk (at \u00B1 5\u2148), and is therefore in the Hardy space.
173 Our benchmarks show that the second integration is performed 22x faster than the first!
174 If you do not understand the structure of your integrand in the complex plane, do performance testing before deployment.
176 Like the trapezoidal quadrature, the tanh-sinh quadrature produces an estimate of the L[sub 1] norm of the integral along with the requested integral.
177 This is to establish a scale against which to measure the tolerance, and to provide an estimate of the condition number of the summation.
178 This can be queried as follows:
180 tanh_sinh<double> integrator;
181 auto f = [](double x) { return 5*x + 7; };
182 double termination = std::sqrt(std::numeric_limits<double>::epsilon());
186 double Q = integrator.integrate(f, 0.0, 1.0, termination, &error, &L1, &levels);
187 double condition_number = L1/std::abs(Q);
189 If the condition number is large, the computed integral is worthless: typically one can assume that Q has lost one digit of precision
190 when the condition number of O(10^Q).
191 The returned error term is not the actual error in the result, but merely an a posteriori error estimate.
192 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.
193 The following table illustrates how the errors and conditioning vary for few sample integrals, in each case the termination condition was set
194 to the square root of epsilon, and all tests were conducted in double precision:
197 [[Integral][Range][Error][Actual measured error][Levels][Condition Number][Comments]]
198 [[`5 * x + 7`][(0,1)][3.5e-15][0][5][1][This trivial case shows just how accurate these methods can be.]]
199 [[`log(x) * log(x)`][0, 1)][0][0][5][1][This is an example of an integral that Gaussian integrators fail to handle.]]
200 [[`exp(-x) / sqrt(x)`][(0,+[infin])][8.0e-10][1.1e-15][5][1][Gaussian integrators typically fail to handle the singularities at the endpoints of this one.]]
201 [[`x * sin(2 * exp(2 * sin(2 * exp(2 * x))))`][(-1,1)][7.2e-16][4.9e-17][9][1.89][This is a truely horrible integral that oscillates wildly and
202 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.]]
203 [[`x == 0 ? 1 : sin(x) / x`][(-[infin], [infin])][3.0e-1][4.0e-1][15][159][This highly oscillatory integral isn't handled at all well by tanh-sinh quadrature: there is so much
204 cancellation in the sum that the result is essentially worthless. The argument transformation of the infinite integral behaves somewhat badly as well, in fact
205 we do ['slightly] better integrating over 2 symmetrical and large finite limits.]]
206 [[`sqrt(x / (1 - x * x))`][(0,1)][1e-8][1e-8][5][1][This an example of an integral that has all its area close to a non-zero endpoint, the problem here is that
207 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:
208 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.]]
209 [[`x < 0.5 ? sqrt(x) / sqrt(1 - x * x) : sqrt(x / ((x + 1) * (xc)))`][(0,1)][0][0][5][1][This is the 2-argument version of the previous integral, the second
210 argument ['xc] is `1-x` in this case, and we use 1-x[super 2] == (1-x)(1+x) to calculate 1-x[super 2] with greater accuracy.]]
213 Although the `tanh-sinh` quadrature can compute integral over infinite domains by variable transformations, these transformations can create a very poorly behaved integrand.
214 For this reason, double-exponential variable transformations have been provided that allow stable computation over infinite domains; these being the exp-sinh and sinh-sinh quadrature.
216 [h4 Complex integrals]
218 The `tanh_sinh` integrator supports integration of functions which return complex results, for example the sine-integral `Si(z)` has the integral representation:
220 [equation sine_integral]
222 Which we can code up directly as:
224 template <class Complex>
225 Complex Si(Complex z)
227 typedef typename Complex::value_type value_type;
228 using std::sin; using std::cos; using std::exp;
229 auto f = [&z](value_type t) { return -exp(-z * cos(t)) * cos(z * sin(t)); };
230 boost::math::quadrature::tanh_sinh<value_type> integrator;
231 return integrator.integrate(f, 0, boost::math::constants::half_pi<value_type>()) + boost::math::constants::half_pi<value_type>();
234 [endsect] [/section:de_tanh_sinh tanh_sinh]
236 [section:de_tanh_sinh_2_arg Handling functions with large features near an endpoint with tanh-sinh quadrature]
238 Tanh-sinh quadrature has a unique feature which makes it well suited to handling integrals with either singularities or large features of interest
239 near one or both endpoints, it turns out that when we calculate and store the abscissa values at which we will be evaluating the function, we can equally
240 well calculate the difference between the abscissa value and the nearest endpoint.
241 This makes it possible to perform quadrature arbitrarily close to an endpoint, without suffering cancellation error.
242 Note however, that we never actually reach the endpoint, so any endpoint singularity will always be excluded from the quadrature.
244 The tanh_sinh integration routine will use this additional information internally when performing range transformation, so that for example,
245 if one end of the range is zero (or infinite) then our transformations will get arbitrarily close to the endpoint without precision loss.
247 However, there are some integrals which may have all of their area near ['both] endpoints, or else near the non-zero endpoint, and in that situation
248 there is a very real risk of loss of precision. For example:
250 tanh_sinh<double> integrator;
251 auto f = [](double x) { return sqrt(x / (1 - x * x); };
252 double Q = integrator.integrate(f, 0.0, 1.0);
254 Results in very low accuracy as all the area of the integral is near 1, and the `1 - x * x` term suffers from cancellation error here.
256 However, both of tanh_sinh's integration routines will automatically handle 2 argument functors: in this case the first argument is the abscissa value as
257 before, while the second is the distance to the nearest endpoint, ie `a - x` or `b - x` if we are integrating over (a,b).
258 You can always differentiate between these 2 cases because the second argument will be negative if we are nearer to the left endpoint.
260 Knowing this, we can rewrite our lambda expression to take advantage of this additional information:
262 tanh_sinh<double> integrator;
263 auto f = [](double x, double xc) { return x <= 0.5 ? sqrt(x) / sqrt(1 - x * x) : sqrt(x / ((x + 1) * (xc))); };
264 double Q = integrator.integrate(f, 0.0, 1.0);
266 Not only is this form accurate to full machine-precision, but it converges to the result faster as well.
268 [endsect] [/section:de_tanh_sinh_2_arg Handling functions with large features near an endpoint with tanh-sinh quadrature]
270 [section:de_sinh_sinh sinh_sinh]
276 sinh_sinh(size_t max_refinements = 9);
279 auto integrate(const F f,
280 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
281 Real* error = nullptr,
283 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;;
286 The sinh-sinh quadrature allows computation over the entire real line, and is called as follows:
288 sinh_sinh<double> integrator;
289 auto f = [](double x) { return exp(-x*x); };
292 double Q = integrator.integrate(f, &error, &L1);
294 Note that the limits of integration are understood to be (-[infin], +[infin]).
296 Complex valued integrands are supported as well, for example the [@https://en.wikipedia.org/wiki/Dirichlet_eta_function Dirichlet Eta function]
297 can be represented via:
299 [equation complex_eta_integral]
301 which we can directly code up as:
303 template <class Complex>
304 Complex eta(Complex s)
306 typedef typename Complex::value_type value_type;
307 using std::pow; using std::exp;
309 value_type pi = boost::math::constants::pi<value_type>();
310 auto f = [&, s, i](value_type t) { return pow(0.5 + i * t, -s) / (exp(pi * t) + exp(-pi * t)); };
311 boost::math::quadrature::sinh_sinh<value_type> integrator;
312 return integrator.integrate(f);
316 [endsect] [/section:de_sinh_sinh sinh_sinh]
318 [section:de_exp_sinh exp_sinh]
324 exp_sinh(size_t max_refinements = 9);
327 auto integrate(const F f, Real a, Real b,
328 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
329 Real* error = nullptr,
331 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
333 auto integrate(const F f,
334 Real tol = sqrt(std::numeric_limits<Real>::epsilon()),
335 Real* error = nullptr,
337 size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
340 For half-infinite intervals, the `exp-sinh` quadrature is provided:
342 exp_sinh<double> integrator;
343 auto f = [](double x) { return exp(-3*x); };
344 double termination = sqrt(std::numeric_limits<double>::epsilon());
347 double Q = integrator.integrate(f, termination, &error, &L1);
350 The native integration range of this integrator is (0, [infin]), but we also support /(a, [infin]), (-[infin], 0)/ and /(-[infin], b)/ via argument transformations.
352 Endpoint singularities and complex-valued integrands are supported by `exp-sinh`.
354 For example, the modified Bessel function K can be represented via:
356 [equation complex_bessel_k_integral]
358 Which we can code up as:
360 template <class Complex>
361 Complex bessel_K(Complex alpha, Complex z)
363 typedef typename Complex::value_type value_type;
364 using std::cosh; using std::exp;
365 auto f = [&, alpha, z](value_type t)
367 value_type ct = cosh(t);
368 if (ct > log(std::numeric_limits<value_type>::max()))
370 return exp(-z * ct) * cosh(alpha * t);
372 boost::math::quadrature::exp_sinh<value_type> integrator;
373 return integrator.integrate(f);
376 The only wrinkle in the above code is the need to check for large `cosh(t)` in which case we assume that
377 `exp(-x cosh(t))` tends to zero faster than `cosh(alpha x)` tends to infinity and return `0`. Without that
378 check we end up with `0 * Infinity` as the result (a NaN).
380 [endsect] [/section:de_exp_sinh exp_sinh]
382 [section:de_tol Setting the Termination Condition for Integration]
384 The integrate method for all three double-exponential quadratures supports ['tolerance] argument that acts as the
385 termination condition for integration.
387 The tolerance is met when two subsequent estimates of the integral have absolute error less than `tolerance*L1`.
389 It is highly recommended that the tolerance be left at the default value of [radic][epsilon], or something similar.
390 Since double exponential quadrature converges exponentially fast for functions in Hardy spaces, then once the routine has *proved* that the error is ~[radic][epsilon],
391 then the error should in fact be ~[epsilon].
393 If you request that the error be ~[epsilon], this tolerance might never be achieved (as the summation is not stabilized ala Kahan),
394 and the routine will simply flounder,
395 dividing the interval in half in order to increase the precision of the integrand, only to be thwarted by floating point roundoff.
397 If for some reason, the default value doesn't quite achieve full precision, then you could try something a little smaller such as
398 [radic][epsilon]/4 or [epsilon][super 2/3].
399 However, more likely, you need to check that your function to be integrated is able to return accurate values, and that there are no other issues with your integration scheme.
401 [endsect] [/section:de_tol Setting the Termination Condition for Integration]
403 [section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
405 The max interval halvings is the maximum number of times the interval can be cut in half before giving up.
406 If the accuracy is not met at that time, the routine returns its best estimate, along with the `error` and `L1`,
407 which allows the user to decide if another quadrature routine should be employed.
409 An example of this is
411 double tol = std::sqrt(std::numeric_limits<double>::epsilon());
412 size_t max_halvings = 12;
413 tanh_sinh<double> integrator(max_halvings);
414 auto f = [](double x) { return 5*x + 7; };
416 double Q = integrator.integrate(f, (double) 0, (double) 1, &error, &L1);
419 Q = some_other_quadrature_method(f, (double) 0, (double) 1);
422 It's important to remember that the number of sample points doubles with each new level, as does the memory footprint
423 of the integrator object. Further, if the integral is smooth, then the precision will be doubling with each new level,
424 so that for example, many integrals can achieve 100 decimal digit precision after just 7 levels. That said, abscissa-weight
425 pairs for new levels are computed only when a new level is actually required (see thread safety), none the less,
426 you should avoid setting the maximum arbitrarily high "just in case" as the time and space requirements for a large
427 number of levels can quickly grow out of control.
429 [endsect] [/section:de_levels Setting the Maximum Interval Halvings and Memory Requirements]
431 [section:de_thread Thread Safety]
433 All three of the double-exponential integrators are thread safe as long as BOOST_MATH_NO_ATOMIC_INT is not set. Since the
434 integrators store a large amount of fairly hard to compute data, it is recommended that these objects are stored and reused
437 Internally all three of the double-exponential integrators use the same caching strategy: they allocate all the vectors needed
438 to store the maximum permitted levels, but only populate the first few levels when constructed. This means a minimal amount of memory
439 is actually allocated when the integrator is first constructed, and already populated levels can be accessed via a lockfree
440 atomic read, and only populating new levels requires a thread lock.
442 In addition, the three built in types (plus `__float128` when available), have the first 7 levels pre-computed: this is generally sufficient for the vast majority
443 of integrals - even at quad precision - and means that integrators for these types are relatively cheap to construct.
445 [endsect] [/section:de_thread Thread Safety]
447 [section:de_caveats Caveats]
449 A few things to keep in mind while using the tanh-sinh, exp-sinh, and sinh-sinh quadratures:
451 These routines are *very* aggressive about approaching the endpoint singularities.
452 This allows lots of significant digits to be extracted, but also has another problem: Roundoff error can cause the function to be evaluated at the endpoints.
453 A few ways to avoid this: Narrow up the bounds of integration to say, [a + [epsilon], b - [epsilon]], make sure (a+b)/2 and (b-a)/2 are representable, and finally,
454 if you think the compromise between accuracy an usability has gone too far in the direction of accuracy, file a ticket.
456 Both exp-sinh and sinh-sinh quadratures evaluate the functions they are passed at *very* large argument.
457 You might understand that x[super 12]exp(-x) is should be zero when x[super 12] overflows, but IEEE floating point arithmetic does not.
458 Hence `std::pow(x, 12)*std::exp(-x)` is an indeterminate form whenever `std::pow(x, 12)` overflows.
459 So make sure your functions have the correct limiting behavior; for example
461 auto f = [](double x) {
470 has the correct behavior for large /x/, but `auto f = [](double x) { return exp(-x)*pow(x, 12); };` does not.
472 Oscillatory integrals, such as the sinc integral, are poorly approximated by double-exponential quadrature.
473 Fortunately the error estimates and L1 norm are massive for these integrals, but nonetheless, oscillatory integrals require different techniques.
475 A special mention should be made about integrating through zero: while our range adaptors preserve precision when one endpoint is zero,
476 things get harder when the origin is neither in the center of the range, nor at an endpoint. Consider integrating:
478 [expression 1 / (1 +x^2)]
480 Over (a, [infin]). As long as `a >= 0` both the tanh_sinh and the exp_sinh integrators will handle this just fine: in fact they provide
481 a rather efficient method for this kind of integral. However, if we have `a < 0` then we are forced to adapt the range in a way that
482 produces abscissa values near zero that have an absolute error of [epsilon], and since all of the area of the integral is near zero
483 both integrators thrash around trying to reach the target accuracy, but never actually get there for `a << 0`. On the other hand, the
484 simple expedient of breaking the integral into two domains: (a, 0) and (0, b) and integrating each seperately using the tanh-sinh
485 integrator, works just fine.
487 Finally, some endpoint singularities are too strong to be handled by `tanh_sinh` or equivalent methods, for example consider integrating
490 double p = some_value;
491 tanh_sinh<double> integrator;
492 auto f = [&](double x){ return pow(tan(x), p); };
493 auto Q = integrator.integrate(f, 0, constants::half_pi<double>());
495 The first problem with this function, is that the singularity is at [pi]/2, so if we're integrating over (0, [pi]/2) then we can never
496 approach closer to the singularity than [epsilon], and for p less than but close to 1, we need to get ['very] close to the singularity
497 to find all the area under the function. If we recall the identity [^tan([pi]/2 - x) == 1/tan(x)] then we can rewrite the function like this:
499 auto f = [&](double x){ return pow(tan(x), -p); };
501 And now the singularity is at the origin and we can get much closer to it when evaluating the integral: all we have done is swap the
502 integral endpoints over.
504 This actually works just fine for p < 0.95, but after that the `tanh_sinh` integrator starts thrashing around and is unable to
505 converge on the integral. The problem is actually a lack of exponent range: if we simply swap type double for something
506 with a greater exponent range (an 80-bit long double or a quad precision type), then we can get to at least p = 0.99. If we want to go
507 beyond that, or stick with type double, then we have to get smart.
509 The easiest method is to notice that for small x, then [^tan(x) [cong] x], and so we are simply integrating x[super -p]. Therefore we can use
510 this approximation over (0, small), and integrate numerically from (small, [pi]/2), using [epsilon] as a suitable crossover point
513 double p = some_value;
514 double crossover = std::numeric_limits<double>::epsilon();
515 tanh_sinh<double> integrator;
516 auto f = [&](double x){ return pow(tan(x), p); };
517 auto Q = integrator.integrate(f, crossover, constants::half_pi<double>()) + pow(crossover, 1 - p) / (1 - p);
519 There is an alternative, more complex method, which is applicable when we are dealing with expressions which can be simplified
520 by evaluating by logs. Let's suppose that as in this case, all the area under the graph is infinitely close to zero, now inagine
521 that we could expand that region out over a much larger range of abscissa values: that's exactly what happens if we perform
522 argument substitution, replacing `x` by `exp(-x)` (note that we must also multiply by the derivative of `exp(-x)`).
523 Now the singularity at zero is moved to +[infin], and the [pi]/2 bound to
524 -log([pi]/2). Initially our argument substituted function looks like:
526 auto f = [&](double x){ return exp(-x) * pow(tan(exp(-x)), -p); };
528 Which is hardly any better, as we still run out of exponent range just as before. However, if we replace `tan(exp(-x))` by `exp(-x)` for suitably
529 small `exp(-x)`, and therefore [^x > -log([epsilon])], we can greatly simplify the expression and evaluate by logs:
531 auto f = [&](double x)
533 static const double crossover = -log(std::numeric_limits<double>::epsilon());
534 return x > crossover ? exp((p - 1) * x) : exp(-x) * pow(tan(exp(-x)), -p);
537 This form integrates just fine over (-log([pi]/2), +[infin]) using either the `tanh_sinh` or `exp_sinh` classes.
539 [endsect] [/section:de_caveats Caveats]
541 [section:de_refes References]
543 * Hidetosi Takahasi and Masatake Mori, ['Double Exponential Formulas for Numerical Integration] Publ. Res. Inst. Math. Sci., 9 (1974), pp. 721-741.
544 * Masatake Mori, ['An IMT-Type Double Exponential Formula for Numerical Integration], Publ RIMS, Kyoto Univ. 14 (1978), 713-729.
545 * David H. Bailey, Karthik Jeyabalan and Xiaoye S. Li ['A Comparison of Three High-Precision Quadrature Schemes] Office of Scientific & Technical Information Technical Reports.
546 * Tanaka, Ken’ichiro, et al. ['Function classes for double exponential integration formulas.] Numerische Mathematik 111.4 (2009): 631-655.
548 [endsect] [/section:de_refes References]
550 [endsect] [/section:double_exponential Double-exponential quadrature]