181 lines · plain
1// Copyright Nick Thompson, 20172// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0.4// (See accompanying file LICENSE_1_0.txt5// 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_HPP16#define BOOST_MATH_QUADRATURE_EXP_SINH_HPP17 18#include <boost/math/tools/config.hpp>19#include <boost/math/quadrature/detail/exp_sinh_detail.hpp>20 21#ifndef BOOST_MATH_HAS_NVRTC22 23#include <cmath>24#include <limits>25#include <memory>26#include <string>27 28namespace boost{ namespace math{ namespace quadrature {29 30template<class Real, class Policy = policies::policy<> >31class exp_sinh32{33public:34 exp_sinh(size_t max_refinements = 9)35 : m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}36 37 template<class F>38 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>()));39 template<class F>40 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>()));41 42private:43 std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;44};45 46template<class Real, class Policy>47template<class F>48auto 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>()))49{50 typedef decltype(f(a)) K;51 static_assert(!std::is_integral<K>::value,52 "The return type cannot be integral, it must be either a real or complex floating point type.");53 using std::abs;54 using boost::math::constants::half;55 using boost::math::quadrature::detail::exp_sinh_detail;56 57 static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";58 59 // Neither limit may be a NaN:60 if((boost::math::isnan)(a) || (boost::math::isnan)(b))61 {62 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()));63 }64 // Right limit is infinite:65 if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))66 {67 // If a = 0, don't use an additional level of indirection:68 if (a == static_cast<Real>(0))69 {70 return m_imp->integrate(f, error, L1, function, tolerance, levels);71 }72 const auto u = [&](Real t)->K { return f(t + a); };73 return m_imp->integrate(u, error, L1, function, tolerance, levels);74 }75 76 if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())77 {78 const auto u = [&](Real t)->K { return f(b-t);};79 return m_imp->integrate(u, error, L1, function, tolerance, levels);80 }81 82 // Infinite limits:83 if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))84 {85 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()));86 }87 // If we get to here then both ends must necessarily be finite:88 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()));89}90 91template<class Real, class Policy>92template<class F>93auto 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>()))94{95 static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";96 using std::abs;97 if (abs(tolerance) > 1) {98 return policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy());99 }100 return m_imp->integrate(f, error, L1, function, tolerance, levels);101}102 103 104}}}105 106#endif // BOOST_MATH_HAS_NVRTC107 108#ifdef BOOST_MATH_ENABLE_CUDA109 110#include <boost/math/tools/type_traits.hpp>111#include <boost/math/tools/cstdint.hpp>112#include <boost/math/tools/precision.hpp>113#include <boost/math/policies/error_handling.hpp>114#include <boost/math/constants/constants.hpp>115 116namespace boost { 117namespace math { 118namespace quadrature {119 120template <class F, class Real, class Policy = policies::policy<> >121__device__ auto exp_sinh_integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, boost::math::size_t* levels)122{123 BOOST_MATH_STD_USING124 125 using K = decltype(f(a));126 static_assert(!boost::math::is_integral<K>::value,127 "The return type cannot be integral, it must be either a real or complex floating point type.");128 129 constexpr auto function = "boost::math::quadrature::exp_sinh<%1%>::integrate";130 131 // Neither limit may be a NaN:132 if((boost::math::isnan)(a) || (boost::math::isnan)(b))133 {134 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()));135 }136 // Right limit is infinite:137 if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))138 {139 // If a = 0, don't use an additional level of indirection:140 if (a == static_cast<Real>(0))141 {142 return detail::exp_sinh_integrate_impl(f, tolerance, error, L1, levels);143 }144 const auto u = [&](Real t)->K { return f(t + a); };145 return detail::exp_sinh_integrate_impl(u, tolerance, error, L1, levels);146 }147 148 if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())149 {150 const auto u = [&](Real t)->K { return f(b-t);};151 return detail::exp_sinh_integrate_impl(u, tolerance, error, L1, levels);152 }153 154 // Infinite limits:155 if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))156 {157 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()));158 }159 // If we get to here then both ends must necessarily be finite:160 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()));161}162 163template <class F, class Real, class Policy = policies::policy<> >164__device__ auto exp_sinh_integrate(const F& f, Real tolerance, Real* error, Real* L1, boost::math::size_t* levels)165{166 BOOST_MATH_STD_USING167 constexpr auto function = "boost::math::quadrature::exp_sinh<%1%>::integrate";168 if (abs(tolerance) > 1) {169 return policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy());170 }171 return detail::exp_sinh_integrate_impl(f, tolerance, error, L1, levels);172}173 174} // namespace quadrature175} // namespace math176} // namespace boost177 178#endif // BOOST_MATH_ENABLE_CUDA179 180#endif // BOOST_MATH_QUADRATURE_EXP_SINH_HPP181