From 60134c06e9713becbd5cd9490a37117005239532 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Tue, 6 Feb 2024 18:26:35 +0000
Subject: [PATCH] Improve Bessel Y derivative coverage. Also fixes a bug in
 bessel_y_derivative_small_z_series corner case.

---
 .../detail/bessel_jy_derivatives_series.hpp      |  5 ++++-
 test/test_bessel_y_prime.cpp                     |  2 +-
 test/test_bessel_y_prime.hpp                     | 16 ++++++++++++++++
 3 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/include/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp b/include/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
index 8914e2c3f..d7bb64f8c 100644
--- a/include/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
+++ b/include/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
@@ -147,6 +147,7 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
    T p = log(x / 2);
    T scale = 1;
    bool need_logs = (v >= boost::math::max_factorial<T>::value) || (boost::math::tools::log_max_value<T>() / v < fabs(p));
+
    if (!need_logs)
    {
       gam = boost::math::tgamma(v, pol);
@@ -204,7 +205,9 @@ inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
    T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 
    result += scale * prefix * b;
-   return result;
+   if(scale * tools::max_value<T>() < result)
+      return boost::math::policies::raise_overflow_error<T>(function, nullptr, pol);
+   return result / scale;
 }
 
 // Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives
diff --git a/test/test_bessel_y_prime.cpp b/test/test_bessel_y_prime.cpp
index 20b4abbee..9848b6335 100644
--- a/test/test_bessel_y_prime.cpp
+++ b/test/test_bessel_y_prime.cpp
@@ -276,7 +276,7 @@ void expected_results()
       ".*",                          // platform
       largest_type,                  // test type(s)
       ".*",                          // test data group
-      ".*", 700, 300);                 // test function
+      ".*", 720, 600);                 // test function
    //
    // One set of float tests has inexact input values, so there is a slight error:
    //
diff --git a/test/test_bessel_y_prime.hpp b/test/test_bessel_y_prime.hpp
index 74600bece..1c0719e1d 100644
--- a/test/test_bessel_y_prime.hpp
+++ b/test/test_bessel_y_prime.hpp
@@ -225,5 +225,21 @@ void test_bessel_prime(T, const char* name)
     BOOST_CHECK_THROW(boost::math::cyl_neumann_prime(T(2.5), T(-1)), std::domain_error);
     BOOST_CHECK_THROW(boost::math::sph_neumann_prime(2, T(0)), std::domain_error);
     BOOST_CHECK_THROW(boost::math::sph_neumann_prime(2, T(-1)), std::domain_error);
+
+    BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity && (std::numeric_limits<T>::min_exponent < -1072))
+    {
+       static const std::array<std::array<T, 3>, 5> yv_prime_coverage_data = { {
+          {{ SC_(170.25), SC_(2), SC_(4.1990285871978876642542582761856953686528755802132926772620e306) }},
+    #if LDBL_MAX_10_EXP > 4936
+          {{ SC_(14.25), ldexp(T(1), -1072), SC_(1.830622575805420777640505999291582751497710200210963689466e4936)}},
+    #else
+          {{ SC_(14.25), ldexp(T(1), -1072), std::numeric_limits<T>::infinity() }},
+    #endif
+          {{ SC_(14.25), ldexp(T(1), -1074), std::numeric_limits<T>::infinity() }},
+          {{ SC_(15.25), ldexp(T(1), -1074), std::numeric_limits<T>::infinity() }},
+          {{ SC_(30.25), ldexp(T(1), -1045), std::numeric_limits<T>::infinity() }},
+       } };
+       do_test_cyl_neumann_y_prime<T>(yv_prime_coverage_data, name, "y': Extra coverage data");
+    }
 }