1 |
|
---|
2 | /* Complex object implementation */
|
---|
3 |
|
---|
4 | /* Borrows heavily from floatobject.c */
|
---|
5 |
|
---|
6 | /* Submitted by Jim Hugunin */
|
---|
7 |
|
---|
8 | #include "Python.h"
|
---|
9 | #include "structmember.h"
|
---|
10 |
|
---|
11 | #ifndef WITHOUT_COMPLEX
|
---|
12 |
|
---|
13 | /* Precisions used by repr() and str(), respectively.
|
---|
14 |
|
---|
15 | The repr() precision (17 significant decimal digits) is the minimal number
|
---|
16 | that is guaranteed to have enough precision so that if the number is read
|
---|
17 | back in the exact same binary value is recreated. This is true for IEEE
|
---|
18 | floating point by design, and also happens to work for all other modern
|
---|
19 | hardware.
|
---|
20 |
|
---|
21 | The str() precision is chosen so that in most cases, the rounding noise
|
---|
22 | created by various operations is suppressed, while giving plenty of
|
---|
23 | precision for practical use.
|
---|
24 | */
|
---|
25 |
|
---|
26 | #define PREC_REPR 17
|
---|
27 | #define PREC_STR 12
|
---|
28 |
|
---|
29 | /* elementary operations on complex numbers */
|
---|
30 |
|
---|
31 | static Py_complex c_1 = {1., 0.};
|
---|
32 |
|
---|
33 | Py_complex
|
---|
34 | c_sum(Py_complex a, Py_complex b)
|
---|
35 | {
|
---|
36 | Py_complex r;
|
---|
37 | r.real = a.real + b.real;
|
---|
38 | r.imag = a.imag + b.imag;
|
---|
39 | return r;
|
---|
40 | }
|
---|
41 |
|
---|
42 | Py_complex
|
---|
43 | c_diff(Py_complex a, Py_complex b)
|
---|
44 | {
|
---|
45 | Py_complex r;
|
---|
46 | r.real = a.real - b.real;
|
---|
47 | r.imag = a.imag - b.imag;
|
---|
48 | return r;
|
---|
49 | }
|
---|
50 |
|
---|
51 | Py_complex
|
---|
52 | c_neg(Py_complex a)
|
---|
53 | {
|
---|
54 | Py_complex r;
|
---|
55 | r.real = -a.real;
|
---|
56 | r.imag = -a.imag;
|
---|
57 | return r;
|
---|
58 | }
|
---|
59 |
|
---|
60 | Py_complex
|
---|
61 | c_prod(Py_complex a, Py_complex b)
|
---|
62 | {
|
---|
63 | Py_complex r;
|
---|
64 | r.real = a.real*b.real - a.imag*b.imag;
|
---|
65 | r.imag = a.real*b.imag + a.imag*b.real;
|
---|
66 | return r;
|
---|
67 | }
|
---|
68 |
|
---|
69 | Py_complex
|
---|
70 | c_quot(Py_complex a, Py_complex b)
|
---|
71 | {
|
---|
72 | /******************************************************************
|
---|
73 | This was the original algorithm. It's grossly prone to spurious
|
---|
74 | overflow and underflow errors. It also merrily divides by 0 despite
|
---|
75 | checking for that(!). The code still serves a doc purpose here, as
|
---|
76 | the algorithm following is a simple by-cases transformation of this
|
---|
77 | one:
|
---|
78 |
|
---|
79 | Py_complex r;
|
---|
80 | double d = b.real*b.real + b.imag*b.imag;
|
---|
81 | if (d == 0.)
|
---|
82 | errno = EDOM;
|
---|
83 | r.real = (a.real*b.real + a.imag*b.imag)/d;
|
---|
84 | r.imag = (a.imag*b.real - a.real*b.imag)/d;
|
---|
85 | return r;
|
---|
86 | ******************************************************************/
|
---|
87 |
|
---|
88 | /* This algorithm is better, and is pretty obvious: first divide the
|
---|
89 | * numerators and denominator by whichever of {b.real, b.imag} has
|
---|
90 | * larger magnitude. The earliest reference I found was to CACM
|
---|
91 | * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
|
---|
92 | * University). As usual, though, we're still ignoring all IEEE
|
---|
93 | * endcases.
|
---|
94 | */
|
---|
95 | Py_complex r; /* the result */
|
---|
96 | const double abs_breal = b.real < 0 ? -b.real : b.real;
|
---|
97 | const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
|
---|
98 |
|
---|
99 | if (abs_breal >= abs_bimag) {
|
---|
100 | /* divide tops and bottom by b.real */
|
---|
101 | if (abs_breal == 0.0) {
|
---|
102 | errno = EDOM;
|
---|
103 | r.real = r.imag = 0.0;
|
---|
104 | }
|
---|
105 | else {
|
---|
106 | const double ratio = b.imag / b.real;
|
---|
107 | const double denom = b.real + b.imag * ratio;
|
---|
108 | r.real = (a.real + a.imag * ratio) / denom;
|
---|
109 | r.imag = (a.imag - a.real * ratio) / denom;
|
---|
110 | }
|
---|
111 | }
|
---|
112 | else {
|
---|
113 | /* divide tops and bottom by b.imag */
|
---|
114 | const double ratio = b.real / b.imag;
|
---|
115 | const double denom = b.real * ratio + b.imag;
|
---|
116 | assert(b.imag != 0.0);
|
---|
117 | r.real = (a.real * ratio + a.imag) / denom;
|
---|
118 | r.imag = (a.imag * ratio - a.real) / denom;
|
---|
119 | }
|
---|
120 | return r;
|
---|
121 | }
|
---|
122 |
|
---|
123 | Py_complex
|
---|
124 | c_pow(Py_complex a, Py_complex b)
|
---|
125 | {
|
---|
126 | Py_complex r;
|
---|
127 | double vabs,len,at,phase;
|
---|
128 | if (b.real == 0. && b.imag == 0.) {
|
---|
129 | r.real = 1.;
|
---|
130 | r.imag = 0.;
|
---|
131 | }
|
---|
132 | else if (a.real == 0. && a.imag == 0.) {
|
---|
133 | if (b.imag != 0. || b.real < 0.)
|
---|
134 | errno = EDOM;
|
---|
135 | r.real = 0.;
|
---|
136 | r.imag = 0.;
|
---|
137 | }
|
---|
138 | else {
|
---|
139 | vabs = hypot(a.real,a.imag);
|
---|
140 | len = pow(vabs,b.real);
|
---|
141 | at = atan2(a.imag, a.real);
|
---|
142 | phase = at*b.real;
|
---|
143 | if (b.imag != 0.0) {
|
---|
144 | len /= exp(at*b.imag);
|
---|
145 | phase += b.imag*log(vabs);
|
---|
146 | }
|
---|
147 | r.real = len*cos(phase);
|
---|
148 | r.imag = len*sin(phase);
|
---|
149 | }
|
---|
150 | return r;
|
---|
151 | }
|
---|
152 |
|
---|
153 | static Py_complex
|
---|
154 | c_powu(Py_complex x, long n)
|
---|
155 | {
|
---|
156 | Py_complex r, p;
|
---|
157 | long mask = 1;
|
---|
158 | r = c_1;
|
---|
159 | p = x;
|
---|
160 | while (mask > 0 && n >= mask) {
|
---|
161 | if (n & mask)
|
---|
162 | r = c_prod(r,p);
|
---|
163 | mask <<= 1;
|
---|
164 | p = c_prod(p,p);
|
---|
165 | }
|
---|
166 | return r;
|
---|
167 | }
|
---|
168 |
|
---|
169 | static Py_complex
|
---|
170 | c_powi(Py_complex x, long n)
|
---|
171 | {
|
---|
172 | Py_complex cn;
|
---|
173 |
|
---|
174 | if (n > 100 || n < -100) {
|
---|
175 | cn.real = (double) n;
|
---|
176 | cn.imag = 0.;
|
---|
177 | return c_pow(x,cn);
|
---|
178 | }
|
---|
179 | else if (n > 0)
|
---|
180 | return c_powu(x,n);
|
---|
181 | else
|
---|
182 | return c_quot(c_1,c_powu(x,-n));
|
---|
183 |
|
---|
184 | }
|
---|
185 |
|
---|
186 | static PyObject *
|
---|
187 | complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
|
---|
188 | {
|
---|
189 | PyObject *op;
|
---|
190 |
|
---|
191 | op = type->tp_alloc(type, 0);
|
---|
192 | if (op != NULL)
|
---|
193 | ((PyComplexObject *)op)->cval = cval;
|
---|
194 | return op;
|
---|
195 | }
|
---|
196 |
|
---|
197 | PyObject *
|
---|
198 | PyComplex_FromCComplex(Py_complex cval)
|
---|
199 | {
|
---|
200 | register PyComplexObject *op;
|
---|
201 |
|
---|
202 | /* Inline PyObject_New */
|
---|
203 | op = (PyComplexObject *) PyObject_MALLOC(sizeof(PyComplexObject));
|
---|
204 | if (op == NULL)
|
---|
205 | return PyErr_NoMemory();
|
---|
206 | PyObject_INIT(op, &PyComplex_Type);
|
---|
207 | op->cval = cval;
|
---|
208 | return (PyObject *) op;
|
---|
209 | }
|
---|
210 |
|
---|
211 | static PyObject *
|
---|
212 | complex_subtype_from_doubles(PyTypeObject *type, double real, double imag)
|
---|
213 | {
|
---|
214 | Py_complex c;
|
---|
215 | c.real = real;
|
---|
216 | c.imag = imag;
|
---|
217 | return complex_subtype_from_c_complex(type, c);
|
---|
218 | }
|
---|
219 |
|
---|
220 | PyObject *
|
---|
221 | PyComplex_FromDoubles(double real, double imag)
|
---|
222 | {
|
---|
223 | Py_complex c;
|
---|
224 | c.real = real;
|
---|
225 | c.imag = imag;
|
---|
226 | return PyComplex_FromCComplex(c);
|
---|
227 | }
|
---|
228 |
|
---|
229 | double
|
---|
230 | PyComplex_RealAsDouble(PyObject *op)
|
---|
231 | {
|
---|
232 | if (PyComplex_Check(op)) {
|
---|
233 | return ((PyComplexObject *)op)->cval.real;
|
---|
234 | }
|
---|
235 | else {
|
---|
236 | return PyFloat_AsDouble(op);
|
---|
237 | }
|
---|
238 | }
|
---|
239 |
|
---|
240 | double
|
---|
241 | PyComplex_ImagAsDouble(PyObject *op)
|
---|
242 | {
|
---|
243 | if (PyComplex_Check(op)) {
|
---|
244 | return ((PyComplexObject *)op)->cval.imag;
|
---|
245 | }
|
---|
246 | else {
|
---|
247 | return 0.0;
|
---|
248 | }
|
---|
249 | }
|
---|
250 |
|
---|
251 | Py_complex
|
---|
252 | PyComplex_AsCComplex(PyObject *op)
|
---|
253 | {
|
---|
254 | Py_complex cv;
|
---|
255 | if (PyComplex_Check(op)) {
|
---|
256 | return ((PyComplexObject *)op)->cval;
|
---|
257 | }
|
---|
258 | else {
|
---|
259 | cv.real = PyFloat_AsDouble(op);
|
---|
260 | cv.imag = 0.;
|
---|
261 | return cv;
|
---|
262 | }
|
---|
263 | }
|
---|
264 |
|
---|
265 | static void
|
---|
266 | complex_dealloc(PyObject *op)
|
---|
267 | {
|
---|
268 | op->ob_type->tp_free(op);
|
---|
269 | }
|
---|
270 |
|
---|
271 |
|
---|
272 | static void
|
---|
273 | complex_to_buf(char *buf, int bufsz, PyComplexObject *v, int precision)
|
---|
274 | {
|
---|
275 | char format[32];
|
---|
276 | if (v->cval.real == 0.) {
|
---|
277 | PyOS_snprintf(format, sizeof(format), "%%.%ig", precision);
|
---|
278 | PyOS_ascii_formatd(buf, bufsz - 1, format, v->cval.imag);
|
---|
279 | strncat(buf, "j", 1);
|
---|
280 | } else {
|
---|
281 | char re[64], im[64];
|
---|
282 | /* Format imaginary part with sign, real part without */
|
---|
283 | PyOS_snprintf(format, sizeof(format), "%%.%ig", precision);
|
---|
284 | PyOS_ascii_formatd(re, sizeof(re), format, v->cval.real);
|
---|
285 | PyOS_snprintf(format, sizeof(format), "%%+.%ig", precision);
|
---|
286 | PyOS_ascii_formatd(im, sizeof(im), format, v->cval.imag);
|
---|
287 | PyOS_snprintf(buf, bufsz, "(%s%sj)", re, im);
|
---|
288 | }
|
---|
289 | }
|
---|
290 |
|
---|
291 | static int
|
---|
292 | complex_print(PyComplexObject *v, FILE *fp, int flags)
|
---|
293 | {
|
---|
294 | char buf[100];
|
---|
295 | complex_to_buf(buf, sizeof(buf), v,
|
---|
296 | (flags & Py_PRINT_RAW) ? PREC_STR : PREC_REPR);
|
---|
297 | fputs(buf, fp);
|
---|
298 | return 0;
|
---|
299 | }
|
---|
300 |
|
---|
301 | static PyObject *
|
---|
302 | complex_repr(PyComplexObject *v)
|
---|
303 | {
|
---|
304 | char buf[100];
|
---|
305 | complex_to_buf(buf, sizeof(buf), v, PREC_REPR);
|
---|
306 | return PyString_FromString(buf);
|
---|
307 | }
|
---|
308 |
|
---|
309 | static PyObject *
|
---|
310 | complex_str(PyComplexObject *v)
|
---|
311 | {
|
---|
312 | char buf[100];
|
---|
313 | complex_to_buf(buf, sizeof(buf), v, PREC_STR);
|
---|
314 | return PyString_FromString(buf);
|
---|
315 | }
|
---|
316 |
|
---|
317 | static long
|
---|
318 | complex_hash(PyComplexObject *v)
|
---|
319 | {
|
---|
320 | long hashreal, hashimag, combined;
|
---|
321 | hashreal = _Py_HashDouble(v->cval.real);
|
---|
322 | if (hashreal == -1)
|
---|
323 | return -1;
|
---|
324 | hashimag = _Py_HashDouble(v->cval.imag);
|
---|
325 | if (hashimag == -1)
|
---|
326 | return -1;
|
---|
327 | /* Note: if the imaginary part is 0, hashimag is 0 now,
|
---|
328 | * so the following returns hashreal unchanged. This is
|
---|
329 | * important because numbers of different types that
|
---|
330 | * compare equal must have the same hash value, so that
|
---|
331 | * hash(x + 0*j) must equal hash(x).
|
---|
332 | */
|
---|
333 | combined = hashreal + 1000003 * hashimag;
|
---|
334 | if (combined == -1)
|
---|
335 | combined = -2;
|
---|
336 | return combined;
|
---|
337 | }
|
---|
338 |
|
---|
339 | static PyObject *
|
---|
340 | complex_add(PyComplexObject *v, PyComplexObject *w)
|
---|
341 | {
|
---|
342 | Py_complex result;
|
---|
343 | PyFPE_START_PROTECT("complex_add", return 0)
|
---|
344 | result = c_sum(v->cval,w->cval);
|
---|
345 | PyFPE_END_PROTECT(result)
|
---|
346 | return PyComplex_FromCComplex(result);
|
---|
347 | }
|
---|
348 |
|
---|
349 | static PyObject *
|
---|
350 | complex_sub(PyComplexObject *v, PyComplexObject *w)
|
---|
351 | {
|
---|
352 | Py_complex result;
|
---|
353 | PyFPE_START_PROTECT("complex_sub", return 0)
|
---|
354 | result = c_diff(v->cval,w->cval);
|
---|
355 | PyFPE_END_PROTECT(result)
|
---|
356 | return PyComplex_FromCComplex(result);
|
---|
357 | }
|
---|
358 |
|
---|
359 | static PyObject *
|
---|
360 | complex_mul(PyComplexObject *v, PyComplexObject *w)
|
---|
361 | {
|
---|
362 | Py_complex result;
|
---|
363 | PyFPE_START_PROTECT("complex_mul", return 0)
|
---|
364 | result = c_prod(v->cval,w->cval);
|
---|
365 | PyFPE_END_PROTECT(result)
|
---|
366 | return PyComplex_FromCComplex(result);
|
---|
367 | }
|
---|
368 |
|
---|
369 | static PyObject *
|
---|
370 | complex_div(PyComplexObject *v, PyComplexObject *w)
|
---|
371 | {
|
---|
372 | Py_complex quot;
|
---|
373 | PyFPE_START_PROTECT("complex_div", return 0)
|
---|
374 | errno = 0;
|
---|
375 | quot = c_quot(v->cval,w->cval);
|
---|
376 | PyFPE_END_PROTECT(quot)
|
---|
377 | if (errno == EDOM) {
|
---|
378 | PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
|
---|
379 | return NULL;
|
---|
380 | }
|
---|
381 | return PyComplex_FromCComplex(quot);
|
---|
382 | }
|
---|
383 |
|
---|
384 | static PyObject *
|
---|
385 | complex_classic_div(PyComplexObject *v, PyComplexObject *w)
|
---|
386 | {
|
---|
387 | Py_complex quot;
|
---|
388 |
|
---|
389 | if (Py_DivisionWarningFlag >= 2 &&
|
---|
390 | PyErr_Warn(PyExc_DeprecationWarning,
|
---|
391 | "classic complex division") < 0)
|
---|
392 | return NULL;
|
---|
393 |
|
---|
394 | PyFPE_START_PROTECT("complex_classic_div", return 0)
|
---|
395 | errno = 0;
|
---|
396 | quot = c_quot(v->cval,w->cval);
|
---|
397 | PyFPE_END_PROTECT(quot)
|
---|
398 | if (errno == EDOM) {
|
---|
399 | PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
|
---|
400 | return NULL;
|
---|
401 | }
|
---|
402 | return PyComplex_FromCComplex(quot);
|
---|
403 | }
|
---|
404 |
|
---|
405 | static PyObject *
|
---|
406 | complex_remainder(PyComplexObject *v, PyComplexObject *w)
|
---|
407 | {
|
---|
408 | Py_complex div, mod;
|
---|
409 |
|
---|
410 | if (PyErr_Warn(PyExc_DeprecationWarning,
|
---|
411 | "complex divmod(), // and % are deprecated") < 0)
|
---|
412 | return NULL;
|
---|
413 |
|
---|
414 | errno = 0;
|
---|
415 | div = c_quot(v->cval,w->cval); /* The raw divisor value. */
|
---|
416 | if (errno == EDOM) {
|
---|
417 | PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
|
---|
418 | return NULL;
|
---|
419 | }
|
---|
420 | div.real = floor(div.real); /* Use the floor of the real part. */
|
---|
421 | div.imag = 0.0;
|
---|
422 | mod = c_diff(v->cval, c_prod(w->cval, div));
|
---|
423 |
|
---|
424 | return PyComplex_FromCComplex(mod);
|
---|
425 | }
|
---|
426 |
|
---|
427 |
|
---|
428 | static PyObject *
|
---|
429 | complex_divmod(PyComplexObject *v, PyComplexObject *w)
|
---|
430 | {
|
---|
431 | Py_complex div, mod;
|
---|
432 | PyObject *d, *m, *z;
|
---|
433 |
|
---|
434 | if (PyErr_Warn(PyExc_DeprecationWarning,
|
---|
435 | "complex divmod(), // and % are deprecated") < 0)
|
---|
436 | return NULL;
|
---|
437 |
|
---|
438 | errno = 0;
|
---|
439 | div = c_quot(v->cval,w->cval); /* The raw divisor value. */
|
---|
440 | if (errno == EDOM) {
|
---|
441 | PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
|
---|
442 | return NULL;
|
---|
443 | }
|
---|
444 | div.real = floor(div.real); /* Use the floor of the real part. */
|
---|
445 | div.imag = 0.0;
|
---|
446 | mod = c_diff(v->cval, c_prod(w->cval, div));
|
---|
447 | d = PyComplex_FromCComplex(div);
|
---|
448 | m = PyComplex_FromCComplex(mod);
|
---|
449 | z = PyTuple_Pack(2, d, m);
|
---|
450 | Py_XDECREF(d);
|
---|
451 | Py_XDECREF(m);
|
---|
452 | return z;
|
---|
453 | }
|
---|
454 |
|
---|
455 | static PyObject *
|
---|
456 | complex_pow(PyComplexObject *v, PyObject *w, PyComplexObject *z)
|
---|
457 | {
|
---|
458 | Py_complex p;
|
---|
459 | Py_complex exponent;
|
---|
460 | long int_exponent;
|
---|
461 |
|
---|
462 | if ((PyObject *)z!=Py_None) {
|
---|
463 | PyErr_SetString(PyExc_ValueError, "complex modulo");
|
---|
464 | return NULL;
|
---|
465 | }
|
---|
466 | PyFPE_START_PROTECT("complex_pow", return 0)
|
---|
467 | errno = 0;
|
---|
468 | exponent = ((PyComplexObject*)w)->cval;
|
---|
469 | int_exponent = (long)exponent.real;
|
---|
470 | if (exponent.imag == 0. && exponent.real == int_exponent)
|
---|
471 | p = c_powi(v->cval,int_exponent);
|
---|
472 | else
|
---|
473 | p = c_pow(v->cval,exponent);
|
---|
474 |
|
---|
475 | PyFPE_END_PROTECT(p)
|
---|
476 | Py_ADJUST_ERANGE2(p.real, p.imag);
|
---|
477 | if (errno == EDOM) {
|
---|
478 | PyErr_SetString(PyExc_ZeroDivisionError,
|
---|
479 | "0.0 to a negative or complex power");
|
---|
480 | return NULL;
|
---|
481 | }
|
---|
482 | else if (errno == ERANGE) {
|
---|
483 | PyErr_SetString(PyExc_OverflowError,
|
---|
484 | "complex exponentiaion");
|
---|
485 | return NULL;
|
---|
486 | }
|
---|
487 | return PyComplex_FromCComplex(p);
|
---|
488 | }
|
---|
489 |
|
---|
490 | static PyObject *
|
---|
491 | complex_int_div(PyComplexObject *v, PyComplexObject *w)
|
---|
492 | {
|
---|
493 | PyObject *t, *r;
|
---|
494 |
|
---|
495 | t = complex_divmod(v, w);
|
---|
496 | if (t != NULL) {
|
---|
497 | r = PyTuple_GET_ITEM(t, 0);
|
---|
498 | Py_INCREF(r);
|
---|
499 | Py_DECREF(t);
|
---|
500 | return r;
|
---|
501 | }
|
---|
502 | return NULL;
|
---|
503 | }
|
---|
504 |
|
---|
505 | static PyObject *
|
---|
506 | complex_neg(PyComplexObject *v)
|
---|
507 | {
|
---|
508 | Py_complex neg;
|
---|
509 | neg.real = -v->cval.real;
|
---|
510 | neg.imag = -v->cval.imag;
|
---|
511 | return PyComplex_FromCComplex(neg);
|
---|
512 | }
|
---|
513 |
|
---|
514 | static PyObject *
|
---|
515 | complex_pos(PyComplexObject *v)
|
---|
516 | {
|
---|
517 | if (PyComplex_CheckExact(v)) {
|
---|
518 | Py_INCREF(v);
|
---|
519 | return (PyObject *)v;
|
---|
520 | }
|
---|
521 | else
|
---|
522 | return PyComplex_FromCComplex(v->cval);
|
---|
523 | }
|
---|
524 |
|
---|
525 | static PyObject *
|
---|
526 | complex_abs(PyComplexObject *v)
|
---|
527 | {
|
---|
528 | double result;
|
---|
529 | PyFPE_START_PROTECT("complex_abs", return 0)
|
---|
530 | result = hypot(v->cval.real,v->cval.imag);
|
---|
531 | PyFPE_END_PROTECT(result)
|
---|
532 | return PyFloat_FromDouble(result);
|
---|
533 | }
|
---|
534 |
|
---|
535 | static int
|
---|
536 | complex_nonzero(PyComplexObject *v)
|
---|
537 | {
|
---|
538 | return v->cval.real != 0.0 || v->cval.imag != 0.0;
|
---|
539 | }
|
---|
540 |
|
---|
541 | static int
|
---|
542 | complex_coerce(PyObject **pv, PyObject **pw)
|
---|
543 | {
|
---|
544 | Py_complex cval;
|
---|
545 | cval.imag = 0.;
|
---|
546 | if (PyInt_Check(*pw)) {
|
---|
547 | cval.real = (double)PyInt_AsLong(*pw);
|
---|
548 | *pw = PyComplex_FromCComplex(cval);
|
---|
549 | Py_INCREF(*pv);
|
---|
550 | return 0;
|
---|
551 | }
|
---|
552 | else if (PyLong_Check(*pw)) {
|
---|
553 | cval.real = PyLong_AsDouble(*pw);
|
---|
554 | if (cval.real == -1.0 && PyErr_Occurred())
|
---|
555 | return -1;
|
---|
556 | *pw = PyComplex_FromCComplex(cval);
|
---|
557 | Py_INCREF(*pv);
|
---|
558 | return 0;
|
---|
559 | }
|
---|
560 | else if (PyFloat_Check(*pw)) {
|
---|
561 | cval.real = PyFloat_AsDouble(*pw);
|
---|
562 | *pw = PyComplex_FromCComplex(cval);
|
---|
563 | Py_INCREF(*pv);
|
---|
564 | return 0;
|
---|
565 | }
|
---|
566 | else if (PyComplex_Check(*pw)) {
|
---|
567 | Py_INCREF(*pv);
|
---|
568 | Py_INCREF(*pw);
|
---|
569 | return 0;
|
---|
570 | }
|
---|
571 | return 1; /* Can't do it */
|
---|
572 | }
|
---|
573 |
|
---|
574 | static PyObject *
|
---|
575 | complex_richcompare(PyObject *v, PyObject *w, int op)
|
---|
576 | {
|
---|
577 | int c;
|
---|
578 | Py_complex i, j;
|
---|
579 | PyObject *res;
|
---|
580 |
|
---|
581 | c = PyNumber_CoerceEx(&v, &w);
|
---|
582 | if (c < 0)
|
---|
583 | return NULL;
|
---|
584 | if (c > 0) {
|
---|
585 | Py_INCREF(Py_NotImplemented);
|
---|
586 | return Py_NotImplemented;
|
---|
587 | }
|
---|
588 | /* Make sure both arguments are complex. */
|
---|
589 | if (!(PyComplex_Check(v) && PyComplex_Check(w))) {
|
---|
590 | Py_DECREF(v);
|
---|
591 | Py_DECREF(w);
|
---|
592 | Py_INCREF(Py_NotImplemented);
|
---|
593 | return Py_NotImplemented;
|
---|
594 | }
|
---|
595 |
|
---|
596 | i = ((PyComplexObject *)v)->cval;
|
---|
597 | j = ((PyComplexObject *)w)->cval;
|
---|
598 | Py_DECREF(v);
|
---|
599 | Py_DECREF(w);
|
---|
600 |
|
---|
601 | if (op != Py_EQ && op != Py_NE) {
|
---|
602 | PyErr_SetString(PyExc_TypeError,
|
---|
603 | "no ordering relation is defined for complex numbers");
|
---|
604 | return NULL;
|
---|
605 | }
|
---|
606 |
|
---|
607 | if ((i.real == j.real && i.imag == j.imag) == (op == Py_EQ))
|
---|
608 | res = Py_True;
|
---|
609 | else
|
---|
610 | res = Py_False;
|
---|
611 |
|
---|
612 | Py_INCREF(res);
|
---|
613 | return res;
|
---|
614 | }
|
---|
615 |
|
---|
616 | static PyObject *
|
---|
617 | complex_int(PyObject *v)
|
---|
618 | {
|
---|
619 | PyErr_SetString(PyExc_TypeError,
|
---|
620 | "can't convert complex to int; use int(abs(z))");
|
---|
621 | return NULL;
|
---|
622 | }
|
---|
623 |
|
---|
624 | static PyObject *
|
---|
625 | complex_long(PyObject *v)
|
---|
626 | {
|
---|
627 | PyErr_SetString(PyExc_TypeError,
|
---|
628 | "can't convert complex to long; use long(abs(z))");
|
---|
629 | return NULL;
|
---|
630 | }
|
---|
631 |
|
---|
632 | static PyObject *
|
---|
633 | complex_float(PyObject *v)
|
---|
634 | {
|
---|
635 | PyErr_SetString(PyExc_TypeError,
|
---|
636 | "can't convert complex to float; use abs(z)");
|
---|
637 | return NULL;
|
---|
638 | }
|
---|
639 |
|
---|
640 | static PyObject *
|
---|
641 | complex_conjugate(PyObject *self)
|
---|
642 | {
|
---|
643 | Py_complex c;
|
---|
644 | c = ((PyComplexObject *)self)->cval;
|
---|
645 | c.imag = -c.imag;
|
---|
646 | return PyComplex_FromCComplex(c);
|
---|
647 | }
|
---|
648 |
|
---|
649 | static PyObject *
|
---|
650 | complex_getnewargs(PyComplexObject *v)
|
---|
651 | {
|
---|
652 | return Py_BuildValue("(D)", &v->cval);
|
---|
653 | }
|
---|
654 |
|
---|
655 | static PyMethodDef complex_methods[] = {
|
---|
656 | {"conjugate", (PyCFunction)complex_conjugate, METH_NOARGS},
|
---|
657 | {"__getnewargs__", (PyCFunction)complex_getnewargs, METH_NOARGS},
|
---|
658 | {NULL, NULL} /* sentinel */
|
---|
659 | };
|
---|
660 |
|
---|
661 | static PyMemberDef complex_members[] = {
|
---|
662 | {"real", T_DOUBLE, offsetof(PyComplexObject, cval.real), READONLY,
|
---|
663 | "the real part of a complex number"},
|
---|
664 | {"imag", T_DOUBLE, offsetof(PyComplexObject, cval.imag), READONLY,
|
---|
665 | "the imaginary part of a complex number"},
|
---|
666 | {0},
|
---|
667 | };
|
---|
668 |
|
---|
669 | static PyObject *
|
---|
670 | complex_subtype_from_string(PyTypeObject *type, PyObject *v)
|
---|
671 | {
|
---|
672 | const char *s, *start;
|
---|
673 | char *end;
|
---|
674 | double x=0.0, y=0.0, z;
|
---|
675 | int got_re=0, got_im=0, done=0;
|
---|
676 | int digit_or_dot;
|
---|
677 | int sw_error=0;
|
---|
678 | int sign;
|
---|
679 | char buffer[256]; /* For errors */
|
---|
680 | #ifdef Py_USING_UNICODE
|
---|
681 | char s_buffer[256];
|
---|
682 | #endif
|
---|
683 | Py_ssize_t len;
|
---|
684 |
|
---|
685 | if (PyString_Check(v)) {
|
---|
686 | s = PyString_AS_STRING(v);
|
---|
687 | len = PyString_GET_SIZE(v);
|
---|
688 | }
|
---|
689 | #ifdef Py_USING_UNICODE
|
---|
690 | else if (PyUnicode_Check(v)) {
|
---|
691 | if (PyUnicode_GET_SIZE(v) >= (Py_ssize_t)sizeof(s_buffer)) {
|
---|
692 | PyErr_SetString(PyExc_ValueError,
|
---|
693 | "complex() literal too large to convert");
|
---|
694 | return NULL;
|
---|
695 | }
|
---|
696 | if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
|
---|
697 | PyUnicode_GET_SIZE(v),
|
---|
698 | s_buffer,
|
---|
699 | NULL))
|
---|
700 | return NULL;
|
---|
701 | s = s_buffer;
|
---|
702 | len = strlen(s);
|
---|
703 | }
|
---|
704 | #endif
|
---|
705 | else if (PyObject_AsCharBuffer(v, &s, &len)) {
|
---|
706 | PyErr_SetString(PyExc_TypeError,
|
---|
707 | "complex() arg is not a string");
|
---|
708 | return NULL;
|
---|
709 | }
|
---|
710 |
|
---|
711 | /* position on first nonblank */
|
---|
712 | start = s;
|
---|
713 | while (*s && isspace(Py_CHARMASK(*s)))
|
---|
714 | s++;
|
---|
715 | if (s[0] == '\0') {
|
---|
716 | PyErr_SetString(PyExc_ValueError,
|
---|
717 | "complex() arg is an empty string");
|
---|
718 | return NULL;
|
---|
719 | }
|
---|
720 |
|
---|
721 | z = -1.0;
|
---|
722 | sign = 1;
|
---|
723 | do {
|
---|
724 |
|
---|
725 | switch (*s) {
|
---|
726 |
|
---|
727 | case '\0':
|
---|
728 | if (s-start != len) {
|
---|
729 | PyErr_SetString(
|
---|
730 | PyExc_ValueError,
|
---|
731 | "complex() arg contains a null byte");
|
---|
732 | return NULL;
|
---|
733 | }
|
---|
734 | if(!done) sw_error=1;
|
---|
735 | break;
|
---|
736 |
|
---|
737 | case '-':
|
---|
738 | sign = -1;
|
---|
739 | /* Fallthrough */
|
---|
740 | case '+':
|
---|
741 | if (done) sw_error=1;
|
---|
742 | s++;
|
---|
743 | if ( *s=='\0'||*s=='+'||*s=='-' ||
|
---|
744 | isspace(Py_CHARMASK(*s)) ) sw_error=1;
|
---|
745 | break;
|
---|
746 |
|
---|
747 | case 'J':
|
---|
748 | case 'j':
|
---|
749 | if (got_im || done) {
|
---|
750 | sw_error = 1;
|
---|
751 | break;
|
---|
752 | }
|
---|
753 | if (z<0.0) {
|
---|
754 | y=sign;
|
---|
755 | }
|
---|
756 | else{
|
---|
757 | y=sign*z;
|
---|
758 | }
|
---|
759 | got_im=1;
|
---|
760 | s++;
|
---|
761 | if (*s!='+' && *s!='-' )
|
---|
762 | done=1;
|
---|
763 | break;
|
---|
764 |
|
---|
765 | default:
|
---|
766 | if (isspace(Py_CHARMASK(*s))) {
|
---|
767 | while (*s && isspace(Py_CHARMASK(*s)))
|
---|
768 | s++;
|
---|
769 | if (s[0] != '\0')
|
---|
770 | sw_error=1;
|
---|
771 | else
|
---|
772 | done = 1;
|
---|
773 | break;
|
---|
774 | }
|
---|
775 | digit_or_dot =
|
---|
776 | (*s=='.' || isdigit(Py_CHARMASK(*s)));
|
---|
777 | if (done||!digit_or_dot) {
|
---|
778 | sw_error=1;
|
---|
779 | break;
|
---|
780 | }
|
---|
781 | errno = 0;
|
---|
782 | PyFPE_START_PROTECT("strtod", return 0)
|
---|
783 | z = PyOS_ascii_strtod(s, &end) ;
|
---|
784 | PyFPE_END_PROTECT(z)
|
---|
785 | if (errno != 0) {
|
---|
786 | PyOS_snprintf(buffer, sizeof(buffer),
|
---|
787 | "float() out of range: %.150s", s);
|
---|
788 | PyErr_SetString(
|
---|
789 | PyExc_ValueError,
|
---|
790 | buffer);
|
---|
791 | return NULL;
|
---|
792 | }
|
---|
793 | s=end;
|
---|
794 | if (*s=='J' || *s=='j') {
|
---|
795 |
|
---|
796 | break;
|
---|
797 | }
|
---|
798 | if (got_re) {
|
---|
799 | sw_error=1;
|
---|
800 | break;
|
---|
801 | }
|
---|
802 |
|
---|
803 | /* accept a real part */
|
---|
804 | x=sign*z;
|
---|
805 | got_re=1;
|
---|
806 | if (got_im) done=1;
|
---|
807 | z = -1.0;
|
---|
808 | sign = 1;
|
---|
809 | break;
|
---|
810 |
|
---|
811 | } /* end of switch */
|
---|
812 |
|
---|
813 | } while (s - start < len && !sw_error);
|
---|
814 |
|
---|
815 | if (sw_error) {
|
---|
816 | PyErr_SetString(PyExc_ValueError,
|
---|
817 | "complex() arg is a malformed string");
|
---|
818 | return NULL;
|
---|
819 | }
|
---|
820 |
|
---|
821 | return complex_subtype_from_doubles(type, x, y);
|
---|
822 | }
|
---|
823 |
|
---|
824 | static PyObject *
|
---|
825 | complex_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
|
---|
826 | {
|
---|
827 | PyObject *r, *i, *tmp, *f;
|
---|
828 | PyNumberMethods *nbr, *nbi = NULL;
|
---|
829 | Py_complex cr, ci;
|
---|
830 | int own_r = 0;
|
---|
831 | static PyObject *complexstr;
|
---|
832 | static char *kwlist[] = {"real", "imag", 0};
|
---|
833 |
|
---|
834 | r = Py_False;
|
---|
835 | i = NULL;
|
---|
836 | if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OO:complex", kwlist,
|
---|
837 | &r, &i))
|
---|
838 | return NULL;
|
---|
839 |
|
---|
840 | /* Special-case for single argument that is already complex */
|
---|
841 | if (PyComplex_CheckExact(r) && i == NULL &&
|
---|
842 | type == &PyComplex_Type) {
|
---|
843 | /* Note that we can't know whether it's safe to return
|
---|
844 | a complex *subclass* instance as-is, hence the restriction
|
---|
845 | to exact complexes here. */
|
---|
846 | Py_INCREF(r);
|
---|
847 | return r;
|
---|
848 | }
|
---|
849 | if (PyString_Check(r) || PyUnicode_Check(r)) {
|
---|
850 | if (i != NULL) {
|
---|
851 | PyErr_SetString(PyExc_TypeError,
|
---|
852 | "complex() can't take second arg"
|
---|
853 | " if first is a string");
|
---|
854 | return NULL;
|
---|
855 | }
|
---|
856 | return complex_subtype_from_string(type, r);
|
---|
857 | }
|
---|
858 | if (i != NULL && (PyString_Check(i) || PyUnicode_Check(i))) {
|
---|
859 | PyErr_SetString(PyExc_TypeError,
|
---|
860 | "complex() second arg can't be a string");
|
---|
861 | return NULL;
|
---|
862 | }
|
---|
863 |
|
---|
864 | /* XXX Hack to support classes with __complex__ method */
|
---|
865 | if (complexstr == NULL) {
|
---|
866 | complexstr = PyString_InternFromString("__complex__");
|
---|
867 | if (complexstr == NULL)
|
---|
868 | return NULL;
|
---|
869 | }
|
---|
870 | f = PyObject_GetAttr(r, complexstr);
|
---|
871 | if (f == NULL)
|
---|
872 | PyErr_Clear();
|
---|
873 | else {
|
---|
874 | PyObject *args = PyTuple_New(0);
|
---|
875 | if (args == NULL)
|
---|
876 | return NULL;
|
---|
877 | r = PyEval_CallObject(f, args);
|
---|
878 | Py_DECREF(args);
|
---|
879 | Py_DECREF(f);
|
---|
880 | if (r == NULL)
|
---|
881 | return NULL;
|
---|
882 | own_r = 1;
|
---|
883 | }
|
---|
884 | nbr = r->ob_type->tp_as_number;
|
---|
885 | if (i != NULL)
|
---|
886 | nbi = i->ob_type->tp_as_number;
|
---|
887 | if (nbr == NULL || nbr->nb_float == NULL ||
|
---|
888 | ((i != NULL) && (nbi == NULL || nbi->nb_float == NULL))) {
|
---|
889 | PyErr_SetString(PyExc_TypeError,
|
---|
890 | "complex() argument must be a string or a number");
|
---|
891 | if (own_r) {
|
---|
892 | Py_DECREF(r);
|
---|
893 | }
|
---|
894 | return NULL;
|
---|
895 | }
|
---|
896 | if (PyComplex_Check(r)) {
|
---|
897 | /* Note that if r is of a complex subtype, we're only
|
---|
898 | retaining its real & imag parts here, and the return
|
---|
899 | value is (properly) of the builtin complex type. */
|
---|
900 | cr = ((PyComplexObject*)r)->cval;
|
---|
901 | if (own_r) {
|
---|
902 | Py_DECREF(r);
|
---|
903 | }
|
---|
904 | }
|
---|
905 | else {
|
---|
906 | tmp = PyNumber_Float(r);
|
---|
907 | if (own_r) {
|
---|
908 | Py_DECREF(r);
|
---|
909 | }
|
---|
910 | if (tmp == NULL)
|
---|
911 | return NULL;
|
---|
912 | if (!PyFloat_Check(tmp)) {
|
---|
913 | PyErr_SetString(PyExc_TypeError,
|
---|
914 | "float(r) didn't return a float");
|
---|
915 | Py_DECREF(tmp);
|
---|
916 | return NULL;
|
---|
917 | }
|
---|
918 | cr.real = PyFloat_AsDouble(tmp);
|
---|
919 | Py_DECREF(tmp);
|
---|
920 | cr.imag = 0.0;
|
---|
921 | }
|
---|
922 | if (i == NULL) {
|
---|
923 | ci.real = 0.0;
|
---|
924 | ci.imag = 0.0;
|
---|
925 | }
|
---|
926 | else if (PyComplex_Check(i))
|
---|
927 | ci = ((PyComplexObject*)i)->cval;
|
---|
928 | else {
|
---|
929 | tmp = (*nbi->nb_float)(i);
|
---|
930 | if (tmp == NULL)
|
---|
931 | return NULL;
|
---|
932 | ci.real = PyFloat_AsDouble(tmp);
|
---|
933 | Py_DECREF(tmp);
|
---|
934 | ci.imag = 0.;
|
---|
935 | }
|
---|
936 | cr.real -= ci.imag;
|
---|
937 | cr.imag += ci.real;
|
---|
938 | return complex_subtype_from_c_complex(type, cr);
|
---|
939 | }
|
---|
940 |
|
---|
941 | PyDoc_STRVAR(complex_doc,
|
---|
942 | "complex(real[, imag]) -> complex number\n"
|
---|
943 | "\n"
|
---|
944 | "Create a complex number from a real part and an optional imaginary part.\n"
|
---|
945 | "This is equivalent to (real + imag*1j) where imag defaults to 0.");
|
---|
946 |
|
---|
947 | static PyNumberMethods complex_as_number = {
|
---|
948 | (binaryfunc)complex_add, /* nb_add */
|
---|
949 | (binaryfunc)complex_sub, /* nb_subtract */
|
---|
950 | (binaryfunc)complex_mul, /* nb_multiply */
|
---|
951 | (binaryfunc)complex_classic_div, /* nb_divide */
|
---|
952 | (binaryfunc)complex_remainder, /* nb_remainder */
|
---|
953 | (binaryfunc)complex_divmod, /* nb_divmod */
|
---|
954 | (ternaryfunc)complex_pow, /* nb_power */
|
---|
955 | (unaryfunc)complex_neg, /* nb_negative */
|
---|
956 | (unaryfunc)complex_pos, /* nb_positive */
|
---|
957 | (unaryfunc)complex_abs, /* nb_absolute */
|
---|
958 | (inquiry)complex_nonzero, /* nb_nonzero */
|
---|
959 | 0, /* nb_invert */
|
---|
960 | 0, /* nb_lshift */
|
---|
961 | 0, /* nb_rshift */
|
---|
962 | 0, /* nb_and */
|
---|
963 | 0, /* nb_xor */
|
---|
964 | 0, /* nb_or */
|
---|
965 | complex_coerce, /* nb_coerce */
|
---|
966 | complex_int, /* nb_int */
|
---|
967 | complex_long, /* nb_long */
|
---|
968 | complex_float, /* nb_float */
|
---|
969 | 0, /* nb_oct */
|
---|
970 | 0, /* nb_hex */
|
---|
971 | 0, /* nb_inplace_add */
|
---|
972 | 0, /* nb_inplace_subtract */
|
---|
973 | 0, /* nb_inplace_multiply*/
|
---|
974 | 0, /* nb_inplace_divide */
|
---|
975 | 0, /* nb_inplace_remainder */
|
---|
976 | 0, /* nb_inplace_power */
|
---|
977 | 0, /* nb_inplace_lshift */
|
---|
978 | 0, /* nb_inplace_rshift */
|
---|
979 | 0, /* nb_inplace_and */
|
---|
980 | 0, /* nb_inplace_xor */
|
---|
981 | 0, /* nb_inplace_or */
|
---|
982 | (binaryfunc)complex_int_div, /* nb_floor_divide */
|
---|
983 | (binaryfunc)complex_div, /* nb_true_divide */
|
---|
984 | 0, /* nb_inplace_floor_divide */
|
---|
985 | 0, /* nb_inplace_true_divide */
|
---|
986 | };
|
---|
987 |
|
---|
988 | PyTypeObject PyComplex_Type = {
|
---|
989 | PyObject_HEAD_INIT(&PyType_Type)
|
---|
990 | 0,
|
---|
991 | "complex",
|
---|
992 | sizeof(PyComplexObject),
|
---|
993 | 0,
|
---|
994 | complex_dealloc, /* tp_dealloc */
|
---|
995 | (printfunc)complex_print, /* tp_print */
|
---|
996 | 0, /* tp_getattr */
|
---|
997 | 0, /* tp_setattr */
|
---|
998 | 0, /* tp_compare */
|
---|
999 | (reprfunc)complex_repr, /* tp_repr */
|
---|
1000 | &complex_as_number, /* tp_as_number */
|
---|
1001 | 0, /* tp_as_sequence */
|
---|
1002 | 0, /* tp_as_mapping */
|
---|
1003 | (hashfunc)complex_hash, /* tp_hash */
|
---|
1004 | 0, /* tp_call */
|
---|
1005 | (reprfunc)complex_str, /* tp_str */
|
---|
1006 | PyObject_GenericGetAttr, /* tp_getattro */
|
---|
1007 | 0, /* tp_setattro */
|
---|
1008 | 0, /* tp_as_buffer */
|
---|
1009 | Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags */
|
---|
1010 | complex_doc, /* tp_doc */
|
---|
1011 | 0, /* tp_traverse */
|
---|
1012 | 0, /* tp_clear */
|
---|
1013 | complex_richcompare, /* tp_richcompare */
|
---|
1014 | 0, /* tp_weaklistoffset */
|
---|
1015 | 0, /* tp_iter */
|
---|
1016 | 0, /* tp_iternext */
|
---|
1017 | complex_methods, /* tp_methods */
|
---|
1018 | complex_members, /* tp_members */
|
---|
1019 | 0, /* tp_getset */
|
---|
1020 | 0, /* tp_base */
|
---|
1021 | 0, /* tp_dict */
|
---|
1022 | 0, /* tp_descr_get */
|
---|
1023 | 0, /* tp_descr_set */
|
---|
1024 | 0, /* tp_dictoffset */
|
---|
1025 | 0, /* tp_init */
|
---|
1026 | PyType_GenericAlloc, /* tp_alloc */
|
---|
1027 | complex_new, /* tp_new */
|
---|
1028 | PyObject_Del, /* tp_free */
|
---|
1029 | };
|
---|
1030 |
|
---|
1031 | #endif
|
---|