Objects/complexobject.c

   1 /* Complex object implementation */
   2
   3 /* Borrows heavily from floatobject.c */
   4
   5 #ifndef WITHOUT_COMPLEX
   6
   7 #include "allobjects.h"
   8 #include "modsupport.h"
   9
  10 #include <errno.h>
  11 #include "mymath.h"
  12
  13 #ifdef i860
  14 /* Cray APP has bogus definition of HUGE_VAL in <math.h> */
  15 #undef HUGE_VAL
  16 #endif
  17
  18 #ifdef HUGE_VAL
  19 #define CHECK(x) if (errno != 0) ; \
  20         else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
  21         else errno = ERANGE
  22 #else
  23 #define CHECK(x) /* Don't know how to check */
  24 #endif
  25
  26 #ifdef HAVE_LIMITS_H
  27 #include <limits.h>
  28 #endif
  29
  30 #ifndef LONG_MAX
  31 #define LONG_MAX 0X7FFFFFFFL
  32 #endif
  33
  34 #ifndef LONG_MIN
  35 #define LONG_MIN (-LONG_MAX-1)
  36 #endif
  37
  38 #ifdef __NeXT__
  39 #ifdef __sparc__
  40 /*
  41  * This works around a bug in the NS/Sparc 3.3 pre-release
  42  * limits.h header file.
  43  * 10-Feb-1995 bwarsaw@cnri.reston.va.us
  44  */
  45 #undef LONG_MIN
  46 #define LONG_MIN (-LONG_MAX-1)
  47 #endif
  48 #endif
  49
  50 #if !defined(__STDC__) && !defined(macintosh)
  51 extern double fmod PROTO((double, double));
  52 extern double pow PROTO((double, double));
  53 #endif
  54
  55
  56 /* elementary operations on complex numbers */
  57
  58 static int c_error;
  59 static complex c_1 = {1., 0.};
  60
  61 complex c_sum(a,b)
  62         complex a,b;
  63 {
  64         complex r;
  65         r.real = a.real + b.real;
  66         r.imag = a.imag + b.imag;
  67         return r;
  68 }
  69
  70 complex c_diff(a,b)
  71         complex a,b;
  72 {
  73         complex r;
  74         r.real = a.real - b.real;
  75         r.imag = a.imag - b.imag;
  76         return r;
  77 }
  78
  79 complex c_neg(a)
  80         complex a;
  81 {
  82         complex r;
  83         r.real = -a.real;
  84         r.imag = -a.imag;
  85         return r;
  86 }
  87
  88 complex c_prod(a,b)
  89         complex a,b;
  90 {
  91         complex r;
  92         r.real = a.real*b.real - a.imag*b.imag;
  93         r.imag = a.real*b.imag + a.imag*b.real;
  94         return r;
  95 }
  96
  97 complex c_quot(a,b)
  98         complex a,b;
  99 {
 100         complex r;
 101         double d = b.real*b.real + b.imag*b.imag;
 102         if (d == 0.)
 103                 c_error = 1;
 104         r.real = (a.real*b.real + a.imag*b.imag)/d;
 105         r.imag = (a.imag*b.real - a.real*b.imag)/d;
 106         return r;
 107 }
 108
 109 complex c_pow(a,b)
 110         complex a,b;
 111 {
 112         complex r;
 113         double vabs,len,at,phase;
 114         if (b.real == 0. && b.imag == 0.) {
 115                 r.real = 1.;
 116                 r.imag = 0.;
 117         }
 118         else if (a.real == 0. && a.imag == 0.) {
 119                 if (b.imag != 0. || b.real < 0.)
 120                         c_error = 2;
 121                 r.real = 0.;
 122                 r.imag = 0.;
 123         }
 124         else {
 125                 vabs = hypot(a.real,a.imag);
 126                 len = pow(vabs,b.real);
 127                 at = atan2(a.imag, a.real);
 128                 phase = at*b.real;
 129                 if (b.imag != 0.0) {
 130                         len /= exp(at*b.imag);
 131                         phase += b.imag*log(vabs);
 132                 }
 133                 r.real = len*cos(phase);
 134                 r.imag = len*sin(phase);
 135         }
 136         return r;
 137 }
 138
 139 static complex c_powu(x, n)
 140         complex x;
 141         long n;
 142 {
 143         complex r = c_1;
 144         complex p = x;
 145         long mask = 1;
 146         while (mask > 0 && n >= mask) {
 147                 if (n & mask)
 148                         r = c_prod(r,p);
 149                 mask <<= 1;
 150                 p = c_prod(p,p);
 151         }
 152         return r;
 153 }
 154
 155 static complex c_powi(x, n)
 156         complex x;
 157         long n;
 158 {
 159         complex cn;
 160
 161         if (n > 100 || n < -100) {
 162                 cn.real = (double) n;
 163                 cn.imag = 0.;
 164                 return c_pow(x,cn);
 165         }
 166         else if (n > 0)
 167                 return c_powu(x,n);
 168         else
 169                 return c_quot(c_1,c_powu(x,-n));
 170
 171 }
 172
 173 PyObject *
 174 PyComplex_FromCComplex(complex cval)
 175 {
 176         register complexobject *op = (complexobject *) malloc(sizeof(complexobject));
 177         if (op == NULL)
 178                 return err_nomem();
 179         op->ob_type = &Complextype;
 180         op->cval = cval;
 181         NEWREF(op);
 182         return (object *) op;
 183 }
 184
 185 PyObject *
 186 PyComplex_FromDoubles(double real, double imag) {
 187   complex c;
 188   c.real = real;
 189   c.imag = imag;
 190   return PyComplex_FromCComplex(c);
 191 }
 192
 193 double
 194 PyComplex_RealAsDouble(PyObject *op) {
 195   if (PyComplex_Check(op)) {
 196     return ((PyComplexObject *)op)->cval.real;
 197   } else {
 198     return PyFloat_AsDouble(op);
 199   }
 200 }
 201
 202 double
 203 PyComplex_ImagAsDouble(PyObject *op) {
 204   if (PyComplex_Check(op)) {
 205     return ((PyComplexObject *)op)->cval.imag;
 206   } else {
 207     return 0.0;
 208   }
 209 }
 210
 211 complex
 212 PyComplex_AsCComplex(PyObject *op) {
 213         complex cv;
 214         if (PyComplex_Check(op)) {
 215                 return ((PyComplexObject *)op)->cval;
 216         } else {
 217                 cv.real = PyFloat_AsDouble(op);
 218                 cv.imag = 0.;
 219                 return cv;
 220         }
 221 }
 222
 223 static void
 224 complex_dealloc(op)
 225         object *op;
 226 {
 227         DEL(op);
 228 }
 229
 230
 231 static void
 232 complex_buf_repr(buf, v)
 233         char *buf;
 234         complexobject *v;
 235 {
 236         if (v->cval.real == 0.)
 237                 sprintf(buf, "%.12gj", v->cval.imag);
 238         else
 239                 sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
 240 }
 241
 242 static int
 243 complex_print(v, fp, flags)
 244         complexobject *v;
 245         FILE *fp;
 246         int flags; /* Not used but required by interface */
 247 {
 248         char buf[100];
 249         complex_buf_repr(buf, v);
 250         fputs(buf, fp);
 251         return 0;
 252 }
 253
 254 static object *
 255 complex_repr(v)
 256         complexobject *v;
 257 {
 258         char buf[100];
 259         complex_buf_repr(buf, v);
 260         return newstringobject(buf);
 261 }
 262
 263 static int
 264 complex_compare(v, w)
 265         complexobject *v, *w;
 266 {
 267 /* Note: "greater" and "smaller" have no meaning for complex numbers,
 268    but Python requires that they be defined nevertheless. */
 269         complex i = v->cval;
 270         complex j = w->cval;
 271         if (i.real == j.real && i.imag == j.imag)
 272            return 0;
 273         else if (i.real != j.real)
 274            return (i.real < j.real) ? -1 : 1;
 275         else
 276            return (i.imag < j.imag) ? -1 : 1;
 277 }
 278
 279 static long
 280 complex_hash(v)
 281         complexobject *v;
 282 {
 283         double intpart, fractpart;
 284         int expo;
 285         long x;
 286         /* This is designed so that Python numbers with the same
 287            value hash to the same value, otherwise comparisons
 288            of mapping keys will turn out weird */
 289
 290 #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
 291 {
 292         extended e;
 293         fractpart = modf(v->cval.real, &e);
 294         intpart = e;
 295 }
 296 #else
 297         fractpart = modf(v->cval.real, &intpart);
 298 #endif
 299
 300         if (fractpart == 0.0) {
 301                 if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
 302                         /* Convert to long int and use its hash... */
 303                         object *w = dnewlongobject(v->cval.real);
 304                         if (w == NULL)
 305                                 return -1;
 306                         x = hashobject(w);
 307                         DECREF(w);
 308                         return x;
 309                 }
 310                 x = (long)intpart;
 311         }
 312         else {
 313                 fractpart = frexp(fractpart, &expo);
 314                 fractpart = fractpart*2147483648.0; /* 2**31 */
 315                 x = (long) (intpart + fractpart) ^ expo; /* Rather arbitrary */
 316         }
 317         if (x == -1)
 318                 x = -2;
 319         return x;
 320 }
 321
 322 static object *
 323 complex_add(v, w)
 324         complexobject *v;
 325         complexobject *w;
 326 {
 327         return newcomplexobject(c_sum(v->cval,w->cval));
 328 }
 329
 330 static object *
 331 complex_sub(v, w)
 332         complexobject *v;
 333         complexobject *w;
 334 {
 335         return newcomplexobject(c_diff(v->cval,w->cval));
 336 }
 337
 338 static object *
 339 complex_mul(v, w)
 340         complexobject *v;
 341         complexobject *w;
 342 {
 343         return newcomplexobject(c_prod(v->cval,w->cval));
 344 }
 345
 346 static object *
 347 complex_div(v, w)
 348         complexobject *v;
 349         complexobject *w;
 350 {
 351         complex quot;
 352         c_error = 0;
 353         quot = c_quot(v->cval,w->cval);
 354         if (c_error == 1) {
 355                 err_setstr(ZeroDivisionError, "float division");
 356                 return NULL;
 357         }
 358         return newcomplexobject(quot);
 359 }
 360
 361
 362 static object *
 363 complex_pow(v, w, z)
 364         complexobject *v;
 365         object *w;
 366         complexobject *z;
 367 {
 368         complex p;
 369         complex exponent;
 370         long int_exponent;
 371
 372         if ((object *)z!=None) {
 373                 err_setstr(ValueError, "complex modulo");
 374                 return NULL;
 375         }
 376
 377         c_error = 0;
 378         exponent = ((complexobject*)w)->cval;
 379         int_exponent = (long)exponent.real;
 380         if (exponent.imag == 0. && exponent.real == int_exponent)
 381                 p = c_powi(v->cval,int_exponent);
 382         else
 383                 p = c_pow(v->cval,exponent);
 384
 385         if (c_error == 2) {
 386                 err_setstr(ValueError, "0.0 to a negative or complex power");
 387                 return NULL;
 388         }
 389
 390         return newcomplexobject(p);
 391 }
 392
 393 static object *
 394 complex_neg(v)
 395         complexobject *v;
 396 {
 397         complex neg;
 398         neg.real = -v->cval.real;
 399         neg.imag = -v->cval.imag;
 400         return newcomplexobject(neg);
 401 }
 402
 403 static object *
 404 complex_pos(v)
 405         complexobject *v;
 406 {
 407         INCREF(v);
 408         return (object *)v;
 409 }
 410
 411 static object *
 412 complex_abs(v)
 413         complexobject *v;
 414 {
 415         return newfloatobject(hypot(v->cval.real,v->cval.imag));
 416 }
 417
 418 static int
 419 complex_nonzero(v)
 420         complexobject *v;
 421 {
 422         return v->cval.real != 0.0 && v->cval.imag != 0.0;
 423 }
 424
 425 static int
 426 complex_coerce(pv, pw)
 427         object **pv;
 428         object **pw;
 429 {
 430         complex cval;
 431         cval.imag = 0.;
 432         if (is_intobject(*pw)) {
 433                 cval.real = (double)getintvalue(*pw);
 434                 *pw = newcomplexobject(cval);
 435                 INCREF(*pv);
 436                 return 0;
 437         }
 438         else if (is_longobject(*pw)) {
 439                 cval.real = dgetlongvalue(*pw);
 440                 *pw = newcomplexobject(cval);
 441                 INCREF(*pv);
 442                 return 0;
 443         }
 444         else if (is_floatobject(*pw)) {
 445                 cval.real = getfloatvalue(*pw);
 446                 *pw = newcomplexobject(cval);
 447                 INCREF(*pv);
 448                 return 0;
 449         }
 450         return 1; /* Can't do it */
 451 }
 452
 453 static object *
 454 complex_int(v)
 455         object *v;
 456 {
 457         double x = ((complexobject *)v)->cval.real;
 458         if (x < 0 ? (x = ceil(x)) < (double)LONG_MIN
 459                   : (x = floor(x)) > (double)LONG_MAX) {
 460                 err_setstr(OverflowError, "float too large to convert");
 461                 return NULL;
 462         }
 463         return newintobject((long)x);
 464 }
 465
 466 static object *
 467 complex_long(v)
 468         object *v;
 469 {
 470         double x = ((complexobject *)v)->cval.real;
 471         return dnewlongobject(x);
 472 }
 473
 474 static object *
 475 complex_float(v)
 476         object *v;
 477 {
 478         double x = ((complexobject *)v)->cval.real;
 479         return newfloatobject(x);
 480 }
 481
 482
 483 static object *
 484 complex_new(self, args)
 485         object *self;
 486         object *args;
 487 {
 488         complex cval;
 489
 490         cval.imag = 0.;
 491         if (!PyArg_ParseTuple(args, "d|d", &cval.real, &cval.imag))
 492                 return NULL;
 493         return newcomplexobject(cval);
 494 }
 495
 496 static object *
 497 complex_conjugate(self)
 498         object *self;
 499 {
 500         complex c = ((complexobject *)self)->cval;
 501         c.imag = -c.imag;
 502         return newcomplexobject(c);
 503 }
 504
 505 static PyMethodDef complex_methods[] = {
 506         {"conjugate",   (PyCFunction)complex_conjugate, 1},
 507         {NULL,          NULL}           /* sentinel */
 508 };
 509
 510
 511 static object *
 512 complex_getattr(self, name)
 513         complexobject *self;
 514         char *name;
 515 {
 516         complex cval;
 517         if (strcmp(name, "real") == 0)
 518                 return (object *)newfloatobject(self->cval.real);
 519         else if (strcmp(name, "imag") == 0)
 520                 return (object *)newfloatobject(self->cval.imag);
 521         else if (strcmp(name, "conj") == 0) {
 522                 cval.real = self->cval.real;
 523                 cval.imag = -self->cval.imag;
 524                 return (object *)newcomplexobject(cval);
 525         }
 526         return findmethod(complex_methods, (object *)self, name);
 527 }
 528
 529 static number_methods complex_as_number = {
 530         (binaryfunc)complex_add, /*nb_add*/
 531         (binaryfunc)complex_sub, /*nb_subtract*/
 532         (binaryfunc)complex_mul, /*nb_multiply*/
 533         (binaryfunc)complex_div, /*nb_divide*/
 534         0,              /*nb_remainder*/
 535         0,              /*nb_divmod*/
 536         (ternaryfunc)complex_pow, /*nb_power*/
 537         (unaryfunc)complex_neg, /*nb_negative*/
 538         (unaryfunc)complex_pos, /*nb_positive*/
 539         (unaryfunc)complex_abs, /*nb_absolute*/
 540         (inquiry)complex_nonzero, /*nb_nonzero*/
 541         0,              /*nb_invert*/
 542         0,              /*nb_lshift*/
 543         0,              /*nb_rshift*/
 544         0,              /*nb_and*/
 545         0,              /*nb_xor*/
 546         0,              /*nb_or*/
 547         (coercion)complex_coerce, /*nb_coerce*/
 548         (unaryfunc)complex_int, /*nb_int*/
 549         (unaryfunc)complex_long, /*nb_long*/
 550         (unaryfunc)complex_float, /*nb_float*/
 551         0,              /*nb_oct*/
 552         0,              /*nb_hex*/
 553 };
 554
 555 typeobject Complextype = {
 556         OB_HEAD_INIT(&Typetype)
 557         0,
 558         "complex",
 559         sizeof(complexobject),
 560         0,
 561         (destructor)complex_dealloc,    /*tp_dealloc*/
 562         (printfunc)complex_print,       /*tp_print*/
 563         (getattrfunc)complex_getattr,   /*tp_getattr*/
 564         0,                              /*tp_setattr*/
 565         (cmpfunc)complex_compare,       /*tp_compare*/
 566         (reprfunc)complex_repr,         /*tp_repr*/
 567         &complex_as_number,             /*tp_as_number*/
 568         0,                              /*tp_as_sequence*/
 569         0,                              /*tp_as_mapping*/
 570         (hashfunc)complex_hash,         /*tp_hash*/
 571 };
 572
 573 #endif