## Python-2.7.3/Modules/mathmodule.c

No issues found

```   1 /* Math module -- standard C math library functions, pi and e */
2
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.
8
9 These are the "spirit of 754" rules:
10
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).
14
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).
18
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.
26
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.
30
31 And these are what Python has historically /tried/ to do (but not
32 always successfully, as platform libm behavior varies a lot):
33
34 For #1, raise OverflowError.
35
36 For #2, return a zero (with the appropriate sign if that happens by
37 accident ;-)).
38
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.
42
43 */
44
45 /*
46    In general, on an IEEE-754 platform the aim is to follow the C99
47    standard, including Annex 'F', whenever possible.  Where the
48    standard recommends raising the 'divide-by-zero' or 'invalid'
49    floating-point exceptions, Python should raise a ValueError.  Where
50    the standard recommends raising 'overflow', Python should raise an
51    OverflowError.  In all other circumstances a value should be
52    returned.
53  */
54
55 #include "Python.h"
56 #include "_math.h"
57
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
62
63 /*
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.
68 */
69
70 static const double pi = 3.141592653589793238462643383279502884197;
71 static const double sqrtpi = 1.772453850905516027298167483341145182798;
72
73 static double
74 sinpi(double x)
75 {
76     double y, r;
77     int n;
78     /* this function should only ever be called for finite arguments */
79     assert(Py_IS_FINITE(x));
80     y = fmod(fabs(x), 2.0);
81     n = (int)round(2.0*y);
82     assert(0 <= n && n <= 4);
83     switch (n) {
84     case 0:
85         r = sin(pi*y);
86         break;
87     case 1:
88         r = cos(pi*(y-0.5));
89         break;
90     case 2:
91         /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
92            -0.0 instead of 0.0 when y == 1.0. */
93         r = sin(pi*(1.0-y));
94         break;
95     case 3:
96         r = -cos(pi*(y-1.5));
97         break;
98     case 4:
99         r = sin(pi*(y-2.0));
100         break;
101     default:
102         assert(0);  /* should never get here */
103         r = -1.23e200; /* silence gcc warning */
104     }
105     return copysign(1.0, x)*r;
106 }
107
108 /* Implementation of the real gamma function.  In extensive but non-exhaustive
109    random tests, this function proved accurate to within <= 10 ulps across the
110    entire float domain.  Note that accuracy may depend on the quality of the
111    system math functions, the pow function in particular.  Special cases
112    follow C99 annex F.  The parameters and method are tailored to platforms
113    whose double format is the IEEE 754 binary64 format.
114
115    Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
116    and g=6.024680040776729583740234375; these parameters are amongst those
117    used by the Boost library.  Following Boost (again), we re-express the
118    Lanczos sum as a rational function, and compute it that way.  The
119    coefficients below were computed independently using MPFR, and have been
120    double-checked against the coefficients in the Boost source code.
121
122    For x < 0.0 we use the reflection formula.
123
124    There's one minor tweak that deserves explanation: Lanczos' formula for
125    Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5).  For many x
126    values, x+g-0.5 can be represented exactly.  However, in cases where it
127    can't be represented exactly the small error in x+g-0.5 can be magnified
128    significantly by the pow and exp calls, especially for large x.  A cheap
129    correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
130    involved in the computation of x+g-0.5 (that is, e = computed value of
131    x+g-0.5 - exact value of x+g-0.5).  Here's the proof:
132
133    Correction factor
134    -----------------
135    Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
136    double, and e is tiny.  Then:
137
138      pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
139      = pow(y, x-0.5)/exp(y) * C,
140
141    where the correction_factor C is given by
142
143      C = pow(1-e/y, x-0.5) * exp(e)
144
145    Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
146
147      C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
148
149    But y-(x-0.5) = g+e, and g+e ~ g.  So we get C ~ 1 + e*g/y, and
150
151      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
153    Note that for accuracy, when computing r*C it's better to do
154
155      r + e*g/y*r;
156
157    than
158
159      r * (1 + e*g/y);
160
161    since the addition in the latter throws away most of the bits of
162    information in e*g/y.
163 */
164
165 #define LANCZOS_N 13
166 static const double lanczos_g = 6.024680040776729583740234375;
167 static const double lanczos_g_minus_half = 5.524680040776729583740234375;
168 static const double lanczos_num_coeffs[LANCZOS_N] = {
169     23531376880.410759688572007674451636754734846804940,
170     42919803642.649098768957899047001988850926355848959,
171     35711959237.355668049440185451547166705960488635843,
172     17921034426.037209699919755754458931112671403265390,
173     6039542586.3520280050642916443072979210699388420708,
174     1439720407.3117216736632230727949123939715485786772,
175     248874557.86205415651146038641322942321632125127801,
176     31426415.585400194380614231628318205362874684987640,
177     2876370.6289353724412254090516208496135991145378768,
178     186056.26539522349504029498971604569928220784236328,
179     8071.6720023658162106380029022722506138218516325024,
180     210.82427775157934587250973392071336271166969580291,
181     2.5066282746310002701649081771338373386264310793408
182 };
183
184 /* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
185 static const double lanczos_den_coeffs[LANCZOS_N] = {
186     0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
187     13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
188
189 /* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
190 #define NGAMMA_INTEGRAL 23
191 static const double gamma_integral[NGAMMA_INTEGRAL] = {
192     1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
193     3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
194     1307674368000.0, 20922789888000.0, 355687428096000.0,
195     6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
196     51090942171709440000.0, 1124000727777607680000.0,
197 };
198
199 /* Lanczos' sum L_g(x), for positive x */
200
201 static double
202 lanczos_sum(double x)
203 {
204     double num = 0.0, den = 0.0;
205     int i;
206     assert(x > 0.0);
207     /* evaluate the rational function lanczos_sum(x).  For large
208        x, the obvious algorithm risks overflow, so we instead
209        rescale the denominator and numerator of the rational
210        function by x**(1-LANCZOS_N) and treat this as a
211        rational function in 1/x.  This also reduces the error for
212        larger x values.  The choice of cutoff point (5.0 below) is
213        somewhat arbitrary; in tests, smaller cutoff values than
214        this resulted in lower accuracy. */
215     if (x < 5.0) {
216         for (i = LANCZOS_N; --i >= 0; ) {
217             num = num * x + lanczos_num_coeffs[i];
218             den = den * x + lanczos_den_coeffs[i];
219         }
220     }
221     else {
222         for (i = 0; i < LANCZOS_N; i++) {
223             num = num / x + lanczos_num_coeffs[i];
224             den = den / x + lanczos_den_coeffs[i];
225         }
226     }
227     return num/den;
228 }
229
230 static double
231 m_tgamma(double x)
232 {
233     double absx, r, y, z, sqrtpow;
234
235     /* special cases */
236     if (!Py_IS_FINITE(x)) {
237         if (Py_IS_NAN(x) || x > 0.0)
238             return x;  /* tgamma(nan) = nan, tgamma(inf) = inf */
239         else {
240             errno = EDOM;
241             return Py_NAN;  /* tgamma(-inf) = nan, invalid */
242         }
243     }
244     if (x == 0.0) {
245         errno = EDOM;
246         return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */
247     }
248
249     /* integer arguments */
250     if (x == floor(x)) {
251         if (x < 0.0) {
252             errno = EDOM;  /* tgamma(n) = nan, invalid for */
253             return Py_NAN; /* negative integers n */
254         }
255         if (x <= NGAMMA_INTEGRAL)
256             return gamma_integral[(int)x - 1];
257     }
258     absx = fabs(x);
259
260     /* tiny arguments:  tgamma(x) ~ 1/x for x near 0 */
261     if (absx < 1e-20) {
262         r = 1.0/x;
263         if (Py_IS_INFINITY(r))
264             errno = ERANGE;
265         return r;
266     }
267
268     /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
269        x > 200, and underflows to +-0.0 for x < -200, not a negative
270        integer. */
271     if (absx > 200.0) {
272         if (x < 0.0) {
273             return 0.0/sinpi(x);
274         }
275         else {
276             errno = ERANGE;
277             return Py_HUGE_VAL;
278         }
279     }
280
281     y = absx + lanczos_g_minus_half;
282     /* compute error in sum */
283     if (absx > lanczos_g_minus_half) {
284         /* note: the correction can be foiled by an optimizing
285            compiler that (incorrectly) thinks that an expression like
286            a + b - a - b can be optimized to 0.0.  This shouldn't
287            happen in a standards-conforming compiler. */
288         double q = y - absx;
289         z = q - lanczos_g_minus_half;
290     }
291     else {
292         double q = y - lanczos_g_minus_half;
293         z = q - absx;
294     }
295     z = z * lanczos_g / y;
296     if (x < 0.0) {
297         r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
298         r -= z * r;
299         if (absx < 140.0) {
300             r /= pow(y, absx - 0.5);
301         }
302         else {
303             sqrtpow = pow(y, absx / 2.0 - 0.25);
304             r /= sqrtpow;
305             r /= sqrtpow;
306         }
307     }
308     else {
309         r = lanczos_sum(absx) / exp(y);
310         r += z * r;
311         if (absx < 140.0) {
312             r *= pow(y, absx - 0.5);
313         }
314         else {
315             sqrtpow = pow(y, absx / 2.0 - 0.25);
316             r *= sqrtpow;
317             r *= sqrtpow;
318         }
319     }
320     if (Py_IS_INFINITY(r))
321         errno = ERANGE;
322     return r;
323 }
324
325 /*
326    lgamma:  natural log of the absolute value of the Gamma function.
327    For large arguments, Lanczos' formula works extremely well here.
328 */
329
330 static double
331 m_lgamma(double x)
332 {
333     double r, absx;
334
335     /* special cases */
336     if (!Py_IS_FINITE(x)) {
337         if (Py_IS_NAN(x))
338             return x;  /* lgamma(nan) = nan */
339         else
340             return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
341     }
342
343     /* integer arguments */
344     if (x == floor(x) && x <= 2.0) {
345         if (x <= 0.0) {
346             errno = EDOM;  /* lgamma(n) = inf, divide-by-zero for */
347             return Py_HUGE_VAL; /* integers n <= 0 */
348         }
349         else {
350             return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
351         }
352     }
353
354     absx = fabs(x);
355     /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
356     if (absx < 1e-20)
357         return -log(absx);
358
359     /* Lanczos' formula */
360     if (x > 0.0) {
361         /* we could save a fraction of a ulp in accuracy by having a
362            second set of numerator coefficients for lanczos_sum that
363            absorbed the exp(-lanczos_g) term, and throwing out the
364            lanczos_g subtraction below; it's probably not worth it. */
365         r = log(lanczos_sum(x)) - lanczos_g +
366             (x-0.5)*(log(x+lanczos_g-0.5)-1);
367     }
368     else {
369         r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -
370             (log(lanczos_sum(absx)) - lanczos_g +
371              (absx-0.5)*(log(absx+lanczos_g-0.5)-1));
372     }
373     if (Py_IS_INFINITY(r))
374         errno = ERANGE;
375     return r;
376 }
377
378 /*
379    Implementations of the error function erf(x) and the complementary error
380    function erfc(x).
381
382    Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,
383    Cambridge University Press), we use a series approximation for erf for
384    small x, and a continued fraction approximation for erfc(x) for larger x;
385    combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
386    this gives us erf(x) and erfc(x) for all x.
387
388    The series expansion used is:
389
390       erf(x) = x*exp(-x*x)/sqrt(pi) * [
391                      2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
392
393    The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
394    This series converges well for smallish x, but slowly for larger x.
395
396    The continued fraction expansion used is:
397
398       erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
399                               3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
400
401    after the first term, the general term has the form:
402
403       k*(k-0.5)/(2*k+0.5 + x**2 - ...).
404
405    This expansion converges fast for larger x, but convergence becomes
406    infinitely slow as x approaches 0.0.  The (somewhat naive) continued
407    fraction evaluation algorithm used below also risks overflow for large x;
408    but for large x, erfc(x) == 0.0 to within machine precision.  (For
409    example, erfc(30.0) is approximately 2.56e-393).
410
411    Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
412    continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
413    ERFC_CONTFRAC_CUTOFF.  ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
414    numbers of terms to use for the relevant expansions.  */
415
416 #define ERF_SERIES_CUTOFF 1.5
417 #define ERF_SERIES_TERMS 25
418 #define ERFC_CONTFRAC_CUTOFF 30.0
419 #define ERFC_CONTFRAC_TERMS 50
420
421 /*
422    Error function, via power series.
423
424    Given a finite float x, return an approximation to erf(x).
425    Converges reasonably fast for small x.
426 */
427
428 static double
429 m_erf_series(double x)
430 {
431     double x2, acc, fk, result;
432     int i, saved_errno;
433
434     x2 = x * x;
435     acc = 0.0;
436     fk = (double)ERF_SERIES_TERMS + 0.5;
437     for (i = 0; i < ERF_SERIES_TERMS; i++) {
438         acc = 2.0 + x2 * acc / fk;
439         fk -= 1.0;
440     }
441     /* Make sure the exp call doesn't affect errno;
442        see m_erfc_contfrac for more. */
443     saved_errno = errno;
444     result = acc * x * exp(-x2) / sqrtpi;
445     errno = saved_errno;
446     return result;
447 }
448
449 /*
450    Complementary error function, via continued fraction expansion.
451
452    Given a positive float x, return an approximation to erfc(x).  Converges
453    reasonably fast for x large (say, x > 2.0), and should be safe from
454    overflow if x and nterms are not too large.  On an IEEE 754 machine, with x
455    <= 30.0, we're safe up to nterms = 100.  For x >= 30.0, erfc(x) is smaller
456    than the smallest representable nonzero float.  */
457
458 static double
459 m_erfc_contfrac(double x)
460 {
461     double x2, a, da, p, p_last, q, q_last, b, result;
462     int i, saved_errno;
463
464     if (x >= ERFC_CONTFRAC_CUTOFF)
465         return 0.0;
466
467     x2 = x*x;
468     a = 0.0;
469     da = 0.5;
470     p = 1.0; p_last = 0.0;
471     q = da + x2; q_last = 1.0;
472     for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
473         double temp;
474         a += da;
475         da += 2.0;
476         b = da + x2;
477         temp = p; p = b*p - a*p_last; p_last = temp;
478         temp = q; q = b*q - a*q_last; q_last = temp;
479     }
480     /* Issue #8986: On some platforms, exp sets errno on underflow to zero;
481        save the current errno value so that we can restore it later. */
482     saved_errno = errno;
483     result = p / q * x * exp(-x2) / sqrtpi;
484     errno = saved_errno;
485     return result;
486 }
487
488 /* Error function erf(x), for general x */
489
490 static double
491 m_erf(double x)
492 {
493     double absx, cf;
494
495     if (Py_IS_NAN(x))
496         return x;
497     absx = fabs(x);
498     if (absx < ERF_SERIES_CUTOFF)
499         return m_erf_series(x);
500     else {
501         cf = m_erfc_contfrac(absx);
502         return x > 0.0 ? 1.0 - cf : cf - 1.0;
503     }
504 }
505
506 /* Complementary error function erfc(x), for general x. */
507
508 static double
509 m_erfc(double x)
510 {
511     double absx, cf;
512
513     if (Py_IS_NAN(x))
514         return x;
515     absx = fabs(x);
516     if (absx < ERF_SERIES_CUTOFF)
517         return 1.0 - m_erf_series(x);
518     else {
519         cf = m_erfc_contfrac(absx);
520         return x > 0.0 ? cf : 2.0 - cf;
521     }
522 }
523
524 /*
525    wrapper for atan2 that deals directly with special cases before
526    delegating to the platform libm for the remaining cases.  This
527    is necessary to get consistent behaviour across platforms.
528    Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
529    always follow C99.
530 */
531
532 static double
533 m_atan2(double y, double x)
534 {
535     if (Py_IS_NAN(x) || Py_IS_NAN(y))
536         return Py_NAN;
537     if (Py_IS_INFINITY(y)) {
538         if (Py_IS_INFINITY(x)) {
539             if (copysign(1., x) == 1.)
540                 /* atan2(+-inf, +inf) == +-pi/4 */
541                 return copysign(0.25*Py_MATH_PI, y);
542             else
543                 /* atan2(+-inf, -inf) == +-pi*3/4 */
544                 return copysign(0.75*Py_MATH_PI, y);
545         }
546         /* atan2(+-inf, x) == +-pi/2 for finite x */
547         return copysign(0.5*Py_MATH_PI, y);
548     }
549     if (Py_IS_INFINITY(x) || y == 0.) {
550         if (copysign(1., x) == 1.)
551             /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
552             return copysign(0., y);
553         else
554             /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
555             return copysign(Py_MATH_PI, y);
556     }
557     return atan2(y, x);
558 }
559
560 /*
561     Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
562     log(-ve), log(NaN).  Here are wrappers for log and log10 that deal with
563     special values directly, passing positive non-special values through to
564     the system log/log10.
565  */
566
567 static double
568 m_log(double x)
569 {
570     if (Py_IS_FINITE(x)) {
571         if (x > 0.0)
572             return log(x);
573         errno = EDOM;
574         if (x == 0.0)
575             return -Py_HUGE_VAL; /* log(0) = -inf */
576         else
577             return Py_NAN; /* log(-ve) = nan */
578     }
579     else if (Py_IS_NAN(x))
580         return x; /* log(nan) = nan */
581     else if (x > 0.0)
582         return x; /* log(inf) = inf */
583     else {
584         errno = EDOM;
585         return Py_NAN; /* log(-inf) = nan */
586     }
587 }
588
589 static double
590 m_log10(double x)
591 {
592     if (Py_IS_FINITE(x)) {
593         if (x > 0.0)
594             return log10(x);
595         errno = EDOM;
596         if (x == 0.0)
597             return -Py_HUGE_VAL; /* log10(0) = -inf */
598         else
599             return Py_NAN; /* log10(-ve) = nan */
600     }
601     else if (Py_IS_NAN(x))
602         return x; /* log10(nan) = nan */
603     else if (x > 0.0)
604         return x; /* log10(inf) = inf */
605     else {
606         errno = EDOM;
607         return Py_NAN; /* log10(-inf) = nan */
608     }
609 }
610
611
612 /* Call is_error when errno != 0, and where x is the result libm
613  * returned.  is_error will usually set up an exception and return
614  * true (1), but may return false (0) without setting up an exception.
615  */
616 static int
617 is_error(double x)
618 {
619     int result = 1;     /* presumption of guilt */
620     assert(errno);      /* non-zero errno is a precondition for calling */
621     if (errno == EDOM)
622         PyErr_SetString(PyExc_ValueError, "math domain error");
623
624     else if (errno == ERANGE) {
625         /* ANSI C generally requires libm functions to set ERANGE
626          * on overflow, but also generally *allows* them to set
627          * ERANGE on underflow too.  There's no consistency about
628          * the latter across platforms.
629          * Alas, C99 never requires that errno be set.
630          * Here we suppress the underflow errors (libm functions
631          * should return a zero on underflow, and +- HUGE_VAL on
632          * overflow, so testing the result for zero suffices to
633          * distinguish the cases).
634          *
635          * On some platforms (Ubuntu/ia64) it seems that errno can be
636          * set to ERANGE for subnormal results that do *not* underflow
637          * to zero.  So to be safe, we'll ignore ERANGE whenever the
638          * function result is less than one in absolute value.
639          */
640         if (fabs(x) < 1.0)
641             result = 0;
642         else
643             PyErr_SetString(PyExc_OverflowError,
644                             "math range error");
645     }
646     else
647         /* Unexpected math error */
648         PyErr_SetFromErrno(PyExc_ValueError);
649     return result;
650 }
651
652 /*
653    math_1 is used to wrap a libm function f that takes a double
654    arguments and returns a double.
655
656    The error reporting follows these rules, which are designed to do
657    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
658    platforms.
659
660    - a NaN result from non-NaN inputs causes ValueError to be raised
661    - an infinite result from finite inputs causes OverflowError to be
662      raised if can_overflow is 1, or raises ValueError if can_overflow
663      is 0.
664    - if the result is finite and errno == EDOM then ValueError is
665      raised
666    - if the result is finite and nonzero and errno == ERANGE then
667      OverflowError is raised
668
669    The last rule is used to catch overflow on platforms which follow
670    C89 but for which HUGE_VAL is not an infinity.
671
672    For the majority of one-argument functions these rules are enough
673    to ensure that Python's functions behave as specified in 'Annex F'
674    of the C99 standard, with the 'invalid' and 'divide-by-zero'
675    floating-point exceptions mapping to Python's ValueError and the
676    'overflow' floating-point exception mapping to OverflowError.
677    math_1 only works for functions that don't have singularities *and*
678    the possibility of overflow; fortunately, that covers everything we
679    care about right now.
680 */
681
682 static PyObject *
683 math_1(PyObject *arg, double (*func) (double), int can_overflow)
684 {
685     double x, r;
686     x = PyFloat_AsDouble(arg);
687     if (x == -1.0 && PyErr_Occurred())
688         return NULL;
689     errno = 0;
690     PyFPE_START_PROTECT("in math_1", return 0);
691     r = (*func)(x);
692     PyFPE_END_PROTECT(r);
693     if (Py_IS_NAN(r)) {
694         if (!Py_IS_NAN(x))
695             errno = EDOM;
696         else
697             errno = 0;
698     }
699     else if (Py_IS_INFINITY(r)) {
700         if (Py_IS_FINITE(x))
701             errno = can_overflow ? ERANGE : EDOM;
702         else
703             errno = 0;
704     }
705     if (errno && is_error(r))
706         return NULL;
707     else
708         return PyFloat_FromDouble(r);
709 }
710
711 /* variant of math_1, to be used when the function being wrapped is known to
712    set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
713    errno = ERANGE for overflow). */
714
715 static PyObject *
716 math_1a(PyObject *arg, double (*func) (double))
717 {
718     double x, r;
719     x = PyFloat_AsDouble(arg);
720     if (x == -1.0 && PyErr_Occurred())
721         return NULL;
722     errno = 0;
723     PyFPE_START_PROTECT("in math_1a", return 0);
724     r = (*func)(x);
725     PyFPE_END_PROTECT(r);
726     if (errno && is_error(r))
727         return NULL;
728     return PyFloat_FromDouble(r);
729 }
730
731 /*
732    math_2 is used to wrap a libm function f that takes two double
733    arguments and returns a double.
734
735    The error reporting follows these rules, which are designed to do
736    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
737    platforms.
738
739    - a NaN result from non-NaN inputs causes ValueError to be raised
740    - an infinite result from finite inputs causes OverflowError to be
741      raised.
742    - if the result is finite and errno == EDOM then ValueError is
743      raised
744    - if the result is finite and nonzero and errno == ERANGE then
745      OverflowError is raised
746
747    The last rule is used to catch overflow on platforms which follow
748    C89 but for which HUGE_VAL is not an infinity.
749
750    For most two-argument functions (copysign, fmod, hypot, atan2)
751    these rules are enough to ensure that Python's functions behave as
752    specified in 'Annex F' of the C99 standard, with the 'invalid' and
753    'divide-by-zero' floating-point exceptions mapping to Python's
754    ValueError and the 'overflow' floating-point exception mapping to
755    OverflowError.
756 */
757
758 static PyObject *
759 math_2(PyObject *args, double (*func) (double, double), char *funcname)
760 {
761     PyObject *ox, *oy;
762     double x, y, r;
763     if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
764         return NULL;
765     x = PyFloat_AsDouble(ox);
766     y = PyFloat_AsDouble(oy);
767     if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
768         return NULL;
769     errno = 0;
770     PyFPE_START_PROTECT("in math_2", return 0);
771     r = (*func)(x, y);
772     PyFPE_END_PROTECT(r);
773     if (Py_IS_NAN(r)) {
774         if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
775             errno = EDOM;
776         else
777             errno = 0;
778     }
779     else if (Py_IS_INFINITY(r)) {
780         if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
781             errno = ERANGE;
782         else
783             errno = 0;
784     }
785     if (errno && is_error(r))
786         return NULL;
787     else
788         return PyFloat_FromDouble(r);
789 }
790
791 #define FUNC1(funcname, func, can_overflow, docstring)                  \
792     static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
793         return math_1(args, func, can_overflow);                            \
794     }\
795     PyDoc_STRVAR(math_##funcname##_doc, docstring);
796
797 #define FUNC1A(funcname, func, docstring)                               \
798     static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
799         return math_1a(args, func);                                     \
800     }\
801     PyDoc_STRVAR(math_##funcname##_doc, docstring);
802
803 #define FUNC2(funcname, func, docstring) \
804     static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
805         return math_2(args, func, #funcname); \
806     }\
807     PyDoc_STRVAR(math_##funcname##_doc, docstring);
808
809 FUNC1(acos, acos, 0,
810       "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
811 FUNC1(acosh, m_acosh, 0,
812       "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
813 FUNC1(asin, asin, 0,
814       "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
815 FUNC1(asinh, m_asinh, 0,
816       "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
817 FUNC1(atan, atan, 0,
818       "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
819 FUNC2(atan2, m_atan2,
820       "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
821       "Unlike atan(y/x), the signs of both x and y are considered.")
822 FUNC1(atanh, m_atanh, 0,
823       "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
824 FUNC1(ceil, ceil, 0,
825       "ceil(x)\n\nReturn the ceiling of x as a float.\n"
826       "This is the smallest integral value >= x.")
827 FUNC2(copysign, copysign,
828       "copysign(x, y)\n\nReturn x with the sign of y.")
829 FUNC1(cos, cos, 0,
830       "cos(x)\n\nReturn the cosine of x (measured in radians).")
831 FUNC1(cosh, cosh, 1,
832       "cosh(x)\n\nReturn the hyperbolic cosine of x.")
833 FUNC1A(erf, m_erf,
834        "erf(x)\n\nError function at x.")
835 FUNC1A(erfc, m_erfc,
836        "erfc(x)\n\nComplementary error function at x.")
837 FUNC1(exp, exp, 1,
838       "exp(x)\n\nReturn e raised to the power of x.")
839 FUNC1(expm1, m_expm1, 1,
840       "expm1(x)\n\nReturn exp(x)-1.\n"
841       "This function avoids the loss of precision involved in the direct "
842       "evaluation of exp(x)-1 for small x.")
843 FUNC1(fabs, fabs, 0,
844       "fabs(x)\n\nReturn the absolute value of the float x.")
845 FUNC1(floor, floor, 0,
846       "floor(x)\n\nReturn the floor of x as a float.\n"
847       "This is the largest integral value <= x.")
848 FUNC1A(gamma, m_tgamma,
849       "gamma(x)\n\nGamma function at x.")
850 FUNC1A(lgamma, m_lgamma,
851       "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
852 FUNC1(log1p, m_log1p, 1,
853       "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
854       "The result is computed in a way which is accurate for x near zero.")
855 FUNC1(sin, sin, 0,
856       "sin(x)\n\nReturn the sine of x (measured in radians).")
857 FUNC1(sinh, sinh, 1,
858       "sinh(x)\n\nReturn the hyperbolic sine of x.")
859 FUNC1(sqrt, sqrt, 0,
860       "sqrt(x)\n\nReturn the square root of x.")
861 FUNC1(tan, tan, 0,
862       "tan(x)\n\nReturn the tangent of x (measured in radians).")
863 FUNC1(tanh, tanh, 0,
864       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
865
866 /* Precision summation function as msum() by Raymond Hettinger in
867    <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
868    enhanced with the exact partials sum and roundoff from Mark
869    Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
870    See those links for more details, proofs and other references.
871
872    Note 1: IEEE 754R floating point semantics are assumed,
873    but the current implementation does not re-establish special
874    value semantics across iterations (i.e. handling -Inf + Inf).
875
876    Note 2:  No provision is made for intermediate overflow handling;
877    therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
878    sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
879    overflow of the first partial sum.
880
881    Note 3: The intermediate values lo, yr, and hi are declared volatile so
882    aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
883    Also, the volatile declaration forces the values to be stored in memory as
884    regular doubles instead of extended long precision (80-bit) values.  This
885    prevents double rounding because any addition or subtraction of two doubles
886    can be resolved exactly into double-sized hi and lo values.  As long as the
887    hi value gets forced into a double before yr and lo are computed, the extra
888    bits in downstream extended precision operations (x87 for example) will be
889    exactly zero and therefore can be losslessly stored back into a double,
890    thereby preventing double rounding.
891
892    Note 4: A similar implementation is in Modules/cmathmodule.c.
893    Be sure to update both when making changes.
894
895    Note 5: The signature of math.fsum() differs from __builtin__.sum()
896    because the start argument doesn't make sense in the context of
897    accurate summation.  Since the partials table is collapsed before
898    returning a result, sum(seq2, start=sum(seq1)) may not equal the
899    accurate result returned by sum(itertools.chain(seq1, seq2)).
900 */
901
902 #define NUM_PARTIALS  32  /* initial partials array size, on stack */
903
904 /* Extend the partials array p[] by doubling its size. */
905 static int                          /* non-zero on error */
906 _fsum_realloc(double **p_ptr, Py_ssize_t  n,
907              double  *ps,    Py_ssize_t *m_ptr)
908 {
909     void *v = NULL;
910     Py_ssize_t m = *m_ptr;
911
912     m += m;  /* double */
913     if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
914         double *p = *p_ptr;
915         if (p == ps) {
916             v = PyMem_Malloc(sizeof(double) * m);
917             if (v != NULL)
918                 memcpy(v, ps, sizeof(double) * n);
919         }
920         else
921             v = PyMem_Realloc(p, sizeof(double) * m);
922     }
923     if (v == NULL) {        /* size overflow or no memory */
924         PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
925         return 1;
926     }
927     *p_ptr = (double*) v;
928     *m_ptr = m;
929     return 0;
930 }
931
932 /* Full precision summation of a sequence of floats.
933
934    def msum(iterable):
935        partials = []  # sorted, non-overlapping partial sums
936        for x in iterable:
937            i = 0
938            for y in partials:
939                if abs(x) < abs(y):
940                    x, y = y, x
941                hi = x + y
942                lo = y - (hi - x)
943                if lo:
944                    partials[i] = lo
945                    i += 1
946                x = hi
947            partials[i:] = [x]
948        return sum_exact(partials)
949
950    Rounded x+y stored in hi with the roundoff stored in lo.  Together hi+lo
951    are exactly equal to x+y.  The inner loop applies hi/lo summation to each
952    partial so that the list of partial sums remains exact.
953
954    Sum_exact() adds the partial sums exactly and correctly rounds the final
955    result (using the round-half-to-even rule).  The items in partials remain
956    non-zero, non-special, non-overlapping and strictly increasing in
957    magnitude, but possibly not all having the same sign.
958
959    Depends on IEEE 754 arithmetic guarantees and half-even rounding.
960 */
961
962 static PyObject*
963 math_fsum(PyObject *self, PyObject *seq)
964 {
965     PyObject *item, *iter, *sum = NULL;
966     Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
967     double x, y, t, ps[NUM_PARTIALS], *p = ps;
968     double xsave, special_sum = 0.0, inf_sum = 0.0;
969     volatile double hi, yr, lo;
970
971     iter = PyObject_GetIter(seq);
972     if (iter == NULL)
973         return NULL;
974
975     PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
976
977     for(;;) {           /* for x in iterable */
978         assert(0 <= n && n <= m);
979         assert((m == NUM_PARTIALS && p == ps) ||
980                (m >  NUM_PARTIALS && p != NULL));
981
982         item = PyIter_Next(iter);
983         if (item == NULL) {
984             if (PyErr_Occurred())
985                 goto _fsum_error;
986             break;
987         }
988         x = PyFloat_AsDouble(item);
989         Py_DECREF(item);
990         if (PyErr_Occurred())
991             goto _fsum_error;
992
993         xsave = x;
994         for (i = j = 0; j < n; j++) {       /* for y in partials */
995             y = p[j];
996             if (fabs(x) < fabs(y)) {
997                 t = x; x = y; y = t;
998             }
999             hi = x + y;
1000             yr = hi - x;
1001             lo = y - yr;
1002             if (lo != 0.0)
1003                 p[i++] = lo;
1004             x = hi;
1005         }
1006
1007         n = i;                              /* ps[i:] = [x] */
1008         if (x != 0.0) {
1009             if (! Py_IS_FINITE(x)) {
1010                 /* a nonfinite x could arise either as
1011                    a result of intermediate overflow, or
1012                    as a result of a nan or inf in the
1013                    summands */
1014                 if (Py_IS_FINITE(xsave)) {
1015                     PyErr_SetString(PyExc_OverflowError,
1016                           "intermediate overflow in fsum");
1017                     goto _fsum_error;
1018                 }
1019                 if (Py_IS_INFINITY(xsave))
1020                     inf_sum += xsave;
1021                 special_sum += xsave;
1022                 /* reset partials */
1023                 n = 0;
1024             }
1025             else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1026                 goto _fsum_error;
1027             else
1028                 p[n++] = x;
1029         }
1030     }
1031
1032     if (special_sum != 0.0) {
1033         if (Py_IS_NAN(inf_sum))
1034             PyErr_SetString(PyExc_ValueError,
1035                             "-inf + inf in fsum");
1036         else
1037             sum = PyFloat_FromDouble(special_sum);
1038         goto _fsum_error;
1039     }
1040
1041     hi = 0.0;
1042     if (n > 0) {
1043         hi = p[--n];
1044         /* sum_exact(ps, hi) from the top, stop when the sum becomes
1045            inexact. */
1046         while (n > 0) {
1047             x = hi;
1048             y = p[--n];
1049             assert(fabs(y) < fabs(x));
1050             hi = x + y;
1051             yr = hi - x;
1052             lo = y - yr;
1053             if (lo != 0.0)
1054                 break;
1055         }
1056         /* Make half-even rounding work across multiple partials.
1057            Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1058            digit to two instead of down to zero (the 1e-16 makes the 1
1059            slightly closer to two).  With a potential 1 ULP rounding
1060            error fixed-up, math.fsum() can guarantee commutativity. */
1061         if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1062                       (lo > 0.0 && p[n-1] > 0.0))) {
1063             y = lo * 2.0;
1064             x = hi + y;
1065             yr = x - hi;
1066             if (y == yr)
1067                 hi = x;
1068         }
1069     }
1070     sum = PyFloat_FromDouble(hi);
1071
1072 _fsum_error:
1073     PyFPE_END_PROTECT(hi)
1074     Py_DECREF(iter);
1075     if (p != ps)
1076         PyMem_Free(p);
1077     return sum;
1078 }
1079
1080 #undef NUM_PARTIALS
1081
1082 PyDoc_STRVAR(math_fsum_doc,
1083 "fsum(iterable)\n\n\
1084 Return an accurate floating point sum of values in the iterable.\n\
1085 Assumes IEEE-754 floating point arithmetic.");
1086
1087 static PyObject *
1088 math_factorial(PyObject *self, PyObject *arg)
1089 {
1090     long i, x;
1091     PyObject *result, *iobj, *newresult;
1092
1093     if (PyFloat_Check(arg)) {
1094         PyObject *lx;
1095         double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1096         if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1097             PyErr_SetString(PyExc_ValueError,
1098                 "factorial() only accepts integral values");
1099             return NULL;
1100         }
1101         lx = PyLong_FromDouble(dx);
1102         if (lx == NULL)
1103             return NULL;
1104         x = PyLong_AsLong(lx);
1105         Py_DECREF(lx);
1106     }
1107     else
1108         x = PyInt_AsLong(arg);
1109
1110     if (x == -1 && PyErr_Occurred())
1111         return NULL;
1112     if (x < 0) {
1113         PyErr_SetString(PyExc_ValueError,
1114             "factorial() not defined for negative values");
1115         return NULL;
1116     }
1117
1118     result = (PyObject *)PyInt_FromLong(1);
1119     if (result == NULL)
1120         return NULL;
1121     for (i=1 ; i<=x ; i++) {
1122         iobj = (PyObject *)PyInt_FromLong(i);
1123         if (iobj == NULL)
1124             goto error;
1125         newresult = PyNumber_Multiply(result, iobj);
1126         Py_DECREF(iobj);
1127         if (newresult == NULL)
1128             goto error;
1129         Py_DECREF(result);
1130         result = newresult;
1131     }
1132     return result;
1133
1134 error:
1135     Py_DECREF(result);
1136     return NULL;
1137 }
1138
1139 PyDoc_STRVAR(math_factorial_doc,
1140 "factorial(x) -> Integral\n"
1141 "\n"
1142 "Find x!. Raise a ValueError if x is negative or non-integral.");
1143
1144 static PyObject *
1145 math_trunc(PyObject *self, PyObject *number)
1146 {
1147     return PyObject_CallMethod(number, "__trunc__", NULL);
1148 }
1149
1150 PyDoc_STRVAR(math_trunc_doc,
1151 "trunc(x:Real) -> Integral\n"
1152 "\n"
1153 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
1154
1155 static PyObject *
1156 math_frexp(PyObject *self, PyObject *arg)
1157 {
1158     int i;
1159     double x = PyFloat_AsDouble(arg);
1160     if (x == -1.0 && PyErr_Occurred())
1161         return NULL;
1162     /* deal with special cases directly, to sidestep platform
1163        differences */
1164     if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
1165         i = 0;
1166     }
1167     else {
1168         PyFPE_START_PROTECT("in math_frexp", return 0);
1169         x = frexp(x, &i);
1170         PyFPE_END_PROTECT(x);
1171     }
1172     return Py_BuildValue("(di)", x, i);
1173 }
1174
1175 PyDoc_STRVAR(math_frexp_doc,
1176 "frexp(x)\n"
1177 "\n"
1178 "Return the mantissa and exponent of x, as pair (m, e).\n"
1179 "m is a float and e is an int, such that x = m * 2.**e.\n"
1180 "If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.");
1181
1182 static PyObject *
1183 math_ldexp(PyObject *self, PyObject *args)
1184 {
1185     double x, r;
1186     PyObject *oexp;
1187     long exp;
1188     int overflow;
1189     if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
1190         return NULL;
1191
1192     if (PyLong_Check(oexp) || PyInt_Check(oexp)) {
1193         /* on overflow, replace exponent with either LONG_MAX
1194            or LONG_MIN, depending on the sign. */
1195         exp = PyLong_AsLongAndOverflow(oexp, &overflow);
1196         if (exp == -1 && PyErr_Occurred())
1197             return NULL;
1198         if (overflow)
1199             exp = overflow < 0 ? LONG_MIN : LONG_MAX;
1200     }
1201     else {
1202         PyErr_SetString(PyExc_TypeError,
1203                         "Expected an int or long as second argument "
1204                         "to ldexp.");
1205         return NULL;
1206     }
1207
1208     if (x == 0. || !Py_IS_FINITE(x)) {
1209         /* NaNs, zeros and infinities are returned unchanged */
1210         r = x;
1211         errno = 0;
1212     } else if (exp > INT_MAX) {
1213         /* overflow */
1214         r = copysign(Py_HUGE_VAL, x);
1215         errno = ERANGE;
1216     } else if (exp < INT_MIN) {
1217         /* underflow to +-0 */
1218         r = copysign(0., x);
1219         errno = 0;
1220     } else {
1221         errno = 0;
1222         PyFPE_START_PROTECT("in math_ldexp", return 0);
1223         r = ldexp(x, (int)exp);
1224         PyFPE_END_PROTECT(r);
1225         if (Py_IS_INFINITY(r))
1226             errno = ERANGE;
1227     }
1228
1229     if (errno && is_error(r))
1230         return NULL;
1231     return PyFloat_FromDouble(r);
1232 }
1233
1234 PyDoc_STRVAR(math_ldexp_doc,
1235 "ldexp(x, i)\n\n\
1236 Return x * (2**i).");
1237
1238 static PyObject *
1239 math_modf(PyObject *self, PyObject *arg)
1240 {
1241     double y, x = PyFloat_AsDouble(arg);
1242     if (x == -1.0 && PyErr_Occurred())
1243         return NULL;
1244     /* some platforms don't do the right thing for NaNs and
1245        infinities, so we take care of special cases directly. */
1246     if (!Py_IS_FINITE(x)) {
1247         if (Py_IS_INFINITY(x))
1248             return Py_BuildValue("(dd)", copysign(0., x), x);
1249         else if (Py_IS_NAN(x))
1250             return Py_BuildValue("(dd)", x, x);
1251     }
1252
1253     errno = 0;
1254     PyFPE_START_PROTECT("in math_modf", return 0);
1255     x = modf(x, &y);
1256     PyFPE_END_PROTECT(x);
1257     return Py_BuildValue("(dd)", x, y);
1258 }
1259
1260 PyDoc_STRVAR(math_modf_doc,
1261 "modf(x)\n"
1262 "\n"
1263 "Return the fractional and integer parts of x.  Both results carry the sign\n"
1264 "of x and are floats.");
1265
1266 /* A decent logarithm is easy to compute even for huge longs, but libm can't
1267    do that by itself -- loghelper can.  func is log or log10, and name is
1268    "log" or "log10".  Note that overflow of the result isn't possible: a long
1269    can contain no more than INT_MAX * SHIFT bits, so has value certainly less
1270    than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
1271    small enough to fit in an IEEE single.  log and log10 are even smaller.
1272    However, intermediate overflow is possible for a long if the number of bits
1273    in that long is larger than PY_SSIZE_T_MAX. */
1274
1275 static PyObject*
1276 loghelper(PyObject* arg, double (*func)(double), char *funcname)
1277 {
1278     /* If it is long, do it ourselves. */
1279     if (PyLong_Check(arg)) {
1280         double x;
1281         Py_ssize_t e;
1282         x = _PyLong_Frexp((PyLongObject *)arg, &e);
1283         if (x == -1.0 && PyErr_Occurred())
1284             return NULL;
1285         if (x <= 0.0) {
1286             PyErr_SetString(PyExc_ValueError,
1287                             "math domain error");
1288             return NULL;
1289         }
1290         /* Special case for log(1), to make sure we get an
1291            exact result there. */
1292         if (e == 1 && x == 0.5)
1293             return PyFloat_FromDouble(0.0);
1294         /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
1295         x = func(x) + func(2.0) * e;
1296         return PyFloat_FromDouble(x);
1297     }
1298
1299     /* Else let libm handle it by itself. */
1300     return math_1(arg, func, 0);
1301 }
1302
1303 static PyObject *
1304 math_log(PyObject *self, PyObject *args)
1305 {
1306     PyObject *arg;
1307     PyObject *base = NULL;
1308     PyObject *num, *den;
1309     PyObject *ans;
1310
1311     if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
1312         return NULL;
1313
1314     num = loghelper(arg, m_log, "log");
1315     if (num == NULL || base == NULL)
1316         return num;
1317
1318     den = loghelper(base, m_log, "log");
1319     if (den == NULL) {
1320         Py_DECREF(num);
1321         return NULL;
1322     }
1323
1324     ans = PyNumber_Divide(num, den);
1325     Py_DECREF(num);
1326     Py_DECREF(den);
1327     return ans;
1328 }
1329
1330 PyDoc_STRVAR(math_log_doc,
1331 "log(x[, base])\n\n\
1332 Return the logarithm of x to the given base.\n\
1333 If the base not specified, returns the natural logarithm (base e) of x.");
1334
1335 static PyObject *
1336 math_log10(PyObject *self, PyObject *arg)
1337 {
1338     return loghelper(arg, m_log10, "log10");
1339 }
1340
1341 PyDoc_STRVAR(math_log10_doc,
1342 "log10(x)\n\nReturn the base 10 logarithm of x.");
1343
1344 static PyObject *
1345 math_fmod(PyObject *self, PyObject *args)
1346 {
1347     PyObject *ox, *oy;
1348     double r, x, y;
1349     if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
1350         return NULL;
1351     x = PyFloat_AsDouble(ox);
1352     y = PyFloat_AsDouble(oy);
1353     if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1354         return NULL;
1355     /* fmod(x, +/-Inf) returns x for finite x. */
1356     if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1357         return PyFloat_FromDouble(x);
1358     errno = 0;
1359     PyFPE_START_PROTECT("in math_fmod", return 0);
1360     r = fmod(x, y);
1361     PyFPE_END_PROTECT(r);
1362     if (Py_IS_NAN(r)) {
1363         if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1364             errno = EDOM;
1365         else
1366             errno = 0;
1367     }
1368     if (errno && is_error(r))
1369         return NULL;
1370     else
1371         return PyFloat_FromDouble(r);
1372 }
1373
1374 PyDoc_STRVAR(math_fmod_doc,
1375 "fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
1376 "  x % y may differ.");
1377
1378 static PyObject *
1379 math_hypot(PyObject *self, PyObject *args)
1380 {
1381     PyObject *ox, *oy;
1382     double r, x, y;
1383     if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
1384         return NULL;
1385     x = PyFloat_AsDouble(ox);
1386     y = PyFloat_AsDouble(oy);
1387     if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1388         return NULL;
1389     /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1390     if (Py_IS_INFINITY(x))
1391         return PyFloat_FromDouble(fabs(x));
1392     if (Py_IS_INFINITY(y))
1393         return PyFloat_FromDouble(fabs(y));
1394     errno = 0;
1395     PyFPE_START_PROTECT("in math_hypot", return 0);
1396     r = hypot(x, y);
1397     PyFPE_END_PROTECT(r);
1398     if (Py_IS_NAN(r)) {
1399         if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1400             errno = EDOM;
1401         else
1402             errno = 0;
1403     }
1404     else if (Py_IS_INFINITY(r)) {
1405         if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1406             errno = ERANGE;
1407         else
1408             errno = 0;
1409     }
1410     if (errno && is_error(r))
1411         return NULL;
1412     else
1413         return PyFloat_FromDouble(r);
1414 }
1415
1416 PyDoc_STRVAR(math_hypot_doc,
1417 "hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
1418
1419 /* pow can't use math_2, but needs its own wrapper: the problem is
1420    that an infinite result can arise either as a result of overflow
1421    (in which case OverflowError should be raised) or as a result of
1422    e.g. 0.**-5. (for which ValueError needs to be raised.)
1423 */
1424
1425 static PyObject *
1426 math_pow(PyObject *self, PyObject *args)
1427 {
1428     PyObject *ox, *oy;
1429     double r, x, y;
1430     int odd_y;
1431
1432     if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
1433         return NULL;
1434     x = PyFloat_AsDouble(ox);
1435     y = PyFloat_AsDouble(oy);
1436     if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1437         return NULL;
1438
1439     /* deal directly with IEEE specials, to cope with problems on various
1440        platforms whose semantics don't exactly match C99 */
1441     r = 0.; /* silence compiler warning */
1442     if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
1443         errno = 0;
1444         if (Py_IS_NAN(x))
1445             r = y == 0. ? 1. : x; /* NaN**0 = 1 */
1446         else if (Py_IS_NAN(y))
1447             r = x == 1. ? 1. : y; /* 1**NaN = 1 */
1448         else if (Py_IS_INFINITY(x)) {
1449             odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1450             if (y > 0.)
1451                 r = odd_y ? x : fabs(x);
1452             else if (y == 0.)
1453                 r = 1.;
1454             else /* y < 0. */
1455                 r = odd_y ? copysign(0., x) : 0.;
1456         }
1457         else if (Py_IS_INFINITY(y)) {
1458             if (fabs(x) == 1.0)
1459                 r = 1.;
1460             else if (y > 0. && fabs(x) > 1.0)
1461                 r = y;
1462             else if (y < 0. && fabs(x) < 1.0) {
1463                 r = -y; /* result is +inf */
1464                 if (x == 0.) /* 0**-inf: divide-by-zero */
1465                     errno = EDOM;
1466             }
1467             else
1468                 r = 0.;
1469         }
1470     }
1471     else {
1472         /* let libm handle finite**finite */
1473         errno = 0;
1474         PyFPE_START_PROTECT("in math_pow", return 0);
1475         r = pow(x, y);
1476         PyFPE_END_PROTECT(r);
1477         /* a NaN result should arise only from (-ve)**(finite
1478            non-integer); in this case we want to raise ValueError. */
1479         if (!Py_IS_FINITE(r)) {
1480             if (Py_IS_NAN(r)) {
1481                 errno = EDOM;
1482             }
1483             /*
1484                an infinite result here arises either from:
1485                (A) (+/-0.)**negative (-> divide-by-zero)
1486                (B) overflow of x**y with x and y finite
1487             */
1488             else if (Py_IS_INFINITY(r)) {
1489                 if (x == 0.)
1490                     errno = EDOM;
1491                 else
1492                     errno = ERANGE;
1493             }
1494         }
1495     }
1496
1497     if (errno && is_error(r))
1498         return NULL;
1499     else
1500         return PyFloat_FromDouble(r);
1501 }
1502
1503 PyDoc_STRVAR(math_pow_doc,
1504 "pow(x, y)\n\nReturn x**y (x to the power of y).");
1505
1506 static const double degToRad = Py_MATH_PI / 180.0;
1507 static const double radToDeg = 180.0 / Py_MATH_PI;
1508
1509 static PyObject *
1510 math_degrees(PyObject *self, PyObject *arg)
1511 {
1512     double x = PyFloat_AsDouble(arg);
1513     if (x == -1.0 && PyErr_Occurred())
1514         return NULL;
1515     return PyFloat_FromDouble(x * radToDeg);
1516 }
1517
1518 PyDoc_STRVAR(math_degrees_doc,
1519 "degrees(x)\n\n\
1520 Convert angle x from radians to degrees.");
1521
1522 static PyObject *
1523 math_radians(PyObject *self, PyObject *arg)
1524 {
1525     double x = PyFloat_AsDouble(arg);
1526     if (x == -1.0 && PyErr_Occurred())
1527         return NULL;
1528     return PyFloat_FromDouble(x * degToRad);
1529 }
1530
1533 Convert angle x from degrees to radians.");
1534
1535 static PyObject *
1536 math_isnan(PyObject *self, PyObject *arg)
1537 {
1538     double x = PyFloat_AsDouble(arg);
1539     if (x == -1.0 && PyErr_Occurred())
1540         return NULL;
1541     return PyBool_FromLong((long)Py_IS_NAN(x));
1542 }
1543
1544 PyDoc_STRVAR(math_isnan_doc,
1545 "isnan(x) -> bool\n\n\
1546 Check if float x is not a number (NaN).");
1547
1548 static PyObject *
1549 math_isinf(PyObject *self, PyObject *arg)
1550 {
1551     double x = PyFloat_AsDouble(arg);
1552     if (x == -1.0 && PyErr_Occurred())
1553         return NULL;
1554     return PyBool_FromLong((long)Py_IS_INFINITY(x));
1555 }
1556
1557 PyDoc_STRVAR(math_isinf_doc,
1558 "isinf(x) -> bool\n\n\
1559 Check if float x is infinite (positive or negative).");
1560
1561 static PyMethodDef math_methods[] = {
1562     {"acos",            math_acos,      METH_O,         math_acos_doc},
1563     {"acosh",           math_acosh,     METH_O,         math_acosh_doc},
1564     {"asin",            math_asin,      METH_O,         math_asin_doc},
1565     {"asinh",           math_asinh,     METH_O,         math_asinh_doc},
1566     {"atan",            math_atan,      METH_O,         math_atan_doc},
1567     {"atan2",           math_atan2,     METH_VARARGS,   math_atan2_doc},
1568     {"atanh",           math_atanh,     METH_O,         math_atanh_doc},
1569     {"ceil",            math_ceil,      METH_O,         math_ceil_doc},
1570     {"copysign",        math_copysign,  METH_VARARGS,   math_copysign_doc},
1571     {"cos",             math_cos,       METH_O,         math_cos_doc},
1572     {"cosh",            math_cosh,      METH_O,         math_cosh_doc},
1573     {"degrees",         math_degrees,   METH_O,         math_degrees_doc},
1574     {"erf",             math_erf,       METH_O,         math_erf_doc},
1575     {"erfc",            math_erfc,      METH_O,         math_erfc_doc},
1576     {"exp",             math_exp,       METH_O,         math_exp_doc},
1577     {"expm1",           math_expm1,     METH_O,         math_expm1_doc},
1578     {"fabs",            math_fabs,      METH_O,         math_fabs_doc},
1579     {"factorial",       math_factorial, METH_O,         math_factorial_doc},
1580     {"floor",           math_floor,     METH_O,         math_floor_doc},
1581     {"fmod",            math_fmod,      METH_VARARGS,   math_fmod_doc},
1582     {"frexp",           math_frexp,     METH_O,         math_frexp_doc},
1583     {"fsum",            math_fsum,      METH_O,         math_fsum_doc},
1584     {"gamma",           math_gamma,     METH_O,         math_gamma_doc},
1585     {"hypot",           math_hypot,     METH_VARARGS,   math_hypot_doc},
1586     {"isinf",           math_isinf,     METH_O,         math_isinf_doc},
1587     {"isnan",           math_isnan,     METH_O,         math_isnan_doc},
1588     {"ldexp",           math_ldexp,     METH_VARARGS,   math_ldexp_doc},
1589     {"lgamma",          math_lgamma,    METH_O,         math_lgamma_doc},
1590     {"log",             math_log,       METH_VARARGS,   math_log_doc},
1591     {"log1p",           math_log1p,     METH_O,         math_log1p_doc},
1592     {"log10",           math_log10,     METH_O,         math_log10_doc},
1593     {"modf",            math_modf,      METH_O,         math_modf_doc},
1594     {"pow",             math_pow,       METH_VARARGS,   math_pow_doc},
1596     {"sin",             math_sin,       METH_O,         math_sin_doc},
1597     {"sinh",            math_sinh,      METH_O,         math_sinh_doc},
1598     {"sqrt",            math_sqrt,      METH_O,         math_sqrt_doc},
1599     {"tan",             math_tan,       METH_O,         math_tan_doc},
1600     {"tanh",            math_tanh,      METH_O,         math_tanh_doc},
1601     {"trunc",           math_trunc,     METH_O,         math_trunc_doc},
1602     {NULL,              NULL}           /* sentinel */
1603 };
1604
1605
1606 PyDoc_STRVAR(module_doc,
1607 "This module is always available.  It provides access to the\n"
1608 "mathematical functions defined by the C standard.");
1609
1610 PyMODINIT_FUNC
1611 initmath(void)
1612 {
1613     PyObject *m;
1614
1615     m = Py_InitModule3("math", math_methods, module_doc);
1616     if (m == NULL)
1617         goto finally;
1618
1619     PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
1620     PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1621
1622     finally:
1623     return;
1624 }
```