Implementation of Fraction.limit_denominator.

Remove Fraction.to_continued_fraction and
Fraction.from_continued_fraction

Doc/library/fractions.rst

@@ -46,6 +46,24 @@ Fraction number class.
    :class:`decimal.Decimal`.
 
 
+.. method:: Fraction.limit_denominator(max_denominator=1000000)
+
+   Finds and returns the closest :class:`Fraction` to ``self`` that
+   has denominator at most max_denominator.  This method is useful for
+   finding rational approximations to a given floating-point number::
+
+      >>> Fraction('3.1415926535897932').limit_denominator(1000)
+      Fraction(355, 113)
+
+   or for recovering a rational number that's represented as a float::
+
+      >>> from math import pi, cos
+      >>> Fraction.from_float(cos(pi/3))
+      Fraction(4503599627370497L, 9007199254740992L)
+      >>> Fraction.from_float(cos(pi/3)).limit_denominator()
+      Fraction(1, 2)
+
+
 .. method:: Fraction.__floor__()
 
    Returns the greatest :class:`int` ``<= self``. Will be accessible

Lib/fractions.py

@@ -140,42 +140,60 @@ class Fraction(Rational):
         else:
             return Fraction(digits, 10 ** -exp)
 
-    @staticmethod
-    def from_continued_fraction(seq):
-        'Build a Fraction from a continued fraction expessed as a sequence'
-        n, d = 1, 0
-        for e in reversed(seq):
-            n, d = d, n
-            n += e * d
-        return Fraction(n, d) if seq else Fraction(0)
-
-    def as_continued_fraction(self):
-        'Return continued fraction expressed as a list'
-        n = self.numerator
-        d = self.denominator
-        cf = []
-        while d:
-            e = int(n // d)
-            cf.append(e)
-            n -= e * d
-            n, d = d, n
-        return cf
-
-    def approximate(self, max_denominator):
-        'Best rational approximation with a denominator <= max_denominator'
-        # XXX First cut at algorithm
-        # Still needs rounding rules as specified at
-        #       http://en.wikipedia.org/wiki/Continued_fraction
+    def limit_denominator(self, max_denominator=1000000):
+        """Closest Fraction to self with denominator at most max_denominator.
+
+        >>> Fraction('3.141592653589793').limit_denominator(10)
+        Fraction(22, 7)
+        >>> Fraction('3.141592653589793').limit_denominator(100)
+        Fraction(311, 99)
+        >>> Fraction(1234, 5678).limit_denominator(10000)
+        Fraction(1234, 5678)
+
+        """
+        # Algorithm notes: For any real number x, define a *best upper
+        # approximation* to x to be a rational number p/q such that:
+        #
+        #   (1) p/q >= x, and
+        #   (2) if p/q > r/s >= x then s > q, for any rational r/s.
+        #
+        # Define *best lower approximation* similarly.  Then it can be
+        # proved that a rational number is a best upper or lower
+        # approximation to x if, and only if, it is a convergent or
+        # semiconvergent of the (unique shortest) continued fraction
+        # associated to x.
+        #
+        # To find a best rational approximation with denominator <= M,
+        # we find the best upper and lower approximations with
+        # denominator <= M and take whichever of these is closer to x.
+        # In the event of a tie, the bound with smaller denominator is
+        # chosen.  If both denominators are equal (which can happen
+        # only when max_denominator == 1 and self is midway between
+        # two integers) the lower bound---i.e., the floor of self, is
+        # taken.
+
+        if max_denominator < 1:
+            raise ValueError("max_denominator should be at least 1")
         if self.denominator <= max_denominator:
-            return self
-        cf = self.as_continued_fraction()
-        result = Fraction(0)
-        for i in range(1, len(cf)):
-            new = self.from_continued_fraction(cf[:i])
-            if new.denominator > max_denominator:
+            return Fraction(self)
+
+        p0, q0, p1, q1 = 0, 1, 1, 0
+        n, d = self.numerator, self.denominator
+        while True:
+            a = n//d
+            q2 = q0+a*q1
+            if q2 > max_denominator:
                 break
-            result = new
-        return result
+            p0, q0, p1, q1 = p1, q1, p0+a*p1, q2
+            n, d = d, n-a*d
+
+        k = (max_denominator-q0)//q1
+        bound1 = Fraction(p0+k*p1, q0+k*q1)
+        bound2 = Fraction(p1, q1)
+        if abs(bound2 - self) <= abs(bound1-self):
+            return bound2
+        else:
+            return bound1
 
     @property
     def numerator(a):

Lib/test/test_fractions.py

@@ -188,28 +188,15 @@ class FractionTest(unittest.TestCase):
             TypeError, "Cannot convert sNaN to Fraction.",
             R.from_decimal, Decimal("snan"))
 
-    def testFromContinuedFraction(self):
-        self.assertRaises(TypeError, R.from_continued_fraction, None)
-        phi = R.from_continued_fraction([1]*100)
-        self.assertEquals(round(phi - (1 + 5 ** 0.5) / 2, 10), 0.0)
-
-        minusphi = R.from_continued_fraction([-1]*100)
-        self.assertEquals(round(minusphi + (1 + 5 ** 0.5) / 2, 10), 0.0)
-
-        self.assertEquals(R.from_continued_fraction([0]), R(0))
-        self.assertEquals(R.from_continued_fraction([]), R(0))
-
-    def testAsContinuedFraction(self):
-        self.assertEqual(R.from_float(math.pi).as_continued_fraction()[:15],
-                         [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3, 3])
-        self.assertEqual(R.from_float(-math.pi).as_continued_fraction()[:16],
-                         [-4, 1, 6, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3, 3])
-        self.assertEqual(R(0).as_continued_fraction(), [0])
-
-    def testApproximateFrom(self):
-        self.assertEqual(R.from_float(math.pi).approximate(10000), R(355, 113))
-        self.assertEqual(R.from_float(-math.pi).approximate(10000), R(-355, 113))
-        self.assertEqual(R.from_float(0.0).approximate(10000), R(0))
+    def testLimitDenominator(self):
+        rpi = R('3.1415926535897932')
+        self.assertEqual(rpi.limit_denominator(10000), R(355, 113))
+        self.assertEqual(-rpi.limit_denominator(10000), R(-355, 113))
+        self.assertEqual(rpi.limit_denominator(113), R(355, 113))
+        self.assertEqual(rpi.limit_denominator(112), R(333, 106))
+        self.assertEqual(R(201, 200).limit_denominator(100), R(1))
+        self.assertEqual(R(201, 200).limit_denominator(101), R(102, 101))
+        self.assertEqual(R(0).limit_denominator(10000), R(0))
 
     def testConversions(self):
         self.assertTypedEquals(-1, math.trunc(R(-11, 10)))