# Implementing the Farey sequence efficiently

There is of course the silly implementation:

FareySequence[n_] := Union[Flatten[Table[j/i, {i, 1, n}, {j, 0, i}]]]


However, there are numerous properties and confinements of Farey sequences (that can be used, potentially, in an indirect manner).

This calls for a very simple, and, very efficient recurring/functional implementation, exhibiting Superiority. But I'm new to Mathematica and can't find the right combination of built-in functions, and pure functions..

Ideas?

• Hmm... Commented May 6, 2013 at 4:23
• I saw the implementation in the link. Isn't there a simpler implementation? Commented May 6, 2013 at 4:34
• Your implementation is certainly simpler. Of course, it's also less efficient. Commented May 6, 2013 at 6:03
• I think the problem with finding a recursive solution using the cited property is that the property is not strictly speaking a recurrence relation, since $D_k$ depends on $F_k$. Indeed for any $F_{k-1}$, there are infinitely many solutions for $F_k$ (depending on the order). Commented May 7, 2013 at 2:35
• Yes, perhaps you should clarify whether you want an efficient way to generate the Farey sequences or to be shown how to implement certain properties. Commented May 7, 2013 at 13:04

Graham, Knuth, and Patashnik in their book Concrete Mathematics (pages 118 and 150) discuss the Farey series. Their recurrence does not require finding Subsets, computing the elements in order starting with $0/1$ and $1/n$. Although very fast, Subsets can use too much memory when very large $n$ are required, as for some PE problems.

FareyIterate[{f1_,f2_},n_Integer]:=
{f2,(#*Numerator[f2]-Numerator[f1])/(#*Denominator[f2]-Denominator[f1])}&
[Floor[(Denominator[f1]+n)/Denominator[f2]]]

FareyLength[n_Integer]:=Total[EulerPhi[Range[n]]] + 1

ConcreteFarey[n_]:=NestList[FareyIterate[#,n]&, {0, 1/n}, FareyLength[n]-1][[All,1]]


A NestWhile formulation is possible to pick out certain values without storing the entire list. Nevertheless, this function is only half as fast as farey2[n] of @Michael E2 and @J.M.

• This is effectively the method I used for plotting Thomae's function (which I linked to in the comments), but was glibly dismissed by the OP. On that note, I'd use NestWhile[] if I wanted to avoid the use of EulerPhi[], and would replace Floor[(Denominator[f1]+n)/Denominator[f2]] with Quotient[Denominator[f1] + n, Denominator[f2]]. Commented May 7, 2013 at 19:22
• Thank you @J.M. for these performance tips. I would also be interested to know your experience with GCD versus CoprimeQ in farey2[n]. I've found GCD faster in some cases... Commented May 7, 2013 at 19:42
• If you use Divide @@@ NestList[FareyIterate[#, n] &, {0, 1, 1, n}, FareyLength[n] - 1][[All, 1 ;; 2]] for ConcreteFarey and alter to FareyIterate[{f1N_, f1D_, f2N_, f2D_}, n_Integer] etc -- J.M. beat me to recommending Quotient -- the speed is about halfway between your original and farey2. (Integers are faster than Rationals.) Commented May 7, 2013 at 19:52
• @Michael E2: Brilliant. Your recommendations bring the timings down to 2/3 of previous values, much closer to those of farey2[n]. Commented May 7, 2013 at 20:13

Here's a way to exploit the mediant property of the Farey series. To calculate the mediant:

med[{a_, b_}] := (Numerator[a] + Numerator[b])/(Denominator[a] + Denominator[b]);


Then the Farey series is:

farey[n_] := farey[n] = DeleteCases[ Riffle[
farey[n - 1], med /@ Partition[farey[n - 1], 2, 1]], _?(Denominator[#] > n &)];


with initial conditions

farey[2]={0,1/2,1}


Now you can get farey[n] for any fixed n straightforwardly.

• med[v_List] := Total[Numerator[v]]/Total[Denominator[v]] is a bit more compact. Commented May 6, 2013 at 8:29
• That's good, but you know, this whole approach doesn't seem any faster than the simple way! Commented May 6, 2013 at 8:33
• Well, the OP did dismiss my proposal so glibly... Commented May 6, 2013 at 8:46
• @Dror Well, you did ask for a functional approach and this certainly provides it. Also, you never indicated that speed was your primary concern and I see no immediate reason that it should be. The primary advantage of this approach that I see is the clarity provided by its immediate connection to the mediant. By removing the DeleteCases step, for example, we essentially recover the Stern-Brocot tree. Commented May 6, 2013 at 13:32
• @dror - It is interesting that the functional answer is so much slower than the "simple" answer. Maybe it's because the very basic commands like Union and Table are quicker than comparatively obscure commands like Numerator and DeleteCases. Commented May 6, 2013 at 16:00

Here's a functional way to use the property (the property, which has been removed from the original question, was $N'/D' = N/D + 1/D'D$ or equivalently $N'D-D'N=1$):

farey1[n_] :=
NestWhileList[
With[{num0 = Numerator[#], den0 = Denominator[#]},
First @ Minimize[{num/den,
num den0 - num0 den == 1 && 1 <= den <= n && 1 <= num <= n},
{num, den}, Integers]] &,
1/n,
# < 1 - 1/n &]


Mighty slow:

foo1 = farey1[15]; // Timing
(* {0.441386, Null} *)


Here's faster way, without using the property:

farey2[n_] := Sort @ Pick[Rational @@@ #, GCD @@ Transpose@#, 1] &@ Subsets[Range[n], {2}];

foo2 = Farey2[15]; // Timing
(* {0.000193, Null} *)


Following J.M.'s comment, this is more succinct:

farey2[n_] := Sort @ Pick[Divide @@@ #, CoprimeQ @@@ #] &@ Subsets[Range[n], {2}];


Both methods give the "interior" of the traditional sequence, omitting 0 and 1:

farey1[5]
(* {1/5, 1/4, 1/3, 2/5, 1/2, 3/5, 2/3, 3/4, 4/5} *)

• Any reason why you use Rational @@@ (* stuff *) instead of Divide @@@ (* stuff *)? Anyway, something for your consideration: Sort[Pick[Divide @@@ #, CoprimeQ @@@ #]] &[Subsets[Range[n], {2}]]. Commented May 7, 2013 at 14:38
• @J.M. No particular reason for Rational -- any reason not to? Thanks for the CoprimeQ tip -- it's a bit faster by a little. Commented May 7, 2013 at 14:50
• "any reason not to?" - not really; I was asking out of curiosity more than anything, since I conventionally use Divide in such situations. Commented May 7, 2013 at 16:19
• @J.M. In my answer I adapted some pre-Pick code from way back, which had Rational. At the time, I may have thought it was more direct since Divide[1, 2] evaluates to Rational[1, 2]. It seems to be about 5% faster, but that hardly matters here. Commented May 7, 2013 at 18:29

And with an updated Mathematica, there is an updated (and far simpler) answer: FareySequence is now built in. For example:

FareySequence[4]
{0, 1/4, 1/3, 1/2, 2/3, 3/4, 1}


There is a very nice, and relatively new (2008) algorithm for the Farey sequence that is extremely efficient. It computes any Farey sequence in just one pass, in order, with very little compute time or overhead.

It can be found on Wikipedia (Farey Sequence) and is cited there as:

Routledge, Norman (March 2008). "Computing Farey series". The Mathematical Gazette. Vol. 92, no. 523. pp. 55–62.

In javascript:

function farey(n) {
let a = 0,
b = 1,
c = 1,
d = n,
result = "0/1";

while (c <= n) {
let k = Math.floor((n + b) / d);
let a_new = a,
b_new = b,
c_new = c,
d_new = d;

a = c_new;
b = d_new;
c = k * c_new - a_new;
d = k * d_new - b_new;
result += , ${a}/${b};
}

return result;
}


farey(6) returns:

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

This is extraordinary because it doesn't need recursion, nor nested loops, nor insertion within a list (only append the next value, in order). It also does not need GCD, as it immediately finds the reduced form of each term with the clever use of k in the algorithm.

I think this might satisfy the OP's request for an algorithm exhibiting "superiority."

A good explanation of how it works is in the Wikipedia article.

• Providing the Mathematic code for this algorithm would be a service to the community. Commented Jan 9, 2023 at 22:57
• @bgodfrey I think the response by Kenny Colnago effectively provides said code. Commented Jan 10, 2023 at 15:25