/* Complex object implementation */
/* Borrows heavily from floatobject.c */
/* Submitted by Jim Hugunin */
#include "Python.h"
#include "structmember.h"
/*[clinic input]
class complex "PyComplexObject *" "&PyComplex_Type"
[clinic start generated code]*/
/*[clinic end generated code: output=da39a3ee5e6b4b0d input=819e057d2d10f5ec]*/
#include "clinic/complexobject.c.h"
/* elementary operations on complex numbers */
static Py_complex c_1 = {1., 0.};
Py_complex
_Py_c_sum(Py_complex a, Py_complex b)
{
Py_complex r;
r.real = a.real + b.real;
r.imag = a.imag + b.imag;
return r;
}
Py_complex
_Py_c_diff(Py_complex a, Py_complex b)
{
Py_complex r;
r.real = a.real - b.real;
r.imag = a.imag - b.imag;
return r;
}
Py_complex
_Py_c_neg(Py_complex a)
{
Py_complex r;
r.real = -a.real;
r.imag = -a.imag;
return r;
}
Py_complex
_Py_c_prod(Py_complex a, Py_complex b)
{
Py_complex r;
r.real = a.real*b.real - a.imag*b.imag;
r.imag = a.real*b.imag + a.imag*b.real;
return r;
}
Py_complex
_Py_c_quot(Py_complex a, Py_complex b)
{
/******************************************************************
This was the original algorithm. It's grossly prone to spurious
overflow and underflow errors. It also merrily divides by 0 despite
checking for that(!). The code still serves a doc purpose here, as
the algorithm following is a simple by-cases transformation of this
one:
Py_complex r;
double d = b.real*b.real + b.imag*b.imag;
if (d == 0.)
errno = EDOM;
r.real = (a.real*b.real + a.imag*b.imag)/d;
r.imag = (a.imag*b.real - a.real*b.imag)/d;
return r;
******************************************************************/
/* This algorithm is better, and is pretty obvious: first divide the
* numerators and denominator by whichever of {b.real, b.imag} has
* larger magnitude. The earliest reference I found was to CACM
* Algorithm 116 (Complex Division, Robert L. Smith, Stanford
* University). As usual, though, we're still ignoring all IEEE
* endcases.
*/
Py_complex r; /* the result */
const double abs_breal = b.real < 0 ? -b.real : b.real;
const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
if (abs_breal >= abs_bimag) {
/* divide tops and bottom by b.real */
if (abs_breal == 0.0) {
errno = EDOM;
r.real = r.imag = 0.0;
}
else {
const double ratio = b.imag / b.real;
const double denom = b.real + b.imag * ratio;
r.real = (a.real + a.imag * ratio) / denom;
r.imag = (a.imag - a.real * ratio) / denom;
}
}
else if (abs_bimag >= abs_breal) {
/* divide tops and bottom by b.imag */
const double ratio = b.real / b.imag;
const double denom = b.real * ratio + b.imag;
assert(b.imag != 0.0);
r.real = (a.real * ratio + a.imag) / denom;
r.imag = (a.imag * ratio - a.real) / denom;
}
else {
/* At least one of b.real or b.imag is a NaN */
r.real = r.imag = Py_NAN;
}
return r;
}
Py_complex
_Py_c_pow(Py_complex a, Py_complex b)
{
Py_complex r;
double vabs,len,at,phase;
if (b.real == 0. && b.imag == 0.) {
r.real = 1.;
r.imag = 0.;
}
else if (a.real == 0. && a.imag == 0.) {
if (b.imag != 0. || b.real < 0.)
errno = EDOM;
r.real = 0.;
r.imag = 0.;
}
else {
vabs = hypot(a.real,a.imag);
len = pow(vabs,b.real);
at = atan2(a.imag, a.real);
phase = at*b.real;
if (b.imag != 0.0) {
len /= exp(at*b.imag);
phase += b.imag*log(vabs);
}
r.real = len*cos(phase);
r.imag = len*sin(phase);
}
return r;
}
static Py_complex
c_powu(Py_complex x, long n)
{
Py_complex r, p;
long mask = 1;
r = c_1;
p = x;
while (mask > 0 && n >= mask) {
if (n & mask)
r = _Py_c_prod(r,p);
mask <<= 1;
p = _Py_c_prod(p,p);
}
return r;
}
static Py_complex
c_powi(Py_complex x, long n)
{
Py_complex cn;
if (n > 100 || n < -100) {
cn.real = (double) n;
cn.imag = 0.;
return _Py_c_pow(x,cn);
}
else if (n > 0)
return c_powu(x,n);
else
return _Py_c_quot(c_1, c_powu(x,-n));
}
double
_Py_c_abs(Py_complex z)
{
/* sets errno = ERANGE on overflow; otherwise errno = 0 */
double result;
if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
/* C99 rules: if either the real or the imaginary part is an
infinity, return infinity, even if the other part is a
NaN. */
if (Py_IS_INFINITY(z.real)) {
result = fabs(z.real);
errno = 0;
return result;
}
if (Py_IS_INFINITY(z.imag)) {
result = fabs(z.imag);
errno = 0;
return result;
}
/* either the real or imaginary part is a NaN,
and neither is infinite. Result should be NaN. */
return Py_NAN;
}
result = hypot(z.real, z.imag);
if (!Py_IS_FINITE(result))
errno = ERANGE;
else
errno = 0;
return result;
}
static PyObject *
complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
{
PyObject *op;
op = type->tp_alloc(type, 0);
if (op != NULL)
((PyComplexObject *)op)->cval = cval;
return op;
}
PyObject *
PyComplex_FromCComplex(Py_complex cval)
{
PyComplexObject *op;
/* Inline PyObject_New */
op = (PyComplexObject *) PyObject_MALLOC(sizeof(PyComplexObject));
if (op == NULL)
return PyErr_NoMemory();
PyObject_INIT(op, &PyComplex_Type);
op->cval = cval;
return (PyObject *) op;
}
static PyObject *
complex_subtype_from_doubles(PyTypeObject *type, double real, double imag)
{
Py_complex c;
c.real = real;
c.imag = imag;
return complex_subtype_from_c_complex(type, c);
}
PyObject *
PyComplex_FromDoubles(double real, double imag)
{
Py_complex c;
c.real = real;
c.imag = imag;
return PyComplex_FromCComplex(c);
}
double
PyComplex_RealAsDouble(PyObject *op)
{
if (PyComplex_Check(op)) {
return ((PyComplexObject *)op)->cval.real;
}
else {
return PyFloat_AsDouble(op);
}
}
double
PyComplex_ImagAsDouble(PyObject *op)
{
if (PyComplex_Check(op)) {
return ((PyComplexObject *)op)->cval.imag;
}
else {
return 0.0;
}
}
static PyObject *
try_complex_special_method(PyObject *op)
{
PyObject *f;
_Py_IDENTIFIER(__complex__);
f = _PyObject_LookupSpecial(op, &PyId___complex__);
if (f) {
PyObject *res = _PyObject_CallNoArg(f);
Py_DECREF(f);
if (!res || PyComplex_CheckExact(res)) {
return res;
}
if (!PyComplex_Check(res)) {
PyErr_Format(PyExc_TypeError,
"__complex__ returned non-complex (type %.200s)",
res->ob_type->tp_name);
Py_DECREF(res);
return NULL;
}
/* Issue #29894: warn if 'res' not of exact type complex. */
if (PyErr_WarnFormat(PyExc_DeprecationWarning, 1,
"__complex__ returned non-complex (type %.200s). "
"The ability to return an instance of a strict subclass of complex "
"is deprecated, and may be removed in a future version of Python.",
res->ob_type->tp_name)) {
Py_DECREF(res);
return NULL;
}
return res;
}
return NULL;
}
Py_complex
PyComplex_AsCComplex(PyObject *op)
{
Py_complex cv;
PyObject *newop = NULL;
assert(op);
/* If op is already of type PyComplex_Type, return its value */
if (PyComplex_Check(op)) {
return ((PyComplexObject *)op)->cval;
}
/* If not, use op's __complex__ method, if it exists */
/* return -1 on failure */
cv.real = -1.;
cv.imag = 0.;
newop = try_complex_special_method(op);
if (newop) {
cv = ((PyComplexObject *)newop)->cval;
Py_DECREF(newop);
return cv;
}
else if (PyErr_Occurred()) {
return cv;
}
/* If neither of the above works, interpret op as a float giving the
real part of the result, and fill in the imaginary part as 0. */
else {
/* PyFloat_AsDouble will return -1 on failure */
cv.real = PyFloat_AsDouble(op);
return cv;
}
}
static void
complex_dealloc(PyObject *op)
{
op->ob_type->tp_free(op);
}
static PyObject *
complex_repr(PyComplexObject *v)
{
int precision = 0;
char format_code = 'r';
PyObject *result = NULL;
/* If these are non-NULL, they'll need to be freed. */
char *pre = NULL;
char *im = NULL;
/* These do not need to be freed. re is either an alias
for pre or a pointer to a constant. lead and tail
are pointers to constants. */
const char *re = NULL;
const char *lead = "";
const char *tail = "";
if (v->cval.real == 0. && copysign(1.0, v->cval.real)==1.0) {
/* Real part is +0: just output the imaginary part and do not
include parens. */
re = "";
im = PyOS_double_to_string(v->cval.imag, format_code,
precision, 0, NULL);
if (!im) {
PyErr_NoMemory();
goto done;
}
} else {
/* Format imaginary part with sign, real part without. Include
parens in the result. */
pre = PyOS_double_to_string(v->cval.real, format_code,
precision, 0, NULL);
if (!pre) {
PyErr_NoMemory();
goto done;
}
re = pre;
im = PyOS_double_to_string(v->cval.imag, format_code,
precision, Py_DTSF_SIGN, NULL);
if (!im) {
PyErr_NoMemory();
goto done;
}
lead = "(";
tail = ")";
}
result = PyUnicode_FromFormat("%s%s%sj%s", lead, re, im, tail);
done:
PyMem_Free(im);
PyMem_Free(pre);
return result;
}
static Py_hash_t
complex_hash(PyComplexObject *v)
{
Py_uhash_t hashreal, hashimag, combined;
hashreal = (Py_uhash_t)_Py_HashDouble(v->cval.real);
if (hashreal == (Py_uhash_t)-1)
return -1;
hashimag = (Py_uhash_t)_Py_HashDouble(v->cval.imag);
if (hashimag == (Py_uhash_t)-1)
return -1;
/* Note: if the imaginary part is 0, hashimag is 0 now,
* so the following returns hashreal unchanged. This is
* important because numbers of different types that
* compare equal must have the same hash value, so that
* hash(x + 0*j) must equal hash(x).
*/
combined = hashreal + _PyHASH_IMAG * hashimag;
if (combined == (Py_uhash_t)-1)
combined = (Py_uhash_t)-2;
return (Py_hash_t)combined;
}
/* This macro may return! */
#define TO_COMPLEX(obj, c) \
if (PyComplex_Check(obj)) \
c = ((PyComplexObject *)(obj))->cval; \
else if (to_complex(&(obj), &(c)) < 0) \
return (obj)
static int
to_complex(PyObject **pobj, Py_complex *pc)
{
PyObject *obj = *pobj;
pc->real = pc->imag = 0.0;
if (PyLong_Check(obj)) {
pc->real = PyLong_AsDouble(obj);
if (pc->real == -1.0 && PyErr_Occurred()) {
*pobj = NULL;
return -1;
}
return 0;
}
if (PyFloat_Check(obj)) {
pc->real = PyFloat_AsDouble(obj);
return 0;
}
Py_INCREF(Py_NotImplemented);
*pobj = Py_NotImplemented;
return -1;
}
static PyObject *
complex_add(PyObject *v, PyObject *w)
{
Py_complex result;
Py_complex a, b;
TO_COMPLEX(v, a);
TO_COMPLEX(w, b);
PyFPE_START_PROTECT("complex_add", return 0)
result = _Py_c_sum(a, b);
PyFPE_END_PROTECT(result)
return PyComplex_FromCComplex(result);
}
static PyObject *
complex_sub(PyObject *v, PyObject *w)
{
Py_complex result;
Py_complex a, b;
TO_COMPLEX(v, a);
TO_COMPLEX(w, b);
PyFPE_START_PROTECT("complex_sub", return 0)
result = _Py_c_diff(a, b);
PyFPE_END_PROTECT(result)
return PyComplex_FromCComplex(result);
}
static PyObject *
complex_mul(PyObject *v, PyObject *w)
{
Py_complex result;
Py_complex a, b;
TO_COMPLEX(v, a);
TO_COMPLEX(w, b);
PyFPE_START_PROTECT("complex_mul", return 0)
result = _Py_c_prod(a, b);
PyFPE_END_PROTECT(result)
return PyComplex_FromCComplex(result);
}
static PyObject *
complex_div(PyObject *v, PyObject *w)
{
Py_complex quot;
Py_complex a, b;
TO_COMPLEX(v, a);
TO_COMPLEX(w, b);
PyFPE_START_PROTECT("complex_div", return 0)
errno = 0;
quot = _Py_c_quot(a, b);
PyFPE_END_PROTECT(quot)
if (errno == EDOM) {
PyErr_SetString(PyExc_ZeroDivisionError, "complex division by zero");
return NULL;
}
return PyComplex_FromCComplex(quot);
}
static PyObject *
complex_remainder(PyObject *v, PyObject *w)
{
PyErr_SetString(PyExc_TypeError,
"can't mod complex numbers.");
return NULL;
}
static PyObject *
complex_divmod(PyObject *v, PyObject *w)
{
PyErr_SetString(PyExc_TypeError,
"can't take floor or mod of complex number.");
return NULL;
}
static PyObject *
complex_pow(PyObject *v, PyObject *w, PyObject *z)
{
Py_complex p;
Py_complex exponent;
long int_exponent;
Py_complex a, b;
TO_COMPLEX(v, a);
TO_COMPLEX(w, b);
if (z != Py_None) {
PyErr_SetString(PyExc_ValueError, "complex modulo");
return NULL;
}
PyFPE_START_PROTECT("complex_pow", return 0)
errno = 0;
exponent = b;
int_exponent = (long)exponent.real;
if (exponent.imag == 0. && exponent.real == int_exponent)
p = c_powi(a, int_exponent);
else
p = _Py_c_pow(a, exponent);
PyFPE_END_PROTECT(p)
Py_ADJUST_ERANGE2(p.real, p.imag);
if (errno == EDOM) {
PyErr_SetString(PyExc_ZeroDivisionError,
"0.0 to a negative or complex power");
return NULL;
}
else if (errno == ERANGE) {
PyErr_SetString(PyExc_OverflowError,
"complex exponentiation");
return NULL;
}
return PyComplex_FromCComplex(p);
}
static PyObject *
complex_int_div(PyObject *v, PyObject *w)
{
PyErr_SetString(PyExc_TypeError,
"can't take floor of complex number.");
return NULL;
}
static PyObject *
complex_neg(PyComplexObject *v)
{
Py_complex neg;
neg.real = -v->cval.real;
neg.imag = -v->cval.imag;
return PyComplex_FromCComplex(neg);
}
static PyObject *
complex_pos(PyComplexObject *v)
{
if (PyComplex_CheckExact(v)) {
Py_INCREF(v);
return (PyObject *)v;
}
else
return PyComplex_FromCComplex(v->cval);
}
static PyObject *
complex_abs(PyComplexObject *v)
{
double result;
PyFPE_START_PROTECT("complex_abs", return 0)
result = _Py_c_abs(v->cval);
PyFPE_END_PROTECT(result)
if (errno == ERANGE) {
PyErr_SetString(PyExc_OverflowError,
"absolute value too large");
return NULL;
}
return PyFloat_FromDouble(result);
}
static int
complex_bool(PyComplexObject *v)
{
return v->cval.real != 0.0 || v->cval.imag != 0.0;
}
static PyObject *
complex_richcompare(PyObject *v, PyObject *w, int op)
{
PyObject *res;
Py_complex i;
int equal;
if (op != Py_EQ && op != Py_NE) {
goto Unimplemented;
}
assert(PyComplex_Check(v));
TO_COMPLEX(v, i);
if (PyLong_Check(w)) {
/* Check for 0.0 imaginary part first to avoid the rich
* comparison when possible.
*/
if (i.imag == 0.0) {
PyObject *j, *sub_res;
j = PyFloat_FromDouble(i.real);
if (j == NULL)
return NULL;
sub_res = PyObject_RichCompare(j, w, op);
Py_DECREF(j);
return sub_res;
}
else {
equal = 0;
}
}
else if (PyFloat_Check(w)) {
equal = (i.real == PyFloat_AsDouble(w) && i.imag == 0.0);
}
else if (PyComplex_Check(w)) {
Py_complex j;
TO_COMPLEX(w, j);
equal = (i.real == j.real && i.imag == j.imag);
}
else {
goto Unimplemented;
}
if (equal == (op == Py_EQ))
res = Py_True;
else
res = Py_False;
Py_INCREF(res);
return res;
Unimplemented:
Py_RETURN_NOTIMPLEMENTED;
}
static PyObject *
complex_int(PyObject *v)
{
PyErr_SetString(PyExc_TypeError,
"can't convert complex to int");
return NULL;
}
static PyObject *
complex_float(PyObject *v)
{
PyErr_SetString(PyExc_TypeError,
"can't convert complex to float");
return NULL;
}
static PyObject *
complex_conjugate(PyObject *self, PyObject *Py_UNUSED(ignored))
{
Py_complex c;
c = ((PyComplexObject *)self)->cval;
c.imag = -c.imag;
return PyComplex_FromCComplex(c);
}
PyDoc_STRVAR(complex_conjugate_doc,
"complex.conjugate() -> complex\n"
"\n"
"Return the complex conjugate of its argument. (3-4j).conjugate() == 3+4j.");
static PyObject *
complex_getnewargs(PyComplexObject *v, PyObject *Py_UNUSED(ignored))
{
Py_complex c = v->cval;
return Py_BuildValue("(dd)", c.real, c.imag);
}
PyDoc_STRVAR(complex__format__doc,
"complex.__format__() -> str\n"
"\n"
"Convert to a string according to format_spec.");
static PyObject *
complex__format__(PyObject* self, PyObject* args)
{
PyObject *format_spec;
_PyUnicodeWriter writer;
int ret;
if (!PyArg_ParseTuple(args, "U:__format__", &format_spec))
return NULL;
_PyUnicodeWriter_Init(&writer);
ret = _PyComplex_FormatAdvancedWriter(
&writer,
self,
format_spec, 0, PyUnicode_GET_LENGTH(format_spec));
if (ret == -1) {
_PyUnicodeWriter_Dealloc(&writer);
return NULL;
}
return _PyUnicodeWriter_Finish(&writer);
}
#if 0
static PyObject *
complex_is_finite(PyObject *self)
{
Py_complex c;
c = ((PyComplexObject *)self)->cval;
return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
Py_IS_FINITE(c.imag)));
}
PyDoc_STRVAR(complex_is_finite_doc,
"complex.is_finite() -> bool\n"
"\n"
"Returns True if the real and the imaginary part is finite.");
#endif
static PyMethodDef complex_methods[] = {
{"conjugate", (PyCFunction)complex_conjugate, METH_NOARGS,
complex_conjugate_doc},
#if 0
{"is_finite", (PyCFunction)complex_is_finite, METH_NOARGS,
complex_is_finite_doc},
#endif
{"__getnewargs__", (PyCFunction)complex_getnewargs, METH_NOARGS},
{"__format__", (PyCFunction)complex__format__,
METH_VARARGS, complex__format__doc},
{NULL, NULL} /* sentinel */
};
static PyMemberDef complex_members[] = {
{"real", T_DOUBLE, offsetof(PyComplexObject, cval.real), READONLY,
"the real part of a complex number"},
{"imag", T_DOUBLE, offsetof(PyComplexObject, cval.imag), READONLY,
"the imaginary part of a complex number"},
{0},
};
static PyObject *
complex_from_string_inner(const char *s, Py_ssize_t len, void *type)
{
double x=0.0, y=0.0, z;
int got_bracket=0;
const char *start;
char *end;
/* position on first nonblank */
start = s;
while (Py_ISSPACE(*s))
s++;
if (*s == '(') {
/* Skip over possible bracket from repr(). */
got_bracket = 1;
s++;
while (Py_ISSPACE(*s))
s++;
}
/* a valid complex string usually takes one of the three forms: