Cython中的复数

Question

在Cython中处理复杂数字的正确方法是什么？

我想用dtype np.complex128的numpy.ndarray编写一个纯C循环.在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版本0.22(询问问题时Pypi的最新版本)和GCC 4.9.2.

知识库

我创建的例子(一个小库hg clone https://bitbucket.org/paugier/test_cython_complex),并用3个目标(一个微小的Makefile make clean,make build,make html),所以很容易测试任何东西.

Answer 1

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

如果在testcplx.pyx我改变

varc128 = varc128 * varf64

至

varc128 = varf64 * varc128

我从失败的情况改为描述为正常的情况.这种情况很有用,因为它允许生成的C代码的直接差异.

TL;博士

乘法的顺序改变了转换,这意味着在失败的版本中,乘法是通过__pyx_t_npy_float64_complex类型尝试的,而在工作版本中,它是通过__pyx_t_double_complex类型完成的.这又引入了typedef行typedef npy_float64 _Complex __pyx_t_npy_float64_complex;,该行无效.

我很确定这是一个cython bug(更新:在这里报道).虽然这是一个非常古老的gcc错误报告,但响应明确指出(事实上,它不是gcc错误,而是用户代码错误):

typedef R _Complex C;

这不是有效的代码; 您不能将_Complex与typedef一起使用,只能与C99中列出的其中一种形式的"float","double"或"long double"一起使用.

他们得出结论,这double _Complex是一个有效的类型说明符,而ArbitraryType _Complex不是. 这个更新的报告具有相同类型的响应 - 尝试_Complex在非基本类型上使用是在规范之外,并且GCC手册指示_Complex只能使用float,double和long double

所以-我们可以破解用Cython生成C代码来测试:更换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代码一样简单,它从哪个类型用于从它遇到的第一个变量中进行乘法.在第一种情况下,它会看到一个浮点数64,因此根据它生成(无效的)C代码,而在第二种情况下,它会看到(double)complex128类型并将其转换基于该类型.这个解释有点波动,如果时间允许,我希望回到它的分析......

基于此的一个说明- 在这里我们看到的是typedef对npy_float64的double,所以在这种特殊情况下,修复可能包括修改的这里的代码使用的double _Complex地方type是npy_float64,但是这能超过的SO回答的范围,并没有呈现一般解决方案

C代码差异结果

工作版

从`varc128 = varf64*varc128行创建此C代码

__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