blob: bb97a2c1023f15b1fd9e975ea7c44ed8245e6087 [file] [log] [blame]
Guido van Rossum85a5fbb1990-10-14 12:07:46 +00001/* Math module -- standard C math library functions, pi and e */
2
Christian Heimes6f341092008-04-18 23:13:07 +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
Neal Norwitz5f95a792008-01-25 08:04:16 +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 Heimes6f341092008-04-18 23:13:07 +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 *
Christian Heimes6f341092008-04-18 23:13:07 +0000129math_1(PyObject *arg, double (*func) (double), int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000130{
Christian Heimes6f341092008-04-18 23:13:07 +0000131 double x, r;
132 x = PyFloat_AsDouble(arg);
Neal Norwitz45e230a2006-11-19 21:26:53 +0000133 if (x == -1.0 && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000134 return NULL;
135 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000136 PyFPE_START_PROTECT("in math_1", return 0);
137 r = (*func)(x);
138 PyFPE_END_PROTECT(r);
139 if (Py_IS_NAN(r)) {
140 if (!Py_IS_NAN(x))
141 errno = EDOM;
142 else
143 errno = 0;
144 }
145 else if (Py_IS_INFINITY(r)) {
146 if (Py_IS_FINITE(x))
147 errno = can_overflow ? ERANGE : EDOM;
148 else
149 errno = 0;
150 }
151 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000152 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000153 else
Christian Heimes6f341092008-04-18 23:13:07 +0000154 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000155}
156
Christian Heimes6f341092008-04-18 23:13:07 +0000157/*
158 math_2 is used to wrap a libm function f that takes two double
159 arguments and returns a double.
160
161 The error reporting follows these rules, which are designed to do
162 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
163 platforms.
164
165 - a NaN result from non-NaN inputs causes ValueError to be raised
166 - an infinite result from finite inputs causes OverflowError to be
167 raised.
168 - if the result is finite and errno == EDOM then ValueError is
169 raised
170 - if the result is finite and nonzero and errno == ERANGE then
171 OverflowError is raised
172
173 The last rule is used to catch overflow on platforms which follow
174 C89 but for which HUGE_VAL is not an infinity.
175
176 For most two-argument functions (copysign, fmod, hypot, atan2)
177 these rules are enough to ensure that Python's functions behave as
178 specified in 'Annex F' of the C99 standard, with the 'invalid' and
179 'divide-by-zero' floating-point exceptions mapping to Python's
180 ValueError and the 'overflow' floating-point exception mapping to
181 OverflowError.
182*/
183
Barry Warsaw8b43b191996-12-09 22:32:36 +0000184static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000185math_2(PyObject *args, double (*func) (double, double), char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000186{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000187 PyObject *ox, *oy;
Christian Heimes6f341092008-04-18 23:13:07 +0000188 double x, y, r;
Neal Norwitz45e230a2006-11-19 21:26:53 +0000189 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
190 return NULL;
191 x = PyFloat_AsDouble(ox);
192 y = PyFloat_AsDouble(oy);
193 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000194 return NULL;
195 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000196 PyFPE_START_PROTECT("in math_2", return 0);
197 r = (*func)(x, y);
198 PyFPE_END_PROTECT(r);
199 if (Py_IS_NAN(r)) {
200 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
201 errno = EDOM;
202 else
203 errno = 0;
204 }
205 else if (Py_IS_INFINITY(r)) {
206 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
207 errno = ERANGE;
208 else
209 errno = 0;
210 }
211 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000212 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000213 else
Christian Heimes6f341092008-04-18 23:13:07 +0000214 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000215}
216
Christian Heimes6f341092008-04-18 23:13:07 +0000217#define FUNC1(funcname, func, can_overflow, docstring) \
Fred Drake40c48682000-07-03 18:11:56 +0000218 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Christian Heimes6f341092008-04-18 23:13:07 +0000219 return math_1(args, func, can_overflow); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000220 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000221 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000222
Fred Drake40c48682000-07-03 18:11:56 +0000223#define FUNC2(funcname, func, docstring) \
224 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Neal Norwitz45e230a2006-11-19 21:26:53 +0000225 return math_2(args, func, #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000226 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000227 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000228
Christian Heimes6f341092008-04-18 23:13:07 +0000229FUNC1(acos, acos, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000230 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000231FUNC1(acosh, acosh, 0,
232 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
233FUNC1(asin, asin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000234 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000235FUNC1(asinh, asinh, 0,
236 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
237FUNC1(atan, atan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000238 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Martin v. Löwis387c5472001-09-06 08:16:17 +0000239FUNC2(atan2, atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000240 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
241 "Unlike atan(y/x), the signs of both x and y are considered.")
Christian Heimes6f341092008-04-18 23:13:07 +0000242FUNC1(atanh, atanh, 0,
243 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
244FUNC1(ceil, ceil, 0,
Jeffrey Yasskin9871d8f2008-01-05 08:47:13 +0000245 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
246 "This is the smallest integral value >= x.")
Christian Heimeseebb79c2008-01-03 22:32:26 +0000247FUNC2(copysign, copysign,
Christian Heimes6f341092008-04-18 23:13:07 +0000248 "copysign(x,y)\n\nReturn x with the sign of y.")
249FUNC1(cos, cos, 0,
250 "cos(x)\n\nReturn the cosine of x (measured in radians).")
251FUNC1(cosh, cosh, 1,
252 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
253FUNC1(exp, exp, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000254 "exp(x)\n\nReturn e raised to the power of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000255FUNC1(fabs, fabs, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000256 "fabs(x)\n\nReturn the absolute value of the float x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000257FUNC1(floor, floor, 0,
Jeffrey Yasskin9871d8f2008-01-05 08:47:13 +0000258 "floor(x)\n\nReturn the floor of x as a float.\n"
259 "This is the largest integral value <= x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000260FUNC1(log1p, log1p, 1,
261 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
262 The result is computed in a way which is accurate for x near zero.")
263FUNC1(sin, sin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000264 "sin(x)\n\nReturn the sine of x (measured in radians).")
Christian Heimes6f341092008-04-18 23:13:07 +0000265FUNC1(sinh, sinh, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000266 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000267FUNC1(sqrt, sqrt, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000268 "sqrt(x)\n\nReturn the square root of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000269FUNC1(tan, tan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000270 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Christian Heimes6f341092008-04-18 23:13:07 +0000271FUNC1(tanh, tanh, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000272 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000273
Barry Warsaw8b43b191996-12-09 22:32:36 +0000274static PyObject *
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000275math_trunc(PyObject *self, PyObject *number)
276{
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000277 return PyObject_CallMethod(number, "__trunc__", NULL);
278}
279
280PyDoc_STRVAR(math_trunc_doc,
281"trunc(x:Real) -> Integral\n"
282"\n"
Raymond Hettingerfe424f72008-02-02 05:24:44 +0000283"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000284
285static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000286math_frexp(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000287{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000288 int i;
Neal Norwitz45e230a2006-11-19 21:26:53 +0000289 double x = PyFloat_AsDouble(arg);
290 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000291 return NULL;
Christian Heimes6f341092008-04-18 23:13:07 +0000292 /* deal with special cases directly, to sidestep platform
293 differences */
294 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
295 i = 0;
296 }
297 else {
298 PyFPE_START_PROTECT("in math_frexp", return 0);
299 x = frexp(x, &i);
300 PyFPE_END_PROTECT(x);
301 }
302 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000303}
304
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000305PyDoc_STRVAR(math_frexp_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000306"frexp(x)\n"
307"\n"
308"Return the mantissa and exponent of x, as pair (m, e).\n"
309"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 +0000310"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 +0000311
Barry Warsaw8b43b191996-12-09 22:32:36 +0000312static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000313math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000314{
Christian Heimes6f341092008-04-18 23:13:07 +0000315 double x, r;
Guido van Rossumc5545052000-05-08 14:29:38 +0000316 int exp;
Fred Drake40c48682000-07-03 18:11:56 +0000317 if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000318 return NULL;
319 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000320 PyFPE_START_PROTECT("in math_ldexp", return 0)
321 r = ldexp(x, exp);
322 PyFPE_END_PROTECT(r)
323 if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
324 errno = ERANGE;
325 /* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
326 we unset it to avoid raising a ValueError here. */
327 if (errno == EDOM)
328 errno = 0;
329 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000330 return NULL;
Guido van Rossumd18ad581991-10-24 14:57:21 +0000331 else
Christian Heimes6f341092008-04-18 23:13:07 +0000332 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000333}
334
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000335PyDoc_STRVAR(math_ldexp_doc,
336"ldexp(x, i) -> x * (2**i)");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000337
Barry Warsaw8b43b191996-12-09 22:32:36 +0000338static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000339math_modf(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000340{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000341 double y, x = PyFloat_AsDouble(arg);
342 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000343 return NULL;
344 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000345 PyFPE_START_PROTECT("in math_modf", return 0);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000346 x = modf(x, &y);
Christian Heimes6f341092008-04-18 23:13:07 +0000347 PyFPE_END_PROTECT(x);
348 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000349}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000350
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000351PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000352"modf(x)\n"
353"\n"
354"Return the fractional and integer parts of x. Both results carry the sign\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000355"of x. The integer part is returned as a real.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000356
Tim Peters78526162001-09-05 00:53:45 +0000357/* A decent logarithm is easy to compute even for huge longs, but libm can't
358 do that by itself -- loghelper can. func is log or log10, and name is
359 "log" or "log10". Note that overflow isn't possible: a long can contain
360 no more than INT_MAX * SHIFT bits, so has value certainly less than
361 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
362 small enough to fit in an IEEE single. log and log10 are even smaller.
363*/
364
365static PyObject*
Neal Norwitz45e230a2006-11-19 21:26:53 +0000366loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +0000367{
Tim Peters78526162001-09-05 00:53:45 +0000368 /* If it is long, do it ourselves. */
369 if (PyLong_Check(arg)) {
370 double x;
371 int e;
372 x = _PyLong_AsScaledDouble(arg, &e);
373 if (x <= 0.0) {
374 PyErr_SetString(PyExc_ValueError,
375 "math domain error");
376 return NULL;
377 }
Christian Heimes543cabc2008-01-25 14:54:23 +0000378 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
379 log(x) + log(2) * e * PyLong_SHIFT.
380 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
Tim Peters78526162001-09-05 00:53:45 +0000381 so force use of double. */
Christian Heimes543cabc2008-01-25 14:54:23 +0000382 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000383 return PyFloat_FromDouble(x);
384 }
385
386 /* Else let libm handle it by itself. */
Christian Heimes6f341092008-04-18 23:13:07 +0000387 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +0000388}
389
390static PyObject *
391math_log(PyObject *self, PyObject *args)
392{
Raymond Hettinger866964c2002-12-14 19:51:34 +0000393 PyObject *arg;
394 PyObject *base = NULL;
395 PyObject *num, *den;
396 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000397
Raymond Hettingerea3fdf42002-12-29 16:33:45 +0000398 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
Raymond Hettinger866964c2002-12-14 19:51:34 +0000399 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000400
Neal Norwitz45e230a2006-11-19 21:26:53 +0000401 num = loghelper(arg, log, "log");
402 if (num == NULL || base == NULL)
403 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000404
Neal Norwitz45e230a2006-11-19 21:26:53 +0000405 den = loghelper(base, log, "log");
Raymond Hettinger866964c2002-12-14 19:51:34 +0000406 if (den == NULL) {
407 Py_DECREF(num);
408 return NULL;
409 }
410
411 ans = PyNumber_Divide(num, den);
412 Py_DECREF(num);
413 Py_DECREF(den);
414 return ans;
Tim Peters78526162001-09-05 00:53:45 +0000415}
416
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000417PyDoc_STRVAR(math_log_doc,
Raymond Hettinger866964c2002-12-14 19:51:34 +0000418"log(x[, base]) -> the logarithm of x to the given base.\n\
419If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +0000420
421static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000422math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +0000423{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000424 return loghelper(arg, log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +0000425}
426
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000427PyDoc_STRVAR(math_log10_doc,
428"log10(x) -> the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +0000429
Christian Heimes6f341092008-04-18 23:13:07 +0000430static PyObject *
431math_fmod(PyObject *self, PyObject *args)
432{
433 PyObject *ox, *oy;
434 double r, x, y;
435 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
436 return NULL;
437 x = PyFloat_AsDouble(ox);
438 y = PyFloat_AsDouble(oy);
439 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
440 return NULL;
441 /* fmod(x, +/-Inf) returns x for finite x. */
442 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
443 return PyFloat_FromDouble(x);
444 errno = 0;
445 PyFPE_START_PROTECT("in math_fmod", return 0);
446 r = fmod(x, y);
447 PyFPE_END_PROTECT(r);
448 if (Py_IS_NAN(r)) {
449 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
450 errno = EDOM;
451 else
452 errno = 0;
453 }
454 if (errno && is_error(r))
455 return NULL;
456 else
457 return PyFloat_FromDouble(r);
458}
459
460PyDoc_STRVAR(math_fmod_doc,
461"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
462" x % y may differ.");
463
464static PyObject *
465math_hypot(PyObject *self, PyObject *args)
466{
467 PyObject *ox, *oy;
468 double r, x, y;
469 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
470 return NULL;
471 x = PyFloat_AsDouble(ox);
472 y = PyFloat_AsDouble(oy);
473 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
474 return NULL;
475 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
476 if (Py_IS_INFINITY(x))
477 return PyFloat_FromDouble(fabs(x));
478 if (Py_IS_INFINITY(y))
479 return PyFloat_FromDouble(fabs(y));
480 errno = 0;
481 PyFPE_START_PROTECT("in math_hypot", return 0);
482 r = hypot(x, y);
483 PyFPE_END_PROTECT(r);
484 if (Py_IS_NAN(r)) {
485 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
486 errno = EDOM;
487 else
488 errno = 0;
489 }
490 else if (Py_IS_INFINITY(r)) {
491 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
492 errno = ERANGE;
493 else
494 errno = 0;
495 }
496 if (errno && is_error(r))
497 return NULL;
498 else
499 return PyFloat_FromDouble(r);
500}
501
502PyDoc_STRVAR(math_hypot_doc,
503"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
504
505/* pow can't use math_2, but needs its own wrapper: the problem is
506 that an infinite result can arise either as a result of overflow
507 (in which case OverflowError should be raised) or as a result of
508 e.g. 0.**-5. (for which ValueError needs to be raised.)
509*/
510
511static PyObject *
512math_pow(PyObject *self, PyObject *args)
513{
514 PyObject *ox, *oy;
515 double r, x, y;
Mark Dickinsone941d972008-04-19 18:51:48 +0000516 int y_is_odd;
Christian Heimes6f341092008-04-18 23:13:07 +0000517
518 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
519 return NULL;
520 x = PyFloat_AsDouble(ox);
521 y = PyFloat_AsDouble(oy);
522 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
523 return NULL;
Mark Dickinsona1293eb2008-04-19 19:41:52 +0000524
525 /* deal directly with various special cases, to cope with problems on
526 various platforms whose semantics don't exactly match C99 */
527
528 /* 1**x, x**0, and (-1)**(+-infinity) return 1., even if x is NaN or
529 an infinity. */
530 if (x == 1. || y == 0. || (x == -1. && Py_IS_INFINITY(y)))
Christian Heimes6f341092008-04-18 23:13:07 +0000531 return PyFloat_FromDouble(1.);
Mark Dickinsona1293eb2008-04-19 19:41:52 +0000532 /* otherwise, return a NaN if either input was a NaN */
533 if (Py_IS_NAN(x))
534 return PyFloat_FromDouble(x);
535 if (Py_IS_NAN(y))
536 return PyFloat_FromDouble(y);
Mark Dickinsone941d972008-04-19 18:51:48 +0000537 /* inf ** (nonzero, non-NaN) is one of +-0, +-infinity */
538 if (Py_IS_INFINITY(x) && !Py_IS_NAN(y)) {
539 y_is_odd = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
540 if (y > 0.)
541 r = y_is_odd ? x : fabs(x);
542 else
543 r = y_is_odd ? copysign(0., x) : 0.;
544 return PyFloat_FromDouble(r);
545 }
546
Christian Heimes6f341092008-04-18 23:13:07 +0000547 errno = 0;
548 PyFPE_START_PROTECT("in math_pow", return 0);
549 r = pow(x, y);
550 PyFPE_END_PROTECT(r);
551 if (Py_IS_NAN(r)) {
Mark Dickinsona1293eb2008-04-19 19:41:52 +0000552 errno = EDOM;
Christian Heimes6f341092008-04-18 23:13:07 +0000553 }
554 /* an infinite result arises either from:
555
556 (A) (+/-0.)**negative,
557 (B) overflow of x**y with both x and y finite (and x nonzero)
558 (C) (+/-inf)**positive, or
559 (D) x**inf with |x| > 1, or x**-inf with |x| < 1.
560
561 In case (A) we want ValueError to be raised. In case (B)
562 OverflowError should be raised. In cases (C) and (D) the infinite
563 result should be returned.
564 */
565 else if (Py_IS_INFINITY(r)) {
566 if (x == 0.)
567 errno = EDOM;
568 else if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
569 errno = ERANGE;
570 else
571 errno = 0;
572 }
573
574 if (errno && is_error(r))
575 return NULL;
576 else
577 return PyFloat_FromDouble(r);
578}
579
580PyDoc_STRVAR(math_pow_doc,
581"pow(x,y)\n\nReturn x**y (x to the power of y).");
582
Christian Heimese2ca4242008-01-03 20:23:15 +0000583static const double degToRad = Py_MATH_PI / 180.0;
584static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000585
586static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000587math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000588{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000589 double x = PyFloat_AsDouble(arg);
590 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000591 return NULL;
Christian Heimese2ca4242008-01-03 20:23:15 +0000592 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000593}
594
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000595PyDoc_STRVAR(math_degrees_doc,
596"degrees(x) -> converts angle x from radians to degrees");
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000597
598static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000599math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000600{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000601 double x = PyFloat_AsDouble(arg);
602 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000603 return NULL;
604 return PyFloat_FromDouble(x * degToRad);
605}
606
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000607PyDoc_STRVAR(math_radians_doc,
608"radians(x) -> converts angle x from degrees to radians");
Tim Peters78526162001-09-05 00:53:45 +0000609
Christian Heimese2ca4242008-01-03 20:23:15 +0000610static PyObject *
611math_isnan(PyObject *self, PyObject *arg)
612{
613 double x = PyFloat_AsDouble(arg);
614 if (x == -1.0 && PyErr_Occurred())
615 return NULL;
616 return PyBool_FromLong((long)Py_IS_NAN(x));
617}
618
619PyDoc_STRVAR(math_isnan_doc,
620"isnan(x) -> bool\n\
621Checks if float x is not a number (NaN)");
622
623static PyObject *
624math_isinf(PyObject *self, PyObject *arg)
625{
626 double x = PyFloat_AsDouble(arg);
627 if (x == -1.0 && PyErr_Occurred())
628 return NULL;
629 return PyBool_FromLong((long)Py_IS_INFINITY(x));
630}
631
632PyDoc_STRVAR(math_isinf_doc,
633"isinf(x) -> bool\n\
634Checks if float x is infinite (positive or negative)");
635
Barry Warsaw8b43b191996-12-09 22:32:36 +0000636static PyMethodDef math_methods[] = {
Neal Norwitz45e230a2006-11-19 21:26:53 +0000637 {"acos", math_acos, METH_O, math_acos_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000638 {"acosh", math_acosh, METH_O, math_acosh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000639 {"asin", math_asin, METH_O, math_asin_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000640 {"asinh", math_asinh, METH_O, math_asinh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000641 {"atan", math_atan, METH_O, math_atan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000642 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000643 {"atanh", math_atanh, METH_O, math_atanh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000644 {"ceil", math_ceil, METH_O, math_ceil_doc},
Christian Heimeseebb79c2008-01-03 22:32:26 +0000645 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000646 {"cos", math_cos, METH_O, math_cos_doc},
647 {"cosh", math_cosh, METH_O, math_cosh_doc},
648 {"degrees", math_degrees, METH_O, math_degrees_doc},
649 {"exp", math_exp, METH_O, math_exp_doc},
650 {"fabs", math_fabs, METH_O, math_fabs_doc},
651 {"floor", math_floor, METH_O, math_floor_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000652 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000653 {"frexp", math_frexp, METH_O, math_frexp_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000654 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Christian Heimese2ca4242008-01-03 20:23:15 +0000655 {"isinf", math_isinf, METH_O, math_isinf_doc},
656 {"isnan", math_isnan, METH_O, math_isnan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000657 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
658 {"log", math_log, METH_VARARGS, math_log_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000659 {"log1p", math_log1p, METH_O, math_log1p_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000660 {"log10", math_log10, METH_O, math_log10_doc},
661 {"modf", math_modf, METH_O, math_modf_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000662 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000663 {"radians", math_radians, METH_O, math_radians_doc},
664 {"sin", math_sin, METH_O, math_sin_doc},
665 {"sinh", math_sinh, METH_O, math_sinh_doc},
666 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
667 {"tan", math_tan, METH_O, math_tan_doc},
668 {"tanh", math_tanh, METH_O, math_tanh_doc},
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000669 {"trunc", math_trunc, METH_O, math_trunc_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000670 {NULL, NULL} /* sentinel */
671};
672
Guido van Rossumc6e22901998-12-04 19:26:43 +0000673
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000674PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000675"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000676"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000677
Mark Hammondfe51c6d2002-08-02 02:27:13 +0000678PyMODINIT_FUNC
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000679initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000680{
Christian Heimes6f341092008-04-18 23:13:07 +0000681 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +0000682
Guido van Rossumc6e22901998-12-04 19:26:43 +0000683 m = Py_InitModule3("math", math_methods, module_doc);
Neal Norwitz1ac754f2006-01-19 06:09:39 +0000684 if (m == NULL)
685 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +0000686
Christian Heimes6f341092008-04-18 23:13:07 +0000687 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
688 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +0000689
Christian Heimes6f341092008-04-18 23:13:07 +0000690 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000691 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000692}