home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / py2s152.zip / Modules / cmathmodule.c < prev    next >
C/C++ Source or Header  |  1999-06-27  |  8KB  |  423 lines
  1. /* Complex math module */
  2.  
  3. /* much code borrowed from mathmodule.c */
  4.  
  5. #include "Python.h"
  6.  
  7. #include "mymath.h"
  8.  
  9. #ifdef i860
  10. /* Cray APP has bogus definition of HUGE_VAL in <math.h> */
  11. #undef HUGE_VAL
  12. #endif
  13.  
  14. #ifdef HUGE_VAL
  15. #define CHECK(x) if (errno != 0) ; \
  16.     else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
  17.     else errno = ERANGE
  18. #else
  19. #define CHECK(x) /* Don't know how to check */
  20. #endif
  21.  
  22. #ifndef M_PI
  23. #define M_PI (3.141592653589793239)
  24. #endif
  25.  
  26. /* First, the C functions that do the real work */
  27.  
  28. /* constants */
  29. static Py_complex c_1 = {1., 0.};
  30. static Py_complex c_half = {0.5, 0.};
  31. static Py_complex c_i = {0., 1.};
  32. static Py_complex c_i2 = {0., 0.5};
  33. #if 0
  34. static Py_complex c_mi = {0., -1.};
  35. static Py_complex c_pi2 = {M_PI/2., 0.};
  36. #endif
  37.  
  38. /* forward declarations */
  39. staticforward Py_complex c_log();
  40. staticforward Py_complex c_prodi();
  41. staticforward Py_complex c_sqrt();
  42.  
  43.  
  44. static Py_complex c_acos(x)
  45.     Py_complex x;
  46. {
  47.     return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
  48.             c_sqrt(c_diff(c_1,c_prod(x,x))))))));
  49. }
  50.  
  51. static char c_acos_doc [] =
  52. "acos(x)\n\
  53. \n\
  54. Return the arc cosine of x.";
  55.  
  56.  
  57. static Py_complex c_acosh(x)
  58.     Py_complex x;
  59. {
  60.     return c_log(c_sum(x,c_prod(c_i,
  61.             c_sqrt(c_diff(c_1,c_prod(x,x))))));
  62. }
  63.  
  64. static char c_acosh_doc [] =
  65. "acosh(x)\n\
  66. \n\
  67. Return the hyperbolic arccosine of x.";
  68.  
  69.  
  70. static Py_complex c_asin(x)
  71.     Py_complex x;
  72. {
  73.     return c_neg(c_prodi(c_log(c_sum(c_prod(c_i,x),
  74.             c_sqrt(c_diff(c_1,c_prod(x,x)))))));
  75. }
  76.  
  77. static char c_asin_doc [] =
  78. "asin(x)\n\
  79. \n\
  80. Return the arc sine of x.";
  81.  
  82.  
  83. static Py_complex c_asinh(x)
  84.     Py_complex x;
  85. {
  86.     /* Break up long expression for WATCOM */
  87.     Py_complex z;
  88.     z = c_sum(c_1,c_prod(x,x));
  89.     z = c_diff(c_sqrt(z),x);
  90.     return c_neg(c_log(z));
  91. }
  92.  
  93. static char c_asinh_doc [] =
  94. "asinh(x)\n\
  95. \n\
  96. Return the hyperbolic arc sine of x.";
  97.  
  98.  
  99. static Py_complex c_atan(x)
  100.     Py_complex x;
  101. {
  102.     return c_prod(c_i2,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
  103. }
  104.  
  105. static char c_atan_doc [] =
  106. "atan(x)\n\
  107. \n\
  108. Return the arc tangent of x.";
  109.  
  110.  
  111. static Py_complex c_atanh(x)
  112.     Py_complex x;
  113. {
  114.     return c_prod(c_half,c_log(c_quot(c_sum(c_1,x),c_diff(c_1,x))));
  115. }
  116.  
  117. static char c_atanh_doc [] =
  118. "atanh(x)\n\
  119. \n\
  120. Return the hyperbolic arc tangent of x.";
  121.  
  122.  
  123. static Py_complex c_cos(x)
  124.     Py_complex x;
  125. {
  126.     Py_complex r;
  127.     r.real = cos(x.real)*cosh(x.imag);
  128.     r.imag = -sin(x.real)*sinh(x.imag);
  129.     return r;
  130. }
  131.  
  132. static char c_cos_doc [] =
  133. "cos(x)\n\
  134. \n\
  135. Return the cosine of x.";
  136.  
  137.  
  138. static Py_complex c_cosh(x)
  139.     Py_complex x;
  140. {
  141.     Py_complex r;
  142.     r.real = cos(x.imag)*cosh(x.real);
  143.     r.imag = sin(x.imag)*sinh(x.real);
  144.     return r;
  145. }
  146.  
  147. static char c_cosh_doc [] =
  148. "cosh(x)\n\
  149. \n\
  150. Return the hyperbolic cosine of x.";
  151.  
  152.  
  153. static Py_complex c_exp(x)
  154.     Py_complex x;
  155. {
  156.     Py_complex r;
  157.     double l = exp(x.real);
  158.     r.real = l*cos(x.imag);
  159.     r.imag = l*sin(x.imag);
  160.     return r;
  161. }
  162.  
  163. static char c_exp_doc [] =
  164. "exp(x)\n\
  165. \n\
  166. Return the exponential value e**x.";
  167.  
  168.  
  169. static Py_complex c_log(x)
  170.     Py_complex x;
  171. {
  172.     Py_complex r;
  173.     double l = hypot(x.real,x.imag);
  174.     r.imag = atan2(x.imag, x.real);
  175.     r.real = log(l);
  176.     return r;
  177. }
  178.  
  179. static char c_log_doc [] =
  180. "log(x)\n\
  181. \n\
  182. Return the natural logarithm of x.";
  183.  
  184.  
  185. static Py_complex c_log10(x)
  186.     Py_complex x;
  187. {
  188.     Py_complex r;
  189.     double l = hypot(x.real,x.imag);
  190.     r.imag = atan2(x.imag, x.real)/log(10.);
  191.     r.real = log10(l);
  192.     return r;
  193. }
  194.  
  195. static char c_log10_doc [] =
  196. "log10(x)\n\
  197. \n\
  198. Return the base-10 logarithm of x.";
  199.  
  200.  
  201. /* internal function not available from Python */
  202. static Py_complex c_prodi(x)
  203.      Py_complex x;
  204. {
  205.     Py_complex r;
  206.     r.real = -x.imag;
  207.     r.imag = x.real;
  208.     return r;
  209. }
  210.  
  211.  
  212. static Py_complex c_sin(x)
  213.     Py_complex x;
  214. {
  215.     Py_complex r;
  216.     r.real = sin(x.real)*cosh(x.imag);
  217.     r.imag = cos(x.real)*sinh(x.imag);
  218.     return r;
  219. }
  220.  
  221. static char c_sin_doc [] =
  222. "sin(x)\n\
  223. \n\
  224. Return the sine of x.";
  225.  
  226.  
  227. static Py_complex c_sinh(x)
  228.     Py_complex x;
  229. {
  230.     Py_complex r;
  231.     r.real = cos(x.imag)*sinh(x.real);
  232.     r.imag = sin(x.imag)*cosh(x.real);
  233.     return r;
  234. }
  235.  
  236. static char c_sinh_doc [] =
  237. "sinh(x)\n\
  238. \n\
  239. Return the hyperbolic sine of x.";
  240.  
  241.  
  242. static Py_complex c_sqrt(x)
  243.     Py_complex x;
  244. {
  245.     Py_complex r;
  246.     double s,d;
  247.     if (x.real == 0. && x.imag == 0.)
  248.         r = x;
  249.     else {
  250.         s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
  251.         d = 0.5*x.imag/s;
  252.         if (x.real > 0.) {
  253.             r.real = s;
  254.             r.imag = d;
  255.         }
  256.         else if (x.imag >= 0.) {
  257.             r.real = d;
  258.             r.imag = s;
  259.         }
  260.         else {
  261.             r.real = -d;
  262.             r.imag = -s;
  263.         }
  264.     }
  265.     return r;
  266. }
  267.  
  268. static char c_sqrt_doc [] =
  269. "sqrt(x)\n\
  270. \n\
  271. Return the square root of x.";
  272.  
  273.  
  274. static Py_complex c_tan(x)
  275.     Py_complex x;
  276. {
  277.     Py_complex r;
  278.     double sr,cr,shi,chi;
  279.     double rs,is,rc,ic;
  280.     double d;
  281.     sr = sin(x.real);
  282.     cr = cos(x.real);
  283.     shi = sinh(x.imag);
  284.     chi = cosh(x.imag);
  285.     rs = sr*chi;
  286.     is = cr*shi;
  287.     rc = cr*chi;
  288.     ic = -sr*shi;
  289.     d = rc*rc + ic*ic;
  290.     r.real = (rs*rc+is*ic)/d;
  291.     r.imag = (is*rc-rs*ic)/d;
  292.     return r;
  293. }
  294.  
  295. static char c_tan_doc [] =
  296. "tan(x)\n\
  297. \n\
  298. Return the tangent of x.";
  299.  
  300.  
  301. static Py_complex c_tanh(x)
  302.     Py_complex x;
  303. {
  304.     Py_complex r;
  305.     double si,ci,shr,chr;
  306.     double rs,is,rc,ic;
  307.     double d;
  308.     si = sin(x.imag);
  309.     ci = cos(x.imag);
  310.     shr = sinh(x.real);
  311.     chr = cosh(x.real);
  312.     rs = ci*shr;
  313.     is = si*chr;
  314.     rc = ci*chr;
  315.     ic = si*shr;
  316.     d = rc*rc + ic*ic;
  317.     r.real = (rs*rc+is*ic)/d;
  318.     r.imag = (is*rc-rs*ic)/d;
  319.     return r;
  320. }
  321.  
  322. static char c_tanh_doc [] =
  323. "tanh(x)\n\
  324. \n\
  325. Return the hyperbolic tangent of x.";
  326.  
  327.  
  328. /* And now the glue to make them available from Python: */
  329.  
  330. static PyObject *
  331. math_error()
  332. {
  333.     if (errno == EDOM)
  334.         PyErr_SetString(PyExc_ValueError, "math domain error");
  335.     else if (errno == ERANGE)
  336.         PyErr_SetString(PyExc_OverflowError, "math range error");
  337.     else    /* Unexpected math error */
  338.         PyErr_SetFromErrno(PyExc_ValueError); 
  339.     return NULL;
  340. }
  341.  
  342. static PyObject *
  343. math_1(args, func)
  344.     PyObject *args;
  345.     Py_complex (*func) Py_FPROTO((Py_complex));
  346. {
  347.     Py_complex x;
  348.     if (!PyArg_ParseTuple(args, "D", &x))
  349.         return NULL;
  350.     errno = 0;
  351.     PyFPE_START_PROTECT("complex function", return 0)
  352.     x = (*func)(x);
  353.     PyFPE_END_PROTECT(x)
  354.     CHECK(x.real);
  355.     CHECK(x.imag);
  356.     if (errno != 0)
  357.         return math_error();
  358.     else
  359.         return PyComplex_FromCComplex(x);
  360. }
  361.  
  362. #define FUNC1(stubname, func) \
  363.     static PyObject * stubname(self, args) PyObject *self, *args; { \
  364.         return math_1(args, func); \
  365.     }
  366.  
  367. FUNC1(cmath_acos, c_acos)
  368. FUNC1(cmath_acosh, c_acosh)
  369. FUNC1(cmath_asin, c_asin)
  370. FUNC1(cmath_asinh, c_asinh)
  371. FUNC1(cmath_atan, c_atan)
  372. FUNC1(cmath_atanh, c_atanh)
  373. FUNC1(cmath_cos, c_cos)
  374. FUNC1(cmath_cosh, c_cosh)
  375. FUNC1(cmath_exp, c_exp)
  376. FUNC1(cmath_log, c_log)
  377. FUNC1(cmath_log10, c_log10)
  378. FUNC1(cmath_sin, c_sin)
  379. FUNC1(cmath_sinh, c_sinh)
  380. FUNC1(cmath_sqrt, c_sqrt)
  381. FUNC1(cmath_tan, c_tan)
  382. FUNC1(cmath_tanh, c_tanh)
  383.  
  384.  
  385. static char module_doc [] =
  386. "This module is always available. It provides access to mathematical\n\
  387. functions for complex numbers.";
  388.  
  389.  
  390. static PyMethodDef cmath_methods[] = {
  391.     {"acos", cmath_acos, 1, c_acos_doc},
  392.     {"acosh", cmath_acosh, 1, c_acosh_doc},
  393.     {"asin", cmath_asin, 1, c_asin_doc},
  394.     {"asinh", cmath_asinh, 1, c_asinh_doc},
  395.     {"atan", cmath_atan, 1, c_atan_doc},
  396.     {"atanh", cmath_atanh, 1, c_atanh_doc},
  397.     {"cos", cmath_cos, 1, c_cos_doc},
  398.     {"cosh", cmath_cosh, 1, c_cosh_doc},
  399.     {"exp", cmath_exp, 1, c_exp_doc},
  400.     {"log", cmath_log, 1, c_log_doc},
  401.     {"log10", cmath_log10, 1, c_log10_doc},
  402.     {"sin", cmath_sin, 1, c_sin_doc},
  403.     {"sinh", cmath_sinh, 1, c_sinh_doc},
  404.     {"sqrt", cmath_sqrt, 1, c_sqrt_doc},
  405.     {"tan", cmath_tan, 1, c_tan_doc},
  406.     {"tanh", cmath_tanh, 1, c_tanh_doc},
  407.     {NULL,        NULL}        /* sentinel */
  408. };
  409.  
  410. DL_EXPORT(void)
  411. initcmath()
  412. {
  413.     PyObject *m, *d, *v;
  414.     
  415.     m = Py_InitModule3("cmath", cmath_methods, module_doc);
  416.     d = PyModule_GetDict(m);
  417.     PyDict_SetItemString(d, "pi",
  418.                  v = PyFloat_FromDouble(atan(1.0) * 4.0));
  419.     Py_DECREF(v);
  420.     PyDict_SetItemString(d, "e", v = PyFloat_FromDouble(exp(1.0)));
  421.     Py_DECREF(v);
  422. }
  423.