python - 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: error: two or more data types in declaration specifiers”，似乎是由于行 typedef npy_float64 _Complex __pyx_t_npy_float64_complex; ) .已经报告了此错误(例如 here )，但我没有找到任何好的解释和/或干净的解决方案。

不包括 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 创建了一个小型存储库，因此可以轻松测试任何内容。

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 错误(更新: reported here)。虽然 this is a very old gcc bug report ，响应明确指出(实际上，它不是 gcc 错误，而是用户代码错误):

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.

他们得出的结论是 double _Complex是有效的类型说明符，而 ArbitraryType _Complex不是。 This more recent report具有相同类型的响应 - 尝试使用 _Complex非基本类型超出规范， GCC manual表示 _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并且代码编译得很好。下一个问题是为什么首先生产这条线......

这似乎是由 this line 产生的在 Cython 源代码中。

为什么乘法的顺序会显着改变代码 - 这样类型 __pyx_t_npy_float64_complex被引入，并以失败的方式引入？

在失败的实例中，实现乘法的代码变成 varf64成 __pyx_t_npy_float64_complex type，对实部和虚部进行乘法，然后重新组合复数。在工作版本中，它直接通过 __pyx_t_double_complex 做产品。使用函数输入 __Pyx_c_prod
我想这就像 cython 代码从它遇到的第一个变量中获取用于乘法的类型的提示一样简单。在第一种情况下，它看到一个浮点数 64，因此基于它生成(无效)C 代码，而在第二种情况下，它看到(双)complex128 类型并以此为基础进行翻译。这个解释有点笨拙，如果时间允许，我希望回到对它的分析......

关于此的说明 - here we see typedef为 npy_float64是 double ，因此在这种特殊情况下，修复可能包括修改 the code here使用 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);

失败版本

从行 varc128 = varc128 * varf64 创建此 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));

这需要这些额外的进口 - 而违规行是说 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