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; use int(abs(z))");

   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; use long(abs(z))");

   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; use abs(z)");

   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 != 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.