Hi, This is my first look at the numpy internals, so I'm trying to help update pytensor to be compatible with numpy 2.0, and we have some structs that inherit from npy_complex64 and npy_complex128, which currently use .real and .imag to access real and imaginary parts. Assuming a 64 bit system, it is easy to update the code using npy_crealf (etc) for the struct that inherits from npy_complex64, and so on. (I'll put an example of the code at the bottom of the post.) Since numpy doesn't assume a 64 bit system, I'm writing some aliases for npy_crealf etc., depending on NPY_SIZEOF_FLOAT etc. I'm wondering if there is a smarter way to do this. Also, the pytensor code redefines complex arithmetic in terms of the standard "math" definitions. For C99 complex types, this can be achieved using #pragma STDC CX_LIMITED_RANGE (in theory, but really depending on the compiler). Is there any way to ask numpy to use this directive? (Googling says this makes complex arithmetic 3-5 times faster. The C99 complex arithmetic tries to avoid overflow, and this directive is only recommended if you know that won't be an issue. I don't know if we can assume that, but this code has been this way since it was added to Theano.) Example code: struct pytensor_complex64 : public npy_complex64 { typedef pytensor_complex64 complex_type; typedef npy_float32 scalar_type; complex_type operator+(const complex_type &y) const { complex_type ret; // ret.real = this->real + y.real; // ret.imag = this->imag + npy_cimagf(y); npy_csetrealf(&ret, npy_crealf(*this) + npy_crealf(y)); npy_csetimagf(&ret, npy_cimagf(*this) + npy_cimagf(y)); return ret; } complex_type operator-() const { complex_type ret; npy_csetrealf(&ret, -npy_crealf(*this)); npy_csetimagf(&ret, -npy_cimagf(*this)); return ret; } bool operator==(const complex_type &y) const { return (npy_crealf(*this) == npy_crealf(y)) && (npy_cimagf(*this) == npy_cimagf(y)); } bool operator==(const scalar_type &y) const { return (npy_crealf(*this) == y) && (npy_crealf(*this) == 0); } complex_type operator-(const complex_type &y) const { complex_type ret; npy_csetrealf(&ret, npy_crealf(*this) - npy_crealf(y)); npy_csetimagf(&ret, npy_cimagf(*this) - npy_cimagf(y)); return ret; } complex_type operator*(const complex_type &y) const { complex_type ret; npy_csetrealf(&ret, npy_crealf(*this) * npy_crealf(y) - npy_cimagf(*this) * npy_cimagf(y)); npy_csetimagf(&ret, npy_crealf(*this) * npy_cimagf(y) + npy_cimagf(*this) * npy_crealf(y)); return ret; } complex_type operator/(const complex_type &y) const { complex_type ret; scalar_type y_norm_square = npy_crealf(y) * npy_crealf(y) + npy_cimagf(y) * npy_cimagf(y); npy_csetrealf(&ret, (npy_crealf(*this) * npy_crealf(y) + npy_cimagf(*this) * npy_cimagf(y)) / y_norm_square); npy_csetimagf(&ret, (npy_crealf(*this) * npy_crealf(y) - npy_cimagf(*this) * npy_cimagf(y)) / y_norm_square); return ret; } template <typename T> complex_type &operator=(const T &y); pytensor_complex64() {} template <typename T> pytensor_complex64(const T &y) { *this = y; } template <typename TR, typename TI> pytensor_complex64(const TR &r, const TI &i) { npy_csetrealf(this, r); npy_csetimagf(this, i); } };