Implementing the Farey sequence efficiently

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.


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.


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} *)