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

 Generated by: LCOV version 1.10