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

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

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

 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 версии 0.22 (последняя версия в Pypi, когда был задан вопрос) и GCC 4.9.2.

вместилище

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

One Solution collect form web for “Комплексные числа в Китоне”

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

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

 varc128 = varc128 * varf64 

в

 varc128 = varf64 * varc128 

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

ТЛ; др

Порядок умножения изменяет перевод, а это означает, что в неудачной версии попытка умножается с помощью типов __pyx_t_npy_float64_complex , тогда как в рабочей версии это делается с помощью типов __pyx_t_double_complex . Это, в свою очередь, вводит строку 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 тип __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 где typenpy_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 
  • Как добавить столбец в массив numpy
  • заменить значения в массиве
  • нанесение маркеров на топоры
  • Python numpy.square vs **
  • многократные формы умножения матрицы
  • Округление до значительных цифр в numpy
  • Кластеризация ~ 100 000 коротких строк в Python
  • Список индексирования массива python в списке
  • Python - лучший язык программирования в мире.