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 *name)
 225 {
 226         PyObject *arg;
 227         char format[16];
 228
 229         /* See whether this is a long. */
 230         format[0] = 'O';
 231         format[1] = ':';
 232         strcpy(format + 2, name);
 233         if (! PyArg_ParseTuple(args, format, &arg))
 234                 return NULL;
 235
 236         /* If it is long, do it ourselves. */
 237         if (PyLong_Check(arg)) {
 238                 double x;
 239                 int e;
 240                 x = _PyLong_AsScaledDouble(arg, &e);
 241                 if (x <= 0.0) {
 242                         PyErr_SetString(PyExc_ValueError,
 243                                         "math domain error");
 244                         return NULL;
 245                 }
 246                 /* Value is ~= x * 2**(e*SHIFT), so the log ~=
 247                    log(x) + log(2) * e * SHIFT.
 248                    CAUTION:  e*SHIFT may overflow using int arithmetic,
 249                    so force use of double. */
 250                 x = func(x) + (e * (double)SHIFT) * func(2.0);
 251                 return PyFloat_FromDouble(x);
 252         }
 253
 254         /* Else let libm handle it by itself. */
 255         format[0] = 'd';
 256         return math_1(args, func, format);
 257 }
 258
 259 static PyObject *
 260 math_log(PyObject *self, PyObject *args)
 261 {
 262         return loghelper(args, log, "log");
 263 }
 264
 265 PyDoc_STRVAR(math_log_doc,
 266 "log(x) -> the natural logarithm (base e) of x.");
 267
 268 static PyObject *
 269 math_log10(PyObject *self, PyObject *args)
 270 {
 271         return loghelper(args, log10, "log10");
 272 }
 273
 274 PyDoc_STRVAR(math_log10_doc,
 275 "log10(x) -> the base 10 logarithm of x.");
 276
 277 static const double degToRad = 3.141592653589793238462643383 / 180.0;
 278
 279 static PyObject *
 280 math_degrees(PyObject *self, PyObject *args)
 281 {
 282         double x;
 283         if (! PyArg_ParseTuple(args, "d:degrees", &x))
 284                 return NULL;
 285         return PyFloat_FromDouble(x / degToRad);
 286 }
 287
 288 PyDoc_STRVAR(math_degrees_doc,
 289 "degrees(x) -> converts angle x from radians to degrees");
 290
 291 static PyObject *
 292 math_radians(PyObject *self, PyObject *args)
 293 {
 294         double x;
 295         if (! PyArg_ParseTuple(args, "d:radians", &x))
 296                 return NULL;
 297         return PyFloat_FromDouble(x * degToRad);
 298 }
 299
 300 PyDoc_STRVAR(math_radians_doc,
 301 "radians(x) -> converts angle x from degrees to radians");
 302
 303 static PyMethodDef math_methods[] = {
 304         {"acos",        math_acos,      METH_VARARGS,   math_acos_doc},
 305         {"asin",        math_asin,      METH_VARARGS,   math_asin_doc},
 306         {"atan",        math_atan,      METH_VARARGS,   math_atan_doc},
 307         {"atan2",       math_atan2,     METH_VARARGS,   math_atan2_doc},
 308         {"ceil",        math_ceil,      METH_VARARGS,   math_ceil_doc},
 309         {"cos",         math_cos,       METH_VARARGS,   math_cos_doc},
 310         {"cosh",        math_cosh,      METH_VARARGS,   math_cosh_doc},
 311         {"degrees",     math_degrees,   METH_VARARGS,   math_degrees_doc},
 312         {"exp",         math_exp,       METH_VARARGS,   math_exp_doc},
 313         {"fabs",        math_fabs,      METH_VARARGS,   math_fabs_doc},
 314         {"floor",       math_floor,     METH_VARARGS,   math_floor_doc},
 315         {"fmod",        math_fmod,      METH_VARARGS,   math_fmod_doc},
 316         {"frexp",       math_frexp,     METH_VARARGS,   math_frexp_doc},
 317         {"hypot",       math_hypot,     METH_VARARGS,   math_hypot_doc},
 318         {"ldexp",       math_ldexp,     METH_VARARGS,   math_ldexp_doc},
 319         {"log",         math_log,       METH_VARARGS,   math_log_doc},
 320         {"log10",       math_log10,     METH_VARARGS,   math_log10_doc},
 321         {"modf",        math_modf,      METH_VARARGS,   math_modf_doc},
 322         {"pow",         math_pow,       METH_VARARGS,   math_pow_doc},
 323         {"radians",     math_radians,   METH_VARARGS,   math_radians_doc},
 324         {"sin",         math_sin,       METH_VARARGS,   math_sin_doc},
 325         {"sinh",        math_sinh,      METH_VARARGS,   math_sinh_doc},
 326         {"sqrt",        math_sqrt,      METH_VARARGS,   math_sqrt_doc},
 327         {"tan",         math_tan,       METH_VARARGS,   math_tan_doc},
 328         {"tanh",        math_tanh,      METH_VARARGS,   math_tanh_doc},
 329         {NULL,          NULL}           /* sentinel */
 330 };
 331
 332
 333 PyDoc_STRVAR(module_doc,
 334 "This module is always available.  It provides access to the\n"
 335 "mathematical functions defined by the C standard.");
 336
 337 PyMODINIT_FUNC
 338 initmath(void)
 339 {
 340         PyObject *m, *d, *v;
 341
 342         m = Py_InitModule3("math", math_methods, module_doc);
 343         d = PyModule_GetDict(m);
 344
 345         if (!(v = PyFloat_FromDouble(atan(1.0) * 4.0)))
 346                 goto finally;
 347         if (PyDict_SetItemString(d, "pi", v) < 0)
 348                 goto finally;
 349         Py_DECREF(v);
 350
 351         if (!(v = PyFloat_FromDouble(exp(1.0))))
 352                 goto finally;
 353         if (PyDict_SetItemString(d, "e", v) < 0)
 354                 goto finally;
 355         Py_DECREF(v);
 356
 357   finally:
 358         return;
 359 }