From a1c17b0a07c50f454be53d4ce36c5881cbe5725d Mon Sep 17 00:00:00 2001 From: SabrinaJewson Date: Sat, 26 Jul 2025 01:36:57 +0100 Subject: [PATCH 1/2] Replace uses of `long` with more portable types MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `long` is unportable, since its width depends a lot on the target platform: it’s 32-bit on 32 bit platforms and Windows, but 64-bit elsewhere. `long k` in `lambertw.h` is replaced with `int64_t`, since it can be any integer and is castable to a float. `long n` in the spherical Bessel functions is replaced with `uint64_t`. This change also fixes three bugs: + Previously, `int`s were used in loops up to `n`, but this is the wrong type. This is changed to use `uint64_t`. + `n + 1` was used liberally before, but this causes UB when `n` is its maximum value for the type. We instead cast to the float type before adding one, thus avoiding UB (or in the case of the new unsigned type, overflow). + `std::pow(-1, n)` was previously used – however, this is incorrect when `n` is a `long`, since it will be implicitly converted to an `int`. We replace these with a simple ternary to control the sign. --- include/xsf/lambertw.h | 6 +-- include/xsf/sph_bessel.h | 84 +++++++++++----------------------------- 2 files changed, 25 insertions(+), 65 deletions(-) diff --git a/include/xsf/lambertw.h b/include/xsf/lambertw.h index eeeca8aef7..bacfeeee73 100644 --- a/include/xsf/lambertw.h +++ b/include/xsf/lambertw.h @@ -54,7 +54,7 @@ namespace detail { return z * cevalpoly(num, 2, z) / cevalpoly(denom, 2, z); } - XSF_HOST_DEVICE inline std::complex lambertw_asy(std::complex z, long k) { + XSF_HOST_DEVICE inline std::complex lambertw_asy(std::complex z, int64_t k) { /* Compute the W function using the first two terms of the * asymptotic series. See 4.20 in [1]. */ @@ -64,7 +64,7 @@ namespace detail { } // namespace detail -XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, long k, double tol) { +XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, int64_t k, double tol) { double absz; std::complex w; std::complex ew, wew, wewz, wn; @@ -142,7 +142,7 @@ XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, lon return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } -XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, long k, float tol) { +XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, int64_t k, float tol) { return static_cast>(lambertw(static_cast>(z), k, static_cast(tol)) ); } diff --git a/include/xsf/sph_bessel.h b/include/xsf/sph_bessel.h index ea1105c4da..fbebc9a1cf 100644 --- a/include/xsf/sph_bessel.h +++ b/include/xsf/sph_bessel.h @@ -34,16 +34,11 @@ Translated to C++ by SciPy developers in 2024. namespace xsf { template -T sph_bessel_j(long n, T x) { +T sph_bessel_j(uint64_t n, T x) { if (std::isnan(x)) { return x; } - if (n < 0) { - set_error("spherical_jn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if ((x == std::numeric_limits::infinity()) || (x == -std::numeric_limits::infinity())) { return 0; } @@ -71,7 +66,7 @@ T sph_bessel_j(long n, T x) { } T sn; - for (int i = 0; i < n - 1; ++i) { + for (uint64_t i = 0; i < n - 1; ++i) { sn = (2 * i + 3) * s1 / x - s0; s0 = s1; s1 = sn; @@ -85,16 +80,11 @@ T sph_bessel_j(long n, T x) { } template -std::complex sph_bessel_j(long n, std::complex z) { +std::complex sph_bessel_j(uint64_t n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } - if (n < 0) { - set_error("spherical_jn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (std::real(z) == std::numeric_limits::infinity() || std::real(z) == -std::numeric_limits::infinity()) { // https://dlmf.nist.gov/10.52.E3 if (std::imag(z) == 0) { @@ -121,7 +111,7 @@ std::complex sph_bessel_j(long n, std::complex z) { } template -T sph_bessel_j_jac(long n, T z) { +T sph_bessel_j_jac(uint64_t n, T z) { if (n == 0) { return -sph_bessel_j(1, z); } @@ -136,25 +126,20 @@ T sph_bessel_j_jac(long n, T z) { } // DLMF 10.51.2 - return sph_bessel_j(n - 1, z) - static_cast(n + 1) * sph_bessel_j(n, z) / z; + return sph_bessel_j(n - 1, z) - (static_cast(n) + 1.0) * sph_bessel_j(n, z) / z; } template -T sph_bessel_y(long n, T x) { +T sph_bessel_y(uint64_t n, T x) { T s0, s1, sn; - int idx; + uint64_t idx; if (std::isnan(x)) { return x; } - if (n < 0) { - set_error("spherical_yn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (x < 0) { - return std::pow(-1, n + 1) * sph_bessel_y(n, -x); + return n % 2 == 0 ? -sph_bessel_y(n, -x) : sph_bessel_y(n, -x); } if (x == std::numeric_limits::infinity() || x == -std::numeric_limits::infinity()) { @@ -188,19 +173,14 @@ T sph_bessel_y(long n, T x) { return sn; } -inline float sph_bessel_y(long n, float x) { return sph_bessel_y(n, static_cast(x)); } +inline float sph_bessel_y(uint64_t n, float x) { return sph_bessel_y(n, static_cast(x)); } template -std::complex sph_bessel_y(long n, std::complex z) { +std::complex sph_bessel_y(uint64_t n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } - if (n < 0) { - set_error("spherical_yn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (std::real(z) == 0 && std::imag(z) == 0) { // https://dlmf.nist.gov/10.52.E2 return std::numeric_limits::quiet_NaN(); @@ -219,25 +199,20 @@ std::complex sph_bessel_y(long n, std::complex z) { } template -T sph_bessel_y_jac(long n, T x) { +T sph_bessel_y_jac(uint64_t n, T x) { if (n == 0) { return -sph_bessel_y(1, x); } - return sph_bessel_y(n - 1, x) - static_cast(n + 1) * sph_bessel_y(n, x) / x; + return sph_bessel_y(n - 1, x) - (static_cast(n) + 1.0) * sph_bessel_y(n, x) / x; } template -T sph_bessel_i(long n, T x) { +T sph_bessel_i(uint64_t n, T x) { if (std::isnan(x)) { return x; } - if (n < 0) { - set_error("spherical_in", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (x == 0) { // https://dlmf.nist.gov/10.52.E1 if (n == 0) { @@ -249,7 +224,7 @@ T sph_bessel_i(long n, T x) { if (std::isinf(x)) { // https://dlmf.nist.gov/10.49.E8 if (x == -std::numeric_limits::infinity()) { - return std::pow(-1, n) * std::numeric_limits::infinity(); + return n % 2 == 0 ? std::numeric_limits::infinity() : -std::numeric_limits::infinity(); } return std::numeric_limits::infinity(); @@ -259,16 +234,11 @@ T sph_bessel_i(long n, T x) { } template -std::complex sph_bessel_i(long n, std::complex z) { +std::complex sph_bessel_i(uint64_t n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } - if (n < 0) { - set_error("spherical_in", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (std::abs(z) == 0) { // https://dlmf.nist.gov/10.52.E1 if (n == 0) { @@ -282,7 +252,7 @@ std::complex sph_bessel_i(long n, std::complex z) { // https://dlmf.nist.gov/10.52.E5 if (std::imag(z) == 0) { if (std::real(z) == -std::numeric_limits::infinity()) { - return std::pow(-1, n) * std::numeric_limits::infinity(); + return n % 2 == 0 ? std::numeric_limits::infinity() : -std::numeric_limits::infinity(); } return std::numeric_limits::infinity(); @@ -295,7 +265,7 @@ std::complex sph_bessel_i(long n, std::complex z) { } template -T sph_bessel_i_jac(long n, T z) { +T sph_bessel_i_jac(uint64_t n, T z) { if (n == 0) { return sph_bessel_i(1, z); } @@ -308,20 +278,15 @@ T sph_bessel_i_jac(long n, T z) { } } - return sph_bessel_i(n - 1, z) - static_cast(n + 1) * sph_bessel_i(n, z) / z; + return sph_bessel_i(n - 1, z) - (static_cast(n) + 1.0) * sph_bessel_i(n, z) / z; } template -T sph_bessel_k(long n, T z) { +T sph_bessel_k(uint64_t n, T z) { if (std::isnan(z)) { return z; } - if (n < 0) { - set_error("spherical_kn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (z == 0) { return std::numeric_limits::infinity(); } @@ -339,16 +304,11 @@ T sph_bessel_k(long n, T z) { } template -std::complex sph_bessel_k(long n, std::complex z) { +std::complex sph_bessel_k(uint64_t n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } - if (n < 0) { - set_error("spherical_kn", SF_ERROR_DOMAIN, nullptr); - return std::numeric_limits::quiet_NaN(); - } - if (std::abs(z) == 0) { return std::numeric_limits::quiet_NaN(); } @@ -370,12 +330,12 @@ std::complex sph_bessel_k(long n, std::complex z) { } template -T sph_bessel_k_jac(long n, T x) { +T sph_bessel_k_jac(uint64_t n, T x) { if (n == 0) { return -sph_bessel_k(1, x); } - return -sph_bessel_k(n - 1, x) - static_cast(n + 1) * sph_bessel_k(n, x) / x; + return -sph_bessel_k(n - 1, x) - (static_cast(n) + 1.0) * sph_bessel_k(n, x) / x; } } // namespace xsf From 7c2348dd51e4bc98a0f80eb57735068438dc516a Mon Sep 17 00:00:00 2001 From: SabrinaJewson Date: Sat, 26 Jul 2025 21:07:28 +0100 Subject: [PATCH 2/2] Be generic over signed `N` integer types Additionally, add checks in `sph_bessel_*_jac` functions to prevent overflow when calculating `n - 1`. --- include/xsf/config.h | 11 ++++ include/xsf/lambertw.h | 8 ++- include/xsf/sph_bessel.h | 115 ++++++++++++++++++++++++++++++--------- 3 files changed, 104 insertions(+), 30 deletions(-) diff --git a/include/xsf/config.h b/include/xsf/config.h index ec7f0bbff8..645bde5c9c 100644 --- a/include/xsf/config.h +++ b/include/xsf/config.h @@ -53,6 +53,17 @@ #define M_SQRT1_2 0.707106781186547524401 #endif +namespace xsf { + +template +constexpr bool is_supported_int_v = + std::is_same_v || std::is_same_v; + +template +using enable_if_supported_int_t = std::enable_if_t, int>; + +} + #ifdef __CUDACC__ #define XSF_HOST_DEVICE __host__ __device__ diff --git a/include/xsf/lambertw.h b/include/xsf/lambertw.h index bacfeeee73..39ca1bab54 100644 --- a/include/xsf/lambertw.h +++ b/include/xsf/lambertw.h @@ -54,7 +54,7 @@ namespace detail { return z * cevalpoly(num, 2, z) / cevalpoly(denom, 2, z); } - XSF_HOST_DEVICE inline std::complex lambertw_asy(std::complex z, int64_t k) { + XSF_HOST_DEVICE inline std::complex lambertw_asy(std::complex z, double k) { /* Compute the W function using the first two terms of the * asymptotic series. See 4.20 in [1]. */ @@ -64,7 +64,8 @@ namespace detail { } // namespace detail -XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, int64_t k, double tol) { +template = 0> +XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, K k, double tol) { double absz; std::complex w; std::complex ew, wew, wewz, wn; @@ -142,7 +143,8 @@ XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, int return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; } -XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, int64_t k, float tol) { +template = 0> +XSF_HOST_DEVICE inline std::complex lambertw(std::complex z, K k, float tol) { return static_cast>(lambertw(static_cast>(z), k, static_cast(tol)) ); } diff --git a/include/xsf/sph_bessel.h b/include/xsf/sph_bessel.h index fbebc9a1cf..40cc1fda60 100644 --- a/include/xsf/sph_bessel.h +++ b/include/xsf/sph_bessel.h @@ -33,12 +33,17 @@ Translated to C++ by SciPy developers in 2024. namespace xsf { -template -T sph_bessel_j(uint64_t n, T x) { +template = 0> +T sph_bessel_j(N n, T x) { if (std::isnan(x)) { return x; } + if (n < 0) { + set_error("spherical_jn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if ((x == std::numeric_limits::infinity()) || (x == -std::numeric_limits::infinity())) { return 0; } @@ -66,7 +71,7 @@ T sph_bessel_j(uint64_t n, T x) { } T sn; - for (uint64_t i = 0; i < n - 1; ++i) { + for (N i = 0; i < n - 1; ++i) { sn = (2 * i + 3) * s1 / x - s0; s0 = s1; s1 = sn; @@ -79,12 +84,17 @@ T sph_bessel_j(uint64_t n, T x) { return sn; } -template -std::complex sph_bessel_j(uint64_t n, std::complex z) { +template = 0> +std::complex sph_bessel_j(N n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } + if (n < 0) { + set_error("spherical_jn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (std::real(z) == std::numeric_limits::infinity() || std::real(z) == -std::numeric_limits::infinity()) { // https://dlmf.nist.gov/10.52.E3 if (std::imag(z) == 0) { @@ -110,8 +120,13 @@ std::complex sph_bessel_j(uint64_t n, std::complex z) { return out; } -template -T sph_bessel_j_jac(uint64_t n, T z) { +template = 0> +T sph_bessel_j_jac(N n, T z) { + if (n < 0) { + set_error("spherical_j_jac", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (n == 0) { return -sph_bessel_j(1, z); } @@ -129,15 +144,20 @@ T sph_bessel_j_jac(uint64_t n, T z) { return sph_bessel_j(n - 1, z) - (static_cast(n) + 1.0) * sph_bessel_j(n, z) / z; } -template -T sph_bessel_y(uint64_t n, T x) { +template = 0> +T sph_bessel_y(N n, T x) { T s0, s1, sn; - uint64_t idx; + N idx; if (std::isnan(x)) { return x; } + if (n < 0) { + set_error("spherical_yn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (x < 0) { return n % 2 == 0 ? -sph_bessel_y(n, -x) : sph_bessel_y(n, -x); } @@ -173,14 +193,20 @@ T sph_bessel_y(uint64_t n, T x) { return sn; } -inline float sph_bessel_y(uint64_t n, float x) { return sph_bessel_y(n, static_cast(x)); } +template = 0> +inline float sph_bessel_y(N n, float x) { return sph_bessel_y(n, static_cast(x)); } -template -std::complex sph_bessel_y(uint64_t n, std::complex z) { +template = 0> +std::complex sph_bessel_y(N n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } + if (n < 0) { + set_error("spherical_yn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (std::real(z) == 0 && std::imag(z) == 0) { // https://dlmf.nist.gov/10.52.E2 return std::numeric_limits::quiet_NaN(); @@ -198,8 +224,13 @@ std::complex sph_bessel_y(uint64_t n, std::complex z) { return std::sqrt(static_cast(M_PI_2) / z) * cyl_bessel_y(n + 1 / static_cast(2), z); } -template -T sph_bessel_y_jac(uint64_t n, T x) { +template = 0> +T sph_bessel_y_jac(N n, T x) { + if (n < 0) { + set_error("spherical_y_jac", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (n == 0) { return -sph_bessel_y(1, x); } @@ -207,12 +238,17 @@ T sph_bessel_y_jac(uint64_t n, T x) { return sph_bessel_y(n - 1, x) - (static_cast(n) + 1.0) * sph_bessel_y(n, x) / x; } -template -T sph_bessel_i(uint64_t n, T x) { +template = 0> +T sph_bessel_i(N n, T x) { if (std::isnan(x)) { return x; } + if (n < 0) { + set_error("spherical_in", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (x == 0) { // https://dlmf.nist.gov/10.52.E1 if (n == 0) { @@ -233,12 +269,17 @@ T sph_bessel_i(uint64_t n, T x) { return sqrt(static_cast(M_PI_2) / x) * cyl_bessel_i(n + 1 / static_cast(2), x); } -template -std::complex sph_bessel_i(uint64_t n, std::complex z) { +template = 0> +std::complex sph_bessel_i(N n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } + if (n < 0) { + set_error("spherical_in", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (std::abs(z) == 0) { // https://dlmf.nist.gov/10.52.E1 if (n == 0) { @@ -264,8 +305,13 @@ std::complex sph_bessel_i(uint64_t n, std::complex z) { return std::sqrt(static_cast(M_PI_2) / z) * cyl_bessel_i(n + 1 / static_cast(2), z); } -template -T sph_bessel_i_jac(uint64_t n, T z) { +template = 0> +T sph_bessel_i_jac(N n, T z) { + if (n < 0) { + set_error("spherical_i_jac", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (n == 0) { return sph_bessel_i(1, z); } @@ -281,12 +327,17 @@ T sph_bessel_i_jac(uint64_t n, T z) { return sph_bessel_i(n - 1, z) - (static_cast(n) + 1.0) * sph_bessel_i(n, z) / z; } -template -T sph_bessel_k(uint64_t n, T z) { +template = 0> +T sph_bessel_k(N n, T z) { if (std::isnan(z)) { return z; } + if (n < 0) { + set_error("spherical_kn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (z == 0) { return std::numeric_limits::infinity(); } @@ -303,12 +354,17 @@ T sph_bessel_k(uint64_t n, T z) { return std::sqrt(M_PI_2 / z) * cyl_bessel_k(n + 1 / static_cast(2), z); } -template -std::complex sph_bessel_k(uint64_t n, std::complex z) { +template = 0> +std::complex sph_bessel_k(N n, std::complex z) { if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) { return z; } + if (n < 0) { + set_error("spherical_kn", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (std::abs(z) == 0) { return std::numeric_limits::quiet_NaN(); } @@ -329,8 +385,13 @@ std::complex sph_bessel_k(uint64_t n, std::complex z) { return std::sqrt(static_cast(M_PI_2) / z) * cyl_bessel_k(n + 1 / static_cast(2), z); } -template -T sph_bessel_k_jac(uint64_t n, T x) { +template = 0> +T sph_bessel_k_jac(N n, T x) { + if (n < 0) { + set_error("spherical_k_jac", SF_ERROR_DOMAIN, nullptr); + return std::numeric_limits::quiet_NaN(); + } + if (n == 0) { return -sph_bessel_k(1, x); }