/*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 next    C
C  fraction in the series is computed from the current fraction using Haros'  C
C  theorem that if h/k and h'/k' are successive fractions in a Farey series,  C
C  k*h' - h*k' = 1.  Lagrange's algorithm is used to solve the Diophantine    C
C  equation k*x - h*y = 1.  The execution time of the algorithm is            C
C  proportional to n*n*log2(n).                                               C
C                                                                             C
C  The "Q" array is used to store partial quotients.  The numbers generating  C
C  the most partial quotients are two successive Fibonacci numbers.  These    C
C  numbers can be defined by the recurrence relation f(i+1) = f(i) + f(i-1),  C
C  i=1,2,3,..., f(0)=f(1)=1.  To hold the partial quotients generated by the  C
C  j-1 and j th Fibonacci numbers, the dimension of the "Q" array must be at  C
C  least j-1.  The sixteenth Fibonacci number is 1597, so an array size of 20 C
C  is more than adequate for computing Farey series of orders less than 1000. C
C  (The "E" and "F" arrays are used to store the convergents and must be as   C
C  large as the "Q" array.)                                                   C
C                                                                             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC*/
      void lagrange2(unsigned int N, unsigned int D, unsigned int O,
		     unsigned int *ND) {
      unsigned int A,B,R,I,J,Q[400],E[400],F[400];
//
// SOLVE D*X - N*Y = 1 USING LAGRANGE'S ALGORITHM;
// FIND QUOTIENTS USING EUCLID'S GREATEST COMMON DIVISOR ALGORITHM
//
      A=N;
      B=D;
      I=0;
L300: I=I+1;
      Q[I-1]=B/A;
      R=B-Q[I-1]*A;
      if(R==0) goto L400;
      B=A;
      A=R;
      goto L300;
//
// COMPUTE CONVERGENTS
//
L400: A=Q[0]-1;
      B=1;
      if(I==1) goto L600;
      A=A+1;
      if(I==2) goto L600;
      E[0]=A;
      F[0]=B;
      A=A*Q[1]+1;
      B=Q[1];
      if(I==3) goto L500;
      E[1]=A;
      F[1]=B;
      for (J=3; J<I; J++) {
	 A=A*Q[J-1]+E[J-3];
	 B=B*Q[J-1]+F[J-3];
	 E[J-1]=A;
	 F[J-1]=B;
	 }
      if(I==(I/2)*2) goto L600;
//
// SOLUTION OF D*X - N*Y = -1; ADJUST SOLUTION
//
L500: A=D-A;
      B=N-B;
//
// COMPUTE NEXT FRACTION;
// FIND R SUCH THAT O - D < A + D*R <= O
//
L600: R=(O-A)/D;
      N=B+N*R;
      D=A+D*R;
      ND[0]=N;
      ND[1]=D;
      return;
      }