Fix duplicate test numbers in extra.decTest
[python.git] / Modules / mathmodule.c
blobb2fee3ff3f83df4aa8046481857889d46141e78d
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
7 exceptions.
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
37 accident ;-)).
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
52 returned.
55 #include "Python.h"
56 #include "longintrepr.h" /* just for SHIFT */
58 #ifdef _OSF_SOURCE
59 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
60 extern double copysign(double, double);
61 #endif
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;
72 static double
73 sinpi(double x)
75 double y, r;
76 int n;
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);
82 switch (n) {
83 case 0:
84 r = sin(pi*y);
85 break;
86 case 1:
87 r = cos(pi*(y-0.5));
88 break;
89 case 2:
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. */
92 r = sin(pi*(1.0-y));
93 break;
94 case 3:
95 r = -cos(pi*(y-1.5));
96 break;
97 case 4:
98 r = sin(pi*(y-2.0));
99 break;
100 default:
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:
132 Correction factor
133 -----------------
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
154 r + e*g/y*r;
156 than
158 r * (1 + e*g/y);
160 since the addition in the latter throws away most of the bits of
161 information in e*g/y.
164 #define LANCZOS_N 13
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 */
200 static double
201 lanczos_sum(double x)
203 double num = 0.0, den = 0.0;
204 int i;
205 assert(x > 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. */
214 if (x < 5.0) {
215 for (i = LANCZOS_N; --i >= 0; ) {
216 num = num * x + lanczos_num_coeffs[i];
217 den = den * x + lanczos_den_coeffs[i];
220 else {
221 for (i = 0; i < LANCZOS_N; i++) {
222 num = num / x + lanczos_num_coeffs[i];
223 den = den / x + lanczos_den_coeffs[i];
226 return num/den;
229 static double
230 m_tgamma(double x)
232 double absx, r, y, z, sqrtpow;
234 /* special cases */
235 if (!Py_IS_FINITE(x)) {
236 if (Py_IS_NAN(x) || x > 0.0)
237 return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
238 else {
239 errno = EDOM;
240 return Py_NAN; /* tgamma(-inf) = nan, invalid */
243 if (x == 0.0) {
244 errno = EDOM;
245 return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */
248 /* integer arguments */
249 if (x == floor(x)) {
250 if (x < 0.0) {
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];
257 absx = fabs(x);
259 /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
260 if (absx < 1e-20) {
261 r = 1.0/x;
262 if (Py_IS_INFINITY(r))
263 errno = ERANGE;
264 return 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
269 integer. */
270 if (absx > 200.0) {
271 if (x < 0.0) {
272 return 0.0/sinpi(x);
274 else {
275 errno = ERANGE;
276 return Py_HUGE_VAL;
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. */
287 double q = y - absx;
288 z = q - lanczos_g_minus_half;
290 else {
291 double q = y - lanczos_g_minus_half;
292 z = q - absx;
294 z = z * lanczos_g / y;
295 if (x < 0.0) {
296 r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
297 r -= z * r;
298 if (absx < 140.0) {
299 r /= pow(y, absx - 0.5);
301 else {
302 sqrtpow = pow(y, absx / 2.0 - 0.25);
303 r /= sqrtpow;
304 r /= sqrtpow;
307 else {
308 r = lanczos_sum(absx) / exp(y);
309 r += z * r;
310 if (absx < 140.0) {
311 r *= pow(y, absx - 0.5);
313 else {
314 sqrtpow = pow(y, absx / 2.0 - 0.25);
315 r *= sqrtpow;
316 r *= sqrtpow;
319 if (Py_IS_INFINITY(r))
320 errno = ERANGE;
321 return r;
325 wrapper for atan2 that deals directly with special cases before
326 delegating to the platform libm for the remaining cases. This
327 is necessary to get consistent behaviour across platforms.
328 Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
329 always follow C99.
332 static double
333 m_atan2(double y, double x)
335 if (Py_IS_NAN(x) || Py_IS_NAN(y))
336 return Py_NAN;
337 if (Py_IS_INFINITY(y)) {
338 if (Py_IS_INFINITY(x)) {
339 if (copysign(1., x) == 1.)
340 /* atan2(+-inf, +inf) == +-pi/4 */
341 return copysign(0.25*Py_MATH_PI, y);
342 else
343 /* atan2(+-inf, -inf) == +-pi*3/4 */
344 return copysign(0.75*Py_MATH_PI, y);
346 /* atan2(+-inf, x) == +-pi/2 for finite x */
347 return copysign(0.5*Py_MATH_PI, y);
349 if (Py_IS_INFINITY(x) || y == 0.) {
350 if (copysign(1., x) == 1.)
351 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
352 return copysign(0., y);
353 else
354 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
355 return copysign(Py_MATH_PI, y);
357 return atan2(y, x);
361 Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
362 log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
363 special values directly, passing positive non-special values through to
364 the system log/log10.
367 static double
368 m_log(double x)
370 if (Py_IS_FINITE(x)) {
371 if (x > 0.0)
372 return log(x);
373 errno = EDOM;
374 if (x == 0.0)
375 return -Py_HUGE_VAL; /* log(0) = -inf */
376 else
377 return Py_NAN; /* log(-ve) = nan */
379 else if (Py_IS_NAN(x))
380 return x; /* log(nan) = nan */
381 else if (x > 0.0)
382 return x; /* log(inf) = inf */
383 else {
384 errno = EDOM;
385 return Py_NAN; /* log(-inf) = nan */
389 static double
390 m_log10(double x)
392 if (Py_IS_FINITE(x)) {
393 if (x > 0.0)
394 return log10(x);
395 errno = EDOM;
396 if (x == 0.0)
397 return -Py_HUGE_VAL; /* log10(0) = -inf */
398 else
399 return Py_NAN; /* log10(-ve) = nan */
401 else if (Py_IS_NAN(x))
402 return x; /* log10(nan) = nan */
403 else if (x > 0.0)
404 return x; /* log10(inf) = inf */
405 else {
406 errno = EDOM;
407 return Py_NAN; /* log10(-inf) = nan */
412 /* Call is_error when errno != 0, and where x is the result libm
413 * returned. is_error will usually set up an exception and return
414 * true (1), but may return false (0) without setting up an exception.
416 static int
417 is_error(double x)
419 int result = 1; /* presumption of guilt */
420 assert(errno); /* non-zero errno is a precondition for calling */
421 if (errno == EDOM)
422 PyErr_SetString(PyExc_ValueError, "math domain error");
424 else if (errno == ERANGE) {
425 /* ANSI C generally requires libm functions to set ERANGE
426 * on overflow, but also generally *allows* them to set
427 * ERANGE on underflow too. There's no consistency about
428 * the latter across platforms.
429 * Alas, C99 never requires that errno be set.
430 * Here we suppress the underflow errors (libm functions
431 * should return a zero on underflow, and +- HUGE_VAL on
432 * overflow, so testing the result for zero suffices to
433 * distinguish the cases).
435 * On some platforms (Ubuntu/ia64) it seems that errno can be
436 * set to ERANGE for subnormal results that do *not* underflow
437 * to zero. So to be safe, we'll ignore ERANGE whenever the
438 * function result is less than one in absolute value.
440 if (fabs(x) < 1.0)
441 result = 0;
442 else
443 PyErr_SetString(PyExc_OverflowError,
444 "math range error");
446 else
447 /* Unexpected math error */
448 PyErr_SetFromErrno(PyExc_ValueError);
449 return result;
453 math_1 is used to wrap a libm function f that takes a double
454 arguments and returns a double.
456 The error reporting follows these rules, which are designed to do
457 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
458 platforms.
460 - a NaN result from non-NaN inputs causes ValueError to be raised
461 - an infinite result from finite inputs causes OverflowError to be
462 raised if can_overflow is 1, or raises ValueError if can_overflow
463 is 0.
464 - if the result is finite and errno == EDOM then ValueError is
465 raised
466 - if the result is finite and nonzero and errno == ERANGE then
467 OverflowError is raised
469 The last rule is used to catch overflow on platforms which follow
470 C89 but for which HUGE_VAL is not an infinity.
472 For the majority of one-argument functions these rules are enough
473 to ensure that Python's functions behave as specified in 'Annex F'
474 of the C99 standard, with the 'invalid' and 'divide-by-zero'
475 floating-point exceptions mapping to Python's ValueError and the
476 'overflow' floating-point exception mapping to OverflowError.
477 math_1 only works for functions that don't have singularities *and*
478 the possibility of overflow; fortunately, that covers everything we
479 care about right now.
482 static PyObject *
483 math_1(PyObject *arg, double (*func) (double), int can_overflow)
485 double x, r;
486 x = PyFloat_AsDouble(arg);
487 if (x == -1.0 && PyErr_Occurred())
488 return NULL;
489 errno = 0;
490 PyFPE_START_PROTECT("in math_1", return 0);
491 r = (*func)(x);
492 PyFPE_END_PROTECT(r);
493 if (Py_IS_NAN(r)) {
494 if (!Py_IS_NAN(x))
495 errno = EDOM;
496 else
497 errno = 0;
499 else if (Py_IS_INFINITY(r)) {
500 if (Py_IS_FINITE(x))
501 errno = can_overflow ? ERANGE : EDOM;
502 else
503 errno = 0;
505 if (errno && is_error(r))
506 return NULL;
507 else
508 return PyFloat_FromDouble(r);
511 /* variant of math_1, to be used when the function being wrapped is known to
512 set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
513 errno = ERANGE for overflow). */
515 static PyObject *
516 math_1a(PyObject *arg, double (*func) (double))
518 double x, r;
519 x = PyFloat_AsDouble(arg);
520 if (x == -1.0 && PyErr_Occurred())
521 return NULL;
522 errno = 0;
523 PyFPE_START_PROTECT("in math_1a", return 0);
524 r = (*func)(x);
525 PyFPE_END_PROTECT(r);
526 if (errno && is_error(r))
527 return NULL;
528 return PyFloat_FromDouble(r);
532 math_2 is used to wrap a libm function f that takes two double
533 arguments and returns a double.
535 The error reporting follows these rules, which are designed to do
536 the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
537 platforms.
539 - a NaN result from non-NaN inputs causes ValueError to be raised
540 - an infinite result from finite inputs causes OverflowError to be
541 raised.
542 - if the result is finite and errno == EDOM then ValueError is
543 raised
544 - if the result is finite and nonzero and errno == ERANGE then
545 OverflowError is raised
547 The last rule is used to catch overflow on platforms which follow
548 C89 but for which HUGE_VAL is not an infinity.
550 For most two-argument functions (copysign, fmod, hypot, atan2)
551 these rules are enough to ensure that Python's functions behave as
552 specified in 'Annex F' of the C99 standard, with the 'invalid' and
553 'divide-by-zero' floating-point exceptions mapping to Python's
554 ValueError and the 'overflow' floating-point exception mapping to
555 OverflowError.
558 static PyObject *
559 math_2(PyObject *args, double (*func) (double, double), char *funcname)
561 PyObject *ox, *oy;
562 double x, y, r;
563 if (! PyArg_UnpackTuple(args, funcname, 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;
569 errno = 0;
570 PyFPE_START_PROTECT("in math_2", return 0);
571 r = (*func)(x, y);
572 PyFPE_END_PROTECT(r);
573 if (Py_IS_NAN(r)) {
574 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
575 errno = EDOM;
576 else
577 errno = 0;
579 else if (Py_IS_INFINITY(r)) {
580 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
581 errno = ERANGE;
582 else
583 errno = 0;
585 if (errno && is_error(r))
586 return NULL;
587 else
588 return PyFloat_FromDouble(r);
591 #define FUNC1(funcname, func, can_overflow, docstring) \
592 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
593 return math_1(args, func, can_overflow); \
595 PyDoc_STRVAR(math_##funcname##_doc, docstring);
597 #define FUNC1A(funcname, func, docstring) \
598 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
599 return math_1a(args, func); \
601 PyDoc_STRVAR(math_##funcname##_doc, docstring);
603 #define FUNC2(funcname, func, docstring) \
604 static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
605 return math_2(args, func, #funcname); \
607 PyDoc_STRVAR(math_##funcname##_doc, docstring);
609 FUNC1(acos, acos, 0,
610 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
611 FUNC1(acosh, acosh, 0,
612 "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
613 FUNC1(asin, asin, 0,
614 "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
615 FUNC1(asinh, asinh, 0,
616 "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
617 FUNC1(atan, atan, 0,
618 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
619 FUNC2(atan2, m_atan2,
620 "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
621 "Unlike atan(y/x), the signs of both x and y are considered.")
622 FUNC1(atanh, atanh, 0,
623 "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
624 FUNC1(ceil, ceil, 0,
625 "ceil(x)\n\nReturn the ceiling of x as a float.\n"
626 "This is the smallest integral value >= x.")
627 FUNC2(copysign, copysign,
628 "copysign(x,y)\n\nReturn x with the sign of y.")
629 FUNC1(cos, cos, 0,
630 "cos(x)\n\nReturn the cosine of x (measured in radians).")
631 FUNC1(cosh, cosh, 1,
632 "cosh(x)\n\nReturn the hyperbolic cosine of x.")
633 FUNC1(exp, exp, 1,
634 "exp(x)\n\nReturn e raised to the power of x.")
635 FUNC1(fabs, fabs, 0,
636 "fabs(x)\n\nReturn the absolute value of the float x.")
637 FUNC1(floor, floor, 0,
638 "floor(x)\n\nReturn the floor of x as a float.\n"
639 "This is the largest integral value <= x.")
640 FUNC1A(gamma, m_tgamma,
641 "gamma(x)\n\nGamma function at x.")
642 FUNC1(log1p, log1p, 1,
643 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
644 The result is computed in a way which is accurate for x near zero.")
645 FUNC1(sin, sin, 0,
646 "sin(x)\n\nReturn the sine of x (measured in radians).")
647 FUNC1(sinh, sinh, 1,
648 "sinh(x)\n\nReturn the hyperbolic sine of x.")
649 FUNC1(sqrt, sqrt, 0,
650 "sqrt(x)\n\nReturn the square root of x.")
651 FUNC1(tan, tan, 0,
652 "tan(x)\n\nReturn the tangent of x (measured in radians).")
653 FUNC1(tanh, tanh, 0,
654 "tanh(x)\n\nReturn the hyperbolic tangent of x.")
656 /* Precision summation function as msum() by Raymond Hettinger in
657 <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
658 enhanced with the exact partials sum and roundoff from Mark
659 Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
660 See those links for more details, proofs and other references.
662 Note 1: IEEE 754R floating point semantics are assumed,
663 but the current implementation does not re-establish special
664 value semantics across iterations (i.e. handling -Inf + Inf).
666 Note 2: No provision is made for intermediate overflow handling;
667 therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
668 sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
669 overflow of the first partial sum.
671 Note 3: The intermediate values lo, yr, and hi are declared volatile so
672 aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
673 Also, the volatile declaration forces the values to be stored in memory as
674 regular doubles instead of extended long precision (80-bit) values. This
675 prevents double rounding because any addition or subtraction of two doubles
676 can be resolved exactly into double-sized hi and lo values. As long as the
677 hi value gets forced into a double before yr and lo are computed, the extra
678 bits in downstream extended precision operations (x87 for example) will be
679 exactly zero and therefore can be losslessly stored back into a double,
680 thereby preventing double rounding.
682 Note 4: A similar implementation is in Modules/cmathmodule.c.
683 Be sure to update both when making changes.
685 Note 5: The signature of math.fsum() differs from __builtin__.sum()
686 because the start argument doesn't make sense in the context of
687 accurate summation. Since the partials table is collapsed before
688 returning a result, sum(seq2, start=sum(seq1)) may not equal the
689 accurate result returned by sum(itertools.chain(seq1, seq2)).
692 #define NUM_PARTIALS 32 /* initial partials array size, on stack */
694 /* Extend the partials array p[] by doubling its size. */
695 static int /* non-zero on error */
696 _fsum_realloc(double **p_ptr, Py_ssize_t n,
697 double *ps, Py_ssize_t *m_ptr)
699 void *v = NULL;
700 Py_ssize_t m = *m_ptr;
702 m += m; /* double */
703 if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
704 double *p = *p_ptr;
705 if (p == ps) {
706 v = PyMem_Malloc(sizeof(double) * m);
707 if (v != NULL)
708 memcpy(v, ps, sizeof(double) * n);
710 else
711 v = PyMem_Realloc(p, sizeof(double) * m);
713 if (v == NULL) { /* size overflow or no memory */
714 PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
715 return 1;
717 *p_ptr = (double*) v;
718 *m_ptr = m;
719 return 0;
722 /* Full precision summation of a sequence of floats.
724 def msum(iterable):
725 partials = [] # sorted, non-overlapping partial sums
726 for x in iterable:
727 i = 0
728 for y in partials:
729 if abs(x) < abs(y):
730 x, y = y, x
731 hi = x + y
732 lo = y - (hi - x)
733 if lo:
734 partials[i] = lo
735 i += 1
736 x = hi
737 partials[i:] = [x]
738 return sum_exact(partials)
740 Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
741 are exactly equal to x+y. The inner loop applies hi/lo summation to each
742 partial so that the list of partial sums remains exact.
744 Sum_exact() adds the partial sums exactly and correctly rounds the final
745 result (using the round-half-to-even rule). The items in partials remain
746 non-zero, non-special, non-overlapping and strictly increasing in
747 magnitude, but possibly not all having the same sign.
749 Depends on IEEE 754 arithmetic guarantees and half-even rounding.
752 static PyObject*
753 math_fsum(PyObject *self, PyObject *seq)
755 PyObject *item, *iter, *sum = NULL;
756 Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
757 double x, y, t, ps[NUM_PARTIALS], *p = ps;
758 double xsave, special_sum = 0.0, inf_sum = 0.0;
759 volatile double hi, yr, lo;
761 iter = PyObject_GetIter(seq);
762 if (iter == NULL)
763 return NULL;
765 PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
767 for(;;) { /* for x in iterable */
768 assert(0 <= n && n <= m);
769 assert((m == NUM_PARTIALS && p == ps) ||
770 (m > NUM_PARTIALS && p != NULL));
772 item = PyIter_Next(iter);
773 if (item == NULL) {
774 if (PyErr_Occurred())
775 goto _fsum_error;
776 break;
778 x = PyFloat_AsDouble(item);
779 Py_DECREF(item);
780 if (PyErr_Occurred())
781 goto _fsum_error;
783 xsave = x;
784 for (i = j = 0; j < n; j++) { /* for y in partials */
785 y = p[j];
786 if (fabs(x) < fabs(y)) {
787 t = x; x = y; y = t;
789 hi = x + y;
790 yr = hi - x;
791 lo = y - yr;
792 if (lo != 0.0)
793 p[i++] = lo;
794 x = hi;
797 n = i; /* ps[i:] = [x] */
798 if (x != 0.0) {
799 if (! Py_IS_FINITE(x)) {
800 /* a nonfinite x could arise either as
801 a result of intermediate overflow, or
802 as a result of a nan or inf in the
803 summands */
804 if (Py_IS_FINITE(xsave)) {
805 PyErr_SetString(PyExc_OverflowError,
806 "intermediate overflow in fsum");
807 goto _fsum_error;
809 if (Py_IS_INFINITY(xsave))
810 inf_sum += xsave;
811 special_sum += xsave;
812 /* reset partials */
813 n = 0;
815 else if (n >= m && _fsum_realloc(&p, n, ps, &m))
816 goto _fsum_error;
817 else
818 p[n++] = x;
822 if (special_sum != 0.0) {
823 if (Py_IS_NAN(inf_sum))
824 PyErr_SetString(PyExc_ValueError,
825 "-inf + inf in fsum");
826 else
827 sum = PyFloat_FromDouble(special_sum);
828 goto _fsum_error;
831 hi = 0.0;
832 if (n > 0) {
833 hi = p[--n];
834 /* sum_exact(ps, hi) from the top, stop when the sum becomes
835 inexact. */
836 while (n > 0) {
837 x = hi;
838 y = p[--n];
839 assert(fabs(y) < fabs(x));
840 hi = x + y;
841 yr = hi - x;
842 lo = y - yr;
843 if (lo != 0.0)
844 break;
846 /* Make half-even rounding work across multiple partials.
847 Needed so that sum([1e-16, 1, 1e16]) will round-up the last
848 digit to two instead of down to zero (the 1e-16 makes the 1
849 slightly closer to two). With a potential 1 ULP rounding
850 error fixed-up, math.fsum() can guarantee commutativity. */
851 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
852 (lo > 0.0 && p[n-1] > 0.0))) {
853 y = lo * 2.0;
854 x = hi + y;
855 yr = x - hi;
856 if (y == yr)
857 hi = x;
860 sum = PyFloat_FromDouble(hi);
862 _fsum_error:
863 PyFPE_END_PROTECT(hi)
864 Py_DECREF(iter);
865 if (p != ps)
866 PyMem_Free(p);
867 return sum;
870 #undef NUM_PARTIALS
872 PyDoc_STRVAR(math_fsum_doc,
873 "sum(iterable)\n\n\
874 Return an accurate floating point sum of values in the iterable.\n\
875 Assumes IEEE-754 floating point arithmetic.");
877 static PyObject *
878 math_factorial(PyObject *self, PyObject *arg)
880 long i, x;
881 PyObject *result, *iobj, *newresult;
883 if (PyFloat_Check(arg)) {
884 double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
885 if (dx != floor(dx)) {
886 PyErr_SetString(PyExc_ValueError,
887 "factorial() only accepts integral values");
888 return NULL;
892 x = PyInt_AsLong(arg);
893 if (x == -1 && PyErr_Occurred())
894 return NULL;
895 if (x < 0) {
896 PyErr_SetString(PyExc_ValueError,
897 "factorial() not defined for negative values");
898 return NULL;
901 result = (PyObject *)PyInt_FromLong(1);
902 if (result == NULL)
903 return NULL;
904 for (i=1 ; i<=x ; i++) {
905 iobj = (PyObject *)PyInt_FromLong(i);
906 if (iobj == NULL)
907 goto error;
908 newresult = PyNumber_Multiply(result, iobj);
909 Py_DECREF(iobj);
910 if (newresult == NULL)
911 goto error;
912 Py_DECREF(result);
913 result = newresult;
915 return result;
917 error:
918 Py_DECREF(result);
919 return NULL;
922 PyDoc_STRVAR(math_factorial_doc,
923 "factorial(x) -> Integral\n"
924 "\n"
925 "Find x!. Raise a ValueError if x is negative or non-integral.");
927 static PyObject *
928 math_trunc(PyObject *self, PyObject *number)
930 return PyObject_CallMethod(number, "__trunc__", NULL);
933 PyDoc_STRVAR(math_trunc_doc,
934 "trunc(x:Real) -> Integral\n"
935 "\n"
936 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
938 static PyObject *
939 math_frexp(PyObject *self, PyObject *arg)
941 int i;
942 double x = PyFloat_AsDouble(arg);
943 if (x == -1.0 && PyErr_Occurred())
944 return NULL;
945 /* deal with special cases directly, to sidestep platform
946 differences */
947 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
948 i = 0;
950 else {
951 PyFPE_START_PROTECT("in math_frexp", return 0);
952 x = frexp(x, &i);
953 PyFPE_END_PROTECT(x);
955 return Py_BuildValue("(di)", x, i);
958 PyDoc_STRVAR(math_frexp_doc,
959 "frexp(x)\n"
960 "\n"
961 "Return the mantissa and exponent of x, as pair (m, e).\n"
962 "m is a float and e is an int, such that x = m * 2.**e.\n"
963 "If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
965 static PyObject *
966 math_ldexp(PyObject *self, PyObject *args)
968 double x, r;
969 PyObject *oexp;
970 long exp;
971 if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
972 return NULL;
974 if (PyLong_Check(oexp)) {
975 /* on overflow, replace exponent with either LONG_MAX
976 or LONG_MIN, depending on the sign. */
977 exp = PyLong_AsLong(oexp);
978 if (exp == -1 && PyErr_Occurred()) {
979 if (PyErr_ExceptionMatches(PyExc_OverflowError)) {
980 if (Py_SIZE(oexp) < 0) {
981 exp = LONG_MIN;
983 else {
984 exp = LONG_MAX;
986 PyErr_Clear();
988 else {
989 /* propagate any unexpected exception */
990 return NULL;
994 else if (PyInt_Check(oexp)) {
995 exp = PyInt_AS_LONG(oexp);
997 else {
998 PyErr_SetString(PyExc_TypeError,
999 "Expected an int or long as second argument "
1000 "to ldexp.");
1001 return NULL;
1004 if (x == 0. || !Py_IS_FINITE(x)) {
1005 /* NaNs, zeros and infinities are returned unchanged */
1006 r = x;
1007 errno = 0;
1008 } else if (exp > INT_MAX) {
1009 /* overflow */
1010 r = copysign(Py_HUGE_VAL, x);
1011 errno = ERANGE;
1012 } else if (exp < INT_MIN) {
1013 /* underflow to +-0 */
1014 r = copysign(0., x);
1015 errno = 0;
1016 } else {
1017 errno = 0;
1018 PyFPE_START_PROTECT("in math_ldexp", return 0);
1019 r = ldexp(x, (int)exp);
1020 PyFPE_END_PROTECT(r);
1021 if (Py_IS_INFINITY(r))
1022 errno = ERANGE;
1025 if (errno && is_error(r))
1026 return NULL;
1027 return PyFloat_FromDouble(r);
1030 PyDoc_STRVAR(math_ldexp_doc,
1031 "ldexp(x, i) -> x * (2**i)");
1033 static PyObject *
1034 math_modf(PyObject *self, PyObject *arg)
1036 double y, x = PyFloat_AsDouble(arg);
1037 if (x == -1.0 && PyErr_Occurred())
1038 return NULL;
1039 /* some platforms don't do the right thing for NaNs and
1040 infinities, so we take care of special cases directly. */
1041 if (!Py_IS_FINITE(x)) {
1042 if (Py_IS_INFINITY(x))
1043 return Py_BuildValue("(dd)", copysign(0., x), x);
1044 else if (Py_IS_NAN(x))
1045 return Py_BuildValue("(dd)", x, x);
1048 errno = 0;
1049 PyFPE_START_PROTECT("in math_modf", return 0);
1050 x = modf(x, &y);
1051 PyFPE_END_PROTECT(x);
1052 return Py_BuildValue("(dd)", x, y);
1055 PyDoc_STRVAR(math_modf_doc,
1056 "modf(x)\n"
1057 "\n"
1058 "Return the fractional and integer parts of x. Both results carry the sign\n"
1059 "of x and are floats.");
1061 /* A decent logarithm is easy to compute even for huge longs, but libm can't
1062 do that by itself -- loghelper can. func is log or log10, and name is
1063 "log" or "log10". Note that overflow isn't possible: a long can contain
1064 no more than INT_MAX * SHIFT bits, so has value certainly less than
1065 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
1066 small enough to fit in an IEEE single. log and log10 are even smaller.
1069 static PyObject*
1070 loghelper(PyObject* arg, double (*func)(double), char *funcname)
1072 /* If it is long, do it ourselves. */
1073 if (PyLong_Check(arg)) {
1074 double x;
1075 int e;
1076 x = _PyLong_AsScaledDouble(arg, &e);
1077 if (x <= 0.0) {
1078 PyErr_SetString(PyExc_ValueError,
1079 "math domain error");
1080 return NULL;
1082 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
1083 log(x) + log(2) * e * PyLong_SHIFT.
1084 CAUTION: e*PyLong_SHIFT may overflow using int arithmetic,
1085 so force use of double. */
1086 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
1087 return PyFloat_FromDouble(x);
1090 /* Else let libm handle it by itself. */
1091 return math_1(arg, func, 0);
1094 static PyObject *
1095 math_log(PyObject *self, PyObject *args)
1097 PyObject *arg;
1098 PyObject *base = NULL;
1099 PyObject *num, *den;
1100 PyObject *ans;
1102 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
1103 return NULL;
1105 num = loghelper(arg, m_log, "log");
1106 if (num == NULL || base == NULL)
1107 return num;
1109 den = loghelper(base, m_log, "log");
1110 if (den == NULL) {
1111 Py_DECREF(num);
1112 return NULL;
1115 ans = PyNumber_Divide(num, den);
1116 Py_DECREF(num);
1117 Py_DECREF(den);
1118 return ans;
1121 PyDoc_STRVAR(math_log_doc,
1122 "log(x[, base]) -> the logarithm of x to the given base.\n\
1123 If the base not specified, returns the natural logarithm (base e) of x.");
1125 static PyObject *
1126 math_log10(PyObject *self, PyObject *arg)
1128 return loghelper(arg, m_log10, "log10");
1131 PyDoc_STRVAR(math_log10_doc,
1132 "log10(x) -> the base 10 logarithm of x.");
1134 static PyObject *
1135 math_fmod(PyObject *self, PyObject *args)
1137 PyObject *ox, *oy;
1138 double r, x, y;
1139 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
1140 return NULL;
1141 x = PyFloat_AsDouble(ox);
1142 y = PyFloat_AsDouble(oy);
1143 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1144 return NULL;
1145 /* fmod(x, +/-Inf) returns x for finite x. */
1146 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1147 return PyFloat_FromDouble(x);
1148 errno = 0;
1149 PyFPE_START_PROTECT("in math_fmod", return 0);
1150 r = fmod(x, y);
1151 PyFPE_END_PROTECT(r);
1152 if (Py_IS_NAN(r)) {
1153 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1154 errno = EDOM;
1155 else
1156 errno = 0;
1158 if (errno && is_error(r))
1159 return NULL;
1160 else
1161 return PyFloat_FromDouble(r);
1164 PyDoc_STRVAR(math_fmod_doc,
1165 "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
1166 " x % y may differ.");
1168 static PyObject *
1169 math_hypot(PyObject *self, PyObject *args)
1171 PyObject *ox, *oy;
1172 double r, x, y;
1173 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
1174 return NULL;
1175 x = PyFloat_AsDouble(ox);
1176 y = PyFloat_AsDouble(oy);
1177 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1178 return NULL;
1179 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1180 if (Py_IS_INFINITY(x))
1181 return PyFloat_FromDouble(fabs(x));
1182 if (Py_IS_INFINITY(y))
1183 return PyFloat_FromDouble(fabs(y));
1184 errno = 0;
1185 PyFPE_START_PROTECT("in math_hypot", return 0);
1186 r = hypot(x, y);
1187 PyFPE_END_PROTECT(r);
1188 if (Py_IS_NAN(r)) {
1189 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1190 errno = EDOM;
1191 else
1192 errno = 0;
1194 else if (Py_IS_INFINITY(r)) {
1195 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1196 errno = ERANGE;
1197 else
1198 errno = 0;
1200 if (errno && is_error(r))
1201 return NULL;
1202 else
1203 return PyFloat_FromDouble(r);
1206 PyDoc_STRVAR(math_hypot_doc,
1207 "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
1209 /* pow can't use math_2, but needs its own wrapper: the problem is
1210 that an infinite result can arise either as a result of overflow
1211 (in which case OverflowError should be raised) or as a result of
1212 e.g. 0.**-5. (for which ValueError needs to be raised.)
1215 static PyObject *
1216 math_pow(PyObject *self, PyObject *args)
1218 PyObject *ox, *oy;
1219 double r, x, y;
1220 int odd_y;
1222 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
1223 return NULL;
1224 x = PyFloat_AsDouble(ox);
1225 y = PyFloat_AsDouble(oy);
1226 if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1227 return NULL;
1229 /* deal directly with IEEE specials, to cope with problems on various
1230 platforms whose semantics don't exactly match C99 */
1231 r = 0.; /* silence compiler warning */
1232 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
1233 errno = 0;
1234 if (Py_IS_NAN(x))
1235 r = y == 0. ? 1. : x; /* NaN**0 = 1 */
1236 else if (Py_IS_NAN(y))
1237 r = x == 1. ? 1. : y; /* 1**NaN = 1 */
1238 else if (Py_IS_INFINITY(x)) {
1239 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1240 if (y > 0.)
1241 r = odd_y ? x : fabs(x);
1242 else if (y == 0.)
1243 r = 1.;
1244 else /* y < 0. */
1245 r = odd_y ? copysign(0., x) : 0.;
1247 else if (Py_IS_INFINITY(y)) {
1248 if (fabs(x) == 1.0)
1249 r = 1.;
1250 else if (y > 0. && fabs(x) > 1.0)
1251 r = y;
1252 else if (y < 0. && fabs(x) < 1.0) {
1253 r = -y; /* result is +inf */
1254 if (x == 0.) /* 0**-inf: divide-by-zero */
1255 errno = EDOM;
1257 else
1258 r = 0.;
1261 else {
1262 /* let libm handle finite**finite */
1263 errno = 0;
1264 PyFPE_START_PROTECT("in math_pow", return 0);
1265 r = pow(x, y);
1266 PyFPE_END_PROTECT(r);
1267 /* a NaN result should arise only from (-ve)**(finite
1268 non-integer); in this case we want to raise ValueError. */
1269 if (!Py_IS_FINITE(r)) {
1270 if (Py_IS_NAN(r)) {
1271 errno = EDOM;
1274 an infinite result here arises either from:
1275 (A) (+/-0.)**negative (-> divide-by-zero)
1276 (B) overflow of x**y with x and y finite
1278 else if (Py_IS_INFINITY(r)) {
1279 if (x == 0.)
1280 errno = EDOM;
1281 else
1282 errno = ERANGE;
1287 if (errno && is_error(r))
1288 return NULL;
1289 else
1290 return PyFloat_FromDouble(r);
1293 PyDoc_STRVAR(math_pow_doc,
1294 "pow(x,y)\n\nReturn x**y (x to the power of y).");
1296 static const double degToRad = Py_MATH_PI / 180.0;
1297 static const double radToDeg = 180.0 / Py_MATH_PI;
1299 static PyObject *
1300 math_degrees(PyObject *self, PyObject *arg)
1302 double x = PyFloat_AsDouble(arg);
1303 if (x == -1.0 && PyErr_Occurred())
1304 return NULL;
1305 return PyFloat_FromDouble(x * radToDeg);
1308 PyDoc_STRVAR(math_degrees_doc,
1309 "degrees(x) -> converts angle x from radians to degrees");
1311 static PyObject *
1312 math_radians(PyObject *self, PyObject *arg)
1314 double x = PyFloat_AsDouble(arg);
1315 if (x == -1.0 && PyErr_Occurred())
1316 return NULL;
1317 return PyFloat_FromDouble(x * degToRad);
1320 PyDoc_STRVAR(math_radians_doc,
1321 "radians(x) -> converts angle x from degrees to radians");
1323 static PyObject *
1324 math_isnan(PyObject *self, PyObject *arg)
1326 double x = PyFloat_AsDouble(arg);
1327 if (x == -1.0 && PyErr_Occurred())
1328 return NULL;
1329 return PyBool_FromLong((long)Py_IS_NAN(x));
1332 PyDoc_STRVAR(math_isnan_doc,
1333 "isnan(x) -> bool\n\
1334 Checks if float x is not a number (NaN)");
1336 static PyObject *
1337 math_isinf(PyObject *self, PyObject *arg)
1339 double x = PyFloat_AsDouble(arg);
1340 if (x == -1.0 && PyErr_Occurred())
1341 return NULL;
1342 return PyBool_FromLong((long)Py_IS_INFINITY(x));
1345 PyDoc_STRVAR(math_isinf_doc,
1346 "isinf(x) -> bool\n\
1347 Checks if float x is infinite (positive or negative)");
1349 static PyMethodDef math_methods[] = {
1350 {"acos", math_acos, METH_O, math_acos_doc},
1351 {"acosh", math_acosh, METH_O, math_acosh_doc},
1352 {"asin", math_asin, METH_O, math_asin_doc},
1353 {"asinh", math_asinh, METH_O, math_asinh_doc},
1354 {"atan", math_atan, METH_O, math_atan_doc},
1355 {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
1356 {"atanh", math_atanh, METH_O, math_atanh_doc},
1357 {"ceil", math_ceil, METH_O, math_ceil_doc},
1358 {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
1359 {"cos", math_cos, METH_O, math_cos_doc},
1360 {"cosh", math_cosh, METH_O, math_cosh_doc},
1361 {"degrees", math_degrees, METH_O, math_degrees_doc},
1362 {"exp", math_exp, METH_O, math_exp_doc},
1363 {"fabs", math_fabs, METH_O, math_fabs_doc},
1364 {"factorial", math_factorial, METH_O, math_factorial_doc},
1365 {"floor", math_floor, METH_O, math_floor_doc},
1366 {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
1367 {"frexp", math_frexp, METH_O, math_frexp_doc},
1368 {"fsum", math_fsum, METH_O, math_fsum_doc},
1369 {"gamma", math_gamma, METH_O, math_gamma_doc},
1370 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
1371 {"isinf", math_isinf, METH_O, math_isinf_doc},
1372 {"isnan", math_isnan, METH_O, math_isnan_doc},
1373 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
1374 {"log", math_log, METH_VARARGS, math_log_doc},
1375 {"log1p", math_log1p, METH_O, math_log1p_doc},
1376 {"log10", math_log10, METH_O, math_log10_doc},
1377 {"modf", math_modf, METH_O, math_modf_doc},
1378 {"pow", math_pow, METH_VARARGS, math_pow_doc},
1379 {"radians", math_radians, METH_O, math_radians_doc},
1380 {"sin", math_sin, METH_O, math_sin_doc},
1381 {"sinh", math_sinh, METH_O, math_sinh_doc},
1382 {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
1383 {"tan", math_tan, METH_O, math_tan_doc},
1384 {"tanh", math_tanh, METH_O, math_tanh_doc},
1385 {"trunc", math_trunc, METH_O, math_trunc_doc},
1386 {NULL, NULL} /* sentinel */
1390 PyDoc_STRVAR(module_doc,
1391 "This module is always available. It provides access to the\n"
1392 "mathematical functions defined by the C standard.");
1394 PyMODINIT_FUNC
1395 initmath(void)
1397 PyObject *m;
1399 m = Py_InitModule3("math", math_methods, module_doc);
1400 if (m == NULL)
1401 goto finally;
1403 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
1404 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1406 finally:
1407 return;