Modules/cmathmodule.c

   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 cosine 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 }