Line data Source code
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 : /*
59 : sin(pi*x), giving accurate results for all finite x (especially x
60 : integral or close to an integer). This is here for use in the
61 : reflection formula for the gamma function. It conforms to IEEE
62 : 754-2008 for finite arguments, but not for infinities or nans.
63 : */
64 :
65 : static const double pi = 3.141592653589793238462643383279502884197;
66 : static const double sqrtpi = 1.772453850905516027298167483341145182798;
67 : static const double logpi = 1.144729885849400174143427351353058711647;
68 :
69 : static double
70 0 : sinpi(double x)
71 : {
72 : double y, r;
73 : int n;
74 : /* this function should only ever be called for finite arguments */
75 : assert(Py_IS_FINITE(x));
76 0 : y = fmod(fabs(x), 2.0);
77 0 : n = (int)round(2.0*y);
78 : assert(0 <= n && n <= 4);
79 0 : switch (n) {
80 : case 0:
81 0 : r = sin(pi*y);
82 0 : break;
83 : case 1:
84 0 : r = cos(pi*(y-0.5));
85 0 : break;
86 : case 2:
87 : /* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give
88 : -0.0 instead of 0.0 when y == 1.0. */
89 0 : r = sin(pi*(1.0-y));
90 0 : break;
91 : case 3:
92 0 : r = -cos(pi*(y-1.5));
93 0 : break;
94 : case 4:
95 0 : r = sin(pi*(y-2.0));
96 0 : break;
97 : default:
98 : assert(0); /* should never get here */
99 0 : r = -1.23e200; /* silence gcc warning */
100 : }
101 0 : return copysign(1.0, x)*r;
102 : }
103 :
104 : /* Implementation of the real gamma function. In extensive but non-exhaustive
105 : random tests, this function proved accurate to within <= 10 ulps across the
106 : entire float domain. Note that accuracy may depend on the quality of the
107 : system math functions, the pow function in particular. Special cases
108 : follow C99 annex F. The parameters and method are tailored to platforms
109 : whose double format is the IEEE 754 binary64 format.
110 :
111 : Method: for x > 0.0 we use the Lanczos approximation with parameters N=13
112 : and g=6.024680040776729583740234375; these parameters are amongst those
113 : used by the Boost library. Following Boost (again), we re-express the
114 : Lanczos sum as a rational function, and compute it that way. The
115 : coefficients below were computed independently using MPFR, and have been
116 : double-checked against the coefficients in the Boost source code.
117 :
118 : For x < 0.0 we use the reflection formula.
119 :
120 : There's one minor tweak that deserves explanation: Lanczos' formula for
121 : Gamma(x) involves computing pow(x+g-0.5, x-0.5) / exp(x+g-0.5). For many x
122 : values, x+g-0.5 can be represented exactly. However, in cases where it
123 : can't be represented exactly the small error in x+g-0.5 can be magnified
124 : significantly by the pow and exp calls, especially for large x. A cheap
125 : correction is to multiply by (1 + e*g/(x+g-0.5)), where e is the error
126 : involved in the computation of x+g-0.5 (that is, e = computed value of
127 : x+g-0.5 - exact value of x+g-0.5). Here's the proof:
128 :
129 : Correction factor
130 : -----------------
131 : Write x+g-0.5 = y-e, where y is exactly representable as an IEEE 754
132 : double, and e is tiny. Then:
133 :
134 : pow(x+g-0.5,x-0.5)/exp(x+g-0.5) = pow(y-e, x-0.5)/exp(y-e)
135 : = pow(y, x-0.5)/exp(y) * C,
136 :
137 : where the correction_factor C is given by
138 :
139 : C = pow(1-e/y, x-0.5) * exp(e)
140 :
141 : Since e is tiny, pow(1-e/y, x-0.5) ~ 1-(x-0.5)*e/y, and exp(x) ~ 1+e, so:
142 :
143 : C ~ (1-(x-0.5)*e/y) * (1+e) ~ 1 + e*(y-(x-0.5))/y
144 :
145 : But y-(x-0.5) = g+e, and g+e ~ g. So we get C ~ 1 + e*g/y, and
146 :
147 : pow(x+g-0.5,x-0.5)/exp(x+g-0.5) ~ pow(y, x-0.5)/exp(y) * (1 + e*g/y),
148 :
149 : Note that for accuracy, when computing r*C it's better to do
150 :
151 : r + e*g/y*r;
152 :
153 : than
154 :
155 : r * (1 + e*g/y);
156 :
157 : since the addition in the latter throws away most of the bits of
158 : information in e*g/y.
159 : */
160 :
161 : #define LANCZOS_N 13
162 : static const double lanczos_g = 6.024680040776729583740234375;
163 : static const double lanczos_g_minus_half = 5.524680040776729583740234375;
164 : static const double lanczos_num_coeffs[LANCZOS_N] = {
165 : 23531376880.410759688572007674451636754734846804940,
166 : 42919803642.649098768957899047001988850926355848959,
167 : 35711959237.355668049440185451547166705960488635843,
168 : 17921034426.037209699919755754458931112671403265390,
169 : 6039542586.3520280050642916443072979210699388420708,
170 : 1439720407.3117216736632230727949123939715485786772,
171 : 248874557.86205415651146038641322942321632125127801,
172 : 31426415.585400194380614231628318205362874684987640,
173 : 2876370.6289353724412254090516208496135991145378768,
174 : 186056.26539522349504029498971604569928220784236328,
175 : 8071.6720023658162106380029022722506138218516325024,
176 : 210.82427775157934587250973392071336271166969580291,
177 : 2.5066282746310002701649081771338373386264310793408
178 : };
179 :
180 : /* denominator is x*(x+1)*...*(x+LANCZOS_N-2) */
181 : static const double lanczos_den_coeffs[LANCZOS_N] = {
182 : 0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
183 : 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};
184 :
185 : /* gamma values for small positive integers, 1 though NGAMMA_INTEGRAL */
186 : #define NGAMMA_INTEGRAL 23
187 : static const double gamma_integral[NGAMMA_INTEGRAL] = {
188 : 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
189 : 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
190 : 1307674368000.0, 20922789888000.0, 355687428096000.0,
191 : 6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
192 : 51090942171709440000.0, 1124000727777607680000.0,
193 : };
194 :
195 : /* Lanczos' sum L_g(x), for positive x */
196 :
197 : static double
198 0 : lanczos_sum(double x)
199 : {
200 0 : double num = 0.0, den = 0.0;
201 : int i;
202 : assert(x > 0.0);
203 : /* evaluate the rational function lanczos_sum(x). For large
204 : x, the obvious algorithm risks overflow, so we instead
205 : rescale the denominator and numerator of the rational
206 : function by x**(1-LANCZOS_N) and treat this as a
207 : rational function in 1/x. This also reduces the error for
208 : larger x values. The choice of cutoff point (5.0 below) is
209 : somewhat arbitrary; in tests, smaller cutoff values than
210 : this resulted in lower accuracy. */
211 0 : if (x < 5.0) {
212 0 : for (i = LANCZOS_N; --i >= 0; ) {
213 0 : num = num * x + lanczos_num_coeffs[i];
214 0 : den = den * x + lanczos_den_coeffs[i];
215 : }
216 : }
217 : else {
218 0 : for (i = 0; i < LANCZOS_N; i++) {
219 0 : num = num / x + lanczos_num_coeffs[i];
220 0 : den = den / x + lanczos_den_coeffs[i];
221 : }
222 : }
223 0 : return num/den;
224 : }
225 :
226 : static double
227 0 : m_tgamma(double x)
228 : {
229 : double absx, r, y, z, sqrtpow;
230 :
231 : /* special cases */
232 0 : if (!Py_IS_FINITE(x)) {
233 0 : if (Py_IS_NAN(x) || x > 0.0)
234 0 : return x; /* tgamma(nan) = nan, tgamma(inf) = inf */
235 : else {
236 0 : errno = EDOM;
237 0 : return Py_NAN; /* tgamma(-inf) = nan, invalid */
238 : }
239 : }
240 0 : if (x == 0.0) {
241 0 : errno = EDOM;
242 : /* tgamma(+-0.0) = +-inf, divide-by-zero */
243 0 : return copysign(Py_HUGE_VAL, x);
244 : }
245 :
246 : /* integer arguments */
247 0 : if (x == floor(x)) {
248 0 : if (x < 0.0) {
249 0 : errno = EDOM; /* tgamma(n) = nan, invalid for */
250 0 : return Py_NAN; /* negative integers n */
251 : }
252 0 : if (x <= NGAMMA_INTEGRAL)
253 0 : return gamma_integral[(int)x - 1];
254 : }
255 0 : absx = fabs(x);
256 :
257 : /* tiny arguments: tgamma(x) ~ 1/x for x near 0 */
258 0 : if (absx < 1e-20) {
259 0 : r = 1.0/x;
260 0 : if (Py_IS_INFINITY(r))
261 0 : errno = ERANGE;
262 0 : return r;
263 : }
264 :
265 : /* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for
266 : x > 200, and underflows to +-0.0 for x < -200, not a negative
267 : integer. */
268 0 : if (absx > 200.0) {
269 0 : if (x < 0.0) {
270 0 : return 0.0/sinpi(x);
271 : }
272 : else {
273 0 : errno = ERANGE;
274 0 : return Py_HUGE_VAL;
275 : }
276 : }
277 :
278 0 : y = absx + lanczos_g_minus_half;
279 : /* compute error in sum */
280 0 : if (absx > lanczos_g_minus_half) {
281 : /* note: the correction can be foiled by an optimizing
282 : compiler that (incorrectly) thinks that an expression like
283 : a + b - a - b can be optimized to 0.0. This shouldn't
284 : happen in a standards-conforming compiler. */
285 0 : double q = y - absx;
286 0 : z = q - lanczos_g_minus_half;
287 : }
288 : else {
289 0 : double q = y - lanczos_g_minus_half;
290 0 : z = q - absx;
291 : }
292 0 : z = z * lanczos_g / y;
293 0 : if (x < 0.0) {
294 0 : r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);
295 0 : r -= z * r;
296 0 : if (absx < 140.0) {
297 0 : r /= pow(y, absx - 0.5);
298 : }
299 : else {
300 0 : sqrtpow = pow(y, absx / 2.0 - 0.25);
301 0 : r /= sqrtpow;
302 0 : r /= sqrtpow;
303 : }
304 : }
305 : else {
306 0 : r = lanczos_sum(absx) / exp(y);
307 0 : r += z * r;
308 0 : if (absx < 140.0) {
309 0 : r *= pow(y, absx - 0.5);
310 : }
311 : else {
312 0 : sqrtpow = pow(y, absx / 2.0 - 0.25);
313 0 : r *= sqrtpow;
314 0 : r *= sqrtpow;
315 : }
316 : }
317 0 : if (Py_IS_INFINITY(r))
318 0 : errno = ERANGE;
319 0 : return r;
320 : }
321 :
322 : /*
323 : lgamma: natural log of the absolute value of the Gamma function.
324 : For large arguments, Lanczos' formula works extremely well here.
325 : */
326 :
327 : static double
328 0 : m_lgamma(double x)
329 : {
330 : double r, absx;
331 :
332 : /* special cases */
333 0 : if (!Py_IS_FINITE(x)) {
334 0 : if (Py_IS_NAN(x))
335 0 : return x; /* lgamma(nan) = nan */
336 : else
337 0 : return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
338 : }
339 :
340 : /* integer arguments */
341 0 : if (x == floor(x) && x <= 2.0) {
342 0 : if (x <= 0.0) {
343 0 : errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
344 0 : return Py_HUGE_VAL; /* integers n <= 0 */
345 : }
346 : else {
347 0 : return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
348 : }
349 : }
350 :
351 0 : absx = fabs(x);
352 : /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
353 0 : if (absx < 1e-20)
354 0 : return -log(absx);
355 :
356 : /* Lanczos' formula. We could save a fraction of a ulp in accuracy by
357 : having a second set of numerator coefficients for lanczos_sum that
358 : absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
359 : subtraction below; it's probably not worth it. */
360 0 : r = log(lanczos_sum(absx)) - lanczos_g;
361 0 : r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
362 0 : if (x < 0.0)
363 : /* Use reflection formula to get value for negative x. */
364 0 : r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;
365 0 : if (Py_IS_INFINITY(r))
366 0 : errno = ERANGE;
367 0 : return r;
368 : }
369 :
370 : /*
371 : Implementations of the error function erf(x) and the complementary error
372 : function erfc(x).
373 :
374 : Method: following 'Numerical Recipes' by Flannery, Press et. al. (2nd ed.,
375 : Cambridge University Press), we use a series approximation for erf for
376 : small x, and a continued fraction approximation for erfc(x) for larger x;
377 : combined with the relations erf(-x) = -erf(x) and erfc(x) = 1.0 - erf(x),
378 : this gives us erf(x) and erfc(x) for all x.
379 :
380 : The series expansion used is:
381 :
382 : erf(x) = x*exp(-x*x)/sqrt(pi) * [
383 : 2/1 + 4/3 x**2 + 8/15 x**4 + 16/105 x**6 + ...]
384 :
385 : The coefficient of x**(2k-2) here is 4**k*factorial(k)/factorial(2*k).
386 : This series converges well for smallish x, but slowly for larger x.
387 :
388 : The continued fraction expansion used is:
389 :
390 : erfc(x) = x*exp(-x*x)/sqrt(pi) * [1/(0.5 + x**2 -) 0.5/(2.5 + x**2 - )
391 : 3.0/(4.5 + x**2 - ) 7.5/(6.5 + x**2 - ) ...]
392 :
393 : after the first term, the general term has the form:
394 :
395 : k*(k-0.5)/(2*k+0.5 + x**2 - ...).
396 :
397 : This expansion converges fast for larger x, but convergence becomes
398 : infinitely slow as x approaches 0.0. The (somewhat naive) continued
399 : fraction evaluation algorithm used below also risks overflow for large x;
400 : but for large x, erfc(x) == 0.0 to within machine precision. (For
401 : example, erfc(30.0) is approximately 2.56e-393).
402 :
403 : Parameters: use series expansion for abs(x) < ERF_SERIES_CUTOFF and
404 : continued fraction expansion for ERF_SERIES_CUTOFF <= abs(x) <
405 : ERFC_CONTFRAC_CUTOFF. ERFC_SERIES_TERMS and ERFC_CONTFRAC_TERMS are the
406 : numbers of terms to use for the relevant expansions. */
407 :
408 : #define ERF_SERIES_CUTOFF 1.5
409 : #define ERF_SERIES_TERMS 25
410 : #define ERFC_CONTFRAC_CUTOFF 30.0
411 : #define ERFC_CONTFRAC_TERMS 50
412 :
413 : /*
414 : Error function, via power series.
415 :
416 : Given a finite float x, return an approximation to erf(x).
417 : Converges reasonably fast for small x.
418 : */
419 :
420 : static double
421 0 : m_erf_series(double x)
422 : {
423 : double x2, acc, fk, result;
424 : int i, saved_errno;
425 :
426 0 : x2 = x * x;
427 0 : acc = 0.0;
428 0 : fk = (double)ERF_SERIES_TERMS + 0.5;
429 0 : for (i = 0; i < ERF_SERIES_TERMS; i++) {
430 0 : acc = 2.0 + x2 * acc / fk;
431 0 : fk -= 1.0;
432 : }
433 : /* Make sure the exp call doesn't affect errno;
434 : see m_erfc_contfrac for more. */
435 0 : saved_errno = errno;
436 0 : result = acc * x * exp(-x2) / sqrtpi;
437 0 : errno = saved_errno;
438 0 : return result;
439 : }
440 :
441 : /*
442 : Complementary error function, via continued fraction expansion.
443 :
444 : Given a positive float x, return an approximation to erfc(x). Converges
445 : reasonably fast for x large (say, x > 2.0), and should be safe from
446 : overflow if x and nterms are not too large. On an IEEE 754 machine, with x
447 : <= 30.0, we're safe up to nterms = 100. For x >= 30.0, erfc(x) is smaller
448 : than the smallest representable nonzero float. */
449 :
450 : static double
451 0 : m_erfc_contfrac(double x)
452 : {
453 : double x2, a, da, p, p_last, q, q_last, b, result;
454 : int i, saved_errno;
455 :
456 0 : if (x >= ERFC_CONTFRAC_CUTOFF)
457 0 : return 0.0;
458 :
459 0 : x2 = x*x;
460 0 : a = 0.0;
461 0 : da = 0.5;
462 0 : p = 1.0; p_last = 0.0;
463 0 : q = da + x2; q_last = 1.0;
464 0 : for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {
465 : double temp;
466 0 : a += da;
467 0 : da += 2.0;
468 0 : b = da + x2;
469 0 : temp = p; p = b*p - a*p_last; p_last = temp;
470 0 : temp = q; q = b*q - a*q_last; q_last = temp;
471 : }
472 : /* Issue #8986: On some platforms, exp sets errno on underflow to zero;
473 : save the current errno value so that we can restore it later. */
474 0 : saved_errno = errno;
475 0 : result = p / q * x * exp(-x2) / sqrtpi;
476 0 : errno = saved_errno;
477 0 : return result;
478 : }
479 :
480 : /* Error function erf(x), for general x */
481 :
482 : static double
483 0 : m_erf(double x)
484 : {
485 : double absx, cf;
486 :
487 0 : if (Py_IS_NAN(x))
488 0 : return x;
489 0 : absx = fabs(x);
490 0 : if (absx < ERF_SERIES_CUTOFF)
491 0 : return m_erf_series(x);
492 : else {
493 0 : cf = m_erfc_contfrac(absx);
494 0 : return x > 0.0 ? 1.0 - cf : cf - 1.0;
495 : }
496 : }
497 :
498 : /* Complementary error function erfc(x), for general x. */
499 :
500 : static double
501 0 : m_erfc(double x)
502 : {
503 : double absx, cf;
504 :
505 0 : if (Py_IS_NAN(x))
506 0 : return x;
507 0 : absx = fabs(x);
508 0 : if (absx < ERF_SERIES_CUTOFF)
509 0 : return 1.0 - m_erf_series(x);
510 : else {
511 0 : cf = m_erfc_contfrac(absx);
512 0 : return x > 0.0 ? cf : 2.0 - cf;
513 : }
514 : }
515 :
516 : /*
517 : wrapper for atan2 that deals directly with special cases before
518 : delegating to the platform libm for the remaining cases. This
519 : is necessary to get consistent behaviour across platforms.
520 : Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
521 : always follow C99.
522 : */
523 :
524 : static double
525 0 : m_atan2(double y, double x)
526 : {
527 0 : if (Py_IS_NAN(x) || Py_IS_NAN(y))
528 0 : return Py_NAN;
529 0 : if (Py_IS_INFINITY(y)) {
530 0 : if (Py_IS_INFINITY(x)) {
531 0 : if (copysign(1., x) == 1.)
532 : /* atan2(+-inf, +inf) == +-pi/4 */
533 0 : return copysign(0.25*Py_MATH_PI, y);
534 : else
535 : /* atan2(+-inf, -inf) == +-pi*3/4 */
536 0 : return copysign(0.75*Py_MATH_PI, y);
537 : }
538 : /* atan2(+-inf, x) == +-pi/2 for finite x */
539 0 : return copysign(0.5*Py_MATH_PI, y);
540 : }
541 0 : if (Py_IS_INFINITY(x) || y == 0.) {
542 0 : if (copysign(1., x) == 1.)
543 : /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
544 0 : return copysign(0., y);
545 : else
546 : /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
547 0 : return copysign(Py_MATH_PI, y);
548 : }
549 0 : return atan2(y, x);
550 : }
551 :
552 : /*
553 : Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
554 : log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
555 : special values directly, passing positive non-special values through to
556 : the system log/log10.
557 : */
558 :
559 : static double
560 0 : m_log(double x)
561 : {
562 0 : if (Py_IS_FINITE(x)) {
563 0 : if (x > 0.0)
564 0 : return log(x);
565 0 : errno = EDOM;
566 0 : if (x == 0.0)
567 0 : return -Py_HUGE_VAL; /* log(0) = -inf */
568 : else
569 0 : return Py_NAN; /* log(-ve) = nan */
570 : }
571 0 : else if (Py_IS_NAN(x))
572 0 : return x; /* log(nan) = nan */
573 0 : else if (x > 0.0)
574 0 : return x; /* log(inf) = inf */
575 : else {
576 0 : errno = EDOM;
577 0 : return Py_NAN; /* log(-inf) = nan */
578 : }
579 : }
580 :
581 : /*
582 : log2: log to base 2.
583 :
584 : Uses an algorithm that should:
585 :
586 : (a) produce exact results for powers of 2, and
587 : (b) give a monotonic log2 (for positive finite floats),
588 : assuming that the system log is monotonic.
589 : */
590 :
591 : static double
592 0 : m_log2(double x)
593 : {
594 0 : if (!Py_IS_FINITE(x)) {
595 0 : if (Py_IS_NAN(x))
596 0 : return x; /* log2(nan) = nan */
597 0 : else if (x > 0.0)
598 0 : return x; /* log2(+inf) = +inf */
599 : else {
600 0 : errno = EDOM;
601 0 : return Py_NAN; /* log2(-inf) = nan, invalid-operation */
602 : }
603 : }
604 :
605 0 : if (x > 0.0) {
606 : #ifdef HAVE_LOG2
607 0 : return log2(x);
608 : #else
609 : double m;
610 : int e;
611 : m = frexp(x, &e);
612 : /* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
613 : * x is just greater than 1.0: in that case e is 1, log(m) is negative,
614 : * and we get significant cancellation error from the addition of
615 : * log(m) / log(2) to e. The slight rewrite of the expression below
616 : * avoids this problem.
617 : */
618 : if (x >= 1.0) {
619 : return log(2.0 * m) / log(2.0) + (e - 1);
620 : }
621 : else {
622 : return log(m) / log(2.0) + e;
623 : }
624 : #endif
625 : }
626 0 : else if (x == 0.0) {
627 0 : errno = EDOM;
628 0 : return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */
629 : }
630 : else {
631 0 : errno = EDOM;
632 0 : return Py_NAN; /* log2(-inf) = nan, invalid-operation */
633 : }
634 : }
635 :
636 : static double
637 0 : m_log10(double x)
638 : {
639 0 : if (Py_IS_FINITE(x)) {
640 0 : if (x > 0.0)
641 0 : return log10(x);
642 0 : errno = EDOM;
643 0 : if (x == 0.0)
644 0 : return -Py_HUGE_VAL; /* log10(0) = -inf */
645 : else
646 0 : return Py_NAN; /* log10(-ve) = nan */
647 : }
648 0 : else if (Py_IS_NAN(x))
649 0 : return x; /* log10(nan) = nan */
650 0 : else if (x > 0.0)
651 0 : return x; /* log10(inf) = inf */
652 : else {
653 0 : errno = EDOM;
654 0 : return Py_NAN; /* log10(-inf) = nan */
655 : }
656 : }
657 :
658 :
659 : /* Call is_error when errno != 0, and where x is the result libm
660 : * returned. is_error will usually set up an exception and return
661 : * true (1), but may return false (0) without setting up an exception.
662 : */
663 : static int
664 0 : is_error(double x)
665 : {
666 0 : int result = 1; /* presumption of guilt */
667 : assert(errno); /* non-zero errno is a precondition for calling */
668 0 : if (errno == EDOM)
669 0 : PyErr_SetString(PyExc_ValueError, "math domain error");
670 :
671 0 : else if (errno == ERANGE) {
672 : /* ANSI C generally requires libm functions to set ERANGE
673 : * on overflow, but also generally *allows* them to set
674 : * ERANGE on underflow too. There's no consistency about
675 : * the latter across platforms.
676 : * Alas, C99 never requires that errno be set.
677 : * Here we suppress the underflow errors (libm functions
678 : * should return a zero on underflow, and +- HUGE_VAL on
679 : * overflow, so testing the result for zero suffices to
680 : * distinguish the cases).
681 : *
682 : * On some platforms (Ubuntu/ia64) it seems that errno can be
683 : * set to ERANGE for subnormal results that do *not* underflow
684 : * to zero. So to be safe, we'll ignore ERANGE whenever the
685 : * function result is less than one in absolute value.
686 : */
687 0 : if (fabs(x) < 1.0)
688 0 : result = 0;
689 : else
690 0 : PyErr_SetString(PyExc_OverflowError,
691 : "math range error");
692 : }
693 : else
694 : /* Unexpected math error */
695 0 : PyErr_SetFromErrno(PyExc_ValueError);
696 0 : return result;
697 : }
698 :
699 : /*
700 : math_1 is used to wrap a libm function f that takes a double
701 : arguments and returns a double.
702 :
703 : The error reporting follows these rules, which are designed to do
704 : the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
705 : platforms.
706 :
707 : - a NaN result from non-NaN inputs causes ValueError to be raised
708 : - an infinite result from finite inputs causes OverflowError to be
709 : raised if can_overflow is 1, or raises ValueError if can_overflow
710 : is 0.
711 : - if the result is finite and errno == EDOM then ValueError is
712 : raised
713 : - if the result is finite and nonzero and errno == ERANGE then
714 : OverflowError is raised
715 :
716 : The last rule is used to catch overflow on platforms which follow
717 : C89 but for which HUGE_VAL is not an infinity.
718 :
719 : For the majority of one-argument functions these rules are enough
720 : to ensure that Python's functions behave as specified in 'Annex F'
721 : of the C99 standard, with the 'invalid' and 'divide-by-zero'
722 : floating-point exceptions mapping to Python's ValueError and the
723 : 'overflow' floating-point exception mapping to OverflowError.
724 : math_1 only works for functions that don't have singularities *and*
725 : the possibility of overflow; fortunately, that covers everything we
726 : care about right now.
727 : */
728 :
729 : static PyObject *
730 0 : math_1_to_whatever(PyObject *arg, double (*func) (double),
731 : PyObject *(*from_double_func) (double),
732 : int can_overflow)
733 : {
734 : double x, r;
735 0 : x = PyFloat_AsDouble(arg);
736 0 : if (x == -1.0 && PyErr_Occurred())
737 0 : return NULL;
738 0 : errno = 0;
739 : PyFPE_START_PROTECT("in math_1", return 0);
740 0 : r = (*func)(x);
741 : PyFPE_END_PROTECT(r);
742 0 : if (Py_IS_NAN(r) && !Py_IS_NAN(x)) {
743 0 : PyErr_SetString(PyExc_ValueError,
744 : "math domain error"); /* invalid arg */
745 0 : return NULL;
746 : }
747 0 : if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {
748 0 : if (can_overflow)
749 0 : PyErr_SetString(PyExc_OverflowError,
750 : "math range error"); /* overflow */
751 : else
752 0 : PyErr_SetString(PyExc_ValueError,
753 : "math domain error"); /* singularity */
754 0 : return NULL;
755 : }
756 0 : if (Py_IS_FINITE(r) && errno && is_error(r))
757 : /* this branch unnecessary on most platforms */
758 0 : return NULL;
759 :
760 0 : return (*from_double_func)(r);
761 : }
762 :
763 : /* variant of math_1, to be used when the function being wrapped is known to
764 : set errno properly (that is, errno = EDOM for invalid or divide-by-zero,
765 : errno = ERANGE for overflow). */
766 :
767 : static PyObject *
768 0 : math_1a(PyObject *arg, double (*func) (double))
769 : {
770 : double x, r;
771 0 : x = PyFloat_AsDouble(arg);
772 0 : if (x == -1.0 && PyErr_Occurred())
773 0 : return NULL;
774 0 : errno = 0;
775 : PyFPE_START_PROTECT("in math_1a", return 0);
776 0 : r = (*func)(x);
777 : PyFPE_END_PROTECT(r);
778 0 : if (errno && is_error(r))
779 0 : return NULL;
780 0 : return PyFloat_FromDouble(r);
781 : }
782 :
783 : /*
784 : math_2 is used to wrap a libm function f that takes two double
785 : arguments and returns a double.
786 :
787 : The error reporting follows these rules, which are designed to do
788 : the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
789 : platforms.
790 :
791 : - a NaN result from non-NaN inputs causes ValueError to be raised
792 : - an infinite result from finite inputs causes OverflowError to be
793 : raised.
794 : - if the result is finite and errno == EDOM then ValueError is
795 : raised
796 : - if the result is finite and nonzero and errno == ERANGE then
797 : OverflowError is raised
798 :
799 : The last rule is used to catch overflow on platforms which follow
800 : C89 but for which HUGE_VAL is not an infinity.
801 :
802 : For most two-argument functions (copysign, fmod, hypot, atan2)
803 : these rules are enough to ensure that Python's functions behave as
804 : specified in 'Annex F' of the C99 standard, with the 'invalid' and
805 : 'divide-by-zero' floating-point exceptions mapping to Python's
806 : ValueError and the 'overflow' floating-point exception mapping to
807 : OverflowError.
808 : */
809 :
810 : static PyObject *
811 0 : math_1(PyObject *arg, double (*func) (double), int can_overflow)
812 : {
813 0 : return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
814 : }
815 :
816 : static PyObject *
817 0 : math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
818 : {
819 0 : return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
820 : }
821 :
822 : static PyObject *
823 0 : math_2(PyObject *args, double (*func) (double, double), char *funcname)
824 : {
825 : PyObject *ox, *oy;
826 : double x, y, r;
827 0 : if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
828 0 : return NULL;
829 0 : x = PyFloat_AsDouble(ox);
830 0 : y = PyFloat_AsDouble(oy);
831 0 : if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
832 0 : return NULL;
833 0 : errno = 0;
834 : PyFPE_START_PROTECT("in math_2", return 0);
835 0 : r = (*func)(x, y);
836 : PyFPE_END_PROTECT(r);
837 0 : if (Py_IS_NAN(r)) {
838 0 : if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
839 0 : errno = EDOM;
840 : else
841 0 : errno = 0;
842 : }
843 0 : else if (Py_IS_INFINITY(r)) {
844 0 : if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
845 0 : errno = ERANGE;
846 : else
847 0 : errno = 0;
848 : }
849 0 : if (errno && is_error(r))
850 0 : return NULL;
851 : else
852 0 : return PyFloat_FromDouble(r);
853 : }
854 :
855 : #define FUNC1(funcname, func, can_overflow, docstring) \
856 : static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
857 : return math_1(args, func, can_overflow); \
858 : }\
859 : PyDoc_STRVAR(math_##funcname##_doc, docstring);
860 :
861 : #define FUNC1A(funcname, func, docstring) \
862 : static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
863 : return math_1a(args, func); \
864 : }\
865 : PyDoc_STRVAR(math_##funcname##_doc, docstring);
866 :
867 : #define FUNC2(funcname, func, docstring) \
868 : static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
869 : return math_2(args, func, #funcname); \
870 : }\
871 : PyDoc_STRVAR(math_##funcname##_doc, docstring);
872 :
873 0 : FUNC1(acos, acos, 0,
874 0 : "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
875 0 : FUNC1(acosh, m_acosh, 0,
876 0 : "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
877 0 : FUNC1(asin, asin, 0,
878 0 : "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
879 0 : FUNC1(asinh, m_asinh, 0,
880 0 : "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
881 0 : FUNC1(atan, atan, 0,
882 0 : "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
883 0 : FUNC2(atan2, m_atan2,
884 : "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
885 0 : "Unlike atan(y/x), the signs of both x and y are considered.")
886 0 : FUNC1(atanh, m_atanh, 0,
887 0 : "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
888 :
889 0 : static PyObject * math_ceil(PyObject *self, PyObject *number) {
890 : _Py_IDENTIFIER(__ceil__);
891 : PyObject *method, *result;
892 :
893 0 : method = _PyObject_LookupSpecial(number, &PyId___ceil__);
894 0 : if (method == NULL) {
895 0 : if (PyErr_Occurred())
896 0 : return NULL;
897 0 : return math_1_to_int(number, ceil, 0);
898 : }
899 0 : result = PyObject_CallFunctionObjArgs(method, NULL);
900 0 : Py_DECREF(method);
901 0 : return result;
902 : }
903 :
904 : PyDoc_STRVAR(math_ceil_doc,
905 : "ceil(x)\n\nReturn the ceiling of x as an int.\n"
906 : "This is the smallest integral value >= x.");
907 :
908 0 : FUNC2(copysign, copysign,
909 0 : "copysign(x, y)\n\nReturn x with the sign of y.")
910 0 : FUNC1(cos, cos, 0,
911 0 : "cos(x)\n\nReturn the cosine of x (measured in radians).")
912 0 : FUNC1(cosh, cosh, 1,
913 0 : "cosh(x)\n\nReturn the hyperbolic cosine of x.")
914 0 : FUNC1A(erf, m_erf,
915 0 : "erf(x)\n\nError function at x.")
916 0 : FUNC1A(erfc, m_erfc,
917 0 : "erfc(x)\n\nComplementary error function at x.")
918 0 : FUNC1(exp, exp, 1,
919 0 : "exp(x)\n\nReturn e raised to the power of x.")
920 0 : FUNC1(expm1, m_expm1, 1,
921 : "expm1(x)\n\nReturn exp(x)-1.\n"
922 : "This function avoids the loss of precision involved in the direct "
923 0 : "evaluation of exp(x)-1 for small x.")
924 0 : FUNC1(fabs, fabs, 0,
925 0 : "fabs(x)\n\nReturn the absolute value of the float x.")
926 :
927 0 : static PyObject * math_floor(PyObject *self, PyObject *number) {
928 : _Py_IDENTIFIER(__floor__);
929 : PyObject *method, *result;
930 :
931 0 : method = _PyObject_LookupSpecial(number, &PyId___floor__);
932 0 : if (method == NULL) {
933 0 : if (PyErr_Occurred())
934 0 : return NULL;
935 0 : return math_1_to_int(number, floor, 0);
936 : }
937 0 : result = PyObject_CallFunctionObjArgs(method, NULL);
938 0 : Py_DECREF(method);
939 0 : return result;
940 : }
941 :
942 : PyDoc_STRVAR(math_floor_doc,
943 : "floor(x)\n\nReturn the floor of x as an int.\n"
944 : "This is the largest integral value <= x.");
945 :
946 0 : FUNC1A(gamma, m_tgamma,
947 0 : "gamma(x)\n\nGamma function at x.")
948 0 : FUNC1A(lgamma, m_lgamma,
949 0 : "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
950 0 : FUNC1(log1p, m_log1p, 0,
951 : "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
952 0 : "The result is computed in a way which is accurate for x near zero.")
953 0 : FUNC1(sin, sin, 0,
954 0 : "sin(x)\n\nReturn the sine of x (measured in radians).")
955 0 : FUNC1(sinh, sinh, 1,
956 0 : "sinh(x)\n\nReturn the hyperbolic sine of x.")
957 0 : FUNC1(sqrt, sqrt, 0,
958 0 : "sqrt(x)\n\nReturn the square root of x.")
959 0 : FUNC1(tan, tan, 0,
960 0 : "tan(x)\n\nReturn the tangent of x (measured in radians).")
961 0 : FUNC1(tanh, tanh, 0,
962 0 : "tanh(x)\n\nReturn the hyperbolic tangent of x.")
963 :
964 : /* Precision summation function as msum() by Raymond Hettinger in
965 : <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
966 : enhanced with the exact partials sum and roundoff from Mark
967 : Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
968 : See those links for more details, proofs and other references.
969 :
970 : Note 1: IEEE 754R floating point semantics are assumed,
971 : but the current implementation does not re-establish special
972 : value semantics across iterations (i.e. handling -Inf + Inf).
973 :
974 : Note 2: No provision is made for intermediate overflow handling;
975 : therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
976 : sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
977 : overflow of the first partial sum.
978 :
979 : Note 3: The intermediate values lo, yr, and hi are declared volatile so
980 : aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
981 : Also, the volatile declaration forces the values to be stored in memory as
982 : regular doubles instead of extended long precision (80-bit) values. This
983 : prevents double rounding because any addition or subtraction of two doubles
984 : can be resolved exactly into double-sized hi and lo values. As long as the
985 : hi value gets forced into a double before yr and lo are computed, the extra
986 : bits in downstream extended precision operations (x87 for example) will be
987 : exactly zero and therefore can be losslessly stored back into a double,
988 : thereby preventing double rounding.
989 :
990 : Note 4: A similar implementation is in Modules/cmathmodule.c.
991 : Be sure to update both when making changes.
992 :
993 : Note 5: The signature of math.fsum() differs from __builtin__.sum()
994 : because the start argument doesn't make sense in the context of
995 : accurate summation. Since the partials table is collapsed before
996 : returning a result, sum(seq2, start=sum(seq1)) may not equal the
997 : accurate result returned by sum(itertools.chain(seq1, seq2)).
998 : */
999 :
1000 : #define NUM_PARTIALS 32 /* initial partials array size, on stack */
1001 :
1002 : /* Extend the partials array p[] by doubling its size. */
1003 : static int /* non-zero on error */
1004 0 : _fsum_realloc(double **p_ptr, Py_ssize_t n,
1005 : double *ps, Py_ssize_t *m_ptr)
1006 : {
1007 0 : void *v = NULL;
1008 0 : Py_ssize_t m = *m_ptr;
1009 :
1010 0 : m += m; /* double */
1011 0 : if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
1012 0 : double *p = *p_ptr;
1013 0 : if (p == ps) {
1014 0 : v = PyMem_Malloc(sizeof(double) * m);
1015 0 : if (v != NULL)
1016 0 : memcpy(v, ps, sizeof(double) * n);
1017 : }
1018 : else
1019 0 : v = PyMem_Realloc(p, sizeof(double) * m);
1020 : }
1021 0 : if (v == NULL) { /* size overflow or no memory */
1022 0 : PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
1023 0 : return 1;
1024 : }
1025 0 : *p_ptr = (double*) v;
1026 0 : *m_ptr = m;
1027 0 : return 0;
1028 : }
1029 :
1030 : /* Full precision summation of a sequence of floats.
1031 :
1032 : def msum(iterable):
1033 : partials = [] # sorted, non-overlapping partial sums
1034 : for x in iterable:
1035 : i = 0
1036 : for y in partials:
1037 : if abs(x) < abs(y):
1038 : x, y = y, x
1039 : hi = x + y
1040 : lo = y - (hi - x)
1041 : if lo:
1042 : partials[i] = lo
1043 : i += 1
1044 : x = hi
1045 : partials[i:] = [x]
1046 : return sum_exact(partials)
1047 :
1048 : Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo
1049 : are exactly equal to x+y. The inner loop applies hi/lo summation to each
1050 : partial so that the list of partial sums remains exact.
1051 :
1052 : Sum_exact() adds the partial sums exactly and correctly rounds the final
1053 : result (using the round-half-to-even rule). The items in partials remain
1054 : non-zero, non-special, non-overlapping and strictly increasing in
1055 : magnitude, but possibly not all having the same sign.
1056 :
1057 : Depends on IEEE 754 arithmetic guarantees and half-even rounding.
1058 : */
1059 :
1060 : static PyObject*
1061 0 : math_fsum(PyObject *self, PyObject *seq)
1062 : {
1063 0 : PyObject *item, *iter, *sum = NULL;
1064 0 : Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
1065 0 : double x, y, t, ps[NUM_PARTIALS], *p = ps;
1066 0 : double xsave, special_sum = 0.0, inf_sum = 0.0;
1067 : volatile double hi, yr, lo;
1068 :
1069 0 : iter = PyObject_GetIter(seq);
1070 0 : if (iter == NULL)
1071 0 : return NULL;
1072 :
1073 : PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
1074 :
1075 : for(;;) { /* for x in iterable */
1076 : assert(0 <= n && n <= m);
1077 : assert((m == NUM_PARTIALS && p == ps) ||
1078 : (m > NUM_PARTIALS && p != NULL));
1079 :
1080 0 : item = PyIter_Next(iter);
1081 0 : if (item == NULL) {
1082 0 : if (PyErr_Occurred())
1083 0 : goto _fsum_error;
1084 0 : break;
1085 : }
1086 0 : x = PyFloat_AsDouble(item);
1087 0 : Py_DECREF(item);
1088 0 : if (PyErr_Occurred())
1089 0 : goto _fsum_error;
1090 :
1091 0 : xsave = x;
1092 0 : for (i = j = 0; j < n; j++) { /* for y in partials */
1093 0 : y = p[j];
1094 0 : if (fabs(x) < fabs(y)) {
1095 0 : t = x; x = y; y = t;
1096 : }
1097 0 : hi = x + y;
1098 0 : yr = hi - x;
1099 0 : lo = y - yr;
1100 0 : if (lo != 0.0)
1101 0 : p[i++] = lo;
1102 0 : x = hi;
1103 : }
1104 :
1105 0 : n = i; /* ps[i:] = [x] */
1106 0 : if (x != 0.0) {
1107 0 : if (! Py_IS_FINITE(x)) {
1108 : /* a nonfinite x could arise either as
1109 : a result of intermediate overflow, or
1110 : as a result of a nan or inf in the
1111 : summands */
1112 0 : if (Py_IS_FINITE(xsave)) {
1113 0 : PyErr_SetString(PyExc_OverflowError,
1114 : "intermediate overflow in fsum");
1115 0 : goto _fsum_error;
1116 : }
1117 0 : if (Py_IS_INFINITY(xsave))
1118 0 : inf_sum += xsave;
1119 0 : special_sum += xsave;
1120 : /* reset partials */
1121 0 : n = 0;
1122 : }
1123 0 : else if (n >= m && _fsum_realloc(&p, n, ps, &m))
1124 : goto _fsum_error;
1125 : else
1126 0 : p[n++] = x;
1127 : }
1128 0 : }
1129 :
1130 0 : if (special_sum != 0.0) {
1131 0 : if (Py_IS_NAN(inf_sum))
1132 0 : PyErr_SetString(PyExc_ValueError,
1133 : "-inf + inf in fsum");
1134 : else
1135 0 : sum = PyFloat_FromDouble(special_sum);
1136 0 : goto _fsum_error;
1137 : }
1138 :
1139 0 : hi = 0.0;
1140 0 : if (n > 0) {
1141 0 : hi = p[--n];
1142 : /* sum_exact(ps, hi) from the top, stop when the sum becomes
1143 : inexact. */
1144 0 : while (n > 0) {
1145 0 : x = hi;
1146 0 : y = p[--n];
1147 : assert(fabs(y) < fabs(x));
1148 0 : hi = x + y;
1149 0 : yr = hi - x;
1150 0 : lo = y - yr;
1151 0 : if (lo != 0.0)
1152 0 : break;
1153 : }
1154 : /* Make half-even rounding work across multiple partials.
1155 : Needed so that sum([1e-16, 1, 1e16]) will round-up the last
1156 : digit to two instead of down to zero (the 1e-16 makes the 1
1157 : slightly closer to two). With a potential 1 ULP rounding
1158 : error fixed-up, math.fsum() can guarantee commutativity. */
1159 0 : if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
1160 0 : (lo > 0.0 && p[n-1] > 0.0))) {
1161 0 : y = lo * 2.0;
1162 0 : x = hi + y;
1163 0 : yr = x - hi;
1164 0 : if (y == yr)
1165 0 : hi = x;
1166 : }
1167 : }
1168 0 : sum = PyFloat_FromDouble(hi);
1169 :
1170 : _fsum_error:
1171 : PyFPE_END_PROTECT(hi)
1172 0 : Py_DECREF(iter);
1173 0 : if (p != ps)
1174 0 : PyMem_Free(p);
1175 0 : return sum;
1176 : }
1177 :
1178 : #undef NUM_PARTIALS
1179 :
1180 : PyDoc_STRVAR(math_fsum_doc,
1181 : "fsum(iterable)\n\n\
1182 : Return an accurate floating point sum of values in the iterable.\n\
1183 : Assumes IEEE-754 floating point arithmetic.");
1184 :
1185 : /* Return the smallest integer k such that n < 2**k, or 0 if n == 0.
1186 : * Equivalent to floor(lg(x))+1. Also equivalent to: bitwidth_of_type -
1187 : * count_leading_zero_bits(x)
1188 : */
1189 :
1190 : /* XXX: This routine does more or less the same thing as
1191 : * bits_in_digit() in Objects/longobject.c. Someday it would be nice to
1192 : * consolidate them. On BSD, there's a library function called fls()
1193 : * that we could use, and GCC provides __builtin_clz().
1194 : */
1195 :
1196 : static unsigned long
1197 0 : bit_length(unsigned long n)
1198 : {
1199 0 : unsigned long len = 0;
1200 0 : while (n != 0) {
1201 0 : ++len;
1202 0 : n >>= 1;
1203 : }
1204 0 : return len;
1205 : }
1206 :
1207 : static unsigned long
1208 0 : count_set_bits(unsigned long n)
1209 : {
1210 0 : unsigned long count = 0;
1211 0 : while (n != 0) {
1212 0 : ++count;
1213 0 : n &= n - 1; /* clear least significant bit */
1214 : }
1215 0 : return count;
1216 : }
1217 :
1218 : /* Divide-and-conquer factorial algorithm
1219 : *
1220 : * Based on the formula and psuedo-code provided at:
1221 : * http://www.luschny.de/math/factorial/binarysplitfact.html
1222 : *
1223 : * Faster algorithms exist, but they're more complicated and depend on
1224 : * a fast prime factorization algorithm.
1225 : *
1226 : * Notes on the algorithm
1227 : * ----------------------
1228 : *
1229 : * factorial(n) is written in the form 2**k * m, with m odd. k and m are
1230 : * computed separately, and then combined using a left shift.
1231 : *
1232 : * The function factorial_odd_part computes the odd part m (i.e., the greatest
1233 : * odd divisor) of factorial(n), using the formula:
1234 : *
1235 : * factorial_odd_part(n) =
1236 : *
1237 : * product_{i >= 0} product_{0 < j <= n / 2**i, j odd} j
1238 : *
1239 : * Example: factorial_odd_part(20) =
1240 : *
1241 : * (1) *
1242 : * (1) *
1243 : * (1 * 3 * 5) *
1244 : * (1 * 3 * 5 * 7 * 9)
1245 : * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1246 : *
1247 : * Here i goes from large to small: the first term corresponds to i=4 (any
1248 : * larger i gives an empty product), and the last term corresponds to i=0.
1249 : * Each term can be computed from the last by multiplying by the extra odd
1250 : * numbers required: e.g., to get from the penultimate term to the last one,
1251 : * we multiply by (11 * 13 * 15 * 17 * 19).
1252 : *
1253 : * To see a hint of why this formula works, here are the same numbers as above
1254 : * but with the even parts (i.e., the appropriate powers of 2) included. For
1255 : * each subterm in the product for i, we multiply that subterm by 2**i:
1256 : *
1257 : * factorial(20) =
1258 : *
1259 : * (16) *
1260 : * (8) *
1261 : * (4 * 12 * 20) *
1262 : * (2 * 6 * 10 * 14 * 18) *
1263 : * (1 * 3 * 5 * 7 * 9 * 11 * 13 * 15 * 17 * 19)
1264 : *
1265 : * The factorial_partial_product function computes the product of all odd j in
1266 : * range(start, stop) for given start and stop. It's used to compute the
1267 : * partial products like (11 * 13 * 15 * 17 * 19) in the example above. It
1268 : * operates recursively, repeatedly splitting the range into two roughly equal
1269 : * pieces until the subranges are small enough to be computed using only C
1270 : * integer arithmetic.
1271 : *
1272 : * The two-valuation k (i.e., the exponent of the largest power of 2 dividing
1273 : * the factorial) is computed independently in the main math_factorial
1274 : * function. By standard results, its value is:
1275 : *
1276 : * two_valuation = n//2 + n//4 + n//8 + ....
1277 : *
1278 : * It can be shown (e.g., by complete induction on n) that two_valuation is
1279 : * equal to n - count_set_bits(n), where count_set_bits(n) gives the number of
1280 : * '1'-bits in the binary expansion of n.
1281 : */
1282 :
1283 : /* factorial_partial_product: Compute product(range(start, stop, 2)) using
1284 : * divide and conquer. Assumes start and stop are odd and stop > start.
1285 : * max_bits must be >= bit_length(stop - 2). */
1286 :
1287 : static PyObject *
1288 0 : factorial_partial_product(unsigned long start, unsigned long stop,
1289 : unsigned long max_bits)
1290 : {
1291 : unsigned long midpoint, num_operands;
1292 0 : PyObject *left = NULL, *right = NULL, *result = NULL;
1293 :
1294 : /* If the return value will fit an unsigned long, then we can
1295 : * multiply in a tight, fast loop where each multiply is O(1).
1296 : * Compute an upper bound on the number of bits required to store
1297 : * the answer.
1298 : *
1299 : * Storing some integer z requires floor(lg(z))+1 bits, which is
1300 : * conveniently the value returned by bit_length(z). The
1301 : * product x*y will require at most
1302 : * bit_length(x) + bit_length(y) bits to store, based
1303 : * on the idea that lg product = lg x + lg y.
1304 : *
1305 : * We know that stop - 2 is the largest number to be multiplied. From
1306 : * there, we have: bit_length(answer) <= num_operands *
1307 : * bit_length(stop - 2)
1308 : */
1309 :
1310 0 : num_operands = (stop - start) / 2;
1311 : /* The "num_operands <= 8 * SIZEOF_LONG" check guards against the
1312 : * unlikely case of an overflow in num_operands * max_bits. */
1313 0 : if (num_operands <= 8 * SIZEOF_LONG &&
1314 0 : num_operands * max_bits <= 8 * SIZEOF_LONG) {
1315 : unsigned long j, total;
1316 0 : for (total = start, j = start + 2; j < stop; j += 2)
1317 0 : total *= j;
1318 0 : return PyLong_FromUnsignedLong(total);
1319 : }
1320 :
1321 : /* find midpoint of range(start, stop), rounded up to next odd number. */
1322 0 : midpoint = (start + num_operands) | 1;
1323 0 : left = factorial_partial_product(start, midpoint,
1324 : bit_length(midpoint - 2));
1325 0 : if (left == NULL)
1326 0 : goto error;
1327 0 : right = factorial_partial_product(midpoint, stop, max_bits);
1328 0 : if (right == NULL)
1329 0 : goto error;
1330 0 : result = PyNumber_Multiply(left, right);
1331 :
1332 : error:
1333 0 : Py_XDECREF(left);
1334 0 : Py_XDECREF(right);
1335 0 : return result;
1336 : }
1337 :
1338 : /* factorial_odd_part: compute the odd part of factorial(n). */
1339 :
1340 : static PyObject *
1341 0 : factorial_odd_part(unsigned long n)
1342 : {
1343 : long i;
1344 : unsigned long v, lower, upper;
1345 : PyObject *partial, *tmp, *inner, *outer;
1346 :
1347 0 : inner = PyLong_FromLong(1);
1348 0 : if (inner == NULL)
1349 0 : return NULL;
1350 0 : outer = inner;
1351 0 : Py_INCREF(outer);
1352 :
1353 0 : upper = 3;
1354 0 : for (i = bit_length(n) - 2; i >= 0; i--) {
1355 0 : v = n >> i;
1356 0 : if (v <= 2)
1357 0 : continue;
1358 0 : lower = upper;
1359 : /* (v + 1) | 1 = least odd integer strictly larger than n / 2**i */
1360 0 : upper = (v + 1) | 1;
1361 : /* Here inner is the product of all odd integers j in the range (0,
1362 : n/2**(i+1)]. The factorial_partial_product call below gives the
1363 : product of all odd integers j in the range (n/2**(i+1), n/2**i]. */
1364 0 : partial = factorial_partial_product(lower, upper, bit_length(upper-2));
1365 : /* inner *= partial */
1366 0 : if (partial == NULL)
1367 0 : goto error;
1368 0 : tmp = PyNumber_Multiply(inner, partial);
1369 0 : Py_DECREF(partial);
1370 0 : if (tmp == NULL)
1371 0 : goto error;
1372 0 : Py_DECREF(inner);
1373 0 : inner = tmp;
1374 : /* Now inner is the product of all odd integers j in the range (0,
1375 : n/2**i], giving the inner product in the formula above. */
1376 :
1377 : /* outer *= inner; */
1378 0 : tmp = PyNumber_Multiply(outer, inner);
1379 0 : if (tmp == NULL)
1380 0 : goto error;
1381 0 : Py_DECREF(outer);
1382 0 : outer = tmp;
1383 : }
1384 :
1385 0 : goto done;
1386 :
1387 : error:
1388 0 : Py_DECREF(outer);
1389 : done:
1390 0 : Py_DECREF(inner);
1391 0 : return outer;
1392 : }
1393 :
1394 : /* Lookup table for small factorial values */
1395 :
1396 : static const unsigned long SmallFactorials[] = {
1397 : 1, 1, 2, 6, 24, 120, 720, 5040, 40320,
1398 : 362880, 3628800, 39916800, 479001600,
1399 : #if SIZEOF_LONG >= 8
1400 : 6227020800, 87178291200, 1307674368000,
1401 : 20922789888000, 355687428096000, 6402373705728000,
1402 : 121645100408832000, 2432902008176640000
1403 : #endif
1404 : };
1405 :
1406 : static PyObject *
1407 0 : math_factorial(PyObject *self, PyObject *arg)
1408 : {
1409 : long x;
1410 : PyObject *result, *odd_part, *two_valuation;
1411 :
1412 0 : if (PyFloat_Check(arg)) {
1413 : PyObject *lx;
1414 0 : double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
1415 0 : if (!(Py_IS_FINITE(dx) && dx == floor(dx))) {
1416 0 : PyErr_SetString(PyExc_ValueError,
1417 : "factorial() only accepts integral values");
1418 0 : return NULL;
1419 : }
1420 0 : lx = PyLong_FromDouble(dx);
1421 0 : if (lx == NULL)
1422 0 : return NULL;
1423 0 : x = PyLong_AsLong(lx);
1424 0 : Py_DECREF(lx);
1425 : }
1426 : else
1427 0 : x = PyLong_AsLong(arg);
1428 :
1429 0 : if (x == -1 && PyErr_Occurred())
1430 0 : return NULL;
1431 0 : if (x < 0) {
1432 0 : PyErr_SetString(PyExc_ValueError,
1433 : "factorial() not defined for negative values");
1434 0 : return NULL;
1435 : }
1436 :
1437 : /* use lookup table if x is small */
1438 0 : if (x < (long)Py_ARRAY_LENGTH(SmallFactorials))
1439 0 : return PyLong_FromUnsignedLong(SmallFactorials[x]);
1440 :
1441 : /* else express in the form odd_part * 2**two_valuation, and compute as
1442 : odd_part << two_valuation. */
1443 0 : odd_part = factorial_odd_part(x);
1444 0 : if (odd_part == NULL)
1445 0 : return NULL;
1446 0 : two_valuation = PyLong_FromLong(x - count_set_bits(x));
1447 0 : if (two_valuation == NULL) {
1448 0 : Py_DECREF(odd_part);
1449 0 : return NULL;
1450 : }
1451 0 : result = PyNumber_Lshift(odd_part, two_valuation);
1452 0 : Py_DECREF(two_valuation);
1453 0 : Py_DECREF(odd_part);
1454 0 : return result;
1455 : }
1456 :
1457 : PyDoc_STRVAR(math_factorial_doc,
1458 : "factorial(x) -> Integral\n"
1459 : "\n"
1460 : "Find x!. Raise a ValueError if x is negative or non-integral.");
1461 :
1462 : static PyObject *
1463 0 : math_trunc(PyObject *self, PyObject *number)
1464 : {
1465 : _Py_IDENTIFIER(__trunc__);
1466 : PyObject *trunc, *result;
1467 :
1468 0 : if (Py_TYPE(number)->tp_dict == NULL) {
1469 0 : if (PyType_Ready(Py_TYPE(number)) < 0)
1470 0 : return NULL;
1471 : }
1472 :
1473 0 : trunc = _PyObject_LookupSpecial(number, &PyId___trunc__);
1474 0 : if (trunc == NULL) {
1475 0 : if (!PyErr_Occurred())
1476 0 : PyErr_Format(PyExc_TypeError,
1477 : "type %.100s doesn't define __trunc__ method",
1478 0 : Py_TYPE(number)->tp_name);
1479 0 : return NULL;
1480 : }
1481 0 : result = PyObject_CallFunctionObjArgs(trunc, NULL);
1482 0 : Py_DECREF(trunc);
1483 0 : return result;
1484 : }
1485 :
1486 : PyDoc_STRVAR(math_trunc_doc,
1487 : "trunc(x:Real) -> Integral\n"
1488 : "\n"
1489 : "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
1490 :
1491 : static PyObject *
1492 0 : math_frexp(PyObject *self, PyObject *arg)
1493 : {
1494 : int i;
1495 0 : double x = PyFloat_AsDouble(arg);
1496 0 : if (x == -1.0 && PyErr_Occurred())
1497 0 : return NULL;
1498 : /* deal with special cases directly, to sidestep platform
1499 : differences */
1500 0 : if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
1501 0 : i = 0;
1502 : }
1503 : else {
1504 : PyFPE_START_PROTECT("in math_frexp", return 0);
1505 0 : x = frexp(x, &i);
1506 : PyFPE_END_PROTECT(x);
1507 : }
1508 0 : return Py_BuildValue("(di)", x, i);
1509 : }
1510 :
1511 : PyDoc_STRVAR(math_frexp_doc,
1512 : "frexp(x)\n"
1513 : "\n"
1514 : "Return the mantissa and exponent of x, as pair (m, e).\n"
1515 : "m is a float and e is an int, such that x = m * 2.**e.\n"
1516 : "If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0.");
1517 :
1518 : static PyObject *
1519 0 : math_ldexp(PyObject *self, PyObject *args)
1520 : {
1521 : double x, r;
1522 : PyObject *oexp;
1523 : long exp;
1524 : int overflow;
1525 0 : if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
1526 0 : return NULL;
1527 :
1528 0 : if (PyLong_Check(oexp)) {
1529 : /* on overflow, replace exponent with either LONG_MAX
1530 : or LONG_MIN, depending on the sign. */
1531 0 : exp = PyLong_AsLongAndOverflow(oexp, &overflow);
1532 0 : if (exp == -1 && PyErr_Occurred())
1533 0 : return NULL;
1534 0 : if (overflow)
1535 0 : exp = overflow < 0 ? LONG_MIN : LONG_MAX;
1536 : }
1537 : else {
1538 0 : PyErr_SetString(PyExc_TypeError,
1539 : "Expected an int or long as second argument "
1540 : "to ldexp.");
1541 0 : return NULL;
1542 : }
1543 :
1544 0 : if (x == 0. || !Py_IS_FINITE(x)) {
1545 : /* NaNs, zeros and infinities are returned unchanged */
1546 0 : r = x;
1547 0 : errno = 0;
1548 : } else if (exp > INT_MAX) {
1549 : /* overflow */
1550 : r = copysign(Py_HUGE_VAL, x);
1551 : errno = ERANGE;
1552 : } else if (exp < INT_MIN) {
1553 : /* underflow to +-0 */
1554 : r = copysign(0., x);
1555 : errno = 0;
1556 : } else {
1557 0 : errno = 0;
1558 : PyFPE_START_PROTECT("in math_ldexp", return 0);
1559 0 : r = ldexp(x, (int)exp);
1560 : PyFPE_END_PROTECT(r);
1561 0 : if (Py_IS_INFINITY(r))
1562 0 : errno = ERANGE;
1563 : }
1564 :
1565 0 : if (errno && is_error(r))
1566 0 : return NULL;
1567 0 : return PyFloat_FromDouble(r);
1568 : }
1569 :
1570 : PyDoc_STRVAR(math_ldexp_doc,
1571 : "ldexp(x, i)\n\n\
1572 : Return x * (2**i).");
1573 :
1574 : static PyObject *
1575 0 : math_modf(PyObject *self, PyObject *arg)
1576 : {
1577 0 : double y, x = PyFloat_AsDouble(arg);
1578 0 : if (x == -1.0 && PyErr_Occurred())
1579 0 : return NULL;
1580 : /* some platforms don't do the right thing for NaNs and
1581 : infinities, so we take care of special cases directly. */
1582 0 : if (!Py_IS_FINITE(x)) {
1583 0 : if (Py_IS_INFINITY(x))
1584 0 : return Py_BuildValue("(dd)", copysign(0., x), x);
1585 0 : else if (Py_IS_NAN(x))
1586 0 : return Py_BuildValue("(dd)", x, x);
1587 : }
1588 :
1589 0 : errno = 0;
1590 : PyFPE_START_PROTECT("in math_modf", return 0);
1591 0 : x = modf(x, &y);
1592 : PyFPE_END_PROTECT(x);
1593 0 : return Py_BuildValue("(dd)", x, y);
1594 : }
1595 :
1596 : PyDoc_STRVAR(math_modf_doc,
1597 : "modf(x)\n"
1598 : "\n"
1599 : "Return the fractional and integer parts of x. Both results carry the sign\n"
1600 : "of x and are floats.");
1601 :
1602 : /* A decent logarithm is easy to compute even for huge longs, but libm can't
1603 : do that by itself -- loghelper can. func is log or log10, and name is
1604 : "log" or "log10". Note that overflow of the result isn't possible: a long
1605 : can contain no more than INT_MAX * SHIFT bits, so has value certainly less
1606 : than 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
1607 : small enough to fit in an IEEE single. log and log10 are even smaller.
1608 : However, intermediate overflow is possible for a long if the number of bits
1609 : in that long is larger than PY_SSIZE_T_MAX. */
1610 :
1611 : static PyObject*
1612 0 : loghelper(PyObject* arg, double (*func)(double), char *funcname)
1613 : {
1614 : /* If it is long, do it ourselves. */
1615 0 : if (PyLong_Check(arg)) {
1616 : double x, result;
1617 : Py_ssize_t e;
1618 :
1619 : /* Negative or zero inputs give a ValueError. */
1620 0 : if (Py_SIZE(arg) <= 0) {
1621 0 : PyErr_SetString(PyExc_ValueError,
1622 : "math domain error");
1623 0 : return NULL;
1624 : }
1625 :
1626 0 : x = PyLong_AsDouble(arg);
1627 0 : if (x == -1.0 && PyErr_Occurred()) {
1628 0 : if (!PyErr_ExceptionMatches(PyExc_OverflowError))
1629 0 : return NULL;
1630 : /* Here the conversion to double overflowed, but it's possible
1631 : to compute the log anyway. Clear the exception and continue. */
1632 0 : PyErr_Clear();
1633 0 : x = _PyLong_Frexp((PyLongObject *)arg, &e);
1634 0 : if (x == -1.0 && PyErr_Occurred())
1635 0 : return NULL;
1636 : /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */
1637 0 : result = func(x) + func(2.0) * e;
1638 : }
1639 : else
1640 : /* Successfully converted x to a double. */
1641 0 : result = func(x);
1642 0 : return PyFloat_FromDouble(result);
1643 : }
1644 :
1645 : /* Else let libm handle it by itself. */
1646 0 : return math_1(arg, func, 0);
1647 : }
1648 :
1649 : static PyObject *
1650 0 : math_log(PyObject *self, PyObject *args)
1651 : {
1652 : PyObject *arg;
1653 0 : PyObject *base = NULL;
1654 : PyObject *num, *den;
1655 : PyObject *ans;
1656 :
1657 0 : if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
1658 0 : return NULL;
1659 :
1660 0 : num = loghelper(arg, m_log, "log");
1661 0 : if (num == NULL || base == NULL)
1662 0 : return num;
1663 :
1664 0 : den = loghelper(base, m_log, "log");
1665 0 : if (den == NULL) {
1666 0 : Py_DECREF(num);
1667 0 : return NULL;
1668 : }
1669 :
1670 0 : ans = PyNumber_TrueDivide(num, den);
1671 0 : Py_DECREF(num);
1672 0 : Py_DECREF(den);
1673 0 : return ans;
1674 : }
1675 :
1676 : PyDoc_STRVAR(math_log_doc,
1677 : "log(x[, base])\n\n\
1678 : Return the logarithm of x to the given base.\n\
1679 : If the base not specified, returns the natural logarithm (base e) of x.");
1680 :
1681 : static PyObject *
1682 0 : math_log2(PyObject *self, PyObject *arg)
1683 : {
1684 0 : return loghelper(arg, m_log2, "log2");
1685 : }
1686 :
1687 : PyDoc_STRVAR(math_log2_doc,
1688 : "log2(x)\n\nReturn the base 2 logarithm of x.");
1689 :
1690 : static PyObject *
1691 0 : math_log10(PyObject *self, PyObject *arg)
1692 : {
1693 0 : return loghelper(arg, m_log10, "log10");
1694 : }
1695 :
1696 : PyDoc_STRVAR(math_log10_doc,
1697 : "log10(x)\n\nReturn the base 10 logarithm of x.");
1698 :
1699 : static PyObject *
1700 0 : math_fmod(PyObject *self, PyObject *args)
1701 : {
1702 : PyObject *ox, *oy;
1703 : double r, x, y;
1704 0 : if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
1705 0 : return NULL;
1706 0 : x = PyFloat_AsDouble(ox);
1707 0 : y = PyFloat_AsDouble(oy);
1708 0 : if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1709 0 : return NULL;
1710 : /* fmod(x, +/-Inf) returns x for finite x. */
1711 0 : if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
1712 0 : return PyFloat_FromDouble(x);
1713 0 : errno = 0;
1714 : PyFPE_START_PROTECT("in math_fmod", return 0);
1715 0 : r = fmod(x, y);
1716 : PyFPE_END_PROTECT(r);
1717 0 : if (Py_IS_NAN(r)) {
1718 0 : if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1719 0 : errno = EDOM;
1720 : else
1721 0 : errno = 0;
1722 : }
1723 0 : if (errno && is_error(r))
1724 0 : return NULL;
1725 : else
1726 0 : return PyFloat_FromDouble(r);
1727 : }
1728 :
1729 : PyDoc_STRVAR(math_fmod_doc,
1730 : "fmod(x, y)\n\nReturn fmod(x, y), according to platform C."
1731 : " x % y may differ.");
1732 :
1733 : static PyObject *
1734 0 : math_hypot(PyObject *self, PyObject *args)
1735 : {
1736 : PyObject *ox, *oy;
1737 : double r, x, y;
1738 0 : if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
1739 0 : return NULL;
1740 0 : x = PyFloat_AsDouble(ox);
1741 0 : y = PyFloat_AsDouble(oy);
1742 0 : if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1743 0 : return NULL;
1744 : /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
1745 0 : if (Py_IS_INFINITY(x))
1746 0 : return PyFloat_FromDouble(fabs(x));
1747 0 : if (Py_IS_INFINITY(y))
1748 0 : return PyFloat_FromDouble(fabs(y));
1749 0 : errno = 0;
1750 : PyFPE_START_PROTECT("in math_hypot", return 0);
1751 0 : r = hypot(x, y);
1752 : PyFPE_END_PROTECT(r);
1753 0 : if (Py_IS_NAN(r)) {
1754 0 : if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
1755 0 : errno = EDOM;
1756 : else
1757 0 : errno = 0;
1758 : }
1759 0 : else if (Py_IS_INFINITY(r)) {
1760 0 : if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
1761 0 : errno = ERANGE;
1762 : else
1763 0 : errno = 0;
1764 : }
1765 0 : if (errno && is_error(r))
1766 0 : return NULL;
1767 : else
1768 0 : return PyFloat_FromDouble(r);
1769 : }
1770 :
1771 : PyDoc_STRVAR(math_hypot_doc,
1772 : "hypot(x, y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
1773 :
1774 : /* pow can't use math_2, but needs its own wrapper: the problem is
1775 : that an infinite result can arise either as a result of overflow
1776 : (in which case OverflowError should be raised) or as a result of
1777 : e.g. 0.**-5. (for which ValueError needs to be raised.)
1778 : */
1779 :
1780 : static PyObject *
1781 0 : math_pow(PyObject *self, PyObject *args)
1782 : {
1783 : PyObject *ox, *oy;
1784 : double r, x, y;
1785 : int odd_y;
1786 :
1787 0 : if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
1788 0 : return NULL;
1789 0 : x = PyFloat_AsDouble(ox);
1790 0 : y = PyFloat_AsDouble(oy);
1791 0 : if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
1792 0 : return NULL;
1793 :
1794 : /* deal directly with IEEE specials, to cope with problems on various
1795 : platforms whose semantics don't exactly match C99 */
1796 0 : r = 0.; /* silence compiler warning */
1797 0 : if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
1798 0 : errno = 0;
1799 0 : if (Py_IS_NAN(x))
1800 0 : r = y == 0. ? 1. : x; /* NaN**0 = 1 */
1801 0 : else if (Py_IS_NAN(y))
1802 0 : r = x == 1. ? 1. : y; /* 1**NaN = 1 */
1803 0 : else if (Py_IS_INFINITY(x)) {
1804 0 : odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
1805 0 : if (y > 0.)
1806 0 : r = odd_y ? x : fabs(x);
1807 0 : else if (y == 0.)
1808 0 : r = 1.;
1809 : else /* y < 0. */
1810 0 : r = odd_y ? copysign(0., x) : 0.;
1811 : }
1812 0 : else if (Py_IS_INFINITY(y)) {
1813 0 : if (fabs(x) == 1.0)
1814 0 : r = 1.;
1815 0 : else if (y > 0. && fabs(x) > 1.0)
1816 0 : r = y;
1817 0 : else if (y < 0. && fabs(x) < 1.0) {
1818 0 : r = -y; /* result is +inf */
1819 0 : if (x == 0.) /* 0**-inf: divide-by-zero */
1820 0 : errno = EDOM;
1821 : }
1822 : else
1823 0 : r = 0.;
1824 : }
1825 : }
1826 : else {
1827 : /* let libm handle finite**finite */
1828 0 : errno = 0;
1829 : PyFPE_START_PROTECT("in math_pow", return 0);
1830 0 : r = pow(x, y);
1831 : PyFPE_END_PROTECT(r);
1832 : /* a NaN result should arise only from (-ve)**(finite
1833 : non-integer); in this case we want to raise ValueError. */
1834 0 : if (!Py_IS_FINITE(r)) {
1835 0 : if (Py_IS_NAN(r)) {
1836 0 : errno = EDOM;
1837 : }
1838 : /*
1839 : an infinite result here arises either from:
1840 : (A) (+/-0.)**negative (-> divide-by-zero)
1841 : (B) overflow of x**y with x and y finite
1842 : */
1843 0 : else if (Py_IS_INFINITY(r)) {
1844 0 : if (x == 0.)
1845 0 : errno = EDOM;
1846 : else
1847 0 : errno = ERANGE;
1848 : }
1849 : }
1850 : }
1851 :
1852 0 : if (errno && is_error(r))
1853 0 : return NULL;
1854 : else
1855 0 : return PyFloat_FromDouble(r);
1856 : }
1857 :
1858 : PyDoc_STRVAR(math_pow_doc,
1859 : "pow(x, y)\n\nReturn x**y (x to the power of y).");
1860 :
1861 : static const double degToRad = Py_MATH_PI / 180.0;
1862 : static const double radToDeg = 180.0 / Py_MATH_PI;
1863 :
1864 : static PyObject *
1865 0 : math_degrees(PyObject *self, PyObject *arg)
1866 : {
1867 0 : double x = PyFloat_AsDouble(arg);
1868 0 : if (x == -1.0 && PyErr_Occurred())
1869 0 : return NULL;
1870 0 : return PyFloat_FromDouble(x * radToDeg);
1871 : }
1872 :
1873 : PyDoc_STRVAR(math_degrees_doc,
1874 : "degrees(x)\n\n\
1875 : Convert angle x from radians to degrees.");
1876 :
1877 : static PyObject *
1878 0 : math_radians(PyObject *self, PyObject *arg)
1879 : {
1880 0 : double x = PyFloat_AsDouble(arg);
1881 0 : if (x == -1.0 && PyErr_Occurred())
1882 0 : return NULL;
1883 0 : return PyFloat_FromDouble(x * degToRad);
1884 : }
1885 :
1886 : PyDoc_STRVAR(math_radians_doc,
1887 : "radians(x)\n\n\
1888 : Convert angle x from degrees to radians.");
1889 :
1890 : static PyObject *
1891 0 : math_isfinite(PyObject *self, PyObject *arg)
1892 : {
1893 0 : double x = PyFloat_AsDouble(arg);
1894 0 : if (x == -1.0 && PyErr_Occurred())
1895 0 : return NULL;
1896 0 : return PyBool_FromLong((long)Py_IS_FINITE(x));
1897 : }
1898 :
1899 : PyDoc_STRVAR(math_isfinite_doc,
1900 : "isfinite(x) -> bool\n\n\
1901 : Return True if x is neither an infinity nor a NaN, and False otherwise.");
1902 :
1903 : static PyObject *
1904 0 : math_isnan(PyObject *self, PyObject *arg)
1905 : {
1906 0 : double x = PyFloat_AsDouble(arg);
1907 0 : if (x == -1.0 && PyErr_Occurred())
1908 0 : return NULL;
1909 0 : return PyBool_FromLong((long)Py_IS_NAN(x));
1910 : }
1911 :
1912 : PyDoc_STRVAR(math_isnan_doc,
1913 : "isnan(x) -> bool\n\n\
1914 : Return True if x is a NaN (not a number), and False otherwise.");
1915 :
1916 : static PyObject *
1917 0 : math_isinf(PyObject *self, PyObject *arg)
1918 : {
1919 0 : double x = PyFloat_AsDouble(arg);
1920 0 : if (x == -1.0 && PyErr_Occurred())
1921 0 : return NULL;
1922 0 : return PyBool_FromLong((long)Py_IS_INFINITY(x));
1923 : }
1924 :
1925 : PyDoc_STRVAR(math_isinf_doc,
1926 : "isinf(x) -> bool\n\n\
1927 : Return True if x is a positive or negative infinity, and False otherwise.");
1928 :
1929 : static PyMethodDef math_methods[] = {
1930 : {"acos", math_acos, METH_O, math_acos_doc},
1931 : {"acosh", math_acosh, METH_O, math_acosh_doc},
1932 : {"asin", math_asin, METH_O, math_asin_doc},
1933 : {"asinh", math_asinh, METH_O, math_asinh_doc},
1934 : {"atan", math_atan, METH_O, math_atan_doc},
1935 : {"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
1936 : {"atanh", math_atanh, METH_O, math_atanh_doc},
1937 : {"ceil", math_ceil, METH_O, math_ceil_doc},
1938 : {"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
1939 : {"cos", math_cos, METH_O, math_cos_doc},
1940 : {"cosh", math_cosh, METH_O, math_cosh_doc},
1941 : {"degrees", math_degrees, METH_O, math_degrees_doc},
1942 : {"erf", math_erf, METH_O, math_erf_doc},
1943 : {"erfc", math_erfc, METH_O, math_erfc_doc},
1944 : {"exp", math_exp, METH_O, math_exp_doc},
1945 : {"expm1", math_expm1, METH_O, math_expm1_doc},
1946 : {"fabs", math_fabs, METH_O, math_fabs_doc},
1947 : {"factorial", math_factorial, METH_O, math_factorial_doc},
1948 : {"floor", math_floor, METH_O, math_floor_doc},
1949 : {"fmod", math_fmod, METH_VARARGS, math_fmod_doc},
1950 : {"frexp", math_frexp, METH_O, math_frexp_doc},
1951 : {"fsum", math_fsum, METH_O, math_fsum_doc},
1952 : {"gamma", math_gamma, METH_O, math_gamma_doc},
1953 : {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
1954 : {"isfinite", math_isfinite, METH_O, math_isfinite_doc},
1955 : {"isinf", math_isinf, METH_O, math_isinf_doc},
1956 : {"isnan", math_isnan, METH_O, math_isnan_doc},
1957 : {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
1958 : {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
1959 : {"log", math_log, METH_VARARGS, math_log_doc},
1960 : {"log1p", math_log1p, METH_O, math_log1p_doc},
1961 : {"log10", math_log10, METH_O, math_log10_doc},
1962 : {"log2", math_log2, METH_O, math_log2_doc},
1963 : {"modf", math_modf, METH_O, math_modf_doc},
1964 : {"pow", math_pow, METH_VARARGS, math_pow_doc},
1965 : {"radians", math_radians, METH_O, math_radians_doc},
1966 : {"sin", math_sin, METH_O, math_sin_doc},
1967 : {"sinh", math_sinh, METH_O, math_sinh_doc},
1968 : {"sqrt", math_sqrt, METH_O, math_sqrt_doc},
1969 : {"tan", math_tan, METH_O, math_tan_doc},
1970 : {"tanh", math_tanh, METH_O, math_tanh_doc},
1971 : {"trunc", math_trunc, METH_O, math_trunc_doc},
1972 : {NULL, NULL} /* sentinel */
1973 : };
1974 :
1975 :
1976 : PyDoc_STRVAR(module_doc,
1977 : "This module is always available. It provides access to the\n"
1978 : "mathematical functions defined by the C standard.");
1979 :
1980 :
1981 : static struct PyModuleDef mathmodule = {
1982 : PyModuleDef_HEAD_INIT,
1983 : "math",
1984 : module_doc,
1985 : -1,
1986 : math_methods,
1987 : NULL,
1988 : NULL,
1989 : NULL,
1990 : NULL
1991 : };
1992 :
1993 : PyMODINIT_FUNC
1994 0 : PyInit_math(void)
1995 : {
1996 : PyObject *m;
1997 :
1998 0 : m = PyModule_Create(&mathmodule);
1999 0 : if (m == NULL)
2000 0 : goto finally;
2001 :
2002 0 : PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
2003 0 : PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
2004 :
2005 : finally:
2006 0 : return m;
2007 : }
|
