import new upstream nautilus stable release 14.2.8

[ceph.git] / ceph / src / boost / boost / math / quadrature / gauss_kronrod.hpp

diff --git a/ceph/src/boost/boost/math/quadrature/gauss_kronrod.hpp b/ceph/src/boost/boost/math/quadrature/gauss_kronrod.hpp
index 2c3a6299c0bce7de0adff496d65bad50c0b50415..fe6ba92c9d06b610a08f9d8ec413cec35370de9c 100644 (file)
--- a/ceph/src/boost/boost/math/quadrature/gauss_kronrod.hpp
+++ b/ceph/src/boost/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
-