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