diff options
Diffstat (limited to 'Objects/floatobject.c')
-rw-r--r-- | Objects/floatobject.c | 407 |
1 files changed, 407 insertions, 0 deletions
diff --git a/Objects/floatobject.c b/Objects/floatobject.c index 45cb905..f70771e 100644 --- a/Objects/floatobject.c +++ b/Objects/floatobject.c @@ -10,6 +10,11 @@ #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> #endif @@ -1109,6 +1114,404 @@ float_float(PyObject *v) return v; } +/* turn ASCII hex characters into integer values and vice versa */ + +static char +char_from_hex(int x) +{ + assert(0 <= x && x < 16); + return "0123456789abcdef"[x]; +} + +static int +hex_from_char(char c) { + int x; + assert(isxdigit(c)); + switch(c) { + case '0': + x = 0; + break; + case '1': + x = 1; + break; + case '2': + x = 2; + break; + case '3': + x = 3; + break; + case '4': + x = 4; + break; + case '5': + x = 5; + break; + case '6': + x = 6; + break; + case '7': + x = 7; + break; + case '8': + x = 8; + break; + case '9': + x = 9; + break; + case 'a': + case 'A': + x = 10; + break; + case 'b': + case 'B': + x = 11; + break; + case 'c': + case 'C': + x = 12; + break; + case 'd': + case 'D': + x = 13; + break; + case 'e': + case 'E': + x = 14; + break; + case 'f': + case 'F': + x = 15; + break; + default: + x = -1; + break; + } + return x; +} + +/* convert a float to a hexadecimal string */ + +/* TOHEX_NBITS is DBL_MANT_DIG rounded up to the next integer + of the form 4k+1. */ +#define TOHEX_NBITS DBL_MANT_DIG + 3 - (DBL_MANT_DIG+2)%4 + +static PyObject * +float_hex(PyObject *v) +{ + double x, m; + int e, shift, i, si, esign; + /* Space for 1+(TOHEX_NBITS-1)/4 digits, a decimal point, and the + trailing NUL byte. */ + char s[(TOHEX_NBITS-1)/4+3]; + + CONVERT_TO_DOUBLE(v, x); + + if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) + return float_str((PyFloatObject *)v); + + if (x == 0.0) { + if(copysign(1.0, x) == -1.0) + return PyString_FromString("-0x0.0p+0"); + else + return PyString_FromString("0x0.0p+0"); + } + + m = frexp(fabs(x), &e); + shift = 1 - MAX(DBL_MIN_EXP - e, 0); + m = ldexp(m, shift); + e -= shift; + + si = 0; + s[si] = char_from_hex((int)m); + si++; + m -= (int)m; + s[si] = '.'; + si++; + for (i=0; i < (TOHEX_NBITS-1)/4; i++) { + m *= 16.0; + s[si] = char_from_hex((int)m); + si++; + m -= (int)m; + } + s[si] = '\0'; + + if (e < 0) { + esign = (int)'-'; + e = -e; + } + else + esign = (int)'+'; + + if (x < 0.0) + return PyString_FromFormat("-0x%sp%c%d", s, esign, e); + else + return PyString_FromFormat("0x%sp%c%d", s, esign, e); +} + +PyDoc_STRVAR(float_hex_doc, +"float.hex() -> string\n\ +\n\ +Return a hexadecimal representation of a floating-point number.\n\ +>>> (-0.1).hex()\n\ +'-0x1.999999999999ap-4'\n\ +>>> 3.14159.hex()\n\ +'0x1.921f9f01b866ep+1'"); + +/* Convert a hexadecimal string to a float. */ + +static PyObject * +float_fromhex(PyObject *cls, PyObject *arg) +{ + PyObject *result_as_float, *result; + double x; + long exp, top_exp, lsb, key_digit; + char *s, *coeff_start, *s_store, *coeff_end, *exp_start, *s_end; + int half_eps, digit, round_up, sign=1; + Py_ssize_t length, ndigits, fdigits, i; + + /* + * For the sake of simplicity and correctness, we impose an artificial + * limit on ndigits, the total number of hex digits in the coefficient + * The limit is chosen to ensure that, writing exp for the exponent, + * + * (1) if exp > LONG_MAX/2 then the value of the hex string is + * guaranteed to overflow (provided it's nonzero) + * + * (2) if exp < LONG_MIN/2 then the value of the hex string is + * guaranteed to underflow to 0. + * + * (3) if LONG_MIN/2 <= exp <= LONG_MAX/2 then there's no danger of + * overflow in the calculation of exp and top_exp below. + * + * More specifically, ndigits is assumed to satisfy the following + * inequalities: + * + * 4*ndigits <= DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2 + * 4*ndigits <= LONG_MAX/2 + 1 - DBL_MAX_EXP + * + * If either of these inequalities is not satisfied, a ValueError is + * raised. Otherwise, write x for the value of the hex string, and + * assume x is nonzero. Then + * + * 2**(exp-4*ndigits) <= |x| < 2**(exp+4*ndigits). + * + * Now if exp > LONG_MAX/2 then: + * + * exp - 4*ndigits >= LONG_MAX/2 + 1 - (LONG_MAX/2 + 1 - DBL_MAX_EXP) + * = DBL_MAX_EXP + * + * so |x| >= 2**DBL_MAX_EXP, which is too large to be stored in C + * double, so overflows. If exp < LONG_MIN/2, then + * + * exp + 4*ndigits <= LONG_MIN/2 - 1 + ( + * DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2) + * = DBL_MIN_EXP - DBL_MANT_DIG - 1 + * + * and so |x| < 2**(DBL_MIN_EXP-DBL_MANT_DIG-1), hence underflows to 0 + * when converted to a C double. + * + * It's easy to show that if LONG_MIN/2 <= exp <= LONG_MAX/2 then both + * exp+4*ndigits and exp-4*ndigits are within the range of a long. + */ + + if (PyString_AsStringAndSize(arg, &s, &length)) + return NULL; + s_end = s + length; + + /******************** + * Parse the string * + ********************/ + + /* leading whitespace and optional sign */ + while (isspace(*s)) + s++; + if (*s == '-') { + s++; + sign = -1; + } + else if (*s == '+') + s++; + + /* infinities and nans */ + if (PyOS_mystrnicmp(s, "nan", 4) == 0) { + x = Py_NAN; + goto finished; + } + if (PyOS_mystrnicmp(s, "inf", 4) == 0 || + PyOS_mystrnicmp(s, "infinity", 9) == 0) { + x = sign*Py_HUGE_VAL; + goto finished; + } + + /* [0x] */ + s_store = s; + if (*s == '0') { + s++; + if (tolower(*s) == (int)'x') + s++; + else + s = s_store; + } + + /* coefficient: <integer> [. <fraction>] */ + coeff_start = s; + while (isxdigit(*s)) + s++; + s_store = s; + if (*s == '.') { + s++; + while (isxdigit(*s)) + s++; + coeff_end = s-1; + } + else + coeff_end = s; + + /* ndigits = total # of hex digits; fdigits = # after point */ + ndigits = coeff_end - coeff_start; + fdigits = coeff_end - s_store; + if (ndigits == 0) + goto parse_error; + if (ndigits > MIN(DBL_MIN_EXP - DBL_MANT_DIG - LONG_MIN/2, + LONG_MAX/2 + 1 - DBL_MAX_EXP)/4) + goto insane_length_error; + + /* [p <exponent>] */ + if (tolower(*s) == (int)'p') { + s++; + exp_start = s; + if (*s == '-' || *s == '+') + s++; + if (!isdigit(*s)) + goto parse_error; + s++; + while (isdigit(*s)) + s++; + exp = strtol(exp_start, NULL, 10); + } + else + exp = 0; + + /* optional trailing whitespace leading to the end of the string */ + while (isspace(*s)) + s++; + if (s != s_end) + goto parse_error; + +/* for 0 <= j < ndigits, HEX_DIGIT(j) gives the jth most significant digit */ +#define HEX_DIGIT(j) hex_from_char(*((j) < fdigits ? \ + coeff_end-(j) : \ + coeff_end-1-(j))) + + /******************************************* + * Compute rounded value of the hex string * + *******************************************/ + + /* Discard leading zeros, and catch extreme overflow and underflow */ + while (ndigits > 0 && HEX_DIGIT(ndigits-1) == 0) + ndigits--; + if (ndigits == 0 || exp < LONG_MIN/2) { + x = sign * 0.0; + goto finished; + } + if (exp > LONG_MAX/2) + goto overflow_error; + + /* Adjust exponent for fractional part. */ + exp = exp - 4*((long)fdigits); + + /* top_exp = 1 more than exponent of most sig. bit of coefficient */ + top_exp = exp + 4*((long)ndigits - 1); + for (digit = HEX_DIGIT(ndigits-1); digit != 0; digit /= 2) + top_exp++; + + /* catch almost all nonextreme cases of overflow and underflow here */ + if (top_exp < DBL_MIN_EXP - DBL_MANT_DIG) { + x = sign * 0.0; + goto finished; + } + if (top_exp > DBL_MAX_EXP) + goto overflow_error; + + /* lsb = exponent of least significant bit of the *rounded* value. + This is top_exp - DBL_MANT_DIG unless result is subnormal. */ + lsb = MAX(top_exp, (long)DBL_MIN_EXP) - DBL_MANT_DIG; + + x = 0.0; + if (exp >= lsb) { + /* no rounding required */ + for (i = ndigits-1; i >= 0; i--) + x = 16.0*x + HEX_DIGIT(i); + x = sign * ldexp(x, (int)(exp)); + goto finished; + } + /* rounding required. key_digit is the index of the hex digit + containing the first bit to be rounded away. */ + half_eps = 1 << (int)((lsb - exp - 1) % 4); + key_digit = (lsb - exp - 1) / 4; + for (i = ndigits-1; i > key_digit; i--) + x = 16.0*x + HEX_DIGIT(i); + digit = HEX_DIGIT(key_digit); + x = 16.0*x + (double)(digit & (16-2*half_eps)); + + /* round-half-even: round up if bit lsb-1 is 1 and at least one of + bits lsb, lsb-2, lsb-3, lsb-4, ... is 1. */ + if ((digit & half_eps) != 0) { + round_up = 0; + if ((digit & (3*half_eps-1)) != 0 || + (half_eps == 8 && (HEX_DIGIT(key_digit+1) & 1) != 0)) + round_up = 1; + else + for (i = key_digit-1; i >= 0; i--) + if (HEX_DIGIT(i) != 0) { + round_up = 1; + break; + } + if (round_up == 1) { + x += 2*half_eps; + if (top_exp == DBL_MAX_EXP && + x == ldexp((double)(2*half_eps), DBL_MANT_DIG)) + /* overflow corner case: pre-rounded value < + 2**DBL_MAX_EXP; rounded=2**DBL_MAX_EXP. */ + goto overflow_error; + } + } + x = sign * ldexp(x, (int)(exp+4*key_digit)); + + finished: + result_as_float = Py_BuildValue("(d)", x); + if (result_as_float == NULL) + return NULL; + result = PyObject_CallObject(cls, result_as_float); + Py_DECREF(result_as_float); + return result; + + overflow_error: + PyErr_SetString(PyExc_OverflowError, + "hexadecimal value too large to represent as a float"); + return NULL; + + parse_error: + PyErr_SetString(PyExc_ValueError, + "invalid hexadecimal floating-point string"); + return NULL; + + insane_length_error: + PyErr_SetString(PyExc_ValueError, + "hexadecimal string too long to convert"); + return NULL; +} + +PyDoc_STRVAR(float_fromhex_doc, +"float.fromhex(string) -> float\n\ +\n\ +Create a floating-point number from a hexadecimal string.\n\ +>>> float.fromhex('0x1.ffffp10')\n\ +2047.984375\n\ +>>> float.fromhex('-0x1p-1074')\n\ +-4.9406564584124654e-324"); + + static PyObject * float_as_integer_ratio(PyObject *v, PyObject *unused) { @@ -1433,6 +1836,10 @@ static PyMethodDef float_methods[] = { "Returns the Integral closest to x between 0 and x."}, {"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS, float_as_integer_ratio_doc}, + {"fromhex", (PyCFunction)float_fromhex, + METH_O|METH_CLASS, float_fromhex_doc}, + {"hex", (PyCFunction)float_hex, + METH_NOARGS, float_hex_doc}, {"is_integer", (PyCFunction)float_is_integer, METH_NOARGS, "Returns True if the float is an integer."}, #if 0 |