python - resumen - Números complejos en Cython

raiz negativa python

Question

¿Cuál es la forma correcta de trabajar con números complejos en Cython?

Me gustaría escribir un ciclo C puro usando un numpy.ndarray de dtype np.complex128. En Cython, el tipo de C asociado se define en Cython/Includes/numpy/__init__.pxd como

ctypedef double complex complex128_t

así que parece que esto es solo un simple complejo doble de C.

Sin embargo, es fácil obtener comportamientos extraños. En particular, con estas definiciones

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.

la línea

varc128 = varc128 * varf64

puede compilarse por Cython, pero gcc no puede compilar el código C producido (el error es "testcplx.c: 663: 25: error: dos o más tipos de datos en los especificadores de declaraciones" y parece deberse a la línea typedef npy_float64 _Complex __pyx_t_npy_float64_complex; ). Este error ya ha sido reportado (por ejemplo here ) pero no encontré ninguna buena explicación y / o solución limpia.

Sin inclusión de complex.h , no hay ningún error (supongo que porque typedef no está incluido).

Sin embargo, todavía hay un problema ya que en el archivo html producido por cython -a testcplx.pyx , la línea varc128 = varc128 * varf64 es amarilla, lo que significa que no se ha traducido a C. puro. El código C correspondiente es:

__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));

y __Pyx_CREAL y __Pyx_CIMAG son de color naranja (llamadas de Python).

Curiosamente, la línea

vardc = vardc * vard

no produce ningún error y se traduce en C puro (solo __pyx_v_8testcplx_vardc = __Pyx_c_prod(__pyx_v_8testcplx_vardc, __pyx_t_double_complex_from_parts(__pyx_v_8testcplx_vard, 0)); ), mientras que es muy similar al primero.

Puedo evitar el error usando variables intermedias (y se traduce en C pura):

vardc = varc128
vard = varf64
varc128 = vardc * vard

o simplemente lanzando (pero no se traduce en C puro):

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

¿Así que lo que ocurre? ¿Cuál es el significado del error de compilación? ¿Hay alguna forma limpia de evitarlo? ¿Por qué la multiplicación de np.complex128_t y np.float64_t parece implicar llamadas de Python?

Versiones

Cython versión 0.22 (versión más reciente en Pypi cuando se hizo la pregunta) y GCC 4.9.2.

Repositorio

Creé un pequeño repositorio con el ejemplo ( hg clone https://bitbucket.org/paugier/test_cython_complex ) y un pequeño Makefile con 3 objetivos ( make clean , make build , make html ) para que sea fácil probar cualquier cosa.

Answer 1

La forma más sencilla que puedo encontrar para solucionar este problema es simplemente cambiar el orden de multiplicación.

Si en testcplx.pyx cambio

varc128 = varc128 * varf64

a

varc128 = varf64 * varc128

Cambio de la situación de falla a la descrita a una que funciona correctamente. Este escenario es útil ya que permite una diferencia directa del código C producido.

tl; dr

El orden de la multiplicación cambia la traducción, lo que significa que en la versión que falla la multiplicación se intenta a través de los tipos __pyx_t_npy_float64_complex , mientras que en la versión de trabajo se realiza a través de los tipos __pyx_t_double_complex . Esto a su vez introduce la línea typedef npy_float64 _Complex __pyx_t_npy_float64_complex; , que no es válido

Estoy bastante seguro de que esto es un error de cython (Actualización: informada aquí ). Aunque este es un informe de error de gcc muy antiguo , la respuesta establece explícitamente (al decir que no es, de hecho, un error de gcc , sino un error de código de usuario):

typedef R _Complex C;

Este no es un código válido; no puede usar _Complex junto con un typedef, solo junto con "float", "double" o "long double" en una de las formas listadas en C99.

Llegan a la conclusión de que double _Complex es un especificador de tipo válido, mientras que ArbitraryType _Complex no lo es. Este informe más reciente tiene el mismo tipo de respuesta: intentar usar _Complex en un tipo no fundamental está fuera de la especificación, y el manual de GCC indica que _Complex solo se puede usar con float , double y long double

Por lo tanto, podemos piratear el código C generado por cython para probarlo: reemplace typedef npy_float64 _Complex __pyx_t_npy_float64_complex; con typedef double _Complex __pyx_t_npy_float64_complex; y verifique que sea efectivamente válido y que pueda compilar el código de salida.

Breve recorrido por el código

Al intercambiar la orden de multiplicación, solo se resalta el problema del que nos informa el compilador. En el primer caso, la línea ofensiva es la que dice typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - Está intentando asignar el tipo npy_float64 y usar la palabra clave _Complex al tipo __pyx_t_npy_float64_complex .

float _Complex o double _Complex es un tipo válido, mientras que npy_float64 _Complex no lo es. Para ver el efecto, puede simplemente eliminar npy_float64 de esa línea, o reemplazarlo por double o float y el código compila bien. La siguiente pregunta es por qué esa línea se produce en primer lugar ...

Esto parece ser producido por esta línea en el código fuente de Cython.

¿Por qué el orden de la multiplicación cambia el código de manera significativa, de modo que se introduce el tipo __pyx_t_npy_float64_complex , y se introduce de una manera que falla?

En la instancia varf64 , el código para implementar la multiplicación convierte varf64 en un tipo __pyx_t_npy_float64_complex , realiza la multiplicación en partes reales e imaginarias y luego vuelve a ensamblar el número complejo. En la versión de trabajo, hace el producto directamente a través del tipo __pyx_t_double_complex con la función __Pyx_c_prod

Supongo que esto es tan simple como que el código de cython tome su señal para qué tipo usar para la multiplicación de la primera variable que encuentra. En el primer caso, ve un float 64, por lo que produce un código C ( inválido ) basado en eso, mientras que en el segundo, ve el tipo (double) complex128 y basa su traducción en eso. Esta explicación es un poco ondulada y espero volver a analizarla si el tiempo lo permite ...

Una nota al respecto: aquí vemos que typedef para npy_float64 es double , por lo que en este caso particular, una corrección podría consistir en modificar el código aquí para usar double _Complex donde type es npy_float64 , pero esto va más allá del alcance de un SO responder y no presenta una solución general.

Resultado del código de C diff

Versión de trabajo

Crea este código C desde la línea `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);

Versión fallida

Crea este código C desde la línea 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));

Lo que necesita estas importaciones adicionales, y la línea ofensiva es la que dice typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - Está intentando asignar el tipo npy_float64 y el tipo _Complex al tipo __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