/*CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                             C
C  GENERATE FAREY SERIES                                                      C
C  06/11/07 (DKC)							      C
C                                                                             C
C  The Farey series Fn of order n is the ascending series of irreducible      C
C  fractions between 0 and 1 whose denominators do not exceed n.  The         C
C  fractions in the series are generated using the theorem that if h/k,       C
C  h'/k', and h''/k'' are three successive fractions in a Farey series, then  C
C  h'/k' = (h + h'')/(k + k'').  The fraction after two successive fractions  C
C  h/k and h'/k' in the series is then (j*h' - h)/(j*k' - k) where j is some  C
C  positive integer.  Using the theorem that the sum of the denominators of   C
C  successive fractions in a Farey series is greater than the order of the    C
C  series gives j = [(n + k)/k'].                                             C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
void add64(unsigned int *A, unsigned int *B);
void sub64(unsigned int *A, unsigned int *B);
void mul64_64(unsigned int a0, unsigned int a2, unsigned int m0,
	      unsigned int m2, unsigned int *product);
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2,
	       unsigned int a3, unsigned int *quotient, unsigned int d2,
	       unsigned int d3);
void haros9(unsigned *N, unsigned int *H, unsigned int *K, unsigned int *HP,
	    unsigned int *KP, unsigned int *D, unsigned int *ND) {
unsigned int HPP[2],KPP[2],M[2],T[4],Q[4];
M[0]=0;
M[1]=1;
//
// FIND FRACTIONS IN FAREY SERIES FOLLOWING H/K, HP/KP
//
L100: T[2]=K[0];
      T[3]=K[1];
      add64(N,&T[2]);
      T[0]=0;
      T[1]=0;
      div128_64(T[0],T[1],T[2],T[3],Q,KP[0],KP[1]);
      mul64_64(Q[2],Q[3],HP[0],HP[1],T);
      sub64(&T[2],H);
      HPP[0]=H[0];
      HPP[1]=H[1];
      mul64_64(Q[2],Q[3],KP[0],KP[1],T);
      sub64(&T[2],K);
      KPP[0]=K[0];
      KPP[1]=K[1];
      H[0]=HP[0];
      H[1]=HP[1];
      K[0]=KP[0];
      K[1]=KP[1];
      HP[0]=HPP[0];
      HP[1]=HPP[1];
      KP[0]=KPP[0];
      KP[1]=KPP[1];
      M[1]=M[1]+1;
      if (M[1]==0)
	 M[0]=M[0]+1;
      if ((KP[0]!=D[0])||(KP[1]!=D[1]))
	 goto L100;
//    J=(N+K)/KP;
//    HPP=J*HP-H;
//    KPP=J*KP-K;
//    H=HP;
//    K=KP;
//    HP=HPP;
//    KP=KPP;
//    if(KP!=D) goto L100;
//
//
ND[0]=M[0];
ND[1]=M[1];
T[2]=K[0];
T[3]=K[1];
add64(N,&T[2]);
T[0]=0;
T[1]=0;
div128_64(T[0],T[1],T[2],T[3],Q,KP[0],KP[1]);
mul64_64(Q[2],Q[3],HP[0],HP[1],T);
sub64(&T[2],H);
ND[2]=H[0];
ND[3]=H[1];
mul64_64(Q[2],Q[3],KP[0],KP[1],T);
sub64(&T[2],K);
ND[4]=K[0];
ND[5]=K[1];
//    J=(N+K)/KP;
//    HPP=J*HP-H;
//    KPP=J*KP-K;
return;
}