Disallow integral return types in quadrature. (#704)

include/boost/math/quadrature/detail/sinh_sinh_detail.hpp

@@ -158,6 +158,8 @@ auto sinh_sinh_detail<Real, Policy>::integrate(const F f, Real tolerance, Real*
     static const char* function = "boost::math::quadrature::sinh_sinh<%1%>::integrate";
 
     typedef decltype(f(Real(0))) K;
+    static_assert(!std::is_integral<K>::value,
+                  "The return type cannot be integral, it must be either a real or complex floating point type.");
     K y_max = f(boost::math::tools::max_value<Real>());
     if(abs(y_max) > boost::math::tools::epsilon<Real>())
     {

include/boost/math/quadrature/exp_sinh.hpp

@@ -44,6 +44,8 @@ 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;
+    static_assert(!std::is_integral<K>::value,
+                  "The return type cannot be integral, it must be either a real or complex floating point type.");
     using std::abs;
     using boost::math::constants::half;
     using boost::math::quadrature::detail::exp_sinh_detail;

include/boost/math/quadrature/gauss.hpp

@@ -1179,6 +1179,8 @@ class gauss : public detail::gauss_detail<Real, N, detail::gauss_constant_catego
      // In many math texts, K represents the field of real or complex numbers.
      // Too bad we can't put blackboard bold into C++ source!
       typedef decltype(f(Real(0))) K;
+      static_assert(!std::is_integral<K>::value,
+                   "The return type cannot be integral, it must be either a real or complex floating point type.");
       using std::abs;
       unsigned non_zero_start = 1;
       K result = Real(0);

include/boost/math/quadrature/gauss_kronrod.hpp

@@ -1871,6 +1871,8 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
    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_assert(!std::is_integral<K>::value,
+                  "The return type cannot be integral, it must be either a real or complex floating point type.");
       static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)";
       if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b))
       {

include/boost/math/quadrature/tanh_sinh.hpp

@@ -68,24 +68,25 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
     static const char* function = "tanh_sinh<%1%>::integrate";
 
     typedef decltype(std::declval<F>()(std::declval<Real>())) result_type;
-
+    static_assert(!std::is_integral<result_type>::value,
+                  "The return type cannot be integral, it must be either a real or complex floating point type.");
     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, const Real& tc)->result_type
-          { 
-             Real t_sq = t*t; 
+          {
+             Real t_sq = t*t;
              Real inv;
              if (t > 0.5f)
                 inv = 1 / ((2 - tc) * tc);
              else if(t < -0.5)
                 inv = 1 / ((2 + tc) * -tc);
              else
                 inv = 1 / (1 - t_sq);
-             return f(t*inv)*(1 + t_sq)*inv*inv; 
+             return f(t*inv)*(1 + t_sq)*inv*inv;
           };
           Real limit = sqrt(tools::min_value<Real>()) * 4;
           return m_imp->integrate(u, error, L1, function, limit, limit, tolerance, levels);
@@ -95,7 +96,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
        if ((boost::math::isfinite)(a) && (b >= tools::max_value<Real>()))
        {
           auto u = [&](const Real& t, const Real& tc)->result_type
-          { 
+          {
              Real z, arg;
              if (t > -0.5f)
                 z = 1 / (t + 1);
@@ -105,7 +106,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
                 arg = 2 * z + a - 1;
              else
                 arg = a + tc / (2 - tc);
-             return f(arg)*z*z; 
+             return f(arg)*z*z;
           };
           Real left_limit = sqrt(tools::min_value<Real>()) * 4;
           result_type Q = Real(2) * m_imp->integrate(u, error, L1, function, left_limit, tools::min_value<Real>(), tolerance, levels);
@@ -124,7 +125,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
        if ((boost::math::isfinite)(b) && (a <= -tools::max_value<Real>()))
        {
           auto v = [&](const Real& t, const Real& tc)->result_type
-          { 
+          {
              Real z;
              if (t > -0.5)
                 z = 1 / (t + 1);
@@ -185,7 +186,7 @@ auto tanh_sinh<Real, Policy>::integrate(const F f, Real a, Real b, Real toleranc
           BOOST_MATH_ASSERT((left_min_complement * diff + a) > a);
           BOOST_MATH_ASSERT((b - right_min_complement * diff) < b);
           auto u = [&](Real z, Real zc)->result_type
-          { 
+          {
              Real position;
              if (z < -0.5)
              {

include/boost/math/quadrature/trapezoidal.hpp

@@ -35,6 +35,8 @@ auto trapezoidal(F f, Real a, Real b, Real tol, std::size_t max_refinements, Rea
     // In many math texts, K represents the field of real or complex numbers.
     // Too bad we can't put blackboard bold into C++ source!
     typedef decltype(f(a)) K;
+    static_assert(!std::is_integral<K>::value,
+                  "The return type cannot be integral, it must be either a real or complex floating point type.");
     if (!(boost::math::isfinite)(a))
     {
        return static_cast<K>(boost::math::policies::raise_domain_error(function, "Left endpoint of integration must be finite for adaptive trapezoidal integration but got a = %1%.\n", a, pol));