From 8aa0882402ededd8f1a82d739c8df655973095cc Mon Sep 17 00:00:00 2001 From: Kevin Modzelewski <kmod@dropbox.com> Date: Thu, 7 May 2015 16:09:09 -0700 Subject: [PATCH] Delete capi/float.cpp and move its contents to runtime/float.cpp I don't remember the original reasoning for having the separate capi/float.cpp, and all the other runtime/ files have their CAPI methods inline. --- src/CMakeLists.txt | 1 - src/capi/float.cpp | 673 ------------------------------------------ src/runtime/float.cpp | 563 +++++++++++++++++++++++++++++++++++ 3 files changed, 563 insertions(+), 674 deletions(-) delete mode 100644 src/capi/float.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6c499955..93d02926 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -30,7 +30,6 @@ add_library(PYSTON_OBJECTS OBJECT ${OPTIONAL_SRCS} capi/codecs.cpp capi/descrobject.cpp capi/errors.cpp - capi/float.cpp capi/modsupport.cpp capi/object.cpp capi/typeobject.cpp diff --git a/src/capi/float.cpp b/src/capi/float.cpp deleted file mode 100644 index 3110b2cb..00000000 --- a/src/capi/float.cpp +++ /dev/null @@ -1,673 +0,0 @@ -// Copyright (c) 2014-2015 Dropbox, Inc. -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// -// http://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. - -// This file is mostly copied from CPython - -#include <cfloat> - -#include "Python.h" - -#include "core/types.h" -#include "runtime/types.h" - -namespace pyston { - -extern "C" { - -typedef enum { unknown_format, ieee_big_endian_format, ieee_little_endian_format } float_format_type; - -static float_format_type _detect_double_format() { - float_format_type format; -#if SIZEOF_DOUBLE == 8 - { - double x = 9006104071832581.0; - if (memcmp(&x, "\x43\x3f\xff\x01\x02\x03\x04\x05", 8) == 0) - format = ieee_big_endian_format; - else if (memcmp(&x, "\x05\x04\x03\x02\x01\xff\x3f\x43", 8) == 0) - format = ieee_little_endian_format; - else - format = unknown_format; - } -#else - format = unknown_format; -#endif - return format; -} - -static float_format_type _detect_float_format() { - float_format_type format; -#if SIZEOF_FLOAT == 4 - { - float y = 16711938.0; - if (memcmp(&y, "\x4b\x7f\x01\x02", 4) == 0) - format = ieee_big_endian_format; - else if (memcmp(&y, "\x02\x01\x7f\x4b", 4) == 0) - format = ieee_little_endian_format; - else - format = unknown_format; - } -#else - format = unknown_format; -#endif - return format; -} - -static float_format_type double_format, float_format; -static float_format_type detected_double_format = _detect_double_format(), - detected_float_format = _detect_float_format(); - -static PyObject* float_getformat(PyTypeObject* v, PyObject* arg) noexcept { - float_format_type r; - - BoxedString* str = static_cast<BoxedString*>(arg); - - if (!PyString_Check(arg)) { - PyErr_Format(PyExc_TypeError, "__getformat__() argument must be string, not %.500s", Py_TYPE(arg)->tp_name); - return NULL; - } - if (str->s == "double") { - r = double_format; - } else if (str->s == "float") { - r = float_format; - } else { - PyErr_SetString(PyExc_ValueError, "__getformat__() argument 1 must be " - "'double' or 'float'"); - return NULL; - } - - switch (r) { - case unknown_format: - return static_cast<PyObject*>(boxStrConstant("unknown")); - case ieee_little_endian_format: - return static_cast<PyObject*>(boxStrConstant("IEEE, little-endian")); - case ieee_big_endian_format: - return static_cast<PyObject*>(boxStrConstant("IEEE, big-endian")); - default: - Py_FatalError("insane float_format or double_format"); - return NULL; - } -} - -PyDoc_STRVAR(float_getformat_doc, "float.__getformat__(typestr) -> string\n" - "\n" - "You probably don't want to use this function. It exists mainly to be\n" - "used in Python's test suite.\n" - "\n" - "typestr must be 'double' or 'float'. This function returns whichever of\n" - "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n" - "format of floating point numbers used by the C type named by typestr."); - -static PyObject* float_setformat(PyTypeObject* v, PyObject* args) noexcept { - char* typestr; - char* format; - float_format_type f; - float_format_type detected; - float_format_type* p; - - if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format)) - return NULL; - - if (strcmp(typestr, "double") == 0) { - p = &double_format; - detected = detected_double_format; - } else if (strcmp(typestr, "float") == 0) { - p = &float_format; - detected = detected_float_format; - } else { - PyErr_SetString(PyExc_ValueError, "__setformat__() argument 1 must " - "be 'double' or 'float'"); - return NULL; - } - - if (strcmp(format, "unknown") == 0) { - f = unknown_format; - } else if (strcmp(format, "IEEE, little-endian") == 0) { - f = ieee_little_endian_format; - } else if (strcmp(format, "IEEE, big-endian") == 0) { - f = ieee_big_endian_format; - } else { - PyErr_SetString(PyExc_ValueError, "__setformat__() argument 2 must be " - "'unknown', 'IEEE, little-endian' or " - "'IEEE, big-endian'"); - return NULL; - } - - if (f != unknown_format && f != detected) { - PyErr_Format(PyExc_ValueError, "can only set %s format to 'unknown' or the " - "detected platform value", - typestr); - return NULL; - } - - *p = f; - Py_RETURN_NONE; -} - -/*---------------------------------------------------------------------------- - * _PyFloat_{Pack,Unpack}{4,8}. See floatobject.h. - */ -extern "C" int _PyFloat_Pack4(double x, unsigned char* p, int le) noexcept { - if (float_format == unknown_format) { - unsigned char sign; - int e; - double f; - unsigned int fbits; - int incr = 1; - - if (le) { - p += 3; - incr = -1; - } - - if (x < 0) { - sign = 1; - x = -x; - } else - sign = 0; - - f = frexp(x, &e); - - /* Normalize f to be in the range [1.0, 2.0) */ - if (0.5 <= f && f < 1.0) { - f *= 2.0; - e--; - } else if (f == 0.0) - e = 0; - else { - PyErr_SetString(PyExc_SystemError, "frexp() result out of range"); - return -1; - } - - if (e >= 128) - goto Overflow; - else if (e < -126) { - /* Gradual underflow */ - f = ldexp(f, 126 + e); - e = 0; - } else if (!(e == 0 && f == 0.0)) { - e += 127; - f -= 1.0; /* Get rid of leading 1 */ - } - - f *= 8388608.0; /* 2**23 */ - fbits = (unsigned int)(f + 0.5); /* Round */ - assert(fbits <= 8388608); - if (fbits >> 23) { - /* The carry propagated out of a string of 23 1 bits. */ - fbits = 0; - ++e; - if (e >= 255) - goto Overflow; - } - - /* First byte */ - *p = (sign << 7) | (e >> 1); - p += incr; - - /* Second byte */ - *p = (char)(((e & 1) << 7) | (fbits >> 16)); - p += incr; - - /* Third byte */ - *p = (fbits >> 8) & 0xFF; - p += incr; - - /* Fourth byte */ - *p = fbits & 0xFF; - - /* Done */ - return 0; - - } else { - float y = (float)x; - const char* s = (char*)&y; - int i, incr = 1; - - if (Py_IS_INFINITY(y) && !Py_IS_INFINITY(x)) - goto Overflow; - - if ((float_format == ieee_little_endian_format && !le) || (float_format == ieee_big_endian_format && le)) { - p += 3; - incr = -1; - } - - for (i = 0; i < 4; i++) { - *p = *s++; - p += incr; - } - return 0; - } -Overflow: - PyErr_SetString(PyExc_OverflowError, "float too large to pack with f format"); - return -1; -} - -extern "C" int _PyFloat_Pack8(double x, unsigned char* p, int le) noexcept { - if (double_format == unknown_format) { - unsigned char sign; - int e; - double f; - unsigned int fhi, flo; - int incr = 1; - - if (le) { - p += 7; - incr = -1; - } - - if (x < 0) { - sign = 1; - x = -x; - } else - sign = 0; - - f = frexp(x, &e); - - /* Normalize f to be in the range [1.0, 2.0) */ - if (0.5 <= f && f < 1.0) { - f *= 2.0; - e--; - } else if (f == 0.0) - e = 0; - else { - PyErr_SetString(PyExc_SystemError, "frexp() result out of range"); - return -1; - } - - if (e >= 1024) - goto Overflow; - else if (e < -1022) { - /* Gradual underflow */ - f = ldexp(f, 1022 + e); - e = 0; - } else if (!(e == 0 && f == 0.0)) { - e += 1023; - f -= 1.0; /* Get rid of leading 1 */ - } - - /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */ - f *= 268435456.0; /* 2**28 */ - fhi = (unsigned int)f; /* Truncate */ - assert(fhi < 268435456); - - f -= (double)fhi; - f *= 16777216.0; /* 2**24 */ - flo = (unsigned int)(f + 0.5); /* Round */ - assert(flo <= 16777216); - if (flo >> 24) { - /* The carry propagated out of a string of 24 1 bits. */ - flo = 0; - ++fhi; - if (fhi >> 28) { - /* And it also progagated out of the next 28 bits. */ - fhi = 0; - ++e; - if (e >= 2047) - goto Overflow; - } - } - - /* First byte */ - *p = (sign << 7) | (e >> 4); - p += incr; - - /* Second byte */ - *p = (unsigned char)(((e & 0xF) << 4) | (fhi >> 24)); - p += incr; - - /* Third byte */ - *p = (fhi >> 16) & 0xFF; - p += incr; - - /* Fourth byte */ - *p = (fhi >> 8) & 0xFF; - p += incr; - - /* Fifth byte */ - *p = fhi & 0xFF; - p += incr; - - /* Sixth byte */ - *p = (flo >> 16) & 0xFF; - p += incr; - - /* Seventh byte */ - *p = (flo >> 8) & 0xFF; - p += incr; - - /* Eighth byte */ - *p = flo & 0xFF; - /* p += incr; Unneeded (for now) */ - - /* Done */ - return 0; - - Overflow: - PyErr_SetString(PyExc_OverflowError, "float too large to pack with d format"); - return -1; - } else { - const char* s = (char*)&x; - int i, incr = 1; - - if ((double_format == ieee_little_endian_format && !le) || (double_format == ieee_big_endian_format && le)) { - p += 7; - incr = -1; - } - - for (i = 0; i < 8; i++) { - *p = *s++; - p += incr; - } - return 0; - } -} - -double _PyFloat_Unpack4(const unsigned char* p, int le) noexcept { - if (float_format == unknown_format) { - unsigned char sign; - int e; - unsigned int f; - double x; - int incr = 1; - - if (le) { - p += 3; - incr = -1; - } - - /* First byte */ - sign = (*p >> 7) & 1; - e = (*p & 0x7F) << 1; - p += incr; - - /* Second byte */ - e |= (*p >> 7) & 1; - f = (*p & 0x7F) << 16; - p += incr; - - if (e == 255) { - PyErr_SetString(PyExc_ValueError, "can't unpack IEEE 754 special value " - "on non-IEEE platform"); - return -1; - } - - /* Third byte */ - f |= *p << 8; - p += incr; - - /* Fourth byte */ - f |= *p; - - x = (double)f / 8388608.0; - - /* XXX This sadly ignores Inf/NaN issues */ - if (e == 0) - e = -126; - else { - x += 1.0; - e -= 127; - } - x = ldexp(x, e); - - if (sign) - x = -x; - - return x; - } else { - float x; - - if ((float_format == ieee_little_endian_format && !le) || (float_format == ieee_big_endian_format && le)) { - char buf[4]; - char* d = &buf[3]; - int i; - - for (i = 0; i < 4; i++) { - *d-- = *p++; - } - memcpy(&x, buf, 4); - } else { - memcpy(&x, p, 4); - } - - return x; - } -} - -double _PyFloat_Unpack8(const unsigned char* p, int le) noexcept { - if (double_format == unknown_format) { - unsigned char sign; - int e; - unsigned int fhi, flo; - double x; - int incr = 1; - - if (le) { - p += 7; - incr = -1; - } - - /* First byte */ - sign = (*p >> 7) & 1; - e = (*p & 0x7F) << 4; - - p += incr; - - /* Second byte */ - e |= (*p >> 4) & 0xF; - fhi = (*p & 0xF) << 24; - p += incr; - - if (e == 2047) { - PyErr_SetString(PyExc_ValueError, "can't unpack IEEE 754 special value " - "on non-IEEE platform"); - return -1.0; - } - - /* Third byte */ - fhi |= *p << 16; - p += incr; - - /* Fourth byte */ - fhi |= *p << 8; - p += incr; - - /* Fifth byte */ - fhi |= *p; - p += incr; - - /* Sixth byte */ - flo = *p << 16; - p += incr; - - /* Seventh byte */ - flo |= *p << 8; - p += incr; - - /* Eighth byte */ - flo |= *p; - - x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */ - x /= 268435456.0; /* 2**28 */ - - if (e == 0) - e = -1022; - else { - x += 1.0; - e -= 1023; - } - x = ldexp(x, e); - - if (sign) - x = -x; - - return x; - } else { - double x; - - if ((double_format == ieee_little_endian_format && !le) || (double_format == ieee_big_endian_format && le)) { - char buf[8]; - char* d = &buf[7]; - int i; - - for (i = 0; i < 8; i++) { - *d-- = *p++; - } - memcpy(&x, buf, 8); - } else { - memcpy(&x, p, 8); - } - - return x; - } -} - -#if DBL_MANT_DIG == 53 -#define FIVE_POW_LIMIT 22 -#else -#error "C doubles do not appear to be IEEE 754 binary64 format" -#endif - -PyObject* _Py_double_round(double x, int ndigits) noexcept { - - double rounded, m; - Py_ssize_t buflen, mybuflen = 100; - char* buf, *buf_end, shortbuf[100], * mybuf = shortbuf; - int decpt, sign, val, halfway_case; - PyObject* result = NULL; - _Py_SET_53BIT_PRECISION_HEADER; - - /* Easy path for the common case ndigits == 0. */ - if (ndigits == 0) { - rounded = round(x); - if (fabs(rounded - x) == 0.5) - /* halfway between two integers; use round-away-from-zero */ - rounded = x + (x > 0.0 ? 0.5 : -0.5); - return PyFloat_FromDouble(rounded); - } - - /* The basic idea is very simple: convert and round the double to a - decimal string using _Py_dg_dtoa, then convert that decimal string - back to a double with _Py_dg_strtod. There's one minor difficulty: - Python 2.x expects round to do round-half-away-from-zero, while - _Py_dg_dtoa does round-half-to-even. So we need some way to detect - and correct the halfway cases. - - Detection: a halfway value has the form k * 0.5 * 10**-ndigits for - some odd integer k. Or in other words, a rational number x is - exactly halfway between two multiples of 10**-ndigits if its - 2-valuation is exactly -ndigits-1 and its 5-valuation is at least - -ndigits. For ndigits >= 0 the latter condition is automatically - satisfied for a binary float x, since any such float has - nonnegative 5-valuation. For 0 > ndigits >= -22, x needs to be an - integral multiple of 5**-ndigits; we can check this using fmod. - For -22 > ndigits, there are no halfway cases: 5**23 takes 54 bits - to represent exactly, so any odd multiple of 0.5 * 10**n for n >= - 23 takes at least 54 bits of precision to represent exactly. - - Correction: a simple strategy for dealing with halfway cases is to - (for the halfway cases only) call _Py_dg_dtoa with an argument of - ndigits+1 instead of ndigits (thus doing an exact conversion to - decimal), round the resulting string manually, and then convert - back using _Py_dg_strtod. - */ - - /* nans, infinities and zeros should have already been dealt - with by the caller (in this case, builtin_round) */ - assert(std::isfinite(x) && x != 0.0); - - /* find 2-valuation val of x */ - m = frexp(x, &val); - while (m != floor(m)) { - m *= 2.0; - val--; - } - - /* determine whether this is a halfway case */ - if (val == -ndigits - 1) { - if (ndigits >= 0) - halfway_case = 1; - else if (ndigits >= -FIVE_POW_LIMIT) { - double five_pow = 1.0; - int i; - for (i = 0; i < -ndigits; i++) - five_pow *= 5.0; - halfway_case = fmod(x, five_pow) == 0.0; - } else - halfway_case = 0; - } else - halfway_case = 0; - - /* round to a decimal string; use an extra place for halfway case */ - _Py_SET_53BIT_PRECISION_START; - buf = _Py_dg_dtoa(x, 3, ndigits + halfway_case, &decpt, &sign, &buf_end); - _Py_SET_53BIT_PRECISION_END; - if (buf == NULL) { - PyErr_NoMemory(); - return NULL; - } - buflen = buf_end - buf; - - /* in halfway case, do the round-half-away-from-zero manually */ - if (halfway_case) { - int i, carry; - /* sanity check: _Py_dg_dtoa should not have stripped - any zeros from the result: there should be exactly - ndigits+1 places following the decimal point, and - the last digit in the buffer should be a '5'.*/ - assert(buflen - decpt == ndigits + 1); - assert(buf[buflen - 1] == '5'); - - /* increment and shift right at the same time. */ - decpt += 1; - carry = 1; - for (i = buflen - 1; i-- > 0;) { - carry += buf[i] - '0'; - buf[i + 1] = carry % 10 + '0'; - carry /= 10; - } - buf[0] = carry + '0'; - } - - /* Get new buffer if shortbuf is too small. Space needed <= buf_end - - buf + 8: (1 extra for '0', 1 for sign, 5 for exp, 1 for '\0'). */ - if (buflen + 8 > mybuflen) { - mybuflen = buflen + 8; - mybuf = (char*)PyMem_Malloc(mybuflen); - if (mybuf == NULL) { - PyErr_NoMemory(); - goto exit; - } - } - /* copy buf to mybuf, adding exponent, sign and leading 0 */ - PyOS_snprintf(mybuf, mybuflen, "%s0%se%d", (sign ? "-" : ""), buf, decpt - (int)buflen); - - /* and convert the resulting string back to a double */ - errno = 0; - _Py_SET_53BIT_PRECISION_START; - rounded = _Py_dg_strtod(mybuf, NULL); - _Py_SET_53BIT_PRECISION_END; - if (errno == ERANGE && fabs(rounded) >= 1.) - PyErr_SetString(PyExc_OverflowError, "rounded value too large to represent"); - else - result = PyFloat_FromDouble(rounded); - - /* done computing value; now clean up */ - if (mybuf != shortbuf) - PyMem_Free(mybuf); -exit: - _Py_dg_freedtoa(buf); - return result; -} -} -} diff --git a/src/runtime/float.cpp b/src/runtime/float.cpp index 503b3651..caf6a098 100644 --- a/src/runtime/float.cpp +++ b/src/runtime/float.cpp @@ -12,6 +12,7 @@ // See the License for the specific language governing permissions and // limitations under the License. +#include <cfloat> #include <cmath> #include <cstring> @@ -833,6 +834,568 @@ const char* floatGetFormatDoc = "float.__getformat__(typestr) -> string\n" "'unknown', 'IEEE, big-endian' or 'IEEE, little-endian' best describes the\n" "format of floating point numbers used by the C type named by typestr."; +static PyObject* float_setformat(PyTypeObject* v, PyObject* args) noexcept { + char* typestr; + char* format; + float_format_type f; + float_format_type detected; + float_format_type* p; + + if (!PyArg_ParseTuple(args, "ss:__setformat__", &typestr, &format)) + return NULL; + + if (strcmp(typestr, "double") == 0) { + p = &double_format; + detected = detected_double_format; + } else if (strcmp(typestr, "float") == 0) { + p = &float_format; + detected = detected_float_format; + } else { + PyErr_SetString(PyExc_ValueError, "__setformat__() argument 1 must " + "be 'double' or 'float'"); + return NULL; + } + + if (strcmp(format, "unknown") == 0) { + f = unknown_format; + } else if (strcmp(format, "IEEE, little-endian") == 0) { + f = ieee_little_endian_format; + } else if (strcmp(format, "IEEE, big-endian") == 0) { + f = ieee_big_endian_format; + } else { + PyErr_SetString(PyExc_ValueError, "__setformat__() argument 2 must be " + "'unknown', 'IEEE, little-endian' or " + "'IEEE, big-endian'"); + return NULL; + } + + if (f != unknown_format && f != detected) { + PyErr_Format(PyExc_ValueError, "can only set %s format to 'unknown' or the " + "detected platform value", + typestr); + return NULL; + } + + *p = f; + Py_RETURN_NONE; +} + +/*---------------------------------------------------------------------------- + * _PyFloat_{Pack,Unpack}{4,8}. See floatobject.h. + */ +extern "C" int _PyFloat_Pack4(double x, unsigned char* p, int le) noexcept { + if (float_format == unknown_format) { + unsigned char sign; + int e; + double f; + unsigned int fbits; + int incr = 1; + + if (le) { + p += 3; + incr = -1; + } + + if (x < 0) { + sign = 1; + x = -x; + } else + sign = 0; + + f = frexp(x, &e); + + /* Normalize f to be in the range [1.0, 2.0) */ + if (0.5 <= f && f < 1.0) { + f *= 2.0; + e--; + } else if (f == 0.0) + e = 0; + else { + PyErr_SetString(PyExc_SystemError, "frexp() result out of range"); + return -1; + } + + if (e >= 128) + goto Overflow; + else if (e < -126) { + /* Gradual underflow */ + f = ldexp(f, 126 + e); + e = 0; + } else if (!(e == 0 && f == 0.0)) { + e += 127; + f -= 1.0; /* Get rid of leading 1 */ + } + + f *= 8388608.0; /* 2**23 */ + fbits = (unsigned int)(f + 0.5); /* Round */ + assert(fbits <= 8388608); + if (fbits >> 23) { + /* The carry propagated out of a string of 23 1 bits. */ + fbits = 0; + ++e; + if (e >= 255) + goto Overflow; + } + + /* First byte */ + *p = (sign << 7) | (e >> 1); + p += incr; + + /* Second byte */ + *p = (char)(((e & 1) << 7) | (fbits >> 16)); + p += incr; + + /* Third byte */ + *p = (fbits >> 8) & 0xFF; + p += incr; + + /* Fourth byte */ + *p = fbits & 0xFF; + + /* Done */ + return 0; + + } else { + float y = (float)x; + const char* s = (char*)&y; + int i, incr = 1; + + if (Py_IS_INFINITY(y) && !Py_IS_INFINITY(x)) + goto Overflow; + + if ((float_format == ieee_little_endian_format && !le) || (float_format == ieee_big_endian_format && le)) { + p += 3; + incr = -1; + } + + for (i = 0; i < 4; i++) { + *p = *s++; + p += incr; + } + return 0; + } +Overflow: + PyErr_SetString(PyExc_OverflowError, "float too large to pack with f format"); + return -1; +} + +extern "C" int _PyFloat_Pack8(double x, unsigned char* p, int le) noexcept { + if (double_format == unknown_format) { + unsigned char sign; + int e; + double f; + unsigned int fhi, flo; + int incr = 1; + + if (le) { + p += 7; + incr = -1; + } + + if (x < 0) { + sign = 1; + x = -x; + } else + sign = 0; + + f = frexp(x, &e); + + /* Normalize f to be in the range [1.0, 2.0) */ + if (0.5 <= f && f < 1.0) { + f *= 2.0; + e--; + } else if (f == 0.0) + e = 0; + else { + PyErr_SetString(PyExc_SystemError, "frexp() result out of range"); + return -1; + } + + if (e >= 1024) + goto Overflow; + else if (e < -1022) { + /* Gradual underflow */ + f = ldexp(f, 1022 + e); + e = 0; + } else if (!(e == 0 && f == 0.0)) { + e += 1023; + f -= 1.0; /* Get rid of leading 1 */ + } + + /* fhi receives the high 28 bits; flo the low 24 bits (== 52 bits) */ + f *= 268435456.0; /* 2**28 */ + fhi = (unsigned int)f; /* Truncate */ + assert(fhi < 268435456); + + f -= (double)fhi; + f *= 16777216.0; /* 2**24 */ + flo = (unsigned int)(f + 0.5); /* Round */ + assert(flo <= 16777216); + if (flo >> 24) { + /* The carry propagated out of a string of 24 1 bits. */ + flo = 0; + ++fhi; + if (fhi >> 28) { + /* And it also progagated out of the next 28 bits. */ + fhi = 0; + ++e; + if (e >= 2047) + goto Overflow; + } + } + + /* First byte */ + *p = (sign << 7) | (e >> 4); + p += incr; + + /* Second byte */ + *p = (unsigned char)(((e & 0xF) << 4) | (fhi >> 24)); + p += incr; + + /* Third byte */ + *p = (fhi >> 16) & 0xFF; + p += incr; + + /* Fourth byte */ + *p = (fhi >> 8) & 0xFF; + p += incr; + + /* Fifth byte */ + *p = fhi & 0xFF; + p += incr; + + /* Sixth byte */ + *p = (flo >> 16) & 0xFF; + p += incr; + + /* Seventh byte */ + *p = (flo >> 8) & 0xFF; + p += incr; + + /* Eighth byte */ + *p = flo & 0xFF; + /* p += incr; Unneeded (for now) */ + + /* Done */ + return 0; + + Overflow: + PyErr_SetString(PyExc_OverflowError, "float too large to pack with d format"); + return -1; + } else { + const char* s = (char*)&x; + int i, incr = 1; + + if ((double_format == ieee_little_endian_format && !le) || (double_format == ieee_big_endian_format && le)) { + p += 7; + incr = -1; + } + + for (i = 0; i < 8; i++) { + *p = *s++; + p += incr; + } + return 0; + } +} + +extern "C" double _PyFloat_Unpack4(const unsigned char* p, int le) noexcept { + if (float_format == unknown_format) { + unsigned char sign; + int e; + unsigned int f; + double x; + int incr = 1; + + if (le) { + p += 3; + incr = -1; + } + + /* First byte */ + sign = (*p >> 7) & 1; + e = (*p & 0x7F) << 1; + p += incr; + + /* Second byte */ + e |= (*p >> 7) & 1; + f = (*p & 0x7F) << 16; + p += incr; + + if (e == 255) { + PyErr_SetString(PyExc_ValueError, "can't unpack IEEE 754 special value " + "on non-IEEE platform"); + return -1; + } + + /* Third byte */ + f |= *p << 8; + p += incr; + + /* Fourth byte */ + f |= *p; + + x = (double)f / 8388608.0; + + /* XXX This sadly ignores Inf/NaN issues */ + if (e == 0) + e = -126; + else { + x += 1.0; + e -= 127; + } + x = ldexp(x, e); + + if (sign) + x = -x; + + return x; + } else { + float x; + + if ((float_format == ieee_little_endian_format && !le) || (float_format == ieee_big_endian_format && le)) { + char buf[4]; + char* d = &buf[3]; + int i; + + for (i = 0; i < 4; i++) { + *d-- = *p++; + } + memcpy(&x, buf, 4); + } else { + memcpy(&x, p, 4); + } + + return x; + } +} + +extern "C" double _PyFloat_Unpack8(const unsigned char* p, int le) noexcept { + if (double_format == unknown_format) { + unsigned char sign; + int e; + unsigned int fhi, flo; + double x; + int incr = 1; + + if (le) { + p += 7; + incr = -1; + } + + /* First byte */ + sign = (*p >> 7) & 1; + e = (*p & 0x7F) << 4; + + p += incr; + + /* Second byte */ + e |= (*p >> 4) & 0xF; + fhi = (*p & 0xF) << 24; + p += incr; + + if (e == 2047) { + PyErr_SetString(PyExc_ValueError, "can't unpack IEEE 754 special value " + "on non-IEEE platform"); + return -1.0; + } + + /* Third byte */ + fhi |= *p << 16; + p += incr; + + /* Fourth byte */ + fhi |= *p << 8; + p += incr; + + /* Fifth byte */ + fhi |= *p; + p += incr; + + /* Sixth byte */ + flo = *p << 16; + p += incr; + + /* Seventh byte */ + flo |= *p << 8; + p += incr; + + /* Eighth byte */ + flo |= *p; + + x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */ + x /= 268435456.0; /* 2**28 */ + + if (e == 0) + e = -1022; + else { + x += 1.0; + e -= 1023; + } + x = ldexp(x, e); + + if (sign) + x = -x; + + return x; + } else { + double x; + + if ((double_format == ieee_little_endian_format && !le) || (double_format == ieee_big_endian_format && le)) { + char buf[8]; + char* d = &buf[7]; + int i; + + for (i = 0; i < 8; i++) { + *d-- = *p++; + } + memcpy(&x, buf, 8); + } else { + memcpy(&x, p, 8); + } + + return x; + } +} + +#if DBL_MANT_DIG == 53 +#define FIVE_POW_LIMIT 22 +#else +#error "C doubles do not appear to be IEEE 754 binary64 format" +#endif + +extern "C" PyObject* _Py_double_round(double x, int ndigits) noexcept { + + double rounded, m; + Py_ssize_t buflen, mybuflen = 100; + char* buf, *buf_end, shortbuf[100], * mybuf = shortbuf; + int decpt, sign, val, halfway_case; + PyObject* result = NULL; + _Py_SET_53BIT_PRECISION_HEADER; + + /* Easy path for the common case ndigits == 0. */ + if (ndigits == 0) { + rounded = round(x); + if (fabs(rounded - x) == 0.5) + /* halfway between two integers; use round-away-from-zero */ + rounded = x + (x > 0.0 ? 0.5 : -0.5); + return PyFloat_FromDouble(rounded); + } + + /* The basic idea is very simple: convert and round the double to a + decimal string using _Py_dg_dtoa, then convert that decimal string + back to a double with _Py_dg_strtod. There's one minor difficulty: + Python 2.x expects round to do round-half-away-from-zero, while + _Py_dg_dtoa does round-half-to-even. So we need some way to detect + and correct the halfway cases. + + Detection: a halfway value has the form k * 0.5 * 10**-ndigits for + some odd integer k. Or in other words, a rational number x is + exactly halfway between two multiples of 10**-ndigits if its + 2-valuation is exactly -ndigits-1 and its 5-valuation is at least + -ndigits. For ndigits >= 0 the latter condition is automatically + satisfied for a binary float x, since any such float has + nonnegative 5-valuation. For 0 > ndigits >= -22, x needs to be an + integral multiple of 5**-ndigits; we can check this using fmod. + For -22 > ndigits, there are no halfway cases: 5**23 takes 54 bits + to represent exactly, so any odd multiple of 0.5 * 10**n for n >= + 23 takes at least 54 bits of precision to represent exactly. + + Correction: a simple strategy for dealing with halfway cases is to + (for the halfway cases only) call _Py_dg_dtoa with an argument of + ndigits+1 instead of ndigits (thus doing an exact conversion to + decimal), round the resulting string manually, and then convert + back using _Py_dg_strtod. + */ + + /* nans, infinities and zeros should have already been dealt + with by the caller (in this case, builtin_round) */ + assert(std::isfinite(x) && x != 0.0); + + /* find 2-valuation val of x */ + m = frexp(x, &val); + while (m != floor(m)) { + m *= 2.0; + val--; + } + + /* determine whether this is a halfway case */ + if (val == -ndigits - 1) { + if (ndigits >= 0) + halfway_case = 1; + else if (ndigits >= -FIVE_POW_LIMIT) { + double five_pow = 1.0; + int i; + for (i = 0; i < -ndigits; i++) + five_pow *= 5.0; + halfway_case = fmod(x, five_pow) == 0.0; + } else + halfway_case = 0; + } else + halfway_case = 0; + + /* round to a decimal string; use an extra place for halfway case */ + _Py_SET_53BIT_PRECISION_START; + buf = _Py_dg_dtoa(x, 3, ndigits + halfway_case, &decpt, &sign, &buf_end); + _Py_SET_53BIT_PRECISION_END; + if (buf == NULL) { + PyErr_NoMemory(); + return NULL; + } + buflen = buf_end - buf; + + /* in halfway case, do the round-half-away-from-zero manually */ + if (halfway_case) { + int i, carry; + /* sanity check: _Py_dg_dtoa should not have stripped + any zeros from the result: there should be exactly + ndigits+1 places following the decimal point, and + the last digit in the buffer should be a '5'.*/ + assert(buflen - decpt == ndigits + 1); + assert(buf[buflen - 1] == '5'); + + /* increment and shift right at the same time. */ + decpt += 1; + carry = 1; + for (i = buflen - 1; i-- > 0;) { + carry += buf[i] - '0'; + buf[i + 1] = carry % 10 + '0'; + carry /= 10; + } + buf[0] = carry + '0'; + } + + /* Get new buffer if shortbuf is too small. Space needed <= buf_end - + buf + 8: (1 extra for '0', 1 for sign, 5 for exp, 1 for '\0'). */ + if (buflen + 8 > mybuflen) { + mybuflen = buflen + 8; + mybuf = (char*)PyMem_Malloc(mybuflen); + if (mybuf == NULL) { + PyErr_NoMemory(); + goto exit; + } + } + /* copy buf to mybuf, adding exponent, sign and leading 0 */ + PyOS_snprintf(mybuf, mybuflen, "%s0%se%d", (sign ? "-" : ""), buf, decpt - (int)buflen); + + /* and convert the resulting string back to a double */ + errno = 0; + _Py_SET_53BIT_PRECISION_START; + rounded = _Py_dg_strtod(mybuf, NULL); + _Py_SET_53BIT_PRECISION_END; + if (errno == ERANGE && fabs(rounded) >= 1.) + PyErr_SetString(PyExc_OverflowError, "rounded value too large to represent"); + else + result = PyFloat_FromDouble(rounded); + + /* done computing value; now clean up */ + if (mybuf != shortbuf) + PyMem_Free(mybuf); +exit: + _Py_dg_freedtoa(buf); + return result; +} + void setupFloat() { _addFunc("__add__", BOXED_FLOAT, (void*)floatAddFloat, (void*)floatAddInt, (void*)floatAdd); float_cls->giveAttr("__radd__", float_cls->getattr("__add__")); -- 2.30.9