Cython中的复数 [英] Complex numbers in Cython
问题描述
在Cython中处理复数的正确方法是什么?
我想使用dtype np.complex128的numpy.ndarray编写纯C循环.在Cython中，相关的C类型在
Cython/Includes/numpy/__init__.pxd为
ctypedef double complex complex128_t
所以这似乎只是一个简单的C double复数.
但是，很容易获得奇怪的行为.特别是这些定义
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).
但是，仍然存在问题，因为在cython -a testcplx.pyx生成的html文件中，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)创建了一个小型存储库，并创建了一个具有3个目标(make clean，make build，make html)的小型Makefile，因此可以轻松地进行任何测试.
找到该问题的最简单方法是简单地切换乘法顺序.
如果在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错误(更新:在此处报告).尽管这是一个非常古老的gcc错误报告，但响应明确指出(表示实际上不是 gcc 错误，而是用户代码错误):
typedef R _Complex C;这是无效的代码；您不能将_Complex与typedef一起使用， 仅与以下形式之一的"float"，"double"或"long double"一起使用 列在C99中.
他们得出的结论是double _Complex是有效的类型说明符，而ArbitraryType _Complex则不是. 此最新报告具有相同类型的响应-尝试在非基本类型超出了规范，并且 GCC手册表示_Complex只能与float，double和long double
所以-我们可以破解cython生成的C代码以进行测试:用typedef double _Complex __pyx_t_npy_float64_complex;替换typedef npy_float64 _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_c_prod
__pyx_t_double_complex类型直接生产产品.
我想这就像cython代码一样简单，它可以从遇到的第一个变量中得知要用于乘法的类型.在第一种情况下，它看到一个浮点数64，因此基于该值产生( invalid )C代码，而在第二种情况下，它看到了(double)complex128类型，并将其翻译基于该代码.这种解释有点费事，我希望在时间允许的情况下可以重新分析它.
对此的注释-在这里我们看到
从npy_float64的typedef是double，因此在这种特殊情况下，修复程序可能包括修改__pyx_v_8testcplx_varc128 = __Pyx_c_prod(__pyx_t_double_complex_from_parts(__pyx_v_8testcplx_varf64, 0), __pyx_v_8testcplx_varc128);
失败版本
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 和 c19>类型为__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
What is the correct way to work with complex numbers in Cython?
I would like to write a pure C loop using a numpy.ndarray of dtype np.complex128. In Cython, the associated C type is defined in
Cython/Includes/numpy/__init__.pxd as
ctypedef double complex complex128_t
so it seems this is just a simple C double complex.
However, it's easy to obtain strange behaviors. In particular, with these definitions
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.
the line
varc128 = varc128 * varf64
can be compiled by Cython but gcc can not compiled the C code produced (the error is "testcplx.c:663:25: error: two or more data types in declaration specifiers" and seems to be due to the line typedef npy_float64 _Complex __pyx_t_npy_float64_complex;). This error has already been reported (for example here) but I didn't find any good explanation and/or clean solution.
Without inclusion of complex.h, there is no error (I guess because the typedef is then not included).
However, there is still a problem since in the html file produced by cython -a testcplx.pyx, the line varc128 = varc128 * varf64 is yellow, meaning that it has not been translated into pure C. The corresponding C code is:
__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));
and the __Pyx_CREAL and __Pyx_CIMAG are orange (Python calls).
Interestingly, the line
vardc = vardc * vard
does not produce any error and is translated into pure C (just __pyx_v_8testcplx_vardc = __Pyx_c_prod(__pyx_v_8testcplx_vardc, __pyx_t_double_complex_from_parts(__pyx_v_8testcplx_vard, 0));), whereas it is very very similar to the first one.
I can avoid the error by using intermediate variables (and it translates into pure C):
vardc = varc128
vard = varf64
varc128 = vardc * vard
or simply by casting (but it does not translate into pure C):
vardc = <double complex>varc128 * <double>varf64
So what happens? What is the meaning of the compilation error? Is there a clean way to avoid it? Why does the multiplication of a np.complex128_t and np.float64_t seem to involve Python calls?
Versions
Cython version 0.22 (most recent version in Pypi when the question was asked) and GCC 4.9.2.
Repository
I created a tiny repository with the example (hg clone https://bitbucket.org/paugier/test_cython_complex) and a tiny Makefile with 3 targets (make clean, make build, make html) so it is easy to test anything.
The simplest way I can find to work around this issue is to simply switch the order of multiplication.
If in testcplx.pyx I change
varc128 = varc128 * varf64
to
varc128 = varf64 * varc128
I change from the failing situation to described to one that works correctly. This scenario is useful as it allows a direct diff of the produced C code.
tl;dr
The order of the multiplication changes the translation, meaning that in the failing version the multiplication is attempted via __pyx_t_npy_float64_complex types, whereas in the working version it is done via __pyx_t_double_complex types. This in turn introduces the typedef line typedef npy_float64 _Complex __pyx_t_npy_float64_complex;, which is invalid.
I am fairly sure this is a cython bug (Update: reported here). Although this is a very old gcc bug report, the response explicitly states (in saying that it is not, in fact, a gcc bug, but user code error):
typedef R _Complex C;This is not valid code; you can't use _Complex together with a typedef, only together with "float", "double" or "long double" in one of the forms listed in C99.
They conclude that double _Complex is a valid type specifier whereas ArbitraryType _Complex is not. This more recent report has the same type of response - trying to use _Complex on a non fundamental type is outside spec, and the GCC manual indicates that _Complex can only be used with float, double and long double
So - we can hack the cython generated C code to test that: replace typedef npy_float64 _Complex __pyx_t_npy_float64_complex; with typedef double _Complex __pyx_t_npy_float64_complex; and verify that it is indeed valid and can make the output code compile.
Short trek through the code
Swapping the multiplication order only highlights the problem that we are told about by the compiler. In the first case, the offending line is the one that says typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - it is trying to assign the type npy_float64 and use the keyword _Complex to the type __pyx_t_npy_float64_complex.
float _Complex or double _Complex is a valid type, whereas npy_float64 _Complex is not. To see the effect, you can just delete npy_float64 from that line, or replace it with double or float and the code compiles fine. The next question is why that line is produced in the first place...
This seems to be produced by this line in the Cython source code.
Why does the order of the multiplication change the code significantly - such that the type __pyx_t_npy_float64_complex is introduced, and introduced in a way that fails?
In the failing instance, the code to implement the multiplication turns varf64 into a __pyx_t_npy_float64_complex type, does the multiplication on real and imaginary parts and then reassembles the complex number. In the working version, it does the product directly via the __pyx_t_double_complex type using the function __Pyx_c_prod
I guess this is as simple as the cython code taking its cue for which type to use for the multiplication from the first variable it encounters. In the first case, it sees a float 64, so produces (invalid) C code based on that, whereas in the second, it sees the (double) complex128 type and bases its translation on that. This explanation is a little hand-wavy and I hope to return to an analysis of it if time allows...
A note on this - here we see that the typedef for npy_float64 is double, so in this particular case, a fix might consist of modifying the code here to use double _Complex where type is npy_float64, but this is getting beyond the scope of a SO answer and doesn't present a general solution.
C code diff result
Working version
Creates this C code from the line `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);
Failing version
Creates this C code from the line 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));
Which necessitates these extra imports - and the offending line is the one that says typedef npy_float64 _Complex __pyx_t_npy_float64_complex; - it is trying to assign the type npy_float64 and the type _Complex to the type __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
这篇关于Cython中的复数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
