Welcome, guest | Sign In | My Account | Store | Cart

farey: Numeric to rational via Farey fractions (Python recipe) by Scott David Daniels
ActiveState Code (http://code.activestate.com/recipes/52317/)

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, 22 lines
Download
 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
Tags: algorithms

5 comments

Scott David Daniels (author) 23 years, 1 month ago  # | flag

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.

Connelly Barnes 18 years, 6 months ago  # | flag

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.

Dick Moores 17 years, 4 months ago  # | flag

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

Peter Nau 15 years, 7 months ago  # | flag

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.

jerfelix 12 years, 3 months ago  # | flag

I tend to agree with Dick Moores, above, as opposed to Peter Nau. The second parameter specifies the limit to the denominator, and you'd think that the integer ratio which best matches the first parameter and keeps the denominator below the second parameter would be the match.

Here's some code to address the accuracy issue. My code uses the 2nd parameter (lim) as a value that the denominator must be strictly less than (as opposed to less than OR equal to).

from __future__ import division

def not_farey(v,lim):
    d= min((abs(int(round(v*d))/d-v),d) for d in xrange(1,lim))[1]
    return int(round(v*d)),d