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