ActiveState Code

Recipe 496821: Farey Sequence

Function provides farey sequence, F(n), for any integer n.

Python 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
def farey(n):
    def gcd(a,b):
        while b: a,b = b,a%b
        return a

    def simplify(a,b):
        g = gcd(a,b)
        return (a/g,b/g)

    fs = dict()
    for i in xrange(1,n+1):
        for i2 in xrange(1,i+1):
            if i2 < n and i != i2:
                r = simplify(i2,i)
                fs[float(i2)/i] = r

    return [fs[k] for k in sorted(fs.keys())]

# Usage:
# print farey(5)
# => [(1, 5), (1, 4), (1, 3), (2, 5), (1, 2), (3, 5), (2, 3), (3, 4), (4, 5)]

Discussion

There's probably a more elegant solution out there, but I couldn't find it. Note: it doesn't prepend (0,1) and append (1,1)...

For more info: http://mathworld.wolfram.com/FareySequence.html

Comments

  1. 1. At 9:49 a.m. on 26 jun 2006, Scott David Daniels said:

    Some simpler ways based on the Farey properties. See also my Recipe 52317 (whose comments section seems mangled now).

    In the following we have a straight-forward recursive generator named "Farey_r" and an equivalent, but more efficient, generator named "farey". "Farey_r" is provided simply to make it 'clear' (and more provable) what the generator is computing.

    def Farey_r(limit, start=(0, 1), stop=(1, 1)):
    '''recursive definition of a Farey sequence generator'''
    n, d = start
    N, D = stop
    mediant_d = d + D
    if mediant_d &lt;= limit:
        mediant = (n + N), mediant_d
        for pair in Farey_r(limit, start, mediant): yield pair
        for pair in Farey_r(limit, mediant, stop): yield pair
    else:
        yield start


def farey(limit):
    '''Fast computation of Farey sequence as a generator'''
    # n, d is the start fraction n/d (0,1) initially
    # N, D is the stop fraction N/D (1,1) initially
    pend = []
    n = 0
    d = N = D = 1
    while True:
        mediant_d = d + D
        if mediant_d &lt;= limit:
            mediant_n = n + N
            pend.append((mediant_n, mediant_d, N, D))
            N = mediant_n
            D = mediant_d
        else:
            yield n, d
            if pend:
                n, d, N, D = pend.pop()
            else:
                break
  2. 2. At 5:56 p.m. on 27 jun 2006, James Kassemi (the author) said:

    Yep. Sorry about not referencing your recipe, I thought I did a search, but I must have missed it.

    For those who didn't see it, take a look, the mediant calculation approach is certainly more effective.

Sign in to comment