complex numbers a la Konrad Hinsen
diff --git a/Include/complexobject.h b/Include/complexobject.h
new file mode 100644
index 0000000..09a8683
--- /dev/null
+++ b/Include/complexobject.h
@@ -0,0 +1,49 @@
+#ifndef COMPLEXOBJECT_H
+#define COMPLEXOBJECT_H
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Complex number structure */
+
+typedef struct {
+ double real;
+ double imag;
+} complex;
+
+/* Operations on complex numbers from complexmodule.c */
+
+extern complex c_sum();
+extern complex c_diff();
+extern complex c_neg();
+extern complex c_prod();
+extern complex c_quot();
+extern complex c_pow();
+
+
+/* Complex object interface */
+
+/*
+PyComplexObject represents a complex number with double-precision
+real and imaginary parts.
+*/
+
+typedef struct {
+ PyObject_HEAD
+ complex cval;
+} PyComplexObject;
+
+extern DL_IMPORT(PyTypeObject) PyComplex_Type;
+
+#define PyComplex_Check(op) ((op)->ob_type == &PyComplex_Type)
+
+extern PyObject *PyComplex_FromCComplex Py_PROTO((complex));
+extern PyObject *PyComplex_FromDoubles Py_PROTO((double real, double imag));
+
+extern double PyComplex_RealAsDouble Py_PROTO((PyObject *op));
+extern double PyComplex_ImagAsDouble Py_PROTO((PyObject *op));
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* !COMPLEXOBJECT_H */
diff --git a/Objects/complexobject.c b/Objects/complexobject.c
new file mode 100644
index 0000000..8299b0b
--- /dev/null
+++ b/Objects/complexobject.c
@@ -0,0 +1,562 @@
+/* Complex object implementation */
+
+/* Borrows heavily from floatobject.c */
+
+#ifndef WITHOUT_COMPLEX
+
+#include "allobjects.h"
+#include "modsupport.h"
+
+#include <errno.h>
+#include "mymath.h"
+
+#ifdef i860
+/* Cray APP has bogus definition of HUGE_VAL in <math.h> */
+#undef HUGE_VAL
+#endif
+
+#ifdef HUGE_VAL
+#define CHECK(x) if (errno != 0) ; \
+ else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
+ else errno = ERANGE
+#else
+#define CHECK(x) /* Don't know how to check */
+#endif
+
+#ifdef HAVE_LIMITS_H
+#include <limits.h>
+#endif
+
+#ifndef LONG_MAX
+#define LONG_MAX 0X7FFFFFFFL
+#endif
+
+#ifndef LONG_MIN
+#define LONG_MIN (-LONG_MAX-1)
+#endif
+
+#ifdef __NeXT__
+#ifdef __sparc__
+/*
+ * This works around a bug in the NS/Sparc 3.3 pre-release
+ * limits.h header file.
+ * 10-Feb-1995 bwarsaw@cnri.reston.va.us
+ */
+#undef LONG_MIN
+#define LONG_MIN (-LONG_MAX-1)
+#endif
+#endif
+
+#if !defined(__STDC__) && !defined(macintosh)
+extern double fmod PROTO((double, double));
+extern double pow PROTO((double, double));
+#endif
+
+
+/* elementary operations on complex numbers */
+
+int c_error;
+static complex c_1 = {1., 0.};
+
+complex c_sum(a,b)
+ complex a,b;
+{
+ complex r;
+ r.real = a.real + b.real;
+ r.imag = a.imag + b.imag;
+ return r;
+}
+
+complex c_diff(a,b)
+ complex a,b;
+{
+ complex r;
+ r.real = a.real - b.real;
+ r.imag = a.imag - b.imag;
+ return r;
+}
+
+complex c_neg(a)
+ complex a;
+{
+ complex r;
+ r.real = -a.real;
+ r.imag = -a.imag;
+ return r;
+}
+
+complex c_prod(a,b)
+ complex a,b;
+{
+ complex r;
+ r.real = a.real*b.real - a.imag*b.imag;
+ r.imag = a.real*b.imag + a.imag*b.real;
+ return r;
+}
+
+complex c_quot(a,b)
+ complex a,b;
+{
+ complex r;
+ double d = b.real*b.real + b.imag*b.imag;
+ if (d == 0.)
+ c_error = 1;
+ r.real = (a.real*b.real + a.imag*b.imag)/d;
+ r.imag = (a.imag*b.real - a.real*b.imag)/d;
+ return r;
+}
+
+complex c_pow(a,b)
+ complex a,b;
+{
+ 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.)
+ c_error = 2;
+ 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;
+}
+
+complex c_powu(x, n)
+ complex x;
+ long n;
+{
+ complex r = c_1;
+ complex p = x;
+ long mask = 1;
+ while (mask > 0 && n >= mask) {
+ if (n & mask)
+ r = c_prod(r,p);
+ mask <<= 1;
+ p = c_prod(p,p);
+ }
+ return r;
+}
+
+complex c_powi(x, n)
+ complex x;
+ long n;
+{
+ 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(complex cval)
+{
+ register complexobject *op = (complexobject *) malloc(sizeof(complexobject));
+ if (op == NULL)
+ return err_nomem();
+ op->ob_type = &Complextype;
+ op->cval = cval;
+ NEWREF(op);
+ return (object *) op;
+}
+
+PyObject *
+PyComplex_FromDoubles(double real, double imag) {
+ 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 void
+complex_dealloc(op)
+ object *op;
+{
+ DEL(op);
+}
+
+
+void
+complex_buf_repr(buf, v)
+ char *buf;
+ complexobject *v;
+{
+ if (v->cval.real == 0.)
+ sprintf(buf, "%.12gi", v->cval.imag);
+ else
+ sprintf(buf, "(%.12g%+.12gi)", v->cval.real, v->cval.imag);
+}
+
+static int
+complex_print(v, fp, flags)
+ complexobject *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 object *
+complex_repr(v)
+ complexobject *v;
+{
+ char buf[100];
+ complex_buf_repr(buf, v);
+ return newstringobject(buf);
+}
+
+static int
+complex_compare(v, w)
+ complexobject *v, *w;
+{
+/* Note: "greater" and "smaller" have no meaning for complex numbers,
+ but Python requires that they be defined nevertheless. */
+ complex i = v->cval;
+ complex 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)
+ complexobject *v;
+{
+ double intpart, fractpart;
+ int expo;
+ long 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) {
+ if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
+ /* Convert to long int and use its hash... */
+ object *w = dnewlongobject(v->cval.real);
+ if (w == NULL)
+ return -1;
+ x = hashobject(w);
+ DECREF(w);
+ return x;
+ }
+ x = (long)intpart;
+ }
+ else {
+ fractpart = frexp(fractpart, &expo);
+ fractpart = fractpart*2147483648.0; /* 2**31 */
+ x = (long) (intpart + fractpart) ^ expo; /* Rather arbitrary */
+ }
+ if (x == -1)
+ x = -2;
+ return x;
+}
+
+static object *
+complex_add(v, w)
+ complexobject *v;
+ complexobject *w;
+{
+ return newcomplexobject(c_sum(v->cval,w->cval));
+}
+
+static object *
+complex_sub(v, w)
+ complexobject *v;
+ complexobject *w;
+{
+ return newcomplexobject(c_diff(v->cval,w->cval));
+}
+
+static object *
+complex_mul(v, w)
+ complexobject *v;
+ complexobject *w;
+{
+ return newcomplexobject(c_prod(v->cval,w->cval));
+}
+
+static object *
+complex_div(v, w)
+ complexobject *v;
+ complexobject *w;
+{
+ complex quot;
+ c_error = 0;
+ quot = c_quot(v->cval,w->cval);
+ if (c_error == 1) {
+ err_setstr(ZeroDivisionError, "float division");
+ return NULL;
+ }
+ return newcomplexobject(quot);
+}
+
+
+static object *
+complex_pow(v, w, z)
+ complexobject *v;
+ object *w;
+ complexobject *z;
+{
+ complex p;
+ complex exponent;
+ long int_exponent;
+
+ if ((object *)z!=None) {
+ err_setstr(ValueError, "complex modulo");
+ return NULL;
+ }
+
+ c_error = 0;
+ exponent = ((complexobject*)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);
+
+ if (c_error == 2) {
+ err_setstr(ValueError, "0.0 to a negative or complex power");
+ return NULL;
+ }
+
+ return newcomplexobject(p);
+}
+
+static object *
+complex_neg(v)
+ complexobject *v;
+{
+ complex neg;
+ neg.real = -v->cval.real;
+ neg.imag = -v->cval.imag;
+ return newcomplexobject(neg);
+}
+
+static object *
+complex_pos(v)
+ complexobject *v;
+{
+ INCREF(v);
+ return (object *)v;
+}
+
+static object *
+complex_abs(v)
+ complexobject *v;
+{
+ return newfloatobject(hypot(v->cval.real,v->cval.imag));
+}
+
+static int
+complex_nonzero(v)
+ complexobject *v;
+{
+ return v->cval.real != 0.0 && v->cval.imag != 0.0;
+}
+
+static int
+complex_coerce(pv, pw)
+ object **pv;
+ object **pw;
+{
+ complex cval;
+ cval.imag = 0.;
+ if (is_intobject(*pw)) {
+ cval.real = (double)getintvalue(*pw);
+ *pw = newcomplexobject(cval);
+ INCREF(*pv);
+ return 0;
+ }
+ else if (is_longobject(*pw)) {
+ cval.real = dgetlongvalue(*pw);
+ *pw = newcomplexobject(cval);
+ INCREF(*pv);
+ return 0;
+ }
+ else if (is_floatobject(*pw)) {
+ cval.real = getfloatvalue(*pw);
+ *pw = newcomplexobject(cval);
+ INCREF(*pv);
+ return 0;
+ }
+ return 1; /* Can't do it */
+}
+
+static object *
+complex_int(v)
+ object *v;
+{
+ double x = ((complexobject *)v)->cval.real;
+ if (x < 0 ? (x = ceil(x)) < (double)LONG_MIN
+ : (x = floor(x)) > (double)LONG_MAX) {
+ err_setstr(OverflowError, "float too large to convert");
+ return NULL;
+ }
+ return newintobject((long)x);
+}
+
+static object *
+complex_long(v)
+ object *v;
+{
+ double x = ((complexobject *)v)->cval.real;
+ return dnewlongobject(x);
+}
+
+static object *
+complex_float(v)
+ object *v;
+{
+ double x = ((complexobject *)v)->cval.real;
+ return newfloatobject(x);
+}
+
+
+static object *
+complex_new(self, args)
+ object *self;
+ object *args;
+{
+ int n;
+ complex cval;
+
+ cval.imag = 0.;
+ if (!PyArg_ParseTuple(args, "d|d", &cval.real, &cval.imag))
+ return NULL;
+ return newcomplexobject(cval);
+}
+
+static object *
+complex_conjugate(self)
+ object *self;
+{
+ complex c = ((complexobject *)self)->cval;
+ c.imag = -c.imag;
+ return newcomplexobject(c);
+}
+
+static PyMethodDef complex_methods[] = {
+ {"conjugate", (PyCFunction)complex_conjugate, 1},
+ {NULL, NULL} /* sentinel */
+};
+
+
+static object *
+complex_getattr(self, name)
+ complexobject *self;
+ char *name;
+{
+ complex cval;
+ if (strcmp(name, "real") == 0)
+ return (object *)newfloatobject(self->cval.real);
+ else if (strcmp(name, "imag") == 0)
+ return (object *)newfloatobject(self->cval.imag);
+ else if (strcmp(name, "conj") == 0) {
+ cval.real = self->cval.real;
+ cval.imag = -self->cval.imag;
+ return (object *)newcomplexobject(cval);
+ }
+ return findmethod(complex_methods, (object *)self, name);
+}
+
+static number_methods complex_as_number = {
+ (binaryfunc)complex_add, /*nb_add*/
+ (binaryfunc)complex_sub, /*nb_subtract*/
+ (binaryfunc)complex_mul, /*nb_multiply*/
+ (binaryfunc)complex_div, /*nb_divide*/
+ 0, /*nb_remainder*/
+ 0, /*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*/
+};
+
+typeobject Complextype = {
+ OB_HEAD_INIT(&Typetype)
+ 0,
+ "complex",
+ sizeof(complexobject),
+ 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