blob: 3cf113383116397f84d015d24b76ca9dd7dfff33 [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/*
Christian Heimese57950f2008-04-21 13:08:03 +000099 wrapper for atan2 that deals directly with special cases before
100 delegating to the platform libm for the remaining cases. This
101 is necessary to get consistent behaviour across platforms.
102 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
103 always follow C99.
104*/
105
106static double
107m_atan2(double y, double x)
108{
109 if (Py_IS_NAN(x) || Py_IS_NAN(y))
110 return Py_NAN;
111 if (Py_IS_INFINITY(y)) {
112 if (Py_IS_INFINITY(x)) {
113 if (copysign(1., x) == 1.)
114 /* atan2(+-inf, +inf) == +-pi/4 */
115 return copysign(0.25*Py_MATH_PI, y);
116 else
117 /* atan2(+-inf, -inf) == +-pi*3/4 */
118 return copysign(0.75*Py_MATH_PI, y);
119 }
120 /* atan2(+-inf, x) == +-pi/2 for finite x */
121 return copysign(0.5*Py_MATH_PI, y);
122 }
123 if (Py_IS_INFINITY(x) || y == 0.) {
124 if (copysign(1., x) == 1.)
125 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
126 return copysign(0., y);
127 else
128 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
129 return copysign(Py_MATH_PI, y);
130 }
131 return atan2(y, x);
132}
133
134/*
Christian Heimes53876d92008-04-19 00:31:39 +0000135 math_1 is used to wrap a libm function f that takes a double
136 arguments and returns a double.
137
138 The error reporting follows these rules, which are designed to do
139 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
140 platforms.
141
142 - a NaN result from non-NaN inputs causes ValueError to be raised
143 - an infinite result from finite inputs causes OverflowError to be
144 raised if can_overflow is 1, or raises ValueError if can_overflow
145 is 0.
146 - if the result is finite and errno == EDOM then ValueError is
147 raised
148 - if the result is finite and nonzero and errno == ERANGE then
149 OverflowError is raised
150
151 The last rule is used to catch overflow on platforms which follow
152 C89 but for which HUGE_VAL is not an infinity.
153
154 For the majority of one-argument functions these rules are enough
155 to ensure that Python's functions behave as specified in 'Annex F'
156 of the C99 standard, with the 'invalid' and 'divide-by-zero'
157 floating-point exceptions mapping to Python's ValueError and the
158 'overflow' floating-point exception mapping to OverflowError.
159 math_1 only works for functions that don't have singularities *and*
160 the possibility of overflow; fortunately, that covers everything we
161 care about right now.
162*/
163
Barry Warsaw8b43b191996-12-09 22:32:36 +0000164static PyObject *
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000165math_1_to_whatever(PyObject *arg, double (*func) (double),
Christian Heimes53876d92008-04-19 00:31:39 +0000166 PyObject *(*from_double_func) (double),
167 int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000168{
Christian Heimes53876d92008-04-19 00:31:39 +0000169 double x, r;
170 x = PyFloat_AsDouble(arg);
Thomas Wouters89f507f2006-12-13 04:49:30 +0000171 if (x == -1.0 && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000172 return NULL;
173 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000174 PyFPE_START_PROTECT("in math_1", return 0);
175 r = (*func)(x);
176 PyFPE_END_PROTECT(r);
177 if (Py_IS_NAN(r)) {
178 if (!Py_IS_NAN(x))
179 errno = EDOM;
180 else
181 errno = 0;
182 }
183 else if (Py_IS_INFINITY(r)) {
184 if (Py_IS_FINITE(x))
185 errno = can_overflow ? ERANGE : EDOM;
186 else
187 errno = 0;
188 }
189 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000190 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000191 else
Christian Heimes53876d92008-04-19 00:31:39 +0000192 return (*from_double_func)(r);
193}
194
195/*
196 math_2 is used to wrap a libm function f that takes two double
197 arguments and returns a double.
198
199 The error reporting follows these rules, which are designed to do
200 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
201 platforms.
202
203 - a NaN result from non-NaN inputs causes ValueError to be raised
204 - an infinite result from finite inputs causes OverflowError to be
205 raised.
206 - if the result is finite and errno == EDOM then ValueError is
207 raised
208 - if the result is finite and nonzero and errno == ERANGE then
209 OverflowError is raised
210
211 The last rule is used to catch overflow on platforms which follow
212 C89 but for which HUGE_VAL is not an infinity.
213
214 For most two-argument functions (copysign, fmod, hypot, atan2)
215 these rules are enough to ensure that Python's functions behave as
216 specified in 'Annex F' of the C99 standard, with the 'invalid' and
217 'divide-by-zero' floating-point exceptions mapping to Python's
218 ValueError and the 'overflow' floating-point exception mapping to
219 OverflowError.
220*/
221
222static PyObject *
223math_1(PyObject *arg, double (*func) (double), int can_overflow)
224{
225 return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000226}
227
228static PyObject *
Christian Heimes53876d92008-04-19 00:31:39 +0000229math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000230{
Christian Heimes53876d92008-04-19 00:31:39 +0000231 return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000232}
233
Barry Warsaw8b43b191996-12-09 22:32:36 +0000234static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000235math_2(PyObject *args, double (*func) (double, double), char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000236{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000237 PyObject *ox, *oy;
Christian Heimes53876d92008-04-19 00:31:39 +0000238 double x, y, r;
Thomas Wouters89f507f2006-12-13 04:49:30 +0000239 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
240 return NULL;
241 x = PyFloat_AsDouble(ox);
242 y = PyFloat_AsDouble(oy);
243 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000244 return NULL;
245 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000246 PyFPE_START_PROTECT("in math_2", return 0);
247 r = (*func)(x, y);
248 PyFPE_END_PROTECT(r);
249 if (Py_IS_NAN(r)) {
250 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
251 errno = EDOM;
252 else
253 errno = 0;
254 }
255 else if (Py_IS_INFINITY(r)) {
256 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
257 errno = ERANGE;
258 else
259 errno = 0;
260 }
261 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000262 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000263 else
Christian Heimes53876d92008-04-19 00:31:39 +0000264 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000265}
266
Christian Heimes53876d92008-04-19 00:31:39 +0000267#define FUNC1(funcname, func, can_overflow, docstring) \
Fred Drake40c48682000-07-03 18:11:56 +0000268 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Christian Heimes53876d92008-04-19 00:31:39 +0000269 return math_1(args, func, can_overflow); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000270 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000271 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000272
Fred Drake40c48682000-07-03 18:11:56 +0000273#define FUNC2(funcname, func, docstring) \
274 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Thomas Wouters89f507f2006-12-13 04:49:30 +0000275 return math_2(args, func, #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000276 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000277 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000278
Christian Heimes53876d92008-04-19 00:31:39 +0000279FUNC1(acos, acos, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000280 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000281FUNC1(acosh, acosh, 0,
282 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
283FUNC1(asin, asin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000284 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000285FUNC1(asinh, asinh, 0,
286 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
287FUNC1(atan, atan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000288 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Christian Heimese57950f2008-04-21 13:08:03 +0000289FUNC2(atan2, m_atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000290 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
291 "Unlike atan(y/x), the signs of both x and y are considered.")
Christian Heimes53876d92008-04-19 00:31:39 +0000292FUNC1(atanh, atanh, 0,
293 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000294
295static PyObject * math_ceil(PyObject *self, PyObject *number) {
296 static PyObject *ceil_str = NULL;
297 PyObject *method;
298
299 if (ceil_str == NULL) {
Christian Heimesfe82e772008-01-28 02:38:20 +0000300 ceil_str = PyUnicode_InternFromString("__ceil__");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000301 if (ceil_str == NULL)
302 return NULL;
303 }
304
Christian Heimes90aa7642007-12-19 02:45:37 +0000305 method = _PyType_Lookup(Py_TYPE(number), ceil_str);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000306 if (method == NULL)
Christian Heimes53876d92008-04-19 00:31:39 +0000307 return math_1_to_int(number, ceil, 0);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000308 else
309 return PyObject_CallFunction(method, "O", number);
310}
311
312PyDoc_STRVAR(math_ceil_doc,
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000313 "ceil(x)\n\nReturn the ceiling of x as an int.\n"
Guido van Rossum13e05de2007-08-23 22:56:55 +0000314 "This is the smallest integral value >= x.");
315
Christian Heimes072c0f12008-01-03 23:01:04 +0000316FUNC2(copysign, copysign,
Christian Heimes53876d92008-04-19 00:31:39 +0000317 "copysign(x,y)\n\nReturn x with the sign of y.")
318FUNC1(cos, cos, 0,
319 "cos(x)\n\nReturn the cosine of x (measured in radians).")
320FUNC1(cosh, cosh, 1,
321 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
322FUNC1(exp, exp, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000323 "exp(x)\n\nReturn e raised to the power of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000324FUNC1(fabs, fabs, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000325 "fabs(x)\n\nReturn the absolute value of the float x.")
Guido van Rossum13e05de2007-08-23 22:56:55 +0000326
327static PyObject * math_floor(PyObject *self, PyObject *number) {
328 static PyObject *floor_str = NULL;
329 PyObject *method;
330
331 if (floor_str == NULL) {
Christian Heimesfe82e772008-01-28 02:38:20 +0000332 floor_str = PyUnicode_InternFromString("__floor__");
Guido van Rossum13e05de2007-08-23 22:56:55 +0000333 if (floor_str == NULL)
334 return NULL;
335 }
336
Christian Heimes90aa7642007-12-19 02:45:37 +0000337 method = _PyType_Lookup(Py_TYPE(number), floor_str);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000338 if (method == NULL)
Christian Heimes53876d92008-04-19 00:31:39 +0000339 return math_1_to_int(number, floor, 0);
Guido van Rossum13e05de2007-08-23 22:56:55 +0000340 else
341 return PyObject_CallFunction(method, "O", number);
342}
343
344PyDoc_STRVAR(math_floor_doc,
Jeffrey Yasskinc2155832008-01-05 20:03:11 +0000345 "floor(x)\n\nReturn the floor of x as an int.\n"
Guido van Rossum13e05de2007-08-23 22:56:55 +0000346 "This is the largest integral value <= x.");
347
Christian Heimes53876d92008-04-19 00:31:39 +0000348FUNC1(log1p, log1p, 1,
349 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
350 The result is computed in a way which is accurate for x near zero.")
351FUNC1(sin, sin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000352 "sin(x)\n\nReturn the sine of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000353FUNC1(sinh, sinh, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000354 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000355FUNC1(sqrt, sqrt, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000356 "sqrt(x)\n\nReturn the square root of x.")
Christian Heimes53876d92008-04-19 00:31:39 +0000357FUNC1(tan, tan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000358 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Christian Heimes53876d92008-04-19 00:31:39 +0000359FUNC1(tanh, tanh, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000360 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000361
Barry Warsaw8b43b191996-12-09 22:32:36 +0000362static PyObject *
Christian Heimes400adb02008-02-01 08:12:03 +0000363math_trunc(PyObject *self, PyObject *number)
364{
365 static PyObject *trunc_str = NULL;
366 PyObject *trunc;
367
368 if (Py_TYPE(number)->tp_dict == NULL) {
369 if (PyType_Ready(Py_TYPE(number)) < 0)
370 return NULL;
371 }
372
373 if (trunc_str == NULL) {
374 trunc_str = PyUnicode_InternFromString("__trunc__");
375 if (trunc_str == NULL)
376 return NULL;
377 }
378
379 trunc = _PyType_Lookup(Py_TYPE(number), trunc_str);
380 if (trunc == NULL) {
381 PyErr_Format(PyExc_TypeError,
382 "type %.100s doesn't define __trunc__ method",
383 Py_TYPE(number)->tp_name);
384 return NULL;
385 }
386 return PyObject_CallFunctionObjArgs(trunc, number, NULL);
387}
388
389PyDoc_STRVAR(math_trunc_doc,
390"trunc(x:Real) -> Integral\n"
391"\n"
Christian Heimes292d3512008-02-03 16:51:08 +0000392"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
Christian Heimes400adb02008-02-01 08:12:03 +0000393
394static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000395math_frexp(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000396{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000397 int i;
Thomas Wouters89f507f2006-12-13 04:49:30 +0000398 double x = PyFloat_AsDouble(arg);
399 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000400 return NULL;
Christian Heimes53876d92008-04-19 00:31:39 +0000401 /* deal with special cases directly, to sidestep platform
402 differences */
403 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
404 i = 0;
405 }
406 else {
407 PyFPE_START_PROTECT("in math_frexp", return 0);
408 x = frexp(x, &i);
409 PyFPE_END_PROTECT(x);
410 }
411 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000412}
413
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000414PyDoc_STRVAR(math_frexp_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000415"frexp(x)\n"
416"\n"
417"Return the mantissa and exponent of x, as pair (m, e).\n"
418"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 +0000419"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 +0000420
Barry Warsaw8b43b191996-12-09 22:32:36 +0000421static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000422math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000423{
Christian Heimes53876d92008-04-19 00:31:39 +0000424 double x, r;
Guido van Rossumc5545052000-05-08 14:29:38 +0000425 int exp;
Fred Drake40c48682000-07-03 18:11:56 +0000426 if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000427 return NULL;
428 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000429 PyFPE_START_PROTECT("in math_ldexp", return 0)
430 r = ldexp(x, exp);
431 PyFPE_END_PROTECT(r)
432 if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
433 errno = ERANGE;
434 /* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
435 we unset it to avoid raising a ValueError here. */
436 if (errno == EDOM)
437 errno = 0;
438 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000439 return NULL;
Guido van Rossumd18ad581991-10-24 14:57:21 +0000440 else
Christian Heimes53876d92008-04-19 00:31:39 +0000441 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000442}
443
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000444PyDoc_STRVAR(math_ldexp_doc,
445"ldexp(x, i) -> x * (2**i)");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000446
Barry Warsaw8b43b191996-12-09 22:32:36 +0000447static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000448math_modf(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000449{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000450 double y, x = PyFloat_AsDouble(arg);
451 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000452 return NULL;
Christian Heimesa342c012008-04-20 21:01:16 +0000453 /* some platforms don't do the right thing for NaNs and
454 infinities, so we take care of special cases directly. */
455 if (!Py_IS_FINITE(x)) {
456 if (Py_IS_INFINITY(x))
457 return Py_BuildValue("(dd)", copysign(0., x), x);
458 else if (Py_IS_NAN(x))
459 return Py_BuildValue("(dd)", x, x);
460 }
461
Guido van Rossumd18ad581991-10-24 14:57:21 +0000462 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000463 PyFPE_START_PROTECT("in math_modf", return 0);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000464 x = modf(x, &y);
Christian Heimes53876d92008-04-19 00:31:39 +0000465 PyFPE_END_PROTECT(x);
466 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000467}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000468
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000469PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000470"modf(x)\n"
471"\n"
472"Return the fractional and integer parts of x. Both results carry the sign\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000473"of x. The integer part is returned as a real.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000474
Tim Peters78526162001-09-05 00:53:45 +0000475/* A decent logarithm is easy to compute even for huge longs, but libm can't
476 do that by itself -- loghelper can. func is log or log10, and name is
477 "log" or "log10". Note that overflow isn't possible: a long can contain
478 no more than INT_MAX * SHIFT bits, so has value certainly less than
479 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
480 small enough to fit in an IEEE single. log and log10 are even smaller.
481*/
482
483static PyObject*
Thomas Wouters89f507f2006-12-13 04:49:30 +0000484loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +0000485{
Tim Peters78526162001-09-05 00:53:45 +0000486 /* If it is long, do it ourselves. */
487 if (PyLong_Check(arg)) {
488 double x;
489 int e;
490 x = _PyLong_AsScaledDouble(arg, &e);
491 if (x <= 0.0) {
492 PyErr_SetString(PyExc_ValueError,
493 "math domain error");
494 return NULL;
495 }
Christian Heimesaf98da12008-01-27 15:18:18 +0000496 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
497 log(x) + log(2) * e * PyLong_SHIFT.
498 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
Tim Peters78526162001-09-05 00:53:45 +0000499 so force use of double. */
Martin v. Löwis9f2e3462007-07-21 17:22:18 +0000500 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000501 return PyFloat_FromDouble(x);
502 }
503
504 /* Else let libm handle it by itself. */
Christian Heimes53876d92008-04-19 00:31:39 +0000505 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +0000506}
507
508static PyObject *
509math_log(PyObject *self, PyObject *args)
510{
Raymond Hettinger866964c2002-12-14 19:51:34 +0000511 PyObject *arg;
512 PyObject *base = NULL;
513 PyObject *num, *den;
514 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000515
Raymond Hettingerea3fdf42002-12-29 16:33:45 +0000516 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
Raymond Hettinger866964c2002-12-14 19:51:34 +0000517 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000518
Thomas Wouters89f507f2006-12-13 04:49:30 +0000519 num = loghelper(arg, log, "log");
520 if (num == NULL || base == NULL)
521 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000522
Thomas Wouters89f507f2006-12-13 04:49:30 +0000523 den = loghelper(base, log, "log");
Raymond Hettinger866964c2002-12-14 19:51:34 +0000524 if (den == NULL) {
525 Py_DECREF(num);
526 return NULL;
527 }
528
Neal Norwitzbcc0db82006-03-24 08:14:36 +0000529 ans = PyNumber_TrueDivide(num, den);
Raymond Hettinger866964c2002-12-14 19:51:34 +0000530 Py_DECREF(num);
531 Py_DECREF(den);
532 return ans;
Tim Peters78526162001-09-05 00:53:45 +0000533}
534
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000535PyDoc_STRVAR(math_log_doc,
Raymond Hettinger866964c2002-12-14 19:51:34 +0000536"log(x[, base]) -> the logarithm of x to the given base.\n\
537If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +0000538
539static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000540math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +0000541{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000542 return loghelper(arg, log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +0000543}
544
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000545PyDoc_STRVAR(math_log10_doc,
546"log10(x) -> the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +0000547
Christian Heimes53876d92008-04-19 00:31:39 +0000548static PyObject *
549math_fmod(PyObject *self, PyObject *args)
550{
551 PyObject *ox, *oy;
552 double r, x, y;
553 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
554 return NULL;
555 x = PyFloat_AsDouble(ox);
556 y = PyFloat_AsDouble(oy);
557 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
558 return NULL;
559 /* fmod(x, +/-Inf) returns x for finite x. */
560 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
561 return PyFloat_FromDouble(x);
562 errno = 0;
563 PyFPE_START_PROTECT("in math_fmod", return 0);
564 r = fmod(x, y);
565 PyFPE_END_PROTECT(r);
566 if (Py_IS_NAN(r)) {
567 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
568 errno = EDOM;
569 else
570 errno = 0;
571 }
572 if (errno && is_error(r))
573 return NULL;
574 else
575 return PyFloat_FromDouble(r);
576}
577
578PyDoc_STRVAR(math_fmod_doc,
579"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
580" x % y may differ.");
581
582static PyObject *
583math_hypot(PyObject *self, PyObject *args)
584{
585 PyObject *ox, *oy;
586 double r, x, y;
587 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
588 return NULL;
589 x = PyFloat_AsDouble(ox);
590 y = PyFloat_AsDouble(oy);
591 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
592 return NULL;
593 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
594 if (Py_IS_INFINITY(x))
595 return PyFloat_FromDouble(fabs(x));
596 if (Py_IS_INFINITY(y))
597 return PyFloat_FromDouble(fabs(y));
598 errno = 0;
599 PyFPE_START_PROTECT("in math_hypot", return 0);
600 r = hypot(x, y);
601 PyFPE_END_PROTECT(r);
602 if (Py_IS_NAN(r)) {
603 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
604 errno = EDOM;
605 else
606 errno = 0;
607 }
608 else if (Py_IS_INFINITY(r)) {
609 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
610 errno = ERANGE;
611 else
612 errno = 0;
613 }
614 if (errno && is_error(r))
615 return NULL;
616 else
617 return PyFloat_FromDouble(r);
618}
619
620PyDoc_STRVAR(math_hypot_doc,
621"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
622
623/* pow can't use math_2, but needs its own wrapper: the problem is
624 that an infinite result can arise either as a result of overflow
625 (in which case OverflowError should be raised) or as a result of
626 e.g. 0.**-5. (for which ValueError needs to be raised.)
627*/
628
629static PyObject *
630math_pow(PyObject *self, PyObject *args)
631{
632 PyObject *ox, *oy;
633 double r, x, y;
Christian Heimesa342c012008-04-20 21:01:16 +0000634 int odd_y;
Christian Heimes53876d92008-04-19 00:31:39 +0000635
636 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
637 return NULL;
638 x = PyFloat_AsDouble(ox);
639 y = PyFloat_AsDouble(oy);
640 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
641 return NULL;
Christian Heimesa342c012008-04-20 21:01:16 +0000642
643 /* deal directly with IEEE specials, to cope with problems on various
644 platforms whose semantics don't exactly match C99 */
645 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
646 errno = 0;
647 if (Py_IS_NAN(x))
648 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
649 else if (Py_IS_NAN(y))
650 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
651 else if (Py_IS_INFINITY(x)) {
652 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
653 if (y > 0.)
654 r = odd_y ? x : fabs(x);
655 else if (y == 0.)
656 r = 1.;
657 else /* y < 0. */
658 r = odd_y ? copysign(0., x) : 0.;
659 }
660 else if (Py_IS_INFINITY(y)) {
661 if (fabs(x) == 1.0)
662 r = 1.;
663 else if (y > 0. && fabs(x) > 1.0)
664 r = y;
665 else if (y < 0. && fabs(x) < 1.0) {
666 r = -y; /* result is +inf */
667 if (x == 0.) /* 0**-inf: divide-by-zero */
668 errno = EDOM;
669 }
670 else
671 r = 0.;
672 }
Christian Heimes53876d92008-04-19 00:31:39 +0000673 }
Christian Heimesa342c012008-04-20 21:01:16 +0000674 else {
675 /* let libm handle finite**finite */
676 errno = 0;
677 PyFPE_START_PROTECT("in math_pow", return 0);
678 r = pow(x, y);
679 PyFPE_END_PROTECT(r);
680 /* a NaN result should arise only from (-ve)**(finite
681 non-integer); in this case we want to raise ValueError. */
682 if (!Py_IS_FINITE(r)) {
683 if (Py_IS_NAN(r)) {
684 errno = EDOM;
685 }
686 /*
687 an infinite result here arises either from:
688 (A) (+/-0.)**negative (-> divide-by-zero)
689 (B) overflow of x**y with x and y finite
690 */
691 else if (Py_IS_INFINITY(r)) {
692 if (x == 0.)
693 errno = EDOM;
694 else
695 errno = ERANGE;
696 }
697 }
Christian Heimes53876d92008-04-19 00:31:39 +0000698 }
699
700 if (errno && is_error(r))
701 return NULL;
702 else
703 return PyFloat_FromDouble(r);
704}
705
706PyDoc_STRVAR(math_pow_doc,
707"pow(x,y)\n\nReturn x**y (x to the power of y).");
708
Christian Heimes072c0f12008-01-03 23:01:04 +0000709static const double degToRad = Py_MATH_PI / 180.0;
710static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000711
712static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000713math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000714{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000715 double x = PyFloat_AsDouble(arg);
716 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000717 return NULL;
Christian Heimes072c0f12008-01-03 23:01:04 +0000718 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000719}
720
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000721PyDoc_STRVAR(math_degrees_doc,
722"degrees(x) -> converts angle x from radians to degrees");
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000723
724static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000725math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000726{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000727 double x = PyFloat_AsDouble(arg);
728 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000729 return NULL;
730 return PyFloat_FromDouble(x * degToRad);
731}
732
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000733PyDoc_STRVAR(math_radians_doc,
734"radians(x) -> converts angle x from degrees to radians");
Tim Peters78526162001-09-05 00:53:45 +0000735
Christian Heimes072c0f12008-01-03 23:01:04 +0000736static PyObject *
737math_isnan(PyObject *self, PyObject *arg)
738{
739 double x = PyFloat_AsDouble(arg);
740 if (x == -1.0 && PyErr_Occurred())
741 return NULL;
742 return PyBool_FromLong((long)Py_IS_NAN(x));
743}
744
745PyDoc_STRVAR(math_isnan_doc,
746"isnan(x) -> bool\n\
747Checks if float x is not a number (NaN)");
748
749static PyObject *
750math_isinf(PyObject *self, PyObject *arg)
751{
752 double x = PyFloat_AsDouble(arg);
753 if (x == -1.0 && PyErr_Occurred())
754 return NULL;
755 return PyBool_FromLong((long)Py_IS_INFINITY(x));
756}
757
758PyDoc_STRVAR(math_isinf_doc,
759"isinf(x) -> bool\n\
760Checks if float x is infinite (positive or negative)");
761
Barry Warsaw8b43b191996-12-09 22:32:36 +0000762static PyMethodDef math_methods[] = {
Thomas Wouters89f507f2006-12-13 04:49:30 +0000763 {"acos", math_acos, METH_O, math_acos_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000764 {"acosh", math_acosh, METH_O, math_acosh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000765 {"asin", math_asin, METH_O, math_asin_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000766 {"asinh", math_asinh, METH_O, math_asinh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000767 {"atan", math_atan, METH_O, math_atan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000768 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000769 {"atanh", math_atanh, METH_O, math_atanh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000770 {"ceil", math_ceil, METH_O, math_ceil_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000771 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000772 {"cos", math_cos, METH_O, math_cos_doc},
773 {"cosh", math_cosh, METH_O, math_cosh_doc},
774 {"degrees", math_degrees, METH_O, math_degrees_doc},
775 {"exp", math_exp, METH_O, math_exp_doc},
776 {"fabs", math_fabs, METH_O, math_fabs_doc},
777 {"floor", math_floor, METH_O, math_floor_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000778 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000779 {"frexp", math_frexp, METH_O, math_frexp_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000780 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000781 {"isinf", math_isinf, METH_O, math_isinf_doc},
782 {"isnan", math_isnan, METH_O, math_isnan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000783 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
784 {"log", math_log, METH_VARARGS, math_log_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000785 {"log1p", math_log1p, METH_O, math_log1p_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000786 {"log10", math_log10, METH_O, math_log10_doc},
787 {"modf", math_modf, METH_O, math_modf_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000788 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000789 {"radians", math_radians, METH_O, math_radians_doc},
790 {"sin", math_sin, METH_O, math_sin_doc},
791 {"sinh", math_sinh, METH_O, math_sinh_doc},
792 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
793 {"tan", math_tan, METH_O, math_tan_doc},
794 {"tanh", math_tanh, METH_O, math_tanh_doc},
Christian Heimes400adb02008-02-01 08:12:03 +0000795 {"trunc", math_trunc, METH_O, math_trunc_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000796 {NULL, NULL} /* sentinel */
797};
798
Guido van Rossumc6e22901998-12-04 19:26:43 +0000799
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000800PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000801"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000802"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000803
Mark Hammondfe51c6d2002-08-02 02:27:13 +0000804PyMODINIT_FUNC
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000805initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000806{
Christian Heimes53876d92008-04-19 00:31:39 +0000807 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +0000808
Guido van Rossumc6e22901998-12-04 19:26:43 +0000809 m = Py_InitModule3("math", math_methods, module_doc);
Neal Norwitz1ac754f2006-01-19 06:09:39 +0000810 if (m == NULL)
811 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +0000812
Christian Heimes53876d92008-04-19 00:31:39 +0000813 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
814 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +0000815
Christian Heimes53876d92008-04-19 00:31:39 +0000816 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000817 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000818}