#include "Python.h"

#include "longintrepr.h" /* just for SHIFT */

#ifdef _OSF_SOURCE

extern double copysign(double, double);

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

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;

}

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);

}

static double

m_log(double x)

{

if (Py_IS_FINITE(x)) {

if (x > 0.0)

return log(x);

errno = EDOM;

if (x == 0.0)

return -Py_HUGE_VAL; /* log(0) = -inf */

else

return Py_NAN; /* log(-ve) = nan */

}

else if (Py_IS_NAN(x))

return x; /* log(nan) = nan */

else if (x > 0.0)

return x; /* log(inf) = inf */

else {

errno = EDOM;

return Py_NAN; /* log(-inf) = nan */

}

}

static double

m_log10(double x)

{

if (Py_IS_FINITE(x)) {

if (x > 0.0)

return log10(x);

errno = EDOM;

if (x == 0.0)

return -Py_HUGE_VAL; /* log10(0) = -inf */

else

return Py_NAN; /* log10(-ve) = nan */

}

else if (Py_IS_NAN(x))

return x; /* log10(nan) = nan */

else if (x > 0.0)

return x; /* log10(inf) = inf */

else {

errno = EDOM;

return Py_NAN; /* log10(-inf) = nan */

}

}

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);

}

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.")

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

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;

}

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

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

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.

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,

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,

"Find x!. Raise a ValueError if x is negative or non-integral.");

static PyObject *

math_trunc(PyObject *self, PyObject *number)

{

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

}

PyDoc_STRVAR(math_trunc_doc,

"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

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,

"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

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

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,

"of x and are floats.");

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 ~=

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, m_log, "log");

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

return num;

den = loghelper(base, m_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,

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

static PyObject *

math_log10(PyObject *self, PyObject *arg)

{

return loghelper(arg, m_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,

" 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).");

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

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

if (!Py_IS_FINITE(r)) {

if (Py_IS_NAN(r)) {

errno = EDOM;

}

/*

else if (Py_IS_INFINITY(r)) {

if (x == 0.)

errno = EDOM;

else

errno = ERANGE;

}

}

}

if (errno && is_error(r))

return NULL;

else