Complex numbers in Cython

What is the correct way to work with complex numbers in Cython?

I would like to write a pure C loop using a numpy.ndarray of dtype np.complex128. In Cython, the associated C type is defined in
``Cython/Includes/numpy/__init__.pxd`` as

ctypedef double complex complex128_t

so it seems this is just a simple C double complex.

However, it's easy to obtain strange behaviors. In particular, with these definitions

cimport numpy as np
import numpy as np
np.import_array()

cdef extern from "complex.h":
pass

cdef:
np.complex128_t varc128 = 1j
np.float64_t varf64 = 1.
double complex vardc = 1j
double vard = 1.

the line

varc128 = varc128 * varf64

can be compiled by Cython but gcc can not compiled the C code produced (the error is "testcplx.c:663:25: error: two or more data types in declaration specifiers" and seems to be due to the line ``typedef npy_float64 _Complex __pyx_t_npy_float64_complex;``). This error has already been reported (for example [here](

[To see links please register here]

)) but I didn't find any good explanation and/or clean solution.

Without inclusion of ``complex.h``, there is no error (I guess because the ``typedef`` is then not included).

However, there is still a problem since in the html file produced by ``cython -a testcplx.pyx``, the line ``varc128 = varc128 * varf64`` is yellow, meaning that it has not been translated into pure C. The corresponding C code is:

__pyx_t_2 = __Pyx_c_prod_npy_float64(__pyx_t_npy_float64_complex_from_parts(__Pyx_CREAL(__pyx_v_8testcplx_varc128), __Pyx_CIMAG(__pyx_v_8testcplx_varc128)), __pyx_t_npy_float64_complex_from_parts(__pyx_v_8testcplx_varf64, 0));
__pyx_v_8testcplx_varc128 = __pyx_t_double_complex_from_parts(__Pyx_CREAL(__pyx_t_2), __Pyx_CIMAG(__pyx_t_2));

and the ``__Pyx_CREAL`` and ``__Pyx_CIMAG`` are orange (Python calls).

Interestingly, the line

vardc = vardc * vard

does not produce any error and is translated into pure C (just ``__pyx_v_8testcplx_vardc = __Pyx_c_prod(__pyx_v_8testcplx_vardc, __pyx_t_double_complex_from_parts(__pyx_v_8testcplx_vard, 0));``), whereas it is very very similar to the first one.

I can avoid the error by using intermediate variables (and it translates into pure C):

vardc = varc128
vard = varf64
varc128 = vardc * vard

or simply by casting (but it does not translate into pure C):

vardc = <double complex>varc128 * <double>varf64

So what happens? What is the meaning of the compilation error? Is there a clean way to avoid it? Why does the multiplication of a np.complex128_t and np.float64_t seem to involve Python calls?

Versions
--------

Cython version 0.22 (most recent version in Pypi when the question was asked) and GCC 4.9.2.

Repository
----------

I created a tiny repository with the example (`hg clone

[To see links please register here]

) and a tiny Makefile with 3 targets (`make clean`, `make build`, `make html`) so it is easy to test anything.

The simplest way I can find to work around this issue is to simply switch the order of multiplication.

If in `testcplx.pyx` I change

varc128 = varc128 * varf64

to

varc128 = varf64 * varc128

I change from the failing situation to described to one that works correctly. This scenario is useful as it allows a direct diff of the produced C code.

###tl;dr

The order of the multiplication changes the translation, meaning that in the failing version the multiplication is attempted via `__pyx_t_npy_float64_complex` types, whereas in the working version it is done via `__pyx_t_double_complex` types. This in turn introduces the typedef line `typedef npy_float64 _Complex __pyx_t_npy_float64_complex;`, which is invalid.

I am fairly sure this is a cython bug (Update: [reported here][1]). Although [this is a very old gcc bug report][2], the response explicitly states (in saying that it is not, in fact, a ***gcc*** bug, but user code error):

> typedef R _Complex C;

> This is not valid code; you can't use _Complex together with a typedef,
only together with "float", "double" or "long double" in one of the forms
listed in C99.

They conclude that `double _Complex` is a valid type specifier whereas `ArbitraryType _Complex` is not. [This more recent report][3] has the same type of response - trying to use `_Complex` on a non fundamental type is outside spec, and the [GCC manual][4] indicates that `_Complex` can only be used with `float`, `double` and `long double`

So - we can hack the cython generated C code to test that: replace `typedef npy_float64 _Complex __pyx_t_npy_float64_complex;` with `typedef double _Complex __pyx_t_npy_float64_complex;` and verify that it is indeed valid and can make the output code compile.

----------

###Short trek through the code

Swapping the multiplication order only highlights the problem that we are told about by the compiler. In the first case, the offending line is the one that says `typedef npy_float64 _Complex __pyx_t_npy_float64_complex;` - it is trying to assign the type `npy_float64` ***and*** use the keyword `_Complex` to the type `__pyx_t_npy_float64_complex`.

`float _Complex` or `double _Complex` is a valid type, whereas `npy_float64 _Complex` is not. To see the effect, you can just delete `npy_float64` from that line, or replace it with `double` or `float` and the code compiles fine. The next question is why that line is produced in the first place...

This seems to be produced by [this line][5] in the Cython source code.

Why does the order of the multiplication change the code significantly - such that the type `__pyx_t_npy_float64_complex` is introduced, and introduced in a way that fails?

In the failing instance, the code to implement the multiplication turns `varf64` into a `__pyx_t_npy_float64_complex` type, does the multiplication on real and imaginary parts and then reassembles the complex number. In the working version, it does the product directly via the `__pyx_t_double_complex` type using the function `__Pyx_c_prod`

I guess this is as simple as the cython code taking its cue for which type to use for the multiplication from the first variable it encounters. In the first case, it sees a float 64, so produces (*invalid*) C code based on that, whereas in the second, it sees the (double) complex128 type and bases its translation on that. This explanation is a little hand-wavy and I hope to return to an analysis of it if time allows...

A note on this - [here we see][6] that the `typedef` for `npy_float64` is `double`, so in this particular case, a fix might consist of modifying [the code here][7] to use `double _Complex` where `type` is `npy_float64`, but this is getting beyond the scope of a SO answer and doesn't present a general solution.

--------

#C code diff result

###Working version

Creates this C code from the line `varc128 = varf64 * varc128

__pyx_v_8testcplx_varc128 = __Pyx_c_prod(__pyx_t_double_complex_from_parts(__pyx_v_8testcplx_varf64, 0), __pyx_v_8testcplx_varc128);

### Failing version

Creates this C code from the line `varc128 = varc128 * varf64`

__pyx_t_2 = __Pyx_c_prod_npy_float64(__pyx_t_npy_float64_complex_from_parts(__Pyx_CREAL(__pyx_v_8testcplx_varc128), __Pyx_CIMAG(__pyx_v_8testcplx_varc128)), __pyx_t_npy_float64_complex_from_parts(__pyx_v_8testcplx_varf64, 0));
__pyx_v_8testcplx_varc128 = __pyx_t_double_complex_from_parts(__Pyx_CREAL(__pyx_t_2), __Pyx_CIMAG(__pyx_t_2));

Which necessitates these extra imports - and the offending line is the one that says `typedef npy_float64 _Complex __pyx_t_npy_float64_complex;` - it is trying to assign the type `npy_float64` ***and*** the type `_Complex` to the type `__pyx_t_npy_float64_complex`

#if CYTHON_CCOMPLEX
#ifdef __cplusplus
typedef ::std::complex< npy_float64 > __pyx_t_npy_float64_complex;
#else
typedef npy_float64 _Complex __pyx_t_npy_float64_complex;
#endif
#else
typedef struct { npy_float64 real, imag; } __pyx_t_npy_float64_complex;
#endif

/*... loads of other stuff the same ... */

static CYTHON_INLINE __pyx_t_npy_float64_complex __pyx_t_npy_float64_complex_from_parts(npy_float64, npy_float64);

#if CYTHON_CCOMPLEX
#define __Pyx_c_eq_npy_float64(a, b) ((a)==(b))
#define __Pyx_c_sum_npy_float64(a, b) ((a)+(b))
#define __Pyx_c_diff_npy_float64(a, b) ((a)-(b))
#define __Pyx_c_prod_npy_float64(a, b) ((a)*(b))
#define __Pyx_c_quot_npy_float64(a, b) ((a)/(b))
#define __Pyx_c_neg_npy_float64(a) (-(a))
#ifdef __cplusplus
#define __Pyx_c_is_zero_npy_float64(z) ((z)==(npy_float64)0)
#define __Pyx_c_conj_npy_float64(z) (::std::conj(z))
#if 1
#define __Pyx_c_abs_npy_float64(z) (::std::abs(z))
#define __Pyx_c_pow_npy_float64(a, b) (::std::pow(a, b))
#endif
#else
#define __Pyx_c_is_zero_npy_float64(z) ((z)==0)
#define __Pyx_c_conj_npy_float64(z) (conj_npy_float64(z))
#if 1
#define __Pyx_c_abs_npy_float64(z) (cabs_npy_float64(z))
#define __Pyx_c_pow_npy_float64(a, b) (cpow_npy_float64(a, b))
#endif
#endif
#else
static CYTHON_INLINE int __Pyx_c_eq_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_sum_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_diff_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_prod_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_quot_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_neg_npy_float64(__pyx_t_npy_float64_complex);
static CYTHON_INLINE int __Pyx_c_is_zero_npy_float64(__pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_conj_npy_float64(__pyx_t_npy_float64_complex);
#if 1
static CYTHON_INLINE npy_float64 __Pyx_c_abs_npy_float64(__pyx_t_npy_float64_complex);
static CYTHON_INLINE __pyx_t_npy_float64_complex __Pyx_c_pow_npy_float64(__pyx_t_npy_float64_complex, __pyx_t_npy_float64_complex);
#endif
#endif


[1]:

[To see links please register here]

[2]:

[To see links please register here]

[3]:

[To see links please register here]

[4]:

[To see links please register here]

[5]:

[To see links please register here]

[6]:

[To see links please register here]

[7]:

[To see links please register here]