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);
64 sin(pi*x), giving accurate results for all finite x (especially x
65 integral or close to an integer). This is here for use in the
66 reflection formula for the gamma function. It conforms to IEEE
67 754-2008 for finite arguments, but not for infinities or nans.
70 static const double pi
= 3.141592653589793238462643383279502884197;
77 /* this function should only ever be called for finite arguments */
78 assert(Py_IS_FINITE(x
));
79 y
= fmod(fabs(x
), 2.0);
80 n
= (int)round(2.0*y
);
81 assert(0 <= n
&& n
<= 4);
90 /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
91 -0.0 instead of 0.0 when y == 1.0. */
101 assert(0); /* should never get here */
102 r
= -1.23e200
; /* silence gcc warning */
104 return copysign(1.0, x
)*r
;
107 /* Implementation of the real gamma function. In extensive but non-exhaustive
108 random tests, this function proved accurate to within <= 10 ulps across the
109 entire float domain. Note that accuracy may depend on the quality of the
110 system math functions, the pow function in particular. Special cases
111 follow C99 annex F. The parameters and method are tailored to platforms
112 whose double format is the IEEE 754 binary64 format.
114 Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
115 and g=6.024680040776729583740234375; these parameters are amongst those
116 used by the Boost library. Following Boost (again), we re-express the
117 Lanczos sum as a rational function, and compute it that way. The
118 coefficients below were computed independently using MPFR, and have been
119 double-checked against the coefficients in the Boost source code.
121 For x < 0.0 we use the reflection formula.
123 There's one minor tweak that deserves explanation: Lanczos' formula for
124 Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
125 values, x+g-0.5 can be represented exactly. However, in cases where it
126 can't be represented exactly the small error in x+g-0.5 can be magnified
127 significantly by the pow and exp calls, especially for large x. A cheap
128 correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
129 involved in the computation of x+g-0.5 (that is, e = computed value of
130 x+g-0.5 - exact value of x+g-0.5). Here's the proof:
134 Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
135 double, and e is tiny. Then:
137 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
138 = pow(y, x-0.5)/exp(y) * C,
140 where the correction_factor C is given by
142 C = pow(1-e/y, x-0.5) * exp(e)
144 Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
146 C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
148 But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
150 pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
152 Note that for accuracy, when computing r*C it's better to do
160 since the addition in the latter throws away most of the bits of
161 information in e*g/y.
165 static const double lanczos_g
= 6.024680040776729583740234375;
166 static const double lanczos_g_minus_half
= 5.524680040776729583740234375;
167 static const double lanczos_num_coeffs
[LANCZOS_N
] = {
168 23531376880.410759688572007674451636754734846804940,
169 42919803642.649098768957899047001988850926355848959,
170 35711959237.355668049440185451547166705960488635843,
171 17921034426.037209699919755754458931112671403265390,
172 6039542586.3520280050642916443072979210699388420708,
173 1439720407.3117216736632230727949123939715485786772,
174 248874557.86205415651146038641322942321632125127801,
175 31426415.585400194380614231628318205362874684987640,
176 2876370.6289353724412254090516208496135991145378768,
177 186056.26539522349504029498971604569928220784236328,
178 8071.6720023658162106380029022722506138218516325024,
179 210.82427775157934587250973392071336271166969580291,
180 2.5066282746310002701649081771338373386264310793408
183 /* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
184 static const double lanczos_den_coeffs
[LANCZOS_N
] = {
185 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
186 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
188 /* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
189 #define NGAMMA_INTEGRAL 23
190 static const double gamma_integral
[NGAMMA_INTEGRAL
] = {
191 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
192 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
193 1307674368000.0, 20922789888000.0, 355687428096000.0,
194 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
195 51090942171709440000.0, 1124000727777607680000.0,
198 /* Lanczos' sum L_g(x), for positive x */
201 lanczos_sum(double x
)
203 double num
= 0.0, den
= 0.0;
206 /* evaluate the rational function lanczos_sum(x). For large
207 x, the obvious algorithm risks overflow, so we instead
208 rescale the denominator and numerator of the rational
209 function by x**(1-LANCZOS_N) and treat this as a
210 rational function in 1/x. This also reduces the error for
211 larger x values. The choice of cutoff point (5.0 below) is
212 somewhat arbitrary; in tests, smaller cutoff values than
213 this resulted in lower accuracy. */
215 for (i
= LANCZOS_N
; --i
>= 0; ) {
216 num
= num
* x
+ lanczos_num_coeffs
[i
];
217 den
= den
* x
+ lanczos_den_coeffs
[i
];
221 for (i
= 0; i
< LANCZOS_N
; i
++) {
222 num
= num
/ x
+ lanczos_num_coeffs
[i
];
223 den
= den
/ x
+ lanczos_den_coeffs
[i
];
232 double absx
, r
, y
, z
, sqrtpow
;
235 if (!Py_IS_FINITE(x
)) {
236 if (Py_IS_NAN(x
) || x
> 0.0)
237 return x
; /* tgamma(nan) = nan, tgamma(inf) = inf */
240 return Py_NAN
; /* tgamma(-inf) = nan, invalid */
245 return 1.0/x
; /* tgamma(+-0.0) = +-inf, divide-by-zero */
248 /* integer arguments */
251 errno
= EDOM
; /* tgamma(n) = nan, invalid for */
252 return Py_NAN
; /* negative integers n */
254 if (x
<= NGAMMA_INTEGRAL
)
255 return gamma_integral
[(int)x
- 1];
259 /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
262 if (Py_IS_INFINITY(r
))
267 /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
268 x > 200, and underflows to +-0.0 for x < -200, not a negative
280 y
= absx
+ lanczos_g_minus_half
;
281 /* compute error in sum */
282 if (absx
> lanczos_g_minus_half
) {
283 /* note: the correction can be foiled by an optimizing
284 compiler that (incorrectly) thinks that an expression like
285 a + b - a - b can be optimized to 0.0. This shouldn't
286 happen in a standards-conforming compiler. */
288 z
= q
- lanczos_g_minus_half
;
291 double q
= y
- lanczos_g_minus_half
;
294 z
= z
* lanczos_g
/ y
;
296 r
= -pi
/ sinpi(absx
) / absx
* exp(y
) / lanczos_sum(absx
);
299 r
/= pow(y
, absx
- 0.5);
302 sqrtpow
= pow(y
, absx
/ 2.0 - 0.25);
308 r
= lanczos_sum(absx
) / exp(y
);
311 r
*= pow(y
, absx
- 0.5);
314 sqrtpow
= pow(y
, absx
/ 2.0 - 0.25);
319 if (Py_IS_INFINITY(r
))
325 lgamma: natural log of the absolute value of the Gamma function.
326 For large arguments, Lanczos' formula works extremely well here.
335 if (!Py_IS_FINITE(x
)) {
337 return x
; /* lgamma(nan) = nan */
339 return Py_HUGE_VAL
; /* lgamma(+-inf) = +inf */
342 /* integer arguments */
343 if (x
== floor(x
) && x
<= 2.0) {
345 errno
= EDOM
; /* lgamma(n) = inf, divide-by-zero for */
346 return Py_HUGE_VAL
; /* integers n <= 0 */
349 return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
354 /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
358 /* Lanczos' formula */
360 /* we could save a fraction of a ulp in accuracy by having a
361 second set of numerator coefficients for lanczos_sum that
362 absorbed the exp(-lanczos_g) term, and throwing out the
363 lanczos_g subtraction below; it's probably not worth it. */
364 r
= log(lanczos_sum(x
)) - lanczos_g
+
365 (x
-0.5)*(log(x
+lanczos_g
-0.5)-1);
368 r
= log(pi
) - log(fabs(sinpi(absx
))) - log(absx
) -
369 (log(lanczos_sum(absx
)) - lanczos_g
+
370 (absx
-0.5)*(log(absx
+lanczos_g
-0.5)-1));
372 if (Py_IS_INFINITY(r
))
379 wrapper for atan2 that deals directly with special cases before
380 delegating to the platform libm for the remaining cases. This
381 is necessary to get consistent behaviour across platforms.
382 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
387 m_atan2(double y
, double x
)
389 if (Py_IS_NAN(x
) || Py_IS_NAN(y
))
391 if (Py_IS_INFINITY(y
)) {
392 if (Py_IS_INFINITY(x
)) {
393 if (copysign(1., x
) == 1.)
394 /* atan2(+-inf, +inf) == +-pi/4 */
395 return copysign(0.25*Py_MATH_PI
, y
);
397 /* atan2(+-inf, -inf) == +-pi*3/4 */
398 return copysign(0.75*Py_MATH_PI
, y
);
400 /* atan2(+-inf, x) == +-pi/2 for finite x */
401 return copysign(0.5*Py_MATH_PI
, y
);
403 if (Py_IS_INFINITY(x
) || y
== 0.) {
404 if (copysign(1., x
) == 1.)
405 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
406 return copysign(0., y
);
408 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
409 return copysign(Py_MATH_PI
, y
);
415 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
416 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
417 special values directly, passing positive non-special values through to
418 the system log/log10.
424 if (Py_IS_FINITE(x
)) {
429 return -Py_HUGE_VAL
; /* log(0) = -inf */
431 return Py_NAN
; /* log(-ve) = nan */
433 else if (Py_IS_NAN(x
))
434 return x
; /* log(nan) = nan */
436 return x
; /* log(inf) = inf */
439 return Py_NAN
; /* log(-inf) = nan */
446 if (Py_IS_FINITE(x
)) {
451 return -Py_HUGE_VAL
; /* log10(0) = -inf */
453 return Py_NAN
; /* log10(-ve) = nan */
455 else if (Py_IS_NAN(x
))
456 return x
; /* log10(nan) = nan */
458 return x
; /* log10(inf) = inf */
461 return Py_NAN
; /* log10(-inf) = nan */
466 /* Call is_error when errno != 0, and where x is the result libm
467 * returned. is_error will usually set up an exception and return
468 * true (1), but may return false (0) without setting up an exception.
473 int result
= 1; /* presumption of guilt */
474 assert(errno
); /* non-zero errno is a precondition for calling */
476 PyErr_SetString(PyExc_ValueError
, "math domain error");
478 else if (errno
== ERANGE
) {
479 /* ANSI C generally requires libm functions to set ERANGE
480 * on overflow, but also generally *allows* them to set
481 * ERANGE on underflow too. There's no consistency about
482 * the latter across platforms.
483 * Alas, C99 never requires that errno be set.
484 * Here we suppress the underflow errors (libm functions
485 * should return a zero on underflow, and +- HUGE_VAL on
486 * overflow, so testing the result for zero suffices to
487 * distinguish the cases).
489 * On some platforms (Ubuntu/ia64) it seems that errno can be
490 * set to ERANGE for subnormal results that do *not* underflow
491 * to zero. So to be safe, we'll ignore ERANGE whenever the
492 * function result is less than one in absolute value.
497 PyErr_SetString(PyExc_OverflowError
,
501 /* Unexpected math error */
502 PyErr_SetFromErrno(PyExc_ValueError
);
507 math_1 is used to wrap a libm function f that takes a double
508 arguments and returns a double.
510 The error reporting follows these rules, which are designed to do
511 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
514 - a NaN result from non-NaN inputs causes ValueError to be raised
515 - an infinite result from finite inputs causes OverflowError to be
516 raised if can_overflow is 1, or raises ValueError if can_overflow
518 - if the result is finite and errno == EDOM then ValueError is
520 - if the result is finite and nonzero and errno == ERANGE then
521 OverflowError is raised
523 The last rule is used to catch overflow on platforms which follow
524 C89 but for which HUGE_VAL is not an infinity.
526 For the majority of one-argument functions these rules are enough
527 to ensure that Python's functions behave as specified in 'Annex F'
528 of the C99 standard, with the 'invalid' and 'divide-by-zero'
529 floating-point exceptions mapping to Python's ValueError and the
530 'overflow' floating-point exception mapping to OverflowError.
531 math_1 only works for functions that don't have singularities *and*
532 the possibility of overflow; fortunately, that covers everything we
533 care about right now.
537 math_1(PyObject
*arg
, double (*func
) (double), int can_overflow
)
540 x
= PyFloat_AsDouble(arg
);
541 if (x
== -1.0 && PyErr_Occurred())
544 PyFPE_START_PROTECT("in math_1", return 0);
546 PyFPE_END_PROTECT(r
);
553 else if (Py_IS_INFINITY(r
)) {
555 errno
= can_overflow
? ERANGE
: EDOM
;
559 if (errno
&& is_error(r
))
562 return PyFloat_FromDouble(r
);
565 /* variant of math_1, to be used when the function being wrapped is known to
566 set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
567 errno = ERANGE for overflow). */
570 math_1a(PyObject
*arg
, double (*func
) (double))
573 x
= PyFloat_AsDouble(arg
);
574 if (x
== -1.0 && PyErr_Occurred())
577 PyFPE_START_PROTECT("in math_1a", return 0);
579 PyFPE_END_PROTECT(r
);
580 if (errno
&& is_error(r
))
582 return PyFloat_FromDouble(r
);
586 math_2 is used to wrap a libm function f that takes two double
587 arguments and returns a double.
589 The error reporting follows these rules, which are designed to do
590 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
593 - a NaN result from non-NaN inputs causes ValueError to be raised
594 - an infinite result from finite inputs causes OverflowError to be
596 - if the result is finite and errno == EDOM then ValueError is
598 - if the result is finite and nonzero and errno == ERANGE then
599 OverflowError is raised
601 The last rule is used to catch overflow on platforms which follow
602 C89 but for which HUGE_VAL is not an infinity.
604 For most two-argument functions (copysign, fmod, hypot, atan2)
605 these rules are enough to ensure that Python's functions behave as
606 specified in 'Annex F' of the C99 standard, with the 'invalid' and
607 'divide-by-zero' floating-point exceptions mapping to Python's
608 ValueError and the 'overflow' floating-point exception mapping to
613 math_2(PyObject
*args
, double (*func
) (double, double), char *funcname
)
617 if (! PyArg_UnpackTuple(args
, funcname
, 2, 2, &ox
, &oy
))
619 x
= PyFloat_AsDouble(ox
);
620 y
= PyFloat_AsDouble(oy
);
621 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
624 PyFPE_START_PROTECT("in math_2", return 0);
626 PyFPE_END_PROTECT(r
);
628 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
633 else if (Py_IS_INFINITY(r
)) {
634 if (Py_IS_FINITE(x
) && Py_IS_FINITE(y
))
639 if (errno
&& is_error(r
))
642 return PyFloat_FromDouble(r
);
645 #define FUNC1(funcname, func, can_overflow, docstring) \
646 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
647 return math_1(args, func, can_overflow); \
649 PyDoc_STRVAR(math_##funcname##_doc, docstring);
651 #define FUNC1A(funcname, func, docstring) \
652 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
653 return math_1a(args, func); \
655 PyDoc_STRVAR(math_##funcname##_doc, docstring);
657 #define FUNC2(funcname, func, docstring) \
658 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
659 return math_2(args, func, #funcname); \
661 PyDoc_STRVAR(math_##funcname##_doc, docstring);
664 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
665 FUNC1(acosh
, acosh
, 0,
666 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
668 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
669 FUNC1(asinh
, asinh
, 0,
670 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
672 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
673 FUNC2(atan2
, m_atan2
,
674 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
675 "Unlike atan(y/x), the signs of both x and y are considered.")
676 FUNC1(atanh
, atanh
, 0,
677 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
679 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
680 "This is the smallest integral value >= x.")
681 FUNC2(copysign
, copysign
,
682 "copysign(x, y)\n\nReturn x with the sign of y.")
684 "cos(x)\n\nReturn the cosine of x (measured in radians).")
686 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
688 "exp(x)\n\nReturn e raised to the power of x.")
690 "fabs(x)\n\nReturn the absolute value of the float x.")
691 FUNC1(floor
, floor
, 0,
692 "floor(x)\n\nReturn the floor of x as a float.\n"
693 "This is the largest integral value <= x.")
694 FUNC1A(gamma
, m_tgamma
,
695 "gamma(x)\n\nGamma function at x.")
696 FUNC1A(lgamma
, m_lgamma
,
697 "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
698 FUNC1(log1p
, log1p
, 1,
699 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
700 "The result is computed in a way which is accurate for x near zero.")
702 "sin(x)\n\nReturn the sine of x (measured in radians).")
704 "sinh(x)\n\nReturn the hyperbolic sine of x.")
706 "sqrt(x)\n\nReturn the square root of x.")
708 "tan(x)\n\nReturn the tangent of x (measured in radians).")
710 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
712 /* Precision summation function as msum() by Raymond Hettinger in
713 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
714 enhanced with the exact partials sum and roundoff from Mark
715 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
716 See those links for more details, proofs and other references.
718 Note 1: IEEE 754R floating point semantics are assumed,
719 but the current implementation does not re-establish special
720 value semantics across iterations (i.e. handling -Inf + Inf).
722 Note 2: No provision is made for intermediate overflow handling;
723 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
724 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
725 overflow of the first partial sum.
727 Note 3: The intermediate values lo, yr, and hi are declared volatile so
728 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
729 Also, the volatile declaration forces the values to be stored in memory as
730 regular doubles instead of extended long precision (80-bit) values. This
731 prevents double rounding because any addition or subtraction of two doubles
732 can be resolved exactly into double-sized hi and lo values. As long as the
733 hi value gets forced into a double before yr and lo are computed, the extra
734 bits in downstream extended precision operations (x87 for example) will be
735 exactly zero and therefore can be losslessly stored back into a double,
736 thereby preventing double rounding.
738 Note 4: A similar implementation is in Modules/cmathmodule.c.
739 Be sure to update both when making changes.
741 Note 5: The signature of math.fsum() differs from __builtin__.sum()
742 because the start argument doesn't make sense in the context of
743 accurate summation. Since the partials table is collapsed before
744 returning a result, sum(seq2, start=sum(seq1)) may not equal the
745 accurate result returned by sum(itertools.chain(seq1, seq2)).
748 #define NUM_PARTIALS 32 /* initial partials array size, on stack */
750 /* Extend the partials array p[] by doubling its size. */
751 static int /* non-zero on error */
752 _fsum_realloc(double **p_ptr
, Py_ssize_t n
,
753 double *ps
, Py_ssize_t
*m_ptr
)
756 Py_ssize_t m
= *m_ptr
;
759 if (n
< m
&& m
< (PY_SSIZE_T_MAX
/ sizeof(double))) {
762 v
= PyMem_Malloc(sizeof(double) * m
);
764 memcpy(v
, ps
, sizeof(double) * n
);
767 v
= PyMem_Realloc(p
, sizeof(double) * m
);
769 if (v
== NULL
) { /* size overflow or no memory */
770 PyErr_SetString(PyExc_MemoryError
, "math.fsum partials");
773 *p_ptr
= (double*) v
;
778 /* Full precision summation of a sequence of floats.
781 partials = [] # sorted, non-overlapping partial sums
794 return sum_exact(partials)
796 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
797 are exactly equal to x+y. The inner loop applies hi/lo summation to each
798 partial so that the list of partial sums remains exact.
800 Sum_exact() adds the partial sums exactly and correctly rounds the final
801 result (using the round-half-to-even rule). The items in partials remain
802 non-zero, non-special, non-overlapping and strictly increasing in
803 magnitude, but possibly not all having the same sign.
805 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
809 math_fsum(PyObject
*self
, PyObject
*seq
)
811 PyObject
*item
, *iter
, *sum
= NULL
;
812 Py_ssize_t i
, j
, n
= 0, m
= NUM_PARTIALS
;
813 double x
, y
, t
, ps
[NUM_PARTIALS
], *p
= ps
;
814 double xsave
, special_sum
= 0.0, inf_sum
= 0.0;
815 volatile double hi
, yr
, lo
;
817 iter
= PyObject_GetIter(seq
);
821 PyFPE_START_PROTECT("fsum", Py_DECREF(iter
); return NULL
)
823 for(;;) { /* for x in iterable */
824 assert(0 <= n
&& n
<= m
);
825 assert((m
== NUM_PARTIALS
&& p
== ps
) ||
826 (m
> NUM_PARTIALS
&& p
!= NULL
));
828 item
= PyIter_Next(iter
);
830 if (PyErr_Occurred())
834 x
= PyFloat_AsDouble(item
);
836 if (PyErr_Occurred())
840 for (i
= j
= 0; j
< n
; j
++) { /* for y in partials */
842 if (fabs(x
) < fabs(y
)) {
853 n
= i
; /* ps[i:] = [x] */
855 if (! Py_IS_FINITE(x
)) {
856 /* a nonfinite x could arise either as
857 a result of intermediate overflow, or
858 as a result of a nan or inf in the
860 if (Py_IS_FINITE(xsave
)) {
861 PyErr_SetString(PyExc_OverflowError
,
862 "intermediate overflow in fsum");
865 if (Py_IS_INFINITY(xsave
))
867 special_sum
+= xsave
;
871 else if (n
>= m
&& _fsum_realloc(&p
, n
, ps
, &m
))
878 if (special_sum
!= 0.0) {
879 if (Py_IS_NAN(inf_sum
))
880 PyErr_SetString(PyExc_ValueError
,
881 "-inf + inf in fsum");
883 sum
= PyFloat_FromDouble(special_sum
);
890 /* sum_exact(ps, hi) from the top, stop when the sum becomes
895 assert(fabs(y
) < fabs(x
));
902 /* Make half-even rounding work across multiple partials.
903 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
904 digit to two instead of down to zero (the 1e-16 makes the 1
905 slightly closer to two). With a potential 1 ULP rounding
906 error fixed-up, math.fsum() can guarantee commutativity. */
907 if (n
> 0 && ((lo
< 0.0 && p
[n
-1] < 0.0) ||
908 (lo
> 0.0 && p
[n
-1] > 0.0))) {
916 sum
= PyFloat_FromDouble(hi
);
919 PyFPE_END_PROTECT(hi
)
928 PyDoc_STRVAR(math_fsum_doc
,
930 Return an accurate floating point sum of values in the iterable.\n\
931 Assumes IEEE-754 floating point arithmetic.");
934 math_factorial(PyObject
*self
, PyObject
*arg
)
937 PyObject
*result
, *iobj
, *newresult
;
939 if (PyFloat_Check(arg
)) {
940 double dx
= PyFloat_AS_DOUBLE((PyFloatObject
*)arg
);
941 if (dx
!= floor(dx
)) {
942 PyErr_SetString(PyExc_ValueError
,
943 "factorial() only accepts integral values");
948 x
= PyInt_AsLong(arg
);
949 if (x
== -1 && PyErr_Occurred())
952 PyErr_SetString(PyExc_ValueError
,
953 "factorial() not defined for negative values");
957 result
= (PyObject
*)PyInt_FromLong(1);
960 for (i
=1 ; i
<=x
; i
++) {
961 iobj
= (PyObject
*)PyInt_FromLong(i
);
964 newresult
= PyNumber_Multiply(result
, iobj
);
966 if (newresult
== NULL
)
978 PyDoc_STRVAR(math_factorial_doc
,
979 "factorial(x) -> Integral\n"
981 "Find x!. Raise a ValueError if x is negative or non-integral.");
984 math_trunc(PyObject
*self
, PyObject
*number
)
986 return PyObject_CallMethod(number
, "__trunc__", NULL
);
989 PyDoc_STRVAR(math_trunc_doc
,
990 "trunc(x:Real) -> Integral\n"
992 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
995 math_frexp(PyObject
*self
, PyObject
*arg
)
998 double x
= PyFloat_AsDouble(arg
);
999 if (x
== -1.0 && PyErr_Occurred())
1001 /* deal with special cases directly, to sidestep platform
1003 if (Py_IS_NAN(x
) || Py_IS_INFINITY(x
) || !x
) {
1007 PyFPE_START_PROTECT("in math_frexp", return 0);
1009 PyFPE_END_PROTECT(x
);
1011 return Py_BuildValue("(di)", x
, i
);
1014 PyDoc_STRVAR(math_frexp_doc
,
1017 "Return the mantissa and exponent of x, as pair (m, e).\n"
1018 "m is a float and e is an int, such that x = m * 2.**e.\n"
1019 "If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
1022 math_ldexp(PyObject
*self
, PyObject
*args
)
1027 if (! PyArg_ParseTuple(args
, "dO:ldexp", &x
, &oexp
))
1030 if (PyLong_Check(oexp
)) {
1031 /* on overflow, replace exponent with either LONG_MAX
1032 or LONG_MIN, depending on the sign. */
1033 exp
= PyLong_AsLong(oexp
);
1034 if (exp
== -1 && PyErr_Occurred()) {
1035 if (PyErr_ExceptionMatches(PyExc_OverflowError
)) {
1036 if (Py_SIZE(oexp
) < 0) {
1045 /* propagate any unexpected exception */
1050 else if (PyInt_Check(oexp
)) {
1051 exp
= PyInt_AS_LONG(oexp
);
1054 PyErr_SetString(PyExc_TypeError
,
1055 "Expected an int or long as second argument "
1060 if (x
== 0. || !Py_IS_FINITE(x
)) {
1061 /* NaNs, zeros and infinities are returned unchanged */
1064 } else if (exp
> INT_MAX
) {
1066 r
= copysign(Py_HUGE_VAL
, x
);
1068 } else if (exp
< INT_MIN
) {
1069 /* underflow to +-0 */
1070 r
= copysign(0., x
);
1074 PyFPE_START_PROTECT("in math_ldexp", return 0);
1075 r
= ldexp(x
, (int)exp
);
1076 PyFPE_END_PROTECT(r
);
1077 if (Py_IS_INFINITY(r
))
1081 if (errno
&& is_error(r
))
1083 return PyFloat_FromDouble(r
);
1086 PyDoc_STRVAR(math_ldexp_doc
,
1088 Return x * (2**i).");
1091 math_modf(PyObject
*self
, PyObject
*arg
)
1093 double y
, x
= PyFloat_AsDouble(arg
);
1094 if (x
== -1.0 && PyErr_Occurred())
1096 /* some platforms don't do the right thing for NaNs and
1097 infinities, so we take care of special cases directly. */
1098 if (!Py_IS_FINITE(x
)) {
1099 if (Py_IS_INFINITY(x
))
1100 return Py_BuildValue("(dd)", copysign(0., x
), x
);
1101 else if (Py_IS_NAN(x
))
1102 return Py_BuildValue("(dd)", x
, x
);
1106 PyFPE_START_PROTECT("in math_modf", return 0);
1108 PyFPE_END_PROTECT(x
);
1109 return Py_BuildValue("(dd)", x
, y
);
1112 PyDoc_STRVAR(math_modf_doc
,
1115 "Return the fractional and integer parts of x. Both results carry the sign\n"
1116 "of x and are floats.");
1118 /* A decent logarithm is easy to compute even for huge longs, but libm can't
1119 do that by itself -- loghelper can. func is log or log10, and name is
1120 "log" or "log10". Note that overflow isn't possible: a long can contain
1121 no more than INT_MAX * SHIFT bits, so has value certainly less than
1122 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
1123 small enough to fit in an IEEE single. log and log10 are even smaller.
1127 loghelper(PyObject
* arg
, double (*func
)(double), char *funcname
)
1129 /* If it is long, do it ourselves. */
1130 if (PyLong_Check(arg
)) {
1133 x
= _PyLong_AsScaledDouble(arg
, &e
);
1135 PyErr_SetString(PyExc_ValueError
,
1136 "math domain error");
1139 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
1140 log(x) + log(2) * e * PyLong_SHIFT.
1141 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
1142 so force use of double. */
1143 x
= func(x
) + (e
* (double)PyLong_SHIFT
) * func(2.0);
1144 return PyFloat_FromDouble(x
);
1147 /* Else let libm handle it by itself. */
1148 return math_1(arg
, func
, 0);
1152 math_log(PyObject
*self
, PyObject
*args
)
1155 PyObject
*base
= NULL
;
1156 PyObject
*num
, *den
;
1159 if (!PyArg_UnpackTuple(args
, "log", 1, 2, &arg
, &base
))
1162 num
= loghelper(arg
, m_log
, "log");
1163 if (num
== NULL
|| base
== NULL
)
1166 den
= loghelper(base
, m_log
, "log");
1172 ans
= PyNumber_Divide(num
, den
);
1178 PyDoc_STRVAR(math_log_doc
,
1179 "log(x[, base])\n\n\
1180 Return the logarithm of x to the given base.\n\
1181 If the base not specified, returns the natural logarithm (base e) of x.");
1184 math_log10(PyObject
*self
, PyObject
*arg
)
1186 return loghelper(arg
, m_log10
, "log10");
1189 PyDoc_STRVAR(math_log10_doc
,
1190 "log10(x)\n\nReturn the base 10 logarithm of x.");
1193 math_fmod(PyObject
*self
, PyObject
*args
)
1197 if (! PyArg_UnpackTuple(args
, "fmod", 2, 2, &ox
, &oy
))
1199 x
= PyFloat_AsDouble(ox
);
1200 y
= PyFloat_AsDouble(oy
);
1201 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
1203 /* fmod(x, +/-Inf) returns x for finite x. */
1204 if (Py_IS_INFINITY(y
) && Py_IS_FINITE(x
))
1205 return PyFloat_FromDouble(x
);
1207 PyFPE_START_PROTECT("in math_fmod", return 0);
1209 PyFPE_END_PROTECT(r
);
1211 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
1216 if (errno
&& is_error(r
))
1219 return PyFloat_FromDouble(r
);
1222 PyDoc_STRVAR(math_fmod_doc
,
1223 "fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
1224 " x % y may differ.");
1227 math_hypot(PyObject
*self
, PyObject
*args
)
1231 if (! PyArg_UnpackTuple(args
, "hypot", 2, 2, &ox
, &oy
))
1233 x
= PyFloat_AsDouble(ox
);
1234 y
= PyFloat_AsDouble(oy
);
1235 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
1237 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1238 if (Py_IS_INFINITY(x
))
1239 return PyFloat_FromDouble(fabs(x
));
1240 if (Py_IS_INFINITY(y
))
1241 return PyFloat_FromDouble(fabs(y
));
1243 PyFPE_START_PROTECT("in math_hypot", return 0);
1245 PyFPE_END_PROTECT(r
);
1247 if (!Py_IS_NAN(x
) && !Py_IS_NAN(y
))
1252 else if (Py_IS_INFINITY(r
)) {
1253 if (Py_IS_FINITE(x
) && Py_IS_FINITE(y
))
1258 if (errno
&& is_error(r
))
1261 return PyFloat_FromDouble(r
);
1264 PyDoc_STRVAR(math_hypot_doc
,
1265 "hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
1267 /* pow can't use math_2, but needs its own wrapper: the problem is
1268 that an infinite result can arise either as a result of overflow
1269 (in which case OverflowError should be raised) or as a result of
1270 e.g. 0.**-5. (for which ValueError needs to be raised.)
1274 math_pow(PyObject
*self
, PyObject
*args
)
1280 if (! PyArg_UnpackTuple(args
, "pow", 2, 2, &ox
, &oy
))
1282 x
= PyFloat_AsDouble(ox
);
1283 y
= PyFloat_AsDouble(oy
);
1284 if ((x
== -1.0 || y
== -1.0) && PyErr_Occurred())
1287 /* deal directly with IEEE specials, to cope with problems on various
1288 platforms whose semantics don't exactly match C99 */
1289 r
= 0.; /* silence compiler warning */
1290 if (!Py_IS_FINITE(x
) || !Py_IS_FINITE(y
)) {
1293 r
= y
== 0. ? 1. : x
; /* NaN**0 = 1 */
1294 else if (Py_IS_NAN(y
))
1295 r
= x
== 1. ? 1. : y
; /* 1**NaN = 1 */
1296 else if (Py_IS_INFINITY(x
)) {
1297 odd_y
= Py_IS_FINITE(y
) && fmod(fabs(y
), 2.0) == 1.0;
1299 r
= odd_y
? x
: fabs(x
);
1303 r
= odd_y
? copysign(0., x
) : 0.;
1305 else if (Py_IS_INFINITY(y
)) {
1308 else if (y
> 0. && fabs(x
) > 1.0)
1310 else if (y
< 0. && fabs(x
) < 1.0) {
1311 r
= -y
; /* result is +inf */
1312 if (x
== 0.) /* 0**-inf: divide-by-zero */
1320 /* let libm handle finite**finite */
1322 PyFPE_START_PROTECT("in math_pow", return 0);
1324 PyFPE_END_PROTECT(r
);
1325 /* a NaN result should arise only from (-ve)**(finite
1326 non-integer); in this case we want to raise ValueError. */
1327 if (!Py_IS_FINITE(r
)) {
1332 an infinite result here arises either from:
1333 (A) (+/-0.)**negative (-> divide-by-zero)
1334 (B) overflow of x**y with x and y finite
1336 else if (Py_IS_INFINITY(r
)) {
1345 if (errno
&& is_error(r
))
1348 return PyFloat_FromDouble(r
);
1351 PyDoc_STRVAR(math_pow_doc
,
1352 "pow(x, y)\n\nReturn x**y (x to the power of y).");
1354 static const double degToRad
= Py_MATH_PI
/ 180.0;
1355 static const double radToDeg
= 180.0 / Py_MATH_PI
;
1358 math_degrees(PyObject
*self
, PyObject
*arg
)
1360 double x
= PyFloat_AsDouble(arg
);
1361 if (x
== -1.0 && PyErr_Occurred())
1363 return PyFloat_FromDouble(x
* radToDeg
);
1366 PyDoc_STRVAR(math_degrees_doc
,
1368 Convert angle x from radians to degrees.");
1371 math_radians(PyObject
*self
, PyObject
*arg
)
1373 double x
= PyFloat_AsDouble(arg
);
1374 if (x
== -1.0 && PyErr_Occurred())
1376 return PyFloat_FromDouble(x
* degToRad
);
1379 PyDoc_STRVAR(math_radians_doc
,
1381 Convert angle x from degrees to radians.");
1384 math_isnan(PyObject
*self
, PyObject
*arg
)
1386 double x
= PyFloat_AsDouble(arg
);
1387 if (x
== -1.0 && PyErr_Occurred())
1389 return PyBool_FromLong((long)Py_IS_NAN(x
));
1392 PyDoc_STRVAR(math_isnan_doc
,
1393 "isnan(x) -> bool\n\n\
1394 Check if float x is not a number (NaN).");
1397 math_isinf(PyObject
*self
, PyObject
*arg
)
1399 double x
= PyFloat_AsDouble(arg
);
1400 if (x
== -1.0 && PyErr_Occurred())
1402 return PyBool_FromLong((long)Py_IS_INFINITY(x
));
1405 PyDoc_STRVAR(math_isinf_doc
,
1406 "isinf(x) -> bool\n\n\
1407 Check if float x is infinite (positive or negative).");
1409 static PyMethodDef math_methods
[] = {
1410 {"acos", math_acos
, METH_O
, math_acos_doc
},
1411 {"acosh", math_acosh
, METH_O
, math_acosh_doc
},
1412 {"asin", math_asin
, METH_O
, math_asin_doc
},
1413 {"asinh", math_asinh
, METH_O
, math_asinh_doc
},
1414 {"atan", math_atan
, METH_O
, math_atan_doc
},
1415 {"atan2", math_atan2
, METH_VARARGS
, math_atan2_doc
},
1416 {"atanh", math_atanh
, METH_O
, math_atanh_doc
},
1417 {"ceil", math_ceil
, METH_O
, math_ceil_doc
},
1418 {"copysign", math_copysign
, METH_VARARGS
, math_copysign_doc
},
1419 {"cos", math_cos
, METH_O
, math_cos_doc
},
1420 {"cosh", math_cosh
, METH_O
, math_cosh_doc
},
1421 {"degrees", math_degrees
, METH_O
, math_degrees_doc
},
1422 {"exp", math_exp
, METH_O
, math_exp_doc
},
1423 {"fabs", math_fabs
, METH_O
, math_fabs_doc
},
1424 {"factorial", math_factorial
, METH_O
, math_factorial_doc
},
1425 {"floor", math_floor
, METH_O
, math_floor_doc
},
1426 {"fmod", math_fmod
, METH_VARARGS
, math_fmod_doc
},
1427 {"frexp", math_frexp
, METH_O
, math_frexp_doc
},
1428 {"fsum", math_fsum
, METH_O
, math_fsum_doc
},
1429 {"gamma", math_gamma
, METH_O
, math_gamma_doc
},
1430 {"hypot", math_hypot
, METH_VARARGS
, math_hypot_doc
},
1431 {"isinf", math_isinf
, METH_O
, math_isinf_doc
},
1432 {"isnan", math_isnan
, METH_O
, math_isnan_doc
},
1433 {"ldexp", math_ldexp
, METH_VARARGS
, math_ldexp_doc
},
1434 {"lgamma", math_lgamma
, METH_O
, math_lgamma_doc
},
1435 {"log", math_log
, METH_VARARGS
, math_log_doc
},
1436 {"log1p", math_log1p
, METH_O
, math_log1p_doc
},
1437 {"log10", math_log10
, METH_O
, math_log10_doc
},
1438 {"modf", math_modf
, METH_O
, math_modf_doc
},
1439 {"pow", math_pow
, METH_VARARGS
, math_pow_doc
},
1440 {"radians", math_radians
, METH_O
, math_radians_doc
},
1441 {"sin", math_sin
, METH_O
, math_sin_doc
},
1442 {"sinh", math_sinh
, METH_O
, math_sinh_doc
},
1443 {"sqrt", math_sqrt
, METH_O
, math_sqrt_doc
},
1444 {"tan", math_tan
, METH_O
, math_tan_doc
},
1445 {"tanh", math_tanh
, METH_O
, math_tanh_doc
},
1446 {"trunc", math_trunc
, METH_O
, math_trunc_doc
},
1447 {NULL
, NULL
} /* sentinel */
1451 PyDoc_STRVAR(module_doc
,
1452 "This module is always available. It provides access to the\n"
1453 "mathematical functions defined by the C standard.");
1460 m
= Py_InitModule3("math", math_methods
, module_doc
);
1464 PyModule_AddObject(m
, "pi", PyFloat_FromDouble(Py_MATH_PI
));
1465 PyModule_AddObject(m
, "e", PyFloat_FromDouble(Py_MATH_E
));