Merge pull request numpy#24085 from lysnikolaou/real-imag-remove

ENH: Replace npy complex structs with native complex types

numpy/core/include/meson.build

@@ -8,6 +8,7 @@ installed_headers = [
   'numpy/ndarrayobject.h',
   'numpy/ndarraytypes.h',
   'numpy/npy_1_7_deprecated_api.h',
+  'numpy/npy_2_complexcompat.h',
   'numpy/npy_3kcompat.h',
   'numpy/npy_common.h',
   'numpy/npy_cpu.h',

numpy/core/include/numpy/npy_2_complexcompat.h

@@ -0,0 +1,28 @@
+/* This header is designed to be copy-pasted into downstream packages, since it provides
+   a compatibility layer between the old C struct complex types and the new native C99
+   complex types. The new macros are in numpy/npy_math.h, which is why it is included here. */
+#ifndef NUMPY_CORE_INCLUDE_NUMPY_NPY_2_COMPLEXCOMPAT_H_
+#define NUMPY_CORE_INCLUDE_NUMPY_NPY_2_COMPLEXCOMPAT_H_
+
+#include <numpy/npy_math.h>
+
+#ifndef NPY_CSETREALF
+#define NPY_CSETREALF(c, r) (c)->real = (r)
+#endif
+#ifndef NPY_CSETIMAGF
+#define NPY_CSETIMAGF(c, i) (c)->imag = (i)
+#endif
+#ifndef NPY_CSETREAL
+#define NPY_CSETREAL(c, r)  (c)->real = (r)
+#endif
+#ifndef NPY_CSETIMAG
+#define NPY_CSETIMAG(c, i)  (c)->imag = (i)
+#endif
+#ifndef NPY_CSETREALL
+#define NPY_CSETREALL(c, r) (c)->real = (r)
+#endif
+#ifndef NPY_CSETIMAGL
+#define NPY_CSETIMAGL(c, i) (c)->imag = (i)
+#endif
+
+#endif

numpy/core/include/numpy/npy_common.h

@@ -333,9 +333,11 @@ typedef unsigned char npy_bool;
  */
 #if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
     #define NPY_LONGDOUBLE_FMT "g"
+    #define longdouble_t double
     typedef double npy_longdouble;
 #else
     #define NPY_LONGDOUBLE_FMT "Lg"
+    #define longdouble_t long double
     typedef long double npy_longdouble;
 #endif
 
@@ -361,49 +363,56 @@ typedef double npy_double;
 typedef Py_hash_t npy_hash_t;
 #define NPY_SIZEOF_HASH_T NPY_SIZEOF_INTP
 
-/*
- * Disabling C99 complex usage: a lot of C code in numpy/scipy rely on being
- * able to do .real/.imag. Will have to convert code first.
- */
-#if 0
-#if defined(NPY_USE_C99_COMPLEX) && defined(NPY_HAVE_COMPLEX_DOUBLE)
-typedef complex npy_cdouble;
-#else
-typedef struct { double real, imag; } npy_cdouble;
-#endif
-
-#if defined(NPY_USE_C99_COMPLEX) && defined(NPY_HAVE_COMPLEX_FLOAT)
-typedef complex float npy_cfloat;
-#else
-typedef struct { float real, imag; } npy_cfloat;
-#endif
-
-#if defined(NPY_USE_C99_COMPLEX) && defined(NPY_HAVE_COMPLEX_LONG_DOUBLE)
-typedef complex long double npy_clongdouble;
-#else
-typedef struct {npy_longdouble real, imag;} npy_clongdouble;
-#endif
-#endif
 #if NPY_SIZEOF_COMPLEX_DOUBLE != 2 * NPY_SIZEOF_DOUBLE
 #error npy_cdouble definition is not compatible with C99 complex definition ! \
         Please contact NumPy maintainers and give detailed information about your \
         compiler and platform
 #endif
-typedef struct { double real, imag; } npy_cdouble;
-
 #if NPY_SIZEOF_COMPLEX_FLOAT != 2 * NPY_SIZEOF_FLOAT
 #error npy_cfloat definition is not compatible with C99 complex definition ! \
         Please contact NumPy maintainers and give detailed information about your \
         compiler and platform
 #endif
-typedef struct { float real, imag; } npy_cfloat;
-
 #if NPY_SIZEOF_COMPLEX_LONGDOUBLE != 2 * NPY_SIZEOF_LONGDOUBLE
 #error npy_clongdouble definition is not compatible with C99 complex definition ! \
         Please contact NumPy maintainers and give detailed information about your \
         compiler and platform
 #endif
-typedef struct { npy_longdouble real, imag; } npy_clongdouble;
+
+
+#if defined(_MSC_VER) && !defined(__INTEL_COMPILER) && defined(__cplusplus)
+typedef struct 
+{
+    double _Val[2];
+} npy_cdouble;
+
+typedef struct
+{
+    float _Val[2];
+} npy_cfloat;
+
+typedef struct
+{
+    long double _Val[2];
+} npy_clongdouble;
+#elif defined(_MSC_VER) && !defined(__INTEL_COMPILER) /* && !defined(__cplusplus) */
+#include <complex.h>
+typedef _Dcomplex npy_cdouble;
+typedef _Fcomplex npy_cfloat;
+typedef _Lcomplex npy_clongdouble;
+#else /* !defined(_MSC_VER) || defined(__INTEL_COMPILER) */
+#ifdef __cplusplus
+extern "C++" {
+#endif
+#include <complex.h>
+#ifdef __cplusplus
+}
+#endif
+#undef complex
+typedef double _Complex npy_cdouble;
+typedef float _Complex npy_cfloat;
+typedef longdouble_t _Complex npy_clongdouble;
+#endif
 
 /*
  * numarray-style bit-width typedefs

numpy/core/include/numpy/npy_math.h

@@ -360,84 +360,96 @@ NPY_INPLACE npy_longdouble npy_heavisidel(npy_longdouble x, npy_longdouble h0);
  * Complex declarations
  */
 
-/*
- * C99 specifies that complex numbers have the same representation as
- * an array of two elements, where the first element is the real part
- * and the second element is the imaginary part.
- */
-#define __NPY_CPACK_IMP(x, y, type, ctype)   \
-    union {                                  \
-        ctype z;                             \
-        type a[2];                           \
-    } z1;                                    \
-                                             \
-    z1.a[0] = (x);                           \
-    z1.a[1] = (y);                           \
-                                             \
-    return z1.z;
+static inline double npy_creal(const npy_cdouble z)
+{
+    return ((double *) &z)[0];
+}
 
-static inline npy_cdouble npy_cpack(double x, double y)
+static inline void npy_csetreal(npy_cdouble *z, const double r)
 {
-    __NPY_CPACK_IMP(x, y, double, npy_cdouble);
+    ((double *) z)[0] = r;
 }
 
-static inline npy_cfloat npy_cpackf(float x, float y)
+static inline double npy_cimag(const npy_cdouble z)
 {
-    __NPY_CPACK_IMP(x, y, float, npy_cfloat);
+    return ((double *) &z)[1];
 }
 
-static inline npy_clongdouble npy_cpackl(npy_longdouble x, npy_longdouble y)
+static inline void npy_csetimag(npy_cdouble *z, const double i)
 {
-    __NPY_CPACK_IMP(x, y, npy_longdouble, npy_clongdouble);
+    ((double *) z)[1] = i;
 }
-#undef __NPY_CPACK_IMP
 
-/*
- * Same remark as above, but in the other direction: extract first/second
- * member of complex number, assuming a C99-compatible representation
- *
- * Those are defineds as static inline, and such as a reasonable compiler would
- * most likely compile this to one or two instructions (on CISC at least)
- */
-#define __NPY_CEXTRACT_IMP(z, index, type, ctype)   \
-    union {                                         \
-        ctype z;                                    \
-        type a[2];                                  \
-    } __z_repr;                                     \
-    __z_repr.z = z;                                 \
-                                                    \
-    return __z_repr.a[index];
-
-static inline double npy_creal(npy_cdouble z)
+static inline float npy_crealf(const npy_cfloat z)
 {
-    __NPY_CEXTRACT_IMP(z, 0, double, npy_cdouble);
+    return ((float *) &z)[0];
 }
 
-static inline double npy_cimag(npy_cdouble z)
+static inline void npy_csetrealf(npy_cfloat *z, const float r)
 {
-    __NPY_CEXTRACT_IMP(z, 1, double, npy_cdouble);
+    ((float *) z)[0] = r;
 }
 
-static inline float npy_crealf(npy_cfloat z)
+static inline float npy_cimagf(const npy_cfloat z)
 {
-    __NPY_CEXTRACT_IMP(z, 0, float, npy_cfloat);
+    return ((float *) &z)[1];
 }
 
-static inline float npy_cimagf(npy_cfloat z)
+static inline void npy_csetimagf(npy_cfloat *z, const float i)
 {
-    __NPY_CEXTRACT_IMP(z, 1, float, npy_cfloat);
+    ((float *) z)[1] = i;
 }
 
-static inline npy_longdouble npy_creall(npy_clongdouble z)
+static inline npy_longdouble npy_creall(const npy_clongdouble z)
 {
-    __NPY_CEXTRACT_IMP(z, 0, npy_longdouble, npy_clongdouble);
+    return ((longdouble_t *) &z)[0];
 }
 
-static inline npy_longdouble npy_cimagl(npy_clongdouble z)
+static inline void npy_csetreall(npy_clongdouble *z, const longdouble_t r)
+{
+    ((longdouble_t *) z)[0] = r;
+}
+
+static inline npy_longdouble npy_cimagl(const npy_clongdouble z)
+{
+    return ((longdouble_t *) &z)[1];
+}
+
+static inline void npy_csetimagl(npy_clongdouble *z, const longdouble_t i)
+{
+    ((longdouble_t *) z)[1] = i;
+}
+
+#define NPY_CSETREAL(z, r) npy_csetreal(z, r)
+#define NPY_CSETIMAG(z, i) npy_csetimag(z, i)
+#define NPY_CSETREALF(z, r) npy_csetrealf(z, r)
+#define NPY_CSETIMAGF(z, i) npy_csetimagf(z, i)
+#define NPY_CSETREALL(z, r) npy_csetreall(z, r)
+#define NPY_CSETIMAGL(z, i) npy_csetimagl(z, i)
+
+static inline npy_cdouble npy_cpack(double x, double y)
+{
+    npy_cdouble z;
+    npy_csetreal(&z, x);
+    npy_csetimag(&z, y);
+    return z;
+}
+
+static inline npy_cfloat npy_cpackf(float x, float y)
+{
+    npy_cfloat z;
+    npy_csetrealf(&z, x);
+    npy_csetimagf(&z, y);
+    return z;
+}
+
+static inline npy_clongdouble npy_cpackl(npy_longdouble x, npy_longdouble y)
 {
-    __NPY_CEXTRACT_IMP(z, 1, npy_longdouble, npy_clongdouble);
+    npy_clongdouble z;
+    npy_csetreall(&z, x);
+    npy_csetimagl(&z, y);
+    return z;
 }
-#undef __NPY_CEXTRACT_IMP
 
 /*
  * Double precision complex functions

numpy/core/src/common/cblasfuncs.c

@@ -9,6 +9,7 @@
 #include <Python.h>
 
 #include "numpy/arrayobject.h"
+#include "numpy/npy_math.h"
 #include "npy_cblas.h"
 #include "arraytypes.h"
 #include "common.h"
@@ -422,8 +423,10 @@ cblas_matrixproduct(int typenum, PyArrayObject *ap1, PyArrayObject *ap2,
                 ptr1 = (npy_cdouble *)PyArray_DATA(ap2);
                 ptr2 = (npy_cdouble *)PyArray_DATA(ap1);
                 res = (npy_cdouble *)PyArray_DATA(out_buf);
-                res->real = ptr1->real * ptr2->real - ptr1->imag * ptr2->imag;
-                res->imag = ptr1->real * ptr2->imag + ptr1->imag * ptr2->real;
+                npy_csetreal(res, npy_creal(*ptr1) * npy_creal(*ptr2)
+                                            - npy_cimag(*ptr1) * npy_cimag(*ptr2));
+                npy_csetimag(res, npy_creal(*ptr1) * npy_cimag(*ptr2)
+                                            + npy_cimag(*ptr1) * npy_creal(*ptr2));
             }
             else if (ap1shape != _matrix) {
                 CBLAS_FUNC(cblas_zaxpy)(l,
@@ -495,8 +498,10 @@ cblas_matrixproduct(int typenum, PyArrayObject *ap1, PyArrayObject *ap2,
                 ptr1 = (npy_cfloat *)PyArray_DATA(ap2);
                 ptr2 = (npy_cfloat *)PyArray_DATA(ap1);
                 res = (npy_cfloat *)PyArray_DATA(out_buf);
-                res->real = ptr1->real * ptr2->real - ptr1->imag * ptr2->imag;
-                res->imag = ptr1->real * ptr2->imag + ptr1->imag * ptr2->real;
+                npy_csetrealf(res, npy_crealf(*ptr1) * npy_crealf(*ptr2)
+                                           - npy_cimagf(*ptr1) * npy_cimagf(*ptr2));
+                npy_csetimagf(res, npy_crealf(*ptr1) * npy_cimagf(*ptr2)
+                                           + npy_cimagf(*ptr1) * npy_crealf(*ptr2));
             }
             else if (ap1shape != _matrix) {
                 CBLAS_FUNC(cblas_caxpy)(l,