blob: 4609f6040f448f5ecba8e3144c33d8d1e0b8656c [file] [log] [blame]
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001/* Math module -- standard C math library functions, pi and e */
2
Barry Warsaw8b43b191996-12-09 22:32:36 +00003#include "Python.h"
Tim Peters78526162001-09-05 00:53:45 +00004#include "longintrepr.h"
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00005
Guido van Rossum7081cf51996-05-23 22:56:19 +00006#ifndef _MSC_VER
Guido van Rossumb9418681995-01-12 11:28:16 +00007#ifndef __STDC__
Fred Drake40c48682000-07-03 18:11:56 +00008extern double fmod (double, double);
9extern double frexp (double, int *);
10extern double ldexp (double, int);
11extern double modf (double, double *);
Guido van Rossum7081cf51996-05-23 22:56:19 +000012#endif /* __STDC__ */
13#endif /* _MSC_VER */
14
Tim Peters1d120612000-10-12 06:10:25 +000015/* 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 */
19static int
20is_error(double x)
Guido van Rossum8832b621991-12-16 15:44:24 +000021{
Tim Peters1d120612000-10-12 06:10:25 +000022 int result = 1; /* presumption of guilt */
Tim Peters2bf405a2000-10-12 19:42:00 +000023 assert(errno); /* non-zero errno is a precondition for calling */
Guido van Rossum8832b621991-12-16 15:44:24 +000024 if (errno == EDOM)
Barry Warsaw8b43b191996-12-09 22:32:36 +000025 PyErr_SetString(PyExc_ValueError, "math domain error");
Tim Petersa40c7932001-09-05 22:36:56 +000026
Tim Peters1d120612000-10-12 06:10:25 +000027 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
Tim Petersa40c7932001-09-05 22:36:56 +000031 * 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).
Tim Peters1d120612000-10-12 06:10:25 +000037 */
38 if (x)
Tim Petersfe71f812001-08-07 22:10:00 +000039 PyErr_SetString(PyExc_OverflowError,
Tim Peters1d120612000-10-12 06:10:25 +000040 "math range error");
41 else
42 result = 0;
43 }
Guido van Rossum8832b621991-12-16 15:44:24 +000044 else
Barry Warsaw8b43b191996-12-09 22:32:36 +000045 /* Unexpected math error */
46 PyErr_SetFromErrno(PyExc_ValueError);
Tim Peters1d120612000-10-12 06:10:25 +000047 return result;
Guido van Rossum8832b621991-12-16 15:44:24 +000048}
49
Barry Warsaw8b43b191996-12-09 22:32:36 +000050static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +000051math_1(PyObject *args, double (*func) (double), char *argsfmt)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000052{
53 double x;
Fred Drake40c48682000-07-03 18:11:56 +000054 if (! PyArg_ParseTuple(args, argsfmt, &x))
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000055 return NULL;
56 errno = 0;
Guido van Rossum52fa3a61997-02-14 22:59:58 +000057 PyFPE_START_PROTECT("in math_1", return 0)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000058 x = (*func)(x);
Guido van Rossum45b83911997-03-14 04:32:50 +000059 PyFPE_END_PROTECT(x)
Tim Petersa40c7932001-09-05 22:36:56 +000060 Py_SET_ERANGE_IF_OVERFLOW(x);
Tim Peters1d120612000-10-12 06:10:25 +000061 if (errno && is_error(x))
62 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000063 else
Barry Warsaw8b43b191996-12-09 22:32:36 +000064 return PyFloat_FromDouble(x);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000065}
66
Barry Warsaw8b43b191996-12-09 22:32:36 +000067static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +000068math_2(PyObject *args, double (*func) (double, double), char *argsfmt)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000069{
70 double x, y;
Fred Drake40c48682000-07-03 18:11:56 +000071 if (! PyArg_ParseTuple(args, argsfmt, &x, &y))
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000072 return NULL;
73 errno = 0;
Guido van Rossum52fa3a61997-02-14 22:59:58 +000074 PyFPE_START_PROTECT("in math_2", return 0)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000075 x = (*func)(x, y);
Guido van Rossum45b83911997-03-14 04:32:50 +000076 PyFPE_END_PROTECT(x)
Tim Petersa40c7932001-09-05 22:36:56 +000077 Py_SET_ERANGE_IF_OVERFLOW(x);
Tim Peters1d120612000-10-12 06:10:25 +000078 if (errno && is_error(x))
79 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000080 else
Barry Warsaw8b43b191996-12-09 22:32:36 +000081 return PyFloat_FromDouble(x);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000082}
83
Fred Drake40c48682000-07-03 18:11:56 +000084#define FUNC1(funcname, func, docstring) \
85 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
86 return math_1(args, func, "d:" #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +000087 }\
Fred Drake40c48682000-07-03 18:11:56 +000088 static char math_##funcname##_doc [] = docstring;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000089
Fred Drake40c48682000-07-03 18:11:56 +000090#define FUNC2(funcname, func, docstring) \
91 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
92 return math_2(args, func, "dd:" #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +000093 }\
Fred Drake40c48682000-07-03 18:11:56 +000094 static char math_##funcname##_doc [] = docstring;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000095
Fred Drake40c48682000-07-03 18:11:56 +000096FUNC1(acos, acos,
Tim Petersfe71f812001-08-07 22:10:00 +000097 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Fred Drake40c48682000-07-03 18:11:56 +000098FUNC1(asin, asin,
Tim Petersfe71f812001-08-07 22:10:00 +000099 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Fred Drake40c48682000-07-03 18:11:56 +0000100FUNC1(atan, atan,
Tim Petersfe71f812001-08-07 22:10:00 +0000101 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Martin v. Löwis387c5472001-09-06 08:16:17 +0000102FUNC2(atan2, atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000103 "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.")
Fred Drake40c48682000-07-03 18:11:56 +0000105FUNC1(ceil, ceil,
Tim Petersfe71f812001-08-07 22:10:00 +0000106 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
107 "This is the smallest integral value >= x.")
Fred Drake40c48682000-07-03 18:11:56 +0000108FUNC1(cos, cos,
Tim Petersfe71f812001-08-07 22:10:00 +0000109 "cos(x)\n\nReturn the cosine of x (measured in radians).")
Fred Drake40c48682000-07-03 18:11:56 +0000110FUNC1(cosh, cosh,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000111 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
Fred Drake40c48682000-07-03 18:11:56 +0000112FUNC1(exp, exp,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000113 "exp(x)\n\nReturn e raised to the power of x.")
Fred Drake40c48682000-07-03 18:11:56 +0000114FUNC1(fabs, fabs,
Tim Petersfe71f812001-08-07 22:10:00 +0000115 "fabs(x)\n\nReturn the absolute value of the float x.")
Fred Drake40c48682000-07-03 18:11:56 +0000116FUNC1(floor, floor,
Tim Petersfe71f812001-08-07 22:10:00 +0000117 "floor(x)\n\nReturn the floor of x as a float.\n"
118 "This is the largest integral value <= x.")
Peter Schneider-Kamp07c19222001-06-10 16:45:08 +0000119FUNC2(fmod, fmod,
Tim Peters78fc0b52000-09-16 03:54:24 +0000120 "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
121 " x % y may differ.")
Fred Drake40c48682000-07-03 18:11:56 +0000122FUNC2(hypot, hypot,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000123 "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).")
Guido van Rossum1492c271991-07-27 21:38:43 +0000124#ifdef MPW_3_1 /* This hack is needed for MPW 3.1 but not for 3.2 ... */
Fred Drake40c48682000-07-03 18:11:56 +0000125FUNC2(pow, power,
Tim Petersfe71f812001-08-07 22:10:00 +0000126 "pow(x,y)\n\nReturn x**y (x to the power of y).")
Guido van Rossum76f2f2e1991-06-24 22:23:10 +0000127#else
Fred Drake40c48682000-07-03 18:11:56 +0000128FUNC2(pow, pow,
Tim Petersfe71f812001-08-07 22:10:00 +0000129 "pow(x,y)\n\nReturn x**y (x to the power of y).")
Guido van Rossum76f2f2e1991-06-24 22:23:10 +0000130#endif
Fred Drake40c48682000-07-03 18:11:56 +0000131FUNC1(sin, sin,
Tim Petersfe71f812001-08-07 22:10:00 +0000132 "sin(x)\n\nReturn the sine of x (measured in radians).")
Fred Drake40c48682000-07-03 18:11:56 +0000133FUNC1(sinh, sinh,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000134 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Fred Drake40c48682000-07-03 18:11:56 +0000135FUNC1(sqrt, sqrt,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000136 "sqrt(x)\n\nReturn the square root of x.")
Fred Drake40c48682000-07-03 18:11:56 +0000137FUNC1(tan, tan,
Tim Petersfe71f812001-08-07 22:10:00 +0000138 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Fred Drake40c48682000-07-03 18:11:56 +0000139FUNC1(tanh, tanh,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000140 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000141
Barry Warsaw8b43b191996-12-09 22:32:36 +0000142static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000143math_frexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000144{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000145 double x;
146 int i;
Fred Drake40c48682000-07-03 18:11:56 +0000147 if (! PyArg_ParseTuple(args, "d:frexp", &x))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000148 return NULL;
149 errno = 0;
150 x = frexp(x, &i);
Tim Petersa40c7932001-09-05 22:36:56 +0000151 Py_SET_ERANGE_IF_OVERFLOW(x);
Tim Peters1d120612000-10-12 06:10:25 +0000152 if (errno && is_error(x))
153 return NULL;
154 else
155 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000156}
157
Guido van Rossumc6e22901998-12-04 19:26:43 +0000158static char math_frexp_doc [] =
Tim Peters63c94532001-09-04 23:17:42 +0000159"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.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000164
Barry Warsaw8b43b191996-12-09 22:32:36 +0000165static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000166math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000167{
Guido van Rossumc5545052000-05-08 14:29:38 +0000168 double x;
169 int exp;
Fred Drake40c48682000-07-03 18:11:56 +0000170 if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000171 return NULL;
172 errno = 0;
Guido van Rossum52fa3a61997-02-14 22:59:58 +0000173 PyFPE_START_PROTECT("ldexp", return 0)
Guido van Rossumc5545052000-05-08 14:29:38 +0000174 x = ldexp(x, exp);
Guido van Rossum45b83911997-03-14 04:32:50 +0000175 PyFPE_END_PROTECT(x)
Tim Petersa40c7932001-09-05 22:36:56 +0000176 Py_SET_ERANGE_IF_OVERFLOW(x);
Tim Peters1d120612000-10-12 06:10:25 +0000177 if (errno && is_error(x))
178 return NULL;
Guido van Rossumd18ad581991-10-24 14:57:21 +0000179 else
Barry Warsaw8b43b191996-12-09 22:32:36 +0000180 return PyFloat_FromDouble(x);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000181}
182
Tim Petersfe71f812001-08-07 22:10:00 +0000183static char math_ldexp_doc [] =
Tim Peters63c94532001-09-04 23:17:42 +0000184"ldexp(x, i) -> x * (2**i)";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000185
Barry Warsaw8b43b191996-12-09 22:32:36 +0000186static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000187math_modf(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000188{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000189 double x, y;
Fred Drake40c48682000-07-03 18:11:56 +0000190 if (! PyArg_ParseTuple(args, "d:modf", &x))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000191 return NULL;
192 errno = 0;
Guido van Rossumb6775db1994-08-01 11:34:53 +0000193#ifdef MPW /* MPW C modf expects pointer to extended as second argument */
Fred Drake40c48682000-07-03 18:11:56 +0000194 {
195 extended e;
196 x = modf(x, &e);
197 y = e;
198 }
Guido van Rossumb6775db1994-08-01 11:34:53 +0000199#else
Guido van Rossumd18ad581991-10-24 14:57:21 +0000200 x = modf(x, &y);
Guido van Rossumb6775db1994-08-01 11:34:53 +0000201#endif
Tim Petersa40c7932001-09-05 22:36:56 +0000202 Py_SET_ERANGE_IF_OVERFLOW(x);
Tim Peters1d120612000-10-12 06:10:25 +0000203 if (errno && is_error(x))
204 return NULL;
205 else
206 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000207}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000208
Guido van Rossumc6e22901998-12-04 19:26:43 +0000209static char math_modf_doc [] =
Tim Peters63c94532001-09-04 23:17:42 +0000210"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.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000214
Tim Peters78526162001-09-05 00:53:45 +0000215/* 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
223static PyObject*
224loghelper(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. */
Tim Peterse5ca6c72001-09-05 04:33:11 +0000250 x = func(x) + (e * (double)SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000251 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
259static PyObject *
260math_log(PyObject *self, PyObject *args)
261{
262 return loghelper(args, log, "log");
263}
264
265static char math_log_doc[] =
266"log(x) -> the natural logarithm (base e) of x.";
267
268static PyObject *
269math_log10(PyObject *self, PyObject *args)
270{
271 return loghelper(args, log10, "log10");
272}
273
274static char math_log10_doc[] =
275"log10(x) -> the base 10 logarithm of x.";
276
277
Barry Warsaw8b43b191996-12-09 22:32:36 +0000278static PyMethodDef math_methods[] = {
Fred Drake40c48682000-07-03 18:11:56 +0000279 {"acos", math_acos, METH_VARARGS, math_acos_doc},
280 {"asin", math_asin, METH_VARARGS, math_asin_doc},
281 {"atan", math_atan, METH_VARARGS, math_atan_doc},
282 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
283 {"ceil", math_ceil, METH_VARARGS, math_ceil_doc},
284 {"cos", math_cos, METH_VARARGS, math_cos_doc},
285 {"cosh", math_cosh, METH_VARARGS, math_cosh_doc},
286 {"exp", math_exp, METH_VARARGS, math_exp_doc},
287 {"fabs", math_fabs, METH_VARARGS, math_fabs_doc},
288 {"floor", math_floor, METH_VARARGS, math_floor_doc},
289 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
290 {"frexp", math_frexp, METH_VARARGS, math_frexp_doc},
291 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
292 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
293 {"log", math_log, METH_VARARGS, math_log_doc},
294 {"log10", math_log10, METH_VARARGS, math_log10_doc},
295 {"modf", math_modf, METH_VARARGS, math_modf_doc},
296 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000297 {"sin", math_sin, METH_VARARGS, math_sin_doc},
298 {"sinh", math_sinh, METH_VARARGS, math_sinh_doc},
299 {"sqrt", math_sqrt, METH_VARARGS, math_sqrt_doc},
300 {"tan", math_tan, METH_VARARGS, math_tan_doc},
301 {"tanh", math_tanh, METH_VARARGS, math_tanh_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000302 {NULL, NULL} /* sentinel */
303};
304
Guido van Rossumc6e22901998-12-04 19:26:43 +0000305
306static char module_doc [] =
Tim Peters63c94532001-09-04 23:17:42 +0000307"This module is always available. It provides access to the\n"
308"mathematical functions defined by the C standard.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000309
Guido van Rossum3886bb61998-12-04 18:50:17 +0000310DL_EXPORT(void)
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000311initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000312{
Barry Warsaw8b43b191996-12-09 22:32:36 +0000313 PyObject *m, *d, *v;
Tim Petersfe71f812001-08-07 22:10:00 +0000314
Guido van Rossumc6e22901998-12-04 19:26:43 +0000315 m = Py_InitModule3("math", math_methods, module_doc);
Barry Warsaw8b43b191996-12-09 22:32:36 +0000316 d = PyModule_GetDict(m);
Barry Warsawfc93f751996-12-17 00:47:03 +0000317
318 if (!(v = PyFloat_FromDouble(atan(1.0) * 4.0)))
319 goto finally;
320 if (PyDict_SetItemString(d, "pi", v) < 0)
321 goto finally;
Barry Warsaw8b43b191996-12-09 22:32:36 +0000322 Py_DECREF(v);
Barry Warsawfc93f751996-12-17 00:47:03 +0000323
324 if (!(v = PyFloat_FromDouble(exp(1.0))))
325 goto finally;
Guido van Rossum4c4cbf31996-12-18 14:12:22 +0000326 if (PyDict_SetItemString(d, "e", v) < 0)
Barry Warsawfc93f751996-12-17 00:47:03 +0000327 goto finally;
Barry Warsaw8b43b191996-12-09 22:32:36 +0000328 Py_DECREF(v);
Barry Warsawfc93f751996-12-17 00:47:03 +0000329
330 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000331 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000332}