home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 15 / AACD15.ISO / AACD / Programming / Python2 / Python20_source / Objects / complexobject.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-10-25  |  10.8 KB  |  521 lines
  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 int
  214. complex_compare(PyComplexObject *v, PyComplexObject *w)
  215. {
  216.     /* Note: "greater" and "smaller" have no meaning for complex numbers,
  217.        but Python requires that they be defined nevertheless. */
  218.     Py_complex i, j;
  219.     i = v->cval;
  220.     j = w->cval;
  221.     if (i.real == j.real && i.imag == j.imag)
  222.         return 0;
  223.     else if (i.real != j.real)
  224.         return (i.real < j.real) ? -1 : 1;
  225.     else
  226.         return (i.imag < j.imag) ? -1 : 1;
  227. }
  228.  
  229. static long
  230. complex_hash(PyComplexObject *v)
  231. {
  232.     long hashreal, hashimag, combined;
  233.     hashreal = _Py_HashDouble(v->cval.real);
  234.     if (hashreal == -1)
  235.         return -1;
  236.     hashimag = _Py_HashDouble(v->cval.imag);
  237.     if (hashimag == -1)
  238.         return -1;
  239.     /* Note:  if the imaginary part is 0, hashimag is 0 now,
  240.      * so the following returns hashreal unchanged.  This is
  241.      * important because numbers of different types that
  242.      * compare equal must have the same hash value, so that
  243.      * hash(x + 0*j) must equal hash(x).
  244.      */
  245.     combined = hashreal + 1000003 * hashimag;
  246.     if (combined == -1)
  247.         combined = -2;
  248.     return combined;
  249. }
  250.  
  251. static PyObject *
  252. complex_add(PyComplexObject *v, PyComplexObject *w)
  253. {
  254.     Py_complex result;
  255.     PyFPE_START_PROTECT("complex_add", return 0)
  256.     result = c_sum(v->cval,w->cval);
  257.     PyFPE_END_PROTECT(result)
  258.     return PyComplex_FromCComplex(result);
  259. }
  260.  
  261. static PyObject *
  262. complex_sub(PyComplexObject *v, PyComplexObject *w)
  263. {
  264.     Py_complex result;
  265.     PyFPE_START_PROTECT("complex_sub", return 0)
  266.     result = c_diff(v->cval,w->cval);
  267.     PyFPE_END_PROTECT(result)
  268.     return PyComplex_FromCComplex(result);
  269. }
  270.  
  271. static PyObject *
  272. complex_mul(PyComplexObject *v, PyComplexObject *w)
  273. {
  274.     Py_complex result;
  275.     PyFPE_START_PROTECT("complex_mul", return 0)
  276.     result = c_prod(v->cval,w->cval);
  277.     PyFPE_END_PROTECT(result)
  278.     return PyComplex_FromCComplex(result);
  279. }
  280.  
  281. static PyObject *
  282. complex_div(PyComplexObject *v, PyComplexObject *w)
  283. {
  284.     Py_complex quot;
  285.     PyFPE_START_PROTECT("complex_div", return 0)
  286.     errno = 0;
  287.     quot = c_quot(v->cval,w->cval);
  288.     PyFPE_END_PROTECT(quot)
  289.     if (errno == EDOM) {
  290.         PyErr_SetString(PyExc_ZeroDivisionError, "complex division");
  291.         return NULL;
  292.     }
  293.     return PyComplex_FromCComplex(quot);
  294. }
  295.  
  296. static PyObject *
  297. complex_remainder(PyComplexObject *v, PyComplexObject *w)
  298. {
  299.         Py_complex div, mod;
  300.     errno = 0;
  301.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  302.     if (errno == EDOM) {
  303.         PyErr_SetString(PyExc_ZeroDivisionError, "complex remainder");
  304.         return NULL;
  305.     }
  306.     div.real = floor(div.real); /* Use the floor of the real part. */
  307.     div.imag = 0.0;
  308.     mod = c_diff(v->cval, c_prod(w->cval, div));
  309.  
  310.     return PyComplex_FromCComplex(mod);
  311. }
  312.  
  313.  
  314. static PyObject *
  315. complex_divmod(PyComplexObject *v, PyComplexObject *w)
  316. {
  317.         Py_complex div, mod;
  318.     PyObject *d, *m, *z;
  319.     errno = 0;
  320.     div = c_quot(v->cval,w->cval); /* The raw divisor value. */
  321.     if (errno == EDOM) {
  322.         PyErr_SetString(PyExc_ZeroDivisionError, "complex divmod()");
  323.         return NULL;
  324.     }
  325.     div.real = floor(div.real); /* Use the floor of the real part. */
  326.     div.imag = 0.0;
  327.     mod = c_diff(v->cval, c_prod(w->cval, div));
  328.     d = PyComplex_FromCComplex(div);
  329.     m = PyComplex_FromCComplex(mod);
  330.     z = Py_BuildValue("(OO)", d, m);
  331.     Py_XDECREF(d);
  332.     Py_XDECREF(m);
  333.     return z;
  334. }
  335.  
  336. static PyObject *
  337. complex_pow(PyComplexObject *v, PyObject *w, PyComplexObject *z)
  338. {
  339.     Py_complex p;
  340.     Py_complex exponent;
  341.     long int_exponent;
  342.  
  343.      if ((PyObject *)z!=Py_None) {
  344.         PyErr_SetString(PyExc_ValueError, "complex modulo");
  345.         return NULL;
  346.     }
  347.     PyFPE_START_PROTECT("complex_pow", return 0)
  348.     errno = 0;
  349.     exponent = ((PyComplexObject*)w)->cval;
  350.     int_exponent = (long)exponent.real;
  351.     if (exponent.imag == 0. && exponent.real == int_exponent)
  352.         p = c_powi(v->cval,int_exponent);
  353.     else
  354.         p = c_pow(v->cval,exponent);
  355.  
  356.     PyFPE_END_PROTECT(p)
  357.     if (errno == ERANGE) {
  358.         PyErr_SetString(PyExc_ValueError,
  359.                 "0.0 to a negative or complex power");
  360.         return NULL;
  361.     }
  362.     return PyComplex_FromCComplex(p);
  363. }
  364.  
  365. static PyObject *
  366. complex_neg(PyComplexObject *v)
  367. {
  368.     Py_complex neg;
  369.     neg.real = -v->cval.real;
  370.     neg.imag = -v->cval.imag;
  371.     return PyComplex_FromCComplex(neg);
  372. }
  373.  
  374. static PyObject *
  375. complex_pos(PyComplexObject *v)
  376. {
  377.     Py_INCREF(v);
  378.     return (PyObject *)v;
  379. }
  380.  
  381. static PyObject *
  382. complex_abs(PyComplexObject *v)
  383. {
  384.     double result;
  385.     PyFPE_START_PROTECT("complex_abs", return 0)
  386.     result = hypot(v->cval.real,v->cval.imag);
  387.     PyFPE_END_PROTECT(result)
  388.     return PyFloat_FromDouble(result);
  389. }
  390.  
  391. static int
  392. complex_nonzero(PyComplexObject *v)
  393. {
  394.     return v->cval.real != 0.0 || v->cval.imag != 0.0;
  395. }
  396.  
  397. static int
  398. complex_coerce(PyObject **pv, PyObject **pw)
  399. {
  400.     Py_complex cval;
  401.     cval.imag = 0.;
  402.     if (PyInt_Check(*pw)) {
  403.         cval.real = (double)PyInt_AsLong(*pw);
  404.         *pw = PyComplex_FromCComplex(cval);
  405.         Py_INCREF(*pv);
  406.         return 0;
  407.     }
  408.     else if (PyLong_Check(*pw)) {
  409.         cval.real = PyLong_AsDouble(*pw);
  410.         *pw = PyComplex_FromCComplex(cval);
  411.         Py_INCREF(*pv);
  412.         return 0;
  413.     }
  414.     else if (PyFloat_Check(*pw)) {
  415.         cval.real = PyFloat_AsDouble(*pw);
  416.         *pw = PyComplex_FromCComplex(cval);
  417.         Py_INCREF(*pv);
  418.         return 0;
  419.     }
  420.     return 1; /* Can't do it */
  421. }
  422.  
  423. static PyObject *
  424. complex_int(PyObject *v)
  425. {
  426.     PyErr_SetString(PyExc_TypeError,
  427.            "can't convert complex to int; use e.g. int(abs(z))");
  428.     return NULL;
  429. }
  430.  
  431. static PyObject *
  432. complex_long(PyObject *v)
  433. {
  434.     PyErr_SetString(PyExc_TypeError,
  435.            "can't convert complex to long; use e.g. long(abs(z))");
  436.     return NULL;
  437. }
  438.  
  439. static PyObject *
  440. complex_float(PyObject *v)
  441. {
  442.     PyErr_SetString(PyExc_TypeError,
  443.            "can't convert complex to float; use e.g. abs(z)");
  444.     return NULL;
  445. }
  446.  
  447. static PyObject *
  448. complex_conjugate(PyObject *self, PyObject *args)
  449. {
  450.     Py_complex c;
  451.     if (!PyArg_ParseTuple(args, ":conjugate"))
  452.         return NULL;
  453.     c = ((PyComplexObject *)self)->cval;
  454.     c.imag = -c.imag;
  455.     return PyComplex_FromCComplex(c);
  456. }
  457.  
  458. static PyMethodDef complex_methods[] = {
  459.     {"conjugate",    complex_conjugate,    1},
  460.     {NULL,        NULL}        /* sentinel */
  461. };
  462.  
  463.  
  464. static PyObject *
  465. complex_getattr(PyComplexObject *self, char *name)
  466. {
  467.     if (strcmp(name, "real") == 0)
  468.         return (PyObject *)PyFloat_FromDouble(self->cval.real);
  469.     else if (strcmp(name, "imag") == 0)
  470.         return (PyObject *)PyFloat_FromDouble(self->cval.imag);
  471.     else if (strcmp(name, "__members__") == 0)
  472.         return Py_BuildValue("[ss]", "imag", "real");
  473.     return Py_FindMethod(complex_methods, (PyObject *)self, name);
  474. }
  475.  
  476. static PyNumberMethods complex_as_number = {
  477.     (binaryfunc)complex_add, /*nb_add*/
  478.     (binaryfunc)complex_sub, /*nb_subtract*/
  479.     (binaryfunc)complex_mul, /*nb_multiply*/
  480.     (binaryfunc)complex_div, /*nb_divide*/
  481.     (binaryfunc)complex_remainder,    /*nb_remainder*/
  482.     (binaryfunc)complex_divmod,    /*nb_divmod*/
  483.     (ternaryfunc)complex_pow, /*nb_power*/
  484.     (unaryfunc)complex_neg, /*nb_negative*/
  485.     (unaryfunc)complex_pos, /*nb_positive*/
  486.     (unaryfunc)complex_abs, /*nb_absolute*/
  487.     (inquiry)complex_nonzero, /*nb_nonzero*/
  488.     0,        /*nb_invert*/
  489.     0,        /*nb_lshift*/
  490.     0,        /*nb_rshift*/
  491.     0,        /*nb_and*/
  492.     0,        /*nb_xor*/
  493.     0,        /*nb_or*/
  494.     (coercion)complex_coerce, /*nb_coerce*/
  495.     (unaryfunc)complex_int, /*nb_int*/
  496.     (unaryfunc)complex_long, /*nb_long*/
  497.     (unaryfunc)complex_float, /*nb_float*/
  498.     0,        /*nb_oct*/
  499.     0,        /*nb_hex*/
  500. };
  501.  
  502. PyTypeObject PyComplex_Type = {
  503.     PyObject_HEAD_INIT(&PyType_Type)
  504.     0,
  505.     "complex",
  506.     sizeof(PyComplexObject),
  507.     0,
  508.     (destructor)complex_dealloc,    /*tp_dealloc*/
  509.     (printfunc)complex_print,    /*tp_print*/
  510.     (getattrfunc)complex_getattr,    /*tp_getattr*/
  511.     0,                /*tp_setattr*/
  512.     (cmpfunc)complex_compare,    /*tp_compare*/
  513.     (reprfunc)complex_repr,        /*tp_repr*/
  514.     &complex_as_number,            /*tp_as_number*/
  515.     0,                /*tp_as_sequence*/
  516.     0,                /*tp_as_mapping*/
  517.     (hashfunc)complex_hash,     /*tp_hash*/
  518. };
  519.  
  520. #endif
  521.