Objects/complexobject.c

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