#include "Python.h"

#include "pycore_call.h" // _PyObject_CallNoArgs()

#include "pycore_complexobject.h" // _PyComplex_FormatAdvancedWriter()

#include "pycore_floatobject.h" // _Py_convert_int_to_double()

#include "pycore_long.h" // _PyLong_GetZero()

#include "pycore_object.h" // _PyObject_Init()

#include "pycore_pymath.h" // _Py_ADJUST_ERANGE2()

#define _PyComplexObject_CAST(op) ((PyComplexObject *)(op))

#include "clinic/complexobject.c.h"

static Py_complex c_1 = {1., 0.};

_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_cr_sum(Py_complex a, double b)

{

Py_complex r = a;

r.real += b;

return r;

}

static inline Py_complex

_Py_rc_sum(double a, Py_complex b)

{

return _Py_cr_sum(b, a);

}

_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_cr_diff(Py_complex a, double b)

{

Py_complex r = a;

r.real -= b;

return r;

}

_Py_rc_diff(double a, Py_complex b)

{

Py_complex r;

r.real = a - b.real;

r.imag = -b.imag;

return r;

}

_Py_c_neg(Py_complex a)

{

Py_complex r;

r.real = -a.real;

r.imag = -a.imag;

return r;

}

_Py_c_prod(Py_complex z, Py_complex w)

{

double a = z.real, b = z.imag, c = w.real, d = w.imag;

double ac = a*c, bd = b*d, ad = a*d, bc = b*c;

Py_complex r = {ac - bd, ad + bc};

/* Recover infinities that computed as nan+nanj. See e.g. the C11,

if (isnan(r.real) && isnan(r.imag)) {

int recalc = 0;

if (isinf(a) || isinf(b)) { /* z is infinite */

/* "Box" the infinity and change nans in the other factor to 0 */

a = copysign(isinf(a) ? 1.0 : 0.0, a);

b = copysign(isinf(b) ? 1.0 : 0.0, b);

if (isnan(c)) {

c = copysign(0.0, c);

}

if (isnan(d)) {

d = copysign(0.0, d);

}

recalc = 1;

}

if (isinf(c) || isinf(d)) { /* w is infinite */

/* "Box" the infinity and change nans in the other factor to 0 */

c = copysign(isinf(c) ? 1.0 : 0.0, c);

d = copysign(isinf(d) ? 1.0 : 0.0, d);

if (isnan(a)) {

a = copysign(0.0, a);

}

if (isnan(b)) {

b = copysign(0.0, b);

}

recalc = 1;

}

if (!recalc && (isinf(ac) || isinf(bd) || isinf(ad) || isinf(bc))) {

/* Recover infinities from overflow by changing nans to 0 */

if (isnan(a)) {

a = copysign(0.0, a);

}

if (isnan(b)) {

b = copysign(0.0, b);

}

if (isnan(c)) {

c = copysign(0.0, c);

}

if (isnan(d)) {

d = copysign(0.0, d);

}

recalc = 1;

}

if (recalc) {

r.real = Py_INFINITY*(a*c - b*d);

r.imag = Py_INFINITY*(a*d + b*c);

}

}

return r;

}

_Py_cr_prod(Py_complex a, double b)

{

Py_complex r = a;

r.real *= b;

r.imag *= b;

return r;

}

static inline Py_complex

_Py_rc_prod(double a, Py_complex b)

{

return _Py_cr_prod(b, a);

}

#ifdef _M_ARM64

#pragma optimize("", off)

_Py_c_quot(Py_complex a, Py_complex b)

{

/******************************************************************

/* This algorithm is better, and is pretty obvious: first divide the

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;

}

/* Recover infinities and zeros that computed as nan+nanj. See e.g.

if (isnan(r.real) && isnan(r.imag)) {

if ((isinf(a.real) || isinf(a.imag))

&& isfinite(b.real) && isfinite(b.imag))

{

const double x = copysign(isinf(a.real) ? 1.0 : 0.0, a.real);

const double y = copysign(isinf(a.imag) ? 1.0 : 0.0, a.imag);

r.real = Py_INFINITY * (x*b.real + y*b.imag);

r.imag = Py_INFINITY * (y*b.real - x*b.imag);

}

else if ((isinf(abs_breal) || isinf(abs_bimag))

&& isfinite(a.real) && isfinite(a.imag))

{

const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);

const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);

r.real = 0.0 * (a.real*x + a.imag*y);

r.imag = 0.0 * (a.imag*x - a.real*y);

}

}

return r;

}

_Py_cr_quot(Py_complex a, double b)

{

Py_complex r = a;

if (b) {

r.real /= b;

r.imag /= b;

}

else {

errno = EDOM;

r.real = r.imag = 0.0;

}

return r;

}

_Py_rc_quot(double a, Py_complex b)

{

Py_complex r;

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) {

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 / denom;

r.imag = (-a * ratio) / denom;

}

}

else if (abs_bimag >= abs_breal) {

const double ratio = b.real / b.imag;

const double denom = b.real * ratio + b.imag;

assert(b.imag != 0.0);

r.real = (a * ratio) / denom;

r.imag = (-a) / denom;

}

else {

r.real = r.imag = Py_NAN;

}

if (isnan(r.real) && isnan(r.imag) && isfinite(a)

&& (isinf(abs_breal) || isinf(abs_bimag)))

{

const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);

const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);

r.real = 0.0 * (a*x);

r.imag = 0.0 * (-a*y);

}

return r;

}

#ifdef _M_ARM64

#pragma optimize("", on)

_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);

_Py_ADJUST_ERANGE2(r.real, r.imag);

}

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)

{

if (n > 0)

return c_powu(x,n);

else

return _Py_c_quot(c_1, c_powu(x,-n));

}

_Py_c_abs(Py_complex z)

{

/* sets errno = ERANGE on overflow; otherwise errno = 0 */

double result;

if (!isfinite(z.real) || !isfinite(z.imag)) {

/* C99 rules: if either the real or the imaginary part is an

if (isinf(z.real)) {

result = fabs(z.real);

errno = 0;

return result;

}

if (isinf(z.imag)) {

result = fabs(z.imag);

errno = 0;

return result;

}

/* either the real or imaginary part is a NaN,

return Py_NAN;

}

result = hypot(z.real, z.imag);

if (!isfinite(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)

{

/* Inline PyObject_New */

PyComplexObject *op = PyObject_Malloc(sizeof(PyComplexObject));

if (op == NULL) {

return PyErr_NoMemory();

}

_PyObject_Init((PyObject*)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);

}

static PyObject * try_complex_special_method(PyObject *);

PyComplex_RealAsDouble(PyObject *op)

{

double real = -1.0;

if (PyComplex_Check(op)) {

real = ((PyComplexObject *)op)->cval.real;

}

else {

PyObject* newop = try_complex_special_method(op);

if (newop) {

real = ((PyComplexObject *)newop)->cval.real;

Py_DECREF(newop);

} else if (!PyErr_Occurred()) {

real = PyFloat_AsDouble(op);

}

}

return real;

}

PyComplex_ImagAsDouble(PyObject *op)

{

double imag = -1.0;

if (PyComplex_Check(op)) {

imag = ((PyComplexObject *)op)->cval.imag;

}

else {

PyObject* newop = try_complex_special_method(op);

if (newop) {

imag = ((PyComplexObject *)newop)->cval.imag;

Py_DECREF(newop);

} else if (!PyErr_Occurred()) {

PyFloat_AsDouble(op);

if (!PyErr_Occurred()) {

imag = 0.0;

}

}

}

return imag;

}

static PyObject *

try_complex_special_method(PyObject *op)

{

PyObject *f;

f = _PyObject_LookupSpecial(op, &_Py_ID(__complex__));

if (f) {

PyObject *res = _PyObject_CallNoArgs(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)",

Py_TYPE(res)->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.",

Py_TYPE(res)->tp_name)) {

Py_DECREF(res);

return NULL;

}

return res;

}

return NULL;

}

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

else {

/* PyFloat_AsDouble will return -1 on failure */

cv.real = PyFloat_AsDouble(op);

return cv;

}

}

static PyObject *

complex_repr(PyObject *op)

{

int precision = 0;

char format_code = 'r';

PyObject *result = NULL;

PyComplexObject *v = _PyComplexObject_CAST(op);

/* 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

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

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

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(PyObject *op)

{

Py_uhash_t hashreal, hashimag, combined;

PyComplexObject *v = _PyComplexObject_CAST(op);

hashreal = (Py_uhash_t)_Py_HashDouble(op, v->cval.real);

if (hashreal == (Py_uhash_t)-1)

return -1;

hashimag = (Py_uhash_t)_Py_HashDouble(op, v->cval.imag);

if (hashimag == (Py_uhash_t)-1)

return -1;

/* Note: if the imaginary part is 0, hashimag is 0 now,

combined = hashreal + _PyHASH_IMAG * hashimag;

if (combined == (Py_uhash_t)-1)

combined = (Py_uhash_t)-2;

return (Py_hash_t)combined;

}

#define TO_COMPLEX(obj, c) \

if (PyComplex_Check(obj)) \

c = ((PyComplexObject *)(obj))->cval; \

else if (real_to_complex(&(obj), &(c)) < 0) \

return (obj)

static int

real_to_double(PyObject **pobj, double *dbl)

{

PyObject *obj = *pobj;

if (PyFloat_Check(obj)) {

*dbl = PyFloat_AS_DOUBLE(obj);

}

else if (_Py_convert_int_to_double(pobj, dbl) < 0) {

return -1;

}

return 0;

}

static int

real_to_complex(PyObject **pobj, Py_complex *pc)

{

pc->imag = 0.0;

return real_to_double(pobj, &(pc->real));

}

#define COMPLEX_BINOP(NAME, FUNC) \

static PyObject * \

complex_##NAME(PyObject *v, PyObject *w) \

{ \

Py_complex a; \

errno = 0; \

if (PyComplex_Check(w)) { \

Py_complex b = ((PyComplexObject *)w)->cval; \

if (PyComplex_Check(v)) { \

a = ((PyComplexObject *)v)->cval; \

a = _Py_c_##FUNC(a, b); \

} \

else if (real_to_double(&v, &a.real) < 0) { \

return v; \

} \

else { \

a = _Py_rc_##FUNC(a.real, b); \

} \

} \

else if (!PyComplex_Check(v)) { \

Py_RETURN_NOTIMPLEMENTED; \

} \

else { \

a = ((PyComplexObject *)v)->cval; \

double b; \

if (real_to_double(&w, &b) < 0) { \

return w; \

} \

a = _Py_cr_##FUNC(a, b); \

} \

if (errno == EDOM) { \

PyErr_SetString(PyExc_ZeroDivisionError, \

"division by zero"); \

return NULL; \

} \

return PyComplex_FromCComplex(a); \

}

COMPLEX_BINOP(add, sum)

COMPLEX_BINOP(mul, prod)

COMPLEX_BINOP(sub, diff)

COMPLEX_BINOP(div, quot)

static PyObject *

complex_pow(PyObject *v, PyObject *w, PyObject *z)

{

Py_complex p;

Py_complex a, b;

TO_COMPLEX(v, a);

TO_COMPLEX(w, b);

if (z != Py_None) {

PyErr_SetString(PyExc_ValueError, "complex modulo");

return NULL;

}

errno = 0;

// Check whether the exponent has a small integer value, and if so use

// a faster and more accurate algorithm.

if (b.imag == 0.0 && b.real == floor(b.real) && fabs(b.real) <= 100.0) {

p = c_powi(a, (long)b.real);

_Py_ADJUST_ERANGE2(p.real, p.imag);

}

else {

p = _Py_c_pow(a, b);

}

if (errno == EDOM) {

PyErr_SetString(PyExc_ZeroDivisionError,

"zero 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_neg(PyObject *op)

{

PyComplexObject *v = _PyComplexObject_CAST(op);

Py_complex neg;

neg.real = -v->cval.real;

neg.imag = -v->cval.imag;

return PyComplex_FromCComplex(neg);

}

static PyObject *

complex_pos(PyObject *op)

{

PyComplexObject *v = _PyComplexObject_CAST(op);

if (PyComplex_CheckExact(v)) {

return Py_NewRef(v);

}

return PyComplex_FromCComplex(v->cval);

}

static PyObject *

complex_abs(PyObject *op)

{

PyComplexObject *v = _PyComplexObject_CAST(op);

double result;

errno = 0;

result = _Py_c_abs(v->cval);

if (errno == ERANGE) {

PyErr_SetString(PyExc_OverflowError,

"absolute value too large");

return NULL;

}

return PyFloat_FromDouble(result);

}

static int

complex_bool(PyObject *op)

{

PyComplexObject *v = _PyComplexObject_CAST(op);

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

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;

return Py_NewRef(res);

Unimplemented:

Py_RETURN_NOTIMPLEMENTED;

}

static PyObject *

complex_conjugate_impl(PyComplexObject *self)

{

Py_complex c = self->cval;

c.imag = -c.imag;

return PyComplex_FromCComplex(c);

}

static PyObject *

complex___getnewargs___impl(PyComplexObject *self)

{

Py_complex c = self->cval;

return Py_BuildValue("(dd)", c.real, c.imag);

}

static PyObject *

complex___format___impl(PyComplexObject *self, PyObject *format_spec)

{

_PyUnicodeWriter writer;

int ret;

_PyUnicodeWriter_Init(&writer);

ret = _PyComplex_FormatAdvancedWriter(

&writer,

(PyObject *)self,

format_spec, 0, PyUnicode_GET_LENGTH(format_spec));

if (ret == -1) {

_PyUnicodeWriter_Dealloc(&writer);

return NULL;

}

return _PyUnicodeWriter_Finish(&writer);

}

static PyObject *

complex___complex___impl(PyComplexObject *self)

{

if (PyComplex_CheckExact(self)) {

return Py_NewRef(self);

}

else {

return PyComplex_FromCComplex(self->cval);

}

}

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:

/* first look for forms starting with <float> */

z = PyOS_string_to_double(s, &end, NULL);

if (z == -1.0 && PyErr_Occurred()) {

if (PyErr_ExceptionMatches(PyExc_ValueError))

PyErr_Clear();

else

return NULL;

}

if (end != s) {

/* all 4 forms starting with <float> land here */

s = end;

if (*s == '+' || *s == '-') {

/* <float><signed-float>j | <float><sign>j */

x = z;

y = PyOS_string_to_double(s, &end, NULL);

if (y == -1.0 && PyErr_Occurred()) {

if (PyErr_ExceptionMatches(PyExc_ValueError))

PyErr_Clear();

else

return NULL;

}

if (end != s)

/* <float><signed-float>j */

s = end;

else {

/* <float><sign>j */

y = *s == '+' ? 1.0 : -1.0;