/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C GENERATE FAREY SERIES C
C 06/10/07 (DKC) C
C C
C This subroutine generates the fractions in a Farey series after two given C
C successive fractions in the series. The numerator and denominator of the C
C current fraction are denoted by "NUM" and "DEN" respectively and the C
C numerator and denominator of the previous fraction are denoted by "OLDNUM" C
C and "OLDDEN" respectively. Denote NUM-OLDNUM and DEN-OLDDEN by "DELNUM" C
C and "DELDEN" respectively. If the denominators are increasing (DELDEN>0), C
C the next K=(N-DEN)/DELDEN fractions in the Farey series correspond to the C
C remaining lattice points on the line between (OLDNUM,OLDDEN) and (NUM,DEN) C
C and the numerator and denominator of the fraction following these C
C fractions are DELNUM and DELDEN respectively. If K<>0, these fractions C
C are stored in the Farey series array (R) and the current fraction set to C
C the last fraction (DELNUM/DELDEN). If K<>0, DELDEN is less than OLDDEN, C
C so that the current fraction has been adjusted so that OLDDEN>DEN. If C
C K=0 and DELDEN is less than OLDDEN, the current fraction is set to the C
C last fraction (DELNUM/DELDEN) and the current fraction stored in the R C
C array. If K=0 and OLDDEN<DEN, then OLDDEN is less than the denominator of C
C the previous fraction in the Farey series so that (OLDNUM,OLDDEN) is an C
C inflection point in the lattice of points corresponding to the fractions C
C in the Farey series. In the loop, OLDNUM and OLDDEN are adjusted so that C
C (OLDNUM,OLDDEN) is an inflection point by taking OLDNUM modulus NUM and C
C OLDDEN modulus DEN (computing OLDNUM/NUM is not necessary since OLDNUM/NUM C
C equals OLDDEN/DEN). NUM and DEN are then adjusted so that (NUM,DEN) is a C
C following inflection point by taking NUM modulus OLDNUM (when OLDNUM>1) C
C and DEN modulus OLDDEN. (When OLDNUM>=2, NUM/OLDNUM equals DEN/OLDDEN and C
C when OLDNUM=1, NUM/OLDNUM-1 equals DEN/OLDDEN.) When OLDNUM<2, NUM is set C
C to NUM-(DEN/OLDDEN)*OLDNUM. The fractions corresponding to the lattice C
C points on the line between (NUM,DEN) (before adjustment) and the adjusted C
C (NUM,DEN) value are pushed on the stack beginning with the fraction C
C corresponding to the adjusted (NUM,DEN) value and ending with the fraction C
C corresponding to the point after (NUM,DEN). If OLDNUM=0, the number of C
C fractions to be pushed on the stack is decreased by one to avoid pushing C
C the fraction 1/0 on the stack. The fractions in the Farey series between C
C NUM/DEN and the adjusted NUM/DEN value are then computed using the C
C unadjusted NUM/DEN value and the fractions pushed on the stack. For each C
C iteration of the loop, (OLDNUM,OLDDEN) becomes an inflection point further C
C back in the Farey series and (NUM,DEN) becomes an inflection point further C
C forward in the Farey series so that (OLDNUM,OLDDEN) eventually becomes the C
C first fraction in the Farey series (0,1) and (NUM,DEN) eventually becomes C
C the last fraction in the Farey series (1,1). Since OLDDEN is set to C
C OLDDEN modulus DEN and DEN is set to DEN modulus OLDDEN, (NUM,DEN) and C
C (OLDNUM,OLDDEN) rapidly converge to their final values. For repeated C
C applications of the functions B=B modulus A and A=A modulus B, the slowest C
C converging initial values (A,B) are two consecutive Fibonacci numbers. C
C For Farey series of orders less than 1000, fewer than 17 iterations of the C
C loop would be required. C
C C
C The size of the "STACK" array must be as large as two times the order of C
C the series (N) minus two. The R array must be large enough to hold the C
C fractions in the Farey series. C
C 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);
void farsec(unsigned int N, unsigned int M, unsigned int *R,
unsigned int NUM, unsigned int DEN, unsigned int OLDNUM,
unsigned int OLDDEN, unsigned int *STACK) {
unsigned int I,J,K,SP,DELNUM,DELDEN;
//
// STORE LAST FRACTION AND CURRENT FRACTION
//
R[2*M]=OLDNUM;
R[2*M+1]=OLDDEN;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
if(DEN==1) goto L500;
if(DEN<OLDDEN) goto L300;
//
// INCREASE IN DENOMINATOR
//
DELNUM=NUM-OLDNUM;
DELDEN=DEN-OLDDEN;
K=(N-DEN)/DELDEN;
if(K!=0) goto L100;
if(DELDEN<OLDDEN) goto L200;
goto L400;
L100: for (I=0; I<K; I++) {
NUM=NUM+DELNUM;
DEN=DEN+DELDEN;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
}
L200: NUM=DELNUM;
DEN=DELDEN;
M=M+1;
R[2*M]=NUM;
R[2*M+1]=DEN;
if(DEN==1) goto L500;
//
// DECREASE IN DENOMINATOR
//
/***************/
/* BEGIN LOOP */
/***************/
L300: J=OLDDEN/DEN;
OLDNUM=OLDNUM-J*NUM;
OLDDEN=OLDDEN-J*DEN;
L400: J=DEN/OLDDEN;
if(OLDNUM==0) J=J-1;
//
// PUSH FRACTIONS ON THE STACK
//
NUM=NUM-J*OLDNUM;
DEN=DEN-J*OLDDEN;
SP=0;
for (I=0; I<J; I++) {
SP=SP+1;
STACK[2*SP]=NUM;
STACK[2*SP+1]=DEN;
NUM=NUM+OLDNUM;
DEN=DEN+OLDDEN;
}
M=frank(N,M,R,NUM,DEN,SP,STACK);
NUM=STACK[0];
DEN=STACK[1];
if(DEN!=1) goto L300;
/***************/
/* END LOOP */
/***************/
L500: return;
}