diff --git a/Objects/floatobject.c b/Objects/floatobject.c
index 120b561..ba37309 100644
--- a/Objects/floatobject.c
+++ b/Objects/floatobject.c
@@ -359,7 +359,7 @@
 	PyFloatObject *w;
 {
 	double vx, wx;
-	double /* div, */ mod;
+	double mod;
 	wx = w->ob_fval;
 	if (wx == 0.0) {
 		PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
@@ -368,10 +368,10 @@
 	PyFPE_START_PROTECT("modulo", return 0)
 	vx = v->ob_fval;
 	mod = fmod(vx, wx);
-	/* div = (vx - mod) / wx; */
-	if (wx*mod < 0) {
+	/* note: checking mod*wx < 0 is incorrect -- underflows to
+	   0 if wx < sqrt(smallest nonzero double) */
+	if (mod && ((wx < 0) != (mod < 0))) {
 		mod += wx;
-		/* div -= 1.0; */
 	}
 	PyFPE_END_PROTECT(mod)
 	return PyFloat_FromDouble(mod);
@@ -383,7 +383,7 @@
 	PyFloatObject *w;
 {
 	double vx, wx;
-	double div, mod;
+	double div, mod, floordiv;
 	wx = w->ob_fval;
 	if (wx == 0.0) {
 		PyErr_SetString(PyExc_ZeroDivisionError, "float divmod()");
@@ -392,13 +392,25 @@
 	PyFPE_START_PROTECT("divmod", return 0)
 	vx = v->ob_fval;
 	mod = fmod(vx, wx);
+	/* fmod is typically exact, so vx-mod is *mathemtically* an
+	   exact multiple of wx.  But this is fp arithmetic, and fp
+	   vx - mod is an approximation; the result is that div may
+	   not be an exact integral value after the division, although
+	   it will always be very close to one.
+	*/
 	div = (vx - mod) / wx;
-	if (wx*mod < 0) {
+	/* note: checking mod*wx < 0 is incorrect -- underflows to
+	   0 if wx < sqrt(smallest nonzero double) */
+	if (mod && ((wx < 0) != (mod < 0))) {
 		mod += wx;
 		div -= 1.0;
 	}
+	/* snap quotient to nearest integral value */
+	floordiv = floor(div);
+	if (div - floordiv > 0.5)
+		floordiv += 1.0;
 	PyFPE_END_PROTECT(div)
-	return Py_BuildValue("(dd)", div, mod);
+	return Py_BuildValue("(dd)", floordiv, mod);
 }
 
 static double powu(x, n)