diff -r ffa851df0825 -r 2fb8b9db1c86 symbian-qemu-0.9.1-12/python-2.6.1/Modules/mathmodule.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/symbian-qemu-0.9.1-12/python-2.6.1/Modules/mathmodule.c Fri Jul 31 15:01:17 2009 +0100 @@ -0,0 +1,1063 @@ +/* Math module -- standard C math library functions, pi and e */ + +/* Here are some comments from Tim Peters, extracted from the + discussion attached to http://bugs.python.org/issue1640. They + describe the general aims of the math module with respect to + special values, IEEE-754 floating-point exceptions, and Python + exceptions. + +These are the "spirit of 754" rules: + +1. If the mathematical result is a real number, but of magnitude too +large to approximate by a machine float, overflow is signaled and the +result is an infinity (with the appropriate sign). + +2. If the mathematical result is a real number, but of magnitude too +small to approximate by a machine float, underflow is signaled and the +result is a zero (with the appropriate sign). + +3. At a singularity (a value x such that the limit of f(y) as y +approaches x exists and is an infinity), "divide by zero" is signaled +and the result is an infinity (with the appropriate sign). This is +complicated a little by that the left-side and right-side limits may +not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 +from the positive or negative directions. In that specific case, the +sign of the zero determines the result of 1/0. + +4. At a point where a function has no defined result in the extended +reals (i.e., the reals plus an infinity or two), invalid operation is +signaled and a NaN is returned. + +And these are what Python has historically /tried/ to do (but not +always successfully, as platform libm behavior varies a lot): + +For #1, raise OverflowError. + +For #2, return a zero (with the appropriate sign if that happens by +accident ;-)). + +For #3 and #4, raise ValueError. It may have made sense to raise +Python's ZeroDivisionError in #3, but historically that's only been +raised for division by zero and mod by zero. + +*/ + +/* + In general, on an IEEE-754 platform the aim is to follow the C99 + standard, including Annex 'F', whenever possible. Where the + standard recommends raising the 'divide-by-zero' or 'invalid' + floating-point exceptions, Python should raise a ValueError. Where + the standard recommends raising 'overflow', Python should raise an + OverflowError. In all other circumstances a value should be + returned. + */ + +#include "Python.h" +#include "longintrepr.h" /* just for SHIFT */ + +#ifdef _OSF_SOURCE +/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */ +extern double copysign(double, double); +#endif + +/* Call is_error when errno != 0, and where x is the result libm + * returned. is_error will usually set up an exception and return + * true (1), but may return false (0) without setting up an exception. + */ +static int +is_error(double x) +{ + int result = 1; /* presumption of guilt */ + assert(errno); /* non-zero errno is a precondition for calling */ + if (errno == EDOM) + PyErr_SetString(PyExc_ValueError, "math domain error"); + + else if (errno == ERANGE) { + /* ANSI C generally requires libm functions to set ERANGE + * on overflow, but also generally *allows* them to set + * ERANGE on underflow too. There's no consistency about + * the latter across platforms. + * Alas, C99 never requires that errno be set. + * Here we suppress the underflow errors (libm functions + * should return a zero on underflow, and +- HUGE_VAL on + * overflow, so testing the result for zero suffices to + * distinguish the cases). + * + * On some platforms (Ubuntu/ia64) it seems that errno can be + * set to ERANGE for subnormal results that do *not* underflow + * to zero. So to be safe, we'll ignore ERANGE whenever the + * function result is less than one in absolute value. + */ + if (fabs(x) < 1.0) + result = 0; + else + PyErr_SetString(PyExc_OverflowError, + "math range error"); + } + else + /* Unexpected math error */ + PyErr_SetFromErrno(PyExc_ValueError); + return result; +} + +/* + wrapper for atan2 that deals directly with special cases before + delegating to the platform libm for the remaining cases. This + is necessary to get consistent behaviour across platforms. + Windows, FreeBSD and alpha Tru64 are amongst platforms that don't + always follow C99. +*/ + +static double +m_atan2(double y, double x) +{ + if (Py_IS_NAN(x) || Py_IS_NAN(y)) + return Py_NAN; + if (Py_IS_INFINITY(y)) { + if (Py_IS_INFINITY(x)) { + if (copysign(1., x) == 1.) + /* atan2(+-inf, +inf) == +-pi/4 */ + return copysign(0.25*Py_MATH_PI, y); + else + /* atan2(+-inf, -inf) == +-pi*3/4 */ + return copysign(0.75*Py_MATH_PI, y); + } + /* atan2(+-inf, x) == +-pi/2 for finite x */ + return copysign(0.5*Py_MATH_PI, y); + } + if (Py_IS_INFINITY(x) || y == 0.) { + if (copysign(1., x) == 1.) + /* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */ + return copysign(0., y); + else + /* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */ + return copysign(Py_MATH_PI, y); + } + return atan2(y, x); +} + +/* + math_1 is used to wrap a libm function f that takes a double + arguments and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised if can_overflow is 1, or raises ValueError if can_overflow + is 0. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For the majority of one-argument functions these rules are enough + to ensure that Python's functions behave as specified in 'Annex F' + of the C99 standard, with the 'invalid' and 'divide-by-zero' + floating-point exceptions mapping to Python's ValueError and the + 'overflow' floating-point exception mapping to OverflowError. + math_1 only works for functions that don't have singularities *and* + the possibility of overflow; fortunately, that covers everything we + care about right now. +*/ + +static PyObject * +math_1(PyObject *arg, double (*func) (double), int can_overflow) +{ + double x, r; + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + errno = 0; + PyFPE_START_PROTECT("in math_1", return 0); + r = (*func)(x); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x)) + errno = can_overflow ? ERANGE : EDOM; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +/* + math_2 is used to wrap a libm function f that takes two double + arguments and returns a double. + + The error reporting follows these rules, which are designed to do + the right thing on C89/C99 platforms and IEEE 754/non IEEE 754 + platforms. + + - a NaN result from non-NaN inputs causes ValueError to be raised + - an infinite result from finite inputs causes OverflowError to be + raised. + - if the result is finite and errno == EDOM then ValueError is + raised + - if the result is finite and nonzero and errno == ERANGE then + OverflowError is raised + + The last rule is used to catch overflow on platforms which follow + C89 but for which HUGE_VAL is not an infinity. + + For most two-argument functions (copysign, fmod, hypot, atan2) + these rules are enough to ensure that Python's functions behave as + specified in 'Annex F' of the C99 standard, with the 'invalid' and + 'divide-by-zero' floating-point exceptions mapping to Python's + ValueError and the 'overflow' floating-point exception mapping to + OverflowError. +*/ + +static PyObject * +math_2(PyObject *args, double (*func) (double, double), char *funcname) +{ + PyObject *ox, *oy; + double x, y, r; + if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + errno = 0; + PyFPE_START_PROTECT("in math_2", return 0); + r = (*func)(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +#define FUNC1(funcname, func, can_overflow, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_1(args, func, can_overflow); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + +#define FUNC2(funcname, func, docstring) \ + static PyObject * math_##funcname(PyObject *self, PyObject *args) { \ + return math_2(args, func, #funcname); \ + }\ + PyDoc_STRVAR(math_##funcname##_doc, docstring); + +FUNC1(acos, acos, 0, + "acos(x)\n\nReturn the arc cosine (measured in radians) of x.") +FUNC1(acosh, acosh, 0, + "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.") +FUNC1(asin, asin, 0, + "asin(x)\n\nReturn the arc sine (measured in radians) of x.") +FUNC1(asinh, asinh, 0, + "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.") +FUNC1(atan, atan, 0, + "atan(x)\n\nReturn the arc tangent (measured in radians) of x.") +FUNC2(atan2, m_atan2, + "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n" + "Unlike atan(y/x), the signs of both x and y are considered.") +FUNC1(atanh, atanh, 0, + "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.") +FUNC1(ceil, ceil, 0, + "ceil(x)\n\nReturn the ceiling of x as a float.\n" + "This is the smallest integral value >= x.") +FUNC2(copysign, copysign, + "copysign(x,y)\n\nReturn x with the sign of y.") +FUNC1(cos, cos, 0, + "cos(x)\n\nReturn the cosine of x (measured in radians).") +FUNC1(cosh, cosh, 1, + "cosh(x)\n\nReturn the hyperbolic cosine of x.") +FUNC1(exp, exp, 1, + "exp(x)\n\nReturn e raised to the power of x.") +FUNC1(fabs, fabs, 0, + "fabs(x)\n\nReturn the absolute value of the float x.") +FUNC1(floor, floor, 0, + "floor(x)\n\nReturn the floor of x as a float.\n" + "This is the largest integral value <= x.") +FUNC1(log1p, log1p, 1, + "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\ + The result is computed in a way which is accurate for x near zero.") +FUNC1(sin, sin, 0, + "sin(x)\n\nReturn the sine of x (measured in radians).") +FUNC1(sinh, sinh, 1, + "sinh(x)\n\nReturn the hyperbolic sine of x.") +FUNC1(sqrt, sqrt, 0, + "sqrt(x)\n\nReturn the square root of x.") +FUNC1(tan, tan, 0, + "tan(x)\n\nReturn the tangent of x (measured in radians).") +FUNC1(tanh, tanh, 0, + "tanh(x)\n\nReturn the hyperbolic tangent of x.") + +/* Precision summation function as msum() by Raymond Hettinger in + , + enhanced with the exact partials sum and roundoff from Mark + Dickinson's post at . + See those links for more details, proofs and other references. + + Note 1: IEEE 754R floating point semantics are assumed, + but the current implementation does not re-establish special + value semantics across iterations (i.e. handling -Inf + Inf). + + Note 2: No provision is made for intermediate overflow handling; + therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while + sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the + overflow of the first partial sum. + + Note 3: The intermediate values lo, yr, and hi are declared volatile so + aggressive compilers won't algebraically reduce lo to always be exactly 0.0. + Also, the volatile declaration forces the values to be stored in memory as + regular doubles instead of extended long precision (80-bit) values. This + prevents double rounding because any addition or subtraction of two doubles + can be resolved exactly into double-sized hi and lo values. As long as the + hi value gets forced into a double before yr and lo are computed, the extra + bits in downstream extended precision operations (x87 for example) will be + exactly zero and therefore can be losslessly stored back into a double, + thereby preventing double rounding. + + Note 4: A similar implementation is in Modules/cmathmodule.c. + Be sure to update both when making changes. + + Note 5: The signature of math.fsum() differs from __builtin__.sum() + because the start argument doesn't make sense in the context of + accurate summation. Since the partials table is collapsed before + returning a result, sum(seq2, start=sum(seq1)) may not equal the + accurate result returned by sum(itertools.chain(seq1, seq2)). +*/ + +#define NUM_PARTIALS 32 /* initial partials array size, on stack */ + +/* Extend the partials array p[] by doubling its size. */ +static int /* non-zero on error */ +_fsum_realloc(double **p_ptr, Py_ssize_t n, + double *ps, Py_ssize_t *m_ptr) +{ + void *v = NULL; + Py_ssize_t m = *m_ptr; + + m += m; /* double */ + if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) { + double *p = *p_ptr; + if (p == ps) { + v = PyMem_Malloc(sizeof(double) * m); + if (v != NULL) + memcpy(v, ps, sizeof(double) * n); + } + else + v = PyMem_Realloc(p, sizeof(double) * m); + } + if (v == NULL) { /* size overflow or no memory */ + PyErr_SetString(PyExc_MemoryError, "math.fsum partials"); + return 1; + } + *p_ptr = (double*) v; + *m_ptr = m; + return 0; +} + +/* Full precision summation of a sequence of floats. + + def msum(iterable): + partials = [] # sorted, non-overlapping partial sums + for x in iterable: + i = 0 + for y in partials: + if abs(x) < abs(y): + x, y = y, x + hi = x + y + lo = y - (hi - x) + if lo: + partials[i] = lo + i += 1 + x = hi + partials[i:] = [x] + return sum_exact(partials) + + Rounded x+y stored in hi with the roundoff stored in lo. Together hi+lo + are exactly equal to x+y. The inner loop applies hi/lo summation to each + partial so that the list of partial sums remains exact. + + Sum_exact() adds the partial sums exactly and correctly rounds the final + result (using the round-half-to-even rule). The items in partials remain + non-zero, non-special, non-overlapping and strictly increasing in + magnitude, but possibly not all having the same sign. + + Depends on IEEE 754 arithmetic guarantees and half-even rounding. +*/ + +static PyObject* +math_fsum(PyObject *self, PyObject *seq) +{ + PyObject *item, *iter, *sum = NULL; + Py_ssize_t i, j, n = 0, m = NUM_PARTIALS; + double x, y, t, ps[NUM_PARTIALS], *p = ps; + double xsave, special_sum = 0.0, inf_sum = 0.0; + volatile double hi, yr, lo; + + iter = PyObject_GetIter(seq); + if (iter == NULL) + return NULL; + + PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL) + + for(;;) { /* for x in iterable */ + assert(0 <= n && n <= m); + assert((m == NUM_PARTIALS && p == ps) || + (m > NUM_PARTIALS && p != NULL)); + + item = PyIter_Next(iter); + if (item == NULL) { + if (PyErr_Occurred()) + goto _fsum_error; + break; + } + x = PyFloat_AsDouble(item); + Py_DECREF(item); + if (PyErr_Occurred()) + goto _fsum_error; + + xsave = x; + for (i = j = 0; j < n; j++) { /* for y in partials */ + y = p[j]; + if (fabs(x) < fabs(y)) { + t = x; x = y; y = t; + } + hi = x + y; + yr = hi - x; + lo = y - yr; + if (lo != 0.0) + p[i++] = lo; + x = hi; + } + + n = i; /* ps[i:] = [x] */ + if (x != 0.0) { + if (! Py_IS_FINITE(x)) { + /* a nonfinite x could arise either as + a result of intermediate overflow, or + as a result of a nan or inf in the + summands */ + if (Py_IS_FINITE(xsave)) { + PyErr_SetString(PyExc_OverflowError, + "intermediate overflow in fsum"); + goto _fsum_error; + } + if (Py_IS_INFINITY(xsave)) + inf_sum += xsave; + special_sum += xsave; + /* reset partials */ + n = 0; + } + else if (n >= m && _fsum_realloc(&p, n, ps, &m)) + goto _fsum_error; + else + p[n++] = x; + } + } + + if (special_sum != 0.0) { + if (Py_IS_NAN(inf_sum)) + PyErr_SetString(PyExc_ValueError, + "-inf + inf in fsum"); + else + sum = PyFloat_FromDouble(special_sum); + goto _fsum_error; + } + + hi = 0.0; + if (n > 0) { + hi = p[--n]; + /* sum_exact(ps, hi) from the top, stop when the sum becomes + inexact. */ + while (n > 0) { + x = hi; + y = p[--n]; + assert(fabs(y) < fabs(x)); + hi = x + y; + yr = hi - x; + lo = y - yr; + if (lo != 0.0) + break; + } + /* Make half-even rounding work across multiple partials. + Needed so that sum([1e-16, 1, 1e16]) will round-up the last + digit to two instead of down to zero (the 1e-16 makes the 1 + slightly closer to two). With a potential 1 ULP rounding + error fixed-up, math.fsum() can guarantee commutativity. */ + if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) || + (lo > 0.0 && p[n-1] > 0.0))) { + y = lo * 2.0; + x = hi + y; + yr = x - hi; + if (y == yr) + hi = x; + } + } + sum = PyFloat_FromDouble(hi); + +_fsum_error: + PyFPE_END_PROTECT(hi) + Py_DECREF(iter); + if (p != ps) + PyMem_Free(p); + return sum; +} + +#undef NUM_PARTIALS + +PyDoc_STRVAR(math_fsum_doc, +"sum(iterable)\n\n\ +Return an accurate floating point sum of values in the iterable.\n\ +Assumes IEEE-754 floating point arithmetic."); + +static PyObject * +math_factorial(PyObject *self, PyObject *arg) +{ + long i, x; + PyObject *result, *iobj, *newresult; + + if (PyFloat_Check(arg)) { + double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg); + if (dx != floor(dx)) { + PyErr_SetString(PyExc_ValueError, + "factorial() only accepts integral values"); + return NULL; + } + } + + x = PyInt_AsLong(arg); + if (x == -1 && PyErr_Occurred()) + return NULL; + if (x < 0) { + PyErr_SetString(PyExc_ValueError, + "factorial() not defined for negative values"); + return NULL; + } + + result = (PyObject *)PyInt_FromLong(1); + if (result == NULL) + return NULL; + for (i=1 ; i<=x ; i++) { + iobj = (PyObject *)PyInt_FromLong(i); + if (iobj == NULL) + goto error; + newresult = PyNumber_Multiply(result, iobj); + Py_DECREF(iobj); + if (newresult == NULL) + goto error; + Py_DECREF(result); + result = newresult; + } + return result; + +error: + Py_DECREF(result); + return NULL; +} + +PyDoc_STRVAR(math_factorial_doc, "Return n!"); + +static PyObject * +math_trunc(PyObject *self, PyObject *number) +{ + return PyObject_CallMethod(number, "__trunc__", NULL); +} + +PyDoc_STRVAR(math_trunc_doc, +"trunc(x:Real) -> Integral\n" +"\n" +"Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method."); + +static PyObject * +math_frexp(PyObject *self, PyObject *arg) +{ + int i; + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + /* deal with special cases directly, to sidestep platform + differences */ + if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) { + i = 0; + } + else { + PyFPE_START_PROTECT("in math_frexp", return 0); + x = frexp(x, &i); + PyFPE_END_PROTECT(x); + } + return Py_BuildValue("(di)", x, i); +} + +PyDoc_STRVAR(math_frexp_doc, +"frexp(x)\n" +"\n" +"Return the mantissa and exponent of x, as pair (m, e).\n" +"m is a float and e is an int, such that x = m * 2.**e.\n" +"If x is 0, m and e are both 0. Else 0.5 <= abs(m) < 1.0."); + +static PyObject * +math_ldexp(PyObject *self, PyObject *args) +{ + double x, r; + PyObject *oexp; + long exp; + if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp)) + return NULL; + + if (PyLong_Check(oexp)) { + /* on overflow, replace exponent with either LONG_MAX + or LONG_MIN, depending on the sign. */ + exp = PyLong_AsLong(oexp); + if (exp == -1 && PyErr_Occurred()) { + if (PyErr_ExceptionMatches(PyExc_OverflowError)) { + if (Py_SIZE(oexp) < 0) { + exp = LONG_MIN; + } + else { + exp = LONG_MAX; + } + PyErr_Clear(); + } + else { + /* propagate any unexpected exception */ + return NULL; + } + } + } + else if (PyInt_Check(oexp)) { + exp = PyInt_AS_LONG(oexp); + } + else { + PyErr_SetString(PyExc_TypeError, + "Expected an int or long as second argument " + "to ldexp."); + return NULL; + } + + if (x == 0. || !Py_IS_FINITE(x)) { + /* NaNs, zeros and infinities are returned unchanged */ + r = x; + errno = 0; + } else if (exp > INT_MAX) { + /* overflow */ + r = copysign(Py_HUGE_VAL, x); + errno = ERANGE; + } else if (exp < INT_MIN) { + /* underflow to +-0 */ + r = copysign(0., x); + errno = 0; + } else { + errno = 0; + PyFPE_START_PROTECT("in math_ldexp", return 0); + r = ldexp(x, (int)exp); + PyFPE_END_PROTECT(r); + if (Py_IS_INFINITY(r)) + errno = ERANGE; + } + + if (errno && is_error(r)) + return NULL; + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_ldexp_doc, +"ldexp(x, i) -> x * (2**i)"); + +static PyObject * +math_modf(PyObject *self, PyObject *arg) +{ + double y, x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + /* some platforms don't do the right thing for NaNs and + infinities, so we take care of special cases directly. */ + if (!Py_IS_FINITE(x)) { + if (Py_IS_INFINITY(x)) + return Py_BuildValue("(dd)", copysign(0., x), x); + else if (Py_IS_NAN(x)) + return Py_BuildValue("(dd)", x, x); + } + + errno = 0; + PyFPE_START_PROTECT("in math_modf", return 0); + x = modf(x, &y); + PyFPE_END_PROTECT(x); + return Py_BuildValue("(dd)", x, y); +} + +PyDoc_STRVAR(math_modf_doc, +"modf(x)\n" +"\n" +"Return the fractional and integer parts of x. Both results carry the sign\n" +"of x. The integer part is returned as a real."); + +/* A decent logarithm is easy to compute even for huge longs, but libm can't + do that by itself -- loghelper can. func is log or log10, and name is + "log" or "log10". Note that overflow isn't possible: a long can contain + no more than INT_MAX * SHIFT bits, so has value certainly less than + 2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is + small enough to fit in an IEEE single. log and log10 are even smaller. +*/ + +static PyObject* +loghelper(PyObject* arg, double (*func)(double), char *funcname) +{ + /* If it is long, do it ourselves. */ + if (PyLong_Check(arg)) { + double x; + int e; + x = _PyLong_AsScaledDouble(arg, &e); + if (x <= 0.0) { + PyErr_SetString(PyExc_ValueError, + "math domain error"); + return NULL; + } + /* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~= + log(x) + log(2) * e * PyLong_SHIFT. + CAUTION: e*PyLong_SHIFT may overflow using int arithmetic, + so force use of double. */ + x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0); + return PyFloat_FromDouble(x); + } + + /* Else let libm handle it by itself. */ + return math_1(arg, func, 0); +} + +static PyObject * +math_log(PyObject *self, PyObject *args) +{ + PyObject *arg; + PyObject *base = NULL; + PyObject *num, *den; + PyObject *ans; + + if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base)) + return NULL; + + num = loghelper(arg, log, "log"); + if (num == NULL || base == NULL) + return num; + + den = loghelper(base, log, "log"); + if (den == NULL) { + Py_DECREF(num); + return NULL; + } + + ans = PyNumber_Divide(num, den); + Py_DECREF(num); + Py_DECREF(den); + return ans; +} + +PyDoc_STRVAR(math_log_doc, +"log(x[, base]) -> the logarithm of x to the given base.\n\ +If the base not specified, returns the natural logarithm (base e) of x."); + +static PyObject * +math_log10(PyObject *self, PyObject *arg) +{ + return loghelper(arg, log10, "log10"); +} + +PyDoc_STRVAR(math_log10_doc, +"log10(x) -> the base 10 logarithm of x."); + +static PyObject * +math_fmod(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + /* fmod(x, +/-Inf) returns x for finite x. */ + if (Py_IS_INFINITY(y) && Py_IS_FINITE(x)) + return PyFloat_FromDouble(x); + errno = 0; + PyFPE_START_PROTECT("in math_fmod", return 0); + r = fmod(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_fmod_doc, +"fmod(x,y)\n\nReturn fmod(x, y), according to platform C." +" x % y may differ."); + +static PyObject * +math_hypot(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + /* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */ + if (Py_IS_INFINITY(x)) + return PyFloat_FromDouble(fabs(x)); + if (Py_IS_INFINITY(y)) + return PyFloat_FromDouble(fabs(y)); + errno = 0; + PyFPE_START_PROTECT("in math_hypot", return 0); + r = hypot(x, y); + PyFPE_END_PROTECT(r); + if (Py_IS_NAN(r)) { + if (!Py_IS_NAN(x) && !Py_IS_NAN(y)) + errno = EDOM; + else + errno = 0; + } + else if (Py_IS_INFINITY(r)) { + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) + errno = ERANGE; + else + errno = 0; + } + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_hypot_doc, +"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y)."); + +/* pow can't use math_2, but needs its own wrapper: the problem is + that an infinite result can arise either as a result of overflow + (in which case OverflowError should be raised) or as a result of + e.g. 0.**-5. (for which ValueError needs to be raised.) +*/ + +static PyObject * +math_pow(PyObject *self, PyObject *args) +{ + PyObject *ox, *oy; + double r, x, y; + int odd_y; + + if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy)) + return NULL; + x = PyFloat_AsDouble(ox); + y = PyFloat_AsDouble(oy); + if ((x == -1.0 || y == -1.0) && PyErr_Occurred()) + return NULL; + + /* deal directly with IEEE specials, to cope with problems on various + platforms whose semantics don't exactly match C99 */ + r = 0.; /* silence compiler warning */ + if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) { + errno = 0; + if (Py_IS_NAN(x)) + r = y == 0. ? 1. : x; /* NaN**0 = 1 */ + else if (Py_IS_NAN(y)) + r = x == 1. ? 1. : y; /* 1**NaN = 1 */ + else if (Py_IS_INFINITY(x)) { + odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0; + if (y > 0.) + r = odd_y ? x : fabs(x); + else if (y == 0.) + r = 1.; + else /* y < 0. */ + r = odd_y ? copysign(0., x) : 0.; + } + else if (Py_IS_INFINITY(y)) { + if (fabs(x) == 1.0) + r = 1.; + else if (y > 0. && fabs(x) > 1.0) + r = y; + else if (y < 0. && fabs(x) < 1.0) { + r = -y; /* result is +inf */ + if (x == 0.) /* 0**-inf: divide-by-zero */ + errno = EDOM; + } + else + r = 0.; + } + } + else { + /* let libm handle finite**finite */ + errno = 0; + PyFPE_START_PROTECT("in math_pow", return 0); + r = pow(x, y); + PyFPE_END_PROTECT(r); + /* a NaN result should arise only from (-ve)**(finite + non-integer); in this case we want to raise ValueError. */ + if (!Py_IS_FINITE(r)) { + if (Py_IS_NAN(r)) { + errno = EDOM; + } + /* + an infinite result here arises either from: + (A) (+/-0.)**negative (-> divide-by-zero) + (B) overflow of x**y with x and y finite + */ + else if (Py_IS_INFINITY(r)) { + if (x == 0.) + errno = EDOM; + else + errno = ERANGE; + } + } + } + + if (errno && is_error(r)) + return NULL; + else + return PyFloat_FromDouble(r); +} + +PyDoc_STRVAR(math_pow_doc, +"pow(x,y)\n\nReturn x**y (x to the power of y)."); + +static const double degToRad = Py_MATH_PI / 180.0; +static const double radToDeg = 180.0 / Py_MATH_PI; + +static PyObject * +math_degrees(PyObject *self, PyObject *arg) +{ + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyFloat_FromDouble(x * radToDeg); +} + +PyDoc_STRVAR(math_degrees_doc, +"degrees(x) -> converts angle x from radians to degrees"); + +static PyObject * +math_radians(PyObject *self, PyObject *arg) +{ + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyFloat_FromDouble(x * degToRad); +} + +PyDoc_STRVAR(math_radians_doc, +"radians(x) -> converts angle x from degrees to radians"); + +static PyObject * +math_isnan(PyObject *self, PyObject *arg) +{ + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_NAN(x)); +} + +PyDoc_STRVAR(math_isnan_doc, +"isnan(x) -> bool\n\ +Checks if float x is not a number (NaN)"); + +static PyObject * +math_isinf(PyObject *self, PyObject *arg) +{ + double x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyBool_FromLong((long)Py_IS_INFINITY(x)); +} + +PyDoc_STRVAR(math_isinf_doc, +"isinf(x) -> bool\n\ +Checks if float x is infinite (positive or negative)"); + +static PyMethodDef math_methods[] = { + {"acos", math_acos, METH_O, math_acos_doc}, + {"acosh", math_acosh, METH_O, math_acosh_doc}, + {"asin", math_asin, METH_O, math_asin_doc}, + {"asinh", math_asinh, METH_O, math_asinh_doc}, + {"atan", math_atan, METH_O, math_atan_doc}, + {"atan2", math_atan2, METH_VARARGS, math_atan2_doc}, + {"atanh", math_atanh, METH_O, math_atanh_doc}, + {"ceil", math_ceil, METH_O, math_ceil_doc}, + {"copysign", math_copysign, METH_VARARGS, math_copysign_doc}, + {"cos", math_cos, METH_O, math_cos_doc}, + {"cosh", math_cosh, METH_O, math_cosh_doc}, + {"degrees", math_degrees, METH_O, math_degrees_doc}, + {"exp", math_exp, METH_O, math_exp_doc}, + {"fabs", math_fabs, METH_O, math_fabs_doc}, + {"factorial", math_factorial, METH_O, math_factorial_doc}, + {"floor", math_floor, METH_O, math_floor_doc}, + {"fmod", math_fmod, METH_VARARGS, math_fmod_doc}, + {"frexp", math_frexp, METH_O, math_frexp_doc}, + {"fsum", math_fsum, METH_O, math_fsum_doc}, + {"hypot", math_hypot, METH_VARARGS, math_hypot_doc}, + {"isinf", math_isinf, METH_O, math_isinf_doc}, + {"isnan", math_isnan, METH_O, math_isnan_doc}, + {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc}, + {"log", math_log, METH_VARARGS, math_log_doc}, + {"log1p", math_log1p, METH_O, math_log1p_doc}, + {"log10", math_log10, METH_O, math_log10_doc}, + {"modf", math_modf, METH_O, math_modf_doc}, + {"pow", math_pow, METH_VARARGS, math_pow_doc}, + {"radians", math_radians, METH_O, math_radians_doc}, + {"sin", math_sin, METH_O, math_sin_doc}, + {"sinh", math_sinh, METH_O, math_sinh_doc}, + {"sqrt", math_sqrt, METH_O, math_sqrt_doc}, + {"tan", math_tan, METH_O, math_tan_doc}, + {"tanh", math_tanh, METH_O, math_tanh_doc}, + {"trunc", math_trunc, METH_O, math_trunc_doc}, + {NULL, NULL} /* sentinel */ +}; + + +PyDoc_STRVAR(module_doc, +"This module is always available. It provides access to the\n" +"mathematical functions defined by the C standard."); + +PyMODINIT_FUNC +initmath(void) +{ + PyObject *m; + + m = Py_InitModule3("math", math_methods, module_doc); + if (m == NULL) + goto finally; + + PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI)); + PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E)); + + finally: + return; +}