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
1531 PyDoc_STRVAR(math_radians_doc,
1532 "radians(x)\n\n\
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},
1595 {"radians", math_radians, METH_O, math_radians_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 }