Replace uses of long with more portable types

`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.

Author: SabrinaJewson

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<double> lambertw_asy(std::complex<double> z, long k) {
+    XSF_HOST_DEVICE inline std::complex<double> lambertw_asy(std::complex<double> 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<double> lambertw(std::complex<double> z, long k, double tol) {
+XSF_HOST_DEVICE inline std::complex<double> lambertw(std::complex<double> z, int64_t k, double tol) {
     double absz;
     std::complex<double> w;
     std::complex<double> ew, wew, wewz, wn;
@@ -142,7 +142,7 @@ XSF_HOST_DEVICE inline std::complex<double> lambertw(std::complex<double> z, lon
     return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
 }
 
-XSF_HOST_DEVICE inline std::complex<float> lambertw(std::complex<float> z, long k, float tol) {
+XSF_HOST_DEVICE inline std::complex<float> lambertw(std::complex<float> z, int64_t k, float tol) {
     return static_cast<std::complex<float>>(lambertw(static_cast<std::complex<double>>(z), k, static_cast<double>(tol))
     );
 }

include/xsf/sph_bessel.h

@@ -34,16 +34,11 @@ Translated to C++ by SciPy developers in 2024.
 namespace xsf {
 
 template <typename T>
-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<T>::quiet_NaN();
-    }
-
     if ((x == std::numeric_limits<T>::infinity()) || (x == -std::numeric_limits<T>::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 <typename T>
-std::complex<T> sph_bessel_j(long n, std::complex<T> z) {
+std::complex<T> sph_bessel_j(uint64_t n, std::complex<T> 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<T>::quiet_NaN();
-    }
-
     if (std::real(z) == std::numeric_limits<T>::infinity() || std::real(z) == -std::numeric_limits<T>::infinity()) {
         // https://dlmf.nist.gov/10.52.E3
         if (std::imag(z) == 0) {
@@ -121,7 +111,7 @@ std::complex<T> sph_bessel_j(long n, std::complex<T> z) {
 }
 
 template <typename T>
-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<T>(n + 1) * sph_bessel_j(n, z) / z;
+    return sph_bessel_j(n - 1, z) - (static_cast<T>(n) + 1) * sph_bessel_j(n, z) / z;
 }
 
 template <typename T>
-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<T>::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<T>::infinity() || x == -std::numeric_limits<T>::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<double>(x)); }
+inline float sph_bessel_y(uint64_t n, float x) { return sph_bessel_y(n, static_cast<double>(x)); }
 
 template <typename T>
-std::complex<T> sph_bessel_y(long n, std::complex<T> z) {
+std::complex<T> sph_bessel_y(uint64_t n, std::complex<T> 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<T>::quiet_NaN();
-    }
-
     if (std::real(z) == 0 && std::imag(z) == 0) {
         // https://dlmf.nist.gov/10.52.E2
         return std::numeric_limits<T>::quiet_NaN();
@@ -219,25 +199,20 @@ std::complex<T> sph_bessel_y(long n, std::complex<T> z) {
 }
 
 template <typename T>
-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<T>(n + 1) * sph_bessel_y(n, x) / x;
+    return sph_bessel_y(n - 1, x) - (static_cast<T>(n) + 1) * sph_bessel_y(n, x) / x;
 }
 
 template <typename T>
-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<T>::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<T>::infinity()) {
-            return std::pow(-1, n) * std::numeric_limits<T>::infinity();
+            return n % 2 == 0 ? std::numeric_limits<T>::infinity() : -std::numeric_limits<T>::infinity();
         }
 
         return std::numeric_limits<T>::infinity();
@@ -259,16 +234,11 @@ T sph_bessel_i(long n, T x) {
 }
 
 template <typename T>
-std::complex<T> sph_bessel_i(long n, std::complex<T> z) {
+std::complex<T> sph_bessel_i(uint64_t n, std::complex<T> 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<T>::quiet_NaN();
-    }
-
     if (std::abs(z) == 0) {
         // https://dlmf.nist.gov/10.52.E1
         if (n == 0) {
@@ -282,7 +252,7 @@ std::complex<T> sph_bessel_i(long n, std::complex<T> z) {
         // https://dlmf.nist.gov/10.52.E5
         if (std::imag(z) == 0) {
             if (std::real(z) == -std::numeric_limits<T>::infinity()) {
-                return std::pow(-1, n) * std::numeric_limits<T>::infinity();
+                return n % 2 == 0 ? std::numeric_limits<T>::infinity() : -std::numeric_limits<T>::infinity();
             }
 
             return std::numeric_limits<T>::infinity();
@@ -295,7 +265,7 @@ std::complex<T> sph_bessel_i(long n, std::complex<T> z) {
 }
 
 template <typename T>
-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<T>(n + 1) * sph_bessel_i(n, z) / z;
+    return sph_bessel_i(n - 1, z) - (static_cast<T>(n) + 1) * sph_bessel_i(n, z) / z;
 }
 
 template <typename T>
-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<T>::quiet_NaN();
-    }
-
     if (z == 0) {
         return std::numeric_limits<T>::infinity();
     }
@@ -339,16 +304,11 @@ T sph_bessel_k(long n, T z) {
 }
 
 template <typename T>
-std::complex<T> sph_bessel_k(long n, std::complex<T> z) {
+std::complex<T> sph_bessel_k(uint64_t n, std::complex<T> 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<T>::quiet_NaN();
-    }
-
     if (std::abs(z) == 0) {
         return std::numeric_limits<T>::quiet_NaN();
     }
@@ -370,12 +330,12 @@ std::complex<T> sph_bessel_k(long n, std::complex<T> z) {
 }
 
 template <typename T>
-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<T>(n + 1) * sph_bessel_k(n, x) / x;
+    return -sph_bessel_k(n - 1, x) - (static_cast<T>(n) + 1) * sph_bessel_k(n, x) / x;
 }
 
 } // namespace xsf