/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
  35static Py_complex c_1 = {1., 0.};
  36
  37Py_complex
  38c_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
  46Py_complex
  47c_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
  55Py_complex
  56c_neg(Py_complex a)
  57{
  58	Py_complex r;
  59	r.real = -a.real;
  60	r.imag = -a.imag;
  61	return r;
  62}
  63
  64Py_complex
  65c_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
  73Py_complex
  74c_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
 127Py_complex
 128c_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
 157static Py_complex
 158c_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
 173static Py_complex
 174c_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
 190double
 191c_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
 222static PyObject *
 223complex_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
 233PyObject *
 234PyComplex_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
 247static PyObject *
 248complex_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
 256PyObject *
 257PyComplex_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
 265double
 266PyComplex_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
 276double
 277PyComplex_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
 287Py_complex
 288PyComplex_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
 349static void
 350complex_dealloc(PyObject *op)
 351{
 352	op->ob_type->tp_free(op);
 353}
 354
 355
 356static void
 357complex_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
 407static int
 408complex_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
 419static PyObject *
 420complex_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
 427static PyObject *
 428complex_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
 435static long
 436complex_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
 464static int
 465to_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
 492static PyObject *
 493complex_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
 502static PyObject *
 503complex_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
 512static PyObject *
 513complex_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
 522static PyObject *
 523complex_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
 538static PyObject *
 539complex_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
 559static PyObject *
 560complex_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
 582static PyObject *
 583complex_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
 609static PyObject *
 610complex_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
 647static PyObject *
 648complex_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
 666static PyObject *
 667complex_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
 675static PyObject *
 676complex_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
 686static PyObject *
 687complex_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
 703static int
 704complex_nonzero(PyComplexObject *v)
 705{
 706	return v->cval.real != 0.0 || v->cval.imag != 0.0;
 707}
 708
 709static int
 710complex_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
 742static PyObject *
 743complex_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
 784static PyObject *
 785complex_int(PyObject *v)
 786{
 787	PyErr_SetString(PyExc_TypeError,
 788		   "can't convert complex to int");
 789	return NULL;
 790}
 791
 792static PyObject *
 793complex_long(PyObject *v)
 794{
 795	PyErr_SetString(PyExc_TypeError,
 796		   "can't convert complex to long");
 797	return NULL;
 798}
 799
 800static PyObject *
 801complex_float(PyObject *v)
 802{
 803	PyErr_SetString(PyExc_TypeError,
 804		   "can't convert complex to float");
 805	return NULL;
 806}
 807
 808static PyObject *
 809complex_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
 817PyDoc_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
 822static PyObject *
 823complex_getnewargs(PyComplexObject *v)
 824{
 825	Py_complex c = v->cval;
 826	return Py_BuildValue("(dd)", c.real, c.imag);
 827}
 828
 829#if 0
 830static PyObject *
 831complex_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
 839PyDoc_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
 845static 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
 856static 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
 864static PyObject *
 865complex_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
1039static PyObject *
1040complex_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
1184PyDoc_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
1190static 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
1231PyTypeObject 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