Fix tanh_sinh integration for edge case,

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

@@ -229,21 +229,21 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
     // of the abscissa value is greater than the smallest permitted (as specified
     // by the function caller):
     //
-    while (max_left_position && fabs(m_abscissas[0][max_left_position]) < left_min_complement)
+    while ((max_left_position > 1) && fabs(m_abscissas[0][max_left_position]) < left_min_complement)
        --max_left_position;
-    while (max_right_position && fabs(m_abscissas[0][max_right_position]) < right_min_complement)
+    while ((max_right_position > 1) && fabs(m_abscissas[0][max_right_position]) < right_min_complement)
        --max_right_position;
     //
     // Check for non-finite values at the end points:
     // 
     result_type yp, ym, tail_tolerance((std::max)(boost::math::tools::epsilon<Real>(), Real(tolerance * tolerance)));
-    do
+    while (max_left_position)
     {
        yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]);
        if ((boost::math::isfinite)(yp))
           break;
        --max_left_position;
-    } while (m_abscissas[0][max_left_position] < 0);
+    }
     //
     // Also remove points which are insignificant or zero:
     //
@@ -257,13 +257,13 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
     //
     // Over again for the right hand side:
     //
-    do
+    while (max_right_position)
     {
        ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]);
        if ((boost::math::isfinite)(ym))
           break;
        --max_right_position;
-    } while (m_abscissas[0][max_right_position] < 0);
+    }
     while (max_right_position > 1)
     {
        if (abs(ym * m_weights[0][max_right_position]) > abs(L1_I0 * tail_tolerance))
@@ -272,7 +272,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
        ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]);
     }
 
-    if ((max_left_position == 0) && (max_right_position == 0))
+    if ((max_left_position == 0) || (max_right_position == 0))
     {
        return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
     }
@@ -349,8 +349,10 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
         // Thus, we filter which abscissa values generate a call to f(x_i), with a single
         // floating point comparison per loop.  Everything else is integer logic.
         //
+        BOOST_MATH_ASSERT(max_left_position);
         max_left_index = max_left_position - 1;
         max_left_position *= 2;
+        BOOST_MATH_ASSERT(max_right_position);
         max_right_index = max_right_position - 1;
         max_right_position *= 2;
         if ((abscissa_row.size() > max_left_index + 1) && (fabs(abscissa_row[max_left_index + 1]) > left_min_complement))
@@ -371,6 +373,10 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
            yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]);
            if ((boost::math::isfinite)(yp))
               break;
+           if(max_left_position <= 2)
+           {
+              return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
+           }
            max_left_position -= 2;
            --max_left_index;
         } while (abscissa_row[max_left_index] < 0);
@@ -382,6 +388,10 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
            ym = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]);
            if ((boost::math::isfinite)(ym))
               break;
+           if (max_right_position <= 2)
+           {
+              return policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", ym, Policy());
+           }
            --max_right_index;
            max_right_position -= 2;
         } while (abscissa_row[max_right_index] < 0);
@@ -499,7 +509,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
             // trajectory at our new end point, and increase our error estimate by the last
             // good value as an estimate for what we may have discarded.
             //
-            if ((max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement))
+            if (max_left_index && (max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement))
             {
                yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]) * weight_row[max_left_index];
                ym = f(-1 - abscissa_row[max_left_index - 1], abscissa_row[max_left_index - 1]) * weight_row[max_left_index - 1];
@@ -508,7 +518,7 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
                   return policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Integration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint.  Please check your function, and consider providing a 2-argument functor.", I1, Policy());
                }
             }
-            if ((max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement))
+            if (max_right_index && (max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement))
             {
                yp = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]) * weight_row[max_right_index];
                ym = f(1 + abscissa_row[max_right_index - 1], -abscissa_row[max_right_index - 1]) * weight_row[max_right_index - 1];
@@ -521,9 +531,9 @@ decltype(std::declval<F>()(std::declval<Real>(), std::declval<Real>())) tanh_sin
             break;
         }
 
-        if (truncate_left)
+        if ((truncate_left) && (max_left_position > 1))
            --max_left_position;
-        if (truncate_right)
+        if ((truncate_right) && (max_right_position > 1))
            --max_right_position;
 
     }

test/Jamfile.v2

@@ -1461,7 +1461,7 @@ test-suite quadrature :
 
    [ compile compile_test/gauss_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax ] [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : <build>no ] ]
    [ compile compile_test/gauss_kronrod_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax ] [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : <build>no ] ]
-
+   [ run git_issue_898.cpp ]
 ;
 
 test-suite autodiff :

test/git_issue_898.cpp

@@ -0,0 +1,47 @@
+// Copyright Matt Borland, 2022
+// 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)
+
+#include <cmath>
+#include <limits>
+#include <boost/math/quadrature/tanh_sinh.hpp>
+#include <boost/math/constants/constants.hpp>
+#include "math_unit_test.hpp"
+
+// numerically evaluate the integral for Stefan-Boltzmann Law from Planck's Law
+// https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law#Derivation_from_Planck's_law
+template<typename Real>
+Real StefanBoltzmann() 
+{
+   constexpr auto h = Real(6.62607015e-34);
+   constexpr auto c = Real(299792458.0);
+   constexpr auto kB = Real(1.380649e-23);
+
+   constexpr auto c1 = Real(2) * boost::math::constants::pi<Real>() * h * c * c;
+   constexpr auto c2 = h * c / kB;
+   constexpr auto T = 1000;
+
+   auto integrand = [&](const Real l)
+   {
+      return c1 / (std::pow(l, Real(5)) * (std::exp(c2 / (T * l)) - Real(1)));
+   };
+
+   auto integrator = boost::math::quadrature::tanh_sinh<Real>(15);
+   return static_cast<Real>(integrator.integrate(integrand, Real(0.0), std::numeric_limits<Real>::infinity()));
+}
+
+int main() {
+
+   constexpr auto sigma = 56703.74419;
+
+   // Double Precision
+   CHECK_MOLLIFIED_CLOSE(sigma, StefanBoltzmann<double>(), 1e-5);
+
+   // integate using single precision. It gets the right answer, but there is integer underflow
+   // if this line is commented out, there is no underflow
+   CHECK_MOLLIFIED_CLOSE(sigma, StefanBoltzmann<float>(), 1e-5f);
+
+   return boost::math::test::report_errors();
+}