Skip to content

Replace uses of long with more portable types #53

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: main
Choose a base branch
from

Conversation

SabrinaJewson
Copy link
Contributor

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, ints 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.

`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.
@izaid
Copy link
Contributor

izaid commented Jul 26, 2025

So there were actually reasons we went with long as opposed to int or size_t. I'm not against changing them, but I think it should be thought through carefully beforehand. I also don't think it's the most important thing that needs to happen.

With regards to size_t, the reason it can be important to keep things signed as that some of these are used in calculation (e.g. as Bessel orders). And we should be able to do -n or something without worrying we'll get a wrong result due to being unsigned.

@SabrinaJewson
Copy link
Contributor Author

SabrinaJewson commented Jul 26, 2025

For context, I’m trying to use xsf for embedding in a different language, so I’m just trying to make things as easy as possible for this purpose. uint64_t makes things much easier than long, because most languages that aren’t C support fixed-width integer types more easily than C’s platform-dependent ones. The alternative is adding more error conditions to my function, which I’m a little averse to, since it makes it harder for users.

With regards to size_t, the reason it can be important to keep things signed as that some of these are used in calculation (e.g. as Bessel orders). And we should be able to do -n or something without worrying we'll get a wrong result due to being unsigned.

I am sympathic to this, but again it makes it a lot more difficult to use from other languages who would expect an unsigned type to be used here – that is, I favour the unsigned version because it provides strictly more functionality than the signed version.

I suppose it comes down to whether the vision for this library is to provide a universally applicable suite of special functions – in which case uint64_t would be preferable – or to be convenient to use from C – in which case a signed type would be preferable. I would of course personally think the former is more valuable, but I understand if it’s considered out of scope.

Edit: Maybe Scipy in particular can try adopting -Wsign-conversion?

@steppi
Copy link
Collaborator

steppi commented Jul 26, 2025

Thanks @SabrinaJewson. I think @izaid just wants to be cautious about changing things without being aware of their downstream impact in SciPy etc. For these particular cases, I don't think the use of long instead of a fixed width type in lambertw and sph_bessel_j was an intentional decision. It's just what originally appeared in the SciPy code that we started with. It was translated to C++ from separate C, Cython, and Fortran codebases, tens of thousands of lines at a time, and there is still work to be done to better unify the xsf codebase. It would be good to be more consistent around integer types.

It's good you pointed this out. The use of long here is actually problematic in SciPy as well. This is because the NumPy default integer changed in version 2.0 to be 32 bit on 32 bit platforms and 64 bit on 64 bit platforms. Previously the NumPy default integer was long. On Windows, a long is 32 bit. This means that on 64 bit Windows, ufuncs in scipy.special with parameters that only take long will error out when given arrays with type the default NumPy integer. There's actually a hack in place that makes things work for lambertw, which is only possible because lambertw is a wrapper around a ufunc _lambertw.

See below that the wrapper converts the array k to type long.

https://github.com/scipy/scipy/blob/7e055c1bf1a47e1d32f40edb6538422a97c28c8b/scipy/special/_lambertw.py#L146-L149

    # TODO: special expert should inspect this
    # interception; better place to do it?
    k = np.asarray(k, dtype=np.dtype("long"))
    return _lambertw(z, k, tol)

sph_bessel_j also goes through a wrapper which converts the integer array input to type long. So really, I think the only reason long remains for these two functions is because a fix was hacked in on the SciPy side.

For functions in scipy.special which are ufuncs instead of wrappers around them, the working solution has been to use type Py_ssize_t on the C side. See scipy/scipy#21401. In C++, this would correspond to std::ptrdiff_t, at least on all platforms SciPy supports. When this is done, 64 bit ints can't be passed to such ufuncs on 32 bit systems though.

This stuff came up recently in SciPy here, scipy/scipy#22795 (comment). I don't think we can just always use int64_t like you have here in this PR, because there seems to be some desire to support using native 32 bit integers on 32 bit systems. My preference is use templates and type traits to conveniently offer overloads for both int32_t and int64_t, so native ints could be used on 32 bit and 64 bit systems. We could do something like

template<typename T>
constexpr bool is_supported_int_v =
    std::is_same_v<T, int32_t> || std::is_same_v<T, int64_t>;

template<typename T>
using enable_if_supported_int_t = std::enable_if_t<is_supported_int_v<T>, int>;

and then

template <typename T, enable_if_supported_int_t<T> = 0>
XSF_HOST_DEVICE inline std::complex<float> lambertw(std::complex<float> z, T k, float tol) {

As for signed vs unsigned integers. Can you elaborate on the problems you're running into wrapping C++ functions which take signed integers? If it was consistently documented that public functions in xsf always take signed integers, even when an argument should only take nonnegative values, is it something that you could work around? In these cases, the SciPy functions take signed integer arrays as input and return NaN (optionally producing a warning) when given a value less than zero. Using unsigned inputs in xsf would require adding wrappers to SciPy to preserve this behavior. Also, in C++, I think it's generally considered good style to avoid unsigned integers unless necessary, see https://www.learncpp.com/cpp-tutorial/unsigned-integers-and-why-to-avoid-them/ for instance. I'm not sure the algorithms we have can even accurately compute Bessel functions of absurdly large order, to the point where the smaller upper bound on signed ints matters.

@lucascolley
Copy link
Member

scipy/scipy#19462 is still open

@SabrinaJewson
Copy link
Contributor Author

SabrinaJewson commented Jul 26, 2025

Thanks for the long response! That seems very sensible.

Can you elaborate on the problems you're running into wrapping C++ functions which take signed integers?

In the language I’m trying to wrap in (Rust), xsf using signed types mean I have three options:

  1. Using signed types in my API, which is considered highly unidiomatic (unsigned types are preferable for nonnegative values).
  2. Crashing the program (i.e. considering it a programmer error) when a user passes in a value greater than 2 ^ (BITS − 1) − 1 to my function. Probably not ideal?
  3. Returning NaN from my function when the user passes in a too large value. This may seem to be quite counterintuitive behaviour. In this case, passing in values that high will probably error anyway, but I’m also interested in the other functions that accept nonnegative integers, e.g. bdtr.

The idiomatic thing to do is just to support the full range of the unsigned type, which is why I have a preference for it.

Using unsigned inputs in xsf would require adding wrappers to SciPy to preserve this behavior.

I suppose I was under the impression that Scipy in general already did this with _nonneg_int_or_fail, but I guess that’s not done for all functions, which makes sense.

Also, in C++, I think it's generally considered good style to avoid unsigned integers unless necessary, see https://www.learncpp.com/cpp-tutorial/unsigned-integers-and-why-to-avoid-them/ for instance.

Right – and that’s where the friction between C/C++ and other languages comes in. This guideline exists because of the syntactic limitations of C/C++, but from my perspective it’s odd that this should leak into public APIs.

In my opinion, -Wsign-conversion combined with UBSan’s checking of unsigned overflow solves the common problems with using unsigned in C/C++, since it brings the semantics in line with Rust’s (and in my experience, in Rust there are very rarely issues caused by using unsigned values). And in fact, once implicit conversion is disabled, unsigned may become better than signed on the basis that overflow isn’t instant UB anymore. I’m happy to make a PR to get things to work with this setup.

By the way, I’m perfectly happy for array lengths to be signed values, since that doesn’t actually affect the range (arrays can’t be that large anyway). So this is just about the integers that are accepted for the purpose of casting to floats.

@steppi
Copy link
Collaborator

steppi commented Jul 26, 2025

I see. Looking into it more. In the case of Spherical Bessel functions, restriction to nonnegative integers is due to limitations in the implementation, not due to inherent properties of Spherical Bessel functions. The order in principle can be a general complex number. This sort of thing comes up a lot, where there are domain restrictions that are due to limitations in an implementation. The convention in this case has been to return NaN and use set_error to signal the reason to the user in a configurable way. SciPy is configured to have set_error use its internal sf_error that you saw in #51. In this particular case, the reason behind the limitation could be better documented.

Since restriction to nonnegative ints is an implementation detail, I'd be hesitant to change the input type to an unsigned type just to match the domain limitation. If you insist on unsigned integers in your wrappers, of the options you gave above for dealing with unsigned integers which exceed the largest signed integer of the same width, 3) would be the most faithful to the conventions used in xsf.

@SabrinaJewson
Copy link
Contributor Author

Looking into it more. In the case of Spherical Bessel functions, restriction to nonnegative integers is due to limitations in the implementation, not due to inherent properties of Spherical Bessel functions. The order in principle can be a general complex number.

Oh, that’s a good point. In that case, yes I agree it should be signed.

What about cases like the binomial distribution mentioned above?

Additionally, add checks in `sph_bessel_*_jac` functions to prevent
overflow when calculating `n - 1`.
@SabrinaJewson
Copy link
Contributor Author

I updated the PR to use the enable_if_supported_int_t approach. I also added checks for n < 0 in the _jac functions, since otherwise we might overflow n.

@SabrinaJewson
Copy link
Contributor Author

Huh. macOS CI is failing and I don’t know why… unfortunately I don’t have a Mac. CUDA is also failing, and I have absolutely zero experience with that… 😅

Any ideas from people who know C++ better than me? it just says this:

macOS build error
/Users/runner/work/xsf/xsf/./tests/scipy_special_tests/test_lambertw.cpp:21:16: error: no matching function for call to 'lambertw'
   21 |     auto out = xsf::lambertw(z, k, eps);
      |                ^~~~~~~~~~~~~
/Users/runner/work/xsf/xsf/./include/xsf/lambertw.h:68:45: note: candidate template ignored: requirement 'is_supported_int_v<long>' was not satisfied [with K = std::tuple_element<1, std::tuple<std::complex<double>, long, double>>::type]
   68 | XSF_HOST_DEVICE inline std::complex<double> lambertw(std::complex<double> z, K k, double tol) {
      |                                             ^
/Users/runner/work/xsf/xsf/./include/xsf/lambertw.h:147:44: note: candidate template ignored: requirement 'is_supported_int_v<long>' was not satisfied [with K = std::tuple_element<1, std::tuple<std::complex<double>, long, double>>::type]
  147 | XSF_HOST_DEVICE inline std::complex<float> lambertw(std::complex<float> z, K k, float tol) {
      |                                            ^
1 error generated.

@steppi
Copy link
Collaborator

steppi commented Jul 26, 2025

What about cases like the binomial distribution mentioned above?

There are analytic continuations there too. I remember this coming up here, scipy/scipy#2414 (comment). One complication is that when treating a function as CDF function of a distribution with compact support, evaluation at points in the domain greater than the support yields result 1, and evaluation at points less than the support gives output 0. But if you were to do extend the domain through analytic continuation, you get a different extension. There's a tension between extension through analytic continuation and extension by general statistical convention. To allow flexibility to libraries that want to wrap xsf from other languages, maybe it makes sense to offer both signed and unsigned overloads.

Huh. macOS CI is failing and I don’t know why… unfortunately I don’t have a Mac. CUDA is also failing, and I have absolutely zero experience with that… 😅

Any ideas from people who know C++ better than me? it just says this:

The CuPy tests are failing because the version of CuPy in CI is still stuck on C++11, but the changes require C++17. I bumped the C++ language standard for special functions in CuPy to C++17 in cupy/cupy#9159, but this hasn't made it into a CuPy release yet. I'll think of something.

I'm not sure what's happening on macOS, but I have a Mac mini I can test this out on.

@izaid
Copy link
Contributor

izaid commented Jul 27, 2025

The CuPy tests are failing because the version of CuPy in CI is still stuck on C++11, but the changes require C++17. I bumped the C++ language standard for special functions in CuPy to C++17 in cupy/cupy#9159, but this hasn't made it into a CuPy release yet. I'll think of something.

I'm not sure what's happening on macOS, but I have a Mac mini I can test this out on.

I believe you can pass in an argument to the CuPy ElementwiseKernel that enables C++17.

@steppi
Copy link
Collaborator

steppi commented Jul 27, 2025

I believe you can pass in an argument to the CuPy ElementwiseKernel that enables C++17.

Exactly. The tests currently use cupy._core.create_ufunc, which doesn't have that option, but we can change them to use ElementwiseKernel directly.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants