Modules/mathmodule.c

   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 "longintrepr.h" /* just for SHIFT */
  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 /* Call is_error when errno != 0, and where x is the result libm
  64  * returned.  is_error will usually set up an exception and return
  65  * true (1), but may return false (0) without setting up an exception.
  66  */
  67 static int
  68 is_error(double x)
  69 {
  70         int result = 1; /* presumption of guilt */
  71         assert(errno);  /* non-zero errno is a precondition for calling */
  72         if (errno == EDOM)
  73                 PyErr_SetString(PyExc_ValueError, "math domain error");
  74
  75         else if (errno == ERANGE) {
  76                 /* ANSI C generally requires libm functions to set ERANGE
  77                  * on overflow, but also generally *allows* them to set
  78                  * ERANGE on underflow too.  There's no consistency about
  79                  * the latter across platforms.
  80                  * Alas, C99 never requires that errno be set.
  81                  * Here we suppress the underflow errors (libm functions
  82                  * should return a zero on underflow, and +- HUGE_VAL on
  83                  * overflow, so testing the result for zero suffices to
  84                  * distinguish the cases).
  85                  *
  86                  * On some platforms (Ubuntu/ia64) it seems that errno can be
  87                  * set to ERANGE for subnormal results that do *not* underflow
  88                  * to zero.  So to be safe, we'll ignore ERANGE whenever the
  89                  * function result is less than one in absolute value.
  90                  */
  91                 if (fabs(x) < 1.0)
  92                         result = 0;
  93                 else
  94                         PyErr_SetString(PyExc_OverflowError,
  95                                         "math range error");
  96         }
  97         else
  98                 /* Unexpected math error */
  99                 PyErr_SetFromErrno(PyExc_ValueError);
 100         return result;
 101 }
 102
 103 /*
 104    wrapper for atan2 that deals directly with special cases before
 105    delegating to the platform libm for the remaining cases.  This
 106    is necessary to get consistent behaviour across platforms.
 107    Windows, FreeBSD and alpha Tru64 are amongst platforms that don't
 108    always follow C99.
 109 */
 110
 111 static double
 112 m_atan2(double y, double x)
 113 {
 114         if (Py_IS_NAN(x) || Py_IS_NAN(y))
 115                 return Py_NAN;
 116         if (Py_IS_INFINITY(y)) {
 117                 if (Py_IS_INFINITY(x)) {
 118                         if (copysign(1., x) == 1.)
 119                                 /* atan2(+-inf, +inf) == +-pi/4 */
 120                                 return copysign(0.25*Py_MATH_PI, y);
 121                         else
 122                                 /* atan2(+-inf, -inf) == +-pi*3/4 */
 123                                 return copysign(0.75*Py_MATH_PI, y);
 124                 }
 125                 /* atan2(+-inf, x) == +-pi/2 for finite x */
 126                 return copysign(0.5*Py_MATH_PI, y);
 127         }
 128         if (Py_IS_INFINITY(x) || y == 0.) {
 129                 if (copysign(1., x) == 1.)
 130                         /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */
 131                         return copysign(0., y);
 132                 else
 133                         /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */
 134                         return copysign(Py_MATH_PI, y);
 135         }
 136         return atan2(y, x);
 137 }
 138
 139 /*
 140     Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
 141     log(-ve), log(NaN).  Here are wrappers for log and log10 that deal with
 142     special values directly, passing positive non-special values through to
 143     the system log/log10.
 144  */
 145
 146 static double
 147 m_log(double x)
 148 {
 149         if (Py_IS_FINITE(x)) {
 150                 if (x > 0.0)
 151                         return log(x);
 152                 errno = EDOM;
 153                 if (x == 0.0)
 154                         return -Py_HUGE_VAL; /* log(0) = -inf */
 155                 else
 156                         return Py_NAN; /* log(-ve) = nan */
 157         }
 158         else if (Py_IS_NAN(x))
 159                 return x; /* log(nan) = nan */
 160         else if (x > 0.0)
 161                 return x; /* log(inf) = inf */
 162         else {
 163                 errno = EDOM;
 164                 return Py_NAN; /* log(-inf) = nan */
 165         }
 166 }
 167
 168 static double
 169 m_log10(double x)
 170 {
 171         if (Py_IS_FINITE(x)) {
 172                 if (x > 0.0)
 173                         return log10(x);
 174                 errno = EDOM;
 175                 if (x == 0.0)
 176                         return -Py_HUGE_VAL; /* log10(0) = -inf */
 177                 else
 178                         return Py_NAN; /* log10(-ve) = nan */
 179         }
 180         else if (Py_IS_NAN(x))
 181                 return x; /* log10(nan) = nan */
 182         else if (x > 0.0)
 183                 return x; /* log10(inf) = inf */
 184         else {
 185                 errno = EDOM;
 186                 return Py_NAN; /* log10(-inf) = nan */
 187         }
 188 }
 189
 190
 191 /*
 192    math_1 is used to wrap a libm function f that takes a double
 193    arguments and returns a double.
 194
 195    The error reporting follows these rules, which are designed to do
 196    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
 197    platforms.
 198
 199    - a NaN result from non-NaN inputs causes ValueError to be raised
 200    - an infinite result from finite inputs causes OverflowError to be
 201      raised if can_overflow is 1, or raises ValueError if can_overflow
 202      is 0.
 203    - if the result is finite and errno == EDOM then ValueError is
 204      raised
 205    - if the result is finite and nonzero and errno == ERANGE then
 206      OverflowError is raised
 207
 208    The last rule is used to catch overflow on platforms which follow
 209    C89 but for which HUGE_VAL is not an infinity.
 210
 211    For the majority of one-argument functions these rules are enough
 212    to ensure that Python's functions behave as specified in 'Annex F'
 213    of the C99 standard, with the 'invalid' and 'divide-by-zero'
 214    floating-point exceptions mapping to Python's ValueError and the
 215    'overflow' floating-point exception mapping to OverflowError.
 216    math_1 only works for functions that don't have singularities *and*
 217    the possibility of overflow; fortunately, that covers everything we
 218    care about right now.
 219 */
 220
 221 static PyObject *
 222 math_1(PyObject *arg, double (*func) (double), int can_overflow)
 223 {
 224         double x, r;
 225         x = PyFloat_AsDouble(arg);
 226         if (x == -1.0 && PyErr_Occurred())
 227                 return NULL;
 228         errno = 0;
 229         PyFPE_START_PROTECT("in math_1", return 0);
 230         r = (*func)(x);
 231         PyFPE_END_PROTECT(r);
 232         if (Py_IS_NAN(r)) {
 233                 if (!Py_IS_NAN(x))
 234                         errno = EDOM;
 235                 else
 236                         errno = 0;
 237         }
 238         else if (Py_IS_INFINITY(r)) {
 239                 if (Py_IS_FINITE(x))
 240                         errno = can_overflow ? ERANGE : EDOM;
 241                 else
 242                         errno = 0;
 243         }
 244         if (errno && is_error(r))
 245                 return NULL;
 246         else
 247                 return PyFloat_FromDouble(r);
 248 }
 249
 250 /*
 251    math_2 is used to wrap a libm function f that takes two double
 252    arguments and returns a double.
 253
 254    The error reporting follows these rules, which are designed to do
 255    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
 256    platforms.
 257
 258    - a NaN result from non-NaN inputs causes ValueError to be raised
 259    - an infinite result from finite inputs causes OverflowError to be
 260      raised.
 261    - if the result is finite and errno == EDOM then ValueError is
 262      raised
 263    - if the result is finite and nonzero and errno == ERANGE then
 264      OverflowError is raised
 265
 266    The last rule is used to catch overflow on platforms which follow
 267    C89 but for which HUGE_VAL is not an infinity.
 268
 269    For most two-argument functions (copysign, fmod, hypot, atan2)
 270    these rules are enough to ensure that Python's functions behave as
 271    specified in 'Annex F' of the C99 standard, with the 'invalid' and
 272    'divide-by-zero' floating-point exceptions mapping to Python's
 273    ValueError and the 'overflow' floating-point exception mapping to
 274    OverflowError.
 275 */
 276
 277 static PyObject *
 278 math_2(PyObject *args, double (*func) (double, double), char *funcname)
 279 {
 280         PyObject *ox, *oy;
 281         double x, y, r;
 282         if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
 283                 return NULL;
 284         x = PyFloat_AsDouble(ox);
 285         y = PyFloat_AsDouble(oy);
 286         if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
 287                 return NULL;
 288         errno = 0;
 289         PyFPE_START_PROTECT("in math_2", return 0);
 290         r = (*func)(x, y);
 291         PyFPE_END_PROTECT(r);
 292         if (Py_IS_NAN(r)) {
 293                 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
 294                         errno = EDOM;
 295                 else
 296                         errno = 0;
 297         }
 298         else if (Py_IS_INFINITY(r)) {
 299                 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
 300                         errno = ERANGE;
 301                 else
 302                         errno = 0;
 303         }
 304         if (errno && is_error(r))
 305                 return NULL;
 306         else
 307                 return PyFloat_FromDouble(r);
 308 }
 309
 310 #define FUNC1(funcname, func, can_overflow, docstring)                  \
 311         static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
 312                 return math_1(args, func, can_overflow);                    \
 313         }\
 314         PyDoc_STRVAR(math_##funcname##_doc, docstring);
 315
 316 #define FUNC2(funcname, func, docstring) \
 317         static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
 318                 return math_2(args, func, #funcname); \
 319         }\
 320         PyDoc_STRVAR(math_##funcname##_doc, docstring);
 321
 322 FUNC1(acos, acos, 0,
 323       "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
 324 FUNC1(acosh, acosh, 0,
 325       "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
 326 FUNC1(asin, asin, 0,
 327       "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
 328 FUNC1(asinh, asinh, 0,
 329       "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
 330 FUNC1(atan, atan, 0,
 331       "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
 332 FUNC2(atan2, m_atan2,
 333       "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
 334       "Unlike atan(y/x), the signs of both x and y are considered.")
 335 FUNC1(atanh, atanh, 0,
 336       "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
 337 FUNC1(ceil, ceil, 0,
 338       "ceil(x)\n\nReturn the ceiling of x as a float.\n"
 339       "This is the smallest integral value >= x.")
 340 FUNC2(copysign, copysign,
 341       "copysign(x,y)\n\nReturn x with the sign of y.")
 342 FUNC1(cos, cos, 0,
 343       "cos(x)\n\nReturn the cosine of x (measured in radians).")
 344 FUNC1(cosh, cosh, 1,
 345       "cosh(x)\n\nReturn the hyperbolic cosine of x.")
 346 FUNC1(exp, exp, 1,
 347       "exp(x)\n\nReturn e raised to the power of x.")
 348 FUNC1(fabs, fabs, 0,
 349       "fabs(x)\n\nReturn the absolute value of the float x.")
 350 FUNC1(floor, floor, 0,
 351       "floor(x)\n\nReturn the floor of x as a float.\n"
 352       "This is the largest integral value <= x.")
 353 FUNC1(log1p, log1p, 1,
 354       "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
 355       The result is computed in a way which is accurate for x near zero.")
 356 FUNC1(sin, sin, 0,
 357       "sin(x)\n\nReturn the sine of x (measured in radians).")
 358 FUNC1(sinh, sinh, 1,
 359       "sinh(x)\n\nReturn the hyperbolic sine of x.")
 360 FUNC1(sqrt, sqrt, 0,
 361       "sqrt(x)\n\nReturn the square root of x.")
 362 FUNC1(tan, tan, 0,
 363       "tan(x)\n\nReturn the tangent of x (measured in radians).")
 364 FUNC1(tanh, tanh, 0,
 365       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
 366
 367 /* Precision summation function as msum() by Raymond Hettinger in
 368    <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,
 369    enhanced with the exact partials sum and roundoff from Mark
 370    Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.
 371    See those links for more details, proofs and other references.
 372
 373    Note 1: IEEE 754R floating point semantics are assumed,
 374    but the current implementation does not re-establish special
 375    value semantics across iterations (i.e. handling -Inf + Inf).
 376
 377    Note 2:  No provision is made for intermediate overflow handling;
 378    therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while
 379    sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the
 380    overflow of the first partial sum.
 381
 382    Note 3: The intermediate values lo, yr, and hi are declared volatile so
 383    aggressive compilers won't algebraically reduce lo to always be exactly 0.0.
 384    Also, the volatile declaration forces the values to be stored in memory as
 385    regular doubles instead of extended long precision (80-bit) values.  This
 386    prevents double rounding because any addition or subtraction of two doubles
 387    can be resolved exactly into double-sized hi and lo values.  As long as the
 388    hi value gets forced into a double before yr and lo are computed, the extra
 389    bits in downstream extended precision operations (x87 for example) will be
 390    exactly zero and therefore can be losslessly stored back into a double,
 391    thereby preventing double rounding.
 392
 393    Note 4: A similar implementation is in Modules/cmathmodule.c.
 394    Be sure to update both when making changes.
 395
 396    Note 5: The signature of math.fsum() differs from __builtin__.sum()
 397    because the start argument doesn't make sense in the context of
 398    accurate summation.  Since the partials table is collapsed before
 399    returning a result, sum(seq2, start=sum(seq1)) may not equal the
 400    accurate result returned by sum(itertools.chain(seq1, seq2)).
 401 */
 402
 403 #define NUM_PARTIALS  32  /* initial partials array size, on stack */
 404
 405 /* Extend the partials array p[] by doubling its size. */
 406 static int                          /* non-zero on error */
 407 _fsum_realloc(double **p_ptr, Py_ssize_t  n,
 408              double  *ps,    Py_ssize_t *m_ptr)
 409 {
 410         void *v = NULL;
 411         Py_ssize_t m = *m_ptr;
 412
 413         m += m;  /* double */
 414         if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {
 415                 double *p = *p_ptr;
 416                 if (p == ps) {
 417                         v = PyMem_Malloc(sizeof(double) * m);
 418                         if (v != NULL)
 419                                 memcpy(v, ps, sizeof(double) * n);
 420                 }
 421                 else
 422                         v = PyMem_Realloc(p, sizeof(double) * m);
 423         }
 424         if (v == NULL) {        /* size overflow or no memory */
 425                 PyErr_SetString(PyExc_MemoryError, "math.fsum partials");
 426                 return 1;
 427         }
 428         *p_ptr = (double*) v;
 429         *m_ptr = m;
 430         return 0;
 431 }
 432
 433 /* Full precision summation of a sequence of floats.
 434
 435    def msum(iterable):
 436        partials = []  # sorted, non-overlapping partial sums
 437        for x in iterable:
 438            i = 0
 439            for y in partials:
 440                if abs(x) < abs(y):
 441                    x, y = y, x
 442                hi = x + y
 443                lo = y - (hi - x)
 444                if lo:
 445                    partials[i] = lo
 446                    i += 1
 447                x = hi
 448            partials[i:] = [x]
 449        return sum_exact(partials)
 450
 451    Rounded x+y stored in hi with the roundoff stored in lo.  Together hi+lo
 452    are exactly equal to x+y.  The inner loop applies hi/lo summation to each
 453    partial so that the list of partial sums remains exact.
 454
 455    Sum_exact() adds the partial sums exactly and correctly rounds the final
 456    result (using the round-half-to-even rule).  The items in partials remain
 457    non-zero, non-special, non-overlapping and strictly increasing in
 458    magnitude, but possibly not all having the same sign.
 459
 460    Depends on IEEE 754 arithmetic guarantees and half-even rounding.
 461 */
 462
 463 static PyObject*
 464 math_fsum(PyObject *self, PyObject *seq)
 465 {
 466         PyObject *item, *iter, *sum = NULL;
 467         Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;
 468         double x, y, t, ps[NUM_PARTIALS], *p = ps;
 469         double xsave, special_sum = 0.0, inf_sum = 0.0;
 470         volatile double hi, yr, lo;
 471
 472         iter = PyObject_GetIter(seq);
 473         if (iter == NULL)
 474                 return NULL;
 475
 476         PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)
 477
 478         for(;;) {           /* for x in iterable */
 479                 assert(0 <= n && n <= m);
 480                 assert((m == NUM_PARTIALS && p == ps) ||
 481                        (m >  NUM_PARTIALS && p != NULL));
 482
 483                 item = PyIter_Next(iter);
 484                 if (item == NULL) {
 485                         if (PyErr_Occurred())
 486                                 goto _fsum_error;
 487                         break;
 488                 }
 489                 x = PyFloat_AsDouble(item);
 490                 Py_DECREF(item);
 491                 if (PyErr_Occurred())
 492                         goto _fsum_error;
 493
 494                 xsave = x;
 495                 for (i = j = 0; j < n; j++) {       /* for y in partials */
 496                         y = p[j];
 497                         if (fabs(x) < fabs(y)) {
 498                                 t = x; x = y; y = t;
 499                         }
 500                         hi = x + y;
 501                         yr = hi - x;
 502                         lo = y - yr;
 503                         if (lo != 0.0)
 504                                 p[i++] = lo;
 505                         x = hi;
 506                 }
 507
 508                 n = i;                              /* ps[i:] = [x] */
 509                 if (x != 0.0) {
 510                         if (! Py_IS_FINITE(x)) {
 511                                 /* a nonfinite x could arise either as
 512                                    a result of intermediate overflow, or
 513                                    as a result of a nan or inf in the
 514                                    summands */
 515                                 if (Py_IS_FINITE(xsave)) {
 516                                         PyErr_SetString(PyExc_OverflowError,
 517                                               "intermediate overflow in fsum");
 518                                         goto _fsum_error;
 519                                 }
 520                                 if (Py_IS_INFINITY(xsave))
 521                                         inf_sum += xsave;
 522                                 special_sum += xsave;
 523                                 /* reset partials */
 524                                 n = 0;
 525                         }
 526                         else if (n >= m && _fsum_realloc(&p, n, ps, &m))
 527                                 goto _fsum_error;
 528                         else
 529                                 p[n++] = x;
 530                 }
 531         }
 532
 533         if (special_sum != 0.0) {
 534                 if (Py_IS_NAN(inf_sum))
 535                         PyErr_SetString(PyExc_ValueError,
 536                                         "-inf + inf in fsum");
 537                 else
 538                         sum = PyFloat_FromDouble(special_sum);
 539                 goto _fsum_error;
 540         }
 541
 542         hi = 0.0;
 543         if (n > 0) {
 544                 hi = p[--n];
 545                 /* sum_exact(ps, hi) from the top, stop when the sum becomes
 546                    inexact. */
 547                 while (n > 0) {
 548                         x = hi;
 549                         y = p[--n];
 550                         assert(fabs(y) < fabs(x));
 551                         hi = x + y;
 552                         yr = hi - x;
 553                         lo = y - yr;
 554                         if (lo != 0.0)
 555                                 break;
 556                 }
 557                 /* Make half-even rounding work across multiple partials.
 558                    Needed so that sum([1e-16, 1, 1e16]) will round-up the last
 559                    digit to two instead of down to zero (the 1e-16 makes the 1
 560                    slightly closer to two).  With a potential 1 ULP rounding
 561                    error fixed-up, math.fsum() can guarantee commutativity. */
 562                 if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||
 563                               (lo > 0.0 && p[n-1] > 0.0))) {
 564                         y = lo * 2.0;
 565                         x = hi + y;
 566                         yr = x - hi;
 567                         if (y == yr)
 568                                 hi = x;
 569                 }
 570         }
 571         sum = PyFloat_FromDouble(hi);
 572
 573 _fsum_error:
 574         PyFPE_END_PROTECT(hi)
 575         Py_DECREF(iter);
 576         if (p != ps)
 577                 PyMem_Free(p);
 578         return sum;
 579 }
 580
 581 #undef NUM_PARTIALS
 582
 583 PyDoc_STRVAR(math_fsum_doc,
 584 "sum(iterable)\n\n\
 585 Return an accurate floating point sum of values in the iterable.\n\
 586 Assumes IEEE-754 floating point arithmetic.");
 587
 588 static PyObject *
 589 math_factorial(PyObject *self, PyObject *arg)
 590 {
 591         long i, x;
 592         PyObject *result, *iobj, *newresult;
 593
 594         if (PyFloat_Check(arg)) {
 595                 double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);
 596                 if (dx != floor(dx)) {
 597                         PyErr_SetString(PyExc_ValueError,
 598                                 "factorial() only accepts integral values");
 599                         return NULL;
 600                 }
 601         }
 602
 603         x = PyInt_AsLong(arg);
 604         if (x == -1 && PyErr_Occurred())
 605                 return NULL;
 606         if (x < 0) {
 607                 PyErr_SetString(PyExc_ValueError,
 608                         "factorial() not defined for negative values");
 609                 return NULL;
 610         }
 611
 612         result = (PyObject *)PyInt_FromLong(1);
 613         if (result == NULL)
 614                 return NULL;
 615         for (i=1 ; i<=x ; i++) {
 616                 iobj = (PyObject *)PyInt_FromLong(i);
 617                 if (iobj == NULL)
 618                         goto error;
 619                 newresult = PyNumber_Multiply(result, iobj);
 620                 Py_DECREF(iobj);
 621                 if (newresult == NULL)
 622                         goto error;
 623                 Py_DECREF(result);
 624                 result = newresult;
 625         }
 626         return result;
 627
 628 error:
 629         Py_DECREF(result);
 630         return NULL;
 631 }
 632
 633 PyDoc_STRVAR(math_factorial_doc,
 634 "factorial(x) -> Integral\n"
 635 "\n"
 636 "Find x!. Raise a ValueError if x is negative or non-integral.");
 637
 638 static PyObject *
 639 math_trunc(PyObject *self, PyObject *number)
 640 {
 641         return PyObject_CallMethod(number, "__trunc__", NULL);
 642 }
 643
 644 PyDoc_STRVAR(math_trunc_doc,
 645 "trunc(x:Real) -> Integral\n"
 646 "\n"
 647 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");
 648
 649 static PyObject *
 650 math_frexp(PyObject *self, PyObject *arg)
 651 {
 652         int i;
 653         double x = PyFloat_AsDouble(arg);
 654         if (x == -1.0 && PyErr_Occurred())
 655                 return NULL;
 656         /* deal with special cases directly, to sidestep platform
 657            differences */
 658         if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
 659                 i = 0;
 660         }
 661         else {
 662                 PyFPE_START_PROTECT("in math_frexp", return 0);
 663                 x = frexp(x, &i);
 664                 PyFPE_END_PROTECT(x);
 665         }
 666         return Py_BuildValue("(di)", x, i);
 667 }
 668
 669 PyDoc_STRVAR(math_frexp_doc,
 670 "frexp(x)\n"
 671 "\n"
 672 "Return the mantissa and exponent of x, as pair (m, e).\n"
 673 "m is a float and e is an int, such that x = m * 2.**e.\n"
 674 "If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.");
 675
 676 static PyObject *
 677 math_ldexp(PyObject *self, PyObject *args)
 678 {
 679         double x, r;
 680         PyObject *oexp;
 681         long exp;
 682         if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))
 683                 return NULL;
 684
 685         if (PyLong_Check(oexp)) {
 686                 /* on overflow, replace exponent with either LONG_MAX
 687                    or LONG_MIN, depending on the sign. */
 688                 exp = PyLong_AsLong(oexp);
 689                 if (exp == -1 && PyErr_Occurred()) {
 690                         if (PyErr_ExceptionMatches(PyExc_OverflowError)) {
 691                                 if (Py_SIZE(oexp) < 0) {
 692                                         exp = LONG_MIN;
 693                                 }
 694                                 else {
 695                                         exp = LONG_MAX;
 696                                 }
 697                                 PyErr_Clear();
 698                         }
 699                         else {
 700                                 /* propagate any unexpected exception */
 701                                 return NULL;
 702                         }
 703                 }
 704         }
 705         else if (PyInt_Check(oexp)) {
 706                 exp = PyInt_AS_LONG(oexp);
 707         }
 708         else {
 709                 PyErr_SetString(PyExc_TypeError,
 710                                 "Expected an int or long as second argument "
 711                                 "to ldexp.");
 712                 return NULL;
 713         }
 714
 715         if (x == 0. || !Py_IS_FINITE(x)) {
 716                 /* NaNs, zeros and infinities are returned unchanged */
 717                 r = x;
 718                 errno = 0;
 719         } else if (exp > INT_MAX) {
 720                 /* overflow */
 721                 r = copysign(Py_HUGE_VAL, x);
 722                 errno = ERANGE;
 723         } else if (exp < INT_MIN) {
 724                 /* underflow to +-0 */
 725                 r = copysign(0., x);
 726                 errno = 0;
 727         } else {
 728                 errno = 0;
 729                 PyFPE_START_PROTECT("in math_ldexp", return 0);
 730                 r = ldexp(x, (int)exp);
 731                 PyFPE_END_PROTECT(r);
 732                 if (Py_IS_INFINITY(r))
 733                         errno = ERANGE;
 734         }
 735
 736         if (errno && is_error(r))
 737                 return NULL;
 738         return PyFloat_FromDouble(r);
 739 }
 740
 741 PyDoc_STRVAR(math_ldexp_doc,
 742 "ldexp(x, i) -> x * (2**i)");
 743
 744 static PyObject *
 745 math_modf(PyObject *self, PyObject *arg)
 746 {
 747         double y, x = PyFloat_AsDouble(arg);
 748         if (x == -1.0 && PyErr_Occurred())
 749                 return NULL;
 750         /* some platforms don't do the right thing for NaNs and
 751            infinities, so we take care of special cases directly. */
 752         if (!Py_IS_FINITE(x)) {
 753                 if (Py_IS_INFINITY(x))
 754                         return Py_BuildValue("(dd)", copysign(0., x), x);
 755                 else if (Py_IS_NAN(x))
 756                         return Py_BuildValue("(dd)", x, x);
 757         }
 758
 759         errno = 0;
 760         PyFPE_START_PROTECT("in math_modf", return 0);
 761         x = modf(x, &y);
 762         PyFPE_END_PROTECT(x);
 763         return Py_BuildValue("(dd)", x, y);
 764 }
 765
 766 PyDoc_STRVAR(math_modf_doc,
 767 "modf(x)\n"
 768 "\n"
 769 "Return the fractional and integer parts of x.  Both results carry the sign\n"
 770 "of x and are floats.");
 771
 772 /* A decent logarithm is easy to compute even for huge longs, but libm can't
 773    do that by itself -- loghelper can.  func is log or log10, and name is
 774    "log" or "log10".  Note that overflow isn't possible:  a long can contain
 775    no more than INT_MAX * SHIFT bits, so has value certainly less than
 776    2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is
 777    small enough to fit in an IEEE single.  log and log10 are even smaller.
 778 */
 779
 780 static PyObject*
 781 loghelper(PyObject* arg, double (*func)(double), char *funcname)
 782 {
 783         /* If it is long, do it ourselves. */
 784         if (PyLong_Check(arg)) {
 785                 double x;
 786                 int e;
 787                 x = _PyLong_AsScaledDouble(arg, &e);
 788                 if (x <= 0.0) {
 789                         PyErr_SetString(PyExc_ValueError,
 790                                         "math domain error");
 791                         return NULL;
 792                 }
 793                 /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=
 794                    log(x) + log(2) * e * PyLong_SHIFT.
 795                    CAUTION:  e*PyLong_SHIFT may overflow using int arithmetic,
 796                    so force use of double. */
 797                 x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);
 798                 return PyFloat_FromDouble(x);
 799         }
 800
 801         /* Else let libm handle it by itself. */
 802         return math_1(arg, func, 0);
 803 }
 804
 805 static PyObject *
 806 math_log(PyObject *self, PyObject *args)
 807 {
 808         PyObject *arg;
 809         PyObject *base = NULL;
 810         PyObject *num, *den;
 811         PyObject *ans;
 812
 813         if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))
 814                 return NULL;
 815
 816         num = loghelper(arg, m_log, "log");
 817         if (num == NULL || base == NULL)
 818                 return num;
 819
 820         den = loghelper(base, m_log, "log");
 821         if (den == NULL) {
 822                 Py_DECREF(num);
 823                 return NULL;
 824         }
 825
 826         ans = PyNumber_Divide(num, den);
 827         Py_DECREF(num);
 828         Py_DECREF(den);
 829         return ans;
 830 }
 831
 832 PyDoc_STRVAR(math_log_doc,
 833 "log(x[, base]) -> the logarithm of x to the given base.\n\
 834 If the base not specified, returns the natural logarithm (base e) of x.");
 835
 836 static PyObject *
 837 math_log10(PyObject *self, PyObject *arg)
 838 {
 839         return loghelper(arg, m_log10, "log10");
 840 }
 841
 842 PyDoc_STRVAR(math_log10_doc,
 843 "log10(x) -> the base 10 logarithm of x.");
 844
 845 static PyObject *
 846 math_fmod(PyObject *self, PyObject *args)
 847 {
 848         PyObject *ox, *oy;
 849         double r, x, y;
 850         if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
 851                 return NULL;
 852         x = PyFloat_AsDouble(ox);
 853         y = PyFloat_AsDouble(oy);
 854         if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
 855                 return NULL;
 856         /* fmod(x, +/-Inf) returns x for finite x. */
 857         if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
 858                 return PyFloat_FromDouble(x);
 859         errno = 0;
 860         PyFPE_START_PROTECT("in math_fmod", return 0);
 861         r = fmod(x, y);
 862         PyFPE_END_PROTECT(r);
 863         if (Py_IS_NAN(r)) {
 864                 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
 865                         errno = EDOM;
 866                 else
 867                         errno = 0;
 868         }
 869         if (errno && is_error(r))
 870                 return NULL;
 871         else
 872                 return PyFloat_FromDouble(r);
 873 }
 874
 875 PyDoc_STRVAR(math_fmod_doc,
 876 "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
 877 "  x % y may differ.");
 878
 879 static PyObject *
 880 math_hypot(PyObject *self, PyObject *args)
 881 {
 882         PyObject *ox, *oy;
 883         double r, x, y;
 884         if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
 885                 return NULL;
 886         x = PyFloat_AsDouble(ox);
 887         y = PyFloat_AsDouble(oy);
 888         if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
 889                 return NULL;
 890         /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
 891         if (Py_IS_INFINITY(x))
 892                 return PyFloat_FromDouble(fabs(x));
 893         if (Py_IS_INFINITY(y))
 894                 return PyFloat_FromDouble(fabs(y));
 895         errno = 0;
 896         PyFPE_START_PROTECT("in math_hypot", return 0);
 897         r = hypot(x, y);
 898         PyFPE_END_PROTECT(r);
 899         if (Py_IS_NAN(r)) {
 900                 if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
 901                         errno = EDOM;
 902                 else
 903                         errno = 0;
 904         }
 905         else if (Py_IS_INFINITY(r)) {
 906                 if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
 907                         errno = ERANGE;
 908                 else
 909                         errno = 0;
 910         }
 911         if (errno && is_error(r))
 912                 return NULL;
 913         else
 914                 return PyFloat_FromDouble(r);
 915 }
 916
 917 PyDoc_STRVAR(math_hypot_doc,
 918 "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
 919
 920 /* pow can't use math_2, but needs its own wrapper: the problem is
 921    that an infinite result can arise either as a result of overflow
 922    (in which case OverflowError should be raised) or as a result of
 923    e.g. 0.**-5. (for which ValueError needs to be raised.)
 924 */
 925
 926 static PyObject *
 927 math_pow(PyObject *self, PyObject *args)
 928 {
 929         PyObject *ox, *oy;
 930         double r, x, y;
 931         int odd_y;
 932
 933         if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
 934                 return NULL;
 935         x = PyFloat_AsDouble(ox);
 936         y = PyFloat_AsDouble(oy);
 937         if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
 938                 return NULL;
 939
 940         /* deal directly with IEEE specials, to cope with problems on various
 941            platforms whose semantics don't exactly match C99 */
 942         r = 0.; /* silence compiler warning */
 943         if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {
 944                 errno = 0;
 945                 if (Py_IS_NAN(x))
 946                         r = y == 0. ? 1. : x; /* NaN**0 = 1 */
 947                 else if (Py_IS_NAN(y))
 948                         r = x == 1. ? 1. : y; /* 1**NaN = 1 */
 949                 else if (Py_IS_INFINITY(x)) {
 950                         odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;
 951                         if (y > 0.)
 952                                 r = odd_y ? x : fabs(x);
 953                         else if (y == 0.)
 954                                 r = 1.;
 955                         else /* y < 0. */
 956                                 r = odd_y ? copysign(0., x) : 0.;
 957                 }
 958                 else if (Py_IS_INFINITY(y)) {
 959                         if (fabs(x) == 1.0)
 960                                 r = 1.;
 961                         else if (y > 0. && fabs(x) > 1.0)
 962                                 r = y;
 963                         else if (y < 0. && fabs(x) < 1.0) {
 964                                 r = -y; /* result is +inf */
 965                                 if (x == 0.) /* 0**-inf: divide-by-zero */
 966                                         errno = EDOM;
 967                         }
 968                         else
 969                                 r = 0.;
 970                 }
 971         }
 972         else {
 973                 /* let libm handle finite**finite */
 974                 errno = 0;
 975                 PyFPE_START_PROTECT("in math_pow", return 0);
 976                 r = pow(x, y);
 977                 PyFPE_END_PROTECT(r);
 978                 /* a NaN result should arise only from (-ve)**(finite
 979                    non-integer); in this case we want to raise ValueError. */
 980                 if (!Py_IS_FINITE(r)) {
 981                         if (Py_IS_NAN(r)) {
 982                                 errno = EDOM;
 983                         }
 984                         /*
 985                            an infinite result here arises either from:
 986                            (A) (+/-0.)**negative (-> divide-by-zero)
 987                            (B) overflow of x**y with x and y finite
 988                         */
 989                         else if (Py_IS_INFINITY(r)) {
 990                                 if (x == 0.)
 991                                         errno = EDOM;
 992                                 else
 993                                         errno = ERANGE;
 994                         }
 995                 }
 996         }
 997
 998         if (errno && is_error(r))
 999                 return NULL;
1000         else
1001                 return PyFloat_FromDouble(r);
1002 }
1003
1004 PyDoc_STRVAR(math_pow_doc,
1005 "pow(x,y)\n\nReturn x**y (x to the power of y).");
1006
1007 static const double degToRad = Py_MATH_PI / 180.0;
1008 static const double radToDeg = 180.0 / Py_MATH_PI;
1009
1010 static PyObject *
1011 math_degrees(PyObject *self, PyObject *arg)
1012 {
1013         double x = PyFloat_AsDouble(arg);
1014         if (x == -1.0 && PyErr_Occurred())
1015                 return NULL;
1016         return PyFloat_FromDouble(x * radToDeg);
1017 }
1018
1019 PyDoc_STRVAR(math_degrees_doc,
1020 "degrees(x) -> converts angle x from radians to degrees");
1021
1022 static PyObject *
1023 math_radians(PyObject *self, PyObject *arg)
1024 {
1025         double x = PyFloat_AsDouble(arg);
1026         if (x == -1.0 && PyErr_Occurred())
1027                 return NULL;
1028         return PyFloat_FromDouble(x * degToRad);
1029 }
1030
1031 PyDoc_STRVAR(math_radians_doc,
1032 "radians(x) -> converts angle x from degrees to radians");
1033
1034 static PyObject *
1035 math_isnan(PyObject *self, PyObject *arg)
1036 {
1037         double x = PyFloat_AsDouble(arg);
1038         if (x == -1.0 && PyErr_Occurred())
1039                 return NULL;
1040         return PyBool_FromLong((long)Py_IS_NAN(x));
1041 }
1042
1043 PyDoc_STRVAR(math_isnan_doc,
1044 "isnan(x) -> bool\n\
1045 Checks if float x is not a number (NaN)");
1046
1047 static PyObject *
1048 math_isinf(PyObject *self, PyObject *arg)
1049 {
1050         double x = PyFloat_AsDouble(arg);
1051         if (x == -1.0 && PyErr_Occurred())
1052                 return NULL;
1053         return PyBool_FromLong((long)Py_IS_INFINITY(x));
1054 }
1055
1056 PyDoc_STRVAR(math_isinf_doc,
1057 "isinf(x) -> bool\n\
1058 Checks if float x is infinite (positive or negative)");
1059
1060 static PyMethodDef math_methods[] = {
1061         {"acos",        math_acos,      METH_O,         math_acos_doc},
1062         {"acosh",       math_acosh,     METH_O,         math_acosh_doc},
1063         {"asin",        math_asin,      METH_O,         math_asin_doc},
1064         {"asinh",       math_asinh,     METH_O,         math_asinh_doc},
1065         {"atan",        math_atan,      METH_O,         math_atan_doc},
1066         {"atan2",       math_atan2,     METH_VARARGS,   math_atan2_doc},
1067         {"atanh",       math_atanh,     METH_O,         math_atanh_doc},
1068         {"ceil",        math_ceil,      METH_O,         math_ceil_doc},
1069         {"copysign",    math_copysign,  METH_VARARGS,   math_copysign_doc},
1070         {"cos",         math_cos,       METH_O,         math_cos_doc},
1071         {"cosh",        math_cosh,      METH_O,         math_cosh_doc},
1072         {"degrees",     math_degrees,   METH_O,         math_degrees_doc},
1073         {"exp",         math_exp,       METH_O,         math_exp_doc},
1074         {"fabs",        math_fabs,      METH_O,         math_fabs_doc},
1075         {"factorial",   math_factorial, METH_O,         math_factorial_doc},
1076         {"floor",       math_floor,     METH_O,         math_floor_doc},
1077         {"fmod",        math_fmod,      METH_VARARGS,   math_fmod_doc},
1078         {"frexp",       math_frexp,     METH_O,         math_frexp_doc},
1079         {"fsum",        math_fsum,      METH_O,         math_fsum_doc},
1080         {"hypot",       math_hypot,     METH_VARARGS,   math_hypot_doc},
1081         {"isinf",       math_isinf,     METH_O,         math_isinf_doc},
1082         {"isnan",       math_isnan,     METH_O,         math_isnan_doc},
1083         {"ldexp",       math_ldexp,     METH_VARARGS,   math_ldexp_doc},
1084         {"log",         math_log,       METH_VARARGS,   math_log_doc},
1085         {"log1p",       math_log1p,     METH_O,         math_log1p_doc},
1086         {"log10",       math_log10,     METH_O,         math_log10_doc},
1087         {"modf",        math_modf,      METH_O,         math_modf_doc},
1088         {"pow",         math_pow,       METH_VARARGS,   math_pow_doc},
1089         {"radians",     math_radians,   METH_O,         math_radians_doc},
1090         {"sin",         math_sin,       METH_O,         math_sin_doc},
1091         {"sinh",        math_sinh,      METH_O,         math_sinh_doc},
1092         {"sqrt",        math_sqrt,      METH_O,         math_sqrt_doc},
1093         {"tan",         math_tan,       METH_O,         math_tan_doc},
1094         {"tanh",        math_tanh,      METH_O,         math_tanh_doc},
1095         {"trunc",       math_trunc,     METH_O,         math_trunc_doc},
1096         {NULL,          NULL}           /* sentinel */
1097 };
1098
1099
1100 PyDoc_STRVAR(module_doc,
1101 "This module is always available.  It provides access to the\n"
1102 "mathematical functions defined by the C standard.");
1103
1104 PyMODINIT_FUNC
1105 initmath(void)
1106 {
1107         PyObject *m;
1108
1109         m = Py_InitModule3("math", math_methods, module_doc);
1110         if (m == NULL)
1111                 goto finally;
1112
1113         PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
1114         PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
1115
1116     finally:
1117         return;
1118 }