• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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
c_sum(Py_complex a,Py_complex b)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
c_diff(Py_complex a,Py_complex b)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
c_neg(Py_complex a)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
c_prod(Py_complex a,Py_complex b)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
c_quot(Py_complex a,Py_complex b)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 if (abs_bimag >= abs_breal) {
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     else {
121         /* At least one of b.real or b.imag is a NaN */
122         r.real = r.imag = Py_NAN;
123     }
124     return r;
125 }
126 
127 Py_complex
c_pow(Py_complex a,Py_complex b)128 c_pow(Py_complex a, Py_complex b)
129 {
130     Py_complex r;
131     double vabs,len,at,phase;
132     if (b.real == 0. && b.imag == 0.) {
133         r.real = 1.;
134         r.imag = 0.;
135     }
136     else if (a.real == 0. && a.imag == 0.) {
137         if (b.imag != 0. || b.real < 0.)
138             errno = EDOM;
139         r.real = 0.;
140         r.imag = 0.;
141     }
142     else {
143         vabs = hypot(a.real,a.imag);
144         len = pow(vabs,b.real);
145         at = atan2(a.imag, a.real);
146         phase = at*b.real;
147         if (b.imag != 0.0) {
148             len /= exp(at*b.imag);
149             phase += b.imag*log(vabs);
150         }
151         r.real = len*cos(phase);
152         r.imag = len*sin(phase);
153     }
154     return r;
155 }
156 
157 static Py_complex
c_powu(Py_complex x,long n)158 c_powu(Py_complex x, long n)
159 {
160     Py_complex r, p;
161     long mask = 1;
162     r = c_1;
163     p = x;
164     while (mask > 0 && n >= mask) {
165         if (n & mask)
166             r = c_prod(r,p);
167         mask <<= 1;
168         p = c_prod(p,p);
169     }
170     return r;
171 }
172 
173 static Py_complex
c_powi(Py_complex x,long n)174 c_powi(Py_complex x, long n)
175 {
176     Py_complex cn;
177 
178     if (n > 100 || n < -100) {
179         cn.real = (double) n;
180         cn.imag = 0.;
181         return c_pow(x,cn);
182     }
183     else if (n > 0)
184         return c_powu(x,n);
185     else
186         return c_quot(c_1,c_powu(x,-n));
187 
188 }
189 
190 double
c_abs(Py_complex z)191 c_abs(Py_complex z)
192 {
193     /* sets errno = ERANGE on overflow;  otherwise errno = 0 */
194     double result;
195 
196     if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
197         /* C99 rules: if either the real or the imaginary part is an
198            infinity, return infinity, even if the other part is a
199            NaN. */
200         if (Py_IS_INFINITY(z.real)) {
201             result = fabs(z.real);
202             errno = 0;
203             return result;
204         }
205         if (Py_IS_INFINITY(z.imag)) {
206             result = fabs(z.imag);
207             errno = 0;
208             return result;
209         }
210         /* either the real or imaginary part is a NaN,
211            and neither is infinite. Result should be NaN. */
212         return Py_NAN;
213     }
214     result = hypot(z.real, z.imag);
215     if (!Py_IS_FINITE(result))
216         errno = ERANGE;
217     else
218         errno = 0;
219     return result;
220 }
221 
222 static PyObject *
complex_subtype_from_c_complex(PyTypeObject * type,Py_complex cval)223 complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
224 {
225     PyObject *op;
226 
227     op = type->tp_alloc(type, 0);
228     if (op != NULL)
229         ((PyComplexObject *)op)->cval = cval;
230     return op;
231 }
232 
233 PyObject *
PyComplex_FromCComplex(Py_complex cval)234 PyComplex_FromCComplex(Py_complex cval)
235 {
236     register PyComplexObject *op;
237 
238     /* Inline PyObject_New */
239     op = (PyComplexObject *) PyObject_MALLOC(sizeof(PyComplexObject));
240     if (op == NULL)
241         return PyErr_NoMemory();
242     (void)PyObject_INIT(op, &PyComplex_Type);
243     op->cval = cval;
244     return (PyObject *) op;
245 }
246 
247 static PyObject *
complex_subtype_from_doubles(PyTypeObject * type,double real,double imag)248 complex_subtype_from_doubles(PyTypeObject *type, double real, double imag)
249 {
250     Py_complex c;
251     c.real = real;
252     c.imag = imag;
253     return complex_subtype_from_c_complex(type, c);
254 }
255 
256 PyObject *
PyComplex_FromDoubles(double real,double imag)257 PyComplex_FromDoubles(double real, double imag)
258 {
259     Py_complex c;
260     c.real = real;
261     c.imag = imag;
262     return PyComplex_FromCComplex(c);
263 }
264 
265 double
PyComplex_RealAsDouble(PyObject * op)266 PyComplex_RealAsDouble(PyObject *op)
267 {
268     if (PyComplex_Check(op)) {
269         return ((PyComplexObject *)op)->cval.real;
270     }
271     else {
272         return PyFloat_AsDouble(op);
273     }
274 }
275 
276 double
PyComplex_ImagAsDouble(PyObject * op)277 PyComplex_ImagAsDouble(PyObject *op)
278 {
279     if (PyComplex_Check(op)) {
280         return ((PyComplexObject *)op)->cval.imag;
281     }
282     else {
283         return 0.0;
284     }
285 }
286 
287 static PyObject *
try_complex_special_method(PyObject * op)288 try_complex_special_method(PyObject *op) {
289     PyObject *f;
290     static PyObject *complexstr;
291 
292     if (complexstr == NULL) {
293         complexstr = PyString_InternFromString("__complex__");
294         if (complexstr == NULL)
295             return NULL;
296     }
297     if (PyInstance_Check(op)) {
298         f = PyObject_GetAttr(op, complexstr);
299         if (f == NULL) {
300             if (PyErr_ExceptionMatches(PyExc_AttributeError))
301                 PyErr_Clear();
302             else
303                 return NULL;
304         }
305     }
306     else {
307         f = _PyObject_LookupSpecial(op, "__complex__", &complexstr);
308         if (f == NULL && PyErr_Occurred())
309             return NULL;
310     }
311     if (f != NULL) {
312         PyObject *res = PyObject_CallFunctionObjArgs(f, NULL);
313         Py_DECREF(f);
314         return res;
315     }
316     return NULL;
317 }
318 
319 Py_complex
PyComplex_AsCComplex(PyObject * op)320 PyComplex_AsCComplex(PyObject *op)
321 {
322     Py_complex cv;
323     PyObject *newop = NULL;
324 
325     assert(op);
326     /* If op is already of type PyComplex_Type, return its value */
327     if (PyComplex_Check(op)) {
328         return ((PyComplexObject *)op)->cval;
329     }
330     /* If not, use op's __complex__  method, if it exists */
331 
332     /* return -1 on failure */
333     cv.real = -1.;
334     cv.imag = 0.;
335 
336     newop = try_complex_special_method(op);
337 
338     if (newop) {
339         if (!PyComplex_Check(newop)) {
340             PyErr_SetString(PyExc_TypeError,
341                 "__complex__ should return a complex object");
342             Py_DECREF(newop);
343             return cv;
344         }
345         cv = ((PyComplexObject *)newop)->cval;
346         Py_DECREF(newop);
347         return cv;
348     }
349     else if (PyErr_Occurred()) {
350         return cv;
351     }
352     /* If neither of the above works, interpret op as a float giving the
353        real part of the result, and fill in the imaginary part as 0. */
354     else {
355         /* PyFloat_AsDouble will return -1 on failure */
356         cv.real = PyFloat_AsDouble(op);
357         return cv;
358     }
359 }
360 
361 static void
complex_dealloc(PyObject * op)362 complex_dealloc(PyObject *op)
363 {
364     op->ob_type->tp_free(op);
365 }
366 
367 
368 static PyObject *
complex_format(PyComplexObject * v,int precision,char format_code)369 complex_format(PyComplexObject *v, int precision, char format_code)
370 {
371     PyObject *result = NULL;
372     Py_ssize_t len;
373 
374     /* If these are non-NULL, they'll need to be freed. */
375     char *pre = NULL;
376     char *im = NULL;
377     char *buf = NULL;
378 
379     /* These do not need to be freed. re is either an alias
380        for pre or a pointer to a constant.  lead and tail
381        are pointers to constants. */
382     char *re = NULL;
383     char *lead = "";
384     char *tail = "";
385 
386     if (v->cval.real == 0. && copysign(1.0, v->cval.real)==1.0) {
387         re = "";
388         im = PyOS_double_to_string(v->cval.imag, format_code,
389                                    precision, 0, NULL);
390         if (!im) {
391             PyErr_NoMemory();
392             goto done;
393         }
394     } else {
395         /* Format imaginary part with sign, real part without */
396         pre = PyOS_double_to_string(v->cval.real, format_code,
397                                     precision, 0, NULL);
398         if (!pre) {
399             PyErr_NoMemory();
400             goto done;
401         }
402         re = pre;
403 
404         im = PyOS_double_to_string(v->cval.imag, format_code,
405                                    precision, Py_DTSF_SIGN, NULL);
406         if (!im) {
407             PyErr_NoMemory();
408             goto done;
409         }
410         lead = "(";
411         tail = ")";
412     }
413     /* Alloc the final buffer. Add one for the "j" in the format string,
414        and one for the trailing zero. */
415     len = strlen(lead) + strlen(re) + strlen(im) + strlen(tail) + 2;
416     buf = PyMem_Malloc(len);
417     if (!buf) {
418         PyErr_NoMemory();
419         goto done;
420     }
421     PyOS_snprintf(buf, len, "%s%s%sj%s", lead, re, im, tail);
422     result = PyString_FromString(buf);
423   done:
424     PyMem_Free(im);
425     PyMem_Free(pre);
426     PyMem_Free(buf);
427 
428     return result;
429 }
430 
431 static int
complex_print(PyComplexObject * v,FILE * fp,int flags)432 complex_print(PyComplexObject *v, FILE *fp, int flags)
433 {
434     PyObject *formatv;
435     char *buf;
436     if (flags & Py_PRINT_RAW)
437         formatv = complex_format(v, PyFloat_STR_PRECISION, 'g');
438     else
439         formatv = complex_format(v, 0, 'r');
440     if (formatv == NULL)
441         return -1;
442     buf = PyString_AS_STRING(formatv);
443     Py_BEGIN_ALLOW_THREADS
444     fputs(buf, fp);
445     Py_END_ALLOW_THREADS
446     Py_DECREF(formatv);
447     return 0;
448 }
449 
450 static PyObject *
complex_repr(PyComplexObject * v)451 complex_repr(PyComplexObject *v)
452 {
453     return complex_format(v, 0, 'r');
454 }
455 
456 static PyObject *
complex_str(PyComplexObject * v)457 complex_str(PyComplexObject *v)
458 {
459     return complex_format(v, PyFloat_STR_PRECISION, 'g');
460 }
461 
462 static long
complex_hash(PyComplexObject * v)463 complex_hash(PyComplexObject *v)
464 {
465     long hashreal, hashimag, combined;
466     hashreal = _Py_HashDouble(v->cval.real);
467     if (hashreal == -1)
468         return -1;
469     hashimag = _Py_HashDouble(v->cval.imag);
470     if (hashimag == -1)
471         return -1;
472     /* Note:  if the imaginary part is 0, hashimag is 0 now,
473      * so the following returns hashreal unchanged.  This is
474      * important because numbers of different types that
475      * compare equal must have the same hash value, so that
476      * hash(x + 0*j) must equal hash(x).
477      */
478     combined = hashreal + 1000003 * hashimag;
479     if (combined == -1)
480         combined = -2;
481     return combined;
482 }
483 
484 /* This macro may return! */
485 #define TO_COMPLEX(obj, c) \
486     if (PyComplex_Check(obj)) \
487         c = ((PyComplexObject *)(obj))->cval; \
488     else if (to_complex(&(obj), &(c)) < 0) \
489         return (obj)
490 
491 static int
to_complex(PyObject ** pobj,Py_complex * pc)492 to_complex(PyObject **pobj, Py_complex *pc)
493 {
494     PyObject *obj = *pobj;
495 
496     pc->real = pc->imag = 0.0;
497     if (PyInt_Check(obj)) {
498         pc->real = PyInt_AS_LONG(obj);
499         return 0;
500     }
501     if (PyLong_Check(obj)) {
502         pc->real = PyLong_AsDouble(obj);
503         if (pc->real == -1.0 && PyErr_Occurred()) {
504             *pobj = NULL;
505             return -1;
506         }
507         return 0;
508     }
509     if (PyFloat_Check(obj)) {
510         pc->real = PyFloat_AsDouble(obj);
511         return 0;
512     }
513     Py_INCREF(Py_NotImplemented);
514     *pobj = Py_NotImplemented;
515     return -1;
516 }
517 
518 
519 static PyObject *
complex_add(PyObject * v,PyObject * w)520 complex_add(PyObject *v, PyObject *w)
521 {
522     Py_complex result;
523     Py_complex a, b;
524     TO_COMPLEX(v, a);
525     TO_COMPLEX(w, b);
526     PyFPE_START_PROTECT("complex_add", return 0)
527     result = c_sum(a, b);
528     PyFPE_END_PROTECT(result)
529     return PyComplex_FromCComplex(result);
530 }
531 
532 static PyObject *
complex_sub(PyObject * v,PyObject * w)533 complex_sub(PyObject *v, PyObject *w)
534 {
535     Py_complex result;
536     Py_complex a, b;
537     TO_COMPLEX(v, a);
538     TO_COMPLEX(w, b);;
539     PyFPE_START_PROTECT("complex_sub", return 0)
540     result = c_diff(a, b);
541     PyFPE_END_PROTECT(result)
542     return PyComplex_FromCComplex(result);
543 }
544 
545 static PyObject *
complex_mul(PyObject * v,PyObject * w)546 complex_mul(PyObject *v, PyObject *w)
547 {
548     Py_complex result;
549     Py_complex a, b;
550     TO_COMPLEX(v, a);
551     TO_COMPLEX(w, b);
552     PyFPE_START_PROTECT("complex_mul", return 0)
553     result = c_prod(a, b);
554     PyFPE_END_PROTECT(result)
555     return PyComplex_FromCComplex(result);
556 }
557 
558 static PyObject *
complex_div(PyObject * v,PyObject * w)559 complex_div(PyObject *v, PyObject *w)
560 {
561     Py_complex quot;
562     Py_complex a, b;
563     TO_COMPLEX(v, a);
564     TO_COMPLEX(w, b);
565     PyFPE_START_PROTECT("complex_div", return 0)
566     errno = 0;
567     quot = c_quot(a, b);
568     PyFPE_END_PROTECT(quot)
569     if (errno == EDOM) {
570         PyErr_SetString(PyExc_ZeroDivisionError, "complex division by zero");
571         return NULL;
572     }
573     return PyComplex_FromCComplex(quot);
574 }
575 
576 static PyObject *
complex_classic_div(PyObject * v,PyObject * w)577 complex_classic_div(PyObject *v, PyObject *w)
578 {
579     Py_complex quot;
580     Py_complex a, b;
581     TO_COMPLEX(v, a);
582     TO_COMPLEX(w, b);
583     if (Py_DivisionWarningFlag >= 2 &&
584         PyErr_Warn(PyExc_DeprecationWarning,
585                    "classic complex division") < 0)
586         return NULL;
587 
588     PyFPE_START_PROTECT("complex_classic_div", return 0)
589     errno = 0;
590     quot = c_quot(a, b);
591     PyFPE_END_PROTECT(quot)
592     if (errno == EDOM) {
593         PyErr_SetString(PyExc_ZeroDivisionError, "complex division by zero");
594         return NULL;
595     }
596     return PyComplex_FromCComplex(quot);
597 }
598 
599 static PyObject *
complex_remainder(PyObject * v,PyObject * w)600 complex_remainder(PyObject *v, PyObject *w)
601 {
602     Py_complex div, mod;
603     Py_complex a, b;
604     TO_COMPLEX(v, a);
605     TO_COMPLEX(w, b);
606     if (PyErr_Warn(PyExc_DeprecationWarning,
607                    "complex divmod(), // and % are deprecated") < 0)
608         return NULL;
609 
610     errno = 0;
611     div = c_quot(a, b); /* The raw divisor value. */
612     if (errno == EDOM) {
613         PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
614         return NULL;
615     }
616     div.real = floor(div.real); /* Use the floor of the real part. */
617     div.imag = 0.0;
618     mod = c_diff(a, c_prod(b, div));
619 
620     return PyComplex_FromCComplex(mod);
621 }
622 
623 
624 static PyObject *
complex_divmod(PyObject * v,PyObject * w)625 complex_divmod(PyObject *v, PyObject *w)
626 {
627     Py_complex div, mod;
628     PyObject *d, *m, *z;
629     Py_complex a, b;
630     TO_COMPLEX(v, a);
631     TO_COMPLEX(w, b);
632     if (PyErr_Warn(PyExc_DeprecationWarning,
633                    "complex divmod(), // and % are deprecated") < 0)
634         return NULL;
635 
636     errno = 0;
637     div = c_quot(a, b); /* The raw divisor value. */
638     if (errno == EDOM) {
639         PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
640         return NULL;
641     }
642     div.real = floor(div.real); /* Use the floor of the real part. */
643     div.imag = 0.0;
644     mod = c_diff(a, c_prod(b, div));
645     d = PyComplex_FromCComplex(div);
646     m = PyComplex_FromCComplex(mod);
647     z = PyTuple_Pack(2, d, m);
648     Py_XDECREF(d);
649     Py_XDECREF(m);
650     return z;
651 }
652 
653 static PyObject *
complex_pow(PyObject * v,PyObject * w,PyObject * z)654 complex_pow(PyObject *v, PyObject *w, PyObject *z)
655 {
656     Py_complex p;
657     Py_complex exponent;
658     long int_exponent;
659     Py_complex a, b;
660     TO_COMPLEX(v, a);
661     TO_COMPLEX(w, b);
662     if (z!=Py_None) {
663         PyErr_SetString(PyExc_ValueError, "complex modulo");
664         return NULL;
665     }
666     PyFPE_START_PROTECT("complex_pow", return 0)
667     errno = 0;
668     exponent = b;
669     int_exponent = (long)exponent.real;
670     if (exponent.imag == 0. && exponent.real == int_exponent)
671         p = c_powi(a,int_exponent);
672     else
673         p = c_pow(a,exponent);
674 
675     PyFPE_END_PROTECT(p)
676     Py_ADJUST_ERANGE2(p.real, p.imag);
677     if (errno == EDOM) {
678         PyErr_SetString(PyExc_ZeroDivisionError,
679                         "0.0 to a negative or complex power");
680         return NULL;
681     }
682     else if (errno == ERANGE) {
683         PyErr_SetString(PyExc_OverflowError,
684                         "complex exponentiation");
685         return NULL;
686     }
687     return PyComplex_FromCComplex(p);
688 }
689 
690 static PyObject *
complex_int_div(PyObject * v,PyObject * w)691 complex_int_div(PyObject *v, PyObject *w)
692 {
693     PyObject *t, *r;
694     Py_complex a, b;
695     TO_COMPLEX(v, a);
696     TO_COMPLEX(w, b);
697     if (PyErr_Warn(PyExc_DeprecationWarning,
698                    "complex divmod(), // and % are deprecated") < 0)
699         return NULL;
700 
701     t = complex_divmod(v, w);
702     if (t != NULL) {
703         r = PyTuple_GET_ITEM(t, 0);
704         Py_INCREF(r);
705         Py_DECREF(t);
706         return r;
707     }
708     return NULL;
709 }
710 
711 static PyObject *
complex_neg(PyComplexObject * v)712 complex_neg(PyComplexObject *v)
713 {
714     Py_complex neg;
715     neg.real = -v->cval.real;
716     neg.imag = -v->cval.imag;
717     return PyComplex_FromCComplex(neg);
718 }
719 
720 static PyObject *
complex_pos(PyComplexObject * v)721 complex_pos(PyComplexObject *v)
722 {
723     if (PyComplex_CheckExact(v)) {
724         Py_INCREF(v);
725         return (PyObject *)v;
726     }
727     else
728         return PyComplex_FromCComplex(v->cval);
729 }
730 
731 static PyObject *
complex_abs(PyComplexObject * v)732 complex_abs(PyComplexObject *v)
733 {
734     double result;
735 
736     PyFPE_START_PROTECT("complex_abs", return 0)
737     result = c_abs(v->cval);
738     PyFPE_END_PROTECT(result)
739 
740     if (errno == ERANGE) {
741         PyErr_SetString(PyExc_OverflowError,
742                         "absolute value too large");
743         return NULL;
744     }
745     return PyFloat_FromDouble(result);
746 }
747 
748 static int
complex_nonzero(PyComplexObject * v)749 complex_nonzero(PyComplexObject *v)
750 {
751     return v->cval.real != 0.0 || v->cval.imag != 0.0;
752 }
753 
754 static int
complex_coerce(PyObject ** pv,PyObject ** pw)755 complex_coerce(PyObject **pv, PyObject **pw)
756 {
757     Py_complex cval;
758     cval.imag = 0.;
759     if (PyInt_Check(*pw)) {
760         cval.real = (double)PyInt_AsLong(*pw);
761         *pw = PyComplex_FromCComplex(cval);
762         Py_INCREF(*pv);
763         return 0;
764     }
765     else if (PyLong_Check(*pw)) {
766         cval.real = PyLong_AsDouble(*pw);
767         if (cval.real == -1.0 && PyErr_Occurred())
768             return -1;
769         *pw = PyComplex_FromCComplex(cval);
770         Py_INCREF(*pv);
771         return 0;
772     }
773     else if (PyFloat_Check(*pw)) {
774         cval.real = PyFloat_AsDouble(*pw);
775         *pw = PyComplex_FromCComplex(cval);
776         Py_INCREF(*pv);
777         return 0;
778     }
779     else if (PyComplex_Check(*pw)) {
780         Py_INCREF(*pv);
781         Py_INCREF(*pw);
782         return 0;
783     }
784     return 1; /* Can't do it */
785 }
786 
787 static PyObject *
complex_richcompare(PyObject * v,PyObject * w,int op)788 complex_richcompare(PyObject *v, PyObject *w, int op)
789 {
790     PyObject *res;
791     Py_complex i;
792     int equal;
793 
794     if (op != Py_EQ && op != Py_NE) {
795         /* for backwards compatibility, comparisons with non-numbers return
796          * NotImplemented.  Only comparisons with core numeric types raise
797          * TypeError.
798          */
799         if (_PyAnyInt_Check(w) ||
800             PyFloat_Check(w) || PyComplex_Check(w)) {
801             PyErr_SetString(PyExc_TypeError,
802                             "no ordering relation is defined "
803                             "for complex numbers");
804             return NULL;
805         }
806         goto Unimplemented;
807     }
808 
809     assert(PyComplex_Check(v));
810     TO_COMPLEX(v, i);
811 
812     if (_PyAnyInt_Check(w)) {
813         /* Check for 0.0 imaginary part first to avoid the rich
814          * comparison when possible.
815          */
816         if (i.imag == 0.0) {
817             PyObject *j, *sub_res;
818             j = PyFloat_FromDouble(i.real);
819             if (j == NULL)
820                 return NULL;
821 
822             sub_res = PyObject_RichCompare(j, w, op);
823             Py_DECREF(j);
824             return sub_res;
825         }
826         else {
827             equal = 0;
828         }
829     }
830     else if (PyFloat_Check(w)) {
831         equal = (i.real == PyFloat_AsDouble(w) && i.imag == 0.0);
832     }
833     else if (PyComplex_Check(w)) {
834         Py_complex j;
835 
836         TO_COMPLEX(w, j);
837         equal = (i.real == j.real && i.imag == j.imag);
838     }
839     else {
840         goto Unimplemented;
841     }
842 
843     if (equal == (op == Py_EQ))
844          res = Py_True;
845     else
846          res = Py_False;
847 
848     Py_INCREF(res);
849     return res;
850 
851   Unimplemented:
852     Py_INCREF(Py_NotImplemented);
853     return Py_NotImplemented;
854 }
855 
856 static PyObject *
complex_int(PyObject * v)857 complex_int(PyObject *v)
858 {
859     PyErr_SetString(PyExc_TypeError,
860                "can't convert complex to int");
861     return NULL;
862 }
863 
864 static PyObject *
complex_long(PyObject * v)865 complex_long(PyObject *v)
866 {
867     PyErr_SetString(PyExc_TypeError,
868                "can't convert complex to long");
869     return NULL;
870 }
871 
872 static PyObject *
complex_float(PyObject * v)873 complex_float(PyObject *v)
874 {
875     PyErr_SetString(PyExc_TypeError,
876                "can't convert complex to float");
877     return NULL;
878 }
879 
880 static PyObject *
complex_conjugate(PyObject * self)881 complex_conjugate(PyObject *self)
882 {
883     Py_complex c;
884     c = ((PyComplexObject *)self)->cval;
885     c.imag = -c.imag;
886     return PyComplex_FromCComplex(c);
887 }
888 
889 PyDoc_STRVAR(complex_conjugate_doc,
890 "complex.conjugate() -> complex\n"
891 "\n"
892 "Return the complex conjugate of its argument. (3-4j).conjugate() == 3+4j.");
893 
894 static PyObject *
complex_getnewargs(PyComplexObject * v)895 complex_getnewargs(PyComplexObject *v)
896 {
897     Py_complex c = v->cval;
898     return Py_BuildValue("(dd)", c.real, c.imag);
899 }
900 
901 PyDoc_STRVAR(complex__format__doc,
902 "complex.__format__() -> str\n"
903 "\n"
904 "Convert to a string according to format_spec.");
905 
906 static PyObject *
complex__format__(PyObject * self,PyObject * args)907 complex__format__(PyObject* self, PyObject* args)
908 {
909     PyObject *format_spec;
910 
911     if (!PyArg_ParseTuple(args, "O:__format__", &format_spec))
912         return NULL;
913     if (PyBytes_Check(format_spec))
914         return _PyComplex_FormatAdvanced(self,
915                                          PyBytes_AS_STRING(format_spec),
916                                          PyBytes_GET_SIZE(format_spec));
917     if (PyUnicode_Check(format_spec)) {
918         /* Convert format_spec to a str */
919         PyObject *result;
920         PyObject *str_spec = PyObject_Str(format_spec);
921 
922         if (str_spec == NULL)
923             return NULL;
924 
925         result = _PyComplex_FormatAdvanced(self,
926                                            PyBytes_AS_STRING(str_spec),
927                                            PyBytes_GET_SIZE(str_spec));
928 
929         Py_DECREF(str_spec);
930         return result;
931     }
932     PyErr_SetString(PyExc_TypeError, "__format__ requires str or unicode");
933     return NULL;
934 }
935 
936 #if 0
937 static PyObject *
938 complex_is_finite(PyObject *self)
939 {
940     Py_complex c;
941     c = ((PyComplexObject *)self)->cval;
942     return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
943                                   Py_IS_FINITE(c.imag)));
944 }
945 
946 PyDoc_STRVAR(complex_is_finite_doc,
947 "complex.is_finite() -> bool\n"
948 "\n"
949 "Returns True if the real and the imaginary part is finite.");
950 #endif
951 
952 static PyMethodDef complex_methods[] = {
953     {"conjugate",       (PyCFunction)complex_conjugate, METH_NOARGS,
954      complex_conjugate_doc},
955 #if 0
956     {"is_finite",       (PyCFunction)complex_is_finite, METH_NOARGS,
957      complex_is_finite_doc},
958 #endif
959     {"__getnewargs__",          (PyCFunction)complex_getnewargs,        METH_NOARGS},
960     {"__format__",          (PyCFunction)complex__format__,
961                                        METH_VARARGS, complex__format__doc},
962     {NULL,              NULL}           /* sentinel */
963 };
964 
965 static PyMemberDef complex_members[] = {
966     {"real", T_DOUBLE, offsetof(PyComplexObject, cval.real), READONLY,
967      "the real part of a complex number"},
968     {"imag", T_DOUBLE, offsetof(PyComplexObject, cval.imag), READONLY,
969      "the imaginary part of a complex number"},
970     {0},
971 };
972 
973 static PyObject *
complex_subtype_from_string(PyTypeObject * type,PyObject * v)974 complex_subtype_from_string(PyTypeObject *type, PyObject *v)
975 {
976     const char *s, *start;
977     char *end;
978     double x=0.0, y=0.0, z;
979     int got_bracket=0;
980 #ifdef Py_USING_UNICODE
981     char *s_buffer = NULL;
982 #endif
983     Py_ssize_t len;
984 
985     if (PyString_Check(v)) {
986         s = PyString_AS_STRING(v);
987         len = PyString_GET_SIZE(v);
988     }
989 #ifdef Py_USING_UNICODE
990     else if (PyUnicode_Check(v)) {
991         s_buffer = (char *)PyMem_MALLOC(PyUnicode_GET_SIZE(v)+1);
992         if (s_buffer == NULL)
993             return PyErr_NoMemory();
994         if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
995                                     PyUnicode_GET_SIZE(v),
996                                     s_buffer,
997                                     NULL))
998             goto error;
999         s = s_buffer;
1000         len = strlen(s);
1001     }
1002 #endif
1003     else {
1004         PyErr_SetString(PyExc_TypeError,
1005                         "complex() arg is not a string");
1006         return NULL;
1007     }
1008 
1009     /* position on first nonblank */
1010     start = s;
1011     while (Py_ISSPACE(*s))
1012         s++;
1013     if (*s == '(') {
1014         /* Skip over possible bracket from repr(). */
1015         got_bracket = 1;
1016         s++;
1017         while (Py_ISSPACE(*s))
1018             s++;
1019     }
1020 
1021     /* a valid complex string usually takes one of the three forms:
1022 
1023          <float>                  - real part only
1024          <float>j                 - imaginary part only
1025          <float><signed-float>j   - real and imaginary parts
1026 
1027        where <float> represents any numeric string that's accepted by the
1028        float constructor (including 'nan', 'inf', 'infinity', etc.), and
1029        <signed-float> is any string of the form <float> whose first
1030        character is '+' or '-'.
1031 
1032        For backwards compatibility, the extra forms
1033 
1034          <float><sign>j
1035          <sign>j
1036          j
1037 
1038        are also accepted, though support for these forms may be removed from
1039        a future version of Python.
1040     */
1041 
1042     /* first look for forms starting with <float> */
1043     z = PyOS_string_to_double(s, &end, NULL);
1044     if (z == -1.0 && PyErr_Occurred()) {
1045         if (PyErr_ExceptionMatches(PyExc_ValueError))
1046             PyErr_Clear();
1047         else
1048             goto error;
1049     }
1050     if (end != s) {
1051         /* all 4 forms starting with <float> land here */
1052         s = end;
1053         if (*s == '+' || *s == '-') {
1054             /* <float><signed-float>j | <float><sign>j */
1055             x = z;
1056             y = PyOS_string_to_double(s, &end, NULL);
1057             if (y == -1.0 && PyErr_Occurred()) {
1058                 if (PyErr_ExceptionMatches(PyExc_ValueError))
1059                     PyErr_Clear();
1060                 else
1061                     goto error;
1062             }
1063             if (end != s)
1064                 /* <float><signed-float>j */
1065                 s = end;
1066             else {
1067                 /* <float><sign>j */
1068                 y = *s == '+' ? 1.0 : -1.0;
1069                 s++;
1070             }
1071             if (!(*s == 'j' || *s == 'J'))
1072                 goto parse_error;
1073             s++;
1074         }
1075         else if (*s == 'j' || *s == 'J') {
1076             /* <float>j */
1077             s++;
1078             y = z;
1079         }
1080         else
1081             /* <float> */
1082             x = z;
1083     }
1084     else {
1085         /* not starting with <float>; must be <sign>j or j */
1086         if (*s == '+' || *s == '-') {
1087             /* <sign>j */
1088             y = *s == '+' ? 1.0 : -1.0;
1089             s++;
1090         }
1091         else
1092             /* j */
1093             y = 1.0;
1094         if (!(*s == 'j' || *s == 'J'))
1095             goto parse_error;
1096         s++;
1097     }
1098 
1099     /* trailing whitespace and closing bracket */
1100     while (Py_ISSPACE(*s))
1101         s++;
1102     if (got_bracket) {
1103         /* if there was an opening parenthesis, then the corresponding
1104            closing parenthesis should be right here */
1105         if (*s != ')')
1106             goto parse_error;
1107         s++;
1108         while (Py_ISSPACE(*s))
1109             s++;
1110     }
1111 
1112     /* we should now be at the end of the string */
1113     if (s-start != len)
1114         goto parse_error;
1115 
1116 
1117 #ifdef Py_USING_UNICODE
1118     if (s_buffer)
1119         PyMem_FREE(s_buffer);
1120 #endif
1121     return complex_subtype_from_doubles(type, x, y);
1122 
1123   parse_error:
1124     PyErr_SetString(PyExc_ValueError,
1125                     "complex() arg is a malformed string");
1126   error:
1127 #ifdef Py_USING_UNICODE
1128     if (s_buffer)
1129         PyMem_FREE(s_buffer);
1130 #endif
1131     return NULL;
1132 }
1133 
1134 static PyObject *
complex_new(PyTypeObject * type,PyObject * args,PyObject * kwds)1135 complex_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1136 {
1137     PyObject *r, *i, *tmp;
1138     PyNumberMethods *nbr, *nbi = NULL;
1139     Py_complex cr, ci;
1140     int own_r = 0;
1141     int cr_is_complex = 0;
1142     int ci_is_complex = 0;
1143     static char *kwlist[] = {"real", "imag", 0};
1144 
1145     r = Py_False;
1146     i = NULL;
1147     if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OO:complex", kwlist,
1148                                      &r, &i))
1149         return NULL;
1150 
1151     /* Special-case for a single argument when type(arg) is complex. */
1152     if (PyComplex_CheckExact(r) && i == NULL &&
1153         type == &PyComplex_Type) {
1154         /* Note that we can't know whether it's safe to return
1155            a complex *subclass* instance as-is, hence the restriction
1156            to exact complexes here.  If either the input or the
1157            output is a complex subclass, it will be handled below
1158            as a non-orthogonal vector.  */
1159         Py_INCREF(r);
1160         return r;
1161     }
1162     if (PyString_Check(r) || PyUnicode_Check(r)) {
1163         if (i != NULL) {
1164             PyErr_SetString(PyExc_TypeError,
1165                             "complex() can't take second arg"
1166                             " if first is a string");
1167             return NULL;
1168         }
1169         return complex_subtype_from_string(type, r);
1170     }
1171     if (i != NULL && (PyString_Check(i) || PyUnicode_Check(i))) {
1172         PyErr_SetString(PyExc_TypeError,
1173                         "complex() second arg can't be a string");
1174         return NULL;
1175     }
1176 
1177     tmp = try_complex_special_method(r);
1178     if (tmp) {
1179         r = tmp;
1180         own_r = 1;
1181     }
1182     else if (PyErr_Occurred()) {
1183         return NULL;
1184     }
1185 
1186     nbr = r->ob_type->tp_as_number;
1187     if (i != NULL)
1188         nbi = i->ob_type->tp_as_number;
1189     if (nbr == NULL || nbr->nb_float == NULL ||
1190         ((i != NULL) && (nbi == NULL || nbi->nb_float == NULL))) {
1191         PyErr_SetString(PyExc_TypeError,
1192                    "complex() argument must be a string or a number");
1193         if (own_r) {
1194             Py_DECREF(r);
1195         }
1196         return NULL;
1197     }
1198 
1199     /* If we get this far, then the "real" and "imag" parts should
1200        both be treated as numbers, and the constructor should return a
1201        complex number equal to (real + imag*1j).
1202 
1203        Note that we do NOT assume the input to already be in canonical
1204        form; the "real" and "imag" parts might themselves be complex
1205        numbers, which slightly complicates the code below. */
1206     if (PyComplex_Check(r)) {
1207         /* Note that if r is of a complex subtype, we're only
1208            retaining its real & imag parts here, and the return
1209            value is (properly) of the builtin complex type. */
1210         cr = ((PyComplexObject*)r)->cval;
1211         cr_is_complex = 1;
1212         if (own_r) {
1213             Py_DECREF(r);
1214         }
1215     }
1216     else {
1217         /* The "real" part really is entirely real, and contributes
1218            nothing in the imaginary direction.
1219            Just treat it as a double. */
1220         tmp = PyNumber_Float(r);
1221         if (own_r) {
1222             /* r was a newly created complex number, rather
1223                than the original "real" argument. */
1224             Py_DECREF(r);
1225         }
1226         if (tmp == NULL)
1227             return NULL;
1228         if (!PyFloat_Check(tmp)) {
1229             PyErr_SetString(PyExc_TypeError,
1230                             "float(r) didn't return a float");
1231             Py_DECREF(tmp);
1232             return NULL;
1233         }
1234         cr.real = PyFloat_AsDouble(tmp);
1235         cr.imag = 0.0;
1236         Py_DECREF(tmp);
1237     }
1238     if (i == NULL) {
1239         ci.real = cr.imag;
1240     }
1241     else if (PyComplex_Check(i)) {
1242         ci = ((PyComplexObject*)i)->cval;
1243         ci_is_complex = 1;
1244     } else {
1245         /* The "imag" part really is entirely imaginary, and
1246            contributes nothing in the real direction.
1247            Just treat it as a double. */
1248         tmp = (*nbi->nb_float)(i);
1249         if (tmp == NULL)
1250             return NULL;
1251         ci.real = PyFloat_AsDouble(tmp);
1252         Py_DECREF(tmp);
1253     }
1254     /*  If the input was in canonical form, then the "real" and "imag"
1255         parts are real numbers, so that ci.imag and cr.imag are zero.
1256         We need this correction in case they were not real numbers. */
1257 
1258     if (ci_is_complex) {
1259         cr.real -= ci.imag;
1260     }
1261     if (cr_is_complex && i != NULL) {
1262         ci.real += cr.imag;
1263     }
1264     return complex_subtype_from_doubles(type, cr.real, ci.real);
1265 }
1266 
1267 PyDoc_STRVAR(complex_doc,
1268 "complex(real[, imag]) -> complex number\n"
1269 "\n"
1270 "Create a complex number from a real part and an optional imaginary part.\n"
1271 "This is equivalent to (real + imag*1j) where imag defaults to 0.");
1272 
1273 static PyNumberMethods complex_as_number = {
1274     (binaryfunc)complex_add,                    /* nb_add */
1275     (binaryfunc)complex_sub,                    /* nb_subtract */
1276     (binaryfunc)complex_mul,                    /* nb_multiply */
1277     (binaryfunc)complex_classic_div,            /* nb_divide */
1278     (binaryfunc)complex_remainder,              /* nb_remainder */
1279     (binaryfunc)complex_divmod,                 /* nb_divmod */
1280     (ternaryfunc)complex_pow,                   /* nb_power */
1281     (unaryfunc)complex_neg,                     /* nb_negative */
1282     (unaryfunc)complex_pos,                     /* nb_positive */
1283     (unaryfunc)complex_abs,                     /* nb_absolute */
1284     (inquiry)complex_nonzero,                   /* nb_nonzero */
1285     0,                                          /* nb_invert */
1286     0,                                          /* nb_lshift */
1287     0,                                          /* nb_rshift */
1288     0,                                          /* nb_and */
1289     0,                                          /* nb_xor */
1290     0,                                          /* nb_or */
1291     complex_coerce,                             /* nb_coerce */
1292     complex_int,                                /* nb_int */
1293     complex_long,                               /* nb_long */
1294     complex_float,                              /* nb_float */
1295     0,                                          /* nb_oct */
1296     0,                                          /* nb_hex */
1297     0,                                          /* nb_inplace_add */
1298     0,                                          /* nb_inplace_subtract */
1299     0,                                          /* nb_inplace_multiply*/
1300     0,                                          /* nb_inplace_divide */
1301     0,                                          /* nb_inplace_remainder */
1302     0,                                          /* nb_inplace_power */
1303     0,                                          /* nb_inplace_lshift */
1304     0,                                          /* nb_inplace_rshift */
1305     0,                                          /* nb_inplace_and */
1306     0,                                          /* nb_inplace_xor */
1307     0,                                          /* nb_inplace_or */
1308     (binaryfunc)complex_int_div,                /* nb_floor_divide */
1309     (binaryfunc)complex_div,                    /* nb_true_divide */
1310     0,                                          /* nb_inplace_floor_divide */
1311     0,                                          /* nb_inplace_true_divide */
1312 };
1313 
1314 PyTypeObject PyComplex_Type = {
1315     PyVarObject_HEAD_INIT(&PyType_Type, 0)
1316     "complex",
1317     sizeof(PyComplexObject),
1318     0,
1319     complex_dealloc,                            /* tp_dealloc */
1320     (printfunc)complex_print,                   /* tp_print */
1321     0,                                          /* tp_getattr */
1322     0,                                          /* tp_setattr */
1323     0,                                          /* tp_compare */
1324     (reprfunc)complex_repr,                     /* tp_repr */
1325     &complex_as_number,                         /* tp_as_number */
1326     0,                                          /* tp_as_sequence */
1327     0,                                          /* tp_as_mapping */
1328     (hashfunc)complex_hash,                     /* tp_hash */
1329     0,                                          /* tp_call */
1330     (reprfunc)complex_str,                      /* tp_str */
1331     PyObject_GenericGetAttr,                    /* tp_getattro */
1332     0,                                          /* tp_setattro */
1333     0,                                          /* tp_as_buffer */
1334     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_CHECKTYPES |
1335         Py_TPFLAGS_BASETYPE,                    /* tp_flags */
1336     complex_doc,                                /* tp_doc */
1337     0,                                          /* tp_traverse */
1338     0,                                          /* tp_clear */
1339     complex_richcompare,                        /* tp_richcompare */
1340     0,                                          /* tp_weaklistoffset */
1341     0,                                          /* tp_iter */
1342     0,                                          /* tp_iternext */
1343     complex_methods,                            /* tp_methods */
1344     complex_members,                            /* tp_members */
1345     0,                                          /* tp_getset */
1346     0,                                          /* tp_base */
1347     0,                                          /* tp_dict */
1348     0,                                          /* tp_descr_get */
1349     0,                                          /* tp_descr_set */
1350     0,                                          /* tp_dictoffset */
1351     0,                                          /* tp_init */
1352     PyType_GenericAlloc,                        /* tp_alloc */
1353     complex_new,                                /* tp_new */
1354     PyObject_Del,                               /* tp_free */
1355 };
1356 
1357 #endif
1358