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