#include "Python.h"

#include "_math.h"

#include "clinic/mathmodule.c.h"

static const double pi = 3.141592653589793238462643383279502884197;

static const double sqrtpi = 1.772453850905516027298167483341145182798;

static const double logpi = 1.144729885849400174143427351353058711647;

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_inf(void)

{

#ifndef PY_NO_SHORT_FLOAT_REPR

return _Py_dg_infinity(0);

return Py_HUGE_VAL;

}

#if !defined(PY_NO_SHORT_FLOAT_REPR) || defined(Py_NAN)

static double

m_nan(void)

{

#ifndef PY_NO_SHORT_FLOAT_REPR

return _Py_dg_stdnan(0);

return Py_NAN;

}

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;

/* tgamma(+-0.0) = +-inf, divide-by-zero */

return copysign(Py_HUGE_VAL, x);

}

/* 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. We could save a fraction of a ulp in accuracy by

r = log(lanczos_sum(absx)) - lanczos_g;

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

if (x < 0.0)

/* Use reflection formula to get value for negative x. */

r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;

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_log2(double x)

{

if (!Py_IS_FINITE(x)) {

if (Py_IS_NAN(x))

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

else if (x > 0.0)

return x; /* log2(+inf) = +inf */

else {

errno = EDOM;

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

}

}

if (x > 0.0) {

#ifdef HAVE_LOG2

return log2(x);

double m;

int e;

m = frexp(x, &e);

/* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when

if (x >= 1.0) {

return log(2.0 * m) / log(2.0) + (e - 1);

}

else {

return log(m) / log(2.0) + e;

}

}

else if (x == 0.0) {

errno = EDOM;

return -Py_HUGE_VAL; /* log2(0) = -inf, divide-by-zero */

}

else {

errno = EDOM;

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

}

}

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_gcd_impl(PyObject *module, PyObject *a, PyObject *b)

{

PyObject *g;

a = PyNumber_Index(a);

if (a == NULL)

return NULL;

b = PyNumber_Index(b);

if (b == NULL) {

Py_DECREF(a);

return NULL;

}

g = _PyLong_GCD(a, b);

Py_DECREF(a);

Py_DECREF(b);

return g;

}

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

PyObject *(*from_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) && !Py_IS_NAN(x)) {

PyErr_SetString(PyExc_ValueError,

"math domain error"); /* invalid arg */

return NULL;

}

if (Py_IS_INFINITY(r) && Py_IS_FINITE(x)) {

if (can_overflow)

PyErr_SetString(PyExc_OverflowError,

"math range error"); /* overflow */

else

PyErr_SetString(PyExc_ValueError,

"math domain error"); /* singularity */

return NULL;

}

if (Py_IS_FINITE(r) && errno && is_error(r))

/* this branch unnecessary on most platforms */

return NULL;

return (*from_double_func)(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_1(PyObject *arg, double (*func) (double), int can_overflow)

{

return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);

}

static PyObject *

math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)

{

return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);

}

static PyObject *

math_2(PyObject *args, double (*func) (double, double), const 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($module, x, /)\n--\n\n"

"Return the arc cosine (measured in radians) of x.")

FUNC1(acosh, m_acosh, 0,

"acosh($module, x, /)\n--\n\n"

"Return the inverse hyperbolic cosine of x.")

FUNC1(asin, asin, 0,

"asin($module, x, /)\n--\n\n"

"Return the arc sine (measured in radians) of x.")

FUNC1(asinh, m_asinh, 0,

"asinh($module, x, /)\n--\n\n"

"Return the inverse hyperbolic sine of x.")

FUNC1(atan, atan, 0,

"atan($module, x, /)\n--\n\n"

"Return the arc tangent (measured in radians) of x.")

FUNC2(atan2, m_atan2,

"atan2($module, y, x, /)\n--\n\n"

"Return the arc tangent (measured in radians) of y/x.\n\n"

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

FUNC1(atanh, m_atanh, 0,

"atanh($module, x, /)\n--\n\n"

"Return the inverse hyperbolic tangent of x.")

static PyObject *

math_ceil(PyObject *module, PyObject *number)

{

_Py_IDENTIFIER(__ceil__);

PyObject *method, *result;

method = _PyObject_LookupSpecial(number, &PyId___ceil__);

if (method == NULL) {

if (PyErr_Occurred())

return NULL;

return math_1_to_int(number, ceil, 0);

}

result = _PyObject_CallNoArg(method);

Py_DECREF(method);

return result;

}

FUNC2(copysign, copysign,

"copysign($module, x, y, /)\n--\n\n"

"Return a float with the magnitude (absolute value) of x but the sign of y.\n\n"

"On platforms that support signed zeros, copysign(1.0, -0.0)\n"

"returns -1.0.\n")

FUNC1(cos, cos, 0,

"cos($module, x, /)\n--\n\n"

"Return the cosine of x (measured in radians).")

FUNC1(cosh, cosh, 1,