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: 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 m_erf_series(double x) 429 { 430 double x2, acc, fk, result; 431 int i, saved_errno; 432 433 x2 = x * x; 434 acc = 0.0; 435 fk = (double)ERF_SERIES_TERMS + 0.5; 436 for (i = 0; i < ERF_SERIES_TERMS; i++) { 437 acc = 2.0 + x2 * acc / fk; 438 fk -= 1.0; 439 } 440 /* Make sure the exp call doesn't affect errno; 441 see m_erfc_contfrac for more. */ 442 saved_errno = errno; 443 result = acc * x * exp(-x2) / sqrtpi; 444 errno = saved_errno; 445 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 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 if (x >= ERFC_CONTFRAC_CUTOFF) 464 return 0.0; 465 466 x2 = x*x; 467 a = 0.0; 468 da = 0.5; 469 p = 1.0; p_last = 0.0; 470 q = da + x2; q_last = 1.0; 471 for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) { 472 double temp; 473 a += da; 474 da += 2.0; 475 b = da + x2; 476 temp = p; p = b*p - a*p_last; p_last = temp; 477 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 saved_errno = errno; 482 result = p / q * x * exp(-x2) / sqrtpi; 483 errno = saved_errno; 484 return result; 485 } 486 487 /* Error function erf(x), for general x */ 488 489 static double 490 m_erf(double x) 491 { 492 double absx, cf; 493 494 if (Py_IS_NAN(x)) 495 return x; 496 absx = fabs(x); 497 if (absx < ERF_SERIES_CUTOFF) 498 return m_erf_series(x); 499 else { 500 cf = m_erfc_contfrac(absx); 501 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 m_erfc(double x) 509 { 510 double absx, cf; 511 512 if (Py_IS_NAN(x)) 513 return x; 514 absx = fabs(x); 515 if (absx < ERF_SERIES_CUTOFF) 516 return 1.0 - m_erf_series(x); 517 else { 518 cf = m_erfc_contfrac(absx); 519 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 m_atan2(double y, double x) 533 { 534 if (Py_IS_NAN(x) || Py_IS_NAN(y)) 535 return Py_NAN; 536 if (Py_IS_INFINITY(y)) { 537 if (Py_IS_INFINITY(x)) { 538 if (copysign(1., x) == 1.) 539 /* atan2(+-inf, +inf) == +-pi/4 */ 540 return copysign(0.25*Py_MATH_PI, y); 541 else 542 /* atan2(+-inf, -inf) == +-pi*3/4 */ 543 return copysign(0.75*Py_MATH_PI, y); 544 } 545 /* atan2(+-inf, x) == +-pi/2 for finite x */ 546 return copysign(0.5*Py_MATH_PI, y); 547 } 548 if (Py_IS_INFINITY(x) || y == 0.) { 549 if (copysign(1., x) == 1.) 550 /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ 551 return copysign(0., y); 552 else 553 /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ 554 return copysign(Py_MATH_PI, y); 555 } 556 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 m_log(double x) 568 { 569 if (Py_IS_FINITE(x)) { 570 if (x > 0.0) 571 return log(x); 572 errno = EDOM; 573 if (x == 0.0) 574 return -Py_HUGE_VAL; /* log(0) = -inf */ 575 else 576 return Py_NAN; /* log(-ve) = nan */ 577 } 578 else if (Py_IS_NAN(x)) 579 return x; /* log(nan) = nan */ 580 else if (x > 0.0) 581 return x; /* log(inf) = inf */ 582 else { 583 errno = EDOM; 584 return Py_NAN; /* log(-inf) = nan */ 585 } 586 } 587 588 static double 589 m_log10(double x) 590 { 591 if (Py_IS_FINITE(x)) { 592 if (x > 0.0) 593 return log10(x); 594 errno = EDOM; 595 if (x == 0.0) 596 return -Py_HUGE_VAL; /* log10(0) = -inf */ 597 else 598 return Py_NAN; /* log10(-ve) = nan */ 599 } 600 else if (Py_IS_NAN(x)) 601 return x; /* log10(nan) = nan */ 602 else if (x > 0.0) 603 return x; /* log10(inf) = inf */ 604 else { 605 errno = EDOM; 606 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 is_error(double x) 617 { 618 int result = 1; /* presumption of guilt */ 619 assert(errno); /* non-zero errno is a precondition for calling */ 620 if (errno == EDOM) 621 PyErr_SetString(PyExc_ValueError, "math domain error"); 622 623 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 if (fabs(x) < 1.0) 640 result = 0; 641 else 642 PyErr_SetString(PyExc_OverflowError, 643 "math range error"); 644 } 645 else 646 /* Unexpected math error */ 647 PyErr_SetFromErrno(PyExc_ValueError); 648 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 math_1(PyObject *arg, double (*func) (double), int can_overflow) 683 { 684 double x, r; 685 x = PyFloat_AsDouble(arg); 686 if (x == -1.0 && PyErr_Occurred()) 687 return NULL; 688 errno = 0; 689 PyFPE_START_PROTECT("in math_1", return 0); 690 r = (*func)(x); 691 PyFPE_END_PROTECT(r); 692 if (Py_IS_NAN(r)) { 693 if (!Py_IS_NAN(x)) 694 errno = EDOM; 695 else 696 errno = 0; 697 } 698 else if (Py_IS_INFINITY(r)) { 699 if (Py_IS_FINITE(x)) 700 errno = can_overflow ? ERANGE : EDOM; 701 else 702 errno = 0; 703 } 704 if (errno && is_error(r)) 705 return NULL; 706 else 707 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 math_1a(PyObject *arg, double (*func) (double)) 716 { 717 double x, r; 718 x = PyFloat_AsDouble(arg); 719 if (x == -1.0 && PyErr_Occurred()) 720 return NULL; 721 errno = 0; 722 PyFPE_START_PROTECT("in math_1a", return 0); 723 r = (*func)(x); 724 PyFPE_END_PROTECT(r); 725 if (errno && is_error(r)) 726 return NULL; 727 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 math_2(PyObject *args, double (*func) (double, double), char *funcname) 759 { 760 PyObject *ox, *oy; 761 double x, y, r; 762 if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy)) 763 return NULL; 764 x = PyFloat_AsDouble(ox); 765 y = PyFloat_AsDouble(oy); 766 if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) 767 return NULL; 768 errno = 0; 769 PyFPE_START_PROTECT("in math_2", return 0); 770 r = (*func)(x, y); 771 PyFPE_END_PROTECT(r); 772 if (Py_IS_NAN(r)) { 773 if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) 774 errno = EDOM; 775 else 776 errno = 0; 777 } 778 else if (Py_IS_INFINITY(r)) { 779 if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) 780 errno = ERANGE; 781 else 782 errno = 0; 783 } 784 if (errno && is_error(r)) 785 return NULL; 786 else 787 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 FUNC1(acos, acos, 0, 809 "acos(x)\n\nReturn the arc cosine (measured in radians) of x.") 810 FUNC1(acosh, m_acosh, 0, 811 "acosh(x)\n\nReturn the inverse hyperbolic cosine of x.") 812 FUNC1(asin, asin, 0, 813 "asin(x)\n\nReturn the arc sine (measured in radians) of x.") 814 FUNC1(asinh, m_asinh, 0, 815 "asinh(x)\n\nReturn the inverse hyperbolic sine of x.") 816 FUNC1(atan, atan, 0, 817 "atan(x)\n\nReturn the arc tangent (measured in radians) of x.") 818 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 FUNC1(atanh, m_atanh, 0, 822 "atanh(x)\n\nReturn the inverse hyperbolic tangent of x.") 823 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 FUNC2(copysign, copysign, 827 "copysign(x, y)\n\nReturn x with the sign of y.") 828 FUNC1(cos, cos, 0, 829 "cos(x)\n\nReturn the cosine of x (measured in radians).") 830 FUNC1(cosh, cosh, 1, 831 "cosh(x)\n\nReturn the hyperbolic cosine of x.") 832 FUNC1A(erf, m_erf, 833 "erf(x)\n\nError function at x.") 834 FUNC1A(erfc, m_erfc, 835 "erfc(x)\n\nComplementary error function at x.") 836 FUNC1(exp, exp, 1, 837 "exp(x)\n\nReturn e raised to the power of x.") 838 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 FUNC1(fabs, fabs, 0, 843 "fabs(x)\n\nReturn the absolute value of the float x.") 844 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 FUNC1A(gamma, m_tgamma, 848 "gamma(x)\n\nGamma function at x.") 849 FUNC1A(lgamma, m_lgamma, 850 "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.") 851 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 FUNC1(sin, sin, 0, 855 "sin(x)\n\nReturn the sine of x (measured in radians).") 856 FUNC1(sinh, sinh, 1, 857 "sinh(x)\n\nReturn the hyperbolic sine of x.") 858 FUNC1(sqrt, sqrt, 0, 859 "sqrt(x)\n\nReturn the square root of x.") 860 FUNC1(tan, tan, 0, 861 "tan(x)\n\nReturn the tangent of x (measured in radians).") 862 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 _fsum_realloc(double **p_ptr, Py_ssize_t n, 906 double *ps, Py_ssize_t *m_ptr) 907 { 908 void *v = NULL; 909 Py_ssize_t m = *m_ptr; 910 911 m += m; /* double */ 912 if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) { 913 double *p = *p_ptr; 914 if (p == ps) { 915 v = PyMem_Malloc(sizeof(double) * m); 916 if (v != NULL) 917 memcpy(v, ps, sizeof(double) * n); 918 } 919 else 920 v = PyMem_Realloc(p, sizeof(double) * m); 921 } 922 if (v == NULL) { /* size overflow or no memory */ 923 PyErr_SetString(PyExc_MemoryError, "math.fsum partials"); 924 return 1; 925 } 926 *p_ptr = (double*) v; 927 *m_ptr = m; 928 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 math_fsum(PyObject *self, PyObject *seq) 963 { 964 PyObject *item, *iter, *sum = NULL; 965 Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; 966 double x, y, t, ps[NUM_PARTIALS], *p = ps; 967 double xsave, special_sum = 0.0, inf_sum = 0.0; 968 volatile double hi, yr, lo; 969 970 iter = PyObject_GetIter(seq); 971 if (iter == NULL) 972 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 item = PyIter_Next(iter); 982 if (item == NULL) { 983 if (PyErr_Occurred()) 984 goto _fsum_error; 985 break; 986 } 987 x = PyFloat_AsDouble(item); 988 Py_DECREF(item); 989 if (PyErr_Occurred()) 990 goto _fsum_error; 991 992 xsave = x; 993 for (i = j = 0; j < n; j++) { /* for y in partials */ 994 y = p[j]; 995 if (fabs(x) < fabs(y)) { 996 t = x; x = y; y = t; 997 } 998 hi = x + y; 999 yr = hi - x; 1000 lo = y - yr; 1001 if (lo != 0.0) 1002 p[i++] = lo; 1003 x = hi; 1004 } 1005 1006 n = i; /* ps[i:] = [x] */ 1007 if (x != 0.0) { 1008 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 if (Py_IS_FINITE(xsave)) { 1014 PyErr_SetString(PyExc_OverflowError, 1015 "intermediate overflow in fsum"); 1016 goto _fsum_error; 1017 } 1018 if (Py_IS_INFINITY(xsave)) 1019 inf_sum += xsave; 1020 special_sum += xsave; 1021 /* reset partials */ 1022 n = 0; 1023 } 1024 else if (n >= m && _fsum_realloc(&p, n, ps, &m)) 1025 goto _fsum_error; 1026 else 1027 p[n++] = x; 1028 } 1029 } 1030 1031 if (special_sum != 0.0) { 1032 if (Py_IS_NAN(inf_sum)) 1033 PyErr_SetString(PyExc_ValueError, 1034 "-inf + inf in fsum"); 1035 else 1036 sum = PyFloat_FromDouble(special_sum); 1037 goto _fsum_error; 1038 } 1039 1040 hi = 0.0; 1041 if (n > 0) { 1042 hi = p[--n]; 1043 /* sum_exact(ps, hi) from the top, stop when the sum becomes 1044 inexact. */ 1045 while (n > 0) { 1046 x = hi; 1047 y = p[--n]; 1048 assert(fabs(y) < fabs(x)); 1049 hi = x + y; 1050 yr = hi - x; 1051 lo = y - yr; 1052 if (lo != 0.0) 1053 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 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || 1061 (lo > 0.0 && p[n-1] > 0.0))) { 1062 y = lo * 2.0; 1063 x = hi + y; 1064 yr = x - hi; 1065 if (y == yr) 1066 hi = x; 1067 } 1068 } 1069 sum = PyFloat_FromDouble(hi); 1070 1071 _fsum_error: 1072 PyFPE_END_PROTECT(hi) 1073 Py_DECREF(iter); 1074 if (p != ps) 1075 PyMem_Free(p); 1076 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 math_factorial(PyObject *self, PyObject *arg) 1088 { 1089 long i, x; 1090 PyObject *result, *iobj, *newresult; 1091 1092 if (PyFloat_Check(arg)) { 1093 PyObject *lx; 1094 double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); 1095 if (!(Py_IS_FINITE(dx) && dx == floor(dx))) { 1096 PyErr_SetString(PyExc_ValueError, 1097 "factorial() only accepts integral values"); 1098 return NULL; 1099 } 1100 lx = PyLong_FromDouble(dx); 1101 if (lx == NULL) 1102 return NULL; 1103 x = PyLong_AsLong(lx); 1104 Py_DECREF(lx); 1105 } 1106 else 1107 x = PyInt_AsLong(arg); 1108 1109 if (x == -1 && PyErr_Occurred()) 1110 return NULL; 1111 if (x < 0) { 1112 PyErr_SetString(PyExc_ValueError, 1113 "factorial() not defined for negative values"); 1114 return NULL; 1115 } 1116 1117 result = (PyObject *)PyInt_FromLong(1); 1118 if (result == NULL) 1119 return NULL; 1120 for (i=1 ; i<=x ; i++) { 1121 iobj = (PyObject *)PyInt_FromLong(i); 1122 if (iobj == NULL) 1123 goto error; 1124 newresult = PyNumber_Multiply(result, iobj); 1125 Py_DECREF(iobj); 1126 if (newresult == NULL) 1127 goto error; 1128 Py_DECREF(result); 1129 result = newresult; 1130 } 1131 return result; 1132 1133 error: 1134 Py_DECREF(result); 1135 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 math_trunc(PyObject *self, PyObject *number) 1145 { 1146 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 math_frexp(PyObject *self, PyObject *arg) 1156 { 1157 int i; 1158 double x = PyFloat_AsDouble(arg); 1159 if (x == -1.0 && PyErr_Occurred()) 1160 return NULL; 1161 /* deal with special cases directly, to sidestep platform 1162 differences */ 1163 if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { 1164 i = 0; 1165 } 1166 else { 1167 PyFPE_START_PROTECT("in math_frexp", return 0); 1168 x = frexp(x, &i); 1169 PyFPE_END_PROTECT(x); 1170 } 1171 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 math_ldexp(PyObject *self, PyObject *args) 1183 { 1184 double x, r; 1185 PyObject *oexp; 1186 long exp; 1187 int overflow; 1188 if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp)) 1189 return NULL; 1190 1191 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 exp = PyLong_AsLongAndOverflow(oexp, &overflow); 1195 if (exp == -1 && PyErr_Occurred()) 1196 return NULL; 1197 if (overflow) 1198 exp = overflow < 0 ? LONG_MIN : LONG_MAX; 1199 } 1200 else { 1201 PyErr_SetString(PyExc_TypeError, 1202 "Expected an int or long as second argument " 1203 "to ldexp."); 1204 return NULL; 1205 } 1206 1207 if (x == 0. || !Py_IS_FINITE(x)) { 1208 /* NaNs, zeros and infinities are returned unchanged */ 1209 r = x; 1210 errno = 0; 1211 } else if (exp > INT_MAX) { 1212 /* overflow */ 1213 r = copysign(Py_HUGE_VAL, x); 1214 errno = ERANGE; 1215 } else if (exp < INT_MIN) { 1216 /* underflow to +-0 */ 1217 r = copysign(0., x); 1218 errno = 0; 1219 } else { 1220 errno = 0; 1221 PyFPE_START_PROTECT("in math_ldexp", return 0); 1222 r = ldexp(x, (int)exp); 1223 PyFPE_END_PROTECT(r); 1224 if (Py_IS_INFINITY(r)) 1225 errno = ERANGE; 1226 } 1227 1228 if (errno && is_error(r)) 1229 return NULL; 1230 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 math_modf(PyObject *self, PyObject *arg) 1239 { 1240 double y, x = PyFloat_AsDouble(arg); 1241 if (x == -1.0 && PyErr_Occurred()) 1242 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 if (!Py_IS_FINITE(x)) { 1246 if (Py_IS_INFINITY(x)) 1247 return Py_BuildValue("(dd)", copysign(0., x), x); 1248 else if (Py_IS_NAN(x)) 1249 return Py_BuildValue("(dd)", x, x); 1250 } 1251 1252 errno = 0; 1253 PyFPE_START_PROTECT("in math_modf", return 0); 1254 x = modf(x, &y); 1255 PyFPE_END_PROTECT(x); 1256 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 loghelper(PyObject* arg, double (*func)(double), char *funcname) 1276 { 1277 /* If it is long, do it ourselves. */ 1278 if (PyLong_Check(arg)) { 1279 double x, result; 1280 Py_ssize_t e; 1281 1282 /* Negative or zero inputs give a ValueError. */ 1283 if (Py_SIZE(arg) <= 0) { 1284 PyErr_SetString(PyExc_ValueError, 1285 "math domain error"); 1286 return NULL; 1287 } 1288 1289 x = PyLong_AsDouble(arg); 1290 if (x == -1.0 && PyErr_Occurred()) { 1291 if (!PyErr_ExceptionMatches(PyExc_OverflowError)) 1292 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 PyErr_Clear(); 1296 x = _PyLong_Frexp((PyLongObject *)arg, &e); 1297 if (x == -1.0 && PyErr_Occurred()) 1298 return NULL; 1299 /* Value is ~= x * 2**e, so the log ~= log(x) + log(2) * e. */ 1300 result = func(x) + func(2.0) * e; 1301 } 1302 else 1303 /* Successfully converted x to a double. */ 1304 result = func(x); 1305 return PyFloat_FromDouble(result); 1306 } 1307 1308 /* Else let libm handle it by itself. */ 1309 return math_1(arg, func, 0); 1310 } 1311 1312 static PyObject * 1313 math_log(PyObject *self, PyObject *args) 1314 { 1315 PyObject *arg; 1316 PyObject *base = NULL; 1317 PyObject *num, *den; 1318 PyObject *ans; 1319 1320 if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base)) 1321 return NULL; 1322 1323 num = loghelper(arg, m_log, "log"); 1324 if (num == NULL || base == NULL) 1325 return num; 1326 1327 den = loghelper(base, m_log, "log"); 1328 if (den == NULL) { 1329 Py_DECREF(num); 1330 return NULL; 1331 } 1332 1333 ans = PyNumber_Divide(num, den); 1334 Py_DECREF(num); 1335 Py_DECREF(den); 1336 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 math_log10(PyObject *self, PyObject *arg) 1346 { 1347 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 math_fmod(PyObject *self, PyObject *args) 1355 { 1356 PyObject *ox, *oy; 1357 double r, x, y; 1358 if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy)) 1359 return NULL; 1360 x = PyFloat_AsDouble(ox); 1361 y = PyFloat_AsDouble(oy); 1362 if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) 1363 return NULL; 1364 /* fmod(x, +/-Inf) returns x for finite x. */ 1365 if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) 1366 return PyFloat_FromDouble(x); 1367 errno = 0; 1368 PyFPE_START_PROTECT("in math_fmod", return 0); 1369 r = fmod(x, y); 1370 PyFPE_END_PROTECT(r); 1371 if (Py_IS_NAN(r)) { 1372 if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) 1373 errno = EDOM; 1374 else 1375 errno = 0; 1376 } 1377 if (errno && is_error(r)) 1378 return NULL; 1379 else 1380 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 math_hypot(PyObject *self, PyObject *args) 1389 { 1390 PyObject *ox, *oy; 1391 double r, x, y; 1392 if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy)) 1393 return NULL; 1394 x = PyFloat_AsDouble(ox); 1395 y = PyFloat_AsDouble(oy); 1396 if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) 1397 return NULL; 1398 /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */ 1399 if (Py_IS_INFINITY(x)) 1400 return PyFloat_FromDouble(fabs(x)); 1401 if (Py_IS_INFINITY(y)) 1402 return PyFloat_FromDouble(fabs(y)); 1403 errno = 0; 1404 PyFPE_START_PROTECT("in math_hypot", return 0); 1405 r = hypot(x, y); 1406 PyFPE_END_PROTECT(r); 1407 if (Py_IS_NAN(r)) { 1408 if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) 1409 errno = EDOM; 1410 else 1411 errno = 0; 1412 } 1413 else if (Py_IS_INFINITY(r)) { 1414 if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) 1415 errno = ERANGE; 1416 else 1417 errno = 0; 1418 } 1419 if (errno && is_error(r)) 1420 return NULL; 1421 else 1422 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 math_pow(PyObject *self, PyObject *args) 1436 { 1437 PyObject *ox, *oy; 1438 double r, x, y; 1439 int odd_y; 1440 1441 if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy)) 1442 return NULL; 1443 x = PyFloat_AsDouble(ox); 1444 y = PyFloat_AsDouble(oy); 1445 if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) 1446 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 r = 0.; /* silence compiler warning */ 1451 if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) { 1452 errno = 0; 1453 if (Py_IS_NAN(x)) 1454 r = y == 0. ? 1. : x; /* NaN**0 = 1 */ 1455 else if (Py_IS_NAN(y)) 1456 r = x == 1. ? 1. : y; /* 1**NaN = 1 */ 1457 else if (Py_IS_INFINITY(x)) { 1458 odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0; 1459 if (y > 0.) 1460 r = odd_y ? x : fabs(x); 1461 else if (y == 0.) 1462 r = 1.; 1463 else /* y < 0. */ 1464 r = odd_y ? copysign(0., x) : 0.; 1465 } 1466 else if (Py_IS_INFINITY(y)) { 1467 if (fabs(x) == 1.0) 1468 r = 1.; 1469 else if (y > 0. && fabs(x) > 1.0) 1470 r = y; 1471 else if (y < 0. && fabs(x) < 1.0) { 1472 r = -y; /* result is +inf */ 1473 if (x == 0.) /* 0**-inf: divide-by-zero */ 1474 errno = EDOM; 1475 } 1476 else 1477 r = 0.; 1478 } 1479 } 1480 else { 1481 /* let libm handle finite**finite */ 1482 errno = 0; 1483 PyFPE_START_PROTECT("in math_pow", return 0); 1484 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 if (!Py_IS_FINITE(r)) { 1489 if (Py_IS_NAN(r)) { 1490 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 else if (Py_IS_INFINITY(r)) { 1498 if (x == 0.) 1499 errno = EDOM; 1500 else 1501 errno = ERANGE; 1502 } 1503 } 1504 } 1505 1506 if (errno && is_error(r)) 1507 return NULL; 1508 else 1509 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 math_degrees(PyObject *self, PyObject *arg) 1520 { 1521 double x = PyFloat_AsDouble(arg); 1522 if (x == -1.0 && PyErr_Occurred()) 1523 return NULL; 1524 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 math_radians(PyObject *self, PyObject *arg) 1533 { 1534 double x = PyFloat_AsDouble(arg); 1535 if (x == -1.0 && PyErr_Occurred()) 1536 return NULL; 1537 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 math_isnan(PyObject *self, PyObject *arg) 1546 { 1547 double x = PyFloat_AsDouble(arg); 1548 if (x == -1.0 && PyErr_Occurred()) 1549 return NULL; 1550 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 math_isinf(PyObject *self, PyObject *arg) 1559 { 1560 double x = PyFloat_AsDouble(arg); 1561 if (x == -1.0 && PyErr_Occurred()) 1562 return NULL; 1563 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 initmath(void) 1621 { 1622 PyObject *m; 1623 1624 m = Py_InitModule3("math", math_methods, module_doc); 1625 if (m == NULL) 1626 goto finally; 1627 1628 PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); 1629 PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); 1630 1631 finally: 1632 return; 1633 } 1634
