blob: d1c4822cb1cebbb3825b78ede6d55915a87d9eb6 [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/*
Mark Dickinson92483cd2008-04-20 21:39:04 +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 Heimes6f341092008-04-18 23:13:07 +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 *
Christian Heimes6f341092008-04-18 23:13:07 +0000165math_1(PyObject *arg, double (*func) (double), int can_overflow)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000166{
Christian Heimes6f341092008-04-18 23:13:07 +0000167 double x, r;
168 x = PyFloat_AsDouble(arg);
Neal Norwitz45e230a2006-11-19 21:26:53 +0000169 if (x == -1.0 && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000170 return NULL;
171 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000172 PyFPE_START_PROTECT("in math_1", return 0);
173 r = (*func)(x);
174 PyFPE_END_PROTECT(r);
175 if (Py_IS_NAN(r)) {
176 if (!Py_IS_NAN(x))
177 errno = EDOM;
178 else
179 errno = 0;
180 }
181 else if (Py_IS_INFINITY(r)) {
182 if (Py_IS_FINITE(x))
183 errno = can_overflow ? ERANGE : EDOM;
184 else
185 errno = 0;
186 }
187 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000188 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000189 else
Christian Heimes6f341092008-04-18 23:13:07 +0000190 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000191}
192
Christian Heimes6f341092008-04-18 23:13:07 +0000193/*
194 math_2 is used to wrap a libm function f that takes two double
195 arguments and returns a double.
196
197 The error reporting follows these rules, which are designed to do
198 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
199 platforms.
200
201 - a NaN result from non-NaN inputs causes ValueError to be raised
202 - an infinite result from finite inputs causes OverflowError to be
203 raised.
204 - if the result is finite and errno == EDOM then ValueError is
205 raised
206 - if the result is finite and nonzero and errno == ERANGE then
207 OverflowError is raised
208
209 The last rule is used to catch overflow on platforms which follow
210 C89 but for which HUGE_VAL is not an infinity.
211
212 For most two-argument functions (copysign, fmod, hypot, atan2)
213 these rules are enough to ensure that Python's functions behave as
214 specified in 'Annex F' of the C99 standard, with the 'invalid' and
215 'divide-by-zero' floating-point exceptions mapping to Python's
216 ValueError and the 'overflow' floating-point exception mapping to
217 OverflowError.
218*/
219
Barry Warsaw8b43b191996-12-09 22:32:36 +0000220static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000221math_2(PyObject *args, double (*func) (double, double), char *funcname)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000222{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000223 PyObject *ox, *oy;
Christian Heimes6f341092008-04-18 23:13:07 +0000224 double x, y, r;
Neal Norwitz45e230a2006-11-19 21:26:53 +0000225 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
226 return NULL;
227 x = PyFloat_AsDouble(ox);
228 y = PyFloat_AsDouble(oy);
229 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000230 return NULL;
231 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000232 PyFPE_START_PROTECT("in math_2", return 0);
233 r = (*func)(x, y);
234 PyFPE_END_PROTECT(r);
235 if (Py_IS_NAN(r)) {
236 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
237 errno = EDOM;
238 else
239 errno = 0;
240 }
241 else if (Py_IS_INFINITY(r)) {
242 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
243 errno = ERANGE;
244 else
245 errno = 0;
246 }
247 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000248 return NULL;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000249 else
Christian Heimes6f341092008-04-18 23:13:07 +0000250 return PyFloat_FromDouble(r);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000251}
252
Christian Heimes6f341092008-04-18 23:13:07 +0000253#define FUNC1(funcname, func, can_overflow, docstring) \
Fred Drake40c48682000-07-03 18:11:56 +0000254 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Christian Heimes6f341092008-04-18 23:13:07 +0000255 return math_1(args, func, can_overflow); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000256 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000257 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000258
Fred Drake40c48682000-07-03 18:11:56 +0000259#define FUNC2(funcname, func, docstring) \
260 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
Neal Norwitz45e230a2006-11-19 21:26:53 +0000261 return math_2(args, func, #funcname); \
Guido van Rossumc6e22901998-12-04 19:26:43 +0000262 }\
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000263 PyDoc_STRVAR(math_##funcname##_doc, docstring);
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000264
Christian Heimes6f341092008-04-18 23:13:07 +0000265FUNC1(acos, acos, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000266 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000267FUNC1(acosh, acosh, 0,
268 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
269FUNC1(asin, asin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000270 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000271FUNC1(asinh, asinh, 0,
272 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
273FUNC1(atan, atan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000274 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
Mark Dickinson92483cd2008-04-20 21:39:04 +0000275FUNC2(atan2, m_atan2,
Tim Petersfe71f812001-08-07 22:10:00 +0000276 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
277 "Unlike atan(y/x), the signs of both x and y are considered.")
Christian Heimes6f341092008-04-18 23:13:07 +0000278FUNC1(atanh, atanh, 0,
279 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
280FUNC1(ceil, ceil, 0,
Jeffrey Yasskin9871d8f2008-01-05 08:47:13 +0000281 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
282 "This is the smallest integral value >= x.")
Christian Heimeseebb79c2008-01-03 22:32:26 +0000283FUNC2(copysign, copysign,
Christian Heimes6f341092008-04-18 23:13:07 +0000284 "copysign(x,y)\n\nReturn x with the sign of y.")
285FUNC1(cos, cos, 0,
286 "cos(x)\n\nReturn the cosine of x (measured in radians).")
287FUNC1(cosh, cosh, 1,
288 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
289FUNC1(exp, exp, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000290 "exp(x)\n\nReturn e raised to the power of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000291FUNC1(fabs, fabs, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000292 "fabs(x)\n\nReturn the absolute value of the float x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000293FUNC1(floor, floor, 0,
Jeffrey Yasskin9871d8f2008-01-05 08:47:13 +0000294 "floor(x)\n\nReturn the floor of x as a float.\n"
295 "This is the largest integral value <= x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000296FUNC1(log1p, log1p, 1,
297 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
298 The result is computed in a way which is accurate for x near zero.")
299FUNC1(sin, sin, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000300 "sin(x)\n\nReturn the sine of x (measured in radians).")
Christian Heimes6f341092008-04-18 23:13:07 +0000301FUNC1(sinh, sinh, 1,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000302 "sinh(x)\n\nReturn the hyperbolic sine of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000303FUNC1(sqrt, sqrt, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000304 "sqrt(x)\n\nReturn the square root of x.")
Christian Heimes6f341092008-04-18 23:13:07 +0000305FUNC1(tan, tan, 0,
Tim Petersfe71f812001-08-07 22:10:00 +0000306 "tan(x)\n\nReturn the tangent of x (measured in radians).")
Christian Heimes6f341092008-04-18 23:13:07 +0000307FUNC1(tanh, tanh, 0,
Guido van Rossumc6e22901998-12-04 19:26:43 +0000308 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000309
Barry Warsaw8b43b191996-12-09 22:32:36 +0000310static PyObject *
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000311math_trunc(PyObject *self, PyObject *number)
312{
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000313 return PyObject_CallMethod(number, "__trunc__", NULL);
314}
315
316PyDoc_STRVAR(math_trunc_doc,
317"trunc(x:Real) -> Integral\n"
318"\n"
Raymond Hettingerfe424f72008-02-02 05:24:44 +0000319"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000320
321static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000322math_frexp(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000323{
Guido van Rossumd18ad581991-10-24 14:57:21 +0000324 int i;
Neal Norwitz45e230a2006-11-19 21:26:53 +0000325 double x = PyFloat_AsDouble(arg);
326 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000327 return NULL;
Christian Heimes6f341092008-04-18 23:13:07 +0000328 /* deal with special cases directly, to sidestep platform
329 differences */
330 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
331 i = 0;
332 }
333 else {
334 PyFPE_START_PROTECT("in math_frexp", return 0);
335 x = frexp(x, &i);
336 PyFPE_END_PROTECT(x);
337 }
338 return Py_BuildValue("(di)", x, i);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000339}
340
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000341PyDoc_STRVAR(math_frexp_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000342"frexp(x)\n"
343"\n"
344"Return the mantissa and exponent of x, as pair (m, e).\n"
345"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 +0000346"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 +0000347
Barry Warsaw8b43b191996-12-09 22:32:36 +0000348static PyObject *
Fred Drake40c48682000-07-03 18:11:56 +0000349math_ldexp(PyObject *self, PyObject *args)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000350{
Christian Heimes6f341092008-04-18 23:13:07 +0000351 double x, r;
Guido van Rossumc5545052000-05-08 14:29:38 +0000352 int exp;
Fred Drake40c48682000-07-03 18:11:56 +0000353 if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
Guido van Rossumd18ad581991-10-24 14:57:21 +0000354 return NULL;
355 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000356 PyFPE_START_PROTECT("in math_ldexp", return 0)
357 r = ldexp(x, exp);
358 PyFPE_END_PROTECT(r)
359 if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
360 errno = ERANGE;
361 /* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
362 we unset it to avoid raising a ValueError here. */
363 if (errno == EDOM)
364 errno = 0;
365 if (errno && is_error(r))
Tim Peters1d120612000-10-12 06:10:25 +0000366 return NULL;
Guido van Rossumd18ad581991-10-24 14:57:21 +0000367 else
Christian Heimes6f341092008-04-18 23:13:07 +0000368 return PyFloat_FromDouble(r);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000369}
370
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000371PyDoc_STRVAR(math_ldexp_doc,
372"ldexp(x, i) -> x * (2**i)");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000373
Barry Warsaw8b43b191996-12-09 22:32:36 +0000374static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000375math_modf(PyObject *self, PyObject *arg)
Guido van Rossumd18ad581991-10-24 14:57:21 +0000376{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000377 double y, x = PyFloat_AsDouble(arg);
378 if (x == -1.0 && PyErr_Occurred())
Guido van Rossumd18ad581991-10-24 14:57:21 +0000379 return NULL;
Mark Dickinsonb2f70902008-04-20 01:39:24 +0000380 /* some platforms don't do the right thing for NaNs and
381 infinities, so we take care of special cases directly. */
382 if (!Py_IS_FINITE(x)) {
383 if (Py_IS_INFINITY(x))
384 return Py_BuildValue("(dd)", copysign(0., x), x);
385 else if (Py_IS_NAN(x))
386 return Py_BuildValue("(dd)", x, x);
387 }
388
Guido van Rossumd18ad581991-10-24 14:57:21 +0000389 errno = 0;
Christian Heimes6f341092008-04-18 23:13:07 +0000390 PyFPE_START_PROTECT("in math_modf", return 0);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000391 x = modf(x, &y);
Christian Heimes6f341092008-04-18 23:13:07 +0000392 PyFPE_END_PROTECT(x);
393 return Py_BuildValue("(dd)", x, y);
Guido van Rossumd18ad581991-10-24 14:57:21 +0000394}
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000395
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000396PyDoc_STRVAR(math_modf_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000397"modf(x)\n"
398"\n"
399"Return the fractional and integer parts of x. Both results carry the sign\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000400"of x. The integer part is returned as a real.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000401
Tim Peters78526162001-09-05 00:53:45 +0000402/* A decent logarithm is easy to compute even for huge longs, but libm can't
403 do that by itself -- loghelper can. func is log or log10, and name is
404 "log" or "log10". Note that overflow isn't possible: a long can contain
405 no more than INT_MAX * SHIFT bits, so has value certainly less than
406 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
407 small enough to fit in an IEEE single. log and log10 are even smaller.
408*/
409
410static PyObject*
Neal Norwitz45e230a2006-11-19 21:26:53 +0000411loghelper(PyObject* arg, double (*func)(double), char *funcname)
Tim Peters78526162001-09-05 00:53:45 +0000412{
Tim Peters78526162001-09-05 00:53:45 +0000413 /* If it is long, do it ourselves. */
414 if (PyLong_Check(arg)) {
415 double x;
416 int e;
417 x = _PyLong_AsScaledDouble(arg, &e);
418 if (x <= 0.0) {
419 PyErr_SetString(PyExc_ValueError,
420 "math domain error");
421 return NULL;
422 }
Christian Heimes543cabc2008-01-25 14:54:23 +0000423 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
424 log(x) + log(2) * e * PyLong_SHIFT.
425 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
Tim Peters78526162001-09-05 00:53:45 +0000426 so force use of double. */
Christian Heimes543cabc2008-01-25 14:54:23 +0000427 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
Tim Peters78526162001-09-05 00:53:45 +0000428 return PyFloat_FromDouble(x);
429 }
430
431 /* Else let libm handle it by itself. */
Christian Heimes6f341092008-04-18 23:13:07 +0000432 return math_1(arg, func, 0);
Tim Peters78526162001-09-05 00:53:45 +0000433}
434
435static PyObject *
436math_log(PyObject *self, PyObject *args)
437{
Raymond Hettinger866964c2002-12-14 19:51:34 +0000438 PyObject *arg;
439 PyObject *base = NULL;
440 PyObject *num, *den;
441 PyObject *ans;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000442
Raymond Hettingerea3fdf42002-12-29 16:33:45 +0000443 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
Raymond Hettinger866964c2002-12-14 19:51:34 +0000444 return NULL;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000445
Neal Norwitz45e230a2006-11-19 21:26:53 +0000446 num = loghelper(arg, log, "log");
447 if (num == NULL || base == NULL)
448 return num;
Raymond Hettinger866964c2002-12-14 19:51:34 +0000449
Neal Norwitz45e230a2006-11-19 21:26:53 +0000450 den = loghelper(base, log, "log");
Raymond Hettinger866964c2002-12-14 19:51:34 +0000451 if (den == NULL) {
452 Py_DECREF(num);
453 return NULL;
454 }
455
456 ans = PyNumber_Divide(num, den);
457 Py_DECREF(num);
458 Py_DECREF(den);
459 return ans;
Tim Peters78526162001-09-05 00:53:45 +0000460}
461
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000462PyDoc_STRVAR(math_log_doc,
Raymond Hettinger866964c2002-12-14 19:51:34 +0000463"log(x[, base]) -> the logarithm of x to the given base.\n\
464If the base not specified, returns the natural logarithm (base e) of x.");
Tim Peters78526162001-09-05 00:53:45 +0000465
466static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000467math_log10(PyObject *self, PyObject *arg)
Tim Peters78526162001-09-05 00:53:45 +0000468{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000469 return loghelper(arg, log10, "log10");
Tim Peters78526162001-09-05 00:53:45 +0000470}
471
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000472PyDoc_STRVAR(math_log10_doc,
473"log10(x) -> the base 10 logarithm of x.");
Tim Peters78526162001-09-05 00:53:45 +0000474
Christian Heimes6f341092008-04-18 23:13:07 +0000475static PyObject *
476math_fmod(PyObject *self, PyObject *args)
477{
478 PyObject *ox, *oy;
479 double r, x, y;
480 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
481 return NULL;
482 x = PyFloat_AsDouble(ox);
483 y = PyFloat_AsDouble(oy);
484 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
485 return NULL;
486 /* fmod(x, +/-Inf) returns x for finite x. */
487 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
488 return PyFloat_FromDouble(x);
489 errno = 0;
490 PyFPE_START_PROTECT("in math_fmod", return 0);
491 r = fmod(x, y);
492 PyFPE_END_PROTECT(r);
493 if (Py_IS_NAN(r)) {
494 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
495 errno = EDOM;
496 else
497 errno = 0;
498 }
499 if (errno && is_error(r))
500 return NULL;
501 else
502 return PyFloat_FromDouble(r);
503}
504
505PyDoc_STRVAR(math_fmod_doc,
506"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
507" x % y may differ.");
508
509static PyObject *
510math_hypot(PyObject *self, PyObject *args)
511{
512 PyObject *ox, *oy;
513 double r, x, y;
514 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
515 return NULL;
516 x = PyFloat_AsDouble(ox);
517 y = PyFloat_AsDouble(oy);
518 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
519 return NULL;
520 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
521 if (Py_IS_INFINITY(x))
522 return PyFloat_FromDouble(fabs(x));
523 if (Py_IS_INFINITY(y))
524 return PyFloat_FromDouble(fabs(y));
525 errno = 0;
526 PyFPE_START_PROTECT("in math_hypot", return 0);
527 r = hypot(x, y);
528 PyFPE_END_PROTECT(r);
529 if (Py_IS_NAN(r)) {
530 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
531 errno = EDOM;
532 else
533 errno = 0;
534 }
535 else if (Py_IS_INFINITY(r)) {
536 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
537 errno = ERANGE;
538 else
539 errno = 0;
540 }
541 if (errno && is_error(r))
542 return NULL;
543 else
544 return PyFloat_FromDouble(r);
545}
546
547PyDoc_STRVAR(math_hypot_doc,
548"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
549
550/* pow can't use math_2, but needs its own wrapper: the problem is
551 that an infinite result can arise either as a result of overflow
552 (in which case OverflowError should be raised) or as a result of
553 e.g. 0.**-5. (for which ValueError needs to be raised.)
554*/
555
556static PyObject *
557math_pow(PyObject *self, PyObject *args)
558{
559 PyObject *ox, *oy;
560 double r, x, y;
Mark Dickinsoncec3f132008-04-20 04:13:13 +0000561 int odd_y;
Christian Heimes6f341092008-04-18 23:13:07 +0000562
563 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
564 return NULL;
565 x = PyFloat_AsDouble(ox);
566 y = PyFloat_AsDouble(oy);
567 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
568 return NULL;
Mark Dickinsona1293eb2008-04-19 19:41:52 +0000569
Mark Dickinsoncec3f132008-04-20 04:13:13 +0000570 /* deal directly with IEEE specials, to cope with problems on various
571 platforms whose semantics don't exactly match C99 */
572 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
573 errno = 0;
574 if (Py_IS_NAN(x))
575 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
576 else if (Py_IS_NAN(y))
577 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
578 else if (Py_IS_INFINITY(x)) {
579 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
580 if (y > 0.)
581 r = odd_y ? x : fabs(x);
582 else if (y == 0.)
583 r = 1.;
584 else /* y < 0. */
585 r = odd_y ? copysign(0., x) : 0.;
586 }
587 else if (Py_IS_INFINITY(y)) {
588 if (fabs(x) == 1.0)
589 r = 1.;
590 else if (y > 0. && fabs(x) > 1.0)
591 r = y;
592 else if (y < 0. && fabs(x) < 1.0) {
593 r = -y; /* result is +inf */
594 if (x == 0.) /* 0**-inf: divide-by-zero */
595 errno = EDOM;
596 }
597 else
598 r = 0.;
599 }
Mark Dickinsone941d972008-04-19 18:51:48 +0000600 }
Mark Dickinsoncec3f132008-04-20 04:13:13 +0000601 else {
602 /* let libm handle finite**finite */
603 errno = 0;
604 PyFPE_START_PROTECT("in math_pow", return 0);
605 r = pow(x, y);
606 PyFPE_END_PROTECT(r);
607 /* a NaN result should arise only from (-ve)**(finite
608 non-integer); in this case we want to raise ValueError. */
609 if (!Py_IS_FINITE(r)) {
610 if (Py_IS_NAN(r)) {
611 errno = EDOM;
612 }
613 /*
614 an infinite result here arises either from:
615 (A) (+/-0.)**negative (-> divide-by-zero)
616 (B) overflow of x**y with x and y finite
617 */
618 else if (Py_IS_INFINITY(r)) {
619 if (x == 0.)
620 errno = EDOM;
621 else
622 errno = ERANGE;
623 }
624 }
Christian Heimes6f341092008-04-18 23:13:07 +0000625 }
626
627 if (errno && is_error(r))
628 return NULL;
629 else
630 return PyFloat_FromDouble(r);
631}
632
633PyDoc_STRVAR(math_pow_doc,
634"pow(x,y)\n\nReturn x**y (x to the power of y).");
635
Christian Heimese2ca4242008-01-03 20:23:15 +0000636static const double degToRad = Py_MATH_PI / 180.0;
637static const double radToDeg = 180.0 / Py_MATH_PI;
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000638
639static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000640math_degrees(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000641{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000642 double x = PyFloat_AsDouble(arg);
643 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000644 return NULL;
Christian Heimese2ca4242008-01-03 20:23:15 +0000645 return PyFloat_FromDouble(x * radToDeg);
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000646}
647
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000648PyDoc_STRVAR(math_degrees_doc,
649"degrees(x) -> converts angle x from radians to degrees");
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000650
651static PyObject *
Neal Norwitz45e230a2006-11-19 21:26:53 +0000652math_radians(PyObject *self, PyObject *arg)
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000653{
Neal Norwitz45e230a2006-11-19 21:26:53 +0000654 double x = PyFloat_AsDouble(arg);
655 if (x == -1.0 && PyErr_Occurred())
Raymond Hettingerd6f22672002-05-13 03:56:10 +0000656 return NULL;
657 return PyFloat_FromDouble(x * degToRad);
658}
659
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000660PyDoc_STRVAR(math_radians_doc,
661"radians(x) -> converts angle x from degrees to radians");
Tim Peters78526162001-09-05 00:53:45 +0000662
Christian Heimese2ca4242008-01-03 20:23:15 +0000663static PyObject *
664math_isnan(PyObject *self, PyObject *arg)
665{
666 double x = PyFloat_AsDouble(arg);
667 if (x == -1.0 && PyErr_Occurred())
668 return NULL;
669 return PyBool_FromLong((long)Py_IS_NAN(x));
670}
671
672PyDoc_STRVAR(math_isnan_doc,
673"isnan(x) -> bool\n\
674Checks if float x is not a number (NaN)");
675
676static PyObject *
677math_isinf(PyObject *self, PyObject *arg)
678{
679 double x = PyFloat_AsDouble(arg);
680 if (x == -1.0 && PyErr_Occurred())
681 return NULL;
682 return PyBool_FromLong((long)Py_IS_INFINITY(x));
683}
684
685PyDoc_STRVAR(math_isinf_doc,
686"isinf(x) -> bool\n\
687Checks if float x is infinite (positive or negative)");
688
Barry Warsaw8b43b191996-12-09 22:32:36 +0000689static PyMethodDef math_methods[] = {
Neal Norwitz45e230a2006-11-19 21:26:53 +0000690 {"acos", math_acos, METH_O, math_acos_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000691 {"acosh", math_acosh, METH_O, math_acosh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000692 {"asin", math_asin, METH_O, math_asin_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000693 {"asinh", math_asinh, METH_O, math_asinh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000694 {"atan", math_atan, METH_O, math_atan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000695 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000696 {"atanh", math_atanh, METH_O, math_atanh_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000697 {"ceil", math_ceil, METH_O, math_ceil_doc},
Christian Heimeseebb79c2008-01-03 22:32:26 +0000698 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000699 {"cos", math_cos, METH_O, math_cos_doc},
700 {"cosh", math_cosh, METH_O, math_cosh_doc},
701 {"degrees", math_degrees, METH_O, math_degrees_doc},
702 {"exp", math_exp, METH_O, math_exp_doc},
703 {"fabs", math_fabs, METH_O, math_fabs_doc},
704 {"floor", math_floor, METH_O, math_floor_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000705 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000706 {"frexp", math_frexp, METH_O, math_frexp_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000707 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
Christian Heimese2ca4242008-01-03 20:23:15 +0000708 {"isinf", math_isinf, METH_O, math_isinf_doc},
709 {"isnan", math_isnan, METH_O, math_isnan_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000710 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
711 {"log", math_log, METH_VARARGS, math_log_doc},
Christian Heimes6f341092008-04-18 23:13:07 +0000712 {"log1p", math_log1p, METH_O, math_log1p_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000713 {"log10", math_log10, METH_O, math_log10_doc},
714 {"modf", math_modf, METH_O, math_modf_doc},
Fred Drake40c48682000-07-03 18:11:56 +0000715 {"pow", math_pow, METH_VARARGS, math_pow_doc},
Neal Norwitz45e230a2006-11-19 21:26:53 +0000716 {"radians", math_radians, METH_O, math_radians_doc},
717 {"sin", math_sin, METH_O, math_sin_doc},
718 {"sinh", math_sinh, METH_O, math_sinh_doc},
719 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
720 {"tan", math_tan, METH_O, math_tan_doc},
721 {"tanh", math_tanh, METH_O, math_tanh_doc},
Jeffrey Yasskinca2b69f2008-02-01 06:22:46 +0000722 {"trunc", math_trunc, METH_O, math_trunc_doc},
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000723 {NULL, NULL} /* sentinel */
724};
725
Guido van Rossumc6e22901998-12-04 19:26:43 +0000726
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000727PyDoc_STRVAR(module_doc,
Tim Peters63c94532001-09-04 23:17:42 +0000728"This module is always available. It provides access to the\n"
Martin v. Löwis14f8b4c2002-06-13 20:33:02 +0000729"mathematical functions defined by the C standard.");
Guido van Rossumc6e22901998-12-04 19:26:43 +0000730
Mark Hammondfe51c6d2002-08-02 02:27:13 +0000731PyMODINIT_FUNC
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000732initmath(void)
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000733{
Christian Heimes6f341092008-04-18 23:13:07 +0000734 PyObject *m;
Tim Petersfe71f812001-08-07 22:10:00 +0000735
Guido van Rossumc6e22901998-12-04 19:26:43 +0000736 m = Py_InitModule3("math", math_methods, module_doc);
Neal Norwitz1ac754f2006-01-19 06:09:39 +0000737 if (m == NULL)
738 goto finally;
Barry Warsawfc93f751996-12-17 00:47:03 +0000739
Christian Heimes6f341092008-04-18 23:13:07 +0000740 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
741 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
Barry Warsawfc93f751996-12-17 00:47:03 +0000742
Christian Heimes6f341092008-04-18 23:13:07 +0000743 finally:
Barry Warsaw9bfd2bf2000-09-01 09:01:32 +0000744 return;
Guido van Rossum85a5fbb1990-10-14 12:07:46 +0000745}