math/quadrature/exp_sinh.hpp

0001 // Copyright Nick Thompson, 2017
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 /*
0008  * This class performs exp-sinh quadrature on half infinite intervals.
0009  *
0010  * References:
0011  *
0012  * 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655.
0013  */
0014
0015 #ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP
0016 #define BOOST_MATH_QUADRATURE_EXP_SINH_HPP
0017
0018 #include <cmath>
0019 #include <limits>
0020 #include <memory>
0021 #include <string>
0022 #include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
0023
0024 namespace boost{ namespace math{ namespace quadrature {
0025
0026 template<class Real, class Policy = policies::policy<> >
0027 class exp_sinh
0028 {
0029 public:
0030    exp_sinh(size_t max_refinements = 9)
0031       : m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}
0032
0033     template<class F>
0034     auto 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 ->decltype(std::declval<F>()(std::declval<Real>()));
0035     template<class F>
0036     auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));
0037
0038 private:
0039     std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
0040 };
0041
0042 template<class Real, class Policy>
0043 template<class F>
0044 auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
0045 {
0046     typedef decltype(f(a)) K;
0047     static_assert(!std::is_integral<K>::value,
0048                   "The return type cannot be integral, it must be either a real or complex floating point type.");
0049     using std::abs;
0050     using boost::math::constants::half;
0051     using boost::math::quadrature::detail::exp_sinh_detail;
0052
0053     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0054
0055     // Neither limit may be a NaN:
0056     if((boost::math::isnan)(a) || (boost::math::isnan)(b))
0057     {
0058        return static_cast<K>(policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy()));
0059      }
0060     // Right limit is infinite:
0061     if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))
0062     {
0063         // If a = 0, don't use an additional level of indirection:
0064         if (a == static_cast<Real>(0))
0065         {
0066             return m_imp->integrate(f, error, L1, function, tolerance, levels);
0067         }
0068         const auto u = [&](Real t)->K { return f(t + a); };
0069         return m_imp->integrate(u, error, L1, function, tolerance, levels);
0070     }
0071
0072     if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())
0073     {
0074         const auto u = [&](Real t)->K { return f(b-t);};
0075         return m_imp->integrate(u, error, L1, function, tolerance, levels);
0076     }
0077
0078     // Infinite limits:
0079     if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
0080     {
0081         return static_cast<K>(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()));
0082     }
0083     // If we get to here then both ends must necessarily be finite:
0084     return static_cast<K>(policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy()));
0085 }
0086
0087 template<class Real, class Policy>
0088 template<class F>
0089 auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
0090 {
0091     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0092     using std::abs;
0093     if (abs(tolerance) > 1) {
0094         std::string msg = std::string(__FILE__) + ":" + std::to_string(__LINE__) + ":" + std::string(function) + ": The tolerance provided is unusually large; did you confuse it with a domain bound?";
0095         throw std::domain_error(msg);
0096     }
0097     return m_imp->integrate(f, error, L1, function, tolerance, levels);
0098 }
0099
0100
0101 }}}
0102 #endif