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    math_1 is used to wrap a libm function f that takes a double

   141    arguments and returns a double.

   142 

   143    The error reporting follows these rules, which are designed to do

   144    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754

   145    platforms.

   146 

   147    - a NaN result from non-NaN inputs causes ValueError to be raised

   148    - an infinite result from finite inputs causes OverflowError to be

   149      raised if can_overflow is 1, or raises ValueError if can_overflow

   150      is 0.

   151    - if the result is finite and errno == EDOM then ValueError is

   152      raised

   153    - if the result is finite and nonzero and errno == ERANGE then

   154      OverflowError is raised

   155 

   156    The last rule is used to catch overflow on platforms which follow

   157    C89 but for which HUGE_VAL is not an infinity.

   158 

   159    For the majority of one-argument functions these rules are enough

   160    to ensure that Python's functions behave as specified in 'Annex F'

   161    of the C99 standard, with the 'invalid' and 'divide-by-zero'

   162    floating-point exceptions mapping to Python's ValueError and the

   163    'overflow' floating-point exception mapping to OverflowError.

   164    math_1 only works for functions that don't have singularities *and*

   165    the possibility of overflow; fortunately, that covers everything we

   166    care about right now.

   167 */

   168 

   169 static PyObject *

   170 math_1(PyObject *arg, double (*func) (double), int can_overflow)

   171 {

   172 	double x, r;

   173 	x = PyFloat_AsDouble(arg);

   174 	if (x == -1.0 && PyErr_Occurred())

   175 		return NULL;

   176 	errno = 0;

   177 	PyFPE_START_PROTECT("in math_1", return 0);

   178 	r = (*func)(x);

   179 	PyFPE_END_PROTECT(r);

   180 	if (Py_IS_NAN(r)) {

   181 		if (!Py_IS_NAN(x))

   182 			errno = EDOM;

   183 		else

   184 			errno = 0;

   185 	}

   186 	else if (Py_IS_INFINITY(r)) {

   187 		if (Py_IS_FINITE(x))

   188 			errno = can_overflow ? ERANGE : EDOM;

   189 		else

   190 			errno = 0;

   191 	}

   192 	if (errno && is_error(r))

   193 		return NULL;

   194 	else

   195 		return PyFloat_FromDouble(r);

   196 }

   197 

   198 /*

   199    math_2 is used to wrap a libm function f that takes two double

   200    arguments and returns a double.

   201 

   202    The error reporting follows these rules, which are designed to do

   203    the right thing on C89/C99 platforms and IEEE 754/non IEEE 754

   204    platforms.

   205 

   206    - a NaN result from non-NaN inputs causes ValueError to be raised

   207    - an infinite result from finite inputs causes OverflowError to be

   208      raised.

   209    - if the result is finite and errno == EDOM then ValueError is

   210      raised

   211    - if the result is finite and nonzero and errno == ERANGE then

   212      OverflowError is raised

   213 

   214    The last rule is used to catch overflow on platforms which follow

   215    C89 but for which HUGE_VAL is not an infinity.

   216 

   217    For most two-argument functions (copysign, fmod, hypot, atan2)

   218    these rules are enough to ensure that Python's functions behave as

   219    specified in 'Annex F' of the C99 standard, with the 'invalid' and

   220    'divide-by-zero' floating-point exceptions mapping to Python's

   221    ValueError and the 'overflow' floating-point exception mapping to

   222    OverflowError.

   223 */

   224 

   225 static PyObject *

   226 math_2(PyObject *args, double (*func) (double, double), char *funcname)

   227 {

   228 	PyObject *ox, *oy;

   229 	double x, y, r;

   230 	if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))

   231 		return NULL;

   232 	x = PyFloat_AsDouble(ox);

   233 	y = PyFloat_AsDouble(oy);

   234 	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())

   235 		return NULL;

   236 	errno = 0;

   237 	PyFPE_START_PROTECT("in math_2", return 0);

   238 	r = (*func)(x, y);

   239 	PyFPE_END_PROTECT(r);

   240 	if (Py_IS_NAN(r)) {

   241 		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))

   242 			errno = EDOM;

   243 		else

   244 			errno = 0;

   245 	}

   246 	else if (Py_IS_INFINITY(r)) {

   247 		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))

   248 			errno = ERANGE;

   249 		else

   250 			errno = 0;

   251 	}

   252 	if (errno && is_error(r))

   253 		return NULL;

   254 	else

   255 		return PyFloat_FromDouble(r);

   256 }

   257 

   258 #define FUNC1(funcname, func, can_overflow, docstring)			\

   259 	static PyObject * math_##funcname(PyObject *self, PyObject *args) { \

   260 		return math_1(args, func, can_overflow);		    \

   261 	}\

   262         PyDoc_STRVAR(math_##funcname##_doc, docstring);

   263 

   264 #define FUNC2(funcname, func, docstring) \

   265 	static PyObject * math_##funcname(PyObject *self, PyObject *args) { \

   266 		return math_2(args, func, #funcname); \

   267 	}\

   268         PyDoc_STRVAR(math_##funcname##_doc, docstring);

   269 

   270 FUNC1(acos, acos, 0,

   271       "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")

   272 FUNC1(acosh, acosh, 0,

   273       "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")

   274 FUNC1(asin, asin, 0,

   275       "asin(x)\n\nReturn the arc sine (measured in radians) of x.")

   276 FUNC1(asinh, asinh, 0,

   277       "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")

   278 FUNC1(atan, atan, 0,

   279       "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")

   280 FUNC2(atan2, m_atan2,

   281       "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"

   282       "Unlike atan(y/x), the signs of both x and y are considered.")

   283 FUNC1(atanh, atanh, 0,

   284       "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")

   285 FUNC1(ceil, ceil, 0,

   286       "ceil(x)\n\nReturn the ceiling of x as a float.\n"

   287       "This is the smallest integral value >= x.")

   288 FUNC2(copysign, copysign,

   289       "copysign(x,y)\n\nReturn x with the sign of y.")

   290 FUNC1(cos, cos, 0,

   291       "cos(x)\n\nReturn the cosine of x (measured in radians).")

   292 FUNC1(cosh, cosh, 1,

   293       "cosh(x)\n\nReturn the hyperbolic cosine of x.")

   294 FUNC1(exp, exp, 1,

   295       "exp(x)\n\nReturn e raised to the power of x.")

   296 FUNC1(fabs, fabs, 0,

   297       "fabs(x)\n\nReturn the absolute value of the float x.")

   298 FUNC1(floor, floor, 0,

   299       "floor(x)\n\nReturn the floor of x as a float.\n"

   300       "This is the largest integral value <= x.")

   301 FUNC1(log1p, log1p, 1,

   302       "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\

   303       The result is computed in a way which is accurate for x near zero.")

   304 FUNC1(sin, sin, 0,

   305       "sin(x)\n\nReturn the sine of x (measured in radians).")

   306 FUNC1(sinh, sinh, 1,

   307       "sinh(x)\n\nReturn the hyperbolic sine of x.")

   308 FUNC1(sqrt, sqrt, 0,

   309       "sqrt(x)\n\nReturn the square root of x.")

   310 FUNC1(tan, tan, 0,

   311       "tan(x)\n\nReturn the tangent of x (measured in radians).")

   312 FUNC1(tanh, tanh, 0,

   313       "tanh(x)\n\nReturn the hyperbolic tangent of x.")

   314 

   315 /* Precision summation function as msum() by Raymond Hettinger in

   316    <http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090>,

   317    enhanced with the exact partials sum and roundoff from Mark

   318    Dickinson's post at <http://bugs.python.org/file10357/msum4.py>.

   319    See those links for more details, proofs and other references.

   320 

   321    Note 1: IEEE 754R floating point semantics are assumed,

   322    but the current implementation does not re-establish special

   323    value semantics across iterations (i.e. handling -Inf + Inf).

   324 

   325    Note 2:  No provision is made for intermediate overflow handling;

   326    therefore, sum([1e+308, 1e-308, 1e+308]) returns 1e+308 while

   327    sum([1e+308, 1e+308, 1e-308]) raises an OverflowError due to the

   328    overflow of the first partial sum.

   329 

   330    Note 3: The intermediate values lo, yr, and hi are declared volatile so

   331    aggressive compilers won't algebraically reduce lo to always be exactly 0.0.

   332    Also, the volatile declaration forces the values to be stored in memory as

   333    regular doubles instead of extended long precision (80-bit) values.  This

   334    prevents double rounding because any addition or subtraction of two doubles

   335    can be resolved exactly into double-sized hi and lo values.  As long as the 

   336    hi value gets forced into a double before yr and lo are computed, the extra

   337    bits in downstream extended precision operations (x87 for example) will be

   338    exactly zero and therefore can be losslessly stored back into a double,

   339    thereby preventing double rounding.

   340 

   341    Note 4: A similar implementation is in Modules/cmathmodule.c.

   342    Be sure to update both when making changes.

   343 

   344    Note 5: The signature of math.fsum() differs from __builtin__.sum()

   345    because the start argument doesn't make sense in the context of

   346    accurate summation.  Since the partials table is collapsed before

   347    returning a result, sum(seq2, start=sum(seq1)) may not equal the

   348    accurate result returned by sum(itertools.chain(seq1, seq2)).

   349 */

   350 

   351 #define NUM_PARTIALS  32  /* initial partials array size, on stack */

   352 

   353 /* Extend the partials array p[] by doubling its size. */

   354 static int                          /* non-zero on error */

   355 _fsum_realloc(double **p_ptr, Py_ssize_t  n,

   356              double  *ps,    Py_ssize_t *m_ptr)

   357 {

   358 	void *v = NULL;

   359 	Py_ssize_t m = *m_ptr;

   360 

   361 	m += m;  /* double */

   362 	if (n < m && m < (PY_SSIZE_T_MAX / sizeof(double))) {

   363 		double *p = *p_ptr;

   364 		if (p == ps) {

   365 			v = PyMem_Malloc(sizeof(double) * m);

   366 			if (v != NULL)

   367 				memcpy(v, ps, sizeof(double) * n);

   368 		}

   369 		else

   370 			v = PyMem_Realloc(p, sizeof(double) * m);

   371 	}

   372 	if (v == NULL) {        /* size overflow or no memory */

   373 		PyErr_SetString(PyExc_MemoryError, "math.fsum partials");

   374 		return 1;

   375 	}

   376 	*p_ptr = (double*) v;

   377 	*m_ptr = m;

   378 	return 0;

   379 }

   380 

   381 /* Full precision summation of a sequence of floats.

   382 

   383    def msum(iterable):

   384        partials = []  # sorted, non-overlapping partial sums

   385        for x in iterable:

   386            i = 0

   387            for y in partials:

   388                if abs(x) < abs(y):

   389                    x, y = y, x

   390                hi = x + y

   391                lo = y - (hi - x)

   392                if lo:

   393                    partials[i] = lo

   394                    i += 1

   395                x = hi

   396            partials[i:] = [x]

   397        return sum_exact(partials)

   398 

   399    Rounded x+y stored in hi with the roundoff stored in lo.  Together hi+lo

   400    are exactly equal to x+y.  The inner loop applies hi/lo summation to each

   401    partial so that the list of partial sums remains exact.

   402 

   403    Sum_exact() adds the partial sums exactly and correctly rounds the final

   404    result (using the round-half-to-even rule).  The items in partials remain

   405    non-zero, non-special, non-overlapping and strictly increasing in

   406    magnitude, but possibly not all having the same sign.

   407 

   408    Depends on IEEE 754 arithmetic guarantees and half-even rounding.

   409 */

   410 

   411 static PyObject*

   412 math_fsum(PyObject *self, PyObject *seq)

   413 {

   414 	PyObject *item, *iter, *sum = NULL;

   415 	Py_ssize_t i, j, n = 0, m = NUM_PARTIALS;

   416 	double x, y, t, ps[NUM_PARTIALS], *p = ps;

   417 	double xsave, special_sum = 0.0, inf_sum = 0.0;

   418 	volatile double hi, yr, lo;

   419 

   420 	iter = PyObject_GetIter(seq);

   421 	if (iter == NULL)

   422 		return NULL;

   423 

   424 	PyFPE_START_PROTECT("fsum", Py_DECREF(iter); return NULL)

   425 

   426 	for(;;) {           /* for x in iterable */

   427 		assert(0 <= n && n <= m);

   428 		assert((m == NUM_PARTIALS && p == ps) ||

   429 		       (m >  NUM_PARTIALS && p != NULL));

   430 

   431 		item = PyIter_Next(iter);

   432 		if (item == NULL) {

   433 			if (PyErr_Occurred())

   434 				goto _fsum_error;

   435 			break;

   436 		}

   437 		x = PyFloat_AsDouble(item);

   438 		Py_DECREF(item);

   439 		if (PyErr_Occurred())

   440 			goto _fsum_error;

   441 

   442 		xsave = x;

   443 		for (i = j = 0; j < n; j++) {       /* for y in partials */

   444 			y = p[j];

   445 			if (fabs(x) < fabs(y)) {

   446 				t = x; x = y; y = t;

   447 			}

   448 			hi = x + y;

   449 			yr = hi - x;

   450 			lo = y - yr;

   451 			if (lo != 0.0)

   452 				p[i++] = lo;

   453 			x = hi;

   454 		}

   455 

   456 		n = i;                              /* ps[i:] = [x] */

   457 		if (x != 0.0) {

   458 			if (! Py_IS_FINITE(x)) {

   459 				/* a nonfinite x could arise either as

   460 				   a result of intermediate overflow, or

   461 				   as a result of a nan or inf in the

   462 				   summands */

   463 				if (Py_IS_FINITE(xsave)) {

   464 					PyErr_SetString(PyExc_OverflowError,

   465 					      "intermediate overflow in fsum");

   466 					goto _fsum_error;

   467 				}

   468 				if (Py_IS_INFINITY(xsave))

   469 					inf_sum += xsave;

   470 				special_sum += xsave;

   471 				/* reset partials */

   472 				n = 0;

   473 			}

   474 			else if (n >= m && _fsum_realloc(&p, n, ps, &m))

   475 				goto _fsum_error;

   476 			else

   477 				p[n++] = x;

   478 		}

   479 	}

   480 

   481 	if (special_sum != 0.0) {

   482 		if (Py_IS_NAN(inf_sum))

   483 			PyErr_SetString(PyExc_ValueError,

   484 					"-inf + inf in fsum");

   485 		else

   486 			sum = PyFloat_FromDouble(special_sum);

   487 		goto _fsum_error;

   488 	}

   489 

   490 	hi = 0.0;

   491 	if (n > 0) {

   492 		hi = p[--n];

   493 		/* sum_exact(ps, hi) from the top, stop when the sum becomes

   494 		   inexact. */

   495 		while (n > 0) {

   496 			x = hi;

   497 			y = p[--n];

   498 			assert(fabs(y) < fabs(x));

   499 			hi = x + y;

   500 			yr = hi - x;

   501 			lo = y - yr;

   502 			if (lo != 0.0)

   503 				break;

   504 		}

   505 		/* Make half-even rounding work across multiple partials.

   506 		   Needed so that sum([1e-16, 1, 1e16]) will round-up the last

   507 		   digit to two instead of down to zero (the 1e-16 makes the 1

   508 		   slightly closer to two).  With a potential 1 ULP rounding

   509 		   error fixed-up, math.fsum() can guarantee commutativity. */

   510 		if (n > 0 && ((lo < 0.0 && p[n-1] < 0.0) ||

   511 			      (lo > 0.0 && p[n-1] > 0.0))) {

   512 			y = lo * 2.0;

   513 			x = hi + y;

   514 			yr = x - hi;

   515 			if (y == yr)

   516 				hi = x;

   517 		}

   518 	}

   519 	sum = PyFloat_FromDouble(hi);

   520 

   521 _fsum_error:

   522 	PyFPE_END_PROTECT(hi)

   523 	Py_DECREF(iter);

   524 	if (p != ps)

   525 		PyMem_Free(p);

   526 	return sum;

   527 }

   528 

   529 #undef NUM_PARTIALS

   530 

   531 PyDoc_STRVAR(math_fsum_doc,

   532 "sum(iterable)\n\n\

   533 Return an accurate floating point sum of values in the iterable.\n\

   534 Assumes IEEE-754 floating point arithmetic.");

   535 

   536 static PyObject *

   537 math_factorial(PyObject *self, PyObject *arg)

   538 {

   539 	long i, x;

   540 	PyObject *result, *iobj, *newresult;

   541 

   542 	if (PyFloat_Check(arg)) {

   543 		double dx = PyFloat_AS_DOUBLE((PyFloatObject *)arg);

   544 		if (dx != floor(dx)) {

   545 			PyErr_SetString(PyExc_ValueError, 

   546 				"factorial() only accepts integral values");

   547 			return NULL;

   548 		}

   549 	}

   550 

   551 	x = PyInt_AsLong(arg);

   552 	if (x == -1 && PyErr_Occurred())

   553 		return NULL;

   554 	if (x < 0) {

   555 		PyErr_SetString(PyExc_ValueError, 

   556 			"factorial() not defined for negative values");

   557 		return NULL;

   558 	}

   559 

   560 	result = (PyObject *)PyInt_FromLong(1);

   561 	if (result == NULL)

   562 		return NULL;

   563 	for (i=1 ; i<=x ; i++) {

   564 		iobj = (PyObject *)PyInt_FromLong(i);

   565 		if (iobj == NULL)

   566 			goto error;

   567 		newresult = PyNumber_Multiply(result, iobj);

   568 		Py_DECREF(iobj);

   569 		if (newresult == NULL)

   570 			goto error;

   571 		Py_DECREF(result);

   572 		result = newresult;

   573 	}

   574 	return result;

   575 

   576 error:

   577 	Py_DECREF(result);

   578 	return NULL;

   579 }

   580 

   581 PyDoc_STRVAR(math_factorial_doc, "Return n!");

   582 

   583 static PyObject *

   584 math_trunc(PyObject *self, PyObject *number)

   585 {

   586 	return PyObject_CallMethod(number, "__trunc__", NULL);

   587 }

   588 

   589 PyDoc_STRVAR(math_trunc_doc,

   590 "trunc(x:Real) -> Integral\n"

   591 "\n"

   592 "Truncates x to the nearest Integral toward 0. Uses the __trunc__ magic method.");

   593 

   594 static PyObject *

   595 math_frexp(PyObject *self, PyObject *arg)

   596 {

   597 	int i;

   598 	double x = PyFloat_AsDouble(arg);

   599 	if (x == -1.0 && PyErr_Occurred())

   600 		return NULL;

   601 	/* deal with special cases directly, to sidestep platform

   602 	   differences */

   603 	if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {

   604 		i = 0;

   605 	}

   606 	else {

   607 		PyFPE_START_PROTECT("in math_frexp", return 0);

   608 		x = frexp(x, &i);

   609 		PyFPE_END_PROTECT(x);

   610 	}

   611 	return Py_BuildValue("(di)", x, i);

   612 }

   613 

   614 PyDoc_STRVAR(math_frexp_doc,

   615 "frexp(x)\n"

   616 "\n"

   617 "Return the mantissa and exponent of x, as pair (m, e).\n"

   618 "m is a float and e is an int, such that x = m * 2.**e.\n"

   619 "If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.");

   620 

   621 static PyObject *

   622 math_ldexp(PyObject *self, PyObject *args)

   623 {

   624 	double x, r;

   625 	PyObject *oexp;

   626 	long exp;

   627 	if (! PyArg_ParseTuple(args, "dO:ldexp", &x, &oexp))

   628 		return NULL;

   629 

   630 	if (PyLong_Check(oexp)) {

   631 		/* on overflow, replace exponent with either LONG_MAX

   632 		   or LONG_MIN, depending on the sign. */

   633 		exp = PyLong_AsLong(oexp);

   634 		if (exp == -1 && PyErr_Occurred()) {

   635 			if (PyErr_ExceptionMatches(PyExc_OverflowError)) {

   636 				if (Py_SIZE(oexp) < 0) {

   637 					exp = LONG_MIN;

   638 				}

   639 				else {

   640 					exp = LONG_MAX;

   641 				}

   642 				PyErr_Clear();

   643 			}

   644 			else {

   645 				/* propagate any unexpected exception */

   646 				return NULL;

   647 			}

   648 		}

   649 	}

   650 	else if (PyInt_Check(oexp)) {

   651 		exp = PyInt_AS_LONG(oexp);

   652 	}

   653 	else {

   654 		PyErr_SetString(PyExc_TypeError,

   655 				"Expected an int or long as second argument "

   656 				"to ldexp.");

   657 		return NULL;

   658 	}

   659 

   660 	if (x == 0. || !Py_IS_FINITE(x)) {

   661 		/* NaNs, zeros and infinities are returned unchanged */

   662 		r = x;

   663 		errno = 0;

   664 	} else if (exp > INT_MAX) {

   665 		/* overflow */

   666 		r = copysign(Py_HUGE_VAL, x);

   667 		errno = ERANGE;

   668 	} else if (exp < INT_MIN) {

   669 		/* underflow to +-0 */

   670 		r = copysign(0., x);

   671 		errno = 0;

   672 	} else {

   673 		errno = 0;

   674 		PyFPE_START_PROTECT("in math_ldexp", return 0);

   675 		r = ldexp(x, (int)exp);

   676 		PyFPE_END_PROTECT(r);

   677 		if (Py_IS_INFINITY(r))

   678 			errno = ERANGE;

   679 	}

   680 

   681 	if (errno && is_error(r))

   682 		return NULL;

   683 	return PyFloat_FromDouble(r);

   684 }

   685 

   686 PyDoc_STRVAR(math_ldexp_doc,

   687 "ldexp(x, i) -> x * (2**i)");

   688 

   689 static PyObject *

   690 math_modf(PyObject *self, PyObject *arg)

   691 {

   692 	double y, x = PyFloat_AsDouble(arg);

   693 	if (x == -1.0 && PyErr_Occurred())

   694 		return NULL;

   695 	/* some platforms don't do the right thing for NaNs and

   696 	   infinities, so we take care of special cases directly. */

   697 	if (!Py_IS_FINITE(x)) {

   698 		if (Py_IS_INFINITY(x))

   699 			return Py_BuildValue("(dd)", copysign(0., x), x);

   700 		else if (Py_IS_NAN(x))

   701 			return Py_BuildValue("(dd)", x, x);

   702 	}          

   703 

   704 	errno = 0;

   705 	PyFPE_START_PROTECT("in math_modf", return 0);

   706 	x = modf(x, &y);

   707 	PyFPE_END_PROTECT(x);

   708 	return Py_BuildValue("(dd)", x, y);

   709 }

   710 

   711 PyDoc_STRVAR(math_modf_doc,

   712 "modf(x)\n"

   713 "\n"

   714 "Return the fractional and integer parts of x.  Both results carry the sign\n"

   715 "of x.  The integer part is returned as a real.");

   716 

   717 /* A decent logarithm is easy to compute even for huge longs, but libm can't

   718    do that by itself -- loghelper can.  func is log or log10, and name is

   719    "log" or "log10".  Note that overflow isn't possible:  a long can contain

   720    no more than INT_MAX * SHIFT bits, so has value certainly less than

   721    2**(2**64 * 2**16) == 2**2**80, and log2 of that is 2**80, which is

   722    small enough to fit in an IEEE single.  log and log10 are even smaller.

   723 */

   724 

   725 static PyObject*

   726 loghelper(PyObject* arg, double (*func)(double), char *funcname)

   727 {

   728 	/* If it is long, do it ourselves. */

   729 	if (PyLong_Check(arg)) {

   730 		double x;

   731 		int e;

   732 		x = _PyLong_AsScaledDouble(arg, &e);

   733 		if (x <= 0.0) {

   734 			PyErr_SetString(PyExc_ValueError,

   735 					"math domain error");

   736 			return NULL;

   737 		}

   738 		/* Value is ~= x * 2**(e*PyLong_SHIFT), so the log ~=

   739 		   log(x) + log(2) * e * PyLong_SHIFT.

   740 		   CAUTION:  e*PyLong_SHIFT may overflow using int arithmetic,

   741 		   so force use of double. */

   742 		x = func(x) + (e * (double)PyLong_SHIFT) * func(2.0);

   743 		return PyFloat_FromDouble(x);

   744 	}

   745 

   746 	/* Else let libm handle it by itself. */

   747 	return math_1(arg, func, 0);

   748 }

   749 

   750 static PyObject *

   751 math_log(PyObject *self, PyObject *args)

   752 {

   753 	PyObject *arg;

   754 	PyObject *base = NULL;

   755 	PyObject *num, *den;

   756 	PyObject *ans;

   757 

   758 	if (!PyArg_UnpackTuple(args, "log", 1, 2, &arg, &base))

   759 		return NULL;

   760 

   761 	num = loghelper(arg, log, "log");

   762 	if (num == NULL || base == NULL)

   763 		return num;

   764 

   765 	den = loghelper(base, log, "log");

   766 	if (den == NULL) {

   767 		Py_DECREF(num);

   768 		return NULL;

   769 	}

   770 

   771 	ans = PyNumber_Divide(num, den);

   772 	Py_DECREF(num);

   773 	Py_DECREF(den);

   774 	return ans;

   775 }

   776 

   777 PyDoc_STRVAR(math_log_doc,

   778 "log(x[, base]) -> the logarithm of x to the given base.\n\

   779 If the base not specified, returns the natural logarithm (base e) of x.");

   780 

   781 static PyObject *

   782 math_log10(PyObject *self, PyObject *arg)

   783 {

   784 	return loghelper(arg, log10, "log10");

   785 }

   786 

   787 PyDoc_STRVAR(math_log10_doc,

   788 "log10(x) -> the base 10 logarithm of x.");

   789 

   790 static PyObject *

   791 math_fmod(PyObject *self, PyObject *args)

   792 {

   793 	PyObject *ox, *oy;

   794 	double r, x, y;

   795 	if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))

   796 		return NULL;

   797 	x = PyFloat_AsDouble(ox);

   798 	y = PyFloat_AsDouble(oy);

   799 	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())

   800 		return NULL;

   801 	/* fmod(x, +/-Inf) returns x for finite x. */

   802 	if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))

   803 		return PyFloat_FromDouble(x);

   804 	errno = 0;

   805 	PyFPE_START_PROTECT("in math_fmod", return 0);

   806 	r = fmod(x, y);

   807 	PyFPE_END_PROTECT(r);

   808 	if (Py_IS_NAN(r)) {

   809 		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))

   810 			errno = EDOM;

   811 		else

   812 			errno = 0;

   813 	}

   814 	if (errno && is_error(r))

   815 		return NULL;

   816 	else

   817 		return PyFloat_FromDouble(r);

   818 }

   819 

   820 PyDoc_STRVAR(math_fmod_doc,

   821 "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."

   822 "  x % y may differ.");

   823 

   824 static PyObject *

   825 math_hypot(PyObject *self, PyObject *args)

   826 {

   827 	PyObject *ox, *oy;

   828 	double r, x, y;

   829 	if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))

   830 		return NULL;

   831 	x = PyFloat_AsDouble(ox);

   832 	y = PyFloat_AsDouble(oy);

   833 	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())

   834 		return NULL;

   835 	/* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */

   836 	if (Py_IS_INFINITY(x))

   837 		return PyFloat_FromDouble(fabs(x));

   838 	if (Py_IS_INFINITY(y))

   839 		return PyFloat_FromDouble(fabs(y));

   840 	errno = 0;

   841 	PyFPE_START_PROTECT("in math_hypot", return 0);

   842 	r = hypot(x, y);

   843 	PyFPE_END_PROTECT(r);

   844 	if (Py_IS_NAN(r)) {

   845 		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))

   846 			errno = EDOM;

   847 		else

   848 			errno = 0;

   849 	}

   850 	else if (Py_IS_INFINITY(r)) {

   851 		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))

   852 			errno = ERANGE;

   853 		else

   854 			errno = 0;

   855 	}

   856 	if (errno && is_error(r))

   857 		return NULL;

   858 	else

   859 		return PyFloat_FromDouble(r);

   860 }

   861 

   862 PyDoc_STRVAR(math_hypot_doc,

   863 "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");

   864 

   865 /* pow can't use math_2, but needs its own wrapper: the problem is

   866    that an infinite result can arise either as a result of overflow

   867    (in which case OverflowError should be raised) or as a result of

   868    e.g. 0.**-5. (for which ValueError needs to be raised.)

   869 */

   870 

   871 static PyObject *

   872 math_pow(PyObject *self, PyObject *args)

   873 {

   874 	PyObject *ox, *oy;

   875 	double r, x, y;

   876 	int odd_y;

   877 

   878 	if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))

   879 		return NULL;

   880 	x = PyFloat_AsDouble(ox);

   881 	y = PyFloat_AsDouble(oy);

   882 	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())

   883 		return NULL;

   884 

   885 	/* deal directly with IEEE specials, to cope with problems on various

   886 	   platforms whose semantics don't exactly match C99 */

   887 	r = 0.; /* silence compiler warning */

   888 	if (!Py_IS_FINITE(x) || !Py_IS_FINITE(y)) {

   889 		errno = 0;

   890 		if (Py_IS_NAN(x))

   891 			r = y == 0. ? 1. : x; /* NaN**0 = 1 */

   892 		else if (Py_IS_NAN(y))

   893 			r = x == 1. ? 1. : y; /* 1**NaN = 1 */

   894 		else if (Py_IS_INFINITY(x)) {

   895 			odd_y = Py_IS_FINITE(y) && fmod(fabs(y), 2.0) == 1.0;

   896 			if (y > 0.)

   897 				r = odd_y ? x : fabs(x);

   898 			else if (y == 0.)

   899 				r = 1.;

   900 			else /* y < 0. */

   901 				r = odd_y ? copysign(0., x) : 0.;

   902 		}

   903 		else if (Py_IS_INFINITY(y)) {

   904 			if (fabs(x) == 1.0)

   905 				r = 1.;

   906 			else if (y > 0. && fabs(x) > 1.0)

   907 				r = y;

   908 			else if (y < 0. && fabs(x) < 1.0) {

   909 				r = -y; /* result is +inf */

   910 				if (x == 0.) /* 0**-inf: divide-by-zero */

   911 					errno = EDOM;

   912 			}

   913 			else

   914 				r = 0.;

   915 		}

   916 	}

   917 	else {

   918 		/* let libm handle finite**finite */

   919 		errno = 0;

   920 		PyFPE_START_PROTECT("in math_pow", return 0);

   921 		r = pow(x, y);

   922 		PyFPE_END_PROTECT(r);

   923 		/* a NaN result should arise only from (-ve)**(finite

   924 		   non-integer); in this case we want to raise ValueError. */

   925 		if (!Py_IS_FINITE(r)) {

   926 			if (Py_IS_NAN(r)) {

   927 				errno = EDOM;

   928 			}

   929 			/* 

   930 			   an infinite result here arises either from:

   931 			   (A) (+/-0.)**negative (-> divide-by-zero)

   932 			   (B) overflow of x**y with x and y finite

   933 			*/

   934 			else if (Py_IS_INFINITY(r)) {

   935 				if (x == 0.)

   936 					errno = EDOM;

   937 				else

   938 					errno = ERANGE;

   939 			}

   940 		}

   941 	}

   942 

   943 	if (errno && is_error(r))

   944 		return NULL;

   945 	else

   946 		return PyFloat_FromDouble(r);

   947 }

   948 

   949 PyDoc_STRVAR(math_pow_doc,

   950 "pow(x,y)\n\nReturn x**y (x to the power of y).");

   951 

   952 static const double degToRad = Py_MATH_PI / 180.0;

   953 static const double radToDeg = 180.0 / Py_MATH_PI;

   954 

   955 static PyObject *

   956 math_degrees(PyObject *self, PyObject *arg)

   957 {

   958 	double x = PyFloat_AsDouble(arg);

   959 	if (x == -1.0 && PyErr_Occurred())

   960 		return NULL;

   961 	return PyFloat_FromDouble(x * radToDeg);

   962 }

   963 

   964 PyDoc_STRVAR(math_degrees_doc,

   965 "degrees(x) -> converts angle x from radians to degrees");

   966 

   967 static PyObject *

   968 math_radians(PyObject *self, PyObject *arg)

   969 {

   970 	double x = PyFloat_AsDouble(arg);

   971 	if (x == -1.0 && PyErr_Occurred())

   972 		return NULL;

   973 	return PyFloat_FromDouble(x * degToRad);

   974 }

   975 

   976 PyDoc_STRVAR(math_radians_doc,

   977 "radians(x) -> converts angle x from degrees to radians");

   978 

   979 static PyObject *

   980 math_isnan(PyObject *self, PyObject *arg)

   981 {

   982 	double x = PyFloat_AsDouble(arg);

   983 	if (x == -1.0 && PyErr_Occurred())

   984 		return NULL;

   985 	return PyBool_FromLong((long)Py_IS_NAN(x));

   986 }

   987 

   988 PyDoc_STRVAR(math_isnan_doc,

   989 "isnan(x) -> bool\n\

   990 Checks if float x is not a number (NaN)");

   991 

   992 static PyObject *

   993 math_isinf(PyObject *self, PyObject *arg)

   994 {

   995 	double x = PyFloat_AsDouble(arg);

   996 	if (x == -1.0 && PyErr_Occurred())

   997 		return NULL;

   998 	return PyBool_FromLong((long)Py_IS_INFINITY(x));

   999 }

  1000 

  1001 PyDoc_STRVAR(math_isinf_doc,

  1002 "isinf(x) -> bool\n\

  1003 Checks if float x is infinite (positive or negative)");

  1004 

  1005 static PyMethodDef math_methods[] = {

  1006 	{"acos",	math_acos,	METH_O,		math_acos_doc},

  1007 	{"acosh",	math_acosh,	METH_O,		math_acosh_doc},

  1008 	{"asin",	math_asin,	METH_O,		math_asin_doc},

  1009 	{"asinh",	math_asinh,	METH_O,		math_asinh_doc},

  1010 	{"atan",	math_atan,	METH_O,		math_atan_doc},

  1011 	{"atan2",	math_atan2,	METH_VARARGS,	math_atan2_doc},

  1012 	{"atanh",	math_atanh,	METH_O,		math_atanh_doc},

  1013 	{"ceil",	math_ceil,	METH_O,		math_ceil_doc},

  1014 	{"copysign",	math_copysign,	METH_VARARGS,	math_copysign_doc},

  1015 	{"cos",		math_cos,	METH_O,		math_cos_doc},

  1016 	{"cosh",	math_cosh,	METH_O,		math_cosh_doc},

  1017 	{"degrees",	math_degrees,	METH_O,		math_degrees_doc},

  1018 	{"exp",		math_exp,	METH_O,		math_exp_doc},

  1019 	{"fabs",	math_fabs,	METH_O,		math_fabs_doc},

  1020 	{"factorial",	math_factorial,	METH_O,		math_factorial_doc},

  1021 	{"floor",	math_floor,	METH_O,		math_floor_doc},

  1022 	{"fmod",	math_fmod,	METH_VARARGS,	math_fmod_doc},

  1023 	{"frexp",	math_frexp,	METH_O,		math_frexp_doc},

  1024 	{"fsum",	math_fsum,	METH_O,		math_fsum_doc},

  1025 	{"hypot",	math_hypot,	METH_VARARGS,	math_hypot_doc},

  1026 	{"isinf",	math_isinf,	METH_O,		math_isinf_doc},

  1027 	{"isnan",	math_isnan,	METH_O,		math_isnan_doc},

  1028 	{"ldexp",	math_ldexp,	METH_VARARGS,	math_ldexp_doc},

  1029 	{"log",		math_log,	METH_VARARGS,	math_log_doc},

  1030 	{"log1p",	math_log1p,	METH_O,		math_log1p_doc},

  1031 	{"log10",	math_log10,	METH_O,		math_log10_doc},

  1032 	{"modf",	math_modf,	METH_O,		math_modf_doc},

  1033 	{"pow",		math_pow,	METH_VARARGS,	math_pow_doc},

  1034 	{"radians",	math_radians,	METH_O,		math_radians_doc},

  1035 	{"sin",		math_sin,	METH_O,		math_sin_doc},

  1036 	{"sinh",	math_sinh,	METH_O,		math_sinh_doc},

  1037 	{"sqrt",	math_sqrt,	METH_O,		math_sqrt_doc},

  1038 	{"tan",		math_tan,	METH_O,		math_tan_doc},

  1039 	{"tanh",	math_tanh,	METH_O,		math_tanh_doc},

  1040 	{"trunc",	math_trunc,	METH_O,		math_trunc_doc},

  1041 	{NULL,		NULL}		/* sentinel */

  1042 };

  1043 

  1044 

  1045 PyDoc_STRVAR(module_doc,

  1046 "This module is always available.  It provides access to the\n"

  1047 "mathematical functions defined by the C standard.");

  1048 

  1049 PyMODINIT_FUNC

  1050 initmath(void)

  1051 {

  1052 	PyObject *m;

  1053 

  1054 	m = Py_InitModule3("math", math_methods, module_doc);

  1055 	if (m == NULL)

  1056 		goto finally;

  1057 

  1058 	PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));

  1059 	PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));

  1060 

  1061     finally:

  1062 	return;

  1063 }