ceph/src/boost/boost/math/quadrature/tanh_sinh.hpp

   1 // Copyright Nick Thompson, 2017
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 /*
   8  * This class performs tanh-sinh quadrature on the real line.
   9  * Tanh-sinh quadrature is exponentially convergent for integrands in Hardy spaces,
  10  * (see https://en.wikipedia.org/wiki/Hardy_space for a formal definition), and is optimal for a random function from that class.
  11  *
  12  * The tanh-sinh quadrature is one of a class of so called "double exponential quadratures"-there is a large family of them,
  13  * but this one seems to be the most commonly used.
  14  *
  15  * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk,
  16  * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not
  17  * require the function to be holomorphic, only differentiable up to some order.
  18  *
  19  * In addition, if you are integrating a periodic function over a period, the trapezoidal rule is better.
  20  *
  21  * References:
  22  *
  23  * 1) Mori, Masatake. "Quadrature formulas obtained by variable transformation and the DE-rule." Journal of Computational and Applied Mathematics 12 (1985): 119-130.
  24  * 2) Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. "A comparison of three high-precision quadrature schemes." Experimental Mathematics 14.3 (2005): 317-329.
  25  * 3) Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.
  26  *
  27  */
  28
  29 #ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP
  30 #define BOOST_MATH_QUADRATURE_TANH_SINH_HPP
  31
  32 #include <cmath>
  33 #include <limits>
  34 #include <memory>
  35 #include <boost/math/quadrature/detail/tanh_sinh_detail.hpp>
  36
  37 namespace boost{ namespace math{ namespace quadrature {
  38
  39 template<class Real, class Policy = policies::policy<> >
  40 class tanh_sinh
  41 {
  42 public:
  43     tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)
  44     : m_imp(std::make_shared<detail::tanh_sinh_detail<Real, Policy>>(max_refinements, min_complement)) {}
  45
  46     template<class F>
  47     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;
  48     template<class F>
  49     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>(), std::declval<Real>())) const;
  50
  51     template<class F>
  52     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;
  53     template<class F>
  54     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>(), std::declval<Real>())) const;
  55
  56 private:
  57     std::shared_ptr<detail::tanh_sinh_detail<Real, Policy>> m_imp;
  58 };
  59
  60 template<class Real, class Policy>
  61 template<class F>
  62 auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>())) const
  63 {
  64     BOOST_MATH_STD_USING
  65     using boost::math::constants::half;
  66     using boost::math::quadrature::detail::tanh_sinh_detail;
  67
  68     static const char* function = "tanh_sinh<%1%>::integrate";
  69
  70     typedef decltype(std::declval<F>()(std::declval<Real>())) result_type;
  71     static_assert(!std::is_integral<result_type>::value,
  72                   "The return type cannot be integral, it must be either a real or complex floating point type.");
  73     if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
  74     {
  75
  76        // Infinite limits:
  77        if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
  78        {
  79           auto u = [&](const Real& t, const Real& tc)->result_type
  80           {
  81              Real t_sq = t*t;
  82              Real inv;
  83              if (t > 0.5f)
  84                 inv = 1 / ((2 - tc) * tc);
  85              else if(t < -0.5)
  86                 inv = 1 / ((2 + tc) * -tc);
  87              else
  88                 inv = 1 / (1 - t_sq);
  89              return f(t*inv)*(1 + t_sq)*inv*inv;
  90           };
  91           Real limit = sqrt(tools::min_value<Real>()) * 4;
  92           return m_imp->integrate(u, error, L1, function, limit, limit, tolerance, levels);
  93        }
  94
  95        // Right limit is infinite:
  96        if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
  97        {
  98           auto u = [&](const Real& t, const Real& tc)->result_type
  99           {
 100              Real z, arg;
 101              if (t > -0.5f)
 102                 z = 1 / (t + 1);
 103              else
 104                 z = -1 / tc;
 105              if (t < 0.5)
 106                 arg = 2 * z + a - 1;
 107              else
 108                 arg = a + tc / (2 - tc);
 109              return f(arg)*z*z;
 110           };
 111           Real left_limit = sqrt(tools::min_value<Real>()) * 4;
 112           result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
 113           if (L1)
 114           {
 115              *L1 *= 2;
 116           }
 117           if (error)
 118           {
 119              *error *= 2;
 120           }
 121
 122           return Q;
 123        }
 124
 125        if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
 126        {
 127           auto v = [&](const Real& t, const Real& tc)->result_type
 128           {
 129              Real z;
 130              if (t > -0.5)
 131                 z = 1 / (t + 1);
 132              else
 133                 z = -1 / tc;
 134              Real arg;
 135              if (t < 0.5)
 136                 arg = 2 * z - 1;
 137              else
 138                 arg = tc / (2 - tc);
 139              return f(b - arg) * z * z;
 140           };
 141
 142           Real left_limit = sqrt(tools::min_value<Real>()) * 4;
 143           result_type Q = Real(2) * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
 144           if (L1)
 145           {
 146              *L1 *= 2;
 147           }
 148           if (error)
 149           {
 150              *error *= 2;
 151           }
 152           return Q;
 153        }
 154
 155        if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
 156        {
 157           if (a == b)
 158           {
 159              return result_type(0);
 160           }
 161           if (b < a)
 162           {
 163              return -this->integrate(f, b, a, tolerance, error, L1, levels);
 164           }
 165           Real avg = (a + b)*half<Real>();
 166           Real diff = (b - a)*half<Real>();
 167           Real avg_over_diff_m1 = a / diff;
 168           Real avg_over_diff_p1 = b / diff;
 169           bool have_small_left = fabs(a) < 0.5f;
 170           bool have_small_right = fabs(b) < 0.5f;
 171           Real left_min_complement = float_next(avg_over_diff_m1) - avg_over_diff_m1;
 172           Real min_complement_limit = (std::max)(tools::min_value<Real>(), Real(tools::min_value<Real>() / diff));
 173           if (left_min_complement < min_complement_limit)
 174              left_min_complement = min_complement_limit;
 175           Real right_min_complement = avg_over_diff_p1 - float_prior(avg_over_diff_p1);
 176           if (right_min_complement < min_complement_limit)
 177              right_min_complement = min_complement_limit;
 178           //
 179           // These asserts will fail only if rounding errors on
 180           // type Real have accumulated so much error that it's
 181           // broken our internal logic.  Should that prove to be
 182           // a persistent issue, we might need to add a bit of fudge
 183           // factor to move left_min_complement and right_min_complement
 184           // further from the end points of the range.
 185           //
 186           BOOST_MATH_ASSERT((left_min_complement * diff + a) > a);
 187           BOOST_MATH_ASSERT((b - right_min_complement * diff) < b);
 188           auto u = [&](Real z, Real zc)->result_type
 189           {
 190              Real position;
 191              if (z < -0.5)
 192              {
 193                 if(have_small_left)
 194                   return f(diff * (avg_over_diff_m1 - zc));
 195                 position = a - diff * zc;
 196              }
 197              else if (z > 0.5)
 198              {
 199                 if(have_small_right)
 200                   return f(diff * (avg_over_diff_p1 - zc));
 201                 position = b - diff * zc;
 202              }
 203              else
 204                 position = avg + diff*z;
 205              BOOST_MATH_ASSERT(position != a);
 206              BOOST_MATH_ASSERT(position != b);
 207              return f(position);
 208           };
 209           result_type Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
 210
 211           if (L1)
 212           {
 213              *L1 *= diff;
 214           }
 215           if (error)
 216           {
 217              *error *= diff;
 218           }
 219           return Q;
 220        }
 221     }
 222     return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
 223 }
 224
 225 template<class Real, class Policy>
 226 template<class F>
 227 auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const
 228 {
 229    BOOST_MATH_STD_USING
 230       using boost::math::constants::half;
 231    using boost::math::quadrature::detail::tanh_sinh_detail;
 232
 233    static const char* function = "tanh_sinh<%1%>::integrate";
 234
 235    if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))
 236    {
 237       if (b <= a)
 238       {
 239          return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());
 240       }
 241       auto u = [&](Real z, Real zc)->Real
 242       {
 243          if (z < 0)
 244             return f((a - b) * zc / 2 + a, (b - a) * zc / 2);
 245          else
 246             return f((a - b) * zc / 2 + b, (b - a) * zc / 2);
 247       };
 248       Real diff = (b - a)*half<Real>();
 249       Real left_min_complement = tools::min_value<Real>() * 4;
 250       Real right_min_complement = tools::min_value<Real>() * 4;
 251       Real Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);
 252
 253       if (L1)
 254       {
 255          *L1 *= diff;
 256       }
 257       if (error)
 258       {
 259          *error *= diff;
 260       }
 261       return Q;
 262    }
 263    return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
 264 }
 265
 266 template<class Real, class Policy>
 267 template<class F>
 268 auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>())) const
 269 {
 270    using boost::math::quadrature::detail::tanh_sinh_detail;
 271    static const char* function = "tanh_sinh<%1%>::integrate";
 272    Real min_complement = tools::epsilon<Real>();
 273    return m_imp->integrate([&](const Real& arg, const Real&) { return f(arg); }, error, L1, function, min_complement, min_complement, tolerance, levels);
 274 }
 275
 276 template<class Real, class Policy>
 277 template<class F>
 278 auto tanh_sinh<Real, Policy>::integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) const
 279 {
 280    using boost::math::quadrature::detail::tanh_sinh_detail;
 281    static const char* function = "tanh_sinh<%1%>::integrate";
 282    Real min_complement = tools::min_value<Real>() * 4;
 283    return m_imp->integrate(f, error, L1, function, min_complement, min_complement, tolerance, levels);
 284 }
 285
 286 }
 287 }
 288 }
 289 #endif