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