653 lines
14 KiB
C
653 lines
14 KiB
C
/***********************************************************
|
|
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
|