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