Objects/complexobject.c

   1
   2 /* Complex object implementation */
   3
   4 /* Borrows heavily from floatobject.c */
   5
   6 /* Submitted by Jim Hugunin */
   7
   8 #ifndef WITHOUT_COMPLEX
   9
  10 #include "Python.h"
  11
  12 #ifdef HAVE_LIMITS_H
  13 #include <limits.h>
  14 #endif
  15
  16
  17 /* elementary operations on complex numbers */
  18
  19 static Py_complex c_1 = {1., 0.};
  20
  21 Py_complex c_sum(Py_complex a, Py_complex b)
  22 {
  23         Py_complex r;
  24         r.real = a.real + b.real;
  25         r.imag = a.imag + b.imag;
  26         return r;
  27 }
  28
  29 Py_complex c_diff(Py_complex a, Py_complex b)
  30 {
  31         Py_complex r;
  32         r.real = a.real - b.real;
  33         r.imag = a.imag - b.imag;
  34         return r;
  35 }
  36
  37 Py_complex c_neg(Py_complex a)
  38 {
  39         Py_complex r;
  40         r.real = -a.real;
  41         r.imag = -a.imag;
  42         return r;
  43 }
  44
  45 Py_complex c_prod(Py_complex a, Py_complex b)
  46 {
  47         Py_complex r;
  48         r.real = a.real*b.real - a.imag*b.imag;
  49         r.imag = a.real*b.imag + a.imag*b.real;
  50         return r;
  51 }
  52
  53 Py_complex c_quot(Py_complex a, Py_complex b)
  54 {
  55         Py_complex r;
  56         double d = b.real*b.real + b.imag*b.imag;
  57         if (d == 0.)
  58                 errno = EDOM;
  59         r.real = (a.real*b.real + a.imag*b.imag)/d;
  60         r.imag = (a.imag*b.real - a.real*b.imag)/d;
  61         return r;
  62 }
  63
  64 Py_complex c_pow(Py_complex a, Py_complex b)
  65 {
  66         Py_complex r;
  67         double vabs,len,at,phase;
  68         if (b.real == 0. && b.imag == 0.) {
  69                 r.real = 1.;
  70                 r.imag = 0.;
  71         }
  72         else if (a.real == 0. && a.imag == 0.) {
  73                 if (b.imag != 0. || b.real < 0.)
  74                         errno = ERANGE;
  75                 r.real = 0.;
  76                 r.imag = 0.;
  77         }
  78         else {
  79                 vabs = hypot(a.real,a.imag);
  80                 len = pow(vabs,b.real);
  81                 at = atan2(a.imag, a.real);
  82                 phase = at*b.real;
  83                 if (b.imag != 0.0) {
  84                         len /= exp(at*b.imag);
  85                         phase += b.imag*log(vabs);
  86                 }
  87                 r.real = len*cos(phase);
  88                 r.imag = len*sin(phase);
  89         }
  90         return r;
  91 }
  92
  93 static Py_complex c_powu(Py_complex x, long n)
  94 {
  95         Py_complex r, p;
  96         long mask = 1;
  97         r = c_1;
  98         p = x;
  99         while (mask > 0 && n >= mask) {
 100                 if (n & mask)
 101                         r = c_prod(r,p);
 102                 mask <<= 1;
 103                 p = c_prod(p,p);
 104         }
 105         return r;
 106 }
 107
 108 static Py_complex c_powi(Py_complex x, long n)
 109 {
 110         Py_complex cn;
 111
 112         if (n > 100 || n < -100) {
 113                 cn.real = (double) n;
 114                 cn.imag = 0.;
 115                 return c_pow(x,cn);
 116         }
 117         else if (n > 0)
 118                 return c_powu(x,n);
 119         else
 120                 return c_quot(c_1,c_powu(x,-n));
 121
 122 }
 123
 124 PyObject *
 125 PyComplex_FromCComplex(Py_complex cval)
 126 {
 127         register PyComplexObject *op;
 128
 129         /* PyObject_New is inlined */
 130         op = (PyComplexObject *) PyObject_MALLOC(sizeof(PyComplexObject));
 131         if (op == NULL)
 132                 return PyErr_NoMemory();
 133         PyObject_INIT(op, &PyComplex_Type);
 134         op->cval = cval;
 135         return (PyObject *) op;
 136 }
 137
 138 PyObject *
 139 PyComplex_FromDoubles(double real, double imag)
 140 {
 141         Py_complex c;
 142         c.real = real;
 143         c.imag = imag;
 144         return PyComplex_FromCComplex(c);
 145 }
 146
 147 double
 148 PyComplex_RealAsDouble(PyObject *op)
 149 {
 150         if (PyComplex_Check(op)) {
 151                 return ((PyComplexObject *)op)->cval.real;
 152         }
 153         else {
 154                 return PyFloat_AsDouble(op);
 155         }
 156 }
 157
 158 double
 159 PyComplex_ImagAsDouble(PyObject *op)
 160 {
 161         if (PyComplex_Check(op)) {
 162                 return ((PyComplexObject *)op)->cval.imag;
 163         }
 164         else {
 165                 return 0.0;
 166         }
 167 }
 168
 169 Py_complex
 170 PyComplex_AsCComplex(PyObject *op)
 171 {
 172         Py_complex cv;
 173         if (PyComplex_Check(op)) {
 174                 return ((PyComplexObject *)op)->cval;
 175         }
 176         else {
 177                 cv.real = PyFloat_AsDouble(op);
 178                 cv.imag = 0.;
 179                 return cv;
 180         }
 181 }
 182
 183 static void
 184 complex_dealloc(PyObject *op)
 185 {
 186         PyObject_DEL(op);
 187 }
 188
 189
 190 static void
 191 complex_buf_repr(char *buf, PyComplexObject *v)
 192 {
 193         if (v->cval.real == 0.)
 194                 sprintf(buf, "%.12gj", v->cval.imag);
 195         else
 196                 sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
 197 }
 198
 199 static int
 200 complex_print(PyComplexObject *v, FILE *fp, int flags)
 201      /* flags -- not used but required by interface */
 202 {
 203         char buf[100];
 204         complex_buf_repr(buf, v);
 205         fputs(buf, fp);
 206         return 0;
 207 }
 208
 209 static PyObject *
 210 complex_repr(PyComplexObject *v)
 211 {
 212         char buf[100];
 213         complex_buf_repr(buf, v);
 214         return PyString_FromString(buf);
 215 }
 216
 217 static int
 218 complex_compare(PyComplexObject *v, PyComplexObject *w)
 219 {
 220         /* Note: "greater" and "smaller" have no meaning for complex numbers,
 221            but Python requires that they be defined nevertheless. */
 222         Py_complex i, j;
 223         i = v->cval;
 224         j = w->cval;
 225         if (i.real == j.real && i.imag == j.imag)
 226                 return 0;
 227         else if (i.real != j.real)
 228                 return (i.real < j.real) ? -1 : 1;
 229         else
 230                 return (i.imag < j.imag) ? -1 : 1;
 231 }
 232
 233 static long
 234 complex_hash(PyComplexObject *v)
 235 {
 236         long hashreal, hashimag, combined;
 237         hashreal = _Py_HashDouble(v->cval.real);
 238         if (hashreal == -1)
 239                 return -1;
 240         hashimag = _Py_HashDouble(v->cval.imag);
 241         if (hashimag == -1)
 242                 return -1;
 243         /* Note:  if the imaginary part is 0, hashimag is 0 now,
 244          * so the following returns hashreal unchanged.  This is
 245          * important because numbers of different types that
 246          * compare equal must have the same hash value, so that
 247          * hash(x + 0*j) must equal hash(x).
 248          */
 249         combined = hashreal + 1000003 * hashimag;
 250         if (combined == -1)
 251                 combined = -2;
 252         return combined;
 253 }
 254
 255 static PyObject *
 256 complex_add(PyComplexObject *v, PyComplexObject *w)
 257 {
 258         Py_complex result;
 259         PyFPE_START_PROTECT("complex_add", return 0)
 260         result = c_sum(v->cval,w->cval);
 261         PyFPE_END_PROTECT(result)
 262         return PyComplex_FromCComplex(result);
 263 }
 264
 265 static PyObject *
 266 complex_sub(PyComplexObject *v, PyComplexObject *w)
 267 {
 268         Py_complex result;
 269         PyFPE_START_PROTECT("complex_sub", return 0)
 270         result = c_diff(v->cval,w->cval);
 271         PyFPE_END_PROTECT(result)
 272         return PyComplex_FromCComplex(result);
 273 }
 274
 275 static PyObject *
 276 complex_mul(PyComplexObject *v, PyComplexObject *w)
 277 {
 278         Py_complex result;
 279         PyFPE_START_PROTECT("complex_mul", return 0)
 280         result = c_prod(v->cval,w->cval);
 281         PyFPE_END_PROTECT(result)
 282         return PyComplex_FromCComplex(result);
 283 }
 284
 285 static PyObject *
 286 complex_div(PyComplexObject *v, PyComplexObject *w)
 287 {
 288         Py_complex quot;
 289         PyFPE_START_PROTECT("complex_div", return 0)
 290         errno = 0;
 291         quot = c_quot(v->cval,w->cval);
 292         PyFPE_END_PROTECT(quot)
 293         if (errno == EDOM) {
 294                 PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
 295                 return NULL;
 296         }
 297         return PyComplex_FromCComplex(quot);
 298 }
 299
 300 static PyObject *
 301 complex_remainder(PyComplexObject *v, PyComplexObject *w)
 302 {
 303         Py_complex div, mod;
 304         errno = 0;
 305         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 306         if (errno == EDOM) {
 307                 PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
 308                 return NULL;
 309         }
 310         div.real = floor(div.real); /* Use the floor of the real part. */
 311         div.imag = 0.0;
 312         mod = c_diff(v->cval, c_prod(w->cval, div));
 313
 314         return PyComplex_FromCComplex(mod);
 315 }
 316
 317
 318 static PyObject *
 319 complex_divmod(PyComplexObject *v, PyComplexObject *w)
 320 {
 321         Py_complex div, mod;
 322         PyObject *d, *m, *z;
 323         errno = 0;
 324         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 325         if (errno == EDOM) {
 326                 PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
 327                 return NULL;
 328         }
 329         div.real = floor(div.real); /* Use the floor of the real part. */
 330         div.imag = 0.0;
 331         mod = c_diff(v->cval, c_prod(w->cval, div));
 332         d = PyComplex_FromCComplex(div);
 333         m = PyComplex_FromCComplex(mod);
 334         z = Py_BuildValue("(OO)", d, m);
 335         Py_XDECREF(d);
 336         Py_XDECREF(m);
 337         return z;
 338 }
 339
 340 static PyObject *
 341 complex_pow(PyComplexObject *v, PyObject *w, PyComplexObject *z)
 342 {
 343         Py_complex p;
 344         Py_complex exponent;
 345         long int_exponent;
 346
 347         if ((PyObject *)z!=Py_None) {
 348                 PyErr_SetString(PyExc_ValueError, "complex modulo");
 349                 return NULL;
 350         }
 351         PyFPE_START_PROTECT("complex_pow", return 0)
 352         errno = 0;
 353         exponent = ((PyComplexObject*)w)->cval;
 354         int_exponent = (long)exponent.real;
 355         if (exponent.imag == 0. && exponent.real == int_exponent)
 356                 p = c_powi(v->cval,int_exponent);
 357         else
 358                 p = c_pow(v->cval,exponent);
 359
 360         PyFPE_END_PROTECT(p)
 361         if (errno == ERANGE) {
 362                 PyErr_SetString(PyExc_ValueError,
 363                                 "0.0 to a negative or complex power");
 364                 return NULL;
 365         }
 366         return PyComplex_FromCComplex(p);
 367 }
 368
 369 static PyObject *
 370 complex_neg(PyComplexObject *v)
 371 {
 372         Py_complex neg;
 373         neg.real = -v->cval.real;
 374         neg.imag = -v->cval.imag;
 375         return PyComplex_FromCComplex(neg);
 376 }
 377
 378 static PyObject *
 379 complex_pos(PyComplexObject *v)
 380 {
 381         Py_INCREF(v);
 382         return (PyObject *)v;
 383 }
 384
 385 static PyObject *
 386 complex_abs(PyComplexObject *v)
 387 {
 388         double result;
 389         PyFPE_START_PROTECT("complex_abs", return 0)
 390         result = hypot(v->cval.real,v->cval.imag);
 391         PyFPE_END_PROTECT(result)
 392         return PyFloat_FromDouble(result);
 393 }
 394
 395 static int
 396 complex_nonzero(PyComplexObject *v)
 397 {
 398         return v->cval.real != 0.0 || v->cval.imag != 0.0;
 399 }
 400
 401 static int
 402 complex_coerce(PyObject **pv, PyObject **pw)
 403 {
 404         Py_complex cval;
 405         cval.imag = 0.;
 406         if (PyInt_Check(*pw)) {
 407                 cval.real = (double)PyInt_AsLong(*pw);
 408                 *pw = PyComplex_FromCComplex(cval);
 409                 Py_INCREF(*pv);
 410                 return 0;
 411         }
 412         else if (PyLong_Check(*pw)) {
 413                 cval.real = PyLong_AsDouble(*pw);
 414                 *pw = PyComplex_FromCComplex(cval);
 415                 Py_INCREF(*pv);
 416                 return 0;
 417         }
 418         else if (PyFloat_Check(*pw)) {
 419                 cval.real = PyFloat_AsDouble(*pw);
 420                 *pw = PyComplex_FromCComplex(cval);
 421                 Py_INCREF(*pv);
 422                 return 0;
 423         }
 424         return 1; /* Can't do it */
 425 }
 426
 427 static PyObject *
 428 complex_int(PyObject *v)
 429 {
 430         PyErr_SetString(PyExc_TypeError,
 431                    "can't convert complex to int; use e.g. int(abs(z))");
 432         return NULL;
 433 }
 434
 435 static PyObject *
 436 complex_long(PyObject *v)
 437 {
 438         PyErr_SetString(PyExc_TypeError,
 439                    "can't convert complex to long; use e.g. long(abs(z))");
 440         return NULL;
 441 }
 442
 443 static PyObject *
 444 complex_float(PyObject *v)
 445 {
 446         PyErr_SetString(PyExc_TypeError,
 447                    "can't convert complex to float; use e.g. abs(z)");
 448         return NULL;
 449 }
 450
 451 static PyObject *
 452 complex_conjugate(PyObject *self, PyObject *args)
 453 {
 454         Py_complex c;
 455         if (!PyArg_ParseTuple(args, ":conjugate"))
 456                 return NULL;
 457         c = ((PyComplexObject *)self)->cval;
 458         c.imag = -c.imag;
 459         return PyComplex_FromCComplex(c);
 460 }
 461
 462 static PyMethodDef complex_methods[] = {
 463         {"conjugate",   complex_conjugate,      1},
 464         {NULL,          NULL}           /* sentinel */
 465 };
 466
 467
 468 static PyObject *
 469 complex_getattr(PyComplexObject *self, char *name)
 470 {
 471         if (strcmp(name, "real") == 0)
 472                 return (PyObject *)PyFloat_FromDouble(self->cval.real);
 473         else if (strcmp(name, "imag") == 0)
 474                 return (PyObject *)PyFloat_FromDouble(self->cval.imag);
 475         else if (strcmp(name, "__members__") == 0)
 476                 return Py_BuildValue("[ss]", "imag", "real");
 477         return Py_FindMethod(complex_methods, (PyObject *)self, name);
 478 }
 479
 480 static PyNumberMethods complex_as_number = {
 481         (binaryfunc)complex_add, /*nb_add*/
 482         (binaryfunc)complex_sub, /*nb_subtract*/
 483         (binaryfunc)complex_mul, /*nb_multiply*/
 484         (binaryfunc)complex_div, /*nb_divide*/
 485         (binaryfunc)complex_remainder,  /*nb_remainder*/
 486         (binaryfunc)complex_divmod,     /*nb_divmod*/
 487         (ternaryfunc)complex_pow, /*nb_power*/
 488         (unaryfunc)complex_neg, /*nb_negative*/
 489         (unaryfunc)complex_pos, /*nb_positive*/
 490         (unaryfunc)complex_abs, /*nb_absolute*/
 491         (inquiry)complex_nonzero, /*nb_nonzero*/
 492         0,              /*nb_invert*/
 493         0,              /*nb_lshift*/
 494         0,              /*nb_rshift*/
 495         0,              /*nb_and*/
 496         0,              /*nb_xor*/
 497         0,              /*nb_or*/
 498         (coercion)complex_coerce, /*nb_coerce*/
 499         (unaryfunc)complex_int, /*nb_int*/
 500         (unaryfunc)complex_long, /*nb_long*/
 501         (unaryfunc)complex_float, /*nb_float*/
 502         0,              /*nb_oct*/
 503         0,              /*nb_hex*/
 504 };
 505
 506 PyTypeObject PyComplex_Type = {
 507         PyObject_HEAD_INIT(&PyType_Type)
 508         0,
 509         "complex",
 510         sizeof(PyComplexObject),
 511         0,
 512         (destructor)complex_dealloc,    /*tp_dealloc*/
 513         (printfunc)complex_print,       /*tp_print*/
 514         (getattrfunc)complex_getattr,   /*tp_getattr*/
 515         0,                              /*tp_setattr*/
 516         (cmpfunc)complex_compare,       /*tp_compare*/
 517         (reprfunc)complex_repr,         /*tp_repr*/
 518         &complex_as_number,             /*tp_as_number*/
 519         0,                              /*tp_as_sequence*/
 520         0,                              /*tp_as_mapping*/
 521         (hashfunc)complex_hash,         /*tp_hash*/
 522 };
 523
 524 #endif