1 /* Math module -- standard C math library functions, pi and e */
3 /* 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
9 These are the "spirit of 754" rules:
11 1. If the mathematical result is a real number, but of magnitude too
12 large to approximate by a machine float, overflow is signaled and the
13 result is an infinity (with the appropriate sign).
15 2. If the mathematical result is a real number, but of magnitude too
16 small to approximate by a machine float, underflow is signaled and the
17 result is a zero (with the appropriate sign).
19 3. At a singularity (a value x such that the limit of f(y) as y
20 approaches x exists and is an infinity), "divide by zero" is signaled
21 and the result is an infinity (with the appropriate sign). This is
22 complicated a little by that the left-side and right-side limits may
23 not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
24 from the positive or negative directions. In that specific case, the
25 sign of the zero determines the result of 1/0.
27 4. At a point where a function has no defined result in the extended
28 reals (i.e., the reals plus an infinity or two), invalid operation is
29 signaled and a NaN is returned.
31 And these are what Python has historically /tried/ to do (but not
32 always successfully, as platform libm behavior varies a lot):
34 For #1, raise OverflowError.
36 For #2, return a zero (with the appropriate sign if that happens by
39 For #3 and #4, raise ValueError. It may have made sense to raise
40 Python's ZeroDivisionError in #3, but historically that's only been
41 raised for division by zero and mod by zero.
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
56 #include "longintrepr.h" /* just for SHIFT */
59 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
60 extern double copysign(double, double);
63 /* 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.
70 int result
= 1; /* presumption of guilt */
71 assert(errno
); /* non-zero errno is a precondition for calling */
73 PyErr_SetString(PyExc_ValueError
, "math domain error");
75 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
79 * 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).
86 * On some platforms (Ubuntu/ia64) it seems that errno can be
87 * set to ERANGE for subnormal results that do *not* underflow
88 * to zero. So to be safe, we'll ignore ERANGE whenever the
89 * function result is less than one in absolute value.
94 PyErr_SetString(PyExc_OverflowError
,
98 /* Unexpected math error */
99 PyErr_SetFromErrno(PyExc_ValueError
);
104 wrapper for atan2 that deals directly with special cases before
105 delegating to the platform libm for the remaining cases. This
106 is necessary to get consistent behaviour across platforms.
107 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
112 m_atan2(double y
, double x
)
114 if (Py_IS_NAN(x
) || Py_IS_NAN(y
))
116 if (Py_IS_INFINITY(y
)) {
117 if (Py_IS_INFINITY(x
)) {
118 if (copysign(1., x
) == 1.)
119 /* atan2(+-inf, +inf) == +-pi/4 */
120 return copysign(0.25*Py_MATH_PI
, y
);
122 /* atan2(+-inf, -inf) == +-pi*3/4 */
123 return copysign(0.75*Py_MATH_PI
, y
);
125 /* atan2(+-inf, x) == +-pi/2 for finite x */
126 return copysign(0.5*Py_MATH_PI
, y
);
128 if (Py_IS_INFINITY(x
) || y
== 0.) {
129 if (copysign(1., x
) == 1.)
130 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
131 return copysign(0., y
);
133 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
134 return copysign(Py_MATH_PI
, y
);
140 math_1 is used to wrap a libm function f that takes a double
141 arguments and returns a double.
143 The error reporting follows these rules, which are designed to do
144 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
147 - a NaN result from non-NaN inputs causes ValueError to be raised
148 - an infinite result from finite inputs causes OverflowError to be
149 raised if can_overflow is 1, or raises ValueError if can_overflow
151 - if the result is finite and errno == EDOM then ValueError is
153 - if the result is finite and nonzero and errno == ERANGE then
154 OverflowError is raised
156 The last rule is used to catch overflow on platforms which follow
157 C89 but for which HUGE_VAL is not an infinity.
159 For the majority of one-argument functions these rules are enough
160 to ensure that Python's functions behave as specified in 'Annex F'
161 of the C99 standard, with the 'invalid' and 'divide-by-zero'
162 floating-point exceptions mapping to Python's ValueError and the
163 'overflow' floating-point exception mapping to OverflowError.
164 math_1 only works for functions that don't have singularities *and*
165 the possibility of overflow; fortunately, that covers everything we
166 care about right now.
170 math_1(PyObject
*arg
, double (*func
) (double), int can_overflow
)
173 x
= PyFloat_AsDouble(arg
);
174 if (x
== -1.0 && PyErr_Occurred())
177 PyFPE_START_PROTECT("in math_1", return 0);
179 PyFPE_END_PROTECT(r
);
186 else if (Py_IS_INFINITY(r
)) {
188 errno
= can_overflow
? ERANGE
: EDOM
;
192 if (errno
&& is_error(r
))
195 return PyFloat_FromDouble(r
);
199 math_2 is used to wrap a libm function f that takes two double
200 arguments and returns a double.
202 The error reporting follows these rules, which are designed to do
203 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
206 - a NaN result from non-NaN inputs causes ValueError to be raised
207 - an infinite result from finite inputs causes OverflowError to be
209 - if the result is finite and errno == EDOM then ValueError is
211 - if the result is finite and nonzero and errno == ERANGE then
212 OverflowError is raised
214 The last rule is used to catch overflow on platforms which follow
215 C89 but for which HUGE_VAL is not an infinity.
217 For most two-argument functions (copysign, fmod, hypot, atan2)
218 these rules are enough to ensure that Python's functions behave as
219 specified in 'Annex F' of the C99 standard, with the 'invalid' and
220 'divide-by-zero' floating-point exceptions mapping to Python's
221 ValueError and the 'overflow' floating-point exception mapping to
226 math_2(PyObject
*args
, double (*func
) (double, double), char *funcname
)
230 if (! PyArg_UnpackTuple(args
, funcname
, 2, 2, &ox
, &oy
))
232 x
= PyFloat_AsDouble(ox
);
233 y
= PyFloat_AsDouble(oy
);
234 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
237 PyFPE_START_PROTECT("in math_2", return 0);
239 PyFPE_END_PROTECT(r
);
241 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
246 else if (Py_IS_INFINITY(r
)) {
247 if (Py_IS_FINITE(x
) && Py_IS_FINITE(y
))
252 if (errno
&& is_error(r
))
255 return PyFloat_FromDouble(r
);
258 #define FUNC1(funcname, func, can_overflow, docstring) \
259 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
260 return math_1(args, func, can_overflow); \
262 PyDoc_STRVAR(math_##funcname##_doc, docstring);
264 #define FUNC2(funcname, func, docstring) \
265 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
266 return math_2(args, func, #funcname); \
268 PyDoc_STRVAR(math_##funcname##_doc, docstring);
271 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
272 FUNC1(acosh
, acosh
, 0,
273 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
275 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
276 FUNC1(asinh
, asinh
, 0,
277 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
279 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
280 FUNC2(atan2
, m_atan2
,
281 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
282 "Unlike atan(y/x), the signs of both x and y are considered.")
283 FUNC1(atanh
, atanh
, 0,
284 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
286 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
287 "This is the smallest integral value >= x.")
288 FUNC2(copysign
, copysign
,
289 "copysign(x,y)\n\nReturn x with the sign of y.")
291 "cos(x)\n\nReturn the cosine of x (measured in radians).")
293 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
295 "exp(x)\n\nReturn e raised to the power of x.")
297 "fabs(x)\n\nReturn the absolute value of the float x.")
298 FUNC1(floor
, floor
, 0,
299 "floor(x)\n\nReturn the floor of x as a float.\n"
300 "This is the largest integral value <= x.")
301 FUNC1(log1p
, log1p
, 1,
302 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
303 The result is computed in a way which is accurate for x near zero.")
305 "sin(x)\n\nReturn the sine of x (measured in radians).")
307 "sinh(x)\n\nReturn the hyperbolic sine of x.")
309 "sqrt(x)\n\nReturn the square root of x.")
311 "tan(x)\n\nReturn the tangent of x (measured in radians).")
313 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
315 /* Precision summation function as msum() by Raymond Hettinger in
316 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
317 enhanced with the exact partials sum and roundoff from Mark
318 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
319 See those links for more details, proofs and other references.
321 Note 1: IEEE 754R floating point semantics are assumed,
322 but the current implementation does not re-establish special
323 value semantics across iterations (i.e. handling -Inf + Inf).
325 Note 2: No provision is made for intermediate overflow handling;
326 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
327 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
328 overflow of the first partial sum.
330 Note 3: The intermediate values lo, yr, and hi are declared volatile so
331 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
332 Also, the volatile declaration forces the values to be stored in memory as
333 regular doubles instead of extended long precision (80-bit) values. This
334 prevents double rounding because any addition or subtraction of two doubles
335 can be resolved exactly into double-sized hi and lo values. As long as the
336 hi value gets forced into a double before yr and lo are computed, the extra
337 bits in downstream extended precision operations (x87 for example) will be
338 exactly zero and therefore can be losslessly stored back into a double,
339 thereby preventing double rounding.
341 Note 4: A similar implementation is in Modules/cmathmodule.c.
342 Be sure to update both when making changes.
344 Note 5: The signature of math.fsum() differs from __builtin__.sum()
345 because the start argument doesn't make sense in the context of
346 accurate summation. Since the partials table is collapsed before
347 returning a result, sum(seq2, start=sum(seq1)) may not equal the
348 accurate result returned by sum(itertools.chain(seq1, seq2)).
351 #define NUM_PARTIALS 32 /* initial partials array size, on stack */
353 /* Extend the partials array p[] by doubling its size. */
354 static int /* non-zero on error */
355 _fsum_realloc(double **p_ptr
, Py_ssize_t n
,
356 double *ps
, Py_ssize_t
*m_ptr
)
359 Py_ssize_t m
= *m_ptr
;
362 if (n
< m
&& m
< (PY_SSIZE_T_MAX
/ sizeof(double))) {
365 v
= PyMem_Malloc(sizeof(double) * m
);
367 memcpy(v
, ps
, sizeof(double) * n
);
370 v
= PyMem_Realloc(p
, sizeof(double) * m
);
372 if (v
== NULL
) { /* size overflow or no memory */
373 PyErr_SetString(PyExc_MemoryError
, "math.fsum partials");
376 *p_ptr
= (double*) v
;
381 /* Full precision summation of a sequence of floats.
384 partials = [] # sorted, non-overlapping partial sums
397 return sum_exact(partials)
399 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
400 are exactly equal to x+y. The inner loop applies hi/lo summation to each
401 partial so that the list of partial sums remains exact.
403 Sum_exact() adds the partial sums exactly and correctly rounds the final
404 result (using the round-half-to-even rule). The items in partials remain
405 non-zero, non-special, non-overlapping and strictly increasing in
406 magnitude, but possibly not all having the same sign.
408 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
412 math_fsum(PyObject
*self
, PyObject
*seq
)
414 PyObject
*item
, *iter
, *sum
= NULL
;
415 Py_ssize_t i
, j
, n
= 0, m
= NUM_PARTIALS
;
416 double x
, y
, t
, ps
[NUM_PARTIALS
], *p
= ps
;
417 double xsave
, special_sum
= 0.0, inf_sum
= 0.0;
418 volatile double hi
, yr
, lo
;
420 iter
= PyObject_GetIter(seq
);
424 PyFPE_START_PROTECT("fsum", Py_DECREF(iter
); return NULL
)
426 for(;;) { /* for x in iterable */
427 assert(0 <= n
&& n
<= m
);
428 assert((m
== NUM_PARTIALS
&& p
== ps
) ||
429 (m
> NUM_PARTIALS
&& p
!= NULL
));
431 item
= PyIter_Next(iter
);
433 if (PyErr_Occurred())
437 x
= PyFloat_AsDouble(item
);
439 if (PyErr_Occurred())
443 for (i
= j
= 0; j
< n
; j
++) { /* for y in partials */
445 if (fabs(x
) < fabs(y
)) {
456 n
= i
; /* ps[i:] = [x] */
458 if (! Py_IS_FINITE(x
)) {
459 /* a nonfinite x could arise either as
460 a result of intermediate overflow, or
461 as a result of a nan or inf in the
463 if (Py_IS_FINITE(xsave
)) {
464 PyErr_SetString(PyExc_OverflowError
,
465 "intermediate overflow in fsum");
468 if (Py_IS_INFINITY(xsave
))
470 special_sum
+= xsave
;
474 else if (n
>= m
&& _fsum_realloc(&p
, n
, ps
, &m
))
481 if (special_sum
!= 0.0) {
482 if (Py_IS_NAN(inf_sum
))
483 PyErr_SetString(PyExc_ValueError
,
484 "-inf + inf in fsum");
486 sum
= PyFloat_FromDouble(special_sum
);
493 /* sum_exact(ps, hi) from the top, stop when the sum becomes
498 assert(fabs(y
) < fabs(x
));
505 /* Make half-even rounding work across multiple partials.
506 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
507 digit to two instead of down to zero (the 1e-16 makes the 1
508 slightly closer to two). With a potential 1 ULP rounding
509 error fixed-up, math.fsum() can guarantee commutativity. */
510 if (n
> 0 && ((lo
< 0.0 && p
[n
-1] < 0.0) ||
511 (lo
> 0.0 && p
[n
-1] > 0.0))) {
519 sum
= PyFloat_FromDouble(hi
);
522 PyFPE_END_PROTECT(hi
)
531 PyDoc_STRVAR(math_fsum_doc
,
533 Return an accurate floating point sum of values in the iterable.\n\
534 Assumes IEEE-754 floating point arithmetic.");
537 math_factorial(PyObject
*self
, PyObject
*arg
)
540 PyObject
*result
, *iobj
, *newresult
;
542 if (PyFloat_Check(arg
)) {
543 double dx
= PyFloat_AS_DOUBLE((PyFloatObject
*)arg
);
544 if (dx
!= floor(dx
)) {
545 PyErr_SetString(PyExc_ValueError
,
546 "factorial() only accepts integral values");
551 x
= PyInt_AsLong(arg
);
552 if (x
== -1 && PyErr_Occurred())
555 PyErr_SetString(PyExc_ValueError
,
556 "factorial() not defined for negative values");
560 result
= (PyObject
*)PyInt_FromLong(1);
563 for (i
=1 ; i
<=x
; i
++) {
564 iobj
= (PyObject
*)PyInt_FromLong(i
);
567 newresult
= PyNumber_Multiply(result
, iobj
);
569 if (newresult
== NULL
)
581 PyDoc_STRVAR(math_factorial_doc
, "Return n!");
584 math_trunc(PyObject
*self
, PyObject
*number
)
586 return PyObject_CallMethod(number
, "__trunc__", NULL
);
589 PyDoc_STRVAR(math_trunc_doc
,
590 "trunc(x:Real) -> Integral\n"
592 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
595 math_frexp(PyObject
*self
, PyObject
*arg
)
598 double x
= PyFloat_AsDouble(arg
);
599 if (x
== -1.0 && PyErr_Occurred())
601 /* deal with special cases directly, to sidestep platform
603 if (Py_IS_NAN(x
) || Py_IS_INFINITY(x
) || !x
) {
607 PyFPE_START_PROTECT("in math_frexp", return 0);
609 PyFPE_END_PROTECT(x
);
611 return Py_BuildValue("(di)", x
, i
);
614 PyDoc_STRVAR(math_frexp_doc
,
617 "Return the mantissa and exponent of x, as pair (m, e).\n"
618 "m is a float and e is an int, such that x = m * 2.**e.\n"
619 "If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
622 math_ldexp(PyObject
*self
, PyObject
*args
)
627 if (! PyArg_ParseTuple(args
, "dO:ldexp", &x
, &oexp
))
630 if (PyLong_Check(oexp
)) {
631 /* on overflow, replace exponent with either LONG_MAX
632 or LONG_MIN, depending on the sign. */
633 exp
= PyLong_AsLong(oexp
);
634 if (exp
== -1 && PyErr_Occurred()) {
635 if (PyErr_ExceptionMatches(PyExc_OverflowError
)) {
636 if (Py_SIZE(oexp
) < 0) {
645 /* propagate any unexpected exception */
650 else if (PyInt_Check(oexp
)) {
651 exp
= PyInt_AS_LONG(oexp
);
654 PyErr_SetString(PyExc_TypeError
,
655 "Expected an int or long as second argument "
660 if (x
== 0. || !Py_IS_FINITE(x
)) {
661 /* NaNs, zeros and infinities are returned unchanged */
664 } else if (exp
> INT_MAX
) {
666 r
= copysign(Py_HUGE_VAL
, x
);
668 } else if (exp
< INT_MIN
) {
669 /* underflow to +-0 */
674 PyFPE_START_PROTECT("in math_ldexp", return 0);
675 r
= ldexp(x
, (int)exp
);
676 PyFPE_END_PROTECT(r
);
677 if (Py_IS_INFINITY(r
))
681 if (errno
&& is_error(r
))
683 return PyFloat_FromDouble(r
);
686 PyDoc_STRVAR(math_ldexp_doc
,
687 "ldexp(x, i) -> x * (2**i)");
690 math_modf(PyObject
*self
, PyObject
*arg
)
692 double y
, x
= PyFloat_AsDouble(arg
);
693 if (x
== -1.0 && PyErr_Occurred())
695 /* some platforms don't do the right thing for NaNs and
696 infinities, so we take care of special cases directly. */
697 if (!Py_IS_FINITE(x
)) {
698 if (Py_IS_INFINITY(x
))
699 return Py_BuildValue("(dd)", copysign(0., x
), x
);
700 else if (Py_IS_NAN(x
))
701 return Py_BuildValue("(dd)", x
, x
);
705 PyFPE_START_PROTECT("in math_modf", return 0);
707 PyFPE_END_PROTECT(x
);
708 return Py_BuildValue("(dd)", x
, y
);
711 PyDoc_STRVAR(math_modf_doc
,
714 "Return the fractional and integer parts of x. Both results carry the sign\n"
715 "of x. The integer part is returned as a real.");
717 /* A decent logarithm is easy to compute even for huge longs, but libm can't
718 do that by itself -- loghelper can. func is log or log10, and name is
719 "log" or "log10". Note that overflow isn't possible: a long can contain
720 no more than INT_MAX * SHIFT bits, so has value certainly less than
721 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
722 small enough to fit in an IEEE single. log and log10 are even smaller.
726 loghelper(PyObject
* arg
, double (*func
)(double), char *funcname
)
728 /* If it is long, do it ourselves. */
729 if (PyLong_Check(arg
)) {
732 x
= _PyLong_AsScaledDouble(arg
, &e
);
734 PyErr_SetString(PyExc_ValueError
,
735 "math domain error");
738 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
739 log(x) + log(2) * e * PyLong_SHIFT.
740 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
741 so force use of double. */
742 x
= func(x
) + (e
* (double)PyLong_SHIFT
) * func(2.0);
743 return PyFloat_FromDouble(x
);
746 /* Else let libm handle it by itself. */
747 return math_1(arg
, func
, 0);
751 math_log(PyObject
*self
, PyObject
*args
)
754 PyObject
*base
= NULL
;
758 if (!PyArg_UnpackTuple(args
, "log", 1, 2, &arg
, &base
))
761 num
= loghelper(arg
, log
, "log");
762 if (num
== NULL
|| base
== NULL
)
765 den
= loghelper(base
, log
, "log");
771 ans
= PyNumber_Divide(num
, den
);
777 PyDoc_STRVAR(math_log_doc
,
778 "log(x[, base]) -> the logarithm of x to the given base.\n\
779 If the base not specified, returns the natural logarithm (base e) of x.");
782 math_log10(PyObject
*self
, PyObject
*arg
)
784 return loghelper(arg
, log10
, "log10");
787 PyDoc_STRVAR(math_log10_doc
,
788 "log10(x) -> the base 10 logarithm of x.");
791 math_fmod(PyObject
*self
, PyObject
*args
)
795 if (! PyArg_UnpackTuple(args
, "fmod", 2, 2, &ox
, &oy
))
797 x
= PyFloat_AsDouble(ox
);
798 y
= PyFloat_AsDouble(oy
);
799 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
801 /* fmod(x, +/-Inf) returns x for finite x. */
802 if (Py_IS_INFINITY(y
) && Py_IS_FINITE(x
))
803 return PyFloat_FromDouble(x
);
805 PyFPE_START_PROTECT("in math_fmod", return 0);
807 PyFPE_END_PROTECT(r
);
809 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
814 if (errno
&& is_error(r
))
817 return PyFloat_FromDouble(r
);
820 PyDoc_STRVAR(math_fmod_doc
,
821 "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
822 " x % y may differ.");
825 math_hypot(PyObject
*self
, PyObject
*args
)
829 if (! PyArg_UnpackTuple(args
, "hypot", 2, 2, &ox
, &oy
))
831 x
= PyFloat_AsDouble(ox
);
832 y
= PyFloat_AsDouble(oy
);
833 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
835 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
836 if (Py_IS_INFINITY(x
))
837 return PyFloat_FromDouble(fabs(x
));
838 if (Py_IS_INFINITY(y
))
839 return PyFloat_FromDouble(fabs(y
));
841 PyFPE_START_PROTECT("in math_hypot", return 0);
843 PyFPE_END_PROTECT(r
);
845 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
850 else if (Py_IS_INFINITY(r
)) {
851 if (Py_IS_FINITE(x
) && Py_IS_FINITE(y
))
856 if (errno
&& is_error(r
))
859 return PyFloat_FromDouble(r
);
862 PyDoc_STRVAR(math_hypot_doc
,
863 "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
865 /* pow can't use math_2, but needs its own wrapper: the problem is
866 that an infinite result can arise either as a result of overflow
867 (in which case OverflowError should be raised) or as a result of
868 e.g. 0.**-5. (for which ValueError needs to be raised.)
872 math_pow(PyObject
*self
, PyObject
*args
)
878 if (! PyArg_UnpackTuple(args
, "pow", 2, 2, &ox
, &oy
))
880 x
= PyFloat_AsDouble(ox
);
881 y
= PyFloat_AsDouble(oy
);
882 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
885 /* deal directly with IEEE specials, to cope with problems on various
886 platforms whose semantics don't exactly match C99 */
887 r
= 0.; /* silence compiler warning */
888 if (!Py_IS_FINITE(x
) || !Py_IS_FINITE(y
)) {
891 r
= y
== 0. ? 1. : x
; /* NaN**0 = 1 */
892 else if (Py_IS_NAN(y
))
893 r
= x
== 1. ? 1. : y
; /* 1**NaN = 1 */
894 else if (Py_IS_INFINITY(x
)) {
895 odd_y
= Py_IS_FINITE(y
) && fmod(fabs(y
), 2.0) == 1.0;
897 r
= odd_y
? x
: fabs(x
);
901 r
= odd_y
? copysign(0., x
) : 0.;
903 else if (Py_IS_INFINITY(y
)) {
906 else if (y
> 0. && fabs(x
) > 1.0)
908 else if (y
< 0. && fabs(x
) < 1.0) {
909 r
= -y
; /* result is +inf */
910 if (x
== 0.) /* 0**-inf: divide-by-zero */
918 /* let libm handle finite**finite */
920 PyFPE_START_PROTECT("in math_pow", return 0);
922 PyFPE_END_PROTECT(r
);
923 /* a NaN result should arise only from (-ve)**(finite
924 non-integer); in this case we want to raise ValueError. */
925 if (!Py_IS_FINITE(r
)) {
930 an infinite result here arises either from:
931 (A) (+/-0.)**negative (-> divide-by-zero)
932 (B) overflow of x**y with x and y finite
934 else if (Py_IS_INFINITY(r
)) {
943 if (errno
&& is_error(r
))
946 return PyFloat_FromDouble(r
);
949 PyDoc_STRVAR(math_pow_doc
,
950 "pow(x,y)\n\nReturn x**y (x to the power of y).");
952 static const double degToRad
= Py_MATH_PI
/ 180.0;
953 static const double radToDeg
= 180.0 / Py_MATH_PI
;
956 math_degrees(PyObject
*self
, PyObject
*arg
)
958 double x
= PyFloat_AsDouble(arg
);
959 if (x
== -1.0 && PyErr_Occurred())
961 return PyFloat_FromDouble(x
* radToDeg
);
964 PyDoc_STRVAR(math_degrees_doc
,
965 "degrees(x) -> converts angle x from radians to degrees");
968 math_radians(PyObject
*self
, PyObject
*arg
)
970 double x
= PyFloat_AsDouble(arg
);
971 if (x
== -1.0 && PyErr_Occurred())
973 return PyFloat_FromDouble(x
* degToRad
);
976 PyDoc_STRVAR(math_radians_doc
,
977 "radians(x) -> converts angle x from degrees to radians");
980 math_isnan(PyObject
*self
, PyObject
*arg
)
982 double x
= PyFloat_AsDouble(arg
);
983 if (x
== -1.0 && PyErr_Occurred())
985 return PyBool_FromLong((long)Py_IS_NAN(x
));
988 PyDoc_STRVAR(math_isnan_doc
,
990 Checks if float x is not a number (NaN)");
993 math_isinf(PyObject
*self
, PyObject
*arg
)
995 double x
= PyFloat_AsDouble(arg
);
996 if (x
== -1.0 && PyErr_Occurred())
998 return PyBool_FromLong((long)Py_IS_INFINITY(x
));
1001 PyDoc_STRVAR(math_isinf_doc
,
1002 "isinf(x) -> bool\n\
1003 Checks if float x is infinite (positive or negative)");
1005 static PyMethodDef math_methods
[] = {
1006 {"acos", math_acos
, METH_O
, math_acos_doc
},
1007 {"acosh", math_acosh
, METH_O
, math_acosh_doc
},
1008 {"asin", math_asin
, METH_O
, math_asin_doc
},
1009 {"asinh", math_asinh
, METH_O
, math_asinh_doc
},
1010 {"atan", math_atan
, METH_O
, math_atan_doc
},
1011 {"atan2", math_atan2
, METH_VARARGS
, math_atan2_doc
},
1012 {"atanh", math_atanh
, METH_O
, math_atanh_doc
},
1013 {"ceil", math_ceil
, METH_O
, math_ceil_doc
},
1014 {"copysign", math_copysign
, METH_VARARGS
, math_copysign_doc
},
1015 {"cos", math_cos
, METH_O
, math_cos_doc
},
1016 {"cosh", math_cosh
, METH_O
, math_cosh_doc
},
1017 {"degrees", math_degrees
, METH_O
, math_degrees_doc
},
1018 {"exp", math_exp
, METH_O
, math_exp_doc
},
1019 {"fabs", math_fabs
, METH_O
, math_fabs_doc
},
1020 {"factorial", math_factorial
, METH_O
, math_factorial_doc
},
1021 {"floor", math_floor
, METH_O
, math_floor_doc
},
1022 {"fmod", math_fmod
, METH_VARARGS
, math_fmod_doc
},
1023 {"frexp", math_frexp
, METH_O
, math_frexp_doc
},
1024 {"fsum", math_fsum
, METH_O
, math_fsum_doc
},
1025 {"hypot", math_hypot
, METH_VARARGS
, math_hypot_doc
},
1026 {"isinf", math_isinf
, METH_O
, math_isinf_doc
},
1027 {"isnan", math_isnan
, METH_O
, math_isnan_doc
},
1028 {"ldexp", math_ldexp
, METH_VARARGS
, math_ldexp_doc
},
1029 {"log", math_log
, METH_VARARGS
, math_log_doc
},
1030 {"log1p", math_log1p
, METH_O
, math_log1p_doc
},
1031 {"log10", math_log10
, METH_O
, math_log10_doc
},
1032 {"modf", math_modf
, METH_O
, math_modf_doc
},
1033 {"pow", math_pow
, METH_VARARGS
, math_pow_doc
},
1034 {"radians", math_radians
, METH_O
, math_radians_doc
},
1035 {"sin", math_sin
, METH_O
, math_sin_doc
},
1036 {"sinh", math_sinh
, METH_O
, math_sinh_doc
},
1037 {"sqrt", math_sqrt
, METH_O
, math_sqrt_doc
},
1038 {"tan", math_tan
, METH_O
, math_tan_doc
},
1039 {"tanh", math_tanh
, METH_O
, math_tanh_doc
},
1040 {"trunc", math_trunc
, METH_O
, math_trunc_doc
},
1041 {NULL
, NULL
} /* sentinel */
1045 PyDoc_STRVAR(module_doc
,
1046 "This module is always available. It provides access to the\n"
1047 "mathematical functions defined by the C standard.");
1054 m
= Py_InitModule3("math", math_methods
, module_doc
);
1058 PyModule_AddObject(m
, "pi", PyFloat_FromDouble(Py_MATH_PI
));
1059 PyModule_AddObject(m
, "e", PyFloat_FromDouble(Py_MATH_E
));