## Implementing the Farey sequence efficiently

6

1

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?

1Hmm... – J. M.'s ennui – 2013-05-06T04:23:51.557

1

I saw the implementation in the link. Isn't there a simpler implementation?

– user76568 – 2013-05-06T04:34:01.837

1Your implementation is certainly simpler. Of course, it's also less efficient. – J. M.'s ennui – 2013-05-06T06:03:52.227

Actually, the OPs code comes from here: http://www.physicsforums.com/showthread.php?t=489620 or here http://demonstrations.wolfram.com/FareySequence/ whichever came first.

– bill s – 2013-05-06T06:21:37.030

@J.M. I'm sure there is a more elegant, more efficient, simpler to read and comprehend, method. Better than all those mentioned, because it is all very simple. I'll give it a go with a rit or 2 later.... :P – user76568 – 2013-05-06T06:34:40.367

Be my guest, then. :) – J. M.'s ennui – 2013-05-06T07:32:59.863

There are many relationships between members of the Farey sequence. For instance, if a/b and c/d are members, then the mediant (a+c)/(b+d) is in the sequence. This property might be easier to exploit than the one you've given. – bill s – 2013-05-06T07:39:39.923

What you said was not correct, but makes you wonder. – user76568 – 2013-05-06T07:45:28.040

The mediant property is well known, see for instance: http://www.cse.iitd.ernet.in/~mcs103480/farey-icip-100525.pdf

– bill s – 2013-05-06T08:04:22.410

@bill, there are restrictions to be imposed on the denominators for the constructed mediant to be a proper member of the $k$th Farey sequence, of course. – J. M.'s ennui – 2013-05-06T08:14:21.177

Right, it holds as long as the denominator is not larger than k -- the DeleteCases in the answer removes those extra ones. – bill s – 2013-05-06T08:15:54.563

@J.M. Some ideas also here http://codegolf.stackexchange.com/questions/1896/code-challenge-farey-sequence-ii

– Dr. belisarius – 2013-05-06T12:03:20.850

@bel, that just counts the number of terms in the Farey sequence. That is most easily done in Mathematica, which has EulerPhi[]. – J. M.'s ennui – 2013-05-06T12:28:12.277

@J.M. Then I misread it. Sorry. – Dr. belisarius – 2013-05-06T12:39:46.767

@J.M. Nevertheless, EulerPhi[] can be handy even for our purposes... Along with another 10 (or so) tricks... (as of now) :) Let's start a most efficient Farey competition! The winner.. uumm.. will be awarded a book about Number Theory.. or... candy. – user76568 – 2013-05-06T12:50:09.023

1I 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 – 2013-05-07T02:35:03.633

@MichaelE2 you are absolutely correct. Your comment is Superior. :) Do you think I should remake the question? So It'll be Superior? This Farey sequence is packed with logic/rules and confinments, up to a point where I think it can be implmented on/by pure logic :) – user76568 – 2013-05-07T07:07:27.680

1Yes, 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 – 2013-05-07T13:04:02.507

6

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.

1This 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.'s ennui – 2013-05-07T19:22:12.270

1Thank 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 – 2013-05-07T19:42:47.310

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 – 2013-05-07T19:52:25.247

@Michael E2: Brilliant. Your recommendations bring the timings down to 2/3 of previous values, much closer to those of farey2[n]. – KennyColnago – 2013-05-07T20:13:37.157

6

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.

1med[v_List] := Total[Numerator[v]]/Total[Denominator[v]] is a bit more compact. – J. M.'s ennui – 2013-05-06T08:29:37.633

That's good, but you know, this whole approach doesn't seem any faster than the simple way! – bill s – 2013-05-06T08:33:25.520

Well, the OP did dismiss my proposal so glibly... – J. M.'s ennui – 2013-05-06T08:46:33.877

I'll do some performance testing on this marble. Though I'm working on one my self.... – user76568 – 2013-05-06T08:57:56.937

This implementation is slowest thus far as far as I can check – user76568 – 2013-05-06T09:25:03.480

2@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 – 2013-05-06T13:32:03.610

@MarkMcClure , I guess you are right. Obviously, the performance is more than "enough" and the recovery of Stern-Brocot tree obtained by merely removing the head is quite nice. Though, I know there are many many mathematical facts that could be harnessed, for the sake of elegance, simplicity, readability (and performance!!!) among others, while still maintaining functional form.. I actually am inclined to implement it... But I don't have time now. – user76568 – 2013-05-06T14:17:21.980

1@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 – 2013-05-06T16:00:52.573

@bills Well, obviously. Regarding that, I don't think It's TOO weird, that initially I was hoping for a step up in performance. By the way, it isn't much slower... It's same'ish-- :) – user76568 – 2013-05-06T16:04:41.810

@bills Aaaaaaaaaand, I do VERY MUCH seriously thank you for the time you spent on my question, and the new tricks I learned from your implementation. :) – user76568 – 2013-05-06T16:11:05.140

5

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


1Any 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.'s ennui – 2013-05-07T14:38:46.007

@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 – 2013-05-07T14:50:26.653

"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.'s ennui – 2013-05-07T16:19:17.650

@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 – 2013-05-07T18:29:24.403

2

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}