ActiveState Code

Recipe 52317: farey: Numeric to rational via Farey fractions

This converts a Numeric to a rational. The result is always in reduced form, but the proof, while possible, is subtle.

farey(math.pi,100) = (22,7)

Python 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
def farey(v, lim):
  '''Named after James Farey, an English surveyor.
  No error checking on args -- lim = max denominator,
  results are (numerator, denominator), (1,0) is infinity
  '''
  if v < 0:
    n,d = farey(-v, lim)
    return -n,d
  z = lim-lim	# get 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

Comments

  1. 1. At 11:40 a.m. on 6 apr 2001, Scott David Daniels (the author) said:

    notes on farey. Note the trickiness with "z" -- it is a zero of the same type as the argument lim. This allows you to use longs as the limit if need be. To print "odds":

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

    This code is ideally suited for re-implementation in a lower level language (say C or assembly) if you have the need or desire for rational or "odds" output. Because this 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 very near 0., 1., or 0.5 (or _very_ simple fractions), you may find it too slow. You may improve its performance in a "continued fraction" style by appending to the first if:

    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 eighteenth 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 least terms (say, A/B, C/D, E/F), the middle one (C/D) is equal to the (A+E)/(B+F).

    I enjoy my own image of 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.

  2. 2. At 4:40 a.m. on 15 oct 2005, Connelly Barnes said:

    Continued fractions. Alternatively, you can use continued fractions: http://mathworld.wolfram.com/ContinuedFraction.html . Note that in the continued fraction, you can use subtractions as well as additions. Thus math.pi == 3+1/(7+1/(16-.003405593314)) ~= 3+1/(7+1/16.) == 355./113.

  3. 3. At 9 a.m. on 4 jan 2007, Dick Moores said:

    farey() often gives incorrect result. For example, farey(0.36, 10) returns (1, 3), whereas (3, 8) is correct.

    For another, farey(0.584115140346, 100) returns (7, 12) whereas (52, 89) is correct.

    Dick Moores rdmoores@gmail.com

  4. 4. At 11:35 p.m. on 1 oct 2008, Peter Nau said:

    Nice work. Thank you: I plan to reuse this code to convert probability floats to fractional representation (similar to the "odds" comment and example).

    I have comments on the comment about farey(...) accuracy:

    The code appears to make an approximation that is increasingly accurate as one enlarges lim.

    Note that farey(0.36, 100) returns (9, 25), which is the exact answer. (3,8) is "less" correct.

    It's instructive to consider the following examples:

    farey(0.584115140346, 100) returns (7, 12)

    farey(0.584115140346, 1000) returns (125, 214)

    farey(0.584115140346, 10000) returns (5479, 9380)

    farey(0.584115140346, 100000) returns (37419, 64061)

    In fact, farey(0.5842696, 100) returns (52, 89), which is "correct" -- to a desired degree of approximation.

    I hope that's helpful.

Sign in to comment