Cython中的复数

在Cython中使用复数的正确方法是什么？

我想写一个纯粹的C循环使用dtype np.complex128的numpy.ndarray。 在Cython中，关联的Ctypes是在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：error：两个或多个数据types在声明说明符中”，似乎是由于typedef npy_float64 _Complex __pyx_t_npy_float64_complex; ）。 这个错误已经被报告（例如这里 ），但我没有find任何好的解释和/或干净的解决scheme。

没有包含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)); ），而它与第一个非常相似。

我可以通过使用中间variables来避免错误（并将其转换为纯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创build了一个小的版本库，所以很容易testing任何东西。

Answer 1

我可以find解决这个问题最简单的方法是简单地切换乘法的顺序。

如果在testcplx.pyx我改变了

 varc128 = varc128 * varf64 

至

 varc128 = varf64 * varc128 

我从失败的情况转变为描述正确的工作。 这个场景很有用，因为它允许直接比较生成的C代码。

TL;博士

乘法的顺序改变了翻译，这意味着在失败的版本中，通过__pyx_t_npy_float64_complextypes来尝试乘法，而在工作版本中则通过__pyx_t_double_complextypes来完成。 这又引入了typedef行typedef npy_float64 _Complex __pyx_t_npy_float64_complex; ，这是无效的。

我相当确定这是一个cython bug（更新： 在这里报道 ）。 虽然这是一个非常古老的海湾合作委员会的错误报告 ，响应明确指出（事实上，它不是，一个海湾合作委员会的错误，但用户代码错误）：

 typedef R _Complex C; 

这不是有效的代码; 您不能使用_Complex与typedef一起使用，只能使用C99中列出的某种forms的“float”，“double”或“long double”。

他们得出结论double _Complex是一个有效的types说明符，而ArbitraryType _Complex不是。 这个更近期的报告具有相同types的响应 – 尝试在非基本types上使用_Complex是超出规范的，而GCC手册则指出_Complex只能用于float ， double和long double

所以 – 我们可以破解cython生成的C代码来testing：replacetypedef npy_float64 _Complex __pyx_t_npy_float64_complex; 用typedef double _Complex __pyx_t_npy_float64_complex; 并validation它确实是有效的，并且可以使输出代码编译。

---

通过代码短暂跋涉

乘法顺序的交换只突出了编译器告诉我们的问题。 在第一种情况下，违规行是说typedef npy_float64 _Complex __pyx_t_npy_float64_complex; – 它试图分配typesnpy_float64 并使用关键字_Complex到types__pyx_t_npy_float64_complex 。

float _Complex或double _Complex是一个有效的types，而npy_float64 _Complex不是。 为了看到这个效果，你可以直接从该行删除npy_float64 ，或者用double或者float代替它，代码编译的很好。 下一个问题是为什么这条线是在第一个地方生产…

这似乎是由Cython源代码中的这一行产生的。

为什么乘法的顺序会显着改变代码 – 这样引入了types__pyx_t_npy_float64_complex ，并且以失败的方式引入？

在失败的实例中，实现乘法的代码将varf64转换为__pyx_t_npy_float64_complextypes，对实部和虚部进行乘法运算，然后重新组合复数。 在工作版本中，它直接通过__pyx_t_double_complextypes使用函数__Pyx_c_prod

我想这和cython代码一样简单，就是从它遇到的第一个variables中select哪种types的乘法。 在第一种情况下，它看到一个浮点数64，因此会产生（ 无效的 ）C代码，而在第二种情况下，它会看到（double）complex128types，并以此为基础进行翻译。 这个解释是有点手摇的，如果时间允许的话，我希望能回到分析。

关于这个的一个注释 – 在这里我们看到 npy_float64的typedef是double ，所以在这种情况下，一个修复可能包括修改这里的代码来使用double _Complex ，其中type是npy_float64 ，但是这超出了SO的范围答案并不提供一个通用的解决scheme。

---

C代码差异的结果

工作版本

从varc128 = varf64 * varc128这一行创build这个C代码

 __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这一行创build这个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)); 

这需要这些额外的import – 有罪的行是说typedef npy_float64 _Complex __pyx_t_npy_float64_complex; – 它试图将typesnpy_float64 和types_Complex分配给types__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