diff options
| -rw-r--r-- | Doc/library/fractions.rst | 18 | ||||
| -rwxr-xr-x | Lib/fractions.py | 86 | ||||
| -rw-r--r-- | Lib/test/test_fractions.py | 31 |
3 files changed, 79 insertions, 56 deletions
diff --git a/Doc/library/fractions.rst b/Doc/library/fractions.rst index 7f6019f..adc70c6 100644 --- a/Doc/library/fractions.rst +++ b/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 diff --git a/Lib/fractions.py b/Lib/fractions.py index 123ecb6..8593e7e 100755 --- a/Lib/fractions.py +++ b/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): diff --git a/Lib/test/test_fractions.py b/Lib/test/test_fractions.py index a79fedd..14dd868 100644 --- a/Lib/test/test_fractions.py +++ b/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))) |