abs(double complex) in Cython

Question

How do I get an absolute value for a double complex variable?

def f():
    cdef double complex aaa = 1 + 2j
    cdef double bbb = abs(aaa)

The second assignment is highlighted yellow in cython -a html output: aaa is converted to python object prior to applying abs().

How do I call abs() on c/cpp level?

PS I understand that

cdef extern from "complex.h":
    double abs(double complex)

would solve it, but I see the following instructions in the generated c/cpp file:

#if CYTHON_CCOMPLEX
        #define __Pyx_c_abs_double(z)     (::std::abs(z))
#endif

and the like, which is supposed to choose the correct header to include (<complex> or "complex.h" or custom code) depending on the compilation flags.

How do I utilize those instructions?

Answer 1

More useful contribution:

The following is semi-tested addition to fix "cython/compiler/Builtin.py". It should be pretty obvious where to add it:

BuiltinFunction('abs',        None,    None,   "__Pyx_c_abs{0}".format(PyrexTypes.c_double_complex_type.funcsuffix),
                #utility_code = UtilityCode.load('Arithmetic', 'Complex.c', PyrexTypes.c_double_complex_type._utility_code_context()),
                func_type = PyrexTypes.CFuncType(
                    PyrexTypes.c_double_type, [
                        PyrexTypes.CFuncTypeArg("arg", PyrexTypes.c_double_complex_type, None)
                        ],
                        is_strict_signature = True)),
BuiltinFunction('abs',        None,    None,   "__Pyx_c_abs{0}".format(PyrexTypes.c_float_complex_type.funcsuffix),
                #utility_code = UtilityCode.load('Arithmetic', 'Complex.c', PyrexTypes.c_float_complex_type._utility_code_context()),
                func_type = PyrexTypes.CFuncType(
                    PyrexTypes.c_float_type, [
                        PyrexTypes.CFuncTypeArg("arg", PyrexTypes.c_float_complex_type, None)
                        ],
                        is_strict_signature = True)),
BuiltinFunction('abs',        None,    None,   "__Pyx_c_abs{0}".format(PyrexTypes.c_longdouble_complex_type.funcsuffix),
                #utility_code = UtilityCode.load('Arithmetic', 'Complex.c', PyrexTypes.c_longdouble_complex_type._utility_code_context()),
                func_type = PyrexTypes.CFuncType(
                    PyrexTypes.c_longdouble_type, [
                        PyrexTypes.CFuncTypeArg("arg", PyrexTypes.c_longdouble_complex_type, None)
                        ],
                        is_strict_signature = True)),

It appears to work. It hasn't been run through the full Cython test-suite. It doesn't yet generate the correct code at the top of your file(but probably shouldn't need to since you already use complex double which does). It doesn't yet work in a nogil block.

I'll probably submit it to Cython github once I've looked at these issues.

---

Original Answer:

The whole thing is made slightly more complicated since Cython attempts to use the "native" layout for the complex type. From the code generated by Cython:

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

In all cases the types should have a compatible memory layout (I think), so the worst case is that a bit of type casting should let you use any implementation's abs function.

To add to the confusion, if you look at the generated Cython code (search through the various #if CYTHON_CCOMPLEX blocks in the generated file) it appears that Cython defines fast versions of abs (and other useful functions) for all these types, but fails to use them intelligently, falling back to the Python implementation.

If you're going through C you need to tell Cython about cabs from complex.h:

cdef extern from "complex.h":
    double cabs(double complex)

def f():
    cdef double complex aaa = 1 + 2j
    cdef double bbb = cabs(aaa)
    return bbb

If you're going through C++ you need to tell Cython about the C++ standard library abs. Unfortunately it's unable to make the link between the libcpp.complex.complex in its predefined wrapper and the double complex type, so you need to tell it about the function yourself:

cdef extern from "<complex>":
    double abs(double complex)

def f():
    cdef complex aaa = 1 + 2j
    cdef double bbb = abs(aaa)
    return bbb

I don't immediately know what to do if CYTHON_CCOMPLEX isn't defined, but I think that's something you'd have to do explicitly anyway.