Modules/cmathmodule.c

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