Send exact -infinity in log functions to the overflow policy

include/boost/math/distributions/extreme_value.hpp

@@ -190,15 +190,16 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const extreme_value_distribution<R
 
    constexpr auto function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
 
-   if((BOOST_MATH_NAMESPACE::isinf)(x))
-      return x < 0 ? 0.0f : 1.0f;
    RealType a = dist.location();
    RealType b = dist.scale();
    RealType result = 0;
    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
       return result;
    if(0 == detail::check_finite(function, a, &result, Policy()))
       return result;
+   // log(0) at the infinite end is an exact infinity, which goes to the overflow policy:
+   if((BOOST_MATH_NAMESPACE::isinf)(x))
+      return x < 0 ? -policies::raise_overflow_error<RealType>(function, nullptr, Policy()) : RealType(0);
    if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
       return result;
 
@@ -265,15 +266,16 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<extreme_v
 
    constexpr auto function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
 
-   if((BOOST_MATH_NAMESPACE::isinf)(c.param))
-      return c.param < 0 ? 1.0f : 0.0f;
    RealType a = c.dist.location();
    RealType b = c.dist.scale();
    RealType result = 0;
    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
       return result;
    if(0 == detail::check_finite(function, a, &result, Policy()))
       return result;
+   // log(0) at the infinite end is an exact infinity, which goes to the overflow policy:
+   if((BOOST_MATH_NAMESPACE::isinf)(c.param))
+      return c.param < 0 ? RealType(0) : -policies::raise_overflow_error<RealType>(function, nullptr, Policy());
    if(0 == detail::check_x(function, c.param, &result, Policy()))
       return result;
 

include/boost/math/distributions/laplace.hpp

@@ -246,14 +246,10 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const laplace_distribution<RealTyp
       return result;
    }
 
-   // Special cdf values:
+   // The log of the cdf at -infinity is an exact -infinity, which goes to the overflow policy:
    if((BOOST_MATH_NAMESPACE::isinf)(x))
    {
-      if(x < 0) 
-      {
-         return 0; // -infinity.
-      }
-      return 1; // + infinity.
+      return x < 0 ? -policies::raise_overflow_error<RealType>(function, nullptr, Policy()) : RealType(0);
    }
 
    if (false == detail::check_x(function, x, &result, Policy())) 
@@ -369,15 +365,10 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<laplace_d
    // Check scale and location.
     if (false == c.dist.check_parameters(function, &result)) return result;
 
-   // Special cdf values.
+   // The log of the complement at +infinity is an exact -infinity, which goes to the overflow policy:
    if((BOOST_MATH_NAMESPACE::isinf)(x))
    {
-     if(x < 0) 
-     { 
-       return 1; // cdf complement -infinity is unity.
-     }
-
-     return 0; // cdf complement +infinity is zero.
+     return x < 0 ? RealType(0) : -policies::raise_overflow_error<RealType>(function, nullptr, Policy());
    }
    if(false == detail::check_x(function, x, &result, Policy()))return result;
 

include/boost/math/distributions/logistic.hpp

@@ -168,13 +168,10 @@ BOOST_MATH_NAMESPACE_BEGIN
           return result;
        }
 
+       // The log of the cdf at -infinity is an exact -infinity, which goes to the overflow policy:
        if((BOOST_MATH_NAMESPACE::isinf)(x))
        {
-          if(x < 0) 
-          {
-            return 0; // -infinity
-          }
-          return 1; // + infinity
+          return x < 0 ? -policies::raise_overflow_error<RealType>(function, nullptr, Policy()) : RealType(0);
        }
 
        if(false == detail::check_x(function, x, &result, Policy()))
@@ -273,10 +270,10 @@ BOOST_MATH_NAMESPACE_BEGIN
        {
           return result;
        }
+       // The log of the complement at +infinity is an exact -infinity, which goes to the overflow policy:
        if((BOOST_MATH_NAMESPACE::isinf)(x))
        {
-          if(x < 0) return 1; // cdf complement -infinity is unity.
-          return 0; // cdf complement +infinity is zero.
+          return x < 0 ? RealType(0) : -policies::raise_overflow_error<RealType>(function, nullptr, Policy());
        }
        if(false == detail::check_x(function, x, &result, Policy()))
        {

include/boost/math/distributions/students_t.hpp

@@ -325,8 +325,13 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const students_t_distribution<Real
    RealType log_tail;
    if (ax > sqrt(tools::max_value<RealType>()))
    {
-      // x^2 would overflow; df / (df + x^2) is then (df / |x|) / |x|.
-      log_tail = BOOST_MATH_NAMESPACE::libeta(df / 2, static_cast<RealType>(0.5), RealType((df / ax) / ax), Policy());
+      // x^2 would overflow; df / (df + x^2) is then z = (df / |x|) / |x|. Where z is below epsilon (it may
+      // underflow), log I_z(a, 1/2) = a log(z) - log(a) - lbeta(a, 1/2) + O(z), with log(z) taken from the logs:
+      RealType z = (df / ax) / ax;
+      RealType a = df / 2;
+      log_tail = z < tools::epsilon<RealType>()
+         ? RealType(a * (log(df) - 2 * log(ax)) - log(a) - BOOST_MATH_NAMESPACE::lbeta(a, static_cast<RealType>(0.5), Policy()))
+         : RealType(BOOST_MATH_NAMESPACE::libeta(a, static_cast<RealType>(0.5), z, Policy()));
    }
    else
    {

include/boost/math/distributions/weibull.hpp

@@ -181,9 +181,10 @@ BOOST_MATH_GPU_ENABLED inline RealType logpdf(const weibull_distribution<RealTyp
       {
          return log(1 / scale);
       }
+      // pdf(0) is zero, so its log is an exact -infinity, which goes to the overflow policy:
       if(shape > 1)
       {
-         return 0;
+         return -policies::raise_overflow_error<RealType>(function, nullptr, Policy());
       }
       return policies::raise_overflow_error<RealType>(function, 0, Policy());
    }

include/boost/math/special_functions/beta.hpp

@@ -2003,7 +2003,7 @@ BOOST_MATH_GPU_ENABLED T log_ibeta_power_terms(T a, T b, T x, T y, const BOOST_M
 // Logarithm of the regularised incomplete beta, or of its complement when invert is true.
 //
 template <class T, class Policy>
-BOOST_MATH_GPU_ENABLED T libeta_imp(T a, T b, T x, const Policy& pol, bool invert, bool* exact_infinity = nullptr)
+BOOST_MATH_GPU_ENABLED T libeta_imp(T a, T b, T x, const Policy& pol, bool invert)
 {
    BOOST_MATH_STD_USING
    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
@@ -2030,17 +2030,8 @@ BOOST_MATH_GPU_ENABLED T libeta_imp(T a, T b, T x, const Policy& pol, bool inver
    bool target_is_zero = invert ? q_is_zero : p_is_zero;
    if(target_is_zero)
    {
-      // -infinity is the exact answer here, not an overflow; the caller returns it unchecked.
-      BOOST_MATH_IF_CONSTEXPR(BOOST_MATH_NAMESPACE::numeric_limits<T>::has_infinity)
-      {
-         if(exact_infinity)
-            *exact_infinity = true;
-         return -BOOST_MATH_NAMESPACE::numeric_limits<T>::infinity();
-      }
-      else
-      {
-         return -policies::raise_overflow_error<T>(function, nullptr, pol);
-      }
+      // -infinity is exact here; as for every other exact infinity, the overflow policy decides:
+      return -policies::raise_overflow_error<T>(function, nullptr, pol);
    }
 
    // The underflow is ours to handle, so stop ibeta_imp from reporting it:
@@ -2200,12 +2191,6 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   // log(beta(a, b)) is exactly -infinity when either (valid) argument is infinite:
-   BOOST_MATH_IF_CONSTEXPR(BOOST_MATH_NAMESPACE::numeric_limits<value_type>::has_infinity)
-   {
-      if((a > 0) && (b > 0) && ((BOOST_MATH_NAMESPACE::isinf)(a) || (BOOST_MATH_NAMESPACE::isinf)(b)))
-         return -BOOST_MATH_NAMESPACE::numeric_limits<result_type>::infinity();
-   }
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lbeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::lbeta<%1%>(%1%,%1%)");
 }
 
@@ -2310,10 +2295,7 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<RT1, RT2, RT3>
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   bool exact_infinity = false;
-   value_type result = detail::libeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, &exact_infinity);
-   if(exact_infinity)
-      return static_cast<result_type>(result);
+   value_type result = detail::libeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false);
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::libeta<%1%>(%1%,%1%,%1%)");
 }
 
@@ -2338,10 +2320,7 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<RT1, RT2, RT3>
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   bool exact_infinity = false;
-   value_type result = detail::libeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, &exact_infinity);
-   if(exact_infinity)
-      return static_cast<result_type>(result);
+   value_type result = detail::libeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true);
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::libetac<%1%>(%1%,%1%,%1%)");
 }
 

include/boost/math/special_functions/gamma.hpp

@@ -1847,6 +1847,9 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_imp(T a, T x, const Policy& pol)
       return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
    if(x < 0)
       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+   // The log of an exact zero; the overflow policy decides, as for ligamma:
+   if((BOOST_MATH_NAMESPACE::isinf)(x))
+      return -policies::raise_overflow_error<T>(function, nullptr, pol);
 
    if (((x > 1000) || (x > tools::log_max_value<T>() - 10)) && (a + 50 < x))
    {
@@ -1913,6 +1916,9 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol
       return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
    if(x < 0)
       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+   // The log of an exact zero; the overflow policy decides, as for ligamma:
+   if(x == 0)
+      return -policies::raise_overflow_error<T>(function, nullptr, pol);
 
    // Taken from conditions on Line 1709. There are also
    // conditions on Line 1368, but didn't implement that one here. 
@@ -1969,7 +1975,7 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol
 // (ACM TOMS 46(1), Algorithm 1006), the prefix x^a e^-x is never formed.
 //
 template <class T, class Policy>
-BOOST_MATH_GPU_ENABLED T ligamma_imp(T a, T x, const Policy& pol, bool upper, bool* exact_infinity = nullptr)
+BOOST_MATH_GPU_ENABLED T ligamma_imp(T a, T x, const Policy& pol, bool upper)
 {
    BOOST_MATH_STD_USING
    const char* function = upper ? "boost::math::ligamma<%1%>(%1%, %1%)" : "boost::math::ligamma_lower<%1%>(%1%, %1%)";
@@ -1986,17 +1992,8 @@ BOOST_MATH_GPU_ENABLED T ligamma_imp(T a, T x, const Policy& pol, bool upper, bo
    bool target_is_zero = upper ? (BOOST_MATH_NAMESPACE::isinf)(x) : (x == 0);
    if(target_is_zero)
    {
-      // -infinity is the exact answer here, not an overflow; the caller returns it unchecked.
-      BOOST_MATH_IF_CONSTEXPR(BOOST_MATH_NAMESPACE::numeric_limits<T>::has_infinity)
-      {
-         if(exact_infinity)
-            *exact_infinity = true;
-         return -BOOST_MATH_NAMESPACE::numeric_limits<T>::infinity();
-      }
-      else
-      {
-         return -policies::raise_overflow_error<T>(function, nullptr, pol);
-      }
+      // -infinity is exact here; as for every other exact infinity, the overflow policy decides:
+      return -policies::raise_overflow_error<T>(function, nullptr, pol);
    }
    if(upper ? (x == 0) : (BOOST_MATH_NAMESPACE::isinf)(x))
       return BOOST_MATH_NAMESPACE::lgamma(a, pol);
@@ -2728,10 +2725,7 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma(T1 a, T2 z,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   bool exact_infinity = false;
-   value_type result = detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), true, &exact_infinity);
-   if(exact_infinity)
-      return static_cast<result_type>(result);
+   value_type result = detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), true);
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::ligamma<%1%>(%1%, %1%)");
 }
 
@@ -2754,10 +2748,7 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma_lower(T1 a,
       policies::discrete_quantile<>,
       policies::assert_undefined<> >::type forwarding_policy;
 
-   bool exact_infinity = false;
-   value_type result = detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), false, &exact_infinity);
-   if(exact_infinity)
-      return static_cast<result_type>(result);
+   value_type result = detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), false);
    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::ligamma_lower<%1%>(%1%, %1%)");
 }
 

include/boost/math/tools/polynomial_roots.hpp

@@ -9,6 +9,7 @@
 #include <boost/math/tools/config.hpp>
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/tools/fft.hpp>
+#include <boost/math/tools/precision.hpp>
 #include <boost/math/tools/detail/complex_from_real.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 #include <boost/math/special_functions/next.hpp>
@@ -251,7 +252,7 @@ detail::real_from_complex_t<Complex> scaled_residual(std::vector<Coefficient> co
       f = f * x + c;
       scale = scale * ax + abs(c);
    }
-   return abs(f) / ((std::max)(scale, (std::numeric_limits<Real>::min)()));
+   return abs(f) / ((std::max)(scale, tools::min_value<Real>()));
 }
 
 template <class Complex>
@@ -318,7 +319,7 @@ std::vector<Complex> aberth(std::vector<Complex> a,
             Complex(cos(angle), sin(angle));
       }
    }
-   Real eps = std::numeric_limits<Real>::epsilon();
+   Real eps = tools::epsilon<Real>();
    // A root whose last correction was negligible is frozen: it still repels the others but is no longer
    // updated. Each sweep then costs O(n * active) instead of O(n^2).
    const bool global = !(BOOST_MATH_NAMESPACE::isfinite)(interest);
@@ -423,7 +424,7 @@ Real inclusion_radius(std::vector<Real> const& a, Complex z)
    using std::abs;
    using complex_type = Complex;
    std::size_t n = a.size() - 1;
-   Real eps = std::numeric_limits<Real>::epsilon();
+   Real eps = tools::epsilon<Real>();
    bool reversed = abs(z) > 1;
    complex_type x = reversed ? complex_type(Real(1) / z) : z;
    Real ax = abs(x);
@@ -453,7 +454,7 @@ namespace polynomial_roots_detail {
 template <class Real>
 std::vector<detail::complex_from_real_t<Real> > piecewise(std::vector<Real> coefficients)
 {
-   const unsigned precision = std::numeric_limits<Real>::digits;
+   const unsigned precision = static_cast<unsigned>(tools::digits<Real>());
    static_assert(!std::is_integral<Real>::value, "polynomial_roots requires a floating-point coefficient type");
    using std::abs;
    using std::log2;
@@ -485,8 +486,8 @@ std::vector<detail::complex_from_real_t<Real> > piecewise(std::vector<Real> coef
    struct candidate { complex_type z; Real radius; Real residual; };
    std::vector<candidate> candidates;
    std::size_t taylor_degree = (std::min)(degree, std::size_t(4) * precision);
-   Real log_min = log2((std::numeric_limits<Real>::min)());
-   Real log_max = log2((std::numeric_limits<Real>::max)());
+   Real log_min = log2(tools::min_value<Real>());
+   Real log_max = log2(tools::max_value<Real>());
 
    for(auto const& piece : pieces)
    {
@@ -553,7 +554,7 @@ std::vector<detail::complex_from_real_t<Real> > piecewise(std::vector<Real> coef
             coefficient_norm += abs(coefficient);
          Real discarded = 0;
          std::size_t effective_size = local[k].size();
-         Real discard_limit = Real(8) * std::numeric_limits<Real>::epsilon() * coefficient_norm;
+         Real discard_limit = Real(8) * tools::epsilon<Real>() * coefficient_norm;
          while(effective_size > 1 &&
                discarded + abs(local[k][effective_size - 1]) <= discard_limit)
          {
@@ -575,16 +576,16 @@ std::vector<detail::complex_from_real_t<Real> > piecewise(std::vector<Real> coef
                local_value = local_value * t + local[k][q];
                local_scale = local_scale * abs(t) + abs(local[k][q]);
             }
-            local_scale = (std::max)(local_scale, (std::numeric_limits<Real>::min)());
+            local_scale = (std::max)(local_scale, tools::min_value<Real>());
             if(!(BOOST_MATH_NAMESPACE::isfinite)(local_value.real()) || !(BOOST_MATH_NAMESPACE::isfinite)(local_value.imag()) ||
-               abs(local_value) > sqrt(std::numeric_limits<Real>::epsilon()) * local_scale)
+               abs(local_value) > sqrt(tools::epsilon<Real>()) * local_scale)
                continue;
             complex_type z = phase * (gamma + rho * t);
             Real modulus = abs(z);
-            if(modulus >= inner * (Real(1) - Real(32) * std::numeric_limits<Real>::epsilon()) &&
-               modulus <= outer * (Real(1) + Real(32) * std::numeric_limits<Real>::epsilon()))
+            if(modulus >= inner * (Real(1) - Real(32) * tools::epsilon<Real>()) &&
+               modulus <= outer * (Real(1) + Real(32) * tools::epsilon<Real>()))
             {
-               Real error = abs(local_value) + Real(8) * std::numeric_limits<Real>::epsilon() * local_scale;
+               Real error = abs(local_value) + Real(8) * tools::epsilon<Real>() * local_scale;
                // A few Newton steps' worth: copies of one root agree to about the local accuracy, and the
                // guaranteed radius (degree times the step) would also merge neighbouring roots. Copies
                // this misses are merged after polishing, with the guaranteed radius on the full polynomial.
@@ -601,7 +602,7 @@ std::vector<detail::complex_from_real_t<Real> > piecewise(std::vector<Real> coef
    // other's radius (sorted by real part, so there is no branch cut to split a pair), keeping the one with
    // the smaller local backward error, so that each root is polished against the full polynomial only
    // once. Using the smaller radius keeps a poor candidate from absorbing a good one for a different root.
-   Real eps = std::numeric_limits<Real>::epsilon();
+   Real eps = tools::epsilon<Real>();
    auto by_real = [](complex_type a, complex_type b) { return a.real() < b.real() || (a.real() == b.real() && a.imag() < b.imag()); };
    candidates.erase(std::remove_if(candidates.begin(), candidates.end(), [](candidate const& c) {
       return !(BOOST_MATH_NAMESPACE::isfinite)(c.z.real()) || !(BOOST_MATH_NAMESPACE::isfinite)(c.z.imag()) || !(c.radius < std::numeric_limits<Real>::infinity()); }),
@@ -880,7 +881,7 @@ class companion
    std::vector<complex_type> eigenvalues()
    {
       std::size_t n = n_;
-      Real eps = std::numeric_limits<Real>::epsilon();
+      Real eps = tools::epsilon<Real>();
       std::size_t hi = n - 1;
       unsigned its = 0, restarts = 0;
       while(hi > 0)

nanobind/tests/test_smoke.py

@@ -347,7 +347,9 @@ def test_integer_results():
 def test_python_scalars_give_python_floats():
     assert type(bm.libeta(1, 2, 0.1)) is float
     assert type(bm.tgamma(np.float32(3))) is np.float32
-    assert bm.libeta(1, 2, 0) == -np.inf
+    # log(0) is an exact infinity, which goes to the overflow policy:
+    with pytest.raises(OverflowError):
+        bm.libeta(1, 2, 0)
 
 
 # --- conventions ------------------------------------------------------------