/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  GENERATE FAREY SERIES                                                      C
C  06/11/07 (DKC)							      C
C                                                                             C
C  This subroutine computes the fractions in a Farey series after two given   C
C  successive fractions.  The algorithm uses the following theorems;          C
C                                                                             C
C  (1)  Suppose h/k and h'/k' are successive fractions in a Farey series of   C
C       order n where k' > k and denote [(n - k')/(k' - k)] by l.  If l <> 0, C
C       the next l fractions in the Farey series correspond to the remaining  C
C       lattice points on the line through (h,k) and (h',k') and are          C
C       (h' + (h' - h)*i)/(k' + (k' - k)*i), i=1,2,3,...,l.                   C
C                                                                             C
C  (2)  If h/k and h'/k' are successive fractions in a Farey series and k' >  C
C       k, the fraction in the Farey series after the fractions corresponding C
C       to the lattice points on the line through (h,k) and (h',k') is        C
C       (h' - h)/(k' - k).                                                    C
C                                                                             C
C  (3)  Suppose h/k and h'/k' are successive fractions in a Farey series of   C
C       order n where k > k' and denote [{(2k' - n - 1)/(k - k')}/2] by l.    C
C       ({A} denotes the smallest integer greater than  A.)  If l > 0, the    C
C       next l fractions in the Farey series correspond to the next l         C
C       lattice points on the line through (h,k) and (h',k') and are          C
C       (h' - (h - h')*i)/(k' - (k - k')*i), i=1,2,3,...,l.                   C
C                                                                             C
C  (4)  The fraction after the successive fractions h/k and h'/k' in a Farey  C
C       series of order n is (j*h' - h)/(j*k' - k) where j=[(n + k)/k'].      C
C                                                                             C
C  If the denominators of the given successive fractions are increasing, the  C
C  first two theorems can be used to find the next pair of successive         C
C  fractions where the denominators are decreasing.  The last two theorems    C
C  can then be used to find the next pair of successive fractions where the   C
C  denominators are increasing, etc.  Similar logic applies when the          C
C  denominators of the given successive fractions are decreasing.             C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
       void harfar(unsigned int N, unsigned int M, unsigned int *R,
		   unsigned int NUM, unsigned int DEN, unsigned int OLDNUM,
		   unsigned int OLDDEN) {
       unsigned int TNUM,TDEN,II,K;
       int I;
//
// STORE CURRENT FRACTIONS
//
       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 L600;
       TNUM=OLDNUM;
       TDEN=OLDDEN;
       if(DEN>OLDDEN) goto L400;
//
// DECREASE IN DENOMINATOR
//
       I=2*DEN-N;
       if(I<0) goto L200;
       OLDNUM=OLDNUM-NUM;
       OLDDEN=OLDDEN-DEN;
L100:  K=(I+OLDDEN-1)/OLDDEN;
       K=K/2;
       if(K==0) goto L200;
       for (II=0; II<K; II++) {
	  NUM=NUM-OLDNUM;
	  DEN=DEN-OLDDEN;
	  M=M+1;
	  R[2*M]=NUM;
	  R[2*M+1]=DEN;
	  }
//
// INCREASE IN DENOMINATOR
//
L200:  TNUM=NUM;
       TDEN=DEN;
L300:  K=(N+OLDDEN)/DEN;
       NUM=K*NUM-OLDNUM;
       DEN=K*DEN-OLDDEN;
       M=M+1;
       R[2*M]=NUM;
       R[2*M+1]=DEN;
//
L400:  OLDNUM=NUM-TNUM;
       OLDDEN=DEN-TDEN;
       K=(N-DEN)/OLDDEN;
       if(K==0) goto L500;
       for (II=0; II<K; II++) {
	  NUM=NUM+OLDNUM;
	  DEN=DEN+OLDDEN;
	  M=M+1;
	  R[2*M]=NUM;
	  R[2*M+1]=DEN;
	  }
L500:  TNUM=OLDNUM;
       TDEN=OLDDEN;
       OLDNUM=NUM-OLDNUM;
       OLDDEN=DEN-OLDDEN;
       NUM=TNUM;
       DEN=TDEN;
       M=M+1;
       R[2*M]=NUM;
       R[2*M+1]=DEN;
       I=2*DEN-N;
       if(I>0) goto L100;
       if(DEN!=1) goto L300;
//
L600:  return;
       }