Was ist der richtige Weg, um mit komplexen Zahlen in Cython zu arbeiten?

Ich möchte eine reine C-Schleife mit einem numpy.ndarray von dtype np.complex128 schreiben. In Cython ist der zugehörige C-Typ in definiert
Cython/Includes/numpy/__init__.pxd wie

ctypedef double complex complex128_t

Es scheint also, dass dies nur ein einfacher C-Doppelkomplex ist.

Es ist jedoch leicht, seltsame Verhaltensweisen zu erhalten. Insbesondere mit diesen Definitionen

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.

die Linie

varc128 = varc128 * varf64

kann von Cython kompiliert werden, aber gcc kann den erzeugten C-Code nicht kompilieren (der Fehler lautet “testcplx.c:663:25: Fehler: Zwei oder mehr Datentypen in Deklarationsbezeichnern” und scheint auf die Zeile zurückzuführen zu sein typedef npy_float64 _Complex __pyx_t_npy_float64_complex;). Dieser Fehler wurde bereits gemeldet (z hier), aber ich habe keine gute Erklärung und/oder saubere Lösung gefunden.

Ohne Einbeziehung von complex.hes gibt keinen Fehler (ich denke, weil die typedef ist dann nicht enthalten).

Es gibt jedoch immer noch ein Problem, da in der von erzeugten HTML-Datei cython -a testcplx.pyxdie Linie varc128 = varc128 * varf64 ist gelb, was bedeutet, dass es nicht in reines C übersetzt wurde. Der entsprechende C-Code lautet:

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

und die __Pyx_CREAL und __Pyx_CIMAG sind orange (Python-Aufrufe).

Interessanterweise die Linie

vardc = vardc * vard

erzeugt keinen Fehler und wird in reines C übersetzt (nur __pyx_v_8testcplx_vardc = __Pyx_c_prod(__pyx_v_8testcplx_vardc, __pyx_t_double_complex_from_parts(__pyx_v_8testcplx_vard, 0));), wobei es dem ersten sehr, sehr ähnlich ist.

Ich kann den Fehler vermeiden, indem ich Zwischenvariablen verwende (und es in reines C übersetzt):

vardc = varc128
vard = varf64
varc128 = vardc * vard

oder einfach durch Casting (aber es übersetzt sich nicht in reines C):

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

Was passiert also? Was bedeutet der Kompilierungsfehler? Gibt es eine saubere Möglichkeit, das zu vermeiden? Warum scheint die Multiplikation von np.complex128_t und np.float64_t Python-Aufrufe zu beinhalten?

Versionen

Cython Version 0.22 (neueste Version in Pypi, als die Frage gestellt wurde) und GCC 4.9.2.

Repository

Ich habe ein kleines Repository mit dem Beispiel erstellt (hg clone https://bitbucket.org/paugier/test_cython_complex) und ein winziges Makefile mit 3 Zielen (make clean, make build, make html) so ist es einfach, alles zu testen.

Die einfachste Möglichkeit, dieses Problem zu umgehen, besteht darin, einfach die Reihenfolge der Multiplikation zu ändern.

Wenn drin testcplx.pyx Ich ändere

varc128 = varc128 * varf64

zu

varc128 = varf64 * varc128

Ich wechsle von der fehlerhaften Situation zur beschriebenen zu einer korrekt funktionierenden Situation. Dieses Szenario ist nützlich, da es einen direkten Vergleich des produzierten C-Codes ermöglicht.

tl;dr

Die Reihenfolge der Multiplikation ändert die Übersetzung, was bedeutet, dass in der fehlgeschlagenen Version die Multiplikation über versucht wird __pyx_t_npy_float64_complex Typen, während dies in der Arbeitsversion über erfolgt __pyx_t_double_complex Typen. Dies wiederum leitet die typedef-Zeile ein typedef npy_float64 _Complex __pyx_t_npy_float64_complex;was ungültig ist.

Ich bin mir ziemlich sicher, dass dies ein Cython-Bug ist (Update: hier berichtet). Obwohl Dies ist ein sehr alter gcc-Fehlerberichtheißt es in der Antwort ausdrücklich (indem sie sagt, dass es sich tatsächlich nicht um a gcc Fehler, aber Benutzercodefehler):

typedef R _Complex C;

Dies ist kein gültiger Code; Sie können _Complex nicht zusammen mit einem Typedef verwenden, sondern nur zusammen mit “float”, “double” oder “long double” in einer der in C99 aufgeführten Formen.

Das schließen sie double _Complex ist ein gültiger Typbezeichner, wohingegen ArbitraryType _Complex ist nicht. Dieser neuere Bericht hat die gleiche Art von Antwort – versucht zu verwenden _Complex auf einem nicht fundamentalen Typ ist außerhalb der Spezifikation, und die GCC-Handbuch zeigt an, dass _Complex kann nur mit verwendet werden float, double und long double

Also – wir können den von Cython generierten C-Code hacken, um das zu testen: Ersetzen typedef npy_float64 _Complex __pyx_t_npy_float64_complex; mit typedef double _Complex __pyx_t_npy_float64_complex; und vergewissern Sie sich, dass es tatsächlich gültig ist und den Ausgabecode kompilieren kann.

Kurze Wanderung durch den Code

Das Vertauschen der Multiplikationsreihenfolge hebt nur das Problem hervor, von dem uns der Compiler berichtet. Im ersten Fall ist die beleidigende Zeile diejenige, die sagt typedef npy_float64 _Complex __pyx_t_npy_float64_complex; – Es wird versucht, den Typ zuzuweisen npy_float64 und benutze das Schlüsselwort _Complex zum Typ __pyx_t_npy_float64_complex.

float _Complex oder double _Complex ist ein gültiger Typ, wohingegen npy_float64 _Complex ist nicht. Um den Effekt zu sehen, können Sie einfach löschen npy_float64 aus dieser Zeile, oder ersetzen Sie es durch double oder float und der Code kompiliert gut. Die nächste Frage ist, warum diese Linie überhaupt hergestellt wird …

Dies scheint von produziert zu sein diese Linie im Cython-Quellcode.

Warum ändert die Reihenfolge der Multiplikation den Code erheblich – so dass der Typ __pyx_t_npy_float64_complex eingeführt wird und auf eine Weise eingeführt wird, die fehlschlägt?

Im fehlgeschlagenen Fall dreht sich der Code zum Implementieren der Multiplikation varf64 in ein __pyx_t_npy_float64_complex Typ, führt die Multiplikation mit Real- und Imaginärteilen durch und setzt dann die komplexe Zahl wieder zusammen. In der Arbeitsversion macht es das Produkt direkt über die __pyx_t_double_complex mit der Funktion eingeben __Pyx_c_prod

Ich denke, das ist so einfach wie der Cython-Code, der anhand der ersten Variablen, auf die er stößt, darauf hinweist, welcher Typ für die Multiplikation verwendet werden soll. Im ersten Fall sieht es einen Float 64, erzeugt also (ungültig) C-Code basiert darauf, während es im zweiten den (doppelten) komplexen Typ128 sieht und seine Übersetzung darauf basiert. Diese Erklärung ist ein wenig handgewellt und ich hoffe, dass ich zu einer Analyse davon zurückkehren kann, wenn es die Zeit erlaubt …

Eine Anmerkung dazu – hier sehen wir dass die typedef zum npy_float64 ist doubleso dass in diesem speziellen Fall eine Lösung aus einer Änderung bestehen könnte der Code hier benutzen double _Complex wo type ist npy_float64aber dies geht über den Rahmen einer SO-Antwort hinaus und stellt keine allgemeine Lösung dar.

C-Code-Diff-Ergebnis

Arbeitsversion

Erstellt diesen C-Code aus der Zeile `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);

Fehlerhafte Version

Erstellt diesen C-Code aus der Zeile 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));

Was diese zusätzlichen Importe erfordert – und die beleidigende Zeile ist diejenige, die besagt typedef npy_float64 _Complex __pyx_t_npy_float64_complex; – Es wird versucht, den Typ zuzuweisen npy_float64 und der Typ _Complex zum Typ __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