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