[Numpy-discussion] Accessing real and imaginary parts of npy_complex(nbits)

[Brendan Murphy via NumPy-Discussion]

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);
  }
};
_______________________________________________
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com