17.17 Converting Numbers to Rationals via Farey Fractions

Credit: Scott David Daniels

17.17.1 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.

17.17.2 Solution

Farey fractions, whose crucial properties were studied by 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 right type" for denominator
    lower, upper = (z, z+1), (z+1, z)
    while 1:
        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, farey(math.pi, 100) == (22, 7).

17.17.3 Discussion

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

Note the trickiness with z. It is a zero of the same type as the lim argument. This lets you use longs as the limit if necessary, without paying a performance price (not even a test) when there's no such need.

To print odds, you can use:

n, d = farey(probability, lim)
print "Odds are %d : %d" % (n, d-n)

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

If you are using this in an environment where you call it with a lot of values near 0.0, 1.0, or 0.5 (or simple fractions), you may find that its convergence is too slow. You can improve its 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 surveyor who wrote a letter to the Journal of Science around the end of the 18th century. In that letter he observed that, while reading a privately published list of the decimal equivalents of fractions, he noticed the following: for any three consecutive fractions in the simplest terms (e.g., A/B, C/D, E/F), the middle one (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.

17.17.4 See Also

Recipe 17.18 for another mathematical evaluation recipe.