Комплексные числа в Китоне

Каков правильный способ работы с комплексными числами в Cython?

Я хотел бы написать чистый цикл C, используя numpy.ndarray из dtype np.complex128. В Cython связанный тип C определен в Cython/Includes/numpy/__init__.pxd как

ctypedef double complex complex128_t

поэтому кажется, что это всего лишь простой C-двойной комплекс.

Однако легко получить странное поведение. В частности, с этими определениями

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.

строка

varc128 = varc128 * varf64

может быть скомпилирован Cython, но gcc не может скомпилировать полученный код C (ошибка: "testcplx.c: 663: 25: ошибка: два или более типов данных в спецификаторах объявлений" и, похоже, связана с линией typedef npy_float64 _Complex __pyx_t_npy_float64_complex;). Эта ошибка уже сообщалась (например, здесь), но я не нашел никакого объяснения и/или чистого решения.

Без включения complex.h ошибки нет (думаю, потому что typedef тогда не включается).

Однако все еще существует проблема, поскольку в html файле, создаваемом cython -a testcplx.pyx, строка varc128 = varc128 * varf64 является желтой, что означает, что она не была переведена в чистый C. Соответствующий код C:

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

а __Pyx_CREAL и __Pyx_CIMAG - оранжевые (вызовы Python).

Интересно, что линия

vardc = vardc * vard

не вызывает ошибок и переводится в чистый C (просто __pyx_v_8testcplx_vardc = __Pyx_c_prod(__pyx_v_8testcplx_vardc, __pyx_t_double_complex_from_parts(__pyx_v_8testcplx_vard, 0));), тогда как он очень похож на первый.

Я могу избежать ошибки, используя промежуточные переменные (и это переводит в чистый C):

vardc = varc128
vard = varf64
varc128 = vardc * vard

или просто путем литья (но он не преобразуется в чистый C):

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

Так что же происходит? В чем смысл ошибки компиляции? Есть ли чистый способ избежать этого? Почему умножение np.complex128_t и np.float64_t, по-видимому, связано с вызовами Python?

Версия

Cython version 0.22 (самая последняя версия в Pypi, когда был задан вопрос) и GCC 4.9.2.

Repository

Я создал крошечный репозиторий с примером (hg clone https://bitbucket.org/paugier/test_cython_complex) и крошечный Makefile с тремя целями (make clean, make build, make html), поэтому легко проверить что-либо.

Ответы

Ответ 1

Самый простой способ найти эту проблему - просто переключить порядок умножения.

Если в testcplx.pyx я меняю

varc128 = varc128 * varf64

к

varc128 = varf64 * varc128

Я перехожу из ситуации сбойной ситуации к описательной, которая работает правильно. Этот сценарий полезен, поскольку он позволяет напрямую различать полученный код C.

TL;DR

Порядок умножения изменяет перевод, а это означает, что в неудачной версии попытка умножается на типы __pyx_t_npy_float64_complex, тогда как в рабочей версии это делается с помощью типов __pyx_t_double_complex. Это, в свою очередь, вводит строку typedef typedef npy_float64 _Complex __pyx_t_npy_float64_complex;, что является недопустимым.

Я уверен, что это ошибка cython (Update: здесь). Хотя это очень старый отчет об ошибке gcc, ответ явно заявляет (говоря, что это, по сути, ошибка gcc но ошибка кода пользователя):

typedef R _Complex C;

Это недопустимый код; вы не можете использовать _Complex вместе с typedef, только вместе с "float", "double" или "long double" в одной из форм перечисленных на C99.

Они заключают, что double _Complex является допустимым спецификатором типа, тогда как ArbitraryType _Complex не является. Этот более свежий отчет имеет тот же тип ответа - попытка использовать _Complex для не фундаментального типа - вне спецификации, а GCC указывает, что _Complex может использоваться только с float, double и long double

Итак - мы можем взломать C-код, сгенерированный на Cython, чтобы проверить: замените typedef npy_float64 _Complex __pyx_t_npy_float64_complex; на typedef double _Complex __pyx_t_npy_float64_complex; и убедитесь, что он действительно действителен и может скомпилировать выходной код.

Короткий переход через код

Обмен порядком умножения только подчеркивает проблему, о которой нам говорит компилятор. В первом случае строка с нарушением - это typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - она ​​пытается присвоить типу npy_float64 и ключевое слово _Complex для типа __pyx_t_npy_float64_complex.

float _Complex или double _Complex является допустимым типом, тогда как npy_float64 _Complex - нет. Чтобы увидеть эффект, вы можете просто удалить npy_float64 из этой строки или заменить его на double или float и компилировать код отлично. Следующий вопрос: почему эта линия создается в первую очередь...

Это похоже на эту строку в исходном коде Cython.

Почему порядок умножения значительно меняет код - такой, что вводится тип __pyx_t_npy_float64_complex и вводится таким образом, что это не удается?

В неудавшемся экземпляре код для реализации умножения превращает varf64 в тип __pyx_t_npy_float64_complex, выполняет ли умножение на действительную и мнимую части, а затем снова собирает комплексное число. В рабочей версии продукт непосредственно через тип __pyx_t_double_complex использует функцию __Pyx_c_prod

Я думаю, это так же просто, как и код cython, который берет свой сигнал, для какого типа использовать для умножения от первой переменной, с которой он сталкивается. В первом случае он видит float 64, поэтому на основе этого генерирует (недействительный) код C, тогда как во втором он видит (double) complex128 и основывает на нем свой перевод. Это объяснение немного ручно-волнистое, и я надеюсь вернуться к его анализу, если позволит время...

Заметка об этом - здесь мы видим, что typedef для npy_float64 - double, поэтому в данном конкретном случае исправление может состоять из изменяя код здесь, чтобы использовать double _Complex, где type есть npy_float64, но это выходит за рамки ответа SO и не представляет общий решение.

Результат кода кода C

Рабочая версия

Создает этот код C из строки `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);

Ошибка версии

Создает этот код C из строки 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));

Это требует этих дополнительных импортов, а строка с нарушением - typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - она ​​пытается присвоить типу npy_float64 и тип _Complex типу __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