diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp
index 6c8944b331..67824ce100 100644
--- a/include/boost/math/distributions/extreme_value.hpp
+++ b/include/boost/math/distributions/extreme_value.hpp
@@ -190,8 +190,6 @@ 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;
@@ -199,6 +197,9 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const extreme_value_distribution<R
       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,8 +266,6 @@ 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;
@@ -274,6 +273,9 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<extreme_v
       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;
 
diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp
index 01a5169fd5..24466bb6d8 100644
--- a/include/boost/math/distributions/laplace.hpp
+++ b/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;
 
diff --git a/include/boost/math/distributions/logistic.hpp b/include/boost/math/distributions/logistic.hpp
index fcecb35be9..807503783b 100644
--- a/include/boost/math/distributions/logistic.hpp
+++ b/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()))
        {
diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp
index 9dd59075a3..466e949cbd 100644
--- a/include/boost/math/distributions/students_t.hpp
+++ b/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
    {
diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp
index 31d21f416d..0b44488eed 100644
--- a/include/boost/math/distributions/weibull.hpp
+++ b/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());
    }
diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index 1b832b3a86..dfe19ee4b4 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/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%)");
 }
 
diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
index efda6b2d0d..de3de96c13 100644
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/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%)");
 }
 
diff --git a/include/boost/math/tools/polynomial_roots.hpp b/include/boost/math/tools/polynomial_roots.hpp
index a29bcd50d0..8d28463276 100644
--- a/include/boost/math/tools/polynomial_roots.hpp
+++ b/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)
diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py
index a1b6d4bf6e..c232eb96d1 100644
--- a/nanobind/tests/test_smoke.py
+++ b/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 ------------------------------------------------------------
diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp
index 2db591905a..f257b4b92f 100644
--- a/test/test_extreme_value.cpp
+++ b/test/test_extreme_value.cpp
@@ -251,8 +251,13 @@ void test_spots(RealType)
       BOOST_CHECK_EQUAL(cdf(extreme_value_distribution<RealType>(), inf), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), -inf)), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), inf)), 0);
-      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), -inf), 0);
-      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), inf), 1);
+      // log(0) is an exact -infinity, which goes to the overflow policy:
+      BOOST_CHECK_THROW(logcdf(extreme_value_distribution<RealType>(), -inf), std::overflow_error);
+      BOOST_CHECK_THROW(logcdf(complement(extreme_value_distribution<RealType>(), inf)), std::overflow_error);
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), inf), 0);
+      BOOST_CHECK_EQUAL(logcdf(complement(extreme_value_distribution<RealType>(), -inf)), 0);
+      typedef boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error> > ignore_overflow;
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType, ignore_overflow>(), -inf), -inf);
 
       // Test NaN policy
       using boost::math::policies::policy;
diff --git a/test/test_lbeta.cpp b/test/test_lbeta.cpp
index 895453a2d1..55a06f3f36 100644
--- a/test/test_lbeta.cpp
+++ b/test/test_lbeta.cpp
@@ -113,13 +113,18 @@ void test_spots(std::size_t ulps, Real absolute_tolerance)
       CHECK_THROW(lbeta(Real(2), nan), std::domain_error);
 #endif
    }
-   // log(beta) is exactly -infinity when either argument is infinite, so that's returned, not an overflow:
+   // log(beta) is exactly -infinity when either argument is infinite. That's an overflow, as for every
+   // exact infinity; ignoring it gives -infinity:
    if (std::numeric_limits<Real>::has_infinity)
    {
+      using ignore_overflow = boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>>;
       const Real infinity = std::numeric_limits<Real>::infinity();
-      CHECK_EQUAL(Real(lbeta(infinity, Real(2.5))), Real(-infinity));
-      CHECK_EQUAL(Real(lbeta(Real(1e-30), infinity)), Real(-infinity));
-      CHECK_EQUAL(Real(lbeta(infinity, infinity)), Real(-infinity));
+      CHECK_EQUAL(Real(lbeta(infinity, Real(2.5), ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(lbeta(Real(1e-30), infinity, ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(lbeta(infinity, infinity, ignore_overflow())), Real(-infinity));
+#ifndef BOOST_NO_EXCEPTIONS
+      CHECK_THROW(lbeta(infinity, Real(2.5)), std::overflow_error);
+#endif
    }
 }
 
diff --git a/test/test_libeta.cpp b/test/test_libeta.cpp
index d142500658..c5d6397e63 100644
--- a/test/test_libeta.cpp
+++ b/test/test_libeta.cpp
@@ -89,14 +89,20 @@ void test_spots(std::size_t ulps, Real absolute_tolerance)
    CHECK_EQUAL(Real(libetac(Real(2), Real(3), Real(0))), Real(0));
    CHECK_EQUAL(Real(libeta(Real(0), Real(3), Real(0))), Real(0));
    CHECK_EQUAL(Real(libetac(Real(3), Real(0), Real(1))), Real(0));
+   // The log of an exact zero is an overflow, as for every exact infinity; ignoring it gives -infinity:
    if (std::numeric_limits<Real>::has_infinity)
    {
+      using ignore_overflow = boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>>;
       const Real minus_infinity = -std::numeric_limits<Real>::infinity();
-      CHECK_EQUAL(Real(libeta(Real(2), Real(3), Real(0))), minus_infinity);
-      CHECK_EQUAL(Real(libetac(Real(2), Real(3), Real(1))), minus_infinity);
-      CHECK_EQUAL(Real(libeta(Real(2), Real(0), Real(0.5))), minus_infinity);
-      CHECK_EQUAL(Real(libetac(Real(0), Real(2), Real(0.5))), minus_infinity);
+      CHECK_EQUAL(Real(libeta(Real(2), Real(3), Real(0), ignore_overflow())), minus_infinity);
+      CHECK_EQUAL(Real(libetac(Real(2), Real(3), Real(1), ignore_overflow())), minus_infinity);
+      CHECK_EQUAL(Real(libeta(Real(2), Real(0), Real(0.5), ignore_overflow())), minus_infinity);
+      CHECK_EQUAL(Real(libetac(Real(0), Real(2), Real(0.5), ignore_overflow())), minus_infinity);
    }
+#ifndef BOOST_NO_EXCEPTIONS
+   CHECK_THROW(libeta(Real(2), Real(3), Real(0)), std::overflow_error);
+   CHECK_THROW(libetac(Real(2), Real(3), Real(1)), std::overflow_error);
+#endif
    // NaN arguments are domain errors; ignoring the error gives NaN:
    if (std::numeric_limits<Real>::has_quiet_NaN)
    {
diff --git a/test/test_ligamma.cpp b/test/test_ligamma.cpp
index 01262ef62c..b3dfe530dd 100644
--- a/test/test_ligamma.cpp
+++ b/test/test_ligamma.cpp
@@ -125,14 +125,26 @@ void test_spots(std::size_t ulps, Real absolute_tolerance)
       }
    }
 
-   // At the ends of the range the functions are lgamma(a) or exactly -infinity, which is returned, not an overflow:
+   // At the ends of the range the functions are lgamma(a) or exactly -infinity. The latter is an overflow,
+   // as for every exact infinity, and so are lgamma_p and lgamma_q there; ignoring it gives -infinity:
    CHECK_EQUAL(Real(ligamma(Real(2.5), Real(0))), Real(boost::math::lgamma(Real(2.5))));
+#ifndef BOOST_NO_EXCEPTIONS
+   CHECK_THROW(ligamma_lower(Real(2.5), Real(0)), std::overflow_error);
+   CHECK_THROW(boost::math::lgamma_p(Real(2.5), Real(0)), std::overflow_error);
+#endif
    if (std::numeric_limits<Real>::has_infinity)
    {
+      using ignore_overflow = boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>>;
       const Real infinity = std::numeric_limits<Real>::infinity();
       CHECK_EQUAL(Real(ligamma_lower(Real(2.5), infinity)), Real(boost::math::lgamma(Real(2.5))));
-      CHECK_EQUAL(Real(ligamma(Real(2.5), infinity)), Real(-infinity));
-      CHECK_EQUAL(Real(ligamma_lower(Real(2.5), Real(0))), Real(-infinity));
+      CHECK_EQUAL(Real(ligamma(Real(2.5), infinity, ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(ligamma_lower(Real(2.5), Real(0), ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(boost::math::lgamma_p(Real(2.5), Real(0), ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(boost::math::lgamma_q(Real(2.5), infinity, ignore_overflow())), Real(-infinity));
+#ifndef BOOST_NO_EXCEPTIONS
+      CHECK_THROW(ligamma(Real(2.5), infinity), std::overflow_error);
+      CHECK_THROW(boost::math::lgamma_q(Real(2.5), infinity), std::overflow_error);
+#endif
    }
    // NaN arguments are domain errors; ignoring the error gives NaN:
    if (std::numeric_limits<Real>::has_quiet_NaN)
diff --git a/test/test_logcdf.cpp b/test/test_logcdf.cpp
index b96cd680c9..fa699f37ad 100644
--- a/test/test_logcdf.cpp
+++ b/test/test_logcdf.cpp
@@ -12,20 +12,24 @@
 #include <boost/math/distributions/beta.hpp>
 #include <boost/math/distributions/binomial.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
+#include <boost/math/distributions/extreme_value.hpp>
 #include <boost/math/distributions/fisher_f.hpp>
 #include <boost/math/distributions/gamma.hpp>
 #include <boost/math/distributions/inverse_chi_squared.hpp>
 #include <boost/math/distributions/inverse_gamma.hpp>
+#include <boost/math/distributions/laplace.hpp>
+#include <boost/math/distributions/logistic.hpp>
 #include <boost/math/distributions/negative_binomial.hpp>
 #include <boost/math/distributions/poisson.hpp>
 #include <boost/math/distributions/students_t.hpp>
+#include <boost/math/distributions/weibull.hpp>
 #include <cmath>
 
 using namespace boost::math;
 
 // Some references (log(1 - tiny)) are below the range of T, where long double is double.
 template <class T>
-void compare(T got, long double expected)
+void compare(T got, long double expected, long double condition = 1)
 {
    if (std::abs(expected) < (std::numeric_limits<T>::min)())
    {
@@ -35,20 +39,22 @@ void compare(T got, long double expected)
    {
       // log(1 - P) = -P for a tiny tail P, so the error is that of gamma_p/ibeta deep in the tail, which in
       // 80 and 128-bit long double reaches ~1000 ulps.
-      CHECK_ULP_CLOSE(static_cast<T>(expected), got, 2048);
+      CHECK_ULP_CLOSE(static_cast<T>(expected), got, static_cast<unsigned>((std::max)(2048.0L, 4 * condition)));
    }
    else
    {
       // ibeta's own error is the limit at large parameters (82 ulps for beta(200, 300) at 0.01).
-      CHECK_ULP_CLOSE(static_cast<T>(expected), got, 128);
+      CHECK_ULP_CLOSE(static_cast<T>(expected), got, static_cast<unsigned>((std::max)(128.0L, 4 * condition)));
    }
 }
 
+// One tail is log1p(-exp(L)) of the other's log L. Where exp(L) is tiny, that multiplies the error in L
+// by |L|, so the tolerance follows that condition number:
 template <class Dist, class T>
 void check(Dist const& dist, T x, long double log_p, long double log_q)
 {
-   compare(logcdf(dist, x), log_p);
-   compare(logcdf(complement(dist, x)), log_q);
+   compare(logcdf(dist, x), log_p, std::fabs(log_p) < 0.5L ? std::fabs(log_q) : 1);
+   compare(logcdf(complement(dist, x)), log_q, std::fabs(log_q) < 0.5L ? std::fabs(log_p) : 1);
 }
 
 template <class T>
@@ -90,12 +96,39 @@ void test_logcdf()
    check(students_t_distribution<T>(5), T(1e70), -9.490167245562357338012859910823596491189e-350L, -803.6545263121025131443400282712354083919L);
    check(students_t_distribution<T>(1e12), T(-40), -804.6084413729540376061252065587134160101L, -3.655895883611449397336474100007682207859e-350L);
    check(students_t_distribution<T>(1e12), T(40), -3.655895883611449397336474100007682207859e-350L, -804.6084413729540376061252065587134160101L);
+   // |x| > sqrt(max) in double, where df / x^2 may underflow:
+   check(students_t_distribution<T>(1), T(-1e200), -461.6617484846585369474748405007362614333L, -3.183098861837906811720139504195578159356e-201L);
+   check(students_t_distribution<T>(2.5), T(-1e160), -921.3634587421731073387873653405855378588L, -7.193397190831722297482216141085036842542e-401L);
+   check(students_t_distribution<T>(7), T(-1e155), -2493.42162960194500325806622545904493395L, -1.320520676354676868206420818690757852191e-1083L);
+   check(students_t_distribution<T>(0.5), T(-1e300), -346.5250100798783758198490653831970251013L, -3.207009754142228919211973798954086743678e-151L);
+   compare(logcdf(students_t_distribution<T>(1e300), T(-1e200)), -1.151292546497022901892142285621235491416e+302L);
    check(students_t_distribution<T>(1e300), T(-40), -804.6084420137537881666068329186099362001L, -3.655893540915029703748985802688283665054e-350L);
 }
 
+// At an infinite x, the log of the tail that vanishes is an exact -infinity, which goes to the overflow
+// policy, and the log of the other tail is zero:
+template <template <class, class> class Distribution, class T>
+void test_infinite_x()
+{
+   using Dist = Distribution<T, boost::math::policies::policy<>>;
+   using Quiet = Distribution<T, boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>>>;
+   const T inf = std::numeric_limits<T>::infinity();
+   CHECK_EQUAL(boost::math::logcdf(Dist(), inf), T(0));
+   CHECK_EQUAL(boost::math::logcdf(boost::math::complement(Dist(), -inf)), T(0));
+   CHECK_EQUAL(boost::math::logcdf(Quiet(), -inf), -inf);
+   CHECK_EQUAL(boost::math::logcdf(boost::math::complement(Quiet(), inf)), -inf);
+   CHECK_THROW(boost::math::logcdf(Dist(), -inf), std::overflow_error);
+   CHECK_THROW(boost::math::logcdf(boost::math::complement(Dist(), inf)), std::overflow_error);
+}
+
 int main()
 {
    test_logcdf<double>();
    test_logcdf<long double>();
+   test_infinite_x<boost::math::extreme_value_distribution, double>();
+   test_infinite_x<boost::math::laplace_distribution, double>();
+   test_infinite_x<boost::math::logistic_distribution, double>();
+   // The Weibull pdf vanishes at 0 for shape > 1:
+   CHECK_THROW(boost::math::logpdf(boost::math::weibull_distribution<double>(2), 0.0), std::overflow_error);
    return boost::math::test::report_errors();
 }