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 <boost/math/tools/config.hpp>
0019 #include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
0020
0021 #ifndef BOOST_MATH_HAS_NVRTC
0022
0023 #include <cmath>
0024 #include <limits>
0025 #include <memory>
0026 #include <string>
0027
0028 namespace boost{ namespace math{ namespace quadrature {
0029
0030 template<class Real, class Policy = policies::policy<> >
0031 class exp_sinh
0032 {
0033 public:
0034    exp_sinh(size_t max_refinements = 9)
0035       : m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}
0036
0037     template<class F>
0038     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>()));
0039     template<class F>
0040     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>()));
0041
0042 private:
0043     std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
0044 };
0045
0046 template<class Real, class Policy>
0047 template<class F>
0048 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>()))
0049 {
0050     typedef decltype(f(a)) K;
0051     static_assert(!std::is_integral<K>::value,
0052                   "The return type cannot be integral, it must be either a real or complex floating point type.");
0053     using std::abs;
0054     using boost::math::constants::half;
0055     using boost::math::quadrature::detail::exp_sinh_detail;
0056
0057     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0058
0059     // Neither limit may be a NaN:
0060     if((boost::math::isnan)(a) || (boost::math::isnan)(b))
0061     {
0062        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()));
0063      }
0064     // Right limit is infinite:
0065     if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))
0066     {
0067         // If a = 0, don't use an additional level of indirection:
0068         if (a == static_cast<Real>(0))
0069         {
0070             return m_imp->integrate(f, error, L1, function, tolerance, levels);
0071         }
0072         const auto u = [&](Real t)->K { return f(t + a); };
0073         return m_imp->integrate(u, error, L1, function, tolerance, levels);
0074     }
0075
0076     if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())
0077     {
0078         const auto u = [&](Real t)->K { return f(b-t);};
0079         return m_imp->integrate(u, error, L1, function, tolerance, levels);
0080     }
0081
0082     // Infinite limits:
0083     if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
0084     {
0085         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()));
0086     }
0087     // If we get to here then both ends must necessarily be finite:
0088     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()));
0089 }
0090
0091 template<class Real, class Policy>
0092 template<class F>
0093 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>()))
0094 {
0095     static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0096     using std::abs;
0097     if (abs(tolerance) > 1) {
0098         return policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy());
0099     }
0100     return m_imp->integrate(f, error, L1, function, tolerance, levels);
0101 }
0102
0103
0104 }}}
0105
0106 #endif // BOOST_MATH_HAS_NVRTC
0107
0108 #ifdef BOOST_MATH_ENABLE_CUDA
0109
0110 #include <boost/math/tools/type_traits.hpp>
0111 #include <boost/math/tools/cstdint.hpp>
0112 #include <boost/math/tools/precision.hpp>
0113 #include <boost/math/policies/error_handling.hpp>
0114 #include <boost/math/constants/constants.hpp>
0115
0116 namespace boost {
0117 namespace math {
0118 namespace quadrature {
0119
0120 template <class F, class Real, class Policy = policies::policy<> >
0121 __device__ auto exp_sinh_integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, boost::math::size_t* levels)
0122 {
0123     BOOST_MATH_STD_USING
0124
0125     using K = decltype(f(a));
0126     static_assert(!boost::math::is_integral<K>::value,
0127                   "The return type cannot be integral, it must be either a real or complex floating point type.");
0128
0129     constexpr auto function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0130
0131     // Neither limit may be a NaN:
0132     if((boost::math::isnan)(a) || (boost::math::isnan)(b))
0133     {
0134        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()));
0135     }
0136     // Right limit is infinite:
0137     if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value<Real>()))
0138     {
0139         // If a = 0, don't use an additional level of indirection:
0140         if (a == static_cast<Real>(0))
0141         {
0142             return detail::exp_sinh_integrate_impl(f, tolerance, error, L1, levels);
0143         }
0144         const auto u = [&](Real t)->K { return f(t + a); };
0145         return detail::exp_sinh_integrate_impl(u, tolerance, error, L1, levels);
0146     }
0147
0148     if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value<Real>())
0149     {
0150         const auto u = [&](Real t)->K { return f(b-t);};
0151         return detail::exp_sinh_integrate_impl(u, tolerance, error, L1, levels);
0152     }
0153
0154     // Infinite limits:
0155     if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
0156     {
0157         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()));
0158     }
0159     // If we get to here then both ends must necessarily be finite:
0160     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()));
0161 }
0162
0163 template <class F, class Real, class Policy = policies::policy<> >
0164 __device__ auto exp_sinh_integrate(const F& f, Real tolerance, Real* error, Real* L1, boost::math::size_t* levels)
0165 {
0166     BOOST_MATH_STD_USING
0167     constexpr auto function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
0168     if (abs(tolerance) > 1) {
0169         return policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy());
0170     }
0171     return detail::exp_sinh_integrate_impl(f, tolerance, error, L1, levels);
0172 }
0173
0174 } // namespace quadrature
0175 } // namespace math
0176 } // namespace boost
0177
0178 #endif // BOOST_MATH_ENABLE_CUDA
0179
0180 #endif // BOOST_MATH_QUADRATURE_EXP_SINH_HPP