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