Double exponential quadrature

I'm trying to lighten the code review load for the maintainers of boost.math, and I was hoping you guys could help me out. I have a pull request which implements tanh-sinh quadrature, which is provably optimal for integrands in Hardy spaces.

Here's my code:

• Documentation

• Interface

• Implementation

• Tests

which is also reproduced below.

I have a few design worries.

• It is a class and not a function. This is a bit confusing; I worry that people will not recognize that the constructor is doing some one-time calculations to make integrations faster.

• It takes a long time to compile. I generated the abscissas and weights to 100 digits, and then they must be lexically cast to their correct type. I could keep fewer levels of abscissas and weights in the .hpp, but then the runtime would longer for complex integrands.

• For integrands in Hardy spaces, the number of correct digits always doubles on each refinement. However, we want to do just the right amount of work to deliver the requested accuracy, which is almost always overshot.

Interface:

```
#ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP
#define BOOST_MATH_QUADRATURE_TANH_SINH_HPP

#include
#include
#include
#include

namespace boost{ namespace math{

template
class tanh_sinh
{
public:
tanh_sinh(Real tol = sqrt(std::numeric_limits::epsilon()), size_t max_refinements = 20);

template
Real integrate(F f, Real a, Real b, Real* error = nullptr);

private:
std::shared_ptr> m_imp;
};

template
tanh_sinh::tanh_sinh(Real tol, size_t max_refinements) : m_imp(std::make_shared>(tol, max_refinements))
{
return;
}

template
template
Real tanh_sinh::integrate(F f, Real a, Real b, Real* error)
{
using std::isfinite;
using boost::math::constants::half;
using boost::math::detail::tanh_sinh_detail;

if (isfinite(a) && isfinite(b))
{
if (b a.\n");
}
Real avg = (a+b)*half();
Real diff = (b-a)*half();
auto u = = { return f(avg + d

Other than possible code nits (like "use member initalization lists", there are some suggestions I could make:

Calculate tables only once

For the cost of creating the object and the storage requirements, you could have static functions to generate your constant tables.

template 
const std::vector>& get_abcissas() {
  static const std::vector> abcissas() = {
    {
        lexical_cast("0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"), // g(0)
        lexical_cast("0.951367964072746945727055362904639667492765811307212865380079106867050650113429723656597692697291568999499"), // g(1)
        lexical_cast("0.999977477192461592864899636308688849285982470957489530009950811164291603926066803141755297920571692976244"), // g(2)
        lexical_cast("0.999999999999957058389441217592223051901253805502353310119906858210171731776098247943305316472169355401371"), // g(3)
        lexical_cast("0.999999999999999999999999999999999999883235110249013906725663510362671044720752808415603434458608954302982"), // g(4)
        lexical_cast("0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999988520"), // g(5)
    },
    {
        lexical_cast("0.674271492248435826080420090632052144267244477154740136576044169121421222924677035505101849981126759183585"), // g(0.5)
        lexical_cast("0.997514856457224386832717419238820368149231215381809295391585321457677448277585202664213469426402672227688"), // g(1.5)
        lexical_cast("0.999999988875664881984668015033322737014902900245095922058323073481945768209599289672119775131473502267885"), // g(2.5)
        lexical_cast("0.999999999999999999999946215084086063112254432391666747077319911504923816981286361121293026490456057993796"), // g(3.5)
        lexical_cast("0.999999999999999999999999999999999999999999999999999999999999920569786807778838966034206923747918174840316"), // g(4.5)
    },
  };
  return abcissas;
};

At least under C++11, this solution would guarantee that the vector would only be created once during the execution of the program.

Also, for efficiency reasons (math code is something that is bound to be in inner loops), you could consider:

• Keeping all the coefficients in one array, instead of an array of arrays. Locality is always good.

• Since you only have static members, you can use an std::array instead of an std::vector

Calculate parameters in a constexpr

You could go radical and calculate the other members you need (m_t_max) in a constexpr function so that you don't need to calculate them at runtime.

With this, you could theoretically make the number of refinements a template parameter (if it can be known at build time) and you wouldn't need to calculate m_t_max all the time.

Maybe you don't need a functor after all

With the above suggestion, you suddenly don't have any more state for your functor. It would be possible to have just a standalone function with Real and n_refinements template parameters.

StackExchange Code Review Q#158137, answer score: 2