WSJT-X/boost/math/quadrature/exp_sinh.hpp
Bill Somerville d361e123c6 Squashed 'boost/' changes from d9443bc48..c27aa31f0
c27aa31f0 Updated Boost to v1.70.0 including iterator range math numeric crc circular_buffer multi_index intrusive

git-subtree-dir: boost
git-subtree-split: c27aa31f06ebf1a91b3fa3ae9df9b5efdf14ec9f
2019-07-02 23:38:24 +01:00

95 lines
3.8 KiB
C++

// 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 <cmath>
#include <limits>
#include <memory>
#include <boost/math/quadrature/detail/exp_sinh_detail.hpp>
namespace boost{ namespace math{ namespace quadrature {
template<class Real, class Policy = policies::policy<> >
class exp_sinh
{
public:
exp_sinh(size_t max_refinements = 9)
: m_imp(std::make_shared<detail::exp_sinh_detail<Real, Policy>>(max_refinements)) {}
template<class F>
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)->decltype(std::declval<F>()(std::declval<Real>())) const;
template<class F>
auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr)->decltype(std::declval<F>()(std::declval<Real>())) const;
private:
std::shared_ptr<detail::exp_sinh_detail<Real, Policy>> m_imp;
};
template<class Real, class Policy>
template<class F>
auto exp_sinh<Real, Policy>::integrate(const F& f, Real a, Real b, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) 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<K>(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<Real>()))
{
// 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)->Real { 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<Real>())
{
const auto u = [&](Real t)->Real { return f(b-t);};
return m_imp->integrate(u, error, L1, function, tolerance, levels);
}
// Infinite limits:
if ((a <= -boost::math::tools::max_value<Real>()) && (b >= boost::math::tools::max_value<Real>()))
{
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()));
}
// If we get to here then both ends must necessarily be finite:
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()));
}
template<class Real, class Policy>
template<class F>
auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error, Real* L1, std::size_t* levels)->decltype(std::declval<F>()(std::declval<Real>())) const
{
static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}
}}}
#endif