Moved Rational._binary_float_to_ratio() to float.as_integer_ratio() because

it's useful outside of rational numbers.

This is my first C code that had to do anything significant. Please be more
careful when looking over it.

Lib/rational.py

@@ -25,60 +25,6 @@ def gcd(a, b):
     return a
 
 
-def _binary_float_to_ratio(x):
-    """x -> (top, bot), a pair of ints s.t. x = top/bot.
-
-    The conversion is done exactly, without rounding.
-    bot > 0 guaranteed.
-    Some form of binary fp is assumed.
-    Pass NaNs or infinities at your own risk.
-
-    >>> _binary_float_to_ratio(10.0)
-    (10, 1)
-    >>> _binary_float_to_ratio(0.0)
-    (0, 1)
-    >>> _binary_float_to_ratio(-.25)
-    (-1, 4)
-    """
-    # XXX Move this to floatobject.c with a name like
-    # float.as_integer_ratio()
-
-    if x == 0:
-        return 0, 1
-    f, e = math.frexp(x)
-    signbit = 1
-    if f < 0:
-        f = -f
-        signbit = -1
-    assert 0.5 <= f < 1.0
-    # x = signbit * f * 2**e exactly
-
-    # Suck up CHUNK bits at a time; 28 is enough so that we suck
-    # up all bits in 2 iterations for all known binary double-
-    # precision formats, and small enough to fit in an int.
-    CHUNK = 28
-    top = 0
-    # invariant: x = signbit * (top + f) * 2**e exactly
-    while f:
-        f = math.ldexp(f, CHUNK)
-        digit = trunc(f)
-        assert digit >> CHUNK == 0
-        top = (top << CHUNK) | digit
-        f = f - digit
-        assert 0.0 <= f < 1.0
-        e = e - CHUNK
-    assert top
-
-    # Add in the sign bit.
-    top = signbit * top
-
-    # now x = top * 2**e exactly; fold in 2**e
-    if e>0:
-        return (top * 2**e, 1)
-    else:
-        return (top, 2 ** -e)
-
-
 _RATIONAL_FORMAT = re.compile(
     r'^\s*(?P<sign>[-+]?)(?P<num>\d+)'
     r'(?:/(?P<denom>\d+)|\.(?P<decimal>\d+))?\s*$')
@@ -163,7 +109,7 @@ def from_float(cls, f):
                             (cls.__name__, f, type(f).__name__))
         if math.isnan(f) or math.isinf(f):
             raise TypeError("Cannot convert %r to %s." % (f, cls.__name__))
-        return cls(*_binary_float_to_ratio(f))
+        return cls(*f.as_integer_ratio())
 
     @classmethod
     def from_decimal(cls, dec):

Lib/test/test_builtin.py

@@ -5,7 +5,7 @@
                               run_unittest, run_with_locale
 from operator import neg
 
-import sys, warnings, cStringIO, random, UserDict
+import sys, warnings, cStringIO, random, rational, UserDict
 warnings.filterwarnings("ignore", "hex../oct.. of negative int",
                         FutureWarning, __name__)
 warnings.filterwarnings("ignore", "integer argument expected",
@@ -688,6 +688,25 @@ def __float__(self):
         self.assertAlmostEqual(float(Foo3(21)), 42.)
         self.assertRaises(TypeError, float, Foo4(42))
 
+    def test_floatasratio(self):
+        R = rational.Rational
+        self.assertEqual(R(0, 1),
+                         R(*float(0.0).as_integer_ratio()))
+        self.assertEqual(R(5, 2),
+                         R(*float(2.5).as_integer_ratio()))
+        self.assertEqual(R(1, 2),
+                         R(*float(0.5).as_integer_ratio()))
+        self.assertEqual(R(4728779608739021, 2251799813685248),
+                         R(*float(2.1).as_integer_ratio()))
+        self.assertEqual(R(-4728779608739021, 2251799813685248),
+                         R(*float(-2.1).as_integer_ratio()))
+        self.assertEqual(R(-2100, 1),
+                         R(*float(-2100.0).as_integer_ratio()))
+
+        self.assertRaises(OverflowError, float('inf').as_integer_ratio)
+        self.assertRaises(OverflowError, float('-inf').as_integer_ratio)
+        self.assertRaises(ValueError, float('nan').as_integer_ratio)
+
     def test_getattr(self):
         import sys
         self.assert_(getattr(sys, 'stdout') is sys.stdout)

Objects/floatobject.c

@@ -1161,6 +1161,163 @@ float_float(PyObject *v)
 	return v;
 }
 
+static PyObject *
+float_as_integer_ratio(PyObject *v)
+{
+	double self;
+	double float_part;
+	int exponent;
+	int is_negative;
+	const int chunk_size = 28;
+	PyObject *prev;
+	PyObject *py_chunk = NULL;
+	PyObject *py_exponent = NULL;
+	PyObject *numerator = NULL;
+	PyObject *denominator = NULL;
+	PyObject *result_pair = NULL;
+	PyNumberMethods *long_methods;
+
+#define INPLACE_UPDATE(obj, call) \
+	prev = obj; \
+	obj = call; \
+	Py_DECREF(prev); \
+
+	CONVERT_TO_DOUBLE(v, self);
+
+	if (Py_IS_INFINITY(self)) {
+	  PyErr_SetString(PyExc_OverflowError,
+			  "Cannot pass infinity to float.as_integer_ratio.");
+	  return NULL;
+	}
+#ifdef Py_NAN
+	if (Py_IS_NAN(self)) {
+	  PyErr_SetString(PyExc_ValueError,
+			  "Cannot pass nan to float.as_integer_ratio.");
+	  return NULL;
+	}
+#endif
+
+	if (self == 0) {
+		numerator = PyInt_FromLong(0);
+		if (numerator == NULL) goto error;
+		denominator = PyInt_FromLong(1);
+		if (denominator == NULL) goto error;
+		result_pair = PyTuple_Pack(2, numerator, denominator);
+		/* Hand ownership over to the tuple. If the tuple
+		   wasn't created successfully, we want to delete the
+		   ints anyway. */
+		Py_DECREF(numerator);
+		Py_DECREF(denominator);
+		return result_pair;
+	}
+
+	/* XXX: Could perhaps handle FLT_RADIX!=2 by using ilogb and
+	   scalbn, but those may not be in C89. */
+	PyFPE_START_PROTECT("as_integer_ratio", goto error);
+	float_part = frexp(self, &exponent);
+	is_negative = 0;
+	if (float_part < 0) {
+		float_part = -float_part;
+		is_negative = 1;
+		/* 0.5 <= float_part < 1.0 */
+	}
+	PyFPE_END_PROTECT(float_part);
+	/* abs(self) == float_part * 2**exponent exactly */
+
+	/* Suck up chunk_size bits at a time; 28 is enough so that we
+	   suck up all bits in 2 iterations for all known binary
+	   double-precision formats, and small enough to fit in a
+	   long. */
+	numerator = PyLong_FromLong(0);
+	if (numerator == NULL) goto error;
+
+	long_methods = PyLong_Type.tp_as_number;
+
+	py_chunk = PyLong_FromLong(chunk_size);
+	if (py_chunk == NULL) goto error;
+
+	while (float_part != 0) {
+		/* invariant: abs(self) ==
+		   (numerator + float_part) * 2**exponent exactly */
+		long digit;
+		PyObject *py_digit;
+
+		PyFPE_START_PROTECT("as_integer_ratio", goto error);
+		/* Pull chunk_size bits out of float_part, into digits. */
+		float_part = ldexp(float_part, chunk_size);
+		digit = (long)float_part;
+		float_part -= digit;
+                /* 0 <= float_part < 1 */
+		exponent -= chunk_size;
+		PyFPE_END_PROTECT(float_part);
+
+		/* Shift digits into numerator. */
+		// numerator <<= chunk_size
+		INPLACE_UPDATE(numerator,
+			       long_methods->nb_lshift(numerator, py_chunk));
+		if (numerator == NULL) goto error;
+
+		// numerator |= digit
+		py_digit = PyLong_FromLong(digit);
+		if (py_digit == NULL) goto error;
+		INPLACE_UPDATE(numerator,
+			       long_methods->nb_or(numerator, py_digit));
+		Py_DECREF(py_digit);
+		if (numerator == NULL) goto error;
+	}
+
+	/* Add in the sign bit. */
+	if (is_negative) {
+		INPLACE_UPDATE(numerator,
+			       long_methods->nb_negative(numerator));
+		if (numerator == NULL) goto error;
+	}
+
+	/* now self = numerator * 2**exponent exactly; fold in 2**exponent */
+	denominator = PyLong_FromLong(1);
+	py_exponent = PyLong_FromLong(labs(exponent));
+	if (py_exponent == NULL) goto error;
+	INPLACE_UPDATE(py_exponent,
+		       long_methods->nb_lshift(denominator, py_exponent));
+	if (py_exponent == NULL) goto error;
+	if (exponent > 0) {
+		INPLACE_UPDATE(numerator,
+			       long_methods->nb_multiply(numerator,
+							 py_exponent));
+		if (numerator == NULL) goto error;
+	}
+	else {
+		Py_DECREF(denominator);
+		denominator = py_exponent;
+		py_exponent = NULL;
+	}
+
+	result_pair = PyTuple_Pack(2, numerator, denominator);
+
+#undef INPLACE_UPDATE
+error:
+	Py_XDECREF(py_exponent);
+	Py_XDECREF(py_chunk);
+	Py_XDECREF(denominator);
+	Py_XDECREF(numerator);
+	return result_pair;
+}
+
+PyDoc_STRVAR(float_as_integer_ratio_doc,
+"float.as_integer_ratio() -> (int, int)\n"
+"\n"
+"Returns a pair of integers, not necessarily in lowest terms, whose\n"
+"ratio is exactly equal to the original float. This method raises an\n"
+"OverflowError on infinities and a ValueError on nans. The resulting\n"
+"denominator will be positive.\n"
+"\n"
+">>> (10.0).as_integer_ratio()\n"
+"(167772160L, 16777216L)\n"
+">>> (0.0).as_integer_ratio()\n"
+"(0, 1)\n"
+">>> (-.25).as_integer_ratio()\n"
+"(-134217728L, 536870912L)");
+
 
 static PyObject *
 float_subtype_new(PyTypeObject *type, PyObject *args, PyObject *kwds);
@@ -1349,6 +1506,8 @@ static PyMethodDef float_methods[] = {
 	 "Returns self, the complex conjugate of any float."},
 	{"__trunc__",	(PyCFunction)float_trunc, METH_NOARGS,
          "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},
 	{"__getnewargs__",	(PyCFunction)float_getnewargs,	METH_NOARGS},
 	{"__getformat__",	(PyCFunction)float_getformat,	
 	 METH_O|METH_CLASS,		float_getformat_doc},