blob: 19ed1b1b97c179afb7ccc35bb7ef15f8ab49800c [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;
Christian Heimesa342c012008-04-20 21:01:16 +0000417 /* some platforms don't do the right thing for NaNs and
418 infinities, so we take care of special cases directly. */
419 if (!Py_IS_FINITE(x)) {
420 if (Py_IS_INFINITY(x))
421 return Py_BuildValue("(dd)", copysign(0., x), x);
422 else if (Py_IS_NAN(x))
423 return Py_BuildValue("(dd)", x, x);
424 }
425
Guido van Rossumd18ad581991-10-24 14:57:21 +0000426 errno = 0;
Christian Heimes53876d92008-04-19 00:31:39 +0000427 PyFPE_START_PROTECT("in math_modf", return 0);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000428 x = modf(x, &y);
Christian Heimes53876d92008-04-19 00:31:39 +0000429 PyFPE_END_PROTECT(x);
430 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000431}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000432
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000433PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000434"modf(x)\n"
435"\n"
436"Return the fractional and integer parts of x. Both results carry the sign\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000437"of x. The integer part is returned as a real.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000438
Tim Peters78526162001-09-05 00:53:45 +0000439/* A decent logarithm is easy to compute even for huge longs, but libm can't
440 do that by itself -- loghelper can. func is log or log10, and name is
441 "log" or "log10". Note that overflow isn't possible: a long can contain
442 no more than INT_MAX * SHIFT bits, so has value certainly less than
443 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
444 small enough to fit in an IEEE single. log and log10 are even smaller.
445*/
446
447static PyObject*
Thomas Wouters89f507f2006-12-13 04:49:30 +0000448loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +0000449{
Tim Peters78526162001-09-05 00:53:45 +0000450 /* If it is long, do it ourselves. */
451 if (PyLong_Check(arg)) {
452 double x;
453 int e;
454 x = _PyLong_AsScaledDouble(arg, &e);
455 if (x <= 0.0) {
456 PyErr_SetString(PyExc_ValueError,
457 "math domain error");
458 return NULL;
459 }
Christian Heimesaf98da12008-01-27 15:18:18 +0000460 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
461 log(x) + log(2) * e * PyLong_SHIFT.
462 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
Tim Peters78526162001-09-05 00:53:45 +0000463 so force use of double. */
Martin v. Löwis9f2e3462007-07-21 17:22:18 +0000464 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000465 return PyFloat_FromDouble(x);
466 }
467
468 /* Else let libm handle it by itself. */
Christian Heimes53876d92008-04-19 00:31:39 +0000469 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +0000470}
471
472static PyObject *
473math_log(PyObject *self, PyObject *args)
474{
Raymond Hettinger866964c2002-12-14 19:51:34 +0000475 PyObject *arg;
476 PyObject *base = NULL;
477 PyObject *num, *den;
478 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000479
Raymond Hettingerea3fdf42002-12-29 16:33:45 +0000480 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
Raymond Hettinger866964c2002-12-14 19:51:34 +0000481 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000482
Thomas Wouters89f507f2006-12-13 04:49:30 +0000483 num = loghelper(arg, log, "log");
484 if (num == NULL || base == NULL)
485 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000486
Thomas Wouters89f507f2006-12-13 04:49:30 +0000487 den = loghelper(base, log, "log");
Raymond Hettinger866964c2002-12-14 19:51:34 +0000488 if (den == NULL) {
489 Py_DECREF(num);
490 return NULL;
491 }
492
Neal Norwitzbcc0db82006-03-24 08:14:36 +0000493 ans = PyNumber_TrueDivide(num, den);
Raymond Hettinger866964c2002-12-14 19:51:34 +0000494 Py_DECREF(num);
495 Py_DECREF(den);
496 return ans;
Tim Peters78526162001-09-05 00:53:45 +0000497}
498
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000499PyDoc_STRVAR(math_log_doc,
Raymond Hettinger866964c2002-12-14 19:51:34 +0000500"log(x[, base]) -> the logarithm of x to the given base.\n\
501If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +0000502
503static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000504math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +0000505{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000506 return loghelper(arg, log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +0000507}
508
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000509PyDoc_STRVAR(math_log10_doc,
510"log10(x) -> the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +0000511
Christian Heimes53876d92008-04-19 00:31:39 +0000512static PyObject *
513math_fmod(PyObject *self, PyObject *args)
514{
515 PyObject *ox, *oy;
516 double r, x, y;
517 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
518 return NULL;
519 x = PyFloat_AsDouble(ox);
520 y = PyFloat_AsDouble(oy);
521 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
522 return NULL;
523 /* fmod(x, +/-Inf) returns x for finite x. */
524 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
525 return PyFloat_FromDouble(x);
526 errno = 0;
527 PyFPE_START_PROTECT("in math_fmod", return 0);
528 r = fmod(x, y);
529 PyFPE_END_PROTECT(r);
530 if (Py_IS_NAN(r)) {
531 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
532 errno = EDOM;
533 else
534 errno = 0;
535 }
536 if (errno && is_error(r))
537 return NULL;
538 else
539 return PyFloat_FromDouble(r);
540}
541
542PyDoc_STRVAR(math_fmod_doc,
543"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
544" x % y may differ.");
545
546static PyObject *
547math_hypot(PyObject *self, PyObject *args)
548{
549 PyObject *ox, *oy;
550 double r, x, y;
551 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
552 return NULL;
553 x = PyFloat_AsDouble(ox);
554 y = PyFloat_AsDouble(oy);
555 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
556 return NULL;
557 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
558 if (Py_IS_INFINITY(x))
559 return PyFloat_FromDouble(fabs(x));
560 if (Py_IS_INFINITY(y))
561 return PyFloat_FromDouble(fabs(y));
562 errno = 0;
563 PyFPE_START_PROTECT("in math_hypot", return 0);
564 r = hypot(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 else if (Py_IS_INFINITY(r)) {
573 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
574 errno = ERANGE;
575 else
576 errno = 0;
577 }
578 if (errno && is_error(r))
579 return NULL;
580 else
581 return PyFloat_FromDouble(r);
582}
583
584PyDoc_STRVAR(math_hypot_doc,
585"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
586
587/* pow can't use math_2, but needs its own wrapper: the problem is
588 that an infinite result can arise either as a result of overflow
589 (in which case OverflowError should be raised) or as a result of
590 e.g. 0.**-5. (for which ValueError needs to be raised.)
591*/
592
593static PyObject *
594math_pow(PyObject *self, PyObject *args)
595{
596 PyObject *ox, *oy;
597 double r, x, y;
Christian Heimesa342c012008-04-20 21:01:16 +0000598 int odd_y;
Christian Heimes53876d92008-04-19 00:31:39 +0000599
600 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
601 return NULL;
602 x = PyFloat_AsDouble(ox);
603 y = PyFloat_AsDouble(oy);
604 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
605 return NULL;
Christian Heimesa342c012008-04-20 21:01:16 +0000606
607 /* deal directly with IEEE specials, to cope with problems on various
608 platforms whose semantics don't exactly match C99 */
609 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
610 errno = 0;
611 if (Py_IS_NAN(x))
612 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
613 else if (Py_IS_NAN(y))
614 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
615 else if (Py_IS_INFINITY(x)) {
616 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
617 if (y > 0.)
618 r = odd_y ? x : fabs(x);
619 else if (y == 0.)
620 r = 1.;
621 else /* y < 0. */
622 r = odd_y ? copysign(0., x) : 0.;
623 }
624 else if (Py_IS_INFINITY(y)) {
625 if (fabs(x) == 1.0)
626 r = 1.;
627 else if (y > 0. && fabs(x) > 1.0)
628 r = y;
629 else if (y < 0. && fabs(x) < 1.0) {
630 r = -y; /* result is +inf */
631 if (x == 0.) /* 0**-inf: divide-by-zero */
632 errno = EDOM;
633 }
634 else
635 r = 0.;
636 }
Christian Heimes53876d92008-04-19 00:31:39 +0000637 }
Christian Heimesa342c012008-04-20 21:01:16 +0000638 else {
639 /* let libm handle finite**finite */
640 errno = 0;
641 PyFPE_START_PROTECT("in math_pow", return 0);
642 r = pow(x, y);
643 PyFPE_END_PROTECT(r);
644 /* a NaN result should arise only from (-ve)**(finite
645 non-integer); in this case we want to raise ValueError. */
646 if (!Py_IS_FINITE(r)) {
647 if (Py_IS_NAN(r)) {
648 errno = EDOM;
649 }
650 /*
651 an infinite result here arises either from:
652 (A) (+/-0.)**negative (-> divide-by-zero)
653 (B) overflow of x**y with x and y finite
654 */
655 else if (Py_IS_INFINITY(r)) {
656 if (x == 0.)
657 errno = EDOM;
658 else
659 errno = ERANGE;
660 }
661 }
Christian Heimes53876d92008-04-19 00:31:39 +0000662 }
663
664 if (errno && is_error(r))
665 return NULL;
666 else
667 return PyFloat_FromDouble(r);
668}
669
670PyDoc_STRVAR(math_pow_doc,
671"pow(x,y)\n\nReturn x**y (x to the power of y).");
672
Christian Heimes072c0f12008-01-03 23:01:04 +0000673static const double degToRad = Py_MATH_PI / 180.0;
674static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000675
676static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000677math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000678{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000679 double x = PyFloat_AsDouble(arg);
680 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000681 return NULL;
Christian Heimes072c0f12008-01-03 23:01:04 +0000682 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000683}
684
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000685PyDoc_STRVAR(math_degrees_doc,
686"degrees(x) -> converts angle x from radians to degrees");
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000687
688static PyObject *
Thomas Wouters89f507f2006-12-13 04:49:30 +0000689math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000690{
Thomas Wouters89f507f2006-12-13 04:49:30 +0000691 double x = PyFloat_AsDouble(arg);
692 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000693 return NULL;
694 return PyFloat_FromDouble(x * degToRad);
695}
696
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000697PyDoc_STRVAR(math_radians_doc,
698"radians(x) -> converts angle x from degrees to radians");
Tim Peters78526162001-09-05 00:53:45 +0000699
Christian Heimes072c0f12008-01-03 23:01:04 +0000700static PyObject *
701math_isnan(PyObject *self, PyObject *arg)
702{
703 double x = PyFloat_AsDouble(arg);
704 if (x == -1.0 && PyErr_Occurred())
705 return NULL;
706 return PyBool_FromLong((long)Py_IS_NAN(x));
707}
708
709PyDoc_STRVAR(math_isnan_doc,
710"isnan(x) -> bool\n\
711Checks if float x is not a number (NaN)");
712
713static PyObject *
714math_isinf(PyObject *self, PyObject *arg)
715{
716 double x = PyFloat_AsDouble(arg);
717 if (x == -1.0 && PyErr_Occurred())
718 return NULL;
719 return PyBool_FromLong((long)Py_IS_INFINITY(x));
720}
721
722PyDoc_STRVAR(math_isinf_doc,
723"isinf(x) -> bool\n\
724Checks if float x is infinite (positive or negative)");
725
Barry Warsaw8b43b191996-12-09 22:32:36 +0000726static PyMethodDef math_methods[] = {
Thomas Wouters89f507f2006-12-13 04:49:30 +0000727 {"acos", math_acos, METH_O, math_acos_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000728 {"acosh", math_acosh, METH_O, math_acosh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000729 {"asin", math_asin, METH_O, math_asin_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000730 {"asinh", math_asinh, METH_O, math_asinh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000731 {"atan", math_atan, METH_O, math_atan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000732 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000733 {"atanh", math_atanh, METH_O, math_atanh_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000734 {"ceil", math_ceil, METH_O, math_ceil_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000735 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000736 {"cos", math_cos, METH_O, math_cos_doc},
737 {"cosh", math_cosh, METH_O, math_cosh_doc},
738 {"degrees", math_degrees, METH_O, math_degrees_doc},
739 {"exp", math_exp, METH_O, math_exp_doc},
740 {"fabs", math_fabs, METH_O, math_fabs_doc},
741 {"floor", math_floor, METH_O, math_floor_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000742 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000743 {"frexp", math_frexp, METH_O, math_frexp_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000744 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Christian Heimes072c0f12008-01-03 23:01:04 +0000745 {"isinf", math_isinf, METH_O, math_isinf_doc},
746 {"isnan", math_isnan, METH_O, math_isnan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000747 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
748 {"log", math_log, METH_VARARGS, math_log_doc},
Christian Heimes53876d92008-04-19 00:31:39 +0000749 {"log1p", math_log1p, METH_O, math_log1p_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000750 {"log10", math_log10, METH_O, math_log10_doc},
751 {"modf", math_modf, METH_O, math_modf_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000752 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Thomas Wouters89f507f2006-12-13 04:49:30 +0000753 {"radians", math_radians, METH_O, math_radians_doc},
754 {"sin", math_sin, METH_O, math_sin_doc},
755 {"sinh", math_sinh, METH_O, math_sinh_doc},
756 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
757 {"tan", math_tan, METH_O, math_tan_doc},
758 {"tanh", math_tanh, METH_O, math_tanh_doc},
Christian Heimes400adb02008-02-01 08:12:03 +0000759 {"trunc", math_trunc, METH_O, math_trunc_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000760 {NULL, NULL} /* sentinel */
761};
762
Guido van Rossumc6e22901998-12-04 19:26:43 +0000763
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000764PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000765"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000766"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000767
Mark Hammondfe51c6d2002-08-02 02:27:13 +0000768PyMODINIT_FUNC
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000769initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000770{
Christian Heimes53876d92008-04-19 00:31:39 +0000771 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +0000772
Guido van Rossumc6e22901998-12-04 19:26:43 +0000773 m = Py_InitModule3("math", math_methods, module_doc);
Neal Norwitz1ac754f2006-01-19 06:09:39 +0000774 if (m == NULL)
775 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +0000776
Christian Heimes53876d92008-04-19 00:31:39 +0000777 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
778 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +0000779
Christian Heimes53876d92008-04-19 00:31:39 +0000780 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000781 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000782}