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?

DavidW · Accepted Answer · 2017-03-28T12:06:18.367

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.

Ah sorry - I didn't see your edit while I was writing this, I think this is telling you what you already know. It was probably useful to the question as originally written. — DavidW, Mar 28 '17 at 11:19
Really cryptic. I had to make sure it wont format my hard drive prior to testing ;) But it works! — Antony Hatchkins, Mar 28 '17 at 16:19
To be honest, such things frighten me a bit. I'm considering cffi as an alternative to cython. Did you try it? For example it can make dlls right out of the box (if you are a windows guy you'll understand), while cython only does half of the work (`cython --embed` produces the c/cpp code and you have to do the compiling and linking manually - I only found a tutorial for exe, not for dll) — Antony Hatchkins, Mar 28 '17 at 16:20
1) This is really a bug fix/improvement to Cython (so you shouldn't have to write such code yourself hopefully!). However, if you find yourself hitting bugs with Cython often then perhaps it's worth looking at something else. 2) I do like `cffi` in principle. As you say it provides a good way of calling C libraries without too much compiling. It will probably run quite a bit slower than the equivalent Cython code. This may or may not matter, but the question suggests you probably were trying to get a speed improvement (since the Python and C `abs` should give the same answer). — DavidW, Mar 28 '17 at 17:21
afaik cffi also has a use case of including pure c/c++ code into python which should run at the same speed as cython (basically, the same c/c++ code is generated). Back to this question. Why would cython inject that code with three versions of abs into every generated c/c++ file if there was no way of using it before this patch was applied? — Antony Hatchkins, Mar 28 '17 at 18:31
It seems to successfully use all the other pure C functions for complex numbers (e.g. arithmetic and `.conjugate`) so I'd guess it's just an omission. Converting back to a Python complex doesn't actually break anything so I doubt it ever got noticed. — DavidW, Mar 28 '17 at 18:43

abs(double complex) in Cython

1 Answers1

Linked