Complex.c 9.8 KB

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