#include "Python.h"

#include "_math.h"

#ifdef _OSF_SOURCE

extern double copysign(double, double);

static const double pi = 3.141592653589793238462643383279502884197;

static const double sqrtpi = 1.772453850905516027298167483341145182798;

static double

sinpi(double x)

{

double y, r;

int n;

/* this function should only ever be called for finite arguments */

assert(Py_IS_FINITE(x));

y = fmod(fabs(x), 2.0);

n = (int)round(2.0*y);

assert(0 <= n && n <= 4);

switch (n) {

case 0:

r = sin(pi*y);

break;

case 1:

r = cos(pi*(y-0.5));

break;

case 2:

/* N.B. -sin(pi*(y-1.0)) is *not* equivalent: it would give

r = sin(pi*(1.0-y));

break;

case 3:

r = -cos(pi*(y-1.5));

break;

case 4:

r = sin(pi*(y-2.0));

break;

default:

assert(0); /* should never get here */

r = -1.23e200; /* silence gcc warning */

}

return copysign(1.0, x)*r;

}

#define LANCZOS_N 13

static const double lanczos_g = 6.024680040776729583740234375;

static const double lanczos_g_minus_half = 5.524680040776729583740234375;

static const double lanczos_num_coeffs[LANCZOS_N] = {

23531376880.410759688572007674451636754734846804940,

42919803642.649098768957899047001988850926355848959,

35711959237.355668049440185451547166705960488635843,

17921034426.037209699919755754458931112671403265390,

6039542586.3520280050642916443072979210699388420708,

1439720407.3117216736632230727949123939715485786772,

248874557.86205415651146038641322942321632125127801,

31426415.585400194380614231628318205362874684987640,

2876370.6289353724412254090516208496135991145378768,

186056.26539522349504029498971604569928220784236328,

8071.6720023658162106380029022722506138218516325024,

210.82427775157934587250973392071336271166969580291,

2.5066282746310002701649081771338373386264310793408

};

static const double lanczos_den_coeffs[LANCZOS_N] = {

0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,

13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0};

#define NGAMMA_INTEGRAL 23

static const double gamma_integral[NGAMMA_INTEGRAL] = {

1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,

3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,

1307674368000.0, 20922789888000.0, 355687428096000.0,

6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,

51090942171709440000.0, 1124000727777607680000.0,

};

static double

lanczos_sum(double x)

{

double num = 0.0, den = 0.0;

int i;

assert(x > 0.0);

/* evaluate the rational function lanczos_sum(x). For large

if (x < 5.0) {

for (i = LANCZOS_N; --i >= 0; ) {

num = num * x + lanczos_num_coeffs[i];

den = den * x + lanczos_den_coeffs[i];

}

}

else {

for (i = 0; i < LANCZOS_N; i++) {

num = num / x + lanczos_num_coeffs[i];

den = den / x + lanczos_den_coeffs[i];

}

}

return num/den;

}

static double

m_tgamma(double x)

{

double absx, r, y, z, sqrtpow;

/* special cases */

if (!Py_IS_FINITE(x)) {

if (Py_IS_NAN(x) || x > 0.0)

return x; /* tgamma(nan) = nan, tgamma(inf) = inf */

else {

errno = EDOM;

return Py_NAN; /* tgamma(-inf) = nan, invalid */

}

}

if (x == 0.0) {

errno = EDOM;

return 1.0/x; /* tgamma(+-0.0) = +-inf, divide-by-zero */

}

/* integer arguments */

if (x == floor(x)) {

if (x < 0.0) {

errno = EDOM; /* tgamma(n) = nan, invalid for */

return Py_NAN; /* negative integers n */

}

if (x <= NGAMMA_INTEGRAL)

return gamma_integral[(int)x - 1];

}

absx = fabs(x);

/* tiny arguments: tgamma(x) ~ 1/x for x near 0 */

if (absx < 1e-20) {

r = 1.0/x;

if (Py_IS_INFINITY(r))

errno = ERANGE;

return r;

}

/* large arguments: assuming IEEE 754 doubles, tgamma(x) overflows for

if (absx > 200.0) {

if (x < 0.0) {

return 0.0/sinpi(x);

}

else {

errno = ERANGE;

return Py_HUGE_VAL;

}

}

y = absx + lanczos_g_minus_half;

/* compute error in sum */

if (absx > lanczos_g_minus_half) {

/* note: the correction can be foiled by an optimizing

double q = y - absx;

z = q - lanczos_g_minus_half;

}

else {

double q = y - lanczos_g_minus_half;

z = q - absx;

}

z = z * lanczos_g / y;

if (x < 0.0) {

r = -pi / sinpi(absx) / absx * exp(y) / lanczos_sum(absx);

r -= z * r;

if (absx < 140.0) {

r /= pow(y, absx - 0.5);

}

else {

sqrtpow = pow(y, absx / 2.0 - 0.25);

r /= sqrtpow;

r /= sqrtpow;

}

}

else {

r = lanczos_sum(absx) / exp(y);

r += z * r;

if (absx < 140.0) {

r *= pow(y, absx - 0.5);

}

else {

sqrtpow = pow(y, absx / 2.0 - 0.25);

r *= sqrtpow;

r *= sqrtpow;

}

}

if (Py_IS_INFINITY(r))

errno = ERANGE;

return r;

}

static double

m_lgamma(double x)

{

double r, absx;

/* special cases */

if (!Py_IS_FINITE(x)) {

if (Py_IS_NAN(x))

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

else

return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */

}

/* integer arguments */

if (x == floor(x) && x <= 2.0) {

if (x <= 0.0) {

errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */

return Py_HUGE_VAL; /* integers n <= 0 */

}

else {

return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */

}

}

absx = fabs(x);

/* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */

if (absx < 1e-20)

return -log(absx);

/* Lanczos' formula */

if (x > 0.0) {

/* we could save a fraction of a ulp in accuracy by having a

r = log(lanczos_sum(x)) - lanczos_g +

(x-0.5)*(log(x+lanczos_g-0.5)-1);

}

else {

r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -

(log(lanczos_sum(absx)) - lanczos_g +

(absx-0.5)*(log(absx+lanczos_g-0.5)-1));

}

if (Py_IS_INFINITY(r))

errno = ERANGE;

return r;

}

#define ERF_SERIES_CUTOFF 1.5

#define ERF_SERIES_TERMS 25

#define ERFC_CONTFRAC_CUTOFF 30.0

#define ERFC_CONTFRAC_TERMS 50

static double

m_erf_series(double x)

{

double x2, acc, fk, result;

int i, saved_errno;

x2 = x * x;

acc = 0.0;

fk = (double)ERF_SERIES_TERMS + 0.5;

for (i = 0; i < ERF_SERIES_TERMS; i++) {

acc = 2.0 + x2 * acc / fk;

fk -= 1.0;

}

/* Make sure the exp call doesn't affect errno;

saved_errno = errno;

result = acc * x * exp(-x2) / sqrtpi;

errno = saved_errno;

return result;

}

static double

m_erfc_contfrac(double x)

{

double x2, a, da, p, p_last, q, q_last, b, result;

int i, saved_errno;

if (x >= ERFC_CONTFRAC_CUTOFF)

return 0.0;

x2 = x*x;

a = 0.0;

da = 0.5;

p = 1.0; p_last = 0.0;

q = da + x2; q_last = 1.0;

for (i = 0; i < ERFC_CONTFRAC_TERMS; i++) {

double temp;

a += da;

da += 2.0;

b = da + x2;

temp = p; p = b*p - a*p_last; p_last = temp;

temp = q; q = b*q - a*q_last; q_last = temp;

}

/* Issue #8986: On some platforms, exp sets errno on underflow to zero;

saved_errno = errno;

result = p / q * x * exp(-x2) / sqrtpi;

errno = saved_errno;

return result;

}

static double

m_erf(double x)

{

double absx, cf;

if (Py_IS_NAN(x))

return x;

absx = fabs(x);

if (absx < ERF_SERIES_CUTOFF)

return m_erf_series(x);

else {

cf = m_erfc_contfrac(absx);

return x > 0.0 ? 1.0 - cf : cf - 1.0;

}

}

static double

m_erfc(double x)

{

double absx, cf;

if (Py_IS_NAN(x))

return x;

absx = fabs(x);

if (absx < ERF_SERIES_CUTOFF)

return 1.0 - m_erf_series(x);

else {

cf = m_erfc_contfrac(absx);

return x > 0.0 ? cf : 2.0 - cf;

}

}

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 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 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_1a(PyObject *arg, double (*func) (double))

{

double x, r;

x = PyFloat_AsDouble(arg);

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

return NULL;

errno = 0;

PyFPE_START_PROTECT("in math_1a", return 0);

r = (*func)(x);

PyFPE_END_PROTECT(r);

if (errno && is_error(r))

return NULL;

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 FUNC1A(funcname, func, docstring) \

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

return math_1a(args, func); \

}\

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

FUNC1A(erf, m_erf,

"erf(x)\n\nError function at x.")

FUNC1A(erfc, m_erfc,

"erfc(x)\n\nComplementary error function at x.")

FUNC1(exp, exp, 1,

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

FUNC1(expm1, m_expm1, 1,

"expm1(x)\n\nReturn exp(x)-1.\n"

"This function avoids the loss of precision involved in the direct "

"evaluation of exp(x)-1 for small 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.")

FUNC1A(gamma, m_tgamma,

"gamma(x)\n\nGamma function at x.")

FUNC1A(lgamma, m_lgamma,

"lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")

FUNC1(log1p, m_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;