home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / py2s152.zip / Objects / complexobject.c < prev    next >
C/C++ Source or Header  |  1999-06-27  |  15KB  |  653 lines
  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(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
  653.