l4/pkg/python/contrib/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 #ifdef HAVE_IEEEFP_H
  12 #include <ieeefp.h>
  13 #endif
  14
  15 #ifndef WITHOUT_COMPLEX
  16
  17 /* Precisions used by repr() and str(), respectively.
  18
  19    The repr() precision (17 significant decimal digits) is the minimal number
  20    that is guaranteed to have enough precision so that if the number is read
  21    back in the exact same binary value is recreated.  This is true for IEEE
  22    floating point by design, and also happens to work for all other modern
  23    hardware.
  24
  25    The str() precision is chosen so that in most cases, the rounding noise
  26    created by various operations is suppressed, while giving plenty of
  27    precision for practical use.
  28 */
  29
  30 #define PREC_REPR       17
  31 #define PREC_STR        12
  32
  33 /* elementary operations on complex numbers */
  34
  35 static Py_complex c_1 = {1., 0.};
  36
  37 Py_complex
  38 c_sum(Py_complex a, Py_complex b)
  39 {
  40         Py_complex r;
  41         r.real = a.real + b.real;
  42         r.imag = a.imag + b.imag;
  43         return r;
  44 }
  45
  46 Py_complex
  47 c_diff(Py_complex a, Py_complex b)
  48 {
  49         Py_complex r;
  50         r.real = a.real - b.real;
  51         r.imag = a.imag - b.imag;
  52         return r;
  53 }
  54
  55 Py_complex
  56 c_neg(Py_complex a)
  57 {
  58         Py_complex r;
  59         r.real = -a.real;
  60         r.imag = -a.imag;
  61         return r;
  62 }
  63
  64 Py_complex
  65 c_prod(Py_complex a, Py_complex b)
  66 {
  67         Py_complex r;
  68         r.real = a.real*b.real - a.imag*b.imag;
  69         r.imag = a.real*b.imag + a.imag*b.real;
  70         return r;
  71 }
  72
  73 Py_complex
  74 c_quot(Py_complex a, Py_complex b)
  75 {
  76         /******************************************************************
  77         This was the original algorithm.  It's grossly prone to spurious
  78         overflow and underflow errors.  It also merrily divides by 0 despite
  79         checking for that(!).  The code still serves a doc purpose here, as
  80         the algorithm following is a simple by-cases transformation of this
  81         one:
  82
  83         Py_complex r;
  84         double d = b.real*b.real + b.imag*b.imag;
  85         if (d == 0.)
  86                 errno = EDOM;
  87         r.real = (a.real*b.real + a.imag*b.imag)/d;
  88         r.imag = (a.imag*b.real - a.real*b.imag)/d;
  89         return r;
  90         ******************************************************************/
  91
  92         /* This algorithm is better, and is pretty obvious:  first divide the
  93          * numerators and denominator by whichever of {b.real, b.imag} has
  94          * larger magnitude.  The earliest reference I found was to CACM
  95          * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
  96          * University).  As usual, though, we're still ignoring all IEEE
  97          * endcases.
  98          */
  99          Py_complex r;  /* the result */
 100          const double abs_breal = b.real < 0 ? -b.real : b.real;
 101          const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
 102
 103          if (abs_breal >= abs_bimag) {
 104                 /* divide tops and bottom by b.real */
 105                 if (abs_breal == 0.0) {
 106                         errno = EDOM;
 107                         r.real = r.imag = 0.0;
 108                 }
 109                 else {
 110                         const double ratio = b.imag / b.real;
 111                         const double denom = b.real + b.imag * ratio;
 112                         r.real = (a.real + a.imag * ratio) / denom;
 113                         r.imag = (a.imag - a.real * ratio) / denom;
 114                 }
 115         }
 116         else {
 117                 /* divide tops and bottom by b.imag */
 118                 const double ratio = b.real / b.imag;
 119                 const double denom = b.real * ratio + b.imag;
 120                 assert(b.imag != 0.0);
 121                 r.real = (a.real * ratio + a.imag) / denom;
 122                 r.imag = (a.imag * ratio - a.real) / denom;
 123         }
 124         return r;
 125 }
 126
 127 Py_complex
 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
 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
 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
 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 *
 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 *
 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         PyObject_INIT(op, &PyComplex_Type);
 243         op->cval = cval;
 244         return (PyObject *) op;
 245 }
 246
 247 static PyObject *
 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 *
 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
 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
 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 Py_complex
 288 PyComplex_AsCComplex(PyObject *op)
 289 {
 290         Py_complex cv;
 291         PyObject *newop = NULL;
 292         static PyObject *complex_str = NULL;
 293
 294         assert(op);
 295         /* If op is already of type PyComplex_Type, return its value */
 296         if (PyComplex_Check(op)) {
 297                 return ((PyComplexObject *)op)->cval;
 298         }
 299         /* If not, use op's __complex__  method, if it exists */
 300
 301         /* return -1 on failure */
 302         cv.real = -1.;
 303         cv.imag = 0.;
 304
 305         if (complex_str == NULL) {
 306                 if (!(complex_str = PyString_InternFromString("__complex__")))
 307                         return cv;
 308         }
 309
 310         if (PyInstance_Check(op)) {
 311                 /* this can go away in python 3000 */
 312                 if (PyObject_HasAttr(op, complex_str)) {
 313                         newop = PyObject_CallMethod(op, "__complex__", NULL);
 314                         if (!newop)
 315                                 return cv;
 316                 }
 317                 /* else try __float__ */
 318         } else {
 319                 PyObject *complexfunc;
 320                 complexfunc = _PyType_Lookup(op->ob_type, complex_str);
 321                 /* complexfunc is a borrowed reference */
 322                 if (complexfunc) {
 323                         newop = PyObject_CallFunctionObjArgs(complexfunc, op, NULL);
 324                         if (!newop)
 325                                 return cv;
 326                 }
 327         }
 328
 329         if (newop) {
 330                 if (!PyComplex_Check(newop)) {
 331                         PyErr_SetString(PyExc_TypeError,
 332                                 "__complex__ should return a complex object");
 333                         Py_DECREF(newop);
 334                         return cv;
 335                 }
 336                 cv = ((PyComplexObject *)newop)->cval;
 337                 Py_DECREF(newop);
 338                 return cv;
 339         }
 340         /* If neither of the above works, interpret op as a float giving the
 341            real part of the result, and fill in the imaginary part as 0. */
 342         else {
 343                 /* PyFloat_AsDouble will return -1 on failure */
 344                 cv.real = PyFloat_AsDouble(op);
 345                 return cv;
 346         }
 347 }
 348
 349 static void
 350 complex_dealloc(PyObject *op)
 351 {
 352         op->ob_type->tp_free(op);
 353 }
 354
 355
 356 static void
 357 complex_to_buf(char *buf, int bufsz, PyComplexObject *v, int precision)
 358 {
 359         char format[32];
 360         if (v->cval.real == 0.) {
 361                 if (!Py_IS_FINITE(v->cval.imag)) {
 362                         if (Py_IS_NAN(v->cval.imag))
 363                                 strncpy(buf, "nan*j", 6);
 364                         else if (copysign(1, v->cval.imag) == 1)
 365                                 strncpy(buf, "inf*j", 6);
 366                         else
 367                                 strncpy(buf, "-inf*j", 7);
 368                 }
 369                 else {
 370                         PyOS_snprintf(format, sizeof(format), "%%.%ig", precision);
 371                         PyOS_ascii_formatd(buf, bufsz - 1, format, v->cval.imag);
 372                         strncat(buf, "j", 1);
 373                 }
 374         } else {
 375                 char re[64], im[64];
 376                 /* Format imaginary part with sign, real part without */
 377                 if (!Py_IS_FINITE(v->cval.real)) {
 378                         if (Py_IS_NAN(v->cval.real))
 379                                 strncpy(re, "nan", 4);
 380                         /* else if (copysign(1, v->cval.real) == 1) */
 381                         else if (v->cval.real > 0)
 382                                 strncpy(re, "inf", 4);
 383                         else
 384                                 strncpy(re, "-inf", 5);
 385                 }
 386                 else {
 387                         PyOS_snprintf(format, sizeof(format), "%%.%ig", precision);
 388                         PyOS_ascii_formatd(re, sizeof(re), format, v->cval.real);
 389                 }
 390                 if (!Py_IS_FINITE(v->cval.imag)) {
 391                         if (Py_IS_NAN(v->cval.imag))
 392                                 strncpy(im, "+nan*", 6);
 393                         /* else if (copysign(1, v->cval.imag) == 1) */
 394                         else if (v->cval.imag > 0)
 395                                 strncpy(im, "+inf*", 6);
 396                         else
 397                                 strncpy(im, "-inf*", 6);
 398                 }
 399                 else {
 400                         PyOS_snprintf(format, sizeof(format), "%%+.%ig", precision);
 401                         PyOS_ascii_formatd(im, sizeof(im), format, v->cval.imag);
 402                 }
 403                 PyOS_snprintf(buf, bufsz, "(%s%sj)", re, im);
 404         }
 405 }
 406
 407 static int
 408 complex_print(PyComplexObject *v, FILE *fp, int flags)
 409 {
 410         char buf[100];
 411         complex_to_buf(buf, sizeof(buf), v,
 412                        (flags & Py_PRINT_RAW) ? PREC_STR : PREC_REPR);
 413         Py_BEGIN_ALLOW_THREADS
 414         fputs(buf, fp);
 415         Py_END_ALLOW_THREADS
 416         return 0;
 417 }
 418
 419 static PyObject *
 420 complex_repr(PyComplexObject *v)
 421 {
 422         char buf[100];
 423         complex_to_buf(buf, sizeof(buf), v, PREC_REPR);
 424         return PyString_FromString(buf);
 425 }
 426
 427 static PyObject *
 428 complex_str(PyComplexObject *v)
 429 {
 430         char buf[100];
 431         complex_to_buf(buf, sizeof(buf), v, PREC_STR);
 432         return PyString_FromString(buf);
 433 }
 434
 435 static long
 436 complex_hash(PyComplexObject *v)
 437 {
 438         long hashreal, hashimag, combined;
 439         hashreal = _Py_HashDouble(v->cval.real);
 440         if (hashreal == -1)
 441                 return -1;
 442         hashimag = _Py_HashDouble(v->cval.imag);
 443         if (hashimag == -1)
 444                 return -1;
 445         /* Note:  if the imaginary part is 0, hashimag is 0 now,
 446          * so the following returns hashreal unchanged.  This is
 447          * important because numbers of different types that
 448          * compare equal must have the same hash value, so that
 449          * hash(x + 0*j) must equal hash(x).
 450          */
 451         combined = hashreal + 1000003 * hashimag;
 452         if (combined == -1)
 453                 combined = -2;
 454         return combined;
 455 }
 456
 457 /* This macro may return! */
 458 #define TO_COMPLEX(obj, c) \
 459         if (PyComplex_Check(obj)) \
 460                 c = ((PyComplexObject *)(obj))->cval; \
 461         else if (to_complex(&(obj), &(c)) < 0) \
 462                 return (obj)
 463
 464 static int
 465 to_complex(PyObject **pobj, Py_complex *pc)
 466 {
 467     PyObject *obj = *pobj;
 468
 469     pc->real = pc->imag = 0.0;
 470     if (PyInt_Check(obj)) {
 471         pc->real = PyInt_AS_LONG(obj);
 472         return 0;
 473     }
 474     if (PyLong_Check(obj)) {
 475         pc->real = PyLong_AsDouble(obj);
 476         if (pc->real == -1.0 && PyErr_Occurred()) {
 477             *pobj = NULL;
 478             return -1;
 479         }
 480         return 0;
 481     }
 482     if (PyFloat_Check(obj)) {
 483         pc->real = PyFloat_AsDouble(obj);
 484         return 0;
 485     }
 486     Py_INCREF(Py_NotImplemented);
 487     *pobj = Py_NotImplemented;
 488     return -1;
 489 }
 490
 491
 492 static PyObject *
 493 complex_add(PyComplexObject *v, PyComplexObject *w)
 494 {
 495         Py_complex result;
 496         PyFPE_START_PROTECT("complex_add", return 0)
 497         result = c_sum(v->cval,w->cval);
 498         PyFPE_END_PROTECT(result)
 499         return PyComplex_FromCComplex(result);
 500 }
 501
 502 static PyObject *
 503 complex_sub(PyComplexObject *v, PyComplexObject *w)
 504 {
 505         Py_complex result;
 506         PyFPE_START_PROTECT("complex_sub", return 0)
 507         result = c_diff(v->cval,w->cval);
 508         PyFPE_END_PROTECT(result)
 509         return PyComplex_FromCComplex(result);
 510 }
 511
 512 static PyObject *
 513 complex_mul(PyComplexObject *v, PyComplexObject *w)
 514 {
 515         Py_complex result;
 516         PyFPE_START_PROTECT("complex_mul", return 0)
 517         result = c_prod(v->cval,w->cval);
 518         PyFPE_END_PROTECT(result)
 519         return PyComplex_FromCComplex(result);
 520 }
 521
 522 static PyObject *
 523 complex_div(PyComplexObject *v, PyComplexObject *w)
 524 {
 525         Py_complex quot;
 526
 527         PyFPE_START_PROTECT("complex_div", return 0)
 528         errno = 0;
 529         quot = c_quot(v->cval,w->cval);
 530         PyFPE_END_PROTECT(quot)
 531         if (errno == EDOM) {
 532                 PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
 533                 return NULL;
 534         }
 535         return PyComplex_FromCComplex(quot);
 536 }
 537
 538 static PyObject *
 539 complex_classic_div(PyComplexObject *v, PyComplexObject *w)
 540 {
 541         Py_complex quot;
 542
 543         if (Py_DivisionWarningFlag >= 2 &&
 544             PyErr_Warn(PyExc_DeprecationWarning,
 545                        "classic complex division") < 0)
 546                 return NULL;
 547
 548         PyFPE_START_PROTECT("complex_classic_div", return 0)
 549         errno = 0;
 550         quot = c_quot(v->cval,w->cval);
 551         PyFPE_END_PROTECT(quot)
 552         if (errno == EDOM) {
 553                 PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
 554                 return NULL;
 555         }
 556         return PyComplex_FromCComplex(quot);
 557 }
 558
 559 static PyObject *
 560 complex_remainder(PyComplexObject *v, PyComplexObject *w)
 561 {
 562         Py_complex div, mod;
 563
 564         if (PyErr_Warn(PyExc_DeprecationWarning,
 565                        "complex divmod(), // and % are deprecated") < 0)
 566                 return NULL;
 567
 568         errno = 0;
 569         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 570         if (errno == EDOM) {
 571                 PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
 572                 return NULL;
 573         }
 574         div.real = floor(div.real); /* Use the floor of the real part. */
 575         div.imag = 0.0;
 576         mod = c_diff(v->cval, c_prod(w->cval, div));
 577
 578         return PyComplex_FromCComplex(mod);
 579 }
 580
 581
 582 static PyObject *
 583 complex_divmod(PyComplexObject *v, PyComplexObject *w)
 584 {
 585         Py_complex div, mod;
 586         PyObject *d, *m, *z;
 587
 588         if (PyErr_Warn(PyExc_DeprecationWarning,
 589                        "complex divmod(), // and % are deprecated") < 0)
 590                 return NULL;
 591
 592         errno = 0;
 593         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 594         if (errno == EDOM) {
 595                 PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
 596                 return NULL;
 597         }
 598         div.real = floor(div.real); /* Use the floor of the real part. */
 599         div.imag = 0.0;
 600         mod = c_diff(v->cval, c_prod(w->cval, div));
 601         d = PyComplex_FromCComplex(div);
 602         m = PyComplex_FromCComplex(mod);
 603         z = PyTuple_Pack(2, d, m);
 604         Py_XDECREF(d);
 605         Py_XDECREF(m);
 606         return z;
 607 }
 608
 609 static PyObject *
 610 complex_pow(PyObject *v, PyObject *w, PyObject *z)
 611 {
 612         Py_complex p;
 613         Py_complex exponent;
 614         long int_exponent;
 615         Py_complex a, b;
 616         TO_COMPLEX(v, a);
 617         TO_COMPLEX(w, b);
 618
 619         if (z!=Py_None) {
 620                 PyErr_SetString(PyExc_ValueError, "complex modulo");
 621                 return NULL;
 622         }
 623         PyFPE_START_PROTECT("complex_pow", return 0)
 624         errno = 0;
 625         exponent = b;
 626         int_exponent = (long)exponent.real;
 627         if (exponent.imag == 0. && exponent.real == int_exponent)
 628                 p = c_powi(a,int_exponent);
 629         else
 630                 p = c_pow(a,exponent);
 631
 632         PyFPE_END_PROTECT(p)
 633         Py_ADJUST_ERANGE2(p.real, p.imag);
 634         if (errno == EDOM) {
 635                 PyErr_SetString(PyExc_ZeroDivisionError,
 636                                 "0.0 to a negative or complex power");
 637                 return NULL;
 638         }
 639         else if (errno == ERANGE) {
 640                 PyErr_SetString(PyExc_OverflowError,
 641                                 "complex exponentiation");
 642                 return NULL;
 643         }
 644         return PyComplex_FromCComplex(p);
 645 }
 646
 647 static PyObject *
 648 complex_int_div(PyComplexObject *v, PyComplexObject *w)
 649 {
 650         PyObject *t, *r;
 651
 652         if (PyErr_Warn(PyExc_DeprecationWarning,
 653                        "complex divmod(), // and % are deprecated") < 0)
 654                 return NULL;
 655
 656         t = complex_divmod(v, w);
 657         if (t != NULL) {
 658                 r = PyTuple_GET_ITEM(t, 0);
 659                 Py_INCREF(r);
 660                 Py_DECREF(t);
 661                 return r;
 662         }
 663         return NULL;
 664 }
 665
 666 static PyObject *
 667 complex_neg(PyComplexObject *v)
 668 {
 669         Py_complex neg;
 670         neg.real = -v->cval.real;
 671         neg.imag = -v->cval.imag;
 672         return PyComplex_FromCComplex(neg);
 673 }
 674
 675 static PyObject *
 676 complex_pos(PyComplexObject *v)
 677 {
 678         if (PyComplex_CheckExact(v)) {
 679                 Py_INCREF(v);
 680                 return (PyObject *)v;
 681         }
 682         else
 683                 return PyComplex_FromCComplex(v->cval);
 684 }
 685
 686 static PyObject *
 687 complex_abs(PyComplexObject *v)
 688 {
 689         double result;
 690
 691         PyFPE_START_PROTECT("complex_abs", return 0)
 692         result = c_abs(v->cval);
 693         PyFPE_END_PROTECT(result)
 694
 695         if (errno == ERANGE) {
 696                 PyErr_SetString(PyExc_OverflowError,
 697                                 "absolute value too large");
 698                 return NULL;
 699         }
 700         return PyFloat_FromDouble(result);
 701 }
 702
 703 static int
 704 complex_nonzero(PyComplexObject *v)
 705 {
 706         return v->cval.real != 0.0 || v->cval.imag != 0.0;
 707 }
 708
 709 static int
 710 complex_coerce(PyObject **pv, PyObject **pw)
 711 {
 712         Py_complex cval;
 713         cval.imag = 0.;
 714         if (PyInt_Check(*pw)) {
 715                 cval.real = (double)PyInt_AsLong(*pw);
 716                 *pw = PyComplex_FromCComplex(cval);
 717                 Py_INCREF(*pv);
 718                 return 0;
 719         }
 720         else if (PyLong_Check(*pw)) {
 721                 cval.real = PyLong_AsDouble(*pw);
 722                 if (cval.real == -1.0 && PyErr_Occurred())
 723                         return -1;
 724                 *pw = PyComplex_FromCComplex(cval);
 725                 Py_INCREF(*pv);
 726                 return 0;
 727         }
 728         else if (PyFloat_Check(*pw)) {
 729                 cval.real = PyFloat_AsDouble(*pw);
 730                 *pw = PyComplex_FromCComplex(cval);
 731                 Py_INCREF(*pv);
 732                 return 0;
 733         }
 734         else if (PyComplex_Check(*pw)) {
 735                 Py_INCREF(*pv);
 736                 Py_INCREF(*pw);
 737                 return 0;
 738         }
 739         return 1; /* Can't do it */
 740 }
 741
 742 static PyObject *
 743 complex_richcompare(PyObject *v, PyObject *w, int op)
 744 {
 745         int c;
 746         Py_complex i, j;
 747         PyObject *res;
 748
 749         c = PyNumber_CoerceEx(&v, &w);
 750         if (c < 0)
 751                 return NULL;
 752         if (c > 0) {
 753                 Py_INCREF(Py_NotImplemented);
 754                 return Py_NotImplemented;
 755         }
 756         /* Make sure both arguments are complex. */
 757         if (!(PyComplex_Check(v) && PyComplex_Check(w))) {
 758                 Py_DECREF(v);
 759                 Py_DECREF(w);
 760                 Py_INCREF(Py_NotImplemented);
 761                 return Py_NotImplemented;
 762         }
 763
 764         i = ((PyComplexObject *)v)->cval;
 765         j = ((PyComplexObject *)w)->cval;
 766         Py_DECREF(v);
 767         Py_DECREF(w);
 768
 769         if (op != Py_EQ && op != Py_NE) {
 770                 PyErr_SetString(PyExc_TypeError,
 771                         "no ordering relation is defined for complex numbers");
 772                 return NULL;
 773         }
 774
 775         if ((i.real == j.real && i.imag == j.imag) == (op == Py_EQ))
 776                 res = Py_True;
 777         else
 778                 res = Py_False;
 779
 780         Py_INCREF(res);
 781         return res;
 782 }
 783
 784 static PyObject *
 785 complex_int(PyObject *v)
 786 {
 787         PyErr_SetString(PyExc_TypeError,
 788                    "can't convert complex to int");
 789         return NULL;
 790 }
 791
 792 static PyObject *
 793 complex_long(PyObject *v)
 794 {
 795         PyErr_SetString(PyExc_TypeError,
 796                    "can't convert complex to long");
 797         return NULL;
 798 }
 799
 800 static PyObject *
 801 complex_float(PyObject *v)
 802 {
 803         PyErr_SetString(PyExc_TypeError,
 804                    "can't convert complex to float");
 805         return NULL;
 806 }
 807
 808 static PyObject *
 809 complex_conjugate(PyObject *self)
 810 {
 811         Py_complex c;
 812         c = ((PyComplexObject *)self)->cval;
 813         c.imag = -c.imag;
 814         return PyComplex_FromCComplex(c);
 815 }
 816
 817 PyDoc_STRVAR(complex_conjugate_doc,
 818 "complex.conjugate() -> complex\n"
 819 "\n"
 820 "Returns the complex conjugate of its argument. (3-4j).conjugate() == 3+4j.");
 821
 822 static PyObject *
 823 complex_getnewargs(PyComplexObject *v)
 824 {
 825         Py_complex c = v->cval;
 826         return Py_BuildValue("(dd)", c.real, c.imag);
 827 }
 828
 829 #if 0
 830 static PyObject *
 831 complex_is_finite(PyObject *self)
 832 {
 833         Py_complex c;
 834         c = ((PyComplexObject *)self)->cval;
 835         return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
 836                                       Py_IS_FINITE(c.imag)));
 837 }
 838
 839 PyDoc_STRVAR(complex_is_finite_doc,
 840 "complex.is_finite() -> bool\n"
 841 "\n"
 842 "Returns True if the real and the imaginary part is finite.");
 843 #endif
 844
 845 static PyMethodDef complex_methods[] = {
 846         {"conjugate",   (PyCFunction)complex_conjugate, METH_NOARGS,
 847          complex_conjugate_doc},
 848 #if 0
 849         {"is_finite",   (PyCFunction)complex_is_finite, METH_NOARGS,
 850          complex_is_finite_doc},
 851 #endif
 852         {"__getnewargs__",      (PyCFunction)complex_getnewargs,        METH_NOARGS},
 853         {NULL,          NULL}           /* sentinel */
 854 };
 855
 856 static PyMemberDef complex_members[] = {
 857         {"real", T_DOUBLE, offsetof(PyComplexObject, cval.real), READONLY,
 858          "the real part of a complex number"},
 859         {"imag", T_DOUBLE, offsetof(PyComplexObject, cval.imag), READONLY,
 860          "the imaginary part of a complex number"},
 861         {0},
 862 };
 863
 864 static PyObject *
 865 complex_subtype_from_string(PyTypeObject *type, PyObject *v)
 866 {
 867         const char *s, *start;
 868         char *end;
 869         double x=0.0, y=0.0, z;
 870         int got_re=0, got_im=0, got_bracket=0, done=0;
 871         int digit_or_dot;
 872         int sw_error=0;
 873         int sign;
 874         char buffer[256]; /* For errors */
 875 #ifdef Py_USING_UNICODE
 876         char s_buffer[256];
 877 #endif
 878         Py_ssize_t len;
 879
 880         if (PyString_Check(v)) {
 881                 s = PyString_AS_STRING(v);
 882                 len = PyString_GET_SIZE(v);
 883         }
 884 #ifdef Py_USING_UNICODE
 885         else if (PyUnicode_Check(v)) {
 886                 if (PyUnicode_GET_SIZE(v) >= (Py_ssize_t)sizeof(s_buffer)) {
 887                         PyErr_SetString(PyExc_ValueError,
 888                                  "complex() literal too large to convert");
 889                         return NULL;
 890                 }
 891                 if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),
 892                                             PyUnicode_GET_SIZE(v),
 893                                             s_buffer,
 894                                             NULL))
 895                         return NULL;
 896                 s = s_buffer;
 897                 len = strlen(s);
 898         }
 899 #endif
 900         else if (PyObject_AsCharBuffer(v, &s, &len)) {
 901                 PyErr_SetString(PyExc_TypeError,
 902                                 "complex() arg is not a string");
 903                 return NULL;
 904         }
 905
 906         /* position on first nonblank */
 907         start = s;
 908         while (*s && isspace(Py_CHARMASK(*s)))
 909                 s++;
 910         if (s[0] == '\0') {
 911                 PyErr_SetString(PyExc_ValueError,
 912                                 "complex() arg is an empty string");
 913                 return NULL;
 914         }
 915         if (s[0] == '(') {
 916                 /* Skip over possible bracket from repr(). */
 917                 got_bracket = 1;
 918                 s++;
 919                 while (*s && isspace(Py_CHARMASK(*s)))
 920                         s++;
 921         }
 922
 923         z = -1.0;
 924         sign = 1;
 925         do {
 926
 927                 switch (*s) {
 928
 929                 case '\0':
 930                         if (s-start != len) {
 931                                 PyErr_SetString(
 932                                         PyExc_ValueError,
 933                                         "complex() arg contains a null byte");
 934                                 return NULL;
 935                         }
 936                         if(!done) sw_error=1;
 937                         break;
 938
 939                 case ')':
 940                         if (!got_bracket || !(got_re || got_im)) {
 941                                 sw_error=1;
 942                                 break;
 943                         }
 944                         got_bracket=0;
 945                         done=1;
 946                         s++;
 947                         while (*s && isspace(Py_CHARMASK(*s)))
 948                                 s++;
 949                         if (*s) sw_error=1;
 950                         break;
 951
 952                 case '-':
 953                         sign = -1;
 954                                 /* Fallthrough */
 955                 case '+':
 956                         if (done)  sw_error=1;
 957                         s++;
 958                         if  (  *s=='\0'||*s=='+'||*s=='-'||*s==')'||
 959                                isspace(Py_CHARMASK(*s))  )  sw_error=1;
 960                         break;
 961
 962                 case 'J':
 963                 case 'j':
 964                         if (got_im || done) {
 965                                 sw_error = 1;
 966                                 break;
 967                         }
 968                         if  (z<0.0) {
 969                                 y=sign;
 970                         }
 971                         else{
 972                                 y=sign*z;
 973                         }
 974                         got_im=1;
 975                         s++;
 976                         if  (*s!='+' && *s!='-' )
 977                                 done=1;
 978                         break;
 979
 980                 default:
 981                         if (isspace(Py_CHARMASK(*s))) {
 982                                 while (*s && isspace(Py_CHARMASK(*s)))
 983                                         s++;
 984                                 if (*s && *s != ')')
 985                                         sw_error=1;
 986                                 else
 987                                         done = 1;
 988                                 break;
 989                         }
 990                         digit_or_dot =
 991                                 (*s=='.' || isdigit(Py_CHARMASK(*s)));
 992                         if  (done||!digit_or_dot) {
 993                                 sw_error=1;
 994                                 break;
 995                         }
 996                         errno = 0;
 997                         PyFPE_START_PROTECT("strtod", return 0)
 998                         z = PyOS_ascii_strtod(s, &end) ;
 999                         PyFPE_END_PROTECT(z)
1000                         if (errno == ERANGE && fabs(z) >= 1.0) {
1001                                 PyOS_snprintf(buffer, sizeof(buffer),
1002                                         "float() out of range: %.150s", s);
1003                                 PyErr_SetString(
1004                                         PyExc_ValueError,
1005                                         buffer);
1006                                 return NULL;
1007                         }
1008                         s=end;
1009                         if  (*s=='J' || *s=='j') {
1010
1011                                 break;
1012                         }
1013                         if  (got_re) {
1014                                 sw_error=1;
1015                                 break;
1016                         }
1017
1018                                 /* accept a real part */
1019                         x=sign*z;
1020                         got_re=1;
1021                         if  (got_im)  done=1;
1022                         z = -1.0;
1023                         sign = 1;
1024                         break;
1025
1026                 }  /* end of switch  */
1027
1028         } while (s - start < len && !sw_error);
1029
1030         if (sw_error || got_bracket) {
1031                 PyErr_SetString(PyExc_ValueError,
1032                                 "complex() arg is a malformed string");
1033                 return NULL;
1034         }
1035
1036         return complex_subtype_from_doubles(type, x, y);
1037 }
1038
1039 static PyObject *
1040 complex_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
1041 {
1042         PyObject *r, *i, *tmp, *f;
1043         PyNumberMethods *nbr, *nbi = NULL;
1044         Py_complex cr, ci;
1045         int own_r = 0;
1046         int cr_is_complex = 0;
1047         int ci_is_complex = 0;
1048         static PyObject *complexstr;
1049         static char *kwlist[] = {"real", "imag", 0};
1050
1051         r = Py_False;
1052         i = NULL;
1053         if (!PyArg_ParseTupleAndKeywords(args, kwds, "|OO:complex", kwlist,
1054                                          &r, &i))
1055                 return NULL;
1056
1057         /* Special-case for a single argument when type(arg) is complex. */
1058         if (PyComplex_CheckExact(r) && i == NULL &&
1059             type == &PyComplex_Type) {
1060                 /* Note that we can't know whether it's safe to return
1061                    a complex *subclass* instance as-is, hence the restriction
1062                    to exact complexes here.  If either the input or the
1063                    output is a complex subclass, it will be handled below
1064                    as a non-orthogonal vector.  */
1065                 Py_INCREF(r);
1066                 return r;
1067         }
1068         if (PyString_Check(r) || PyUnicode_Check(r)) {
1069                 if (i != NULL) {
1070                         PyErr_SetString(PyExc_TypeError,
1071                                         "complex() can't take second arg"
1072                                         " if first is a string");
1073                         return NULL;
1074                 }
1075                 return complex_subtype_from_string(type, r);
1076         }
1077         if (i != NULL && (PyString_Check(i) || PyUnicode_Check(i))) {
1078                 PyErr_SetString(PyExc_TypeError,
1079                                 "complex() second arg can't be a string");
1080                 return NULL;
1081         }
1082
1083         /* XXX Hack to support classes with __complex__ method */
1084         if (complexstr == NULL) {
1085                 complexstr = PyString_InternFromString("__complex__");
1086                 if (complexstr == NULL)
1087                         return NULL;
1088         }
1089         f = PyObject_GetAttr(r, complexstr);
1090         if (f == NULL)
1091                 PyErr_Clear();
1092         else {
1093                 PyObject *args = PyTuple_New(0);
1094                 if (args == NULL)
1095                         return NULL;
1096                 r = PyEval_CallObject(f, args);
1097                 Py_DECREF(args);
1098                 Py_DECREF(f);
1099                 if (r == NULL)
1100                         return NULL;
1101                 own_r = 1;
1102         }
1103         nbr = r->ob_type->tp_as_number;
1104         if (i != NULL)
1105                 nbi = i->ob_type->tp_as_number;
1106         if (nbr == NULL || nbr->nb_float == NULL ||
1107             ((i != NULL) && (nbi == NULL || nbi->nb_float == NULL))) {
1108                 PyErr_SetString(PyExc_TypeError,
1109                            "complex() argument must be a string or a number");
1110                 if (own_r) {
1111                         Py_DECREF(r);
1112                 }
1113                 return NULL;
1114         }
1115
1116         /* If we get this far, then the "real" and "imag" parts should
1117            both be treated as numbers, and the constructor should return a
1118            complex number equal to (real + imag*1j).
1119
1120            Note that we do NOT assume the input to already be in canonical
1121            form; the "real" and "imag" parts might themselves be complex
1122            numbers, which slightly complicates the code below. */
1123         if (PyComplex_Check(r)) {
1124                 /* Note that if r is of a complex subtype, we're only
1125                    retaining its real & imag parts here, and the return
1126                    value is (properly) of the builtin complex type. */
1127                 cr = ((PyComplexObject*)r)->cval;
1128                 cr_is_complex = 1;
1129                 if (own_r) {
1130                         Py_DECREF(r);
1131                 }
1132         }
1133         else {
1134                 /* The "real" part really is entirely real, and contributes
1135                    nothing in the imaginary direction.
1136                    Just treat it as a double. */
1137                 tmp = PyNumber_Float(r);
1138                 if (own_r) {
1139                         /* r was a newly created complex number, rather
1140                            than the original "real" argument. */
1141                         Py_DECREF(r);
1142                 }
1143                 if (tmp == NULL)
1144                         return NULL;
1145                 if (!PyFloat_Check(tmp)) {
1146                         PyErr_SetString(PyExc_TypeError,
1147                                         "float(r) didn't return a float");
1148                         Py_DECREF(tmp);
1149                         return NULL;
1150                 }
1151                 cr.real = PyFloat_AsDouble(tmp);
1152                 cr.imag = 0.0; /* Shut up compiler warning */
1153                 Py_DECREF(tmp);
1154         }
1155         if (i == NULL) {
1156                 ci.real = 0.0;
1157         }
1158         else if (PyComplex_Check(i)) {
1159                 ci = ((PyComplexObject*)i)->cval;
1160                 ci_is_complex = 1;
1161         } else {
1162                 /* The "imag" part really is entirely imaginary, and
1163                    contributes nothing in the real direction.
1164                    Just treat it as a double. */
1165                 tmp = (*nbi->nb_float)(i);
1166                 if (tmp == NULL)
1167                         return NULL;
1168                 ci.real = PyFloat_AsDouble(tmp);
1169                 Py_DECREF(tmp);
1170         }
1171         /*  If the input was in canonical form, then the "real" and "imag"
1172             parts are real numbers, so that ci.imag and cr.imag are zero.
1173             We need this correction in case they were not real numbers. */
1174
1175         if (ci_is_complex) {
1176                 cr.real -= ci.imag;
1177         }
1178         if (cr_is_complex) {
1179                 ci.real += cr.imag;
1180         }
1181         return complex_subtype_from_doubles(type, cr.real, ci.real);
1182 }
1183
1184 PyDoc_STRVAR(complex_doc,
1185 "complex(real[, imag]) -> complex number\n"
1186 "\n"
1187 "Create a complex number from a real part and an optional imaginary part.\n"
1188 "This is equivalent to (real + imag*1j) where imag defaults to 0.");
1189
1190 static PyNumberMethods complex_as_number = {
1191         (binaryfunc)complex_add,                /* nb_add */
1192         (binaryfunc)complex_sub,                /* nb_subtract */
1193         (binaryfunc)complex_mul,                /* nb_multiply */
1194         (binaryfunc)complex_classic_div,        /* nb_divide */
1195         (binaryfunc)complex_remainder,          /* nb_remainder */
1196         (binaryfunc)complex_divmod,             /* nb_divmod */
1197         (ternaryfunc)complex_pow,               /* nb_power */
1198         (unaryfunc)complex_neg,                 /* nb_negative */
1199         (unaryfunc)complex_pos,                 /* nb_positive */
1200         (unaryfunc)complex_abs,                 /* nb_absolute */
1201         (inquiry)complex_nonzero,               /* nb_nonzero */
1202         0,                                      /* nb_invert */
1203         0,                                      /* nb_lshift */
1204         0,                                      /* nb_rshift */
1205         0,                                      /* nb_and */
1206         0,                                      /* nb_xor */
1207         0,                                      /* nb_or */
1208         complex_coerce,                         /* nb_coerce */
1209         complex_int,                            /* nb_int */
1210         complex_long,                           /* nb_long */
1211         complex_float,                          /* nb_float */
1212         0,                                      /* nb_oct */
1213         0,                                      /* nb_hex */
1214         0,                                      /* nb_inplace_add */
1215         0,                                      /* nb_inplace_subtract */
1216         0,                                      /* nb_inplace_multiply*/
1217         0,                                      /* nb_inplace_divide */
1218         0,                                      /* nb_inplace_remainder */
1219         0,                                      /* nb_inplace_power */
1220         0,                                      /* nb_inplace_lshift */
1221         0,                                      /* nb_inplace_rshift */
1222         0,                                      /* nb_inplace_and */
1223         0,                                      /* nb_inplace_xor */
1224         0,                                      /* nb_inplace_or */
1225         (binaryfunc)complex_int_div,            /* nb_floor_divide */
1226         (binaryfunc)complex_div,                /* nb_true_divide */
1227         0,                                      /* nb_inplace_floor_divide */
1228         0,                                      /* nb_inplace_true_divide */
1229 };
1230
1231 PyTypeObject PyComplex_Type = {
1232         PyVarObject_HEAD_INIT(&PyType_Type, 0)
1233         "complex",
1234         sizeof(PyComplexObject),
1235         0,
1236         complex_dealloc,                        /* tp_dealloc */
1237         (printfunc)complex_print,               /* tp_print */
1238         0,                                      /* tp_getattr */
1239         0,                                      /* tp_setattr */
1240         0,                                      /* tp_compare */
1241         (reprfunc)complex_repr,                 /* tp_repr */
1242         &complex_as_number,                     /* tp_as_number */
1243         0,                                      /* tp_as_sequence */
1244         0,                                      /* tp_as_mapping */
1245         (hashfunc)complex_hash,                 /* tp_hash */
1246         0,                                      /* tp_call */
1247         (reprfunc)complex_str,                  /* tp_str */
1248         PyObject_GenericGetAttr,                /* tp_getattro */
1249         0,                                      /* tp_setattro */
1250         0,                                      /* tp_as_buffer */
1251         Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /* tp_flags */
1252         complex_doc,                            /* tp_doc */
1253         0,                                      /* tp_traverse */
1254         0,                                      /* tp_clear */
1255         complex_richcompare,                    /* tp_richcompare */
1256         0,                                      /* tp_weaklistoffset */
1257         0,                                      /* tp_iter */
1258         0,                                      /* tp_iternext */
1259         complex_methods,                        /* tp_methods */
1260         complex_members,                        /* tp_members */
1261         0,                                      /* tp_getset */
1262         0,                                      /* tp_base */
1263         0,                                      /* tp_dict */
1264         0,                                      /* tp_descr_get */
1265         0,                                      /* tp_descr_set */
1266         0,                                      /* tp_dictoffset */
1267         0,                                      /* tp_init */
1268         PyType_GenericAlloc,                    /* tp_alloc */
1269         complex_new,                            /* tp_new */
1270         PyObject_Del,                           /* tp_free */
1271 };
1272
1273 #endif