From 1ec331ffd241dec5d40c82d30d2c342ce9686aab Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Wed, 29 Oct 2025 20:11:53 -0400 Subject: [PATCH 1/4] Short-circuit return on zero input in cyl_bessel_y() --- include/xsf/bessel.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/xsf/bessel.h b/include/xsf/bessel.h index 74473e8dbb..02947b00ab 100644 --- a/include/xsf/bessel.h +++ b/include/xsf/bessel.h @@ -948,6 +948,7 @@ inline std::complex cyl_bessel_y(double v, std::complex z) { cy_y.real(-INFINITY); cy_y.imag(0); set_error("yv", SF_ERROR_OVERFLOW, NULL); + return cy_y; } else { nz = amos::besy(z, v, kode, n, &cy_y, &ierr); set_error_and_nan("yv:", ierr_to_sferr(nz, ierr), cy_y); From 9b37f3819454f4667247094879f1399a205c32b0 Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Wed, 29 Oct 2025 20:30:19 -0400 Subject: [PATCH 2/4] Fix sign of INFINITY. --- include/xsf/bessel.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/xsf/bessel.h b/include/xsf/bessel.h index 02947b00ab..f2257c9e67 100644 --- a/include/xsf/bessel.h +++ b/include/xsf/bessel.h @@ -945,7 +945,7 @@ inline std::complex cyl_bessel_y(double v, std::complex z) { if (z.real() == 0 && z.imag() == 0) { /* overflow */ - cy_y.real(-INFINITY); + cy_y.real(-sign*INFINITY); cy_y.imag(0); set_error("yv", SF_ERROR_OVERFLOW, NULL); return cy_y; From dc432e3d944ba5b0ab2bb5365d336f86b39985cf Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Wed, 29 Oct 2025 20:32:18 -0400 Subject: [PATCH 3/4] Fix whitespace. --- include/xsf/bessel.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/xsf/bessel.h b/include/xsf/bessel.h index f2257c9e67..f830c96246 100644 --- a/include/xsf/bessel.h +++ b/include/xsf/bessel.h @@ -945,7 +945,7 @@ inline std::complex cyl_bessel_y(double v, std::complex z) { if (z.real() == 0 && z.imag() == 0) { /* overflow */ - cy_y.real(-sign*INFINITY); + cy_y.real(-sign * INFINITY); cy_y.imag(0); set_error("yv", SF_ERROR_OVERFLOW, NULL); return cy_y; From 0f94f1637f66c671c5ffafe3050e83952b9bfed3 Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Wed, 29 Oct 2025 21:44:30 -0400 Subject: [PATCH 4/4] Fix sign of INFINITY when v < 0 and z = 0 --- include/xsf/bessel.h | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/include/xsf/bessel.h b/include/xsf/bessel.h index f830c96246..e7f4e16831 100644 --- a/include/xsf/bessel.h +++ b/include/xsf/bessel.h @@ -945,7 +945,29 @@ inline std::complex cyl_bessel_y(double v, std::complex z) { if (z.real() == 0 && z.imag() == 0) { /* overflow */ - cy_y.real(-sign * INFINITY); + double re; + if (sign == 1 || v < 0.5) { + re = -INFINITY; + } else { + if (v >= std::pow(2.0, 53)) { + // The original v (before sign flip) is an even negative integer. + re = -INFINITY; + } else { + double s = v - 0.5; + if (s != v && std::rint(s) == s) { + re = 0.0; + } else { + double si = std::floor(s); + double half_si = si / 2; + if (std::floor(half_si) == half_si) { + re = INFINITY; + } else { + re = -INFINITY; + } + } + } + } + cy_y.real(re); cy_y.imag(0); set_error("yv", SF_ERROR_OVERFLOW, NULL); return cy_y;