﻿ algorithm for generating a Farey series
```/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  GENERATE FAREY SERIES
C  06/10/07 (DKC)
C
C  This subroutine generates the fractions in a Farey series segment after
C  a given initial fraction.  The numerator and denominator of a current
C  fraction ("NUM" and "DEN") and the numerators and denominators of fractions
C  on a stack are maintained with the numerator and denominator of the
C  fraction on the top of the stack being denoted by "STACK(SP,1)" and
C  "STACK(SP,2)" respectively (where "SP" denotes the stack pointer).  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.  An adjacent pair of fractions
C  initially on the stack, say n/d and n'/d', must satisfy d*n' - n*d' = 1.
C  Also, initially DEN*STACK(SP,1) - NUM*STACK(SP,2) must equal 1.  If
C  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 then set to these values.)  If J > 1, NEWNUM and NEWDEN are
C  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 last
C  fraction has been reached.  If not, another iteration of the algorithm is
C  performed.  The algorithm can be iterated since the stack is always
C  maintained so that two adjacent fractions on the stack, say n/d and n'/d',
C  satisfy d*n' - n*d' = 1.  The execution time of the algorithm is
C  proportional to the number of fractions in the Farey series since all the
C  operations result in a successive fraction being generated (and the number
C  of operations required to generate a fraction is less than some constant)
C  or a fraction being pushed on the stack (not all the fractions are pushed
C  on the stack and if they are, only once).  Other than a relatively small
C  array of (n - 1)*2 elements, the algorithm requires no more memory than
C  that required to hold the fractions.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
unsigned int frank(unsigned int N, unsigned int M, unsigned int *R,
unsigned int NUM, unsigned int DEN, unsigned int SP,
unsigned int *STACK) {
unsigned int I,J,K,DELNUM,DELDEN;
//
// CHECK FOR AN INCREASE IN DENOMINATOR
//
L100: J=(N-STACK[2*SP+1])/DEN;
if(J==0) goto L500;
//
// INCREASE IN DENOMINATOR
//
if(J==1) goto L300;
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 L200;
//
// PUSH FRACTIONS ON THE STACK
//
K=J-2;
for (I=0; I<K; I++) {
SP=SP+1;
STACK[2*SP]=NUM;
STACK[2*SP+1]=DEN;
NUM=NUM+DELNUM;
DEN=DEN+DELDEN;
}
//
// DECREASE IN DENOMINATOR
//
L200: M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
goto L100;
//
// INCREASE IN DENOMINATOR
//
L300: DELNUM=STACK[2*SP];
DELDEN=STACK[2*SP+1];
NUM=DELNUM+NUM;
DEN=DELDEN+DEN;
K=(N-DEN)/DELDEN;
if(K==0) goto L400;
for (I=0; I<K; I++) {
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
NUM=NUM+DELNUM;
DEN=DEN+DELDEN;
}
L400: M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
//
// DECREASE IN DENOMINATOR, POP A FRACTION OFF THE STACK
//
L500: 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;
//
STACK[0]=NUM;
STACK[1]=DEN;
return(M);
}
```