/usr/include/boost/math/quadrature
Edit: /usr/include/boost/math/quadrature/exp_sinh.hpp (3934B)
// Copyright Nick Thompson, 2017
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
/*
* This class performs exp-sinh quadrature on half infinite intervals.
*
* References:
*
* 1) Tanaka, Ken'ichiro, et al. "Function classes for double exponential integration formulas." Numerische Mathematik 111.4 (2009): 631-655.
*/
#ifndef BOOST_MATH_QUADRATURE_EXP_SINH_HPP
#define BOOST_MATH_QUADRATURE_EXP_SINH_HPP
#include
#include
#include
#include
namespace boost{ namespace math{ namespace quadrature {
template >
class exp_sinh
{
public:
exp_sinh(size_t max_refinements = 9)
: m_imp(std::make_shared>(max_refinements)) {}
template
auto integrate(const F& f, Real a, Real b, Real tol = boost::math::tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval()(std::declval())) const;
template
auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval()(std::declval())) const;
private:
std::shared_ptr> m_imp;
};
template
template
auto exp_sinh::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval()(std::declval())) const
{
typedef decltype(f(a)) K;
using std::abs;
using boost::math::constants::half;
using boost::math::quadrature::detail::exp_sinh_detail;
static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
// Neither limit may be a NaN:
if((boost::math::isnan)(a) || (boost::math::isnan)(b))
{
return static_cast(policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy()));
}
// Right limit is infinite:
if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value()))
{
// If a = 0, don't use an additional level of indirection:
if (a == (Real) 0)
{
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}
const auto u = [&](Real t)->K { return f(t + a); };
return m_imp->integrate(u, error, L1, function, tolerance, levels);
}
if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value())
{
const auto u = [&](Real t)->K { return f(b-t);};
return m_imp->integrate(u, error, L1, function, tolerance, levels);
}
// Infinite limits:
if ((a <= -boost::math::tools::max_value()) && (b >= boost::math::tools::max_value()))
{
return static_cast(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()));
}
// If we get to here then both ends must necessarily be finite:
return static_cast(policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy()));
}
template
template
auto exp_sinh::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval()(std::declval())) const
{
static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}
}}}
#endif