Modules/mathmodule.c

   1 /* Math module -- standard C math library functions, pi and e */
   2
   3 #include "Python.h"
   4 #include "longintrepr.h"
   5
   6 #ifndef _MSC_VER
   7 #ifndef __STDC__
   8 extern double fmod (double, double);
   9 extern double frexp (double, int *);
  10 extern double ldexp (double, int);
  11 extern double modf (double, double *);
  12 #endif /* __STDC__ */
  13 #endif /* _MSC_VER */
  14
  15 /* Call is_error when errno != 0, and where x is the result libm
  16  * returned.  is_error will usually set up an exception and return
  17  * true (1), but may return false (0) without setting up an exception.
  18  */
  19 static int
  20 is_error(double x)
  21 {
  22         int result = 1; /* presumption of guilt */
  23         assert(errno);  /* non-zero errno is a precondition for calling */
  24         if (errno == EDOM)
  25                 PyErr_SetString(PyExc_ValueError, "math domain error");
  26
  27         else if (errno == ERANGE) {
  28                 /* ANSI C generally requires libm functions to set ERANGE
  29                  * on overflow, but also generally *allows* them to set
  30                  * ERANGE on underflow too.  There's no consistency about
  31                  * the latter across platforms.
  32                  * Alas, C99 never requires that errno be set.
  33                  * Here we suppress the underflow errors (libm functions
  34                  * should return a zero on underflow, and +- HUGE_VAL on
  35                  * overflow, so testing the result for zero suffices to
  36                  * distinguish the cases).
  37                  */
  38                 if (x)
  39                         PyErr_SetString(PyExc_OverflowError,
  40                                         "math range error");
  41                 else
  42                         result = 0;
  43         }
  44         else
  45                 /* Unexpected math error */
  46                 PyErr_SetFromErrno(PyExc_ValueError);
  47         return result;
  48 }
  49
  50 static PyObject *
  51 math_1(PyObject *args, double (*func) (double), char *argsfmt)
  52 {
  53         double x;
  54         if (!  PyArg_ParseTuple(args, argsfmt, &x))
  55                 return NULL;
  56         errno = 0;
  57         PyFPE_START_PROTECT("in math_1", return 0)
  58         x = (*func)(x);
  59         PyFPE_END_PROTECT(x)
  60         Py_SET_ERANGE_IF_OVERFLOW(x);
  61         if (errno && is_error(x))
  62                 return NULL;
  63         else
  64                 return PyFloat_FromDouble(x);
  65 }
  66
  67 static PyObject *
  68 math_2(PyObject *args, double (*func) (double, double), char *argsfmt)
  69 {
  70         double x, y;
  71         if (! PyArg_ParseTuple(args, argsfmt, &x, &y))
  72                 return NULL;
  73         errno = 0;
  74         PyFPE_START_PROTECT("in math_2", return 0)
  75         x = (*func)(x, y);
  76         PyFPE_END_PROTECT(x)
  77         Py_SET_ERANGE_IF_OVERFLOW(x);
  78         if (errno && is_error(x))
  79                 return NULL;
  80         else
  81                 return PyFloat_FromDouble(x);
  82 }
  83
  84 #define FUNC1(funcname, func, docstring) \
  85         static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
  86                 return math_1(args, func, "d:" #funcname); \
  87         }\
  88         PyDoc_STRVAR(math_##funcname##_doc, docstring);
  89
  90 #define FUNC2(funcname, func, docstring) \
  91         static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
  92                 return math_2(args, func, "dd:" #funcname); \
  93         }\
  94         PyDoc_STRVAR(math_##funcname##_doc, docstring);
  95
  96 FUNC1(acos, acos,
  97       "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
  98 FUNC1(asin, asin,
  99       "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
 100 FUNC1(atan, atan,
 101       "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
 102 FUNC2(atan2, atan2,
 103       "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
 104       "Unlike atan(y/x), the signs of both x and y are considered.")
 105 FUNC1(ceil, ceil,
 106       "ceil(x)\n\nReturn the ceiling of x as a float.\n"
 107       "This is the smallest integral value >= x.")
 108 FUNC1(cos, cos,
 109       "cos(x)\n\nReturn the cosine of x (measured in radians).")
 110 FUNC1(cosh, cosh,
 111       "cosh(x)\n\nReturn the hyperbolic cosine of x.")
 112 FUNC1(exp, exp,
 113       "exp(x)\n\nReturn e raised to the power of x.")
 114 FUNC1(fabs, fabs,
 115       "fabs(x)\n\nReturn the absolute value of the float x.")
 116 FUNC1(floor, floor,
 117       "floor(x)\n\nReturn the floor of x as a float.\n"
 118       "This is the largest integral value <= x.")
 119 FUNC2(fmod, fmod,
 120       "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
 121       "  x % y may differ.")
 122 FUNC2(hypot, hypot,
 123       "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).")
 124 #ifdef MPW_3_1 /* This hack is needed for MPW 3.1 but not for 3.2 ... */
 125 FUNC2(pow, power,
 126       "pow(x,y)\n\nReturn x**y (x to the power of y).")
 127 #else
 128 FUNC2(pow, pow,
 129       "pow(x,y)\n\nReturn x**y (x to the power of y).")
 130 #endif
 131 FUNC1(sin, sin,
 132       "sin(x)\n\nReturn the sine of x (measured in radians).")
 133 FUNC1(sinh, sinh,
 134       "sinh(x)\n\nReturn the hyperbolic sine of x.")
 135 FUNC1(sqrt, sqrt,
 136       "sqrt(x)\n\nReturn the square root of x.")
 137 FUNC1(tan, tan,
 138       "tan(x)\n\nReturn the tangent of x (measured in radians).")
 139 FUNC1(tanh, tanh,
 140       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
 141
 142 static PyObject *
 143 math_frexp(PyObject *self, PyObject *args)
 144 {
 145         double x;
 146         int i;
 147         if (! PyArg_ParseTuple(args, "d:frexp", &x))
 148                 return NULL;
 149         errno = 0;
 150         x = frexp(x, &i);
 151         Py_SET_ERANGE_IF_OVERFLOW(x);
 152         if (errno && is_error(x))
 153                 return NULL;
 154         else
 155                 return Py_BuildValue("(di)", x, i);
 156 }
 157
 158 PyDoc_STRVAR(math_frexp_doc,
 159 "frexp(x)\n"
 160 "\n"
 161 "Return the mantissa and exponent of x, as pair (m, e).\n"
 162 "m is a float and e is an int, such that x = m * 2.**e.\n"
 163 "If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.");
 164
 165 static PyObject *
 166 math_ldexp(PyObject *self, PyObject *args)
 167 {
 168         double x;
 169         int exp;
 170         if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
 171                 return NULL;
 172         errno = 0;
 173         PyFPE_START_PROTECT("ldexp", return 0)
 174         x = ldexp(x, exp);
 175         PyFPE_END_PROTECT(x)
 176         Py_SET_ERANGE_IF_OVERFLOW(x);
 177         if (errno && is_error(x))
 178                 return NULL;
 179         else
 180                 return PyFloat_FromDouble(x);
 181 }
 182
 183 PyDoc_STRVAR(math_ldexp_doc,
 184 "ldexp(x, i) -> x * (2**i)");
 185
 186 static PyObject *
 187 math_modf(PyObject *self, PyObject *args)
 188 {
 189         double x, y;
 190         if (! PyArg_ParseTuple(args, "d:modf", &x))
 191                 return NULL;
 192         errno = 0;
 193 #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
 194         {
 195                 extended e;
 196                 x = modf(x, &e);
 197                 y = e;
 198         }
 199 #else
 200         x = modf(x, &y);
 201 #endif
 202         Py_SET_ERANGE_IF_OVERFLOW(x);
 203         if (errno && is_error(x))
 204                 return NULL;
 205         else
 206                 return Py_BuildValue("(dd)", x, y);
 207 }
 208
 209 PyDoc_STRVAR(math_modf_doc,
 210 "modf(x)\n"
 211 "\n"
 212 "Return the fractional and integer parts of x.  Both results carry the sign\n"
 213 "of x.  The integer part is returned as a real.");
 214
 215 /* A decent logarithm is easy to compute even for huge longs, but libm can't
 216    do that by itself -- loghelper can.  func is log or log10, and name is
 217    "log" or "log10".  Note that overflow isn't possible:  a long can contain
 218    no more than INT_MAX * SHIFT bits, so has value certainly less than
 219    2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
 220    small enough to fit in an IEEE single.  log and log10 are even smaller.
 221 */
 222
 223 static PyObject*
 224 loghelper(PyObject* args, double (*func)(double), char *format, PyObject *arg)
 225 {
 226         /* If it is long, do it ourselves. */
 227         if (PyLong_Check(arg)) {
 228                 double x;
 229                 int e;
 230                 x = _PyLong_AsScaledDouble(arg, &e);
 231                 if (x <= 0.0) {
 232                         PyErr_SetString(PyExc_ValueError,
 233                                         "math domain error");
 234                         return NULL;
 235                 }
 236                 /* Value is ~= x * 2**(e*SHIFT), so the log ~=
 237                    log(x) + log(2) * e * SHIFT.
 238                    CAUTION:  e*SHIFT may overflow using int arithmetic,
 239                    so force use of double. */
 240                 x = func(x) + (e * (double)SHIFT) * func(2.0);
 241                 return PyFloat_FromDouble(x);
 242         }
 243
 244         /* Else let libm handle it by itself. */
 245         return math_1(args, func, format);
 246 }
 247
 248 static PyObject *
 249 math_log(PyObject *self, PyObject *args)
 250 {
 251         PyObject *arg;
 252         PyObject *base = NULL;
 253         PyObject *num, *den;
 254         PyObject *ans;
 255         PyObject *newargs;
 256
 257         if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
 258                 return NULL;
 259         if (base == NULL)
 260                 return loghelper(args, log, "d:log", arg);
 261
 262         newargs = PyTuple_New(1);
 263         if (newargs == NULL)
 264                 return NULL;
 265         Py_INCREF(arg);
 266         PyTuple_SET_ITEM(newargs, 0, arg);
 267         num = loghelper(newargs, log, "d:log", arg);
 268         Py_DECREF(newargs);
 269         if (num == NULL)
 270                 return NULL;
 271
 272         newargs = PyTuple_New(1);
 273         if (newargs == NULL) {
 274                 Py_DECREF(num);
 275                 return NULL;
 276         }
 277         Py_INCREF(base);
 278         PyTuple_SET_ITEM(newargs, 0, base);
 279         den = loghelper(newargs, log, "d:log", base);
 280         Py_DECREF(newargs);
 281         if (den == NULL) {
 282                 Py_DECREF(num);
 283                 return NULL;
 284         }
 285
 286         ans = PyNumber_Divide(num, den);
 287         Py_DECREF(num);
 288         Py_DECREF(den);
 289         return ans;
 290 }
 291
 292 PyDoc_STRVAR(math_log_doc,
 293 "log(x[, base]) -> the logarithm of x to the given base.\n\
 294 If the base not specified, returns the natural logarithm (base e) of x.");
 295
 296 static PyObject *
 297 math_log10(PyObject *self, PyObject *args)
 298 {
 299         PyObject *arg;
 300
 301         if (!PyArg_UnpackTuple(args, "log10", 1, 1, &arg))
 302                 return NULL;
 303         return loghelper(args, log10, "d:log10", arg);
 304 }
 305
 306 PyDoc_STRVAR(math_log10_doc,
 307 "log10(x) -> the base 10 logarithm of x.");
 308
 309 static const double degToRad = 3.141592653589793238462643383 / 180.0;
 310
 311 static PyObject *
 312 math_degrees(PyObject *self, PyObject *args)
 313 {
 314         double x;
 315         if (! PyArg_ParseTuple(args, "d:degrees", &x))
 316                 return NULL;
 317         return PyFloat_FromDouble(x / degToRad);
 318 }
 319
 320 PyDoc_STRVAR(math_degrees_doc,
 321 "degrees(x) -> converts angle x from radians to degrees");
 322
 323 static PyObject *
 324 math_radians(PyObject *self, PyObject *args)
 325 {
 326         double x;
 327         if (! PyArg_ParseTuple(args, "d:radians", &x))
 328                 return NULL;
 329         return PyFloat_FromDouble(x * degToRad);
 330 }
 331
 332 PyDoc_STRVAR(math_radians_doc,
 333 "radians(x) -> converts angle x from degrees to radians");
 334
 335 static PyMethodDef math_methods[] = {
 336         {"acos",        math_acos,      METH_VARARGS,   math_acos_doc},
 337         {"asin",        math_asin,      METH_VARARGS,   math_asin_doc},
 338         {"atan",        math_atan,      METH_VARARGS,   math_atan_doc},
 339         {"atan2",       math_atan2,     METH_VARARGS,   math_atan2_doc},
 340         {"ceil",        math_ceil,      METH_VARARGS,   math_ceil_doc},
 341         {"cos",         math_cos,       METH_VARARGS,   math_cos_doc},
 342         {"cosh",        math_cosh,      METH_VARARGS,   math_cosh_doc},
 343         {"degrees",     math_degrees,   METH_VARARGS,   math_degrees_doc},
 344         {"exp",         math_exp,       METH_VARARGS,   math_exp_doc},
 345         {"fabs",        math_fabs,      METH_VARARGS,   math_fabs_doc},
 346         {"floor",       math_floor,     METH_VARARGS,   math_floor_doc},
 347         {"fmod",        math_fmod,      METH_VARARGS,   math_fmod_doc},
 348         {"frexp",       math_frexp,     METH_VARARGS,   math_frexp_doc},
 349         {"hypot",       math_hypot,     METH_VARARGS,   math_hypot_doc},
 350         {"ldexp",       math_ldexp,     METH_VARARGS,   math_ldexp_doc},
 351         {"log",         math_log,       METH_VARARGS,   math_log_doc},
 352         {"log10",       math_log10,     METH_VARARGS,   math_log10_doc},
 353         {"modf",        math_modf,      METH_VARARGS,   math_modf_doc},
 354         {"pow",         math_pow,       METH_VARARGS,   math_pow_doc},
 355         {"radians",     math_radians,   METH_VARARGS,   math_radians_doc},
 356         {"sin",         math_sin,       METH_VARARGS,   math_sin_doc},
 357         {"sinh",        math_sinh,      METH_VARARGS,   math_sinh_doc},
 358         {"sqrt",        math_sqrt,      METH_VARARGS,   math_sqrt_doc},
 359         {"tan",         math_tan,       METH_VARARGS,   math_tan_doc},
 360         {"tanh",        math_tanh,      METH_VARARGS,   math_tanh_doc},
 361         {NULL,          NULL}           /* sentinel */
 362 };
 363
 364
 365 PyDoc_STRVAR(module_doc,
 366 "This module is always available.  It provides access to the\n"
 367 "mathematical functions defined by the C standard.");
 368
 369 PyMODINIT_FUNC
 370 initmath(void)
 371 {
 372         PyObject *m, *d, *v;
 373
 374         m = Py_InitModule3("math", math_methods, module_doc);
 375         d = PyModule_GetDict(m);
 376
 377         if (!(v = PyFloat_FromDouble(atan(1.0) * 4.0)))
 378                 goto finally;
 379         if (PyDict_SetItemString(d, "pi", v) < 0)
 380                 goto finally;
 381         Py_DECREF(v);
 382
 383         if (!(v = PyFloat_FromDouble(exp(1.0))))
 384                 goto finally;
 385         if (PyDict_SetItemString(d, "e", v) < 0)
 386                 goto finally;
 387         Py_DECREF(v);
 388
 389   finally:
 390         return;
 391 }