Objects/complexobject.c

   1 /***********************************************************
   2 Copyright 1991-1995 by Stichting Mathematisch Centrum, Amsterdam,
   3 The Netherlands.
   4
   5                         All Rights Reserved
   6
   7 Permission to use, copy, modify, and distribute this software and its
   8 documentation for any purpose and without fee is hereby granted,
   9 provided that the above copyright notice appear in all copies and that
  10 both that copyright notice and this permission notice appear in
  11 supporting documentation, and that the names of Stichting Mathematisch
  12 Centrum or CWI or Corporation for National Research Initiatives or
  13 CNRI not be used in advertising or publicity pertaining to
  14 distribution of the software without specific, written prior
  15 permission.
  16
  17 While CWI is the initial source for this software, a modified version
  18 is made available by the Corporation for National Research Initiatives
  19 (CNRI) at the Internet address ftp://ftp.python.org.
  20
  21 STICHTING MATHEMATISCH CENTRUM AND CNRI DISCLAIM ALL WARRANTIES WITH
  22 REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
  23 MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH
  24 CENTRUM OR CNRI BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
  25 DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
  26 PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
  27 TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
  28 PERFORMANCE OF THIS SOFTWARE.
  29
  30 ******************************************************************/
  31
  32 /* Complex object implementation */
  33
  34 /* Borrows heavily from floatobject.c */
  35
  36 /* Submitted by Jim Hugunin */
  37
  38 #ifndef WITHOUT_COMPLEX
  39
  40 #include "Python.h"
  41 #include "mymath.h"
  42
  43 #ifdef HAVE_LIMITS_H
  44 #include <limits.h>
  45 #endif
  46
  47
  48 /* elementary operations on complex numbers */
  49
  50 static Py_complex c_1 = {1., 0.};
  51
  52 Py_complex c_sum(a,b)
  53         Py_complex a,b;
  54 {
  55         Py_complex r;
  56         r.real = a.real + b.real;
  57         r.imag = a.imag + b.imag;
  58         return r;
  59 }
  60
  61 Py_complex c_diff(a,b)
  62         Py_complex a,b;
  63 {
  64         Py_complex r;
  65         r.real = a.real - b.real;
  66         r.imag = a.imag - b.imag;
  67         return r;
  68 }
  69
  70 Py_complex c_neg(a)
  71         Py_complex a;
  72 {
  73         Py_complex r;
  74         r.real = -a.real;
  75         r.imag = -a.imag;
  76         return r;
  77 }
  78
  79 Py_complex c_prod(a,b)
  80         Py_complex a,b;
  81 {
  82         Py_complex r;
  83         r.real = a.real*b.real - a.imag*b.imag;
  84         r.imag = a.real*b.imag + a.imag*b.real;
  85         return r;
  86 }
  87
  88 Py_complex c_quot(a,b)
  89         Py_complex a,b;
  90 {
  91         Py_complex r;
  92         double d = b.real*b.real + b.imag*b.imag;
  93         if (d == 0.)
  94                 errno = EDOM;
  95         r.real = (a.real*b.real + a.imag*b.imag)/d;
  96         r.imag = (a.imag*b.real - a.real*b.imag)/d;
  97         return r;
  98 }
  99
 100 Py_complex c_pow(a,b)
 101         Py_complex a,b;
 102 {
 103         Py_complex r;
 104         double vabs,len,at,phase;
 105         if (b.real == 0. && b.imag == 0.) {
 106                 r.real = 1.;
 107                 r.imag = 0.;
 108         }
 109         else if (a.real == 0. && a.imag == 0.) {
 110                 if (b.imag != 0. || b.real < 0.)
 111                         errno = ERANGE;
 112                 r.real = 0.;
 113                 r.imag = 0.;
 114         }
 115         else {
 116                 vabs = hypot(a.real,a.imag);
 117                 len = pow(vabs,b.real);
 118                 at = atan2(a.imag, a.real);
 119                 phase = at*b.real;
 120                 if (b.imag != 0.0) {
 121                         len /= exp(at*b.imag);
 122                         phase += b.imag*log(vabs);
 123                 }
 124                 r.real = len*cos(phase);
 125                 r.imag = len*sin(phase);
 126         }
 127         return r;
 128 }
 129
 130 static Py_complex c_powu(x, n)
 131         Py_complex x;
 132         long n;
 133 {
 134         Py_complex r, p;
 135         long mask = 1;
 136         r = c_1;
 137         p = x;
 138         while (mask > 0 && n >= mask) {
 139                 if (n & mask)
 140                         r = c_prod(r,p);
 141                 mask <<= 1;
 142                 p = c_prod(p,p);
 143         }
 144         return r;
 145 }
 146
 147 static Py_complex c_powi(x, n)
 148         Py_complex x;
 149         long n;
 150 {
 151         Py_complex cn;
 152
 153         if (n > 100 || n < -100) {
 154                 cn.real = (double) n;
 155                 cn.imag = 0.;
 156                 return c_pow(x,cn);
 157         }
 158         else if (n > 0)
 159                 return c_powu(x,n);
 160         else
 161                 return c_quot(c_1,c_powu(x,-n));
 162
 163 }
 164
 165 PyObject *
 166 PyComplex_FromCComplex(cval)
 167         Py_complex cval;
 168 {
 169         register PyComplexObject *op =
 170                 (PyComplexObject *) malloc(sizeof(PyComplexObject));
 171         if (op == NULL)
 172                 return PyErr_NoMemory();
 173         op->ob_type = &PyComplex_Type;
 174         op->cval = cval;
 175         _Py_NewReference((PyObject *)op);
 176         return (PyObject *) op;
 177 }
 178
 179 PyObject *
 180 PyComplex_FromDoubles(real, imag)
 181         double real, imag;
 182 {
 183         Py_complex c;
 184         c.real = real;
 185         c.imag = imag;
 186         return PyComplex_FromCComplex(c);
 187 }
 188
 189 double
 190 PyComplex_RealAsDouble(op)
 191         PyObject *op;
 192 {
 193         if (PyComplex_Check(op)) {
 194                 return ((PyComplexObject *)op)->cval.real;
 195         } else {
 196                 return PyFloat_AsDouble(op);
 197         }
 198 }
 199
 200 double
 201 PyComplex_ImagAsDouble(op)
 202         PyObject *op;
 203 {
 204         if (PyComplex_Check(op)) {
 205                 return ((PyComplexObject *)op)->cval.imag;
 206         } else {
 207                 return 0.0;
 208         }
 209 }
 210
 211 Py_complex
 212 PyComplex_AsCComplex(op)
 213         PyObject *op;
 214 {
 215         Py_complex cv;
 216         if (PyComplex_Check(op)) {
 217                 return ((PyComplexObject *)op)->cval;
 218         } else {
 219                 cv.real = PyFloat_AsDouble(op);
 220                 cv.imag = 0.;
 221                 return cv;
 222         }
 223 }
 224
 225 static void
 226 complex_dealloc(op)
 227         PyObject *op;
 228 {
 229         PyMem_DEL(op);
 230 }
 231
 232
 233 static void
 234 complex_buf_repr(buf, v)
 235         char *buf;
 236         PyComplexObject *v;
 237 {
 238         if (v->cval.real == 0.)
 239                 sprintf(buf, "%.12gj", v->cval.imag);
 240         else
 241                 sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
 242 }
 243
 244 static int
 245 complex_print(v, fp, flags)
 246         PyComplexObject *v;
 247         FILE *fp;
 248         int flags; /* Not used but required by interface */
 249 {
 250         char buf[100];
 251         complex_buf_repr(buf, v);
 252         fputs(buf, fp);
 253         return 0;
 254 }
 255
 256 static PyObject *
 257 complex_repr(v)
 258         PyComplexObject *v;
 259 {
 260         char buf[100];
 261         complex_buf_repr(buf, v);
 262         return PyString_FromString(buf);
 263 }
 264
 265 static int
 266 complex_compare(v, w)
 267         PyComplexObject *v, *w;
 268 {
 269 /* Note: "greater" and "smaller" have no meaning for complex numbers,
 270    but Python requires that they be defined nevertheless. */
 271         Py_complex i, j;
 272         i = v->cval;
 273         j = w->cval;
 274         if (i.real == j.real && i.imag == j.imag)
 275            return 0;
 276         else if (i.real != j.real)
 277            return (i.real < j.real) ? -1 : 1;
 278         else
 279            return (i.imag < j.imag) ? -1 : 1;
 280 }
 281
 282 static long
 283 complex_hash(v)
 284         PyComplexObject *v;
 285 {
 286         double intpart, fractpart;
 287         int expo;
 288         long hipart, x;
 289         /* This is designed so that Python numbers with the same
 290            value hash to the same value, otherwise comparisons
 291            of mapping keys will turn out weird */
 292
 293 #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
 294 {
 295         extended e;
 296         fractpart = modf(v->cval.real, &e);
 297         intpart = e;
 298 }
 299 #else
 300         fractpart = modf(v->cval.real, &intpart);
 301 #endif
 302
 303         if (fractpart == 0.0 && v->cval.imag == 0.0) {
 304                 if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
 305                         /* Convert to long int and use its hash... */
 306                         PyObject *w = PyLong_FromDouble(v->cval.real);
 307                         if (w == NULL)
 308                                 return -1;
 309                         x = PyObject_Hash(w);
 310                         Py_DECREF(w);
 311                         return x;
 312                 }
 313                 x = (long)intpart;
 314         }
 315         else {
 316                 fractpart = frexp(fractpart, &expo);
 317                 fractpart = fractpart * 2147483648.0; /* 2**31 */
 318                 hipart = (long)fractpart; /* Take the top 32 bits */
 319                 fractpart = (fractpart - (double)hipart) * 2147483648.0;
 320                                                 /* Get the next 32 bits */
 321                 x = hipart + (long)fractpart + (long)intpart + (expo << 15);
 322                                                 /* Combine everything */
 323
 324                 if (v->cval.imag != 0.0) { /* Hash the imaginary part */
 325                         /* XXX Note that this hashes complex(x, y)
 326                            to the same value as complex(y, x).
 327                            Still better than it used to be :-) */
 328 #ifdef MPW
 329                         {
 330                                 extended e;
 331                                 fractpart = modf(v->cval.imag, &e);
 332                                 intpart = e;
 333                         }
 334 #else
 335                         fractpart = modf(v->cval.imag, &intpart);
 336 #endif
 337                         fractpart = frexp(fractpart, &expo);
 338                         fractpart = fractpart * 2147483648.0; /* 2**31 */
 339                         hipart = (long)fractpart; /* Take the top 32 bits */
 340                         fractpart =
 341                                 (fractpart - (double)hipart) * 2147483648.0;
 342                                                 /* Get the next 32 bits */
 343                         x ^= hipart + (long)fractpart +
 344                                 (long)intpart + (expo << 15);
 345                                                 /* Combine everything */
 346                 }
 347         }
 348         if (x == -1)
 349                 x = -2;
 350         return x;
 351 }
 352
 353 static PyObject *
 354 complex_add(v, w)
 355         PyComplexObject *v;
 356         PyComplexObject *w;
 357 {
 358         Py_complex result;
 359         PyFPE_START_PROTECT("complex_add", return 0)
 360         result = c_sum(v->cval,w->cval);
 361         PyFPE_END_PROTECT(result)
 362         return PyComplex_FromCComplex(result);
 363 }
 364
 365 static PyObject *
 366 complex_sub(v, w)
 367         PyComplexObject *v;
 368         PyComplexObject *w;
 369 {
 370         Py_complex result;
 371         PyFPE_START_PROTECT("complex_sub", return 0)
 372         result = c_diff(v->cval,w->cval);
 373         PyFPE_END_PROTECT(result)
 374         return PyComplex_FromCComplex(result);
 375 }
 376
 377 static PyObject *
 378 complex_mul(v, w)
 379         PyComplexObject *v;
 380         PyComplexObject *w;
 381 {
 382         Py_complex result;
 383         PyFPE_START_PROTECT("complex_mul", return 0)
 384         result = c_prod(v->cval,w->cval);
 385         PyFPE_END_PROTECT(result)
 386         return PyComplex_FromCComplex(result);
 387 }
 388
 389 static PyObject *
 390 complex_div(v, w)
 391         PyComplexObject *v;
 392         PyComplexObject *w;
 393 {
 394         Py_complex quot;
 395         PyFPE_START_PROTECT("complex_div", return 0)
 396         errno = 0;
 397         quot = c_quot(v->cval,w->cval);
 398         PyFPE_END_PROTECT(quot)
 399         if (errno == EDOM) {
 400                 PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
 401                 return NULL;
 402         }
 403         return PyComplex_FromCComplex(quot);
 404 }
 405
 406 static PyObject *
 407 complex_remainder(v, w)
 408         PyComplexObject *v;
 409         PyComplexObject *w;
 410 {
 411         Py_complex div, mod;
 412         errno = 0;
 413         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 414         if (errno == EDOM) {
 415                 PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
 416                 return NULL;
 417         }
 418         div.real = floor(div.real); /* Use the floor of the real part. */
 419         div.imag = 0.0;
 420         mod = c_diff(v->cval, c_prod(w->cval, div));
 421
 422         return PyComplex_FromCComplex(mod);
 423 }
 424
 425
 426 static PyObject *
 427 complex_divmod(v, w)
 428         PyComplexObject *v;
 429         PyComplexObject *w;
 430 {
 431         Py_complex div, mod;
 432         PyObject *d, *m, *z;
 433         errno = 0;
 434         div = c_quot(v->cval,w->cval); /* The raw divisor value. */
 435         if (errno == EDOM) {
 436                 PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
 437                 return NULL;
 438         }
 439         div.real = floor(div.real); /* Use the floor of the real part. */
 440         div.imag = 0.0;
 441         mod = c_diff(v->cval, c_prod(w->cval, div));
 442         d = PyComplex_FromCComplex(div);
 443         m = PyComplex_FromCComplex(mod);
 444         z = Py_BuildValue("(OO)", d, m);
 445         Py_XDECREF(d);
 446         Py_XDECREF(m);
 447         return z;
 448 }
 449
 450 static PyObject *
 451 complex_pow(v, w, z)
 452         PyComplexObject *v;
 453         PyObject *w;
 454         PyComplexObject *z;
 455 {
 456         Py_complex p;
 457         Py_complex exponent;
 458         long int_exponent;
 459
 460         if ((PyObject *)z!=Py_None) {
 461                 PyErr_SetString(PyExc_ValueError, "complex modulo");
 462                 return NULL;
 463         }
 464
 465         PyFPE_START_PROTECT("complex_pow", return 0)
 466         errno = 0;
 467         exponent = ((PyComplexObject*)w)->cval;
 468         int_exponent = (long)exponent.real;
 469         if (exponent.imag == 0. && exponent.real == int_exponent)
 470                 p = c_powi(v->cval,int_exponent);
 471         else
 472                 p = c_pow(v->cval,exponent);
 473
 474         PyFPE_END_PROTECT(p)
 475         if (errno == ERANGE) {
 476                 PyErr_SetString(PyExc_ValueError,
 477                                 "0.0 to a negative or complex power");
 478                 return NULL;
 479         }
 480
 481         return PyComplex_FromCComplex(p);
 482 }
 483
 484 static PyObject *
 485 complex_neg(v)
 486         PyComplexObject *v;
 487 {
 488         Py_complex neg;
 489         neg.real = -v->cval.real;
 490         neg.imag = -v->cval.imag;
 491         return PyComplex_FromCComplex(neg);
 492 }
 493
 494 static PyObject *
 495 complex_pos(v)
 496         PyComplexObject *v;
 497 {
 498         Py_INCREF(v);
 499         return (PyObject *)v;
 500 }
 501
 502 static PyObject *
 503 complex_abs(v)
 504         PyComplexObject *v;
 505 {
 506         double result;
 507         PyFPE_START_PROTECT("complex_abs", return 0)
 508         result = hypot(v->cval.real,v->cval.imag);
 509         PyFPE_END_PROTECT(result)
 510         return PyFloat_FromDouble(result);
 511 }
 512
 513 static int
 514 complex_nonzero(v)
 515         PyComplexObject *v;
 516 {
 517         return v->cval.real != 0.0 || v->cval.imag != 0.0;
 518 }
 519
 520 static int
 521 complex_coerce(pv, pw)
 522         PyObject **pv;
 523         PyObject **pw;
 524 {
 525         Py_complex cval;
 526         cval.imag = 0.;
 527         if (PyInt_Check(*pw)) {
 528                 cval.real = (double)PyInt_AsLong(*pw);
 529                 *pw = PyComplex_FromCComplex(cval);
 530                 Py_INCREF(*pv);
 531                 return 0;
 532         }
 533         else if (PyLong_Check(*pw)) {
 534                 cval.real = PyLong_AsDouble(*pw);
 535                 *pw = PyComplex_FromCComplex(cval);
 536                 Py_INCREF(*pv);
 537                 return 0;
 538         }
 539         else if (PyFloat_Check(*pw)) {
 540                 cval.real = PyFloat_AsDouble(*pw);
 541                 *pw = PyComplex_FromCComplex(cval);
 542                 Py_INCREF(*pv);
 543                 return 0;
 544         }
 545         return 1; /* Can't do it */
 546 }
 547
 548 static PyObject *
 549 complex_int(v)
 550         PyObject *v;
 551 {
 552         PyErr_SetString(PyExc_TypeError,
 553                    "can't convert complex to int; use e.g. int(abs(z))");
 554         return NULL;
 555 }
 556
 557 static PyObject *
 558 complex_long(v)
 559         PyObject *v;
 560 {
 561         PyErr_SetString(PyExc_TypeError,
 562                    "can't convert complex to long; use e.g. long(abs(z))");
 563         return NULL;
 564 }
 565
 566 static PyObject *
 567 complex_float(v)
 568         PyObject *v;
 569 {
 570         PyErr_SetString(PyExc_TypeError,
 571                    "can't convert complex to float; use e.g. abs(z)");
 572         return NULL;
 573 }
 574
 575 static PyObject *
 576 complex_conjugate(self, args)
 577         PyObject *self;
 578         PyObject *args;
 579 {
 580         Py_complex c;
 581         if (!PyArg_ParseTuple(args, ""))
 582                 return NULL;
 583         c = ((PyComplexObject *)self)->cval;
 584         c.imag = -c.imag;
 585         return PyComplex_FromCComplex(c);
 586 }
 587
 588 static PyMethodDef complex_methods[] = {
 589         {"conjugate",   complex_conjugate,      1},
 590         {NULL,          NULL}           /* sentinel */
 591 };
 592
 593
 594 static PyObject *
 595 complex_getattr(self, name)
 596         PyComplexObject *self;
 597         char *name;
 598 {
 599         if (strcmp(name, "real") == 0)
 600                 return (PyObject *)PyFloat_FromDouble(self->cval.real);
 601         else if (strcmp(name, "imag") == 0)
 602                 return (PyObject *)PyFloat_FromDouble(self->cval.imag);
 603         else if (strcmp(name, "__members__") == 0)
 604                 return Py_BuildValue("[ss]", "imag", "real");
 605         return Py_FindMethod(complex_methods, (PyObject *)self, name);
 606 }
 607
 608 static PyNumberMethods complex_as_number = {
 609         (binaryfunc)complex_add, /*nb_add*/
 610         (binaryfunc)complex_sub, /*nb_subtract*/
 611         (binaryfunc)complex_mul, /*nb_multiply*/
 612         (binaryfunc)complex_div, /*nb_divide*/
 613         (binaryfunc)complex_remainder,  /*nb_remainder*/
 614         (binaryfunc)complex_divmod,     /*nb_divmod*/
 615         (ternaryfunc)complex_pow, /*nb_power*/
 616         (unaryfunc)complex_neg, /*nb_negative*/
 617         (unaryfunc)complex_pos, /*nb_positive*/
 618         (unaryfunc)complex_abs, /*nb_absolute*/
 619         (inquiry)complex_nonzero, /*nb_nonzero*/
 620         0,              /*nb_invert*/
 621         0,              /*nb_lshift*/
 622         0,              /*nb_rshift*/
 623         0,              /*nb_and*/
 624         0,              /*nb_xor*/
 625         0,              /*nb_or*/
 626         (coercion)complex_coerce, /*nb_coerce*/
 627         (unaryfunc)complex_int, /*nb_int*/
 628         (unaryfunc)complex_long, /*nb_long*/
 629         (unaryfunc)complex_float, /*nb_float*/
 630         0,              /*nb_oct*/
 631         0,              /*nb_hex*/
 632 };
 633
 634 PyTypeObject PyComplex_Type = {
 635         PyObject_HEAD_INIT(&PyType_Type)
 636         0,
 637         "complex",
 638         sizeof(PyComplexObject),
 639         0,
 640         (destructor)complex_dealloc,    /*tp_dealloc*/
 641         (printfunc)complex_print,       /*tp_print*/
 642         (getattrfunc)complex_getattr,   /*tp_getattr*/
 643         0,                              /*tp_setattr*/
 644         (cmpfunc)complex_compare,       /*tp_compare*/
 645         (reprfunc)complex_repr,         /*tp_repr*/
 646         &complex_as_number,             /*tp_as_number*/
 647         0,                              /*tp_as_sequence*/
 648         0,                              /*tp_as_mapping*/
 649         (hashfunc)complex_hash,         /*tp_hash*/
 650 };
 651
 652 #endif