Merged revisions 64974 via svnmerge from

svn+ssh://pythondev@svn.python.org/python/trunk

........
  r64974 | mark.dickinson | 2008-07-15 20:08:33 +0100 (Tue, 15 Jul 2008) | 3 lines

  Issue #3008: add instance method float.hex and class method float.fromhex
  to convert floats to and from hexadecimal strings respectively.
........
This commit is contained in:
Mark Dickinson 2008-07-16 11:30:51 +00:00
parent 0c474d01a1
commit 65fe25e597
5 changed files with 867 additions and 1 deletions

View file

@ -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
@ -1015,6 +1020,405 @@ 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 PyUnicode_FromString("-0x0.0p+0");
else
return PyUnicode_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 PyUnicode_FromFormat("-0x%sp%c%d", s, esign, e);
else
return PyUnicode_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.
*/
s = PyUnicode_AsStringAndSize(arg, &length);
if (s == NULL)
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)
{
@ -1326,6 +1730,10 @@ static PyMethodDef float_methods[] = {
"When an argument is passed, works like built-in round(x, ndigits)."},
{"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