diff -r 25a2aeecb34b Include/floatobject.h --- a/Include/floatobject.h Thu Mar 31 02:05:54 2011 +0200 +++ b/Include/floatobject.h Fri Apr 01 11:56:56 2011 -0700 @@ -74,9 +74,9 @@ * happens in such cases is partly accidental (alas). */ -/* The pack routines write 4 or 8 bytes, starting at p. le is a bool +/* The pack routines write 2, 4 or 8 bytes, starting at p. le is a bool * argument, true if you want the string in little-endian format (exponent - * last, at p+3 or p+7), false if you want big-endian format (exponent + * last, at p+1, p+3 or p+7), false if you want big-endian format (exponent * first, at p). * Return value: 0 if all is OK, -1 if error (and an exception is * set, most likely OverflowError). @@ -84,6 +84,7 @@ * 1): What this does is undefined if x is a NaN or infinity. * 2): -0.0 and +0.0 produce the same string. */ +PyAPI_FUNC(int) _PyFloat_Pack2(double x, unsigned char *p, int le); PyAPI_FUNC(int) _PyFloat_Pack4(double x, unsigned char *p, int le); PyAPI_FUNC(int) _PyFloat_Pack8(double x, unsigned char *p, int le); @@ -96,14 +97,15 @@ PyAPI_FUNC(int) _PyFloat_Digits(char *buf, double v, int *signum); PyAPI_FUNC(void) _PyFloat_DigitsInit(void); -/* The unpack routines read 4 or 8 bytes, starting at p. le is a bool +/* The unpack routines read 2, 4 or 8 bytes, starting at p. le is a bool * argument, true if the string is in little-endian format (exponent - * last, at p+3 or p+7), false if big-endian (exponent first, at p). + * last, at p+1, p+3 or p+7), false if big-endian (exponent first, at p). * Return value: The unpacked double. On error, this is -1.0 and * PyErr_Occurred() is true (and an exception is set, most likely * OverflowError). Note that on a non-IEEE platform this will refuse * to unpack a string that represents a NaN or infinity. */ +PyAPI_FUNC(double) _PyFloat_Unpack2(const unsigned char *p, int le); PyAPI_FUNC(double) _PyFloat_Unpack4(const unsigned char *p, int le); PyAPI_FUNC(double) _PyFloat_Unpack8(const unsigned char *p, int le); diff -r 25a2aeecb34b Lib/test/test_struct.py --- a/Lib/test/test_struct.py Thu Mar 31 02:05:54 2011 +0200 +++ b/Lib/test/test_struct.py Fri Apr 01 11:56:56 2011 -0700 @@ -556,6 +556,58 @@ s = struct.Struct('i') s.__init__('ii') + def test_half_float(self): + import math + + # http://en.wikipedia.org/wiki/Half_precision_floating-point_format + bitsValue_list = [ + ('>e', b'\x3c\x00', 1.0), + ('>e', b'\xc0\x00', -2.0), + ('>e', b'\x7b\xff', 6.5504 * 10**4), # (max half precision) + ('>e', b'\x04\x00', 2**-14), # ~= 6.10352 * 10**-5 (minimum positive normal) + ('>e', b'\x00\x01', 2**-24), # ~= 5.96046 * 10**-8 (minimum strictly positive subnormal) + ('>e', b'\x00\x00', 0.0), + ('>e', b'\x80\x00', -0.0), + ('>e', b'\x7c\x00', float('+inf')), + ('>e', b'\xfc\x00', float('-inf')), + ('>e', b'\x35\x55', 0.333251953125), # ~= 1/3 + + ('e', b'\xfc\x01'), + ] + + for formatcode, bits in bitsValue_list: + self.assertTrue(math.isnan(struct.unpack(formatcode, memoryview(bits))[0])) + + # + # + #data1 = array.array('B', b'\x12\x34\x56\x78') + #data2 = memoryview(b'\x12\x34\x56\x78') # XXX b'......XXXX......', 6, 4 + #for data in [data1, data2]: + # value, = struct.unpack('>e', data) + # self.assertEqual(value, 0x12345678) + + + def test_main(): run_unittest(StructTest) diff -r 25a2aeecb34b Modules/_struct.c --- a/Modules/_struct.c Thu Mar 31 02:05:54 2011 +0200 +++ b/Modules/_struct.c Fri Apr 01 11:56:56 2011 -0700 @@ -218,6 +218,18 @@ /* Floating point helpers */ +static PyObject * +unpack_halffloat(const char *p, /* start of 4-byte string */ + int le) /* true for little-endian, false for big-endian */ +{ + double x; + + x = _PyFloat_Unpack2((unsigned char *)p, le); + if (x == -1.0 && PyErr_Occurred()) + return NULL; + return PyFloat_FromDouble(x); +} + static PyObject * unpack_float(const char *p, /* start of 4-byte string */ @@ -406,6 +418,13 @@ static PyObject * +nu_halffloat(const char *p, const formatdef *f) +{ + int le = 1; + return PyFloat_FromDouble(_PyFloat_Unpack2(p, (int)*(char*)&le)); +} + +static PyObject * nu_float(const char *p, const formatdef *f) { float x; @@ -596,6 +615,20 @@ } static int +np_halffloat(char *p, PyObject *v, const formatdef *f) +{ + double x = PyFloat_AsDouble(v); + if (x == -1 && PyErr_Occurred()) { + PyErr_SetString(StructError, + "required argument is not a float"); + return -1; + } + int le = 1; + _PyFloat_Pack2(x, p, (int)*(char*)&le); + return 0; +} + +static int np_float(char *p, PyObject *v, const formatdef *f) { float x = (float)PyFloat_AsDouble(v); @@ -656,6 +689,7 @@ {'Q', sizeof(PY_LONG_LONG), LONG_LONG_ALIGN, nu_ulonglong,np_ulonglong}, #endif {'?', sizeof(BOOL_TYPE), BOOL_ALIGN, nu_bool, np_bool}, + {'e', sizeof(short), FLOAT_ALIGN, nu_halffloat, np_halffloat}, {'f', sizeof(float), FLOAT_ALIGN, nu_float, np_float}, {'d', sizeof(double), DOUBLE_ALIGN, nu_double, np_double}, {'P', sizeof(void *), VOID_P_ALIGN, nu_void_p, np_void_p}, @@ -739,6 +773,12 @@ } static PyObject * +bu_halffloat(const char *p, const formatdef *f) +{ + return unpack_halffloat(p, 0); +} + +static PyObject * bu_float(const char *p, const formatdef *f) { return unpack_float(p, 0); @@ -835,6 +875,18 @@ } static int +bp_halffloat(char *p, PyObject *v, const formatdef *f) +{ + double x = PyFloat_AsDouble(v); + if (x == -1 && PyErr_Occurred()) { + PyErr_SetString(StructError, + "required argument is not a float"); + return -1; + } + return _PyFloat_Pack2(x, (unsigned char *)p, 0); +} + +static int bp_float(char *p, PyObject *v, const formatdef *f) { double x = PyFloat_AsDouble(v); @@ -885,6 +937,7 @@ {'q', 8, 0, bu_longlong, bp_longlong}, {'Q', 8, 0, bu_ulonglong, bp_ulonglong}, {'?', 1, 0, bu_bool, bp_bool}, + {'e', 2, 0, bu_halffloat, bp_halffloat}, {'f', 4, 0, bu_float, bp_float}, {'d', 8, 0, bu_double, bp_double}, {0} @@ -967,6 +1020,12 @@ } static PyObject * +lu_halffloat(const char *p, const formatdef *f) +{ + return unpack_halffloat(p, 1); +} + +static PyObject * lu_float(const char *p, const formatdef *f) { return unpack_float(p, 1); @@ -1055,6 +1114,18 @@ } static int +lp_halffloat(char *p, PyObject *v, const formatdef *f) +{ + double x = PyFloat_AsDouble(v); + if (x == -1 && PyErr_Occurred()) { + PyErr_SetString(StructError, + "required argument is not a float"); + return -1; + } + return _PyFloat_Pack2(x, (unsigned char *)p, 1); +} + +static int lp_float(char *p, PyObject *v, const formatdef *f) { double x = PyFloat_AsDouble(v); @@ -1095,6 +1166,7 @@ {'Q', 8, 0, lu_ulonglong, lp_ulonglong}, {'?', 1, 0, bu_bool, bp_bool}, /* Std rep not endian dep, but potentially different from native rep -- reuse bx_bool funcs. */ + {'e', 2, 0, lu_halffloat, lp_halffloat}, {'f', 4, 0, lu_float, lp_float}, {'d', 8, 0, lu_double, lp_double}, {0} @@ -1948,7 +2020,7 @@ x: pad byte (no data); c:char; b:signed byte; B:unsigned byte;\n\ ?: _Bool (requires C99; if not available, char is used instead)\n\ h:short; H:unsigned short; i:int; I:unsigned int;\n\ - l:long; L:unsigned long; f:float; d:double.\n\ + l:long; L:unsigned long; f:float; d:double; e:half-float.\n\ Special cases (preceding decimal count indicates length):\n\ s:string (array of char); p: pascal string (with count byte).\n\ Special case (only available in native format):\n\ diff -r 25a2aeecb34b Objects/floatobject.c --- a/Objects/floatobject.c Thu Mar 31 02:05:54 2011 +0200 +++ b/Objects/floatobject.c Fri Apr 01 11:56:56 2011 -0700 @@ -2032,8 +2032,272 @@ } /*---------------------------------------------------------------------------- - * _PyFloat_{Pack,Unpack}{4,8}. See floatobject.h. + * _PyFloat_{Pack,Unpack}{2,4,8}. See floatobject.h. + * {Pack,Unpack}2 adapted from + * https://github.com/numpy/numpy/blob/master/numpy/core/src/npymath/halffloat.c + * as if NPY_HALF_ROUND_TIES_TO_EVEN and NPY_HALF_GENERATE_OVERFLOW were set to + * 1, but NPY_HALF_GENERATE_UNDERFLOW set to 0. */ + +int +_PyFloat_Pack2(double x, unsigned char *p, int le) +{ + int incr = 1; + if (le) { + p += 1; + incr = -1; + } + + typedef unsigned short npy_uint16; + npy_uint16 h_sgn, h_exp, h_sig, bits; + +#if SIZEOF_LONG == 8 || SIZEOF_LONG_LONG == 8 +#if SIZEOF_LONG == 8 + typedef unsigned long npy_uint64; +#else +#if SIZEOF_LONG_LONG == 8 + typedef unsigned long long npy_uint64; +#endif +#endif + + union { double d; npy_uint64 i; } di; + + di.d = x; + npy_uint64 d = di.i; + npy_uint64 d_exp, d_sig; + + h_sgn = (d&0x8000000000000000ULL) >> 48; + d_exp = (d&0x7ff0000000000000ULL); + + /* Exponent overflow/NaN converts to signed inf/NaN */ + if (d_exp >= 0x40f0000000000000ULL) { + if (d_exp == 0x7ff0000000000000ULL) { + /* Inf or NaN */ + d_sig = (d&0x000fffffffffffffULL); + if (d_sig != 0) { + /* NaN - propagate the flag in the significand... */ + npy_uint16 ret = (npy_uint16) (0x7c00u + (d_sig >> 42)); + /* ...but make sure it stays a NaN */ + if (ret == 0x7c00u) { + ret++; + } + bits = h_sgn + ret; + } else { + /* signed inf */ + bits = h_sgn + 0x7c00u; + } + } else { + /* overflow to signed inf */ + goto Overflow; + bits = h_sgn + 0x7c00u; + } + } + + + /* Exponent underflow converts to subnormal half or signed zero */ + else if (d_exp <= 0x3f00000000000000ULL) { + /* + * Signed zeros, subnormal floats, and floats with small + * exponents all convert to signed zero halfs. + */ + if (d_exp < 0x3e60000000000000ULL) { +//#if NPY_HALF_GENERATE_UNDERFLOW +// /* If d != 0, it underflowed to 0 */ +// if ((d&0x7fffffffffffffffULL) != 0) { +// npy_set_floatstatus_underflow(); +// } +//#endif + bits = h_sgn; + } + else { + /* Make the subnormal significand */ + d_exp >>= 52; + d_sig = (0x0010000000000000ULL + (d&0x000fffffffffffffULL)); + //#if NPY_HALF_GENERATE_UNDERFLOW + // /* If it's not exactly represented, it underflowed */ + // if ((d_sig&(((npy_uint64)1 << (1051 - d_exp)) - 1)) != 0) { + // npy_set_floatstatus_underflow(); + // } + //#endif + d_sig >>= (1009 - d_exp); + /* Handle rounding by adding 1 to the bit beyond half precision */ + + /* + * If the last bit in the half significand is 0 (already even), and + * the remaining bit pattern is 1000...0, then we do not add one + * to the bit after the half significand. In all other cases, we do. + */ + if ((d_sig&0x000007ffffffffffULL) != 0x0000020000000000ULL) { + d_sig += 0x0000020000000000ULL; + } + h_sig = (npy_uint16) (d_sig >> 42); + /* + * If the rounding causes a bit to spill into h_exp, it will + * increment h_exp from zero to one and h_sig will be zero. + * This is the correct result. + */ + bits = h_sgn + h_sig; + } + } + else { + /* Regular case with no overflow or underflow */ + h_exp = (npy_uint16) ((d_exp - 0x3f00000000000000ULL) >> 42); + /* Handle rounding by adding 1 to the bit beyond half precision */ + d_sig = (d&0x000fffffffffffffULL); + + /* + * If the last bit in the half significand is 0 (already even), and + * the remaining bit pattern is 1000...0, then we do not add one + * to the bit after the half significand. In all other cases, we do. + */ + if ((d_sig&0x000007ffffffffffULL) != 0x0000020000000000ULL) { + d_sig += 0x0000020000000000ULL; + } + h_sig = (npy_uint16) (d_sig >> 42); + + /* + * If the rounding causes a bit to spill into h_exp, it will + * increment h_exp by one and h_sig will be zero. This is the + * correct result. h_exp may increment to 15, at greatest, in + * which case the result overflows to a signed inf. + */ + h_sig += h_exp; + if (h_sig == 0x7c00u) { + goto Overflow; + } + bits = h_sgn + h_sig; + } + + +#else + typedef unsigned int npy_uint32; + + union { float f; npy_uint32 i; } fi; + + fi.f = x; + npy_uint32 f = fi.i; + npy_uint32 f_exp, f_sig; + + h_sgn = (npy_uint16) ((f&0x80000000u) >> 16); + f_exp = (f&0x7f800000u); + + /* Exponent overflow/NaN converts to signed inf/NaN */ + if (f_exp >= 0x47800000u) { + if (f_exp == 0x7f800000u) { + /* Inf or NaN */ + f_sig = (f&0x007fffffu); + if (f_sig != 0) { + /* NaN - propagate the flag in the significand... */ + npy_uint16 ret = (npy_uint16) (0x7c00u + (f_sig >> 13)); + /* ...but make sure it stays a NaN */ + if (ret == 0x7c00u) { + ret++; + } + bits = h_sgn + ret; + } else { + /* signed inf */ + bits = (npy_uint16) (h_sgn + 0x7c00u); + } + } else { + goto Overflow; + /* overflow to signed inf */ + bits = (npy_uint16) (h_sgn + 0x7c00u); + } + } + + /* Exponent underflow converts to a subnormal half or signed zero */ + else if (f_exp <= 0x38000000u) { + /* + * Signed zeros, subnormal floats, and floats with small + * exponents all convert to signed zero halfs. + */ + if (f_exp < 0x33000000u) { +//#if NPY_HALF_GENERATE_UNDERFLOW +// /* If f != 0, it underflowed to 0 */ +// if ((f&0x7fffffff) != 0) { +// npy_set_floatstatus_underflow(); +// } +//#endif + bits = h_sgn; + } + else { + /* Make the subnormal significand */ + f_exp >>= 23; + f_sig = (0x00800000u + (f&0x007fffffu)); + //#if NPY_HALF_GENERATE_UNDERFLOW + // /* If it's not exactly represented, it underflowed */ + // if ((f_sig&(((npy_uint32)1 << (126 - f_exp)) - 1)) != 0) { + // npy_set_floatstatus_underflow(); + // } + //#endif + f_sig >>= (113 - f_exp); + /* Handle rounding by adding 1 to the bit beyond half precision */ + + /* + * If the last bit in the half significand is 0 (already even), and + * the remaining bit pattern is 1000...0, then we do not add one + * to the bit after the half significand. In all other cases, we do. + */ + if ((f_sig&0x00003fffu) != 0x00001000u) { + f_sig += 0x00001000u; + } + + h_sig = (npy_uint16) (f_sig >> 13); + /* + * If the rounding causes a bit to spill into h_exp, it will + * increment h_exp from zero to one and h_sig will be zero. + * This is the correct result. + */ + bits = (npy_uint16) (h_sgn + h_sig); + } + } + else { + + /* Regular case with no overflow or underflow */ + h_exp = (npy_uint16) ((f_exp - 0x38000000u) >> 13); + /* Handle rounding by adding 1 to the bit beyond half precision */ + f_sig = (f&0x007fffffu); + + /* + * If the last bit in the half significand is 0 (already even), and + * the remaining bit pattern is 1000...0, then we do not add one + * to the bit after the half significand. In all other cases, we do. + */ + if ((f_sig&0x00003fffu) != 0x00001000u) { + f_sig += 0x00001000u; + } + h_sig = (npy_uint16) (f_sig >> 13); + /* + * If the rounding causes a bit to spill into h_exp, it will + * increment h_exp by one and h_sig will be zero. This is the + * correct result. h_exp may increment to 15, at greatest, in + * which case the result overflows to a signed inf. + */ + + h_sig += h_exp; + if (h_sig == 0x7c00u) { + //npy_set_floatstatus_overflow(); + goto Overflow; + } + bits = h_sgn + h_sig; + } +#endif + + /* First byte */ + *p = (bits >> 8) & 0xff; + p += incr; + + /* Second byte */ + *p = bits & 0xFF; + + return 0; + + Overflow: + PyErr_SetString(PyExc_OverflowError, + "float too large to pack with e format"); + return -1; +} + int _PyFloat_Pack4(double x, unsigned char *p, int le) { @@ -2269,6 +2533,61 @@ } double +_PyFloat_Unpack2(const unsigned char *p, int le) +{ + typedef unsigned int npy_uint32; + typedef unsigned short npy_uint16; + + union { float f; int i; } fi; + + int incr = 1; + if (le) { + p += 1; + incr = -1; + } + + npy_uint16 h = *p; + h <<= 8; + p += incr; + h += *p; + + npy_uint16 h_exp, h_sig; + npy_uint32 f_sgn, f_exp, f_sig; + + h_exp = (h&0x7c00u); + f_sgn = ((npy_uint32)h&0x8000u) << 16; + switch (h_exp) { + case 0x0000u: /* 0 or subnormal */ + h_sig = (h&0x03ffu); + /* Signed zero */ + if (h_sig == 0) { + fi.i = f_sgn; + break; + } + /* Subnormal */ + h_sig <<= 1; + while ((h_sig&0x0400u) == 0) { + h_sig <<= 1; + h_exp++; + } + f_exp = ((npy_uint32)(127 - 15 - h_exp)) << 23; + f_sig = ((npy_uint32)(h_sig&0x03ffu)) << 13; + fi.i = f_sgn + f_exp + f_sig; + break; + case 0x7c00u: /* inf or NaN */ + /* All-ones exponent and a copy of the significand */ + fi.i = f_sgn + 0x7f800000u + (((npy_uint32)(h&0x03ffu)) << 13); + break; + default: /* normalized */ + /* Just need to adjust the exponent and shift */ + fi.i = f_sgn + (((npy_uint32)(h&0x7fffu) + 0x1c000u) << 13); + break; + } + + return fi.f; +} + +double _PyFloat_Unpack4(const unsigned char *p, int le) { if (float_format == unknown_format) {