Lagrange's Algorithm for Solving the Diophantine Equation kx-hy=1
Created | Updated Dec 29, 2008
To find integers p1 and q1 satisfying pq1-qp1=±1 where p and q are relatively prime, Lagrange reduced p/q to a continued fraction. A series of fractions converging towards p/q, alternately less than and greater than p/q, is obtained. If p1 is taken to equal the numerator and q1 taken to equal the denominator of the convergent immediately preceding p/q, then pq1-qp1=1 if p1/q1<p/q, or -1 if p1/q1>p/q. Euclid's greatest common divisor algorithm is used to compute the partial quotients of the continued fraction. An embellishment to the algorithm is that a solution of pq1-qp1=1 is given if the number of partial quotients is even. Note that if pq1-qp1=-1, then p(q-q1)-q(p-p1)=1 so that a solution of pq1-qp1=1 is given by setting q1 to q-q1 and p1 to p-p1. Few details of Lagrange's algorithm are given in the mathematical literature, so the algorithm is reproduced here. The following C subroutine generates a Farey series using Lagrange's algorithm;
/* GENERATE FAREY SERIES
The Farey series Fn of order n is the ascending series of irreducible fractions between 0 and 1 whose denominators do not exceed n. The next fraction in the series is computed from the current fraction using Haros' theorem that if h/k and h'/k' are successive fractions in a Farey series, then kh'-hk'=1. Lagrange's algorithm is used to solve the Diophantine equation kx-hy=1. The execution time of the algorithm is proportional to n2log2(n). The "Q" array is used to store the partial quotients. The numbers generating the most partial quotients are two successive Fibonacci numbers. These numbers can be defined by the recurrence relation f(i+1)=f(i)+f(i-1), i=1, 2, 3, ..., f(0)=f(1)=1. To hold the partial quotients generated by the (j-1)th and jth Fibonacci numbers, the dimension of the Q array must be at least j-1. The sixteenth Fibonacci number is 1597, so an array size of 20 is more than adequate for computing Farey series of orders less than 1000. The "E" and "F" arrays are used to store the convergents and must be as large as the Q array. The "S" array must be large enough to hold the Farey series. The order of the Farey series is denoted by "O". */
unsigned int farey(unsigned int O, unsigned int *S) {
unsigned int M,N,D,A,B,R,I,J,Q[20],E[20],F[20];
//
// INITIAL FRACTION
//
M=0; // output array index
S[2*M]=0; // store fraction
S[2*M+1]=1;
N=1; // numerator
D=O; // denominator
M=1;
S[2*M]=N;
S[2*M+1]=D;
if (D==1) goto L800;
if (D==2) goto L700;
//
// SOLVE DX-NY=1 USING LAGRANGE'S ALGORITHM
// FIND QUOTIENTS USING EUCLID'S GREATEST COMMON DIVISOR ALGORITHM
//
L200: 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;
M=M+1;
S[2*M]=N;
S[2*M+1]=D;
if (D!=2) goto L200;
//
// COMPUTE THE FRACTIONS AFTER 1/2
//
L700: J=M;
for (I=1; I≤J; I++) {
S[2*(M+I)]=S[2*(M-I)+1]-S[2*(M-I)];
S[2*(M+I)+1]=S[2*(M-I)+1];
}
M=M+J;
//
L800: return(M+1);
}
Mém, Acad. Berlin, 23, année 1767, 1769, 7; Oeuvres, 2, 1868, 386-8.
/* GENERATE FAREY SERIES
The Farey series Fn of order n is the ascending series of irreducible fractions between 0 and 1 whose denominators do not exceed n. The next fraction in the series is computed from the current fraction using Haros' theorem that if h/k and h'/k' are successive fractions in a Farey series, then kh'-hk'=1. Lagrange's algorithm is used to solve the Diophantine equation kx-hy=1. The execution time of the algorithm is proportional to n2log2(n). The "Q" array is used to store the partial quotients. The numbers generating the most partial quotients are two successive Fibonacci numbers. These numbers can be defined by the recurrence relation f(i+1)=f(i)+f(i-1), i=1, 2, 3, ..., f(0)=f(1)=1. To hold the partial quotients generated by the (j-1)th and jth Fibonacci numbers, the dimension of the Q array must be at least j-1. The sixteenth Fibonacci number is 1597, so an array size of 20 is more than adequate for computing Farey series of orders less than 1000. The "E" and "F" arrays are used to store the convergents and must be as large as the Q array. The "S" array must be large enough to hold the Farey series. The order of the Farey series is denoted by "O". */
unsigned int farey(unsigned int O, unsigned int *S) {
unsigned int M,N,D,A,B,R,I,J,Q[20],E[20],F[20];
//
// INITIAL FRACTION
//
M=0; // output array index
S[2*M]=0; // store fraction
S[2*M+1]=1;
N=1; // numerator
D=O; // denominator
M=1;
S[2*M]=N;
S[2*M+1]=D;
if (D==1) goto L800;
if (D==2) goto L700;
//
// SOLVE DX-NY=1 USING LAGRANGE'S ALGORITHM
// FIND QUOTIENTS USING EUCLID'S GREATEST COMMON DIVISOR ALGORITHM
//
L200: 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;
M=M+1;
S[2*M]=N;
S[2*M+1]=D;
if (D!=2) goto L200;
//
// COMPUTE THE FRACTIONS AFTER 1/2
//
L700: J=M;
for (I=1; I≤J; I++) {
S[2*(M+I)]=S[2*(M-I)+1]-S[2*(M-I)];
S[2*(M+I)+1]=S[2*(M-I)+1];
}
M=M+J;
//
L800: return(M+1);
}
References
J. L. LagrangeMém, Acad. Berlin, 23, année 1767, 1769, 7; Oeuvres, 2, 1868, 386-8.
