#include "Python.h"

#include "structseq.h"

#include <ctype.h>

#include <float.h>

#undef MAX

#undef MIN

#define MAX(x, y) ((x) < (y) ? (y) : (x))

#define MIN(x, y) ((x) < (y) ? (x) : (y))

#ifdef HAVE_IEEEFP_H

#include <ieeefp.h>

#ifdef _OSF_SOURCE

extern int finite(double);

#define BLOCK_SIZE 1000 /* 1K less typical malloc overhead */

#define BHEAD_SIZE 8 /* Enough for a 64-bit pointer */

#define N_FLOATOBJECTS ((BLOCK_SIZE - BHEAD_SIZE) / sizeof(PyFloatObject))

struct _floatblock {

struct _floatblock *next;

PyFloatObject objects[N_FLOATOBJECTS];

};

typedef struct _floatblock PyFloatBlock;

static PyFloatBlock *block_list = NULL;

static PyFloatObject *free_list = NULL;

static PyFloatObject *

fill_free_list(void)

{

PyFloatObject *p, *q;

/* XXX Float blocks escape the object heap. Use PyObject_MALLOC ??? */

p = (PyFloatObject *) PyMem_MALLOC(sizeof(PyFloatBlock));

if (p == NULL)

return (PyFloatObject *) PyErr_NoMemory();

((PyFloatBlock *)p)->next = block_list;

block_list = (PyFloatBlock *)p;

p = &((PyFloatBlock *)p)->objects[0];

q = p + N_FLOATOBJECTS;

while (--q > p)

Py_TYPE(q) = (struct _typeobject *)(q-1);

Py_TYPE(q) = NULL;

return p + N_FLOATOBJECTS - 1;

}

PyFloat_GetMax(void)

{

return DBL_MAX;

}

PyFloat_GetMin(void)

{

return DBL_MIN;

}

static PyTypeObject FloatInfoType = {0, 0, 0, 0, 0, 0};

PyDoc_STRVAR(floatinfo__doc__,

your system's :file:`float.h` for more information.");

static PyStructSequence_Field floatinfo_fields[] = {

{"max", "DBL_MAX -- maximum representable finite float"},

{"max_exp", "DBL_MAX_EXP -- maximum int e such that radix**(e-1) "

"is representable"},

{"max_10_exp", "DBL_MAX_10_EXP -- maximum int e such that 10**e "

"is representable"},

{"min", "DBL_MIN -- Minimum positive normalizer float"},

{"min_exp", "DBL_MIN_EXP -- minimum int e such that radix**(e-1) "

"is a normalized float"},

{"min_10_exp", "DBL_MIN_10_EXP -- minimum int e such that 10**e is "

"a normalized"},

{"dig", "DBL_DIG -- digits"},

{"mant_dig", "DBL_MANT_DIG -- mantissa digits"},

{"epsilon", "DBL_EPSILON -- Difference between 1 and the next "

"representable float"},

{"radix", "FLT_RADIX -- radix of exponent"},

{"rounds", "FLT_ROUNDS -- addition rounds"},

{0}

};

static PyStructSequence_Desc floatinfo_desc = {

"sys.floatinfo", /* name */

floatinfo__doc__, /* doc */

floatinfo_fields, /* fields */

11

};

PyObject *

PyFloat_GetInfo(void)

{

PyObject* floatinfo;

int pos = 0;

floatinfo = PyStructSequence_New(&FloatInfoType);

if (floatinfo == NULL) {

return NULL;

}

#define SetIntFlag(flag) \

PyStructSequence_SET_ITEM(floatinfo, pos++, PyInt_FromLong(flag))

#define SetDblFlag(flag) \

PyStructSequence_SET_ITEM(floatinfo, pos++, PyFloat_FromDouble(flag))

SetDblFlag(DBL_MAX);

SetIntFlag(DBL_MAX_EXP);

SetIntFlag(DBL_MAX_10_EXP);

SetDblFlag(DBL_MIN);

SetIntFlag(DBL_MIN_EXP);

SetIntFlag(DBL_MIN_10_EXP);

SetIntFlag(DBL_DIG);

SetIntFlag(DBL_MANT_DIG);

SetDblFlag(DBL_EPSILON);

SetIntFlag(FLT_RADIX);

SetIntFlag(FLT_ROUNDS);

#undef SetIntFlag

#undef SetDblFlag

if (PyErr_Occurred()) {

Py_CLEAR(floatinfo);

return NULL;

}

return floatinfo;

}

PyObject *

PyFloat_FromDouble(double fval)

{

register PyFloatObject *op;

if (free_list == NULL) {

if ((free_list = fill_free_list()) == NULL)

return NULL;

}

/* Inline PyObject_New */

op = free_list;

free_list = (PyFloatObject *)Py_TYPE(op);

PyObject_INIT(op, &PyFloat_Type);

op->ob_fval = fval;

return (PyObject *) op;

}

PyObject *

PyFloat_FromString(PyObject *v, char **pend)

{

const char *s, *last, *end, *sp;

double x;

char buffer[256]; /* for errors */

#ifdef Py_USING_UNICODE

char s_buffer[256]; /* for objects convertible to a char buffer */

Py_ssize_t len;

if (pend)

*pend = NULL;

if (PyString_Check(v)) {

s = PyString_AS_STRING(v);

len = PyString_GET_SIZE(v);

}

#ifdef Py_USING_UNICODE

else if (PyUnicode_Check(v)) {

if (PyUnicode_GET_SIZE(v) >= (Py_ssize_t)sizeof(s_buffer)) {

PyErr_SetString(PyExc_ValueError,

"Unicode float() literal too long to convert");

return NULL;

}

if (PyUnicode_EncodeDecimal(PyUnicode_AS_UNICODE(v),

PyUnicode_GET_SIZE(v),

s_buffer,

NULL))

return NULL;

s = s_buffer;

len = strlen(s);

}

else if (PyObject_AsCharBuffer(v, &s, &len)) {

PyErr_SetString(PyExc_TypeError,

"float() argument must be a string or a number");

return NULL;

}

last = s + len;

while (*s && isspace(Py_CHARMASK(*s)))

s++;

if (*s == '\0') {

PyErr_SetString(PyExc_ValueError, "empty string for float()");

return NULL;

}

sp = s;

/* We don't care about overflow or underflow. If the platform supports

PyFPE_START_PROTECT("strtod", return NULL)

x = PyOS_ascii_strtod(s, (char **)&end);

PyFPE_END_PROTECT(x)

errno = 0;

/* Believe it or not, Solaris 2.6 can move end *beyond* the null

if (end > last)

end = last;

/* Check for inf and nan. This is done late because it rarely happens. */

if (end == s) {

char *p = (char*)sp;

int sign = 1;

if (*p == '-') {

sign = -1;

p++;

}

if (*p == '+') {

p++;

}

if (PyOS_strnicmp(p, "inf", 4) == 0) {

Py_RETURN_INF(sign);

}

if (PyOS_strnicmp(p, "infinity", 9) == 0) {

Py_RETURN_INF(sign);

}

#ifdef Py_NAN

if(PyOS_strnicmp(p, "nan", 4) == 0) {

Py_RETURN_NAN;

}

PyOS_snprintf(buffer, sizeof(buffer),

"invalid literal for float(): %.200s", s);

PyErr_SetString(PyExc_ValueError, buffer);

return NULL;

}

/* Since end != s, the platform made *some* kind of sense out

while (*end && isspace(Py_CHARMASK(*end)))

end++;

if (*end != '\0') {

PyOS_snprintf(buffer, sizeof(buffer),

"invalid literal for float(): %.200s", s);

PyErr_SetString(PyExc_ValueError, buffer);

return NULL;

}

else if (end != last) {

PyErr_SetString(PyExc_ValueError,

"null byte in argument for float()");

return NULL;

}

if (x == 0.0) {

/* See above -- may have been strtod being anal

PyFPE_START_PROTECT("atof", return NULL)

x = PyOS_ascii_atof(s);

PyFPE_END_PROTECT(x)

errno = 0; /* whether atof ever set errno is undefined */

}

return PyFloat_FromDouble(x);

}

static void

float_dealloc(PyFloatObject *op)

{

if (PyFloat_CheckExact(op)) {

Py_TYPE(op) = (struct _typeobject *)free_list;

free_list = op;

}

else

Py_TYPE(op)->tp_free((PyObject *)op);

}

PyFloat_AsDouble(PyObject *op)

{

PyNumberMethods *nb;

PyFloatObject *fo;

double val;

if (op && PyFloat_Check(op))

return PyFloat_AS_DOUBLE((PyFloatObject*) op);

if (op == NULL) {

PyErr_BadArgument();

return -1;

}

if ((nb = Py_TYPE(op)->tp_as_number) == NULL || nb->nb_float == NULL) {

PyErr_SetString(PyExc_TypeError, "a float is required");

return -1;

}

fo = (PyFloatObject*) (*nb->nb_float) (op);

if (fo == NULL)

return -1;

if (!PyFloat_Check(fo)) {

PyErr_SetString(PyExc_TypeError,

"nb_float should return float object");

return -1;

}

val = PyFloat_AS_DOUBLE(fo);

Py_DECREF(fo);

return val;

}

static void

format_float(char *buf, size_t buflen, PyFloatObject *v, int precision)

{

register char *cp;

char format[32];

int i;

/* Subroutine for float_repr and float_print.

assert(PyFloat_Check(v));

PyOS_snprintf(format, 32, "%%.%ig", precision);

PyOS_ascii_formatd(buf, buflen, format, v->ob_fval);

cp = buf;

if (*cp == '-')

cp++;

for (; *cp != '\0'; cp++) {

/* Any non-digit means it's not an integer;

if (!isdigit(Py_CHARMASK(*cp)))

break;

}

if (*cp == '\0') {

*cp++ = '.';

*cp++ = '0';

*cp++ = '\0';

return;

}

/* Checking the next three chars should be more than enough to

for (i=0; *cp != '\0' && i<3; cp++, i++) {

if (isdigit(Py_CHARMASK(*cp)) || *cp == '.')

continue;

/* found something that is neither a digit nor point

#ifdef Py_NAN

if (Py_IS_NAN(v->ob_fval)) {

strcpy(buf, "nan");

}

else

if (Py_IS_INFINITY(v->ob_fval)) {

cp = buf;

if (*cp == '-')

cp++;

strcpy(cp, "inf");

}

break;

}

}

PyFloat_AsStringEx(char *buf, PyFloatObject *v, int precision)

{

format_float(buf, 100, v, precision);

}

#define CONVERT_TO_DOUBLE(obj, dbl) \

if (PyFloat_Check(obj)) \

dbl = PyFloat_AS_DOUBLE(obj); \

else if (convert_to_double(&(obj), &(dbl)) < 0) \

return obj;

static int

convert_to_double(PyObject **v, double *dbl)

{

register PyObject *obj = *v;

if (PyInt_Check(obj)) {

*dbl = (double)PyInt_AS_LONG(obj);

}

else if (PyLong_Check(obj)) {

*dbl = PyLong_AsDouble(obj);

if (*dbl == -1.0 && PyErr_Occurred()) {

*v = NULL;

return -1;

}

}

else {

Py_INCREF(Py_NotImplemented);

*v = Py_NotImplemented;

return -1;

}

return 0;

}

#define PREC_REPR 17

#define PREC_STR 12

PyFloat_AsString(char *buf, PyFloatObject *v)

{

format_float(buf, 100, v, PREC_STR);

}

PyFloat_AsReprString(char *buf, PyFloatObject *v)

{

format_float(buf, 100, v, PREC_REPR);

}

static int

float_print(PyFloatObject *v, FILE *fp, int flags)

{

char buf[100];

format_float(buf, sizeof(buf), v,

(flags & Py_PRINT_RAW) ? PREC_STR : PREC_REPR);

Py_BEGIN_ALLOW_THREADS

fputs(buf, fp);

Py_END_ALLOW_THREADS

return 0;

}

static PyObject *

float_repr(PyFloatObject *v)

{

char buf[100];

format_float(buf, sizeof(buf), v, PREC_REPR);

return PyString_FromString(buf);

}

static PyObject *

float_str(PyFloatObject *v)

{

char buf[100];

format_float(buf, sizeof(buf), v, PREC_STR);

return PyString_FromString(buf);

}

static PyObject*

float_richcompare(PyObject *v, PyObject *w, int op)

{

double i, j;

int r = 0;

assert(PyFloat_Check(v));

i = PyFloat_AS_DOUBLE(v);

/* Switch on the type of w. Set i and j to doubles to be compared,

if (PyFloat_Check(w))

j = PyFloat_AS_DOUBLE(w);

else if (!Py_IS_FINITE(i)) {

if (PyInt_Check(w) || PyLong_Check(w))

/* If i is an infinity, its magnitude exceeds any

j = 0.0;

else

goto Unimplemented;

}

else if (PyInt_Check(w)) {

long jj = PyInt_AS_LONG(w);

/* In the worst realistic case I can imagine, C double is a

#if SIZEOF_LONG > 6

unsigned long abs = (unsigned long)(jj < 0 ? -jj : jj);

if (abs >> 48) {

/* Needs more than 48 bits. Make it take the

PyObject *result;

PyObject *ww = PyLong_FromLong(jj);

if (ww == NULL)

return NULL;

result = float_richcompare(v, ww, op);

Py_DECREF(ww);

return result;

}

j = (double)jj;

assert((long)j == jj);

}

else if (PyLong_Check(w)) {

int vsign = i == 0.0 ? 0 : i < 0.0 ? -1 : 1;

int wsign = _PyLong_Sign(w);

size_t nbits;

int exponent;

if (vsign != wsign) {

/* Magnitudes are irrelevant -- the signs alone

i = (double)vsign;

j = (double)wsign;

goto Compare;

}

/* The signs are the same. */

/* Convert w to a double if it fits. In particular, 0 fits. */

nbits = _PyLong_NumBits(w);

if (nbits == (size_t)-1 && PyErr_Occurred()) {

/* This long is so large that size_t isn't big enough

PyErr_Clear();

i = (double)vsign;

assert(wsign != 0);

j = wsign * 2.0;

goto Compare;

}

if (nbits <= 48) {

j = PyLong_AsDouble(w);

/* It's impossible that <= 48 bits overflowed. */

assert(j != -1.0 || ! PyErr_Occurred());

goto Compare;

}

assert(wsign != 0); /* else nbits was 0 */

assert(vsign != 0); /* if vsign were 0, then since wsign is

/* We want to work with non-negative numbers. */

if (vsign < 0) {

/* "Multiply both sides" by -1; this also swaps the

i = -i;

op = _Py_SwappedOp[op];

}

assert(i > 0.0);

(void) frexp(i, &exponent);

/* exponent is the # of bits in v before the radix point;

if (exponent < 0 || (size_t)exponent < nbits) {

i = 1.0;

j = 2.0;

goto Compare;

}

if ((size_t)exponent > nbits) {

i = 2.0;

j = 1.0;

goto Compare;

}

/* v and w have the same number of bits before the radix

{

double fracpart;

double intpart;

PyObject *result = NULL;

PyObject *one = NULL;

PyObject *vv = NULL;

PyObject *ww = w;

if (wsign < 0) {

ww = PyNumber_Negative(w);

if (ww == NULL)

goto Error;

}

else

Py_INCREF(ww);

fracpart = modf(i, &intpart);

vv = PyLong_FromDouble(intpart);

if (vv == NULL)

goto Error;

if (fracpart != 0.0) {

/* Shift left, and or a 1 bit into vv

PyObject *temp;

one = PyInt_FromLong(1);

if (one == NULL)

goto Error;

temp = PyNumber_Lshift(ww, one);

if (temp == NULL)

goto Error;

Py_DECREF(ww);

ww = temp;

temp = PyNumber_Lshift(vv, one);

if (temp == NULL)

goto Error;

Py_DECREF(vv);

vv = temp;

temp = PyNumber_Or(vv, one);

if (temp == NULL)

goto Error;

Py_DECREF(vv);

vv = temp;

}

r = PyObject_RichCompareBool(vv, ww, op);

if (r < 0)

goto Error;

result = PyBool_FromLong(r);

Error:

Py_XDECREF(vv);

Py_XDECREF(ww);

Py_XDECREF(one);

return result;

}

} /* else if (PyLong_Check(w)) */

else /* w isn't float, int, or long */

goto Unimplemented;

Compare:

PyFPE_START_PROTECT("richcompare", return NULL)

switch (op) {

case Py_EQ:

r = i == j;

break;

case Py_NE:

r = i != j;

break;

case Py_LE:

r = i <= j;

break;

case Py_GE:

r = i >= j;

break;

case Py_LT:

r = i < j;

break;

case Py_GT:

r = i > j;

break;

}

PyFPE_END_PROTECT(r)

return PyBool_FromLong(r);

Unimplemented:

Py_INCREF(Py_NotImplemented);

return Py_NotImplemented;

}

static long

float_hash(PyFloatObject *v)

{

return _Py_HashDouble(v->ob_fval);

}

static PyObject *

float_add(PyObject *v, PyObject *w)

{

double a,b;

CONVERT_TO_DOUBLE(v, a);

CONVERT_TO_DOUBLE(w, b);

PyFPE_START_PROTECT("add", return 0)

a = a + b;

PyFPE_END_PROTECT(a)

return PyFloat_FromDouble(a);

}

static PyObject *

float_sub(PyObject *v, PyObject *w)

{

double a,b;

CONVERT_TO_DOUBLE(v, a);

CONVERT_TO_DOUBLE(w, b);

PyFPE_START_PROTECT("subtract", return 0)

a = a - b;

PyFPE_END_PROTECT(a)

return PyFloat_FromDouble(a);

}

static PyObject *

float_mul(PyObject *v, PyObject *w)

{

double a,b;

CONVERT_TO_DOUBLE(v, a);

CONVERT_TO_DOUBLE(w, b);

PyFPE_START_PROTECT("multiply", return 0)

a = a * b;

PyFPE_END_PROTECT(a)

return PyFloat_FromDouble(a);

}

static PyObject *

float_div(PyObject *v, PyObject *w)

{

double a,b;

CONVERT_TO_DOUBLE(v, a);

CONVERT_TO_DOUBLE(w, b);

#ifdef Py_NAN

if (b == 0.0) {

PyErr_SetString(PyExc_ZeroDivisionError,

"float division");

return NULL;

}

PyFPE_START_PROTECT("divide", return 0)

a = a / b;

PyFPE_END_PROTECT(a)

return PyFloat_FromDouble(a);

}

static PyObject *

float_classic_div(PyObject *v, PyObject *w)

{

double a,b;

CONVERT_TO_DOUBLE(v, a);

CONVERT_TO_DOUBLE(w, b);

if (Py_DivisionWarningFlag >= 2 &&

PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)

return NULL;

#ifdef Py_NAN

if (b == 0.0) {

PyErr_SetString(PyExc_ZeroDivisionError,

"float division");

return NULL;

}

PyFPE_START_PROTECT("divide", return 0)

a = a / b;

PyFPE_END_PROTECT(a)

return PyFloat_FromDouble(a);

}

static PyObject *

float_rem(PyObject *v, PyObject *w)

{

double vx, wx;

double mod;

CONVERT_TO_DOUBLE(v, vx);

CONVERT_TO_DOUBLE(w, wx);

#ifdef Py_NAN

if (wx == 0.0) {

PyErr_SetString(PyExc_ZeroDivisionError,

"float modulo");

return NULL;

}

PyFPE_START_PROTECT("modulo", return 0)

mod = fmod(vx, wx);

/* note: checking mod*wx < 0 is incorrect -- underflows to

if (mod && ((wx < 0) != (mod < 0))) {

mod += wx;

}

PyFPE_END_PROTECT(mod)

return PyFloat_FromDouble(mod);

}

static PyObject *

float_divmod(PyObject *v, PyObject *w)

{

double vx, wx;

double div, mod, floordiv;

CONVERT_TO_DOUBLE(v, vx);

CONVERT_TO_DOUBLE(w, wx);

if (wx == 0.0) {

PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");

return NULL;

}

PyFPE_START_PROTECT("divmod", return 0)

mod = fmod(vx, wx);

/* fmod is typically exact, so vx-mod is *mathematically* an

div = (vx - mod) / wx;

if (mod) {

/* ensure the remainder has the same sign as the denominator */

if ((wx < 0) != (mod < 0)) {

mod += wx;

div -= 1.0;

}

}

else {

/* the remainder is zero, and in the presence of signed zeroes

mod *= mod; /* hide "mod = +0" from optimizer */

if (wx < 0.0)

mod = -mod;

}

/* snap quotient to nearest integral value */

if (div) {

floordiv = floor(div);

if (div - floordiv > 0.5)

floordiv += 1.0;

}

else {

/* div is zero - get the same sign as the true quotient */

div *= div; /* hide "div = +0" from optimizers */

floordiv = div * vx / wx; /* zero w/ sign of vx/wx */

}

PyFPE_END_PROTECT(floordiv)

return Py_BuildValue("(dd)", floordiv, mod);

}

static PyObject *

float_floor_div(PyObject *v, PyObject *w)

{

PyObject *t, *r;

t = float_divmod(v, w);

if (t == NULL || t == Py_NotImplemented)

return t;

assert(PyTuple_CheckExact(t));

r = PyTuple_GET_ITEM(t, 0);

Py_INCREF(r);

Py_DECREF(t);

return r;

}

static PyObject *

float_pow(PyObject *v, PyObject *w, PyObject *z)

{

double iv, iw, ix;

if ((PyObject *)z != Py_None) {

PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not "

"allowed unless all arguments are integers");

return NULL;

}

CONVERT_TO_DOUBLE(v, iv);

CONVERT_TO_DOUBLE(w, iw);

/* Sort out special cases here instead of relying on pow() */

if (iw == 0) { /* v**0 is 1, even 0**0 */

return PyFloat_FromDouble(1.0);

}

if (iv == 0.0) { /* 0**w is error if w<0, else 1 */

if (iw < 0.0) {

PyErr_SetString(PyExc_ZeroDivisionError,

"0.0 cannot be raised to a negative power");

return NULL;

}

return PyFloat_FromDouble(0.0);

}

if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */

return PyFloat_FromDouble(1.0);

}

if (iv < 0.0) {

/* Whether this is an error is a mess, and bumps into libm

if (iw != floor(iw)) {

PyErr_SetString(PyExc_ValueError, "negative number "

"cannot be raised to a fractional power");

return NULL;

}

/* iw is an exact integer, albeit perhaps a very large one.

if (iv == -1.0 && Py_IS_FINITE(iw)) {

/* Return 1 if iw is even, -1 if iw is odd; there's

ix = floor(iw * 0.5) * 2.0;

return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0);

}

/* Else iv != -1.0, and overflow or underflow are possible.

}

errno = 0;

PyFPE_START_PROTECT("pow", return NULL)

ix = pow(iv, iw);

PyFPE_END_PROTECT(ix)

Py_ADJUST_ERANGE1(ix);

if (errno != 0) {

/* We don't expect any errno value other than ERANGE, but

PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :

PyExc_ValueError);

return NULL;

}

return PyFloat_FromDouble(ix);

}

static PyObject *

float_neg(PyFloatObject *v)

{

return PyFloat_FromDouble(-v->ob_fval);

}

static PyObject *

float_abs(PyFloatObject *v)

{

return PyFloat_FromDouble(fabs(v->ob_fval));

}

static int

float_nonzero(PyFloatObject *v)

{