Python Cookbook 2Nd Edition Jun 1002005 [Electronic resources] نسخه متنی

Recipe 18.13. Converting Numbers to Rationals via Farey Fractions

Credit: Scott David Daniels

Problem

You have a number
v (of almost any type) and need to find a
rational number (in reduced form) that is as close to
v as possible but with a denominator no
larger than a prescribed value.

Solution

Farey fractions, whose crucial properties were studied by Augustin
Louis Cauchy, are an excellent way to find rational approximations of
floating-point values:

def farey(v, lim):
"" No error checking on args.  lim = maximum denominator.
Results are (numerator, denominator); (1, 0) is "infinity".
""
if v < 0:
n, d = farey(-v, lim)
return -n, d
z = lim - lim    # Get a "0 of the right type" for denominator
lower, upper = (z, z+1), (z+1, z)
while True:
mediant = (lower[0] + upper[0]), (lower[1] + upper[1])
if v * mediant[1] > mediant[0]:
if lim < mediant[1]:
return upper
lower = mediant
elif v * mediant[1] == mediant[0]:
if lim >= mediant[1]:
return mediant
if lower[1] < upper[1]:
return lower
return upper
else:
if lim < mediant[1]:
return lower
upper = mediant

For example:

import math
print farey(math.pi, 100)
# emits: (22, 7)

Discussion

The rationals resulting from the algorithm shown in this recipe are
in reduced form (meaning that numerator and denominator are mutually
prime), but the proof, which was given by Cauchy, is rather subtle
(see http://www.cut-the-knot.com/blue/Fareyl).

You can use farey to compute odds from a
probability, such as:

probability = 0.333
n, d = farey(probability, 100)
print "Odds are %d : %d" % (n, d-n)
# emits: Odds are 1 : 2

This recipe's algorithm is ideally suited for
reimplementation in a lower-level language (e.g., C, or assembly, or,
maybe best, Pyrex) if you use it heavily. Since the code uses only
multiplication and addition, it can play optimally to hardware
strengths.

If you are using this recipe in an environment where you call it with
a lot of values near 0.0, 1.0, or 0.5 (or other simple fractions),
you may find that the algorithm's convergence is too
slow. You can improve convergence in a continued fraction style, by
appending to the first if in the
farey function:

if v < 0:
...
elif v < 0.5:
n, d = farey((v-v+1)/v, lim)     # lim is wrong; decide what you want
return d, n
elif v > 1:
intpart = floor(v)
n, d = farey(v-intpart)
return n+intpart*d, d
...

James Farey was an English geologist and surveyor who wrote a letter
to the Journal of Science in 1816. In that
letter he observed that, while reading a privately published list of
the decimal equivalents of fractions, he had noticed an interesting
fact. Consider the set of all the fractions with values between 0 and
1, reduced to the lowest terms, with denominators not exceeding some
integer N. Arrange the set in order of
magnitude to get a sequence. For example, for
N equal to 5, the Farey sequence is:

0/1, 1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5, 1/1

For any three consecutive fractions in this sequence (e.g., A/B, C/D,
E/F), the middle fraction (C/D), called the
mediant, is equal to the ratio (A + E)/(B + F).
I enjoy envisioning Mr. Farey sitting up late on a rainy English
night, reading tables of decimal expansions of fractions by an oil
lamp. Calculation has come a long way since his day, and
I'm pleased to be able to benefit from his work.

See Also

Library Reference and Python in a
Nutshell docs for built-in types int
and long; http://www.cut-the-knot.org/blue/Farey.shtml
for more information about the Farey Series.