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

我想使用 dtype np.complex128 的 numpy.ndarray 编写一个纯 C 循环。在 Cython 中，关联的 C 类型定义在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.

the line

varc128 = varc128 * varf64

可以由 Cython 编译，但 gcc 无法编译生成的 C 代码（错误为“testcplx.c:663:25：错误：声明说明符中的两个或多个数据类型”，似乎是由于该行typedef npy_float64 _Complex __pyx_t_npy_float64_complex;）。这个错误已经被报告过（例如here http://comments.gmane.org/gmane.comp.python.cython.user/10659）但我没有找到任何好的解释和/或干净的解决方案。

不包括complex.h，没有错误（我猜是因为typedef则不包括在内）。

然而，仍然存在一个问题，因为在由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));

and the __Pyx_CREAL and __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 调用？

Versions

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）所以测试任何东西都很容易。

我能找到解决此问题的最简单方法是简单地切换乘法顺序。

If in testcplx.pyx我改变

varc128 = varc128 * varf64

to

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 错误（更新：在这里报道 http://trac.cython.org/ticket/850#ticket）。虽然这是一个非常古老的 gcc 错误报告 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=19514，响应明确指出（事实上，它不是gccbug，但用户代码错误）：

typedef R _Complex C;

这不是有效的代码；你不能将 _Complex 与 typedef 一起使用， 仅与其中一种形式的“float”、“double”或“long double”一起使用 列于C99。

他们的结论是double _Complex是一个有效的类型说明符，而ArbitraryType _Complex不是。这份最新报告 https://gcc.gnu.org/bugzilla/show_bug.cgi?id=36692具有相同类型的响应 - 尝试使用_Complex非基本类型超出规范，并且海湾合作委员会手册 http://www.gnu.org/software/gnu-c-manual/gnu-c-manual.html#Standard-Complex-Number-Types表明_Complex只能与float, double and long double

所以 - 我们可以破解 cython 生成的 C 代码来测试：替换typedef npy_float64 _Complex __pyx_t_npy_float64_complex; with typedef double _Complex __pyx_t_npy_float64_complex;并验证它确实有效并且可以使输出代码编译。

短途浏览代码

交换乘法顺序只会突出编译器告诉我们的问题。在第一种情况下，有问题的行是这样的行：typedef npy_float64 _Complex __pyx_t_npy_float64_complex;- 它正在尝试分配类型npy_float64 and使用关键字_Complex到类型__pyx_t_npy_float64_complex.

float _Complex or double _Complex是一个有效的类型，而npy_float64 _Complex不是。删除即可看到效果npy_float64从该行开始，或将其替换为double or float并且代码编译良好。下一个问题是为什么首先要生产这条线......

这似乎是由这条线 https://github.com/cython/cython/blob/bb4d9c2de71b7c7e1e02d9dfeae53f4547fa9d7d/Cython/Compiler/PyrexTypes.py#L1949在 Cython 源代码中。

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

在失败的实例中，实现乘法的代码变成varf64 into a __pyx_t_npy_float64_complex类型，对实部和虚部进行乘法，然后重新组合复数。在工作版本中，它直接通过__pyx_t_double_complex使用函数输入__Pyx_c_prod

我想这就像 cython 代码从它遇到的第一个变量中获取用于乘法的类型的提示一样简单。在第一种情况下，它看到浮点数 64，因此产生 (invalid）基于此的 C 代码，而在第二个中，它看到（双）complex128 类型并以此为基础进行翻译。这个解释有点夸张，如果时间允许，我希望能回到对它的分析......

关于此的注释 -在这里我们看到 https://github.com/cython/cython/blob/e3f5343f3b648fc0033bdaf0def3268abab7b9ea/Cython/Includes/numpy/__init__.pxd#L331认为typedef for npy_float64 is double，因此在这种特殊情况下，修复可能包括修改代码在这里 https://github.com/cython/cython/blob/bb4d9c2de71b7c7e1e02d9dfeae53f4547fa9d7d/Cython/Compiler/Nodes.py#L1011 to use double _Complex where type is 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 and方式_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