Implementing the Farey sequence efficiently

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

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}

Answer 5

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.