#ifndef Py_BUILD_CORE_BUILTIN

# define Py_BUILD_CORE_MODULE 1

#include "Python.h"

#include "pycore_abstract.h" // _PyNumber_Index()

#include "pycore_bitutils.h" // _Py_bit_length()

#include "pycore_call.h" // _PyObject_CallNoArgs()

#include "pycore_import.h" // _PyImport_SetModuleString()

#include "pycore_long.h" // _PyLong_GetZero()

#include "pycore_moduleobject.h" // _PyModule_GetState()

#include "pycore_object.h" // _PyObject_LookupSpecial()

#include "pycore_pymath.h" // _PY_SHORT_FLOAT_REPR

#include <float.h>

#include "_math.h"

#include <stdbool.h>

#include "clinic/mathmodule.c.h"

typedef struct{ double hi; double lo; } DoubleLength;

static DoubleLength

dl_fast_sum(double a, double b)

{

/* Algorithm 1.1. Compensated summation of two floating-point numbers. */

assert(fabs(a) >= fabs(b));

double x = a + b;

double y = (a - x) + b;

return (DoubleLength) {x, y};

}

static DoubleLength

dl_sum(double a, double b)

{

/* Algorithm 3.1 Error-free transformation of the sum */

double x = a + b;

double z = x - a;

double y = (a - (x - z)) + (b - z);

return (DoubleLength) {x, y};

}

#ifndef UNRELIABLE_FMA

static DoubleLength

dl_mul(double x, double y)

{

/* Algorithm 3.5. Error-free transformation of a product */

double z = x * y;

double zz = fma(x, y, -z);

return (DoubleLength) {z, zz};

}

static DoubleLength

dl_split(double x) {

// Dekker (5.5) and (5.6).

double t = x * 134217729.0; // Veltkamp constant = 2.0 ** 27 + 1

double hi = t - (t - x);

double lo = x - hi;

return (DoubleLength) {hi, lo};

}

static DoubleLength

dl_mul(double x, double y)

{

// Dekker (5.12) and mul12()

DoubleLength xx = dl_split(x);

DoubleLength yy = dl_split(y);

double p = xx.hi * yy.hi;

double q = xx.hi * yy.lo + xx.lo * yy.hi;

double z = p + q;

double zz = p - z + q + xx.lo * yy.lo;

return (DoubleLength) {z, zz};

}

typedef struct { double hi; double lo; double tiny; } TripleLength;

static const TripleLength tl_zero = {0.0, 0.0, 0.0};

static TripleLength

tl_fma(double x, double y, TripleLength total)

{

/* Algorithm 5.10 with SumKVert for K=3 */

DoubleLength pr = dl_mul(x, y);

DoubleLength sm = dl_sum(total.hi, pr.hi);

DoubleLength r1 = dl_sum(total.lo, pr.lo);

DoubleLength r2 = dl_sum(r1.hi, sm.lo);

return (TripleLength) {sm.hi, r2.hi, total.tiny + r1.lo + r2.lo};

}

static double

tl_to_d(TripleLength total)

{

DoubleLength last = dl_sum(total.lo, total.hi);

return total.tiny + last.lo + last.hi;

}

static const double pi = 3.141592653589793238462643383279502884197;

static const double logpi = 1.144729885849400174143427351353058711647;

#define ASSIGN_DOUBLE(target_var, obj, error_label) \

if (PyFloat_CheckExact(obj)) { \

target_var = PyFloat_AS_DOUBLE(obj); \

} \

else if (PyLong_CheckExact(obj)) { \

target_var = PyLong_AsDouble(obj); \

if (target_var == -1.0 && PyErr_Occurred()) { \

goto error_label; \

} \

} \

else { \

target_var = PyFloat_AsDouble(obj); \

if (target_var == -1.0 && PyErr_Occurred()) { \

goto error_label; \

} \

}

static double

m_sinpi(double x)

{

double y, r;

int n;

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

assert(isfinite(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:

Py_UNREACHABLE();

}

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 (!isfinite(x)) {

if (isnan(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(INFINITY, 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 (isinf(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/m_sinpi(x);

}

else {

errno = ERANGE;

return INFINITY;

}

}

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 / m_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 (isinf(r))

errno = ERANGE;

return r;

}

static double

m_lgamma(double x)

{

double r;

double absx;

/* special cases */

if (!isfinite(x)) {

if (isnan(x))

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

else

return INFINITY; /* 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 INFINITY; /* 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(m_sinpi(absx))) - log(absx) - r;

if (isinf(r))

errno = ERANGE;

return r;

}

static double

m_remainder(double x, double y)

{

/* Deal with most common case first. */

if (isfinite(x) && isfinite(y)) {

double absx, absy, c, m, r;

if (y == 0.0) {

return Py_NAN;

}

absx = fabs(x);

absy = fabs(y);

m = fmod(absx, absy);

/*

c = absy - m;

if (m < c) {

r = m;

}

else if (m > c) {

r = -c;

}

else {

/*

assert(m == c);

r = m - 2.0 * fmod(0.5 * (absx - m), absy);

}

return copysign(1.0, x) * r;

}

/* Special values. */

if (isnan(x)) {

return x;

}

if (isnan(y)) {

return y;

}

if (isinf(x)) {

return Py_NAN;

}

assert(isinf(y));

return x;

}

static double

m_log(double x)

{

if (isfinite(x)) {

if (x > 0.0)

return log(x);

errno = EDOM;

if (x == 0.0)

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

else

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

}

else if (isnan(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 (!isfinite(x)) {

if (isnan(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) {

return log2(x);

}

else if (x == 0.0) {

errno = EDOM;

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

}

else {

errno = EDOM;

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

}

}

static double

m_log10(double x)

{

if (isfinite(x)) {

if (x > 0.0)

return log10(x);

errno = EDOM;

if (x == 0.0)

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

else

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

}

else if (isnan(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 raise_edom)

{

int result = 1; /* presumption of guilt */

assert(errno); /* non-zero errno is a precondition for calling */

if (errno == EDOM) {

if (raise_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.5)

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,

const char *err_msg)

{

double x, r;

x = PyFloat_AsDouble(arg);

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

return NULL;

errno = 0;

r = (*func)(x);

if (isnan(r) && !isnan(x))

goto domain_err; /* domain error */

if (isinf(r) && isfinite(x)) {

if (can_overflow)

PyErr_SetString(PyExc_OverflowError,

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

else

goto domain_err; /* singularity */

return NULL;

}

if (isfinite(r) && errno && is_error(r, 1))

/* this branch unnecessary on most platforms */

return NULL;

return PyFloat_FromDouble(r);

domain_err:

if (err_msg) {

char *buf = PyOS_double_to_string(x, 'r', 0, Py_DTSF_ADD_DOT_0, NULL);

if (buf) {

PyErr_Format(PyExc_ValueError, err_msg, buf);

PyMem_Free(buf);

}

}

else {

PyErr_SetString(PyExc_ValueError, "math domain error");

}

return NULL;

}

static PyObject *

math_1a(PyObject *arg, double (*func) (double), const char *err_msg)

{

double x, r;

x = PyFloat_AsDouble(arg);

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

return NULL;

errno = 0;

r = (*func)(x);

if (errno && is_error(r, err_msg ? 0 : 1)) {

if (err_msg && errno == EDOM) {

assert(!PyErr_Occurred()); /* exception is not set by is_error() */

char *buf = PyOS_double_to_string(x, 'r', 0, Py_DTSF_ADD_DOT_0, NULL);

if (buf) {

PyErr_Format(PyExc_ValueError, err_msg, buf);

PyMem_Free(buf);

}

}

return NULL;

}

return PyFloat_FromDouble(r);

}

static PyObject *

math_2(PyObject *const *args, Py_ssize_t nargs,

double (*func) (double, double), const char *funcname)

{

double x, y, r;

if (!_PyArg_CheckPositional(funcname, nargs, 2, 2))

return NULL;

x = PyFloat_AsDouble(args[0]);

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

return NULL;

}

y = PyFloat_AsDouble(args[1]);

if (y == -1.0 && PyErr_Occurred()) {

return NULL;

}

errno = 0;

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

if (isnan(r)) {

if (!isnan(x) && !isnan(y))

errno = EDOM;

else

errno = 0;

}

else if (isinf(r)) {

if (isfinite(x) && isfinite(y))

errno = ERANGE;

else

errno = 0;

}

if (errno && is_error(r, 1))

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, NULL); \

}\

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

#define FUNC1D(funcname, func, can_overflow, docstring, err_msg) \

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

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

}\

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

#define FUNC1A(funcname, func, docstring) \

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

return math_1a(args, func, NULL); \

}\

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

#define FUNC1AD(funcname, func, docstring, err_msg) \

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

return math_1a(args, func, err_msg); \

}\

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

#define FUNC2(funcname, func, docstring) \

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

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

}\

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

FUNC1D(acos, acos, 0,

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

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

"The result is between 0 and pi.",

"expected a number in range from -1 up to 1, got %s")

FUNC1D(acosh, acosh, 0,

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

"Return the inverse hyperbolic cosine of x.",

"expected argument value not less than 1, got %s")

FUNC1D(asin, asin, 0,

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

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

"The result is between -pi/2 and pi/2.",

"expected a number in range from -1 up to 1, got %s")

FUNC1(asinh, 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.\n\n"

"The result is between -pi/2 and pi/2.")

FUNC2(atan2, 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.")

FUNC1D(atanh, atanh, 0,

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

"Return the inverse hyperbolic tangent of x.",

"expected a number between -1 and 1, got %s")

FUNC1(cbrt, cbrt, 0,

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

"Return the cube root of x.")

static PyObject *

math_ceil(PyObject *module, PyObject *number)

{

double x;

if (PyFloat_CheckExact(number)) {

x = PyFloat_AS_DOUBLE(number);

}

else {

PyObject *result = _PyObject_MaybeCallSpecialNoArgs(number, &_Py_ID(__ceil__));