ceph/src/boost/boost/math/quadrature/exp_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 exp-sinh quadrature on half infinite intervals.
   9  *
  10  * References:
  11  *
  12  * 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655.
  13  */
  14
  15 #ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP
  16 #define BOOST_MATH_QUADRATURE_EXP_SINH_HPP
  17
  18 #include <cmath>
  19 #include <limits>
  20 #include <memory>
  21 #include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
  22
  23 namespace boost{ namespace math{ namespace quadrature {
  24
  25 template<class Real, class Policy = policies::policy<> >
  26 class exp_sinh
  27 {
  28 public:
  29    exp_sinh(size_t max_refinements = 9)
  30       : m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}
  31
  32     template<class F>
  33     Real integrate(const F& f, Real a, Real b, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const;
  34     template<class F>
  35     Real integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const;
  36
  37 private:
  38     std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
  39 };
  40
  41 template<class Real, class Policy>
  42 template<class F>
  43 Real exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const
  44 {
  45     using std::abs;
  46     using boost::math::constants::half;
  47     using boost::math::quadrature::detail::exp_sinh_detail;
  48
  49     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
  50
  51     // Neither limit may be a NaN:
  52     if((boost::math::isnan)(a) || (boost::math::isnan)(b))
  53        return policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy());
  54     // Right limit is infinite:
  55     if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))
  56     {
  57         // If a = 0, don't use an additional level of indirection:
  58         if (a == (Real) 0)
  59         {
  60             return m_imp->integrate(f, error, L1, function, tolerance, levels);
  61         }
  62         const auto u = [&](Real t)->Real { return f(t + a); };
  63         return m_imp->integrate(u, error, L1, function, tolerance, levels);
  64     }
  65
  66     if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())
  67     {
  68         const auto u = [&](Real t)->Real { return f(b-t);};
  69         return m_imp->integrate(u, error, L1, function, tolerance, levels);
  70     }
  71
  72     // Infinite limits:
  73     if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
  74     {
  75         return policies::raise_domain_error(function, "Use sinh_sinh quadrature for integration over the whole real line; exp_sinh is for half infinite integrals.", a, Policy());
  76     }
  77     // If we get to here then both ends must necessarily be finite:
  78     return policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy());
  79 }
  80
  81 template<class Real, class Policy>
  82 template<class F>
  83 Real exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const
  84 {
  85     using std::abs;
  86     using boost::math::constants::half;
  87     using boost::math::quadrature::detail::exp_sinh_detail;
  88
  89     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
  90
  91     return m_imp->integrate(f, error, L1, function, tolerance, levels);
  92 }
  93
  94
  95 }}}
  96 #endif