blob: 8c483160995a54e390ba9b65be1990b0d97107d9 [file] [log] [blame]
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001/* Math module -- standard C math library functions, pi and e */
2
Christian Heimes53876d92008-04-19 00:31:39 +00003/* Here are some comments from Tim Peters, extracted from the
4 discussion attached to http://bugs.python.org/issue1640. They
5 describe the general aims of the math module with respect to
6 special values, IEEE-754 floating-point exceptions, and Python
7 exceptions.
8
9These are the "spirit of 754" rules:
10
111. If the mathematical result is a real number, but of magnitude too
12large to approximate by a machine float, overflow is signaled and the
13result is an infinity (with the appropriate sign).
14
152. If the mathematical result is a real number, but of magnitude too
16small to approximate by a machine float, underflow is signaled and the
17result is a zero (with the appropriate sign).
18
193. At a singularity (a value x such that the limit of f(y) as y
20approaches x exists and is an infinity), "divide by zero" is signaled
21and the result is an infinity (with the appropriate sign). This is
22complicated a little by that the left-side and right-side limits may
23not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24from the positive or negative directions. In that specific case, the
25sign of the zero determines the result of 1/0.
26
274. At a point where a function has no defined result in the extended
28reals (i.e., the reals plus an infinity or two), invalid operation is
29signaled and a NaN is returned.
30
31And these are what Python has historically /tried/ to do (but not
32always successfully, as platform libm behavior varies a lot):
33
34For #1, raise OverflowError.
35
36For #2, return a zero (with the appropriate sign if that happens by
37accident ;-)).
38
39For #3 and #4, raise ValueError. It may have made sense to raise
40Python's ZeroDivisionError in #3, but historically that's only been
41raised for division by zero and mod by zero.
42
43*/
44
45/*
46 In general, on an IEEE-754 platform the aim is to follow the C99
47 standard, including Annex 'F', whenever possible. Where the
48 standard recommends raising the 'divide-by-zero' or 'invalid'
49 floating-point exceptions, Python should raise a ValueError. Where
50 the standard recommends raising 'overflow', Python should raise an
51 OverflowError. In all other circumstances a value should be
52 returned.
53 */
54
Barry Warsaw8b43b191996-12-09 22:32:36 +000055#include "Python.h"
Michael W. Hudson9ef852c2005-04-06 13:05:18 +000056#include "longintrepr.h" /* just for SHIFT */
Guido van Rossum85a5fbb1990-10-14 12:07:46 +000057
Christian Heimes969fe572008-01-25 11:23:10 +000058#ifdef _OSF_SOURCE
59/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
60extern double copysign(double, double);
61#endif
62
Tim Peters1d120612000-10-12 06:10:25 +000063/* Call is_error when errno != 0, and where x is the result libm
64 * returned. is_error will usually set up an exception and return
65 * true (1), but may return false (0) without setting up an exception.
66 */
67static int
68is_error(double x)
Guido van Rossum8832b621991-12-16 15:44:24 +000069{
Tim Peters1d120612000-10-12 06:10:25 +000070 int result = 1; /* presumption of guilt */
Tim Peters2bf405a2000-10-12 19:42:00 +000071 assert(errno); /* non-zero errno is a precondition for calling */
Guido van Rossum8832b621991-12-16 15:44:24 +000072 if (errno == EDOM)
Barry Warsaw8b43b191996-12-09 22:32:36 +000073 PyErr_SetString(PyExc_ValueError, "math domain error");
Tim Petersa40c7932001-09-05 22:36:56 +000074
Tim Peters1d120612000-10-12 06:10:25 +000075 else if (errno == ERANGE) {
76 /* ANSI C generally requires libm functions to set ERANGE
77 * on overflow, but also generally *allows* them to set
78 * ERANGE on underflow too. There's no consistency about
Tim Petersa40c7932001-09-05 22:36:56 +000079 * the latter across platforms.
80 * Alas, C99 never requires that errno be set.
81 * Here we suppress the underflow errors (libm functions
82 * should return a zero on underflow, and +- HUGE_VAL on
83 * overflow, so testing the result for zero suffices to
84 * distinguish the cases).
Tim Peters1d120612000-10-12 06:10:25 +000085 */
86 if (x)
Tim Petersfe71f812001-08-07 22:10:00 +000087 PyErr_SetString(PyExc_OverflowError,
Tim Peters1d120612000-10-12 06:10:25 +000088 "math range error");
89 else
90 result = 0;
91 }
Guido van Rossum8832b621991-12-16 15:44:24 +000092 else
Barry Warsaw8b43b191996-12-09 22:32:36 +000093 /* Unexpected math error */
94 PyErr_SetFromErrno(PyExc_ValueError);
Tim Peters1d120612000-10-12 06:10:25 +000095 return result;
Guido van Rossum8832b621991-12-16 15:44:24 +000096}
97
Christian Heimes53876d92008-04-19 00:31:39 +000098/*
99 math_1 is used to wrap a libm function f that takes a double
100 arguments and returns a double.
101
102 The error reporting follows these rules, which are designed to do
103 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
104 platforms.
105
106 - a NaN result from non-NaN inputs causes ValueError to be raised
107 - an infinite result from finite inputs causes OverflowError to be
108 raised if can_overflow is 1, or raises ValueError if can_overflow
109 is 0.
110 - if the result is finite and errno == EDOM then ValueError is
111 raised
112 - if the result is finite and nonzero and errno == ERANGE then
113 OverflowError is raised
114
115 The last rule is used to catch overflow on platforms which follow
116 C89 but for which HUGE_VAL is not an infinity.
117
118 For the majority of one-argument functions these rules are enough
119 to ensure that Python's functions behave as specified in 'Annex F'
120 of the C99 standard, with the 'invalid' and 'divide-by-zero'
121 floating-point exceptions mapping to Python's ValueError and the
122 'overflow' floating-point exception mapping to OverflowError.
123 math_1 only works for functions that don't have singularities *and*
124 the possibility of overflow; fortunately, that covers everything we
125 care about right now.
126*/
127
Barry Warsaw8b43b191996-12-09 22:32:36 +0000128static PyObject *
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000129math_1_to_whatever(PyObject *arg, double (*func) (double),
Christian Heimes53876d92008-04-19 00:31:39 +0000130 PyObject *(*from_double_func) (double),
131 int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000132{
Christian Heimes53876d92008-04-19 00:31:39 +0000133 double x, r;
134 x = PyFloat_AsDouble(arg);
Thomas Wouters89f507f2006-12-13 04:49:30 +0000135 if (x == -1.0 && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000136 return NULL;
137 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000138 PyFPE_START_PROTECT("in math_1", return 0);
139 r = (*func)(x);
140 PyFPE_END_PROTECT(r);
141 if (Py_IS_NAN(r)) {
142 if (!Py_IS_NAN(x))
143 errno = EDOM;
144 else
145 errno = 0;
146 }
147 else if (Py_IS_INFINITY(r)) {
148 if (Py_IS_FINITE(x))
149 errno = can_overflow ? ERANGE : EDOM;
150 else
151 errno = 0;
152 }
153 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000154 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000155 else
Christian Heimes53876d92008-04-19 00:31:39 +0000156 return (*from_double_func)(r);
157}
158
159/*
160 math_2 is used to wrap a libm function f that takes two double
161 arguments and returns a double.
162
163 The error reporting follows these rules, which are designed to do
164 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
165 platforms.
166
167 - a NaN result from non-NaN inputs causes ValueError to be raised
168 - an infinite result from finite inputs causes OverflowError to be
169 raised.
170 - if the result is finite and errno == EDOM then ValueError is
171 raised
172 - if the result is finite and nonzero and errno == ERANGE then
173 OverflowError is raised
174
175 The last rule is used to catch overflow on platforms which follow
176 C89 but for which HUGE_VAL is not an infinity.
177
178 For most two-argument functions (copysign, fmod, hypot, atan2)
179 these rules are enough to ensure that Python's functions behave as
180 specified in 'Annex F' of the C99 standard, with the 'invalid' and
181 'divide-by-zero' floating-point exceptions mapping to Python's
182 ValueError and the 'overflow' floating-point exception mapping to
183 OverflowError.
184*/
185
186static PyObject *
187math_1(PyObject *arg, double (*func) (double), int can_overflow)
188{
189 return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000190}
191
192static PyObject *
Christian Heimes53876d92008-04-19 00:31:39 +0000193math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000194{
Christian Heimes53876d92008-04-19 00:31:39 +0000195 return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000196}
197
Barry Warsaw8b43b191996-12-09 22:32:36 +0000198static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000199math_2(PyObject *args, double (*func) (double, double), char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000200{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000201 PyObject *ox, *oy;
Christian Heimes53876d92008-04-19 00:31:39 +0000202 double x, y, r;
Thomas Wouters89f507f2006-12-13 04:49:30 +0000203 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
204 return NULL;
205 x = PyFloat_AsDouble(ox);
206 y = PyFloat_AsDouble(oy);
207 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000208 return NULL;
209 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000210 PyFPE_START_PROTECT("in math_2", return 0);
211 r = (*func)(x, y);
212 PyFPE_END_PROTECT(r);
213 if (Py_IS_NAN(r)) {
214 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
215 errno = EDOM;
216 else
217 errno = 0;
218 }
219 else if (Py_IS_INFINITY(r)) {
220 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
221 errno = ERANGE;
222 else
223 errno = 0;
224 }
225 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000226 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000227 else
Christian Heimes53876d92008-04-19 00:31:39 +0000228 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000229}
230
Christian Heimes53876d92008-04-19 00:31:39 +0000231#define FUNC1(funcname, func, can_overflow, docstring) \
Fred Drake40c48682000-07-03 18:11:56 +0000232 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Christian Heimes53876d92008-04-19 00:31:39 +0000233 return math_1(args, func, can_overflow); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000234 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000235 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000236
Fred Drake40c48682000-07-03 18:11:56 +0000237#define FUNC2(funcname, func, docstring) \
238 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Thomas Wouters89f507f2006-12-13 04:49:30 +0000239 return math_2(args, func, #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000240 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000241 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000242
Christian Heimes53876d92008-04-19 00:31:39 +0000243FUNC1(acos, acos, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000244 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000245FUNC1(acosh, acosh, 0,
246 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
247FUNC1(asin, asin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000248 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000249FUNC1(asinh, asinh, 0,
250 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
251FUNC1(atan, atan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000252 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Martin v. Löwis387c5472001-09-06 08:16:17 +0000253FUNC2(atan2, atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000254 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
255 "Unlike atan(y/x), the signs of both x and y are considered.")
Christian Heimes53876d92008-04-19 00:31:39 +0000256FUNC1(atanh, atanh, 0,
257 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000258
259static PyObject * math_ceil(PyObject *self, PyObject *number) {
260 static PyObject *ceil_str = NULL;
261 PyObject *method;
262
263 if (ceil_str == NULL) {
Christian Heimesfe82e772008-01-28 02:38:20 +0000264 ceil_str = PyUnicode_InternFromString("__ceil__");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000265 if (ceil_str == NULL)
266 return NULL;
267 }
268
Christian Heimes90aa7642007-12-19 02:45:37 +0000269 method = _PyType_Lookup(Py_TYPE(number), ceil_str);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000270 if (method == NULL)
Christian Heimes53876d92008-04-19 00:31:39 +0000271 return math_1_to_int(number, ceil, 0);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000272 else
273 return PyObject_CallFunction(method, "O", number);
274}
275
276PyDoc_STRVAR(math_ceil_doc,
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000277 "ceil(x)\n\nReturn the ceiling of x as an int.\n"
Guido van Rossum13e05de2007-08-23 22:56:55 +0000278 "This is the smallest integral value >= x.");
279
Christian Heimes072c0f12008-01-03 23:01:04 +0000280FUNC2(copysign, copysign,
Christian Heimes53876d92008-04-19 00:31:39 +0000281 "copysign(x,y)\n\nReturn x with the sign of y.")
282FUNC1(cos, cos, 0,
283 "cos(x)\n\nReturn the cosine of x (measured in radians).")
284FUNC1(cosh, cosh, 1,
285 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
286FUNC1(exp, exp, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000287 "exp(x)\n\nReturn e raised to the power of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000288FUNC1(fabs, fabs, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000289 "fabs(x)\n\nReturn the absolute value of the float x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000290
291static PyObject * math_floor(PyObject *self, PyObject *number) {
292 static PyObject *floor_str = NULL;
293 PyObject *method;
294
295 if (floor_str == NULL) {
Christian Heimesfe82e772008-01-28 02:38:20 +0000296 floor_str = PyUnicode_InternFromString("__floor__");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000297 if (floor_str == NULL)
298 return NULL;
299 }
300
Christian Heimes90aa7642007-12-19 02:45:37 +0000301 method = _PyType_Lookup(Py_TYPE(number), floor_str);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000302 if (method == NULL)
Christian Heimes53876d92008-04-19 00:31:39 +0000303 return math_1_to_int(number, floor, 0);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000304 else
305 return PyObject_CallFunction(method, "O", number);
306}
307
308PyDoc_STRVAR(math_floor_doc,
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000309 "floor(x)\n\nReturn the floor of x as an int.\n"
Guido van Rossum13e05de2007-08-23 22:56:55 +0000310 "This is the largest integral value <= x.");
311
Christian Heimes53876d92008-04-19 00:31:39 +0000312FUNC1(log1p, log1p, 1,
313 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
314 The result is computed in a way which is accurate for x near zero.")
315FUNC1(sin, sin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000316 "sin(x)\n\nReturn the sine of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000317FUNC1(sinh, sinh, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000318 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000319FUNC1(sqrt, sqrt, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000320 "sqrt(x)\n\nReturn the square root of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000321FUNC1(tan, tan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000322 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000323FUNC1(tanh, tanh, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000324 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000325
Barry Warsaw8b43b191996-12-09 22:32:36 +0000326static PyObject *
Christian Heimes400adb02008-02-01 08:12:03 +0000327math_trunc(PyObject *self, PyObject *number)
328{
329 static PyObject *trunc_str = NULL;
330 PyObject *trunc;
331
332 if (Py_TYPE(number)->tp_dict == NULL) {
333 if (PyType_Ready(Py_TYPE(number)) < 0)
334 return NULL;
335 }
336
337 if (trunc_str == NULL) {
338 trunc_str = PyUnicode_InternFromString("__trunc__");
339 if (trunc_str == NULL)
340 return NULL;
341 }
342
343 trunc = _PyType_Lookup(Py_TYPE(number), trunc_str);
344 if (trunc == NULL) {
345 PyErr_Format(PyExc_TypeError,
346 "type %.100s doesn't define __trunc__ method",
347 Py_TYPE(number)->tp_name);
348 return NULL;
349 }
350 return PyObject_CallFunctionObjArgs(trunc, number, NULL);
351}
352
353PyDoc_STRVAR(math_trunc_doc,
354"trunc(x:Real) -> Integral\n"
355"\n"
Christian Heimes292d3512008-02-03 16:51:08 +0000356"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
Christian Heimes400adb02008-02-01 08:12:03 +0000357
358static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000359math_frexp(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000360{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000361 int i;
Thomas Wouters89f507f2006-12-13 04:49:30 +0000362 double x = PyFloat_AsDouble(arg);
363 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000364 return NULL;
Christian Heimes53876d92008-04-19 00:31:39 +0000365 /* deal with special cases directly, to sidestep platform
366 differences */
367 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
368 i = 0;
369 }
370 else {
371 PyFPE_START_PROTECT("in math_frexp", return 0);
372 x = frexp(x, &i);
373 PyFPE_END_PROTECT(x);
374 }
375 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000376}
377
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000378PyDoc_STRVAR(math_frexp_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000379"frexp(x)\n"
380"\n"
381"Return the mantissa and exponent of x, as pair (m, e).\n"
382"m is a float and e is an int, such that x = m * 2.**e.\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000383"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 +0000384
Barry Warsaw8b43b191996-12-09 22:32:36 +0000385static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000386math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000387{
Christian Heimes53876d92008-04-19 00:31:39 +0000388 double x, r;
Guido van Rossumc5545052000-05-08 14:29:38 +0000389 int exp;
Fred Drake40c48682000-07-03 18:11:56 +0000390 if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000391 return NULL;
392 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000393 PyFPE_START_PROTECT("in math_ldexp", return 0)
394 r = ldexp(x, exp);
395 PyFPE_END_PROTECT(r)
396 if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
397 errno = ERANGE;
398 /* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
399 we unset it to avoid raising a ValueError here. */
400 if (errno == EDOM)
401 errno = 0;
402 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000403 return NULL;
Guido van Rossumd18ad581991-10-24 14:57:21 +0000404 else
Christian Heimes53876d92008-04-19 00:31:39 +0000405 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000406}
407
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000408PyDoc_STRVAR(math_ldexp_doc,
409"ldexp(x, i) -> x * (2**i)");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000410
Barry Warsaw8b43b191996-12-09 22:32:36 +0000411static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000412math_modf(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000413{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000414 double y, x = PyFloat_AsDouble(arg);
415 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000416 return NULL;
417 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000418 PyFPE_START_PROTECT("in math_modf", return 0);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000419 x = modf(x, &y);
Christian Heimes53876d92008-04-19 00:31:39 +0000420 PyFPE_END_PROTECT(x);
421 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000422}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000423
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000424PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000425"modf(x)\n"
426"\n"
427"Return the fractional and integer parts of x. Both results carry the sign\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000428"of x. The integer part is returned as a real.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000429
Tim Peters78526162001-09-05 00:53:45 +0000430/* A decent logarithm is easy to compute even for huge longs, but libm can't
431 do that by itself -- loghelper can. func is log or log10, and name is
432 "log" or "log10". Note that overflow isn't possible: a long can contain
433 no more than INT_MAX * SHIFT bits, so has value certainly less than
434 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
435 small enough to fit in an IEEE single. log and log10 are even smaller.
436*/
437
438static PyObject*
Thomas Wouters89f507f2006-12-13 04:49:30 +0000439loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +0000440{
Tim Peters78526162001-09-05 00:53:45 +0000441 /* If it is long, do it ourselves. */
442 if (PyLong_Check(arg)) {
443 double x;
444 int e;
445 x = _PyLong_AsScaledDouble(arg, &e);
446 if (x <= 0.0) {
447 PyErr_SetString(PyExc_ValueError,
448 "math domain error");
449 return NULL;
450 }
Christian Heimesaf98da12008-01-27 15:18:18 +0000451 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
452 log(x) + log(2) * e * PyLong_SHIFT.
453 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
Tim Peters78526162001-09-05 00:53:45 +0000454 so force use of double. */
Martin v. Löwis9f2e3462007-07-21 17:22:18 +0000455 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000456 return PyFloat_FromDouble(x);
457 }
458
459 /* Else let libm handle it by itself. */
Christian Heimes53876d92008-04-19 00:31:39 +0000460 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +0000461}
462
463static PyObject *
464math_log(PyObject *self, PyObject *args)
465{
Raymond Hettinger866964c2002-12-14 19:51:34 +0000466 PyObject *arg;
467 PyObject *base = NULL;
468 PyObject *num, *den;
469 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000470
Raymond Hettingerea3fdf42002-12-29 16:33:45 +0000471 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
Raymond Hettinger866964c2002-12-14 19:51:34 +0000472 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000473
Thomas Wouters89f507f2006-12-13 04:49:30 +0000474 num = loghelper(arg, log, "log");
475 if (num == NULL || base == NULL)
476 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000477
Thomas Wouters89f507f2006-12-13 04:49:30 +0000478 den = loghelper(base, log, "log");
Raymond Hettinger866964c2002-12-14 19:51:34 +0000479 if (den == NULL) {
480 Py_DECREF(num);
481 return NULL;
482 }
483
Neal Norwitzbcc0db82006-03-24 08:14:36 +0000484 ans = PyNumber_TrueDivide(num, den);
Raymond Hettinger866964c2002-12-14 19:51:34 +0000485 Py_DECREF(num);
486 Py_DECREF(den);
487 return ans;
Tim Peters78526162001-09-05 00:53:45 +0000488}
489
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000490PyDoc_STRVAR(math_log_doc,
Raymond Hettinger866964c2002-12-14 19:51:34 +0000491"log(x[, base]) -> the logarithm of x to the given base.\n\
492If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +0000493
494static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000495math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +0000496{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000497 return loghelper(arg, log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +0000498}
499
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000500PyDoc_STRVAR(math_log10_doc,
501"log10(x) -> the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +0000502
Christian Heimes53876d92008-04-19 00:31:39 +0000503static PyObject *
504math_fmod(PyObject *self, PyObject *args)
505{
506 PyObject *ox, *oy;
507 double r, x, y;
508 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
509 return NULL;
510 x = PyFloat_AsDouble(ox);
511 y = PyFloat_AsDouble(oy);
512 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
513 return NULL;
514 /* fmod(x, +/-Inf) returns x for finite x. */
515 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
516 return PyFloat_FromDouble(x);
517 errno = 0;
518 PyFPE_START_PROTECT("in math_fmod", return 0);
519 r = fmod(x, y);
520 PyFPE_END_PROTECT(r);
521 if (Py_IS_NAN(r)) {
522 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
523 errno = EDOM;
524 else
525 errno = 0;
526 }
527 if (errno && is_error(r))
528 return NULL;
529 else
530 return PyFloat_FromDouble(r);
531}
532
533PyDoc_STRVAR(math_fmod_doc,
534"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
535" x % y may differ.");
536
537static PyObject *
538math_hypot(PyObject *self, PyObject *args)
539{
540 PyObject *ox, *oy;
541 double r, x, y;
542 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
543 return NULL;
544 x = PyFloat_AsDouble(ox);
545 y = PyFloat_AsDouble(oy);
546 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
547 return NULL;
548 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
549 if (Py_IS_INFINITY(x))
550 return PyFloat_FromDouble(fabs(x));
551 if (Py_IS_INFINITY(y))
552 return PyFloat_FromDouble(fabs(y));
553 errno = 0;
554 PyFPE_START_PROTECT("in math_hypot", return 0);
555 r = hypot(x, y);
556 PyFPE_END_PROTECT(r);
557 if (Py_IS_NAN(r)) {
558 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
559 errno = EDOM;
560 else
561 errno = 0;
562 }
563 else if (Py_IS_INFINITY(r)) {
564 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
565 errno = ERANGE;
566 else
567 errno = 0;
568 }
569 if (errno && is_error(r))
570 return NULL;
571 else
572 return PyFloat_FromDouble(r);
573}
574
575PyDoc_STRVAR(math_hypot_doc,
576"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
577
578/* pow can't use math_2, but needs its own wrapper: the problem is
579 that an infinite result can arise either as a result of overflow
580 (in which case OverflowError should be raised) or as a result of
581 e.g. 0.**-5. (for which ValueError needs to be raised.)
582*/
583
584static PyObject *
585math_pow(PyObject *self, PyObject *args)
586{
587 PyObject *ox, *oy;
588 double r, x, y;
589
590 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
591 return NULL;
592 x = PyFloat_AsDouble(ox);
593 y = PyFloat_AsDouble(oy);
594 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
595 return NULL;
596 /* 1**x and x**0 return 1., even if x is a NaN or infinity. */
597 if (x == 1.0 || y == 0.0)
598 return PyFloat_FromDouble(1.);
599 errno = 0;
600 PyFPE_START_PROTECT("in math_pow", return 0);
601 r = pow(x, y);
602 PyFPE_END_PROTECT(r);
603 if (Py_IS_NAN(r)) {
604 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
605 errno = EDOM;
606 else
607 errno = 0;
608 }
609 /* an infinite result arises either from:
610
611 (A) (+/-0.)**negative,
612 (B) overflow of x**y with both x and y finite (and x nonzero)
613 (C) (+/-inf)**positive, or
614 (D) x**inf with |x| > 1, or x**-inf with |x| < 1.
615
616 In case (A) we want ValueError to be raised. In case (B)
617 OverflowError should be raised. In cases (C) and (D) the infinite
618 result should be returned.
619 */
620 else if (Py_IS_INFINITY(r)) {
621 if (x == 0.)
622 errno = EDOM;
623 else if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
624 errno = ERANGE;
625 else
626 errno = 0;
627 }
628
629 if (errno && is_error(r))
630 return NULL;
631 else
632 return PyFloat_FromDouble(r);
633}
634
635PyDoc_STRVAR(math_pow_doc,
636"pow(x,y)\n\nReturn x**y (x to the power of y).");
637
Christian Heimes072c0f12008-01-03 23:01:04 +0000638static const double degToRad = Py_MATH_PI / 180.0;
639static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000640
641static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000642math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000643{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000644 double x = PyFloat_AsDouble(arg);
645 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000646 return NULL;
Christian Heimes072c0f12008-01-03 23:01:04 +0000647 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000648}
649
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000650PyDoc_STRVAR(math_degrees_doc,
651"degrees(x) -> converts angle x from radians to degrees");
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000652
653static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000654math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000655{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000656 double x = PyFloat_AsDouble(arg);
657 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000658 return NULL;
659 return PyFloat_FromDouble(x * degToRad);
660}
661
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000662PyDoc_STRVAR(math_radians_doc,
663"radians(x) -> converts angle x from degrees to radians");
Tim Peters78526162001-09-05 00:53:45 +0000664
Christian Heimes072c0f12008-01-03 23:01:04 +0000665static PyObject *
666math_isnan(PyObject *self, PyObject *arg)
667{
668 double x = PyFloat_AsDouble(arg);
669 if (x == -1.0 && PyErr_Occurred())
670 return NULL;
671 return PyBool_FromLong((long)Py_IS_NAN(x));
672}
673
674PyDoc_STRVAR(math_isnan_doc,
675"isnan(x) -> bool\n\
676Checks if float x is not a number (NaN)");
677
678static PyObject *
679math_isinf(PyObject *self, PyObject *arg)
680{
681 double x = PyFloat_AsDouble(arg);
682 if (x == -1.0 && PyErr_Occurred())
683 return NULL;
684 return PyBool_FromLong((long)Py_IS_INFINITY(x));
685}
686
687PyDoc_STRVAR(math_isinf_doc,
688"isinf(x) -> bool\n\
689Checks if float x is infinite (positive or negative)");
690
Barry Warsaw8b43b191996-12-09 22:32:36 +0000691static PyMethodDef math_methods[] = {
Thomas Wouters89f507f2006-12-13 04:49:30 +0000692 {"acos", math_acos, METH_O, math_acos_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000693 {"acosh", math_acosh, METH_O, math_acosh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000694 {"asin", math_asin, METH_O, math_asin_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000695 {"asinh", math_asinh, METH_O, math_asinh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000696 {"atan", math_atan, METH_O, math_atan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000697 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000698 {"atanh", math_atanh, METH_O, math_atanh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000699 {"ceil", math_ceil, METH_O, math_ceil_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000700 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000701 {"cos", math_cos, METH_O, math_cos_doc},
702 {"cosh", math_cosh, METH_O, math_cosh_doc},
703 {"degrees", math_degrees, METH_O, math_degrees_doc},
704 {"exp", math_exp, METH_O, math_exp_doc},
705 {"fabs", math_fabs, METH_O, math_fabs_doc},
706 {"floor", math_floor, METH_O, math_floor_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000707 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000708 {"frexp", math_frexp, METH_O, math_frexp_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000709 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000710 {"isinf", math_isinf, METH_O, math_isinf_doc},
711 {"isnan", math_isnan, METH_O, math_isnan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000712 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
713 {"log", math_log, METH_VARARGS, math_log_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000714 {"log1p", math_log1p, METH_O, math_log1p_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000715 {"log10", math_log10, METH_O, math_log10_doc},
716 {"modf", math_modf, METH_O, math_modf_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000717 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000718 {"radians", math_radians, METH_O, math_radians_doc},
719 {"sin", math_sin, METH_O, math_sin_doc},
720 {"sinh", math_sinh, METH_O, math_sinh_doc},
721 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
722 {"tan", math_tan, METH_O, math_tan_doc},
723 {"tanh", math_tanh, METH_O, math_tanh_doc},
Christian Heimes400adb02008-02-01 08:12:03 +0000724 {"trunc", math_trunc, METH_O, math_trunc_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000725 {NULL, NULL} /* sentinel */
726};
727
Guido van Rossumc6e22901998-12-04 19:26:43 +0000728
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000729PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000730"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000731"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000732
Mark Hammondfe51c6d2002-08-02 02:27:13 +0000733PyMODINIT_FUNC
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000734initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000735{
Christian Heimes53876d92008-04-19 00:31:39 +0000736 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +0000737
Guido van Rossumc6e22901998-12-04 19:26:43 +0000738 m = Py_InitModule3("math", math_methods, module_doc);
Neal Norwitz1ac754f2006-01-19 06:09:39 +0000739 if (m == NULL)
740 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +0000741
Christian Heimes53876d92008-04-19 00:31:39 +0000742 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
743 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +0000744
Christian Heimes53876d92008-04-19 00:31:39 +0000745 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000746 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000747}