New private API functions _PyFloat_{Pack,Unpack}(4,8}. This is a

refactoring to get all the duplicates of this delicate code out of the
cPickle and struct modules.

Include/floatobject.h

@@ -47,6 +47,48 @@ PyAPI_FUNC(void) PyFloat_AsReprString(char*, PyFloatObject *v);
    preserve precision across conversions. */
 PyAPI_FUNC(void) PyFloat_AsString(char*, PyFloatObject *v);
 
+/* _PyFloat_{Pack,Unpack}{4,8}
+ *
+ * The struct and pickle (at least) modules need an efficient platform-
+ * independent way to store floating-point values as byte strings.
+ * The Pack routines produce a string from a C double, and the Unpack
+ * routines produce a C double from such a string.  The suffix (4 or 8)
+ * specifies the number of bytes in the string.
+ *
+ * Excepting NaNs and infinities (which aren't handled correctly), the 4-
+ * byte format is identical to the IEEE-754 single precision format, and
+ * the 8-byte format to the IEEE-754 double precision format.  On non-
+ * IEEE platforms with more precision, or larger dynamic range, than
+ * 754 supports, not all values can be packed; on non-IEEE platforms with
+ * less precision, or smaller dynamic range, not all values can be
+ * unpacked.  What happens in such cases is partly accidental (alas).
+ */
+
+/* The pack routines write 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
+ * first, at p).
+ * Return value:  0 if all is OK, -1 if error (and an exception is
+ * set, most likely OverflowError).
+ * Bug:  What this does is undefined if x is a NaN or infinity.
+ * Bug:  -0.0 and +0.0 produce the same string.
+ */
+PyAPI_FUNC(int) _PyFloat_Pack4(double x, unsigned char *p, int le);
+PyAPI_FUNC(int) _PyFloat_Pack8(double x, unsigned char *p, int le);
+
+/* The unpack routines read 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).
+ * Return value:  The unpacked double.  On error, this is -1.0 and
+ * PyErr_Occurred() is true (and an exception is set, most likely
+ * OverflowError).
+ * Bug:  What this does is undefined if the string represents a NaN or
+ * infinity.
+ */
+PyAPI_FUNC(double) _PyFloat_Unpack4(const unsigned char *p, int le);
+PyAPI_FUNC(double) _PyFloat_Unpack8(const unsigned char *p, int le);
+
+
 #ifdef __cplusplus
 }
 #endif

Modules/cPickle.c

@@ -1124,106 +1124,10 @@ save_float(Picklerobject *self, PyObject *args)
 	double x = PyFloat_AS_DOUBLE((PyFloatObject *)args);
 
 	if (self->bin) {
-		int s, e;
-		double f;
-		long fhi, flo;
 		char str[9];
-		unsigned char *p = (unsigned char *)str;
+		str[0] = BINFLOAT;
-
+		if (_PyFloat_Pack8(x, (unsigned char *)&str[1], 0) < 0)
-		*p = BINFLOAT;
-		p++;
-
-		if (x < 0) {
-			s = 1;
-			x = -x;
-		}
-		else
-			s = 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 = (long) floor(f); /* Truncate */
-		assert(fhi < 268435456);
-
-		f -= (double)fhi;
-		f *= 16777216.0; /* 2**24 */
-		flo = (long) floor(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 = (s<<7) | (e>>4);
-		p++;
-
-		/* Second byte */
-		*p = (unsigned char) (((e&0xF)<<4) | (fhi>>24));
-		p++;
-
-		/* Third byte */
-		*p = (unsigned char) ((fhi>>16) & 0xFF);
-		p++;
-
-		/* Fourth byte */
-		*p = (unsigned char) ((fhi>>8) & 0xFF);
-		p++;
-
-		/* Fifth byte */
-		*p = (unsigned char) (fhi & 0xFF);
-		p++;
-
-		/* Sixth byte */
-		*p = (unsigned char) ((flo>>16) & 0xFF);
-		p++;
-
-		/* Seventh byte */
-		*p = (unsigned char) ((flo>>8) & 0xFF);
-		p++;
-
-		/* Eighth byte */
-		*p = (unsigned char) (flo & 0xFF);
-
 		if (self->write_func(self, str, 9) < 0)
 			return -1;
 	}
@@ -1237,11 +1141,6 @@ save_float(Picklerobject *self, PyObject *args)
 	}
 
 	return 0;
-
- Overflow:
-	PyErr_SetString(PyExc_OverflowError,
-			"float too large to pack with d format");
-	return -1;
 }
 
 
@@ -3372,64 +3271,20 @@ load_float(Unpicklerobject *self)
 static int
 load_binfloat(Unpicklerobject *self)
 {
-	PyObject *py_float = 0;
+	PyObject *py_float;
-	int s, e;
-	long fhi, flo;
 	double x;
 	char *p;
 
 	if (self->read_func(self, &p, 8) < 0)
 		return -1;
 
-	/* First byte */
+	x = _PyFloat_Unpack8((unsigned char *)p, 0);
-	s = (*p>>7) & 1;
+	if (x == -1.0 && PyErr_Occurred())
-	e = (*p & 0x7F) << 4;
+		return -1;
-	p++;
 
-	/* Second byte */
+	py_float = PyFloat_FromDouble(x);
-	e |= (*p>>4) & 0xF;
+	if (py_float == NULL)
-	fhi = (*p & 0xF) << 24;
+		return -1;
-	p++;
-
-	/* Third byte */
-	fhi |= (*p & 0xFF) << 16;
-	p++;
-
-	/* Fourth byte */
-	fhi |= (*p & 0xFF) << 8;
-	p++;
-
-	/* Fifth byte */
-	fhi |= *p & 0xFF;
-	p++;
-
-	/* Sixth byte */
-	flo = (*p & 0xFF) << 16;
-	p++;
-
-	/* Seventh byte */
-	flo |= (*p & 0xFF) << 8;
-	p++;
-
-	/* Eighth byte */
-	flo |= *p & 0xFF;
-
-	x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
-	x /= 268435456.0; /* 2**28 */
-
-	/* XXX This sadly ignores Inf/NaN */
-	if (e == 0)
-		e = -1022;
-	else {
-		x += 1.0;
-		e -= 1023;
-	}
-	x = ldexp(x, e);
-
-	if (s)
-		x = -x;
-
-	if (!( py_float = PyFloat_FromDouble(x)))  return -1;
 
 	PDATA_PUSH(self->stack, py_float, -1);
 	return 0;

Modules/structmodule.c

@@ -185,301 +185,27 @@ get_ulonglong(PyObject *v, unsigned LONG_LONG *p)
 
 /* Floating point helpers */
 
-/* These use ANSI/IEEE Standard 754-1985 (Standard for Binary Floating
-   Point Arithmetic).  See the following URL:
-   http://www.psc.edu/general/software/packages/ieee/ieee.html */
-
-/* XXX Inf/NaN are not handled quite right (but underflow is!) */
-
-static int
-pack_float(double x, /* The number to pack */
-           char *p,  /* Where to pack the high order byte */
-           int incr) /* 1 for big-endian; -1 for little-endian */
-{
-	int s;
-	int e;
-	double f;
-	long fbits;
-
-	if (x < 0) {
-		s = 1;
-		x = -x;
-	}
-	else
-		s = 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 = (long) floor(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 = (s<<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;
-
- Overflow:
-	PyErr_SetString(PyExc_OverflowError,
-			"float too large to pack with f format");
-	return -1;
-}
-
-static int
-pack_double(double x, /* The number to pack */
-            char *p,  /* Where to pack the high order byte */
-            int incr) /* 1 for big-endian; -1 for little-endian */
-{
-	int s;
-	int e;
-	double f;
-	long fhi, flo;
-
-	if (x < 0) {
-		s = 1;
-		x = -x;
-	}
-	else
-		s = 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 = (long) floor(f); /* Truncate */
-	assert(fhi < 268435456);
-
-	f -= (double)fhi;
-	f *= 16777216.0; /* 2**24 */
-	flo = (long) floor(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 = (s<<7) | (e>>4);
-	p += incr;
-
-	/* Second byte */
-	*p = (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;
-
-	/* Done */
-	return 0;
-
- Overflow:
-	PyErr_SetString(PyExc_OverflowError,
-			"float too large to pack with d format");
-	return -1;
-}
-
 static PyObject *
-unpack_float(const char *p,  /* Where the high order byte is */
+unpack_float(const char *p,  /* start of 4-byte string */
-             int incr)       /* 1 for big-endian; -1 for little-endian */
+             int le)	     /* true for little-endian, false for big-endian */
 {
-	int s;
-	int e;
-	long f;
 	double x;
 
-	/* First byte */
+	x = _PyFloat_Unpack4((unsigned char *)p, le);
-	s = (*p>>7) & 1;
+	if (x == -1.0 && PyErr_Occurred())
-	e = (*p & 0x7F) << 1;
+		return NULL;
-	p += incr;
-
-	/* Second byte */
-	e |= (*p>>7) & 1;
-	f = (*p & 0x7F) << 16;
-	p += incr;
-
-	/* Third byte */
-	f |= (*p & 0xFF) << 8;
-	p += incr;
-
-	/* Fourth byte */
-	f |= *p & 0xFF;
-
-	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 (s)
-		x = -x;
-
 	return PyFloat_FromDouble(x);
 }
 
 static PyObject *
-unpack_double(const char *p,  /* Where the high order byte is */
+unpack_double(const char *p,  /* start of 8-byte string */
-              int incr)       /* 1 for big-endian; -1 for little-endian */
+              int le)         /* true for little-endian, false for big-endian */
 {
-	int s;
-	int e;
-	long fhi, flo;
 	double x;
 
-	/* First byte */
+	x = _PyFloat_Unpack8((unsigned char *)p, le);
-	s = (*p>>7) & 1;
+	if (x == -1.0 && PyErr_Occurred())
-	e = (*p & 0x7F) << 4;
+		return NULL;
-	p += incr;
-
-	/* Second byte */
-	e |= (*p>>4) & 0xF;
-	fhi = (*p & 0xF) << 24;
-	p += incr;
-
-	/* Third byte */
-	fhi |= (*p & 0xFF) << 16;
-	p += incr;
-
-	/* Fourth byte */
-	fhi |= (*p & 0xFF) << 8;
-	p += incr;
-
-	/* Fifth byte */
-	fhi |= *p & 0xFF;
-	p += incr;
-
-	/* Sixth byte */
-	flo = (*p & 0xFF) << 16;
-	p += incr;
-
-	/* Seventh byte */
-	flo |= (*p & 0xFF) << 8;
-	p += incr;
-
-	/* Eighth byte */
-	flo |= *p & 0xFF;
-	p += incr;
-
-	x = (double)fhi + (double)flo / 16777216.0; /* 2**24 */
-	x /= 268435456.0; /* 2**28 */
-
-	/* XXX This sadly ignores Inf/NaN */
-	if (e == 0)
-		e = -1022;
-	else {
-		x += 1.0;
-		e -= 1023;
-	}
-	x = ldexp(x, e);
-
-	if (s)
-		x = -x;
-
 	return PyFloat_FromDouble(x);
 }
 
@@ -887,13 +613,13 @@ bu_ulonglong(const char *p, const formatdef *f)
 static PyObject *
 bu_float(const char *p, const formatdef *f)
 {
-	return unpack_float(p, 1);
+	return unpack_float(p, 0);
 }
 
 static PyObject *
 bu_double(const char *p, const formatdef *f)
 {
-	return unpack_double(p, 1);
+	return unpack_double(p, 0);
 }
 
 static int
@@ -967,7 +693,7 @@ bp_float(char *p, PyObject *v, const formatdef *f)
 				"required argument is not a float");
 		return -1;
 	}
-	return pack_float(x, p, 1);
+	return _PyFloat_Pack4(x, (unsigned char *)p, 0);
 }
 
 static int
@@ -979,7 +705,7 @@ bp_double(char *p, PyObject *v, const formatdef *f)
 				"required argument is not a float");
 		return -1;
 	}
-	return pack_double(x, p, 1);
+	return _PyFloat_Pack8(x, (unsigned char *)p, 0);
 }
 
 static formatdef bigendian_table[] = {
@@ -1053,13 +779,13 @@ lu_ulonglong(const char *p, const formatdef *f)
 static PyObject *
 lu_float(const char *p, const formatdef *f)
 {
-	return unpack_float(p+3, -1);
+	return unpack_float(p, 1);
 }
 
 static PyObject *
 lu_double(const char *p, const formatdef *f)
 {
-	return unpack_double(p+7, -1);
+	return unpack_double(p, 1);
 }
 
 static int
@@ -1133,7 +859,7 @@ lp_float(char *p, PyObject *v, const formatdef *f)
 				"required argument is not a float");
 		return -1;
 	}
-	return pack_float(x, p+3, -1);
+	return _PyFloat_Pack4(x, (unsigned char *)p, 1);
 }
 
 static int
@@ -1145,7 +871,7 @@ lp_double(char *p, PyObject *v, const formatdef *f)
 				"required argument is not a float");
 		return -1;
 	}
-	return pack_double(x, p+7, -1);
+	return _PyFloat_Pack8(x, (unsigned char *)p, 1);
 }
 
 static formatdef lilendian_table[] = {

Objects/floatobject.c

@@ -904,3 +904,316 @@ PyFloat_Fini(void)
 		}
 	}
 }
+
+/*----------------------------------------------------------------------------
+ * _PyFloat_{Pack,Unpack}{4,8}.  See floatobject.h.
+ *
+ * TODO:  On platforms that use the standard IEEE-754 single and double
+ * formats natively, these routines could simply copy the bytes.
+ */
+int
+_PyFloat_Pack4(double x, unsigned char *p, int le)
+{
+	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 = (long) floor(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;
+
+ Overflow:
+	PyErr_SetString(PyExc_OverflowError,
+			"float too large to pack with f format");
+	return -1;
+}
+
+int
+_PyFloat_Pack8(double x, unsigned char *p, int le)
+{
+	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;
+
+	/* Done */
+	return 0;
+
+ Overflow:
+	PyErr_SetString(PyExc_OverflowError,
+			"float too large to pack with d format");
+	return -1;
+}
+
+double
+_PyFloat_Unpack4(const unsigned char *p, int le)
+{
+	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;
+
+	/* 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;
+}
+
+double
+_PyFloat_Unpack8(const unsigned char *p, int le)
+{
+	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;
+
+	/* 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 */
+
+	/* XXX This sadly ignores Inf/NaN */
+	if (e == 0)
+		e = -1022;
+	else {
+		x += 1.0;
+		e -= 1023;
+	}
+	x = ldexp(x, e);
+
+	if (sign)
+		x = -x;
+
+	return x;
+}