#ifndef Py_BUILD_CORE_BUILTIN

# define Py_BUILD_CORE_MODULE 1

#include "Python.h"

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

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

#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 {

PyObject *str___ceil__;

PyObject *str___floor__;

PyObject *str___trunc__;

} math_module_state;

static inline math_module_state*

get_math_module_state(PyObject *module)

{

void *state = _PyModule_GetState(module);

assert(state != NULL);

return (math_module_state *)state;

}

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

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 (!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_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 (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/m_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 / 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 (Py_IS_INFINITY(r))

errno = ERANGE;

return r;

}

static double

m_lgamma(double x)

{

double r;

double 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(m_sinpi(absx))) - log(absx) - r;

if (Py_IS_INFINITY(r))

errno = ERANGE;

return r;

}

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_remainder(double x, double y)

{

/* Deal with most common case first. */

if (Py_IS_FINITE(x) && Py_IS_FINITE(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 (Py_IS_NAN(x)) {

return x;

}

if (Py_IS_NAN(y)) {

return y;

}

if (Py_IS_INFINITY(x)) {

return Py_NAN;

}

assert(Py_IS_INFINITY(y));

return 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) {

return log2(x);

}

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(PyObject *module, PyObject * const *args, Py_ssize_t nargs)

{

PyObject *res, *x;

Py_ssize_t i;

if (nargs == 0) {

return PyLong_FromLong(0);

}

res = PyNumber_Index(args[0]);

if (res == NULL) {

return NULL;

}

if (nargs == 1) {

Py_SETREF(res, PyNumber_Absolute(res));

return res;

}

PyObject *one = _PyLong_GetOne(); // borrowed ref

for (i = 1; i < nargs; i++) {

x = _PyNumber_Index(args[i]);

if (x == NULL) {

Py_DECREF(res);

return NULL;

}

if (res == one) {

/* Fast path: just check arguments.

Py_DECREF(x);

continue;

}

Py_SETREF(res, _PyLong_GCD(res, x));

Py_DECREF(x);

if (res == NULL) {

return NULL;

}

}

return res;

}

PyDoc_STRVAR(math_gcd_doc,

"Greatest Common Divisor.");

static PyObject *

long_lcm(PyObject *a, PyObject *b)

{

PyObject *g, *m, *f, *ab;

if (_PyLong_IsZero((PyLongObject *)a) || _PyLong_IsZero((PyLongObject *)b)) {

return PyLong_FromLong(0);

}

g = _PyLong_GCD(a, b);

if (g == NULL) {

return NULL;

}

f = PyNumber_FloorDivide(a, g);

Py_DECREF(g);

if (f == NULL) {

return NULL;

}

m = PyNumber_Multiply(f, b);

Py_DECREF(f);

if (m == NULL) {

return NULL;

}

ab = PyNumber_Absolute(m);

Py_DECREF(m);

return ab;

}

static PyObject *

math_lcm(PyObject *module, PyObject * const *args, Py_ssize_t nargs)

{

PyObject *res, *x;

Py_ssize_t i;

if (nargs == 0) {

return PyLong_FromLong(1);

}

res = PyNumber_Index(args[0]);

if (res == NULL) {

return NULL;

}

if (nargs == 1) {

Py_SETREF(res, PyNumber_Absolute(res));

return res;

}

PyObject *zero = _PyLong_GetZero(); // borrowed ref

for (i = 1; i < nargs; i++) {

x = PyNumber_Index(args[i]);

if (x == NULL) {

Py_DECREF(res);

return NULL;

}

if (res == zero) {

/* Fast path: just check arguments.

Py_DECREF(x);

continue;

}

Py_SETREF(res, long_lcm(res, x));

Py_DECREF(x);

if (res == NULL) {

return NULL;

}

}

return res;

}

PyDoc_STRVAR(math_lcm_doc,

"Least Common Multiple.");

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

{

double x, r;

x = PyFloat_AsDouble(arg);

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

return NULL;

errno = 0;

r = (*func)(x);

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

r = (*func)(x);

if (errno && is_error(r))

return NULL;

return PyFloat_FromDouble(r);