/***********************************************************
Copyright 1991-1995 by Stichting Mathematisch Centrum, Amsterdam,
The Netherlands.

                        All Rights Reserved

Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that copyright notice and this permission notice appear in
supporting documentation, and that the names of Stichting Mathematisch
Centrum or CWI or Corporation for National Research Initiatives or
CNRI not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.

While CWI is the initial source for this software, a modified version
is made available by the Corporation for National Research Initiatives
(CNRI) at the Internet address ftp://ftp.python.org.

STICHTING MATHEMATISCH CENTRUM AND CNRI DISCLAIM ALL WARRANTIES WITH
REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH
CENTRUM OR CNRI BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.

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

/* Complex object implementation */

/* Borrows heavily from floatobject.c */

/* Submitted by Jim Hugunin */

#ifndef WITHOUT_COMPLEX

#include "Python.h"
#include "mymath.h"

#ifdef HAVE_LIMITS_H
#include <limits.h>
#endif


/* elementary operations on complex numbers */

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

Py_complex c_sum(a,b)
	Py_complex a,b;
{
	Py_complex r;
	r.real = a.real + b.real;
	r.imag = a.imag + b.imag;
	return r;
}

Py_complex c_diff(a,b)
	Py_complex a,b;
{
	Py_complex r;
	r.real = a.real - b.real;
	r.imag = a.imag - b.imag;
	return r;
}

Py_complex c_neg(a)
	Py_complex a;
{
	Py_complex r;
	r.real = -a.real;
	r.imag = -a.imag;
	return r;
}

Py_complex c_prod(a,b)
	Py_complex a,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 c_quot(a,b)
	Py_complex a,b;
{
	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;
}

Py_complex c_pow(a,b)
	Py_complex a,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 = ERANGE;
		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(x, n)
	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 = c_prod(r,p);
		mask <<= 1;
		p = c_prod(p,p);
	}
	return r;
}

static Py_complex c_powi(x, n)
	Py_complex x;
	long n;
{
	Py_complex cn;

	if (n > 100 || n < -100) {
		cn.real = (double) n;
		cn.imag = 0.;
		return c_pow(x,cn);
	}
	else if (n > 0)
		return c_powu(x,n);
	else
		return c_quot(c_1,c_powu(x,-n));

}

PyObject *
PyComplex_FromCComplex(cval)
	Py_complex cval;
{
	register PyComplexObject *op =
		(PyComplexObject *) malloc(sizeof(PyComplexObject));
	if (op == NULL)
		return PyErr_NoMemory();
	op->ob_type = &PyComplex_Type;
	op->cval = cval;
	_Py_NewReference(op);
	return (PyObject *) op;
}

PyObject *
PyComplex_FromDoubles(real, imag)
	double real, imag;
{
	Py_complex c;
	c.real = real;
	c.imag = imag;
	return PyComplex_FromCComplex(c);
}

double
PyComplex_RealAsDouble(op)
	PyObject *op;
{
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval.real;
	} else {
		return PyFloat_AsDouble(op);
	}
}

double
PyComplex_ImagAsDouble(op)
	PyObject *op;
{
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval.imag;
	} else {
		return 0.0;
	}
}

Py_complex
PyComplex_AsCComplex(op)
	PyObject *op;
{
	Py_complex cv;
	if (PyComplex_Check(op)) {
		return ((PyComplexObject *)op)->cval;
	} else {
		cv.real = PyFloat_AsDouble(op);
		cv.imag = 0.;
		return cv;
	}   
}

static void
complex_dealloc(op)
	PyObject *op;
{
	PyMem_DEL(op);
}


static void
complex_buf_repr(buf, v)
	char *buf;
	PyComplexObject *v;
{
	if (v->cval.real == 0.)
		sprintf(buf, "%.12gj", v->cval.imag);
	else
		sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
}

static int
complex_print(v, fp, flags)
	PyComplexObject *v;
	FILE *fp;
	int flags; /* Not used but required by interface */
{
	char buf[100];
	complex_buf_repr(buf, v);
	fputs(buf, fp);
	return 0;
}

static PyObject *
complex_repr(v)
	PyComplexObject *v;
{
	char buf[100];
	complex_buf_repr(buf, v);
	return PyString_FromString(buf);
}

static int
complex_compare(v, w)
	PyComplexObject *v, *w;
{
/* Note: "greater" and "smaller" have no meaning for complex numbers,
   but Python requires that they be defined nevertheless. */
	Py_complex i, j;
	i = v->cval;
	j = w->cval;
	if (i.real == j.real && i.imag == j.imag)
	   return 0;
	else if (i.real != j.real)
	   return (i.real < j.real) ? -1 : 1;
	else
	   return (i.imag < j.imag) ? -1 : 1;
}

static long
complex_hash(v)
	PyComplexObject *v;
{
	double intpart, fractpart;
	int expo;
	long hipart, x;
	/* This is designed so that Python numbers with the same
	   value hash to the same value, otherwise comparisons
	   of mapping keys will turn out weird */

#ifdef MPW /* MPW C modf expects pointer to extended as second argument */
{
	extended e;
	fractpart = modf(v->cval.real, &e);
	intpart = e;
}
#else
	fractpart = modf(v->cval.real, &intpart);
#endif

	if (fractpart == 0.0 && v->cval.imag == 0.0) {
		if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
			/* Convert to long int and use its hash... */
			PyObject *w = PyLong_FromDouble(v->cval.real);
			if (w == NULL)
				return -1;
			x = PyObject_Hash(w);
			Py_DECREF(w);
			return x;
		}
		x = (long)intpart;
	}
	else {
		fractpart = frexp(fractpart, &expo);
		fractpart = fractpart * 2147483648.0; /* 2**31 */
		hipart = (long)fractpart; /* Take the top 32 bits */
		fractpart = (fractpart - (double)hipart) * 2147483648.0;
						/* Get the next 32 bits */
		x = hipart + (long)fractpart + (long)intpart + (expo << 15);
						/* Combine everything */

		if (v->cval.imag != 0.0) { /* Hash the imaginary part */
			/* XXX Note that this hashes complex(x, y)
			   to the same value as complex(y, x).
			   Still better than it used to be :-) */
#ifdef MPW
			{
				extended e;
				fractpart = modf(v->cval.imag, &e);
				intpart = e;
			}
#else
			fractpart = modf(v->cval.imag, &intpart);
#endif
			fractpart = frexp(fractpart, &expo);
			fractpart = fractpart * 2147483648.0; /* 2**31 */
			hipart = (long)fractpart; /* Take the top 32 bits */
			fractpart =
				(fractpart - (double)hipart) * 2147483648.0;
						/* Get the next 32 bits */
			x ^= hipart + (long)fractpart +
				(long)intpart + (expo << 15);
						/* Combine everything */
		}
	}
	if (x == -1)
		x = -2;
	return x;
}

static PyObject *
complex_add(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
	Py_complex result;
	PyFPE_START_PROTECT("complex_add", return 0)
	result = c_sum(v->cval,w->cval);
	PyFPE_END_PROTECT(result)
	return PyComplex_FromCComplex(result);
}

static PyObject *
complex_sub(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
	Py_complex result;
	PyFPE_START_PROTECT("complex_sub", return 0)
	result = c_diff(v->cval,w->cval);
	PyFPE_END_PROTECT(result)
	return PyComplex_FromCComplex(result);
}

static PyObject *
complex_mul(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
	Py_complex result;
	PyFPE_START_PROTECT("complex_mul", return 0)
	result = c_prod(v->cval,w->cval);
	PyFPE_END_PROTECT(result)
	return PyComplex_FromCComplex(result);
}

static PyObject *
complex_div(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
	Py_complex quot;
	PyFPE_START_PROTECT("complex_div", return 0)
	errno = 0;
	quot = c_quot(v->cval,w->cval);
	PyFPE_END_PROTECT(quot)
	if (errno == EDOM) {
		PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
		return NULL;
	}
	return PyComplex_FromCComplex(quot);
}

static PyObject *
complex_remainder(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
        Py_complex div, mod;
	errno = 0;
	div = c_quot(v->cval,w->cval); /* The raw divisor value. */
	if (errno == EDOM) {
		PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
		return NULL;
	}
	div.real = floor(div.real); /* Use the floor of the real part. */
	div.imag = 0.0;
	mod = c_diff(v->cval, c_prod(w->cval, div));

	return PyComplex_FromCComplex(mod);
}


static PyObject *
complex_divmod(v, w)
	PyComplexObject *v;
	PyComplexObject *w;
{
        Py_complex div, mod;
	PyObject *d, *m, *z;
	errno = 0;
	div = c_quot(v->cval,w->cval); /* The raw divisor value. */
	if (errno == EDOM) {
		PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
		return NULL;
	}
	div.real = floor(div.real); /* Use the floor of the real part. */
	div.imag = 0.0;
	mod = c_diff(v->cval, c_prod(w->cval, div));
	d = PyComplex_FromCComplex(div);
	m = PyComplex_FromCComplex(mod);
	z = Py_BuildValue("(OO)", d, m);
	Py_XDECREF(d);
	Py_XDECREF(m);
	return z;
}

static PyObject *
complex_pow(v, w, z)
	PyComplexObject *v;
	PyObject *w;
	PyComplexObject *z;
{
	Py_complex p;
	Py_complex exponent;
	long int_exponent;

 	if ((PyObject *)z!=Py_None) {
		PyErr_SetString(PyExc_ValueError, "complex modulo");
		return NULL;
	}

	PyFPE_START_PROTECT("complex_pow", return 0)
	errno = 0;
	exponent = ((PyComplexObject*)w)->cval;
	int_exponent = (long)exponent.real;
	if (exponent.imag == 0. && exponent.real == int_exponent)
		p = c_powi(v->cval,int_exponent);
	else
		p = c_pow(v->cval,exponent);

	PyFPE_END_PROTECT(p)
	if (errno == ERANGE) {
		PyErr_SetString(PyExc_ValueError,
				"0.0 to a negative or complex power");
		return NULL;
	}

	return PyComplex_FromCComplex(p);
}

static PyObject *
complex_neg(v)
	PyComplexObject *v;
{
	Py_complex neg;
	neg.real = -v->cval.real;
	neg.imag = -v->cval.imag;
	return PyComplex_FromCComplex(neg);
}

static PyObject *
complex_pos(v)
	PyComplexObject *v;
{
	Py_INCREF(v);
	return (PyObject *)v;
}

static PyObject *
complex_abs(v)
	PyComplexObject *v;
{
	double result;
	PyFPE_START_PROTECT("complex_abs", return 0)
	result = hypot(v->cval.real,v->cval.imag);
	PyFPE_END_PROTECT(result)
	return PyFloat_FromDouble(result);
}

static int
complex_nonzero(v)
	PyComplexObject *v;
{
	return v->cval.real != 0.0 || v->cval.imag != 0.0;
}

static int
complex_coerce(pv, pw)
	PyObject **pv;
	PyObject **pw;
{
	Py_complex cval;
	cval.imag = 0.;
	if (PyInt_Check(*pw)) {
		cval.real = (double)PyInt_AsLong(*pw);
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	else if (PyLong_Check(*pw)) {
		cval.real = PyLong_AsDouble(*pw);
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	else if (PyFloat_Check(*pw)) {
		cval.real = PyFloat_AsDouble(*pw);
		*pw = PyComplex_FromCComplex(cval);
		Py_INCREF(*pv);
		return 0;
	}
	return 1; /* Can't do it */
}

static PyObject *
complex_int(v)
	PyObject *v;
{
	PyErr_SetString(PyExc_TypeError,
		   "can't convert complex to int; use e.g. int(abs(z))");
	return NULL;
}

static PyObject *
complex_long(v)
	PyObject *v;
{
	PyErr_SetString(PyExc_TypeError,
		   "can't convert complex to long; use e.g. long(abs(z))");
	return NULL;
}

static PyObject *
complex_float(v)
	PyObject *v;
{
	PyErr_SetString(PyExc_TypeError,
		   "can't convert complex to float; use e.g. abs(z)");
	return NULL;
}

static PyObject *
complex_conjugate(self, args)
	PyObject *self;
	PyObject *args;
{
	Py_complex c;
	if (!PyArg_ParseTuple(args, ""))
		return NULL;
	c = ((PyComplexObject *)self)->cval;
	c.imag = -c.imag;
	return PyComplex_FromCComplex(c);
}

static PyMethodDef complex_methods[] = {
	{"conjugate",	complex_conjugate,	1},
	{NULL,		NULL}		/* sentinel */
};


static PyObject *
complex_getattr(self, name)
	PyComplexObject *self;
	char *name;
{
	if (strcmp(name, "real") == 0)
		return (PyObject *)PyFloat_FromDouble(self->cval.real);
	else if (strcmp(name, "imag") == 0)
		return (PyObject *)PyFloat_FromDouble(self->cval.imag);
	else if (strcmp(name, "__members__") == 0)
		return Py_BuildValue("[ss]", "imag", "real");
	return Py_FindMethod(complex_methods, (PyObject *)self, name);
}

static PyNumberMethods complex_as_number = {
	(binaryfunc)complex_add, /*nb_add*/
	(binaryfunc)complex_sub, /*nb_subtract*/
	(binaryfunc)complex_mul, /*nb_multiply*/
	(binaryfunc)complex_div, /*nb_divide*/
	(binaryfunc)complex_remainder,	/*nb_remainder*/
	(binaryfunc)complex_divmod,	/*nb_divmod*/
	(ternaryfunc)complex_pow, /*nb_power*/
	(unaryfunc)complex_neg, /*nb_negative*/
	(unaryfunc)complex_pos, /*nb_positive*/
	(unaryfunc)complex_abs, /*nb_absolute*/
	(inquiry)complex_nonzero, /*nb_nonzero*/
	0,		/*nb_invert*/
	0,		/*nb_lshift*/
	0,		/*nb_rshift*/
	0,		/*nb_and*/
	0,		/*nb_xor*/
	0,		/*nb_or*/
	(coercion)complex_coerce, /*nb_coerce*/
	(unaryfunc)complex_int, /*nb_int*/
	(unaryfunc)complex_long, /*nb_long*/
	(unaryfunc)complex_float, /*nb_float*/
	0,		/*nb_oct*/
	0,		/*nb_hex*/
};

PyTypeObject PyComplex_Type = {
	PyObject_HEAD_INIT(&PyType_Type)
	0,
	"complex",
	sizeof(PyComplexObject),
	0,
	(destructor)complex_dealloc,	/*tp_dealloc*/
	(printfunc)complex_print,	/*tp_print*/
	(getattrfunc)complex_getattr,	/*tp_getattr*/
	0,				/*tp_setattr*/
	(cmpfunc)complex_compare,	/*tp_compare*/
	(reprfunc)complex_repr,		/*tp_repr*/
	&complex_as_number,    		/*tp_as_number*/
	0,				/*tp_as_sequence*/
	0,				/*tp_as_mapping*/
	(hashfunc)complex_hash, 	/*tp_hash*/
};

#endif