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