summaryrefslogtreecommitdiff
path: root/boost/math/quadrature/gauss_kronrod.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/quadrature/gauss_kronrod.hpp')
-rw-r--r--boost/math/quadrature/gauss_kronrod.hpp61
1 files changed, 35 insertions, 26 deletions
diff --git a/boost/math/quadrature/gauss_kronrod.hpp b/boost/math/quadrature/gauss_kronrod.hpp
index 2c3a6299c0..fe6ba92c9d 100644
--- a/boost/math/quadrature/gauss_kronrod.hpp
+++ b/boost/math/quadrature/gauss_kronrod.hpp
@@ -1769,18 +1769,19 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
{
typedef detail::gauss_kronrod_detail<Real, N, detail::gauss_constant_category<Real>::value> base;
public:
- typedef Real value_type;
+ typedef Real value_type;
private:
template <class F>
- static value_type integrate_non_adaptive_m1_1(F f, Real* error = nullptr, Real* pL1 = nullptr)
+ static auto integrate_non_adaptive_m1_1(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
{
- using std::fabs;
+ typedef decltype(f(Real(0))) K;
+ using std::abs;
unsigned gauss_start = 2;
unsigned kronrod_start = 1;
unsigned gauss_order = (N - 1) / 2;
- value_type kronrod_result = 0;
- value_type gauss_result = 0;
- value_type fp, fm;
+ K kronrod_result = 0;
+ K gauss_result = 0;
+ K fp, fm;
if (gauss_order & 1)
{
fp = f(value_type(0));
@@ -1794,13 +1795,13 @@ private:
gauss_start = 1;
kronrod_start = 2;
}
- value_type L1 = fabs(kronrod_result);
+ Real L1 = abs(kronrod_result);
for (unsigned i = gauss_start; i < base::abscissa().size(); i += 2)
{
fp = f(base::abscissa()[i]);
fm = f(-base::abscissa()[i]);
kronrod_result += (fp + fm) * base::weights()[i];
- L1 += (fabs(fp) + fabs(fm)) * base::weights()[i];
+ L1 += (abs(fp) + abs(fm)) * base::weights()[i];
gauss_result += (fp + fm) * gauss<Real, (N - 1) / 2>::weights()[i / 2];
}
for (unsigned i = kronrod_start; i < base::abscissa().size(); i += 2)
@@ -1808,12 +1809,12 @@ private:
fp = f(base::abscissa()[i]);
fm = f(-base::abscissa()[i]);
kronrod_result += (fp + fm) * base::weights()[i];
- L1 += (fabs(fp) + fabs(fm)) * base::weights()[i];
+ L1 += (abs(fp) + abs(fm)) * base::weights()[i];
}
if (pL1)
*pL1 = L1;
if (error)
- *error = (std::max)(static_cast<Real>(fabs(kronrod_result - gauss_result)), static_cast<Real>(fabs(kronrod_result * tools::epsilon<Real>() * 2)));
+ *error = (std::max)(static_cast<Real>(abs(kronrod_result - gauss_result)), static_cast<Real>(abs(kronrod_result * tools::epsilon<Real>() * Real(2))));
return kronrod_result;
}
@@ -1825,19 +1826,22 @@ private:
};
template <class F>
- static value_type recursive_adaptive_integrate(const recursive_info<F>* info, Real a, Real b, unsigned max_levels, Real abs_tol, Real* error, Real* L1)
+ static auto recursive_adaptive_integrate(const recursive_info<F>* info, Real a, Real b, unsigned max_levels, Real abs_tol, Real* error, Real* L1)->decltype(std::declval<F>()(std::declval<Real>()))
{
- using std::fabs;
+ typedef decltype(info->f(Real(a))) K;
+ using std::abs;
Real error_local;
Real mean = (b + a) / 2;
Real scale = (b - a) / 2;
- auto ff = [&](const Real& x)->Real
+ auto ff = [&](const Real& x)->K
{
return info->f(scale * x + mean);
};
- Real estimate = scale * integrate_non_adaptive_m1_1(ff, &error_local, L1);
+ K r1 = integrate_non_adaptive_m1_1(ff, &error_local, L1);
+ K estimate = scale * r1;
- Real abs_tol1 = fabs(estimate * info->tol);
+ K tmp = estimate * info->tol;
+ Real abs_tol1 = abs(tmp);
if (abs_tol == 0)
abs_tol = abs_tol1;
@@ -1862,35 +1866,41 @@ private:
public:
template <class F>
- static value_type integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr)
+ static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
{
+ typedef decltype(f(a)) K;
static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)";
if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
{
// Infinite limits:
if ((a <= -tools::max_value<Real>()) && (b >= tools::max_value<Real>()))
{
- auto u = [&](const Real& t)->Real
+ auto u = [&](const Real& t)->K
{
Real t_sq = t*t;
Real inv = 1 / (1 - t_sq);
- return f(t*inv)*(1 + t_sq)*inv*inv;
+ Real w = (1 + t_sq)*inv*inv;
+ Real arg = t*inv;
+ K res = f(arg)*w;
+ return res;
};
recursive_info<decltype(u)> info = { u, tol };
- return recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+ K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+ return res;
}
// Right limit is infinite:
if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
{
- auto u = [&](const Real& t)->Real
+ auto u = [&](const Real& t)->K
{
Real z = 1 / (t + 1);
Real arg = 2 * z + a - 1;
- return f(arg)*z*z;
+ K res = f(arg)*z*z;
+ return res;
};
recursive_info<decltype(u)> info = { u, tol };
- Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+ K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
if (pL1)
{
*pL1 *= 2;
@@ -1900,14 +1910,14 @@ public:
if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
{
- auto v = [&](const Real& t)->Real
+ auto v = [&](const Real& t)->K
{
Real z = 1 / (t + 1);
Real arg = 2 * z - 1;
return f(b - arg) * z * z;
};
recursive_info<decltype(v)> info = { v, tol };
- Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+ K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
if (pL1)
{
*pL1 *= 2;
@@ -1925,7 +1935,7 @@ public:
return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1);
}
}
- return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy());
+ return static_cast<K>(policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()));
}
};
@@ -1938,4 +1948,3 @@ public:
#endif
#endif // BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP
-