# 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... – J. M. is away May 6 '13 at 4:23
• I saw the implementation in the link. Isn't there a simpler implementation? – user76568 May 6 '13 at 4:34
• Your implementation is certainly simpler. Of course, it's also less efficient. – J. M. is away May 6 '13 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). – Michael E2 May 7 '13 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. – Michael E2 May 7 '13 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]]. – J. M. is away May 7 '13 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... – KennyColnago May 7 '13 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.) – Michael E2 May 7 '13 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]. – KennyColnago May 7 '13 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. – J. M. is away May 6 '13 at 8:29
• That's good, but you know, this whole approach doesn't seem any faster than the simple way! – bill s May 6 '13 at 8:33
• Well, the OP did dismiss my proposal so glibly... – J. M. is away May 6 '13 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. – Mark McClure May 6 '13 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. – bill s May 6 '13 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}]]. – J. M. is away May 7 '13 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. – Michael E2 May 7 '13 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. – J. M. is away May 7 '13 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. – Michael E2 May 7 '13 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}