﻿ algorithm for generating a Farey series
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  GENERATE FAREY SERIES
C  06/10/07 (Darrell Cox)
C
C  This C program computes the first half of a Farey series of order n (the
C  fractions in a Farey series are complementary around the fraction 1/2,
C  that is, if 1/2 is the i th fraction in the Farey series, then F(i+1) =
C  1 - F(i-1), F(i+2) = 1 - F(i-2), etc.).  Temporary results are pushed on
C  the stack ("STACK").  The size of this array must be at least n+1.  The
C  numerators and denominators of the fractions in the first half of the
C  Farey series are output in the array "R".  The size of this array must
C  be as large as the number of fractions in half a Farey series (about
C  .1520*n*n).
C
C  This algorithm for generating a Farey series uses Haros' theorem
C  that if h/k and h'/k' are successive fractions in a Farey series, then
C  kh' - hk' = 1.  The numerator and denominator of a current fraction
C  ("NUM" and "DEN") and the numerators and denominators of fractions on a
C  stack are maintained with the numerator and denominator of the fraction
C  on the top of the stack being denoted by "STACK(SP,1)" and "STACK(SP,2)"
C  respectively (where "SP" denotes the stack pointer).  The initial values
C  of "NUM" and "DEN" are 1 and N respectively where N is the order of the
C  Farey series.  The initial values of the numerators and denominators of
C  the fractions on the stack are 1 for the numerators and N-1, N-2, N-3,...,
C  2 for the denominators (with N-1 being on the top of the stack).  As the
C  numerators and denominators of the successive fractions in the Farey series
C  are computed, they are stored into the array "R(M,1)" and "R(M,2)" where
C  "M" is the current index into the array.  (Actually, as implemented here,
C  the numerators and denominators of the fractions 1/N, 1/(N-1), 1/(N-2),...,
C  1/J where J=N-[N/2] are first stored in the array "R", the index "M" is set
C  to [N/2]+1, "NUM" and "DEN" are set to 1 and "J" respectively, and the
C  numerators and denominators of the rest of the fractions with a numerator
C  of 1 are pushed on the stack.  The only reason for doing this is that the
C  processing for the first few fractions can be simplified and thus a little
C  execution time saved.)  Since the numerators of the fractions initially on
C  the stack are 1 and the denominators are numerically consecutive and
C  decreasing, an adjacent pair of fractions on the stack, say n/d and n'/d',
C  satisfies d*n' - n*d' = 1.  Also, initially DEN*STACK(SP,1)-NUM*STACK(SP,2)
C  = 1.  If also DEN + STACK(SP,2) > N (using the theorem that the sum of the
C  denominators of two successive fractions in a Farey series is greater than
C  the order of the Farey series), then by Haros' theorem (that two successive
C  fractions h/k and h'/k' in a Farey series satisfy k*h' - h*k' = 1),
C  "STACK(SP,1)" and "STACK(SP,2)" are the numerator and denominator of the
C  next fraction in the Farey series.  (There are infinitely many solutions
C  (x,y) of k*x - h*y = 1, (h,k)=1, but only one solution where x <= N,
C  y <= N, and y + k > N.)  The modulus function is used to find a linear
C  combination of "DEN" and STACK(SP,2)", denoted by "NEWDEN", such that
C  N - NEWDEN < DEN.  If "NEWDEN" is set to "STACK(SP,2) + [(N - STACK(SP,2))/
C  DEN]*DEN", then N - NEWDEN < DEN (since "N - NEWDEN" is "(N - STACK(SP,2))
C  modulus DEN)".  Denote "[(N - STACK(SP,2))/DEN]" (where the brackets denote
C  the greatest integer function) by "J" and denote "STACK(SP,1) + NUM*J"
C  by "NEWNUM".  Adding "DEN*NUM*J - NUM*DEN*J" (= 0) to both sides of DEN*
C  STACK(SP,1) - NUM*STACK(SP,2) = 1 gives DEN*NEWNUM - NUM*NEWDEN = 1 and
C  hence "NEWNUM/NEWDEN" is the next fraction in the Farey series.  (If J=0,
C  "NEWNUM" and "NEWDEN" are "STACK(SP,1)" and "STACK(SP,2)".  These values
C  are popped off the stack and stored in the "R" array.  "NUM" and "DEN"
C  are set to these values.  It is not possible that the stack [even after it
C  has been popped] be empty when "J" equals 0 [this will be elaborated later],
C  so checking the stack pointer is not necessary.)  If J > 1, "NEWNUM" and
C  "NEWDEN" are stored in the "R" array.  If J > 1, [STACK(SP,2) + DEN*I]*
C  [STACK(SP,1) + NUM*(I-1)] - [STACK(SP,1) + NUM*I]*[STACK(SP,2) + DEN*
C  (I-1)] = 1, (I=J,J-1,J-2,...,2), (since DEN*STACK(SP,1) - NUM*STACK(SP,2)
C  = 1).  Let "OLDSP" equal "SP".  If J > 1, the numerators and denominators
C  "STACK(OLDSP,1) + NUM*I" and "STACK(OLDSP,2) + DEN*I", (I=1,2,3,...,(J-1)),
C  are pushed on the stack (so that the fraction "[STACK(OLDSP,1) + NUM*(J-1)]/
C  [STACK(OLDSP,2) + DEN*(J-1)]" is on the top of the stack).  [STACK(OLDSP,2)
C  + DEN]*STACK(OLDSP,1) - [STACK(OLDSP,1) + NUM]*STACK(OLDSP,2) = 1, so any
C  adjacent pair of fractions on the stack, say n/d and n'/d', still satisfies
C  d*n' - n*d' = 1.  If J > 1, DEN <= DEN*(J-1) and hence STACK(OLDSP,2) +
C  DEN*J + DEN < 2*STACK(OLDSP,2) + DEN*J + DEN <= 2*STACK(OLDSP,2) + DEN*J +
C  DEN*(J-1).  Since N - [STACK(OLDSP,2) + DEN*J] < DEN, N < STACK(OLDSP,2) +
C  DEN*J + DEN, and hence N < 2*STACK(OLDSP,2) + DEN*J + DEN*(J-1), that is,
C  N - [STACK(OLDSP,2) + DEN*(J-1)] < STACK(OLDSP,2) + DEN*J.  Therefore if J >
C  1, NEWDEN*STACK(SP,1) - NEWNUM*STACK(SP,2) = 1 and N - STACK(SP,2) < NEWDEN,
C  and hence the fraction on the top of the stack is the next fraction in the
C  Farey series.  (In practice, it is not necessary to push this fraction on
C  the stack in the first place since it is known that it will be the next
C  fraction in the series.  If J=2, it is not necessary to push any fractions
C  on the stack.)  If J=1, NEWNUM = STACK(SP,1) + NUM and NEWDEN = STACK(SP,2)
C  + DEN.  "NEWNUM" and "NEWDEN" are then the numerator and denominator of the
C  next fraction in the Farey series and are stored in the "R" array.  Denote
C  the increase in the denominator by "DELDEN".  If J=1, "DELDEN" equals
C  "STACK(SP,2)" (since NEWDEN - DEN = STACK(SP,2)). Let K = [(N - NEWDEN)/
C  DELDEN] ("K" is the maximum number of times "NEWDEN" can be increased by
C  "DELDEN" and still remain less than or equal to "N").  If J=1 and K > 0,
C  [STACK(SP,2) + DEN*I]*[STACK(SP,1) + NUM*(I+1)] - [STACK(SP,1) + NUM*I]*
C  [STACK(SP,2) + DEN*(I+1)] = 1, (I=1,2,3,...,K), (since DEN*STACK(SP,1) -
C  NUM*STACK(SP,2) = 1).  Since STACK(SP,2) < DEN and N - [STACK(SP,2) + DEN] <
C  DEN, N - [STACK(SP,2) + DEN*(I+1)] < [STACK(SP,2) + DEN*I], (I=1,2,3,...,K),
C  and hence (if K>0) "[STACK(SP,1) + NUM*(I+1)]" and "[STACK(SP,2) + DEN*
C  (I+1)]", (I=1,2,3,...,K) are the numerators and denominators of the next "K"
C  successive fractions in the Farey series.  Now denote "DEN + STACK(SP,2)*
C  (K+1)" by "NEWDEN" and "NUM + STACK(SP,1)*(K+1)" by "NEWNUM".  (Note that if
C  K=0, the definition of "NEWDEN" and "NEWNUM" is the same as before.)  Adding
C  "STACK(SP,2)*(K+1)*STACK(SP,1) - STACK(SP,1)*(K+1)*STACK(SP,2)" (=0) to both
C  sides of DEN*STACK(SP,1) - NUM*STACK(SP,2) = 1 gives NEWDEN*STACK(SP,1) -
C  NEWNUM*STACK(SP,2) = 1.  N - STACK(SP,2) = DEN + [N - (STACK(SP,2) + DEN)] <
C  DEN + STACK(SP,2)*{[(N - (STACK(SP,2) + DEN))/STACK(SP,2)] + 1} = DEN +
C  STACK(SP,2)*(K+1) and hence N - STACK(SP,2) < NEWDEN.  Therefore if J=1, the
C  numerator and denominator of the next fraction in the Farey series are
C  "STACK(SP,1)" and "STACK(SP,2)".  "NUM" is set to "STACK(SP,1)" and "DEN"
C  is set to "STACK(SP,2)" and these values are popped off the stack and stored
C  in the "R" array.  The stack pointer is then checked to see if the fraction
C  1/2 has been reached (the fraction 1/2 can only be reached when J=1, so
C  checking the stack pointer in other instances is unnecessary).  If not,
C  another iteration of the algorithm is performed.  The algorithm can be
C  iterated since the stack is always maintained so that two adjacent fractions
C  on the stack, say n/d and n'/d', satisfy d*n' - n*d' = 1.  The execution
C  time of the algorithm is proportional to the number of fractions in the
C  Farey series since all the operations result in a successive fraction being
C  generated (and the number of operations required to generate a fraction is
C  less than some constant) or a fraction being pushed on the stack (not all
C  the fractions are pushed on the stack and if they are, only once).  Other
C  than a relatively small array of [(n - 1)/2]*2 elements (where the brackets
C  denote the greatest integer function), the algorithm requires no more
C  memory than that required to hold the fractions.  This technique can also
C  be used to compute the fractions in a Farey series between the fractions
C  1/a and 1/b where a>b.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
unsigned int farey(unsigned int N, unsigned int *R,
unsigned int *STACK) {
unsigned int M,NUM,DEN,DELNUM,DELDEN,I,J,K,SP;
//
// STORE THE FIRST FRACTION
//
R[0]=0;
R[1]=1;
//
// STORE THE FRACTIONS BEFORE THE FIRST INCREASE IN DENOMINATOR
//
K=(N/2)+2;
J=N;
for (I=1; I<K; I++) {
R[2*I]=1;
R[2*I+1]=J;
J=J-1;
}
M=K-1;
NUM=R[2*M];
DEN=R[2*M+1];
//
// PUSH THE REMAINING FRACTIONS WITH A NUMERATOR OF 1 ON THE STACK
//
K=N-K+1;
if(K<2) goto L600;
SP=0;
for (I=2; I<=K; I++) {
SP=SP+1;
STACK[2*SP]=1;
STACK[2*SP+1]=I;
}
/***************/
/* BEGIN LOOP  */
/***************/
//
// CHECK FOR AN INCREASE IN DENOMINATOR
//
L100: J=(N-STACK[2*SP+1])/DEN;
if(J!=0) goto L200;
//
// DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK
//
NUM=STACK[2*SP];
DEN=STACK[2*SP+1];
SP=SP-1;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
goto L100;
//
// INCREASE IN DENOMINATOR
//
L200: if(J==1) goto L400;
DELNUM=NUM;
DELDEN=DEN;
NUM=STACK[2*SP]+NUM*J;
DEN=STACK[2*SP+1]+DEN*J;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
NUM=STACK[2*SP]+DELNUM;
DEN=STACK[2*SP+1]+DELDEN;
if(J==2) goto L300;
//
// PUSH FRACTIONS ON THE STACK
//
K=J-2;
for (I=1; I<=K; I++) {
SP=SP+1;
STACK[2*SP]=NUM;
STACK[2*SP+1]=DEN;
NUM=NUM+DELNUM;
DEN=DEN+DELDEN;
}
//
// DECREASE IN DENOMINATOR
//
L300: M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
goto L100;
//
// INCREASE IN DENOMINATOR
//
L400: DELNUM=STACK[2*SP];
DELDEN=STACK[2*SP+1];
NUM=DELNUM+NUM;
DEN=DELDEN+DEN;
K=(N-DEN)/DELDEN;
if(K==0) goto L500;
for (I=1; I<=K; I++) {
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
NUM=NUM+DELNUM;
DEN=DEN+DELDEN;
}
L500: M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
//
// DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK
//
NUM=STACK[2*SP];
DEN=STACK[2*SP+1];
SP=SP-1;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
if(SP>0) goto L100;
/***************/
/* END LOOP    */
/***************/
L600: return(M+1);
}
```