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