Cython中的复数

Question

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

我想用努比·恩达雷数据类型的np.复合物128. 在Cython中，关联的C类型定义在 Cython/Includes/numpy/__init__.pxd作为

ctypedef double complex complex128_t

看来这只是一个简单的C双复数。在

然而，很容易获得奇怪的行为。特别是，有了这些定义

^{pr2}$

线

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和{}是橙色的（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

那会发生什么？编译错误的含义是什么？有没有一个干净的方法来避免它？为什么a的乘法np.复合物128以及np.float64似乎涉及Python调用？在

版本

Cython版本0.22（当被问到这个问题时，Pypi中的最新版本）和gcc4.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

到

^{pr2}$

我从失败的情况改为正确的描述。这个场景很有用，因为它允许直接区分生成的C代码。在

tl；dr

乘法的顺序改变了转换，这意味着在失败的版本中，乘法是通过__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表示{}只能与float、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;-它试图将类型^{和使用关键字_Complex给类型__pyx_t_npy_float64_complex。在

float _Complex或double _Complex是有效类型，而npy_float64 _Complex不是。要查看效果，只需从该行中删除npy_float64，或者将其替换为double或{}，代码编译良好。下一个问题是为什么这条生产线首先。。。在

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

为什么乘法的顺序会显著地改变代码-以至于引入__pyx_t_npy_float64_complex类型，并且以失败的方式引入？在

在失败的实例中，实现乘法的代码将varf64转换为__pyx_t_npy_float64_complex类型，对实部和虚部进行乘法，然后重新组合复数。在工作版本中，它使用函数__Pyx_c_prod直接通过__pyx_t_double_complex类型生成产品

我想这很简单，就像cython代码从遇到的第一个变量中获取用于乘法的类型的提示一样简单。在第一种情况下，它看到一个float 64，因此基于它生成（invalid）C代码，而在第二种情况下，它看到（double）complex128类型，并基于该类型进行转换。这个解释有点手忙脚乱，如果时间允许，我希望回到对它的分析上。。。在

关于这个-here we see的一个注释是npy_float64的{}是{}，因此在这个特定的情况下，修复可能包括修改the code here以使用double _Complex，其中type是{}，但这超出了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));

这是一个>>>>>>>>>>><

#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