Modules/cmathmodule.c

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