290 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 tanh-sinh quadrature on the real line.9 * Tanh-sinh quadrature is exponentially convergent for integrands in Hardy spaces,10 * (see https://en.wikipedia.org/wiki/Hardy_space for a formal definition), and is optimal for a random function from that class.11 *12 * The tanh-sinh quadrature is one of a class of so called "double exponential quadratures"-there is a large family of them,13 * but this one seems to be the most commonly used.14 *15 * As always, there are caveats: For instance, if the function you want to integrate is not holomorphic on the unit disk,16 * then the rapid convergence will be spoiled. In this case, a more appropriate quadrature is (say) Romberg, which does not17 * require the function to be holomorphic, only differentiable up to some order.18 *19 * In addition, if you are integrating a periodic function over a period, the trapezoidal rule is better.20 *21 * References:22 *23 * 1) Mori, Masatake. "Quadrature formulas obtained by variable transformation and the DE-rule." Journal of Computational and Applied Mathematics 12 (1985): 119-130.24 * 2) Bailey, David H., Karthik Jeyabalan, and Xiaoye S. Li. "A comparison of three high-precision quadrature schemes." Experimental Mathematics 14.3 (2005): 317-329.25 * 3) Press, William H., et al. "Numerical recipes third edition: the art of scientific computing." Cambridge University Press 32 (2007): 10013-2473.26 *27 */28 29#ifndef BOOST_MATH_QUADRATURE_TANH_SINH_HPP30#define BOOST_MATH_QUADRATURE_TANH_SINH_HPP31 32#include <cmath>33#include <limits>34#include <memory>35#include <boost/math/quadrature/detail/tanh_sinh_detail.hpp>36 37namespace boost{ namespace math{ namespace quadrature {38 39template<class Real, class Policy = policies::policy<> >40class tanh_sinh41{42public:43 tanh_sinh(size_t max_refinements = 15, const Real& min_complement = tools::min_value<Real>() * 4)44 : m_imp(std::make_shared<detail::tanh_sinh_detail<Real, Policy>>(max_refinements, min_complement)) {}45 46 template<class F>47 auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));48 template<class F>49 auto integrate(const F f, Real a, Real b, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));50 51 template<class F>52 auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>()));53 template<class F>54 auto integrate(const F f, Real tolerance = tools::root_epsilon<Real>(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>()));55 56private:57 std::shared_ptr<detail::tanh_sinh_detail<Real, Policy>> m_imp;58};59 60template<class Real, class Policy>61template<class F>62auto tanh_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>()))63{64 BOOST_MATH_STD_USING65 using boost::math::constants::half;66 using boost::math::quadrature::detail::tanh_sinh_detail;67 68 static const char* function = "tanh_sinh<%1%>::integrate";69 70 typedef decltype(std::declval<F>()(std::declval<Real>())) result_type;71 static_assert(!std::is_integral<result_type>::value,72 "The return type cannot be integral, it must be either a real or complex floating point type.");73 if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))74 {75 76 // Infinite limits:77 if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))78 {79 auto u = [&](const Real& t, const Real& tc)->result_type80 {81 Real t_sq = t*t;82 Real inv;83 if (t > 0.5f)84 inv = 1 / ((2 - tc) * tc);85 else if(t < -0.5)86 inv = 1 / ((2 + tc) * -tc);87 else88 inv = 1 / (1 - t_sq);89 return f(t*inv)*(1 + t_sq)*inv*inv;90 };91 Real limit = sqrt(tools::min_value<Real>()) * 4;92 return m_imp->integrate(u, error, L1, function, limit, limit, tolerance, levels);93 }94 95 // Right limit is infinite:96 if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))97 {98 auto u = [&](const Real& t, const Real& tc)->result_type99 {100 Real z, arg;101 if (t > -0.5f)102 z = 1 / (t + 1);103 else104 z = -1 / tc;105 if (t < 0.5)106 arg = 2 * z + a - 1;107 else108 arg = a + tc / (2 - tc);109 return f(arg)*z*z;110 };111 Real left_limit = sqrt(tools::min_value<Real>()) * 4;112 result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);113 if (L1)114 {115 *L1 *= 2;116 }117 if (error)118 {119 *error *= 2;120 }121 122 return Q;123 }124 125 if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))126 {127 auto v = [&](const Real& t, const Real& tc)->result_type128 {129 Real z;130 if (t > -0.5)131 z = 1 / (t + 1);132 else133 z = -1 / tc;134 Real arg;135 if (t < 0.5)136 arg = 2 * z - 1;137 else138 arg = tc / (2 - tc);139 return f(b - arg) * z * z;140 };141 142 Real left_limit = sqrt(tools::min_value<Real>()) * 4;143 result_type Q = Real(2) * m_imp->integrate(v, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);144 if (L1)145 {146 *L1 *= 2;147 }148 if (error)149 {150 *error *= 2;151 }152 return Q;153 }154 155 if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))156 {157 if (a == b)158 {159 return result_type(0);160 }161 if (b < a)162 {163 return -this->integrate(f, b, a, tolerance, error, L1, levels);164 }165 Real avg = (a + b)*half<Real>();166 Real diff = (b - a)*half<Real>();167 Real avg_over_diff_m1 = a / diff;168 Real avg_over_diff_p1 = b / diff;169 bool have_small_left = fabs(a) < 0.5f;170 bool have_small_right = fabs(b) < 0.5f;171 Real left_min_complement = float_next(avg_over_diff_m1) - avg_over_diff_m1;172 Real min_complement_limit = (std::max)(tools::min_value<Real>(), float_next(Real(tools::min_value<Real>() / diff)));173 if (left_min_complement < min_complement_limit)174 left_min_complement = min_complement_limit;175 Real right_min_complement = avg_over_diff_p1 - float_prior(avg_over_diff_p1);176 if (right_min_complement < min_complement_limit)177 right_min_complement = min_complement_limit;178 //179 // These asserts will fail only if rounding errors on180 // type Real have accumulated so much error that it's181 // broken our internal logic. Should that prove to be182 // a persistent issue, we might need to add a bit of fudge183 // factor to move left_min_complement and right_min_complement184 // further from the end points of the range.185 //186 BOOST_MATH_ASSERT((left_min_complement * diff + a) > a);187 BOOST_MATH_ASSERT((b - right_min_complement * diff) < b);188 auto u = [&](Real z, Real zc)->result_type189 {190 Real position;191 if (z < -0.5)192 {193 if(have_small_left)194 return f(diff * (avg_over_diff_m1 - zc));195 position = a - diff * zc;196 }197 else if (z > 0.5)198 {199 if(have_small_right)200 return f(diff * (avg_over_diff_p1 - zc));201 position = b - diff * zc;202 }203 else204 position = avg + diff*z;205 BOOST_MATH_ASSERT(position != a);206 BOOST_MATH_ASSERT(position != b);207 return f(position);208 };209 result_type Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);210 211 if (L1)212 {213 *L1 *= diff;214 }215 if (error)216 {217 *error *= diff;218 }219 return Q;220 }221 }222 return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());223}224 225template<class Real, class Policy>226template<class F>227auto tanh_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>(), std::declval<Real>()))228{229 BOOST_MATH_STD_USING230 using boost::math::constants::half;231 using boost::math::quadrature::detail::tanh_sinh_detail;232 233 static const char* function = "tanh_sinh<%1%>::integrate";234 235 if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b))236 {237 if (b <= a)238 {239 return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy());240 }241 auto u = [&](Real z, Real zc)->Real242 {243 if (z < 0)244 return f((a - b) * zc / 2 + a, (b - a) * zc / 2);245 else246 return f((a - b) * zc / 2 + b, (b - a) * zc / 2);247 };248 Real diff = (b - a)*half<Real>();249 Real left_min_complement = tools::min_value<Real>() * 4;250 Real right_min_complement = tools::min_value<Real>() * 4;251 Real Q = diff*m_imp->integrate(u, error, L1, function, left_min_complement, right_min_complement, tolerance, levels);252 253 if (L1)254 {255 *L1 *= diff;256 }257 if (error)258 {259 *error *= diff;260 }261 return Q;262 }263 return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());264}265 266template<class Real, class Policy>267template<class F>268auto tanh_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>()))269{270 using boost::math::quadrature::detail::tanh_sinh_detail;271 static const char* function = "tanh_sinh<%1%>::integrate";272 Real min_complement = tools::epsilon<Real>();273 return m_imp->integrate([&](const Real& arg, const Real&) { return f(arg); }, error, L1, function, min_complement, min_complement, tolerance, levels);274}275 276template<class Real, class Policy>277template<class F>278auto tanh_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>(), std::declval<Real>()))279{280 using boost::math::quadrature::detail::tanh_sinh_detail;281 static const char* function = "tanh_sinh<%1%>::integrate";282 Real min_complement = tools::min_value<Real>() * 4;283 return m_imp->integrate(f, error, L1, function, min_complement, min_complement, tolerance, levels);284}285 286}287}288}289#endif290