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