Complex.c 9.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292
  1. /////////////// Header.proto ///////////////
  2. //@proto_block: h_code
  3. #if !defined(CYTHON_CCOMPLEX)
  4. #if defined(__cplusplus)
  5. #define CYTHON_CCOMPLEX 1
  6. #elif (defined(_Complex_I) && !defined(_MSC_VER))
  7. // MSVC defines "_Complex_I" but not "_Complex". See https://github.com/cython/cython/issues/5512
  8. #define CYTHON_CCOMPLEX 1
  9. #else
  10. #define CYTHON_CCOMPLEX 0
  11. #endif
  12. #endif
  13. #if CYTHON_CCOMPLEX
  14. #ifdef __cplusplus
  15. #include <complex>
  16. #else
  17. #include <complex.h>
  18. #endif
  19. #endif
  20. #if CYTHON_CCOMPLEX && !defined(__cplusplus) && defined(__sun__) && defined(__GNUC__)
  21. #undef _Complex_I
  22. #define _Complex_I 1.0fj
  23. #endif
  24. /////////////// RealImag.proto ///////////////
  25. #if CYTHON_CCOMPLEX
  26. #ifdef __cplusplus
  27. #define __Pyx_CREAL(z) ((z).real())
  28. #define __Pyx_CIMAG(z) ((z).imag())
  29. #else
  30. #define __Pyx_CREAL(z) (__real__(z))
  31. #define __Pyx_CIMAG(z) (__imag__(z))
  32. #endif
  33. #else
  34. #define __Pyx_CREAL(z) ((z).real)
  35. #define __Pyx_CIMAG(z) ((z).imag)
  36. #endif
  37. #if defined(__cplusplus) && CYTHON_CCOMPLEX \
  38. && (defined(_WIN32) || defined(__clang__) || (defined(__GNUC__) && (__GNUC__ >= 5 || __GNUC__ == 4 && __GNUC_MINOR__ >= 4 )) || __cplusplus >= 201103)
  39. #define __Pyx_SET_CREAL(z,x) ((z).real(x))
  40. #define __Pyx_SET_CIMAG(z,y) ((z).imag(y))
  41. #else
  42. #define __Pyx_SET_CREAL(z,x) __Pyx_CREAL(z) = (x)
  43. #define __Pyx_SET_CIMAG(z,y) __Pyx_CIMAG(z) = (y)
  44. #endif
  45. /////////////// Declarations.proto ///////////////
  46. //@proto_block: complex_type_declarations
  47. #if CYTHON_CCOMPLEX
  48. #ifdef __cplusplus
  49. typedef ::std::complex< {{real_type}} > {{type_name}};
  50. #else
  51. typedef {{real_type}} _Complex {{type_name}};
  52. #endif
  53. #else
  54. typedef struct { {{real_type}} real, imag; } {{type_name}};
  55. #endif
  56. static CYTHON_INLINE {{type}} {{type_name}}_from_parts({{real_type}}, {{real_type}});
  57. /////////////// Declarations ///////////////
  58. #if CYTHON_CCOMPLEX
  59. #ifdef __cplusplus
  60. static CYTHON_INLINE {{type}} {{type_name}}_from_parts({{real_type}} x, {{real_type}} y) {
  61. return ::std::complex< {{real_type}} >(x, y);
  62. }
  63. #else
  64. static CYTHON_INLINE {{type}} {{type_name}}_from_parts({{real_type}} x, {{real_type}} y) {
  65. return x + y*({{type}})_Complex_I;
  66. }
  67. #endif
  68. #else
  69. static CYTHON_INLINE {{type}} {{type_name}}_from_parts({{real_type}} x, {{real_type}} y) {
  70. {{type}} z;
  71. z.real = x;
  72. z.imag = y;
  73. return z;
  74. }
  75. #endif
  76. /////////////// ToPy.proto ///////////////
  77. #define __pyx_PyComplex_FromComplex(z) \
  78. PyComplex_FromDoubles((double)__Pyx_CREAL(z), \
  79. (double)__Pyx_CIMAG(z))
  80. /////////////// FromPy.proto ///////////////
  81. static {{type}} __Pyx_PyComplex_As_{{type_name}}(PyObject*);
  82. /////////////// FromPy ///////////////
  83. static {{type}} __Pyx_PyComplex_As_{{type_name}}(PyObject* o) {
  84. Py_complex cval;
  85. #if !CYTHON_COMPILING_IN_PYPY
  86. if (PyComplex_CheckExact(o))
  87. cval = ((PyComplexObject *)o)->cval;
  88. else
  89. #endif
  90. cval = PyComplex_AsCComplex(o);
  91. return {{type_name}}_from_parts(
  92. ({{real_type}})cval.real,
  93. ({{real_type}})cval.imag);
  94. }
  95. /////////////// Arithmetic.proto ///////////////
  96. #if CYTHON_CCOMPLEX
  97. #define __Pyx_c_eq{{func_suffix}}(a, b) ((a)==(b))
  98. #define __Pyx_c_sum{{func_suffix}}(a, b) ((a)+(b))
  99. #define __Pyx_c_diff{{func_suffix}}(a, b) ((a)-(b))
  100. #define __Pyx_c_prod{{func_suffix}}(a, b) ((a)*(b))
  101. #define __Pyx_c_quot{{func_suffix}}(a, b) ((a)/(b))
  102. #define __Pyx_c_neg{{func_suffix}}(a) (-(a))
  103. #ifdef __cplusplus
  104. #define __Pyx_c_is_zero{{func_suffix}}(z) ((z)==({{real_type}})0)
  105. #define __Pyx_c_conj{{func_suffix}}(z) (::std::conj(z))
  106. #if {{is_float}}
  107. #define __Pyx_c_abs{{func_suffix}}(z) (::std::abs(z))
  108. #define __Pyx_c_pow{{func_suffix}}(a, b) (::std::pow(a, b))
  109. #endif
  110. #else
  111. #define __Pyx_c_is_zero{{func_suffix}}(z) ((z)==0)
  112. #define __Pyx_c_conj{{func_suffix}}(z) (conj{{m}}(z))
  113. #if {{is_float}}
  114. #define __Pyx_c_abs{{func_suffix}}(z) (cabs{{m}}(z))
  115. #define __Pyx_c_pow{{func_suffix}}(a, b) (cpow{{m}}(a, b))
  116. #endif
  117. #endif
  118. #else
  119. static CYTHON_INLINE int __Pyx_c_eq{{func_suffix}}({{type}}, {{type}});
  120. static CYTHON_INLINE {{type}} __Pyx_c_sum{{func_suffix}}({{type}}, {{type}});
  121. static CYTHON_INLINE {{type}} __Pyx_c_diff{{func_suffix}}({{type}}, {{type}});
  122. static CYTHON_INLINE {{type}} __Pyx_c_prod{{func_suffix}}({{type}}, {{type}});
  123. static CYTHON_INLINE {{type}} __Pyx_c_quot{{func_suffix}}({{type}}, {{type}});
  124. static CYTHON_INLINE {{type}} __Pyx_c_neg{{func_suffix}}({{type}});
  125. static CYTHON_INLINE int __Pyx_c_is_zero{{func_suffix}}({{type}});
  126. static CYTHON_INLINE {{type}} __Pyx_c_conj{{func_suffix}}({{type}});
  127. #if {{is_float}}
  128. static CYTHON_INLINE {{real_type}} __Pyx_c_abs{{func_suffix}}({{type}});
  129. static CYTHON_INLINE {{type}} __Pyx_c_pow{{func_suffix}}({{type}}, {{type}});
  130. #endif
  131. #endif
  132. /////////////// Arithmetic ///////////////
  133. #if CYTHON_CCOMPLEX
  134. #else
  135. static CYTHON_INLINE int __Pyx_c_eq{{func_suffix}}({{type}} a, {{type}} b) {
  136. return (a.real == b.real) && (a.imag == b.imag);
  137. }
  138. static CYTHON_INLINE {{type}} __Pyx_c_sum{{func_suffix}}({{type}} a, {{type}} b) {
  139. {{type}} z;
  140. z.real = a.real + b.real;
  141. z.imag = a.imag + b.imag;
  142. return z;
  143. }
  144. static CYTHON_INLINE {{type}} __Pyx_c_diff{{func_suffix}}({{type}} a, {{type}} b) {
  145. {{type}} z;
  146. z.real = a.real - b.real;
  147. z.imag = a.imag - b.imag;
  148. return z;
  149. }
  150. static CYTHON_INLINE {{type}} __Pyx_c_prod{{func_suffix}}({{type}} a, {{type}} b) {
  151. {{type}} z;
  152. z.real = a.real * b.real - a.imag * b.imag;
  153. z.imag = a.real * b.imag + a.imag * b.real;
  154. return z;
  155. }
  156. #if {{is_float}}
  157. static CYTHON_INLINE {{type}} __Pyx_c_quot{{func_suffix}}({{type}} a, {{type}} b) {
  158. if (b.imag == 0) {
  159. return {{type_name}}_from_parts(a.real / b.real, a.imag / b.real);
  160. } else if (fabs{{m}}(b.real) >= fabs{{m}}(b.imag)) {
  161. if (b.real == 0 && b.imag == 0) {
  162. return {{type_name}}_from_parts(a.real / b.real, a.imag / b.imag);
  163. } else {
  164. {{real_type}} r = b.imag / b.real;
  165. {{real_type}} s = ({{real_type}})(1.0) / (b.real + b.imag * r);
  166. return {{type_name}}_from_parts(
  167. (a.real + a.imag * r) * s, (a.imag - a.real * r) * s);
  168. }
  169. } else {
  170. {{real_type}} r = b.real / b.imag;
  171. {{real_type}} s = ({{real_type}})(1.0) / (b.imag + b.real * r);
  172. return {{type_name}}_from_parts(
  173. (a.real * r + a.imag) * s, (a.imag * r - a.real) * s);
  174. }
  175. }
  176. #else
  177. static CYTHON_INLINE {{type}} __Pyx_c_quot{{func_suffix}}({{type}} a, {{type}} b) {
  178. if (b.imag == 0) {
  179. return {{type_name}}_from_parts(a.real / b.real, a.imag / b.real);
  180. } else {
  181. {{real_type}} denom = b.real * b.real + b.imag * b.imag;
  182. return {{type_name}}_from_parts(
  183. (a.real * b.real + a.imag * b.imag) / denom,
  184. (a.imag * b.real - a.real * b.imag) / denom);
  185. }
  186. }
  187. #endif
  188. static CYTHON_INLINE {{type}} __Pyx_c_neg{{func_suffix}}({{type}} a) {
  189. {{type}} z;
  190. z.real = -a.real;
  191. z.imag = -a.imag;
  192. return z;
  193. }
  194. static CYTHON_INLINE int __Pyx_c_is_zero{{func_suffix}}({{type}} a) {
  195. return (a.real == 0) && (a.imag == 0);
  196. }
  197. static CYTHON_INLINE {{type}} __Pyx_c_conj{{func_suffix}}({{type}} a) {
  198. {{type}} z;
  199. z.real = a.real;
  200. z.imag = -a.imag;
  201. return z;
  202. }
  203. #if {{is_float}}
  204. static CYTHON_INLINE {{real_type}} __Pyx_c_abs{{func_suffix}}({{type}} z) {
  205. #if !defined(HAVE_HYPOT) || defined(_MSC_VER)
  206. return sqrt{{m}}(z.real*z.real + z.imag*z.imag);
  207. #else
  208. return hypot{{m}}(z.real, z.imag);
  209. #endif
  210. }
  211. static CYTHON_INLINE {{type}} __Pyx_c_pow{{func_suffix}}({{type}} a, {{type}} b) {
  212. {{type}} z;
  213. {{real_type}} r, lnr, theta, z_r, z_theta;
  214. if (b.imag == 0 && b.real == (int)b.real) {
  215. if (b.real < 0) {
  216. {{real_type}} denom = a.real * a.real + a.imag * a.imag;
  217. a.real = a.real / denom;
  218. a.imag = -a.imag / denom;
  219. b.real = -b.real;
  220. }
  221. switch ((int)b.real) {
  222. case 0:
  223. z.real = 1;
  224. z.imag = 0;
  225. return z;
  226. case 1:
  227. return a;
  228. case 2:
  229. return __Pyx_c_prod{{func_suffix}}(a, a);
  230. case 3:
  231. z = __Pyx_c_prod{{func_suffix}}(a, a);
  232. return __Pyx_c_prod{{func_suffix}}(z, a);
  233. case 4:
  234. z = __Pyx_c_prod{{func_suffix}}(a, a);
  235. return __Pyx_c_prod{{func_suffix}}(z, z);
  236. }
  237. }
  238. if (a.imag == 0) {
  239. if (a.real == 0) {
  240. return a;
  241. } else if ((b.imag == 0) && (a.real >= 0)) {
  242. z.real = pow{{m}}(a.real, b.real);
  243. z.imag = 0;
  244. return z;
  245. } else if (a.real > 0) {
  246. r = a.real;
  247. theta = 0;
  248. } else {
  249. r = -a.real;
  250. theta = atan2{{m}}(0.0, -1.0);
  251. }
  252. } else {
  253. r = __Pyx_c_abs{{func_suffix}}(a);
  254. theta = atan2{{m}}(a.imag, a.real);
  255. }
  256. lnr = log{{m}}(r);
  257. z_r = exp{{m}}(lnr * b.real - theta * b.imag);
  258. z_theta = theta * b.real + lnr * b.imag;
  259. z.real = z_r * cos{{m}}(z_theta);
  260. z.imag = z_r * sin{{m}}(z_theta);
  261. return z;
  262. }
  263. #endif
  264. #endif