The purpose of this article is to show how to solve the Diophantine Equation Ax^{2} + Bxy + Cy^{2} + Dx + Ey + F = 0. The term Diophantine Equation means that the solutions (x, y) should be integer numbers. For example, the equation 4y^{2}  25y + 10 = 0 has solutions given by the horizontal line y = 2.5, but since 2.5 is not an integer number, we will say that the equation has no solutions.
There are several cases that depend on the values of A, B and C. The names are taken from the figures represented by the equation in the plane xy: a line, an ellipse, a parabola or a hyperbola (or two lines). These figures are the set of real solutions. In our situation, the set of solutions are represented by isolated point/s in the plane xy.
If F is multiple of g we can divide all three coefficients by g thus obtaining:
dx + ey = f (where d=D/g, e=E/g and f=F/g).
We will use now the Extended Euclid's Algorithm that can be used to find integers u' and v' such that uu'+vv' = ±gcd(u, v) (where the sign depends on the sign of d and e and the exact implementation of the algorithm).
We can let u = d, v = e. Once the values of u' and v' are found so that du'+ev' = ±1 (since gcd(d, e)=1) we can multiply the equation by f to obtain:
du'+ev' = ±1 => d(±fu')+e(±fv') = f => d(±fu')+det+e(±fv')det = f => d(et±fu')+e(dt±fv') = f
So the general solution set is:
x = et ± fu' y = dt ± fv' t = any integer 

gcd(D, E) = gcd(10, 84) = 2. Since the constant term is also multiple of 2, we will divide the equation by this gcd.
The equation is now: 5x + 42y + 8 = 0.
Now we must apply the Generalized Euclidean algorithm:
Multiplying the last equation by F = 8 we obtain:
(136) * 5 + 16 * 42 = 8
Adding and subtracting de t = 5 * 42 t we obtain:
(136 + 42 t) * 5 + (16  5 t) * 42 = 8
So, the solution is given by the set:
x = 136 + 42 t y = 16  5 t 

where t is any integer number.
Bxy + Dx + Ey + F  =  0 

Bxy + Dx + Ey  =  F 
B^{2}xy + BDx + BEy  =  BF 
B^{2}xy + BDx + BEy + DE  =  DE  BF 
(Bx + E) (By + D)  =  DE  BF 
There are two cases: DE  BF = 0 (two lines parallel to x and y axes respectively) and DE  BF 0 (a hyperbola whose asymptotes are parallel to x and y axes).
In the first case a necessary condition to have solutions occurs when one of the parentheses equal zero, i.e., Bx + E = 0 or By + D = 0. Since B 0, we have solutions for:


In the second case the values of x and y are found by finding all divisors of DE  BF. Let d_{1}, d_{2}, ..., d_{n} be the set of divisors of DE  BF.So,
Bx + E  =  d_{i}  By + D  =  (DE  BF) / d_{i}  

Bx  =  d_{i}  E  By  =  (DE  BF) / d_{i}  D 


Example 2: Solve 2xy + 5x + 56y + 7 = 0.
In this case the divisors of DE  BF = 5*56  2*7 = 266 are: ±1, ±2, ±7, ±14, ±19, ±38, ±133, ±266.
Since (2x + 56) (2y + 5) = 266 we obtain:
The only 8 solutions to the requested equations are marked above in red.
Operating with the original quadratic equation:
Ax^{2} + Bxy + Cy^{2} + Dx + Ey + F = 0
Cy^{2} + (Bx + E)y + (Ax^{2} + Dx + F) = 0
y =  (Bx + E) ± sqrt[(Bx + E)^{2}  4C(Ax^{2} + Dx + F)] 2C  (*) 

For any value of x there will be two values of y except at the left and right extremes of the ellipse. In this case there will be only one value of y. To determine the location of the left and right extremes we should equal the square root to zero, so the previous expression returns only one value of y.
(Bx + E)^{2}  4C(Ax^{2} + Dx + F) = 0
(B^{2}  4AC)x^{2} + 2(BE  2CD)x + (E^{2}  4CF) = 0
So the values of x should be between the roots of this equation. If the roots are not real, there will be no solutions to the original equation, else, all integer values of x should be replaced in equation (*) in order to find an integer value of y.
Example 3: Solve 42x^{2} + 8xy + 15y^{2} + 23x + 17y  4915 = 0.
Since B^{2}  4AC = 8^{2}  4*42*15 = 2456 < 0 the equation is elliptical.
The values of x should be between the roots of (B^{2}  4AC)x^{2} + 2(BE  2CD)x + (E^{2}  4CF) = 2456x^{2}  1108x + 295189 = 0. The roots equal 11.19... and 10.74..., so we have to check the values of x from 11 to 10.
The only value of x that replaced in (*) makes y integer occurs for x = 11, where y = 1, therefore this is the only solution to this problem.
Since b^{2} = 4ac is positive, we can choose g with the same sign of A. In this way a and c will be positive (or one of them zero).
The expression b^{2}  4ac = 0 implies that b^{2}/4 = ac. Since gcd(a,c) = 1, both a and c are perfect squares.
Multiplying the original equation by :
g(ax^{2} + bxy + cy^{2}) + Dx + Ey + F = 0
g(x + y)^{2} + Dx + Ey + F = 0
where for the sign of B/A is taken.
Adding and subtracting Dy:
g(x + y)^{2} + D(x + y)  Dy + Ey + F = 0
Let u = x + y: (i)
gu^{2} + Du + (E  D)y + F = 0
(D  E)y = gu^{2} + Du + F (ii)
There are two cases: D  E = 0 (two parallel lines) or D  E 0 (a parabola).
In the first case, D  E = 0.
From (ii): gu^{2} + Du + F = 0
Since x and y should be integer numbers, the equation (i) implies that the number u (the root of the above equation) should be also integer. Let u_{1} and u_{2} be the roots of the above equation.
From (i) we have: x + y  u_{1} = 0 and x + y  u_{2} = 0 which can be solved with the methods for the linear equation.
In the second case, gu^{2} + Du + F should be multiple of D  E.
Let u_{0}, u_{1},... the values of u in the range 0 <= u < D  E for which the above condition holds.
So u = u_{i} + (D  E)t, where t is any integer number. (iii)
Replacing (iii) in (ii):
(D  E)y = g[u_{i} + (D  E)t]^{2} + D[u_{i} + (D  E)t] + F
y = g(D  E)t^{2} + (D + 2gu_{i})t +  gu_{i}^{2} + Du_{i} + F D  E 
From (i) and (iii):
u = x + y = u_{i} + (D  E)t
x  =  g(E  D)t^{2} + (D  E  2gu_{i}  D)t +  
+ 

x  =  g(E  D)t^{2} + ( E  2gu_{i})t +  
+ 

x  =  g(E  D)t^{2} + ( E  2gu_{i})t   
 


Example 4: Find the solutions for 8 x^{2}  24 xy + 18 y^{2} + 5x + 7y + 16 = 0
We have to calculate the values g, a, c, , , D  E and gu^{2} + Du + F.
g = gcd(8, 18) = 2
a = 8/2 = 4
c = 18/2 = 9
= 2
= 3 (since B/A = 24/8 < 0)
D  E = 3 * 5  2 * 7 = 29 (second case)
gu^{2} + Du + F = 4u^{2} + 5u + 32
We have to determine the values of u in the range 0 <= u < 29 for which 4u^{2} + 5u + 32 is multiple of 29.
The values of u are: u_{0} = 2 and u_{1} = 4.
For u_{0} = 2:
x = 174 t^{2}  17 t  2 y = 116 t^{2}  21 t  2 
For u_{0} = 4:
x = 174 t^{2}  41 t  4 y = 116 t^{2}  37 t  4 
Ax^{2} + Bxy + Cy^{2} = F
Multiplying by 4A:
4A^{2}x^{2} + 4ABxy + 4ACy^{2} = 4AF
4A^{2}x^{2} + 4ABxy + B^{2}y^{2}  B^{2}y^{2} + 4ACy^{2} = 4AF
(2Ax + By)^{2}  (B^{2}  4AC)y^{2} = 4AF
This can be interpreted as a difference of squares:
(2Ax + By + sqrt(B^{2}  4AC) y) (2Ax + By  sqrt(B^{2}  4AC) y) = 4AF
(2Ax + (B + sqrt(B^{2}  4AC))y) (2Ax + (B  sqrt(B^{2}  4AC))y) = 4AF
Since 4AF = 0, the condition to have more solutions is that B^{2}  4AC should be a perfect square.
Now the same method used for the linear equation (since the equation are represented by two lines in the plane xy intersecting at the point (0, 0)) can be used in order to find the solutions.
If F 0 and B^{2}  4AC = k^{2} for some integer k, the parentheses in the equation above should be factors of 4AF.
Let u_{1}, u_{2},... be the positive and negative divisors of 4AF.
Then we have the following set of two linear equations in two unknowns:
2Ax + (B+k)y = u_{i}
2Ax + (Bk)y = 4AF/u_{i}
So we have:
y = (u_{i} + 4AF/u_{i}) / (2k)
x = (u_{i}  (B+k)y)
We should discard the values of u_{i} that makes x or y noninteger.
If F is not a multiple of gcd(A, B, C), the equation has no solutions, otherwise we can divide all coefficients of the equation by this gcd.
If 4 F^{2} < B^{2}  4AC, the solutions of the equation will be amongst the convergents of the continued fraction of the roots of the equation At^{2} + Bt + C = 0.
The continued fraction expansion of a quadratic irrationality is periodic. Since B^{2}  4AC is not a perfect square the number of solutions will be infinite or none.
In the other hand, if 4 F^{2} >= B^{2}  4AC solutions can be obtained as follows:
Let G = gcd(x,y), x = Gu and y = Gv.
The original equation is then: AG^{2}u^{2} + BG^{2}uv + CG^{2}v^{2} + F = 0, so F will be multiple of G^{2}.
Dividing the equation by G^{2}:
Au^{2} + Buv + Cv^{2} + F/G^{2} = 0 (1).
Once the values of u and v are found, we can easily determine x = Gu and y = Gv.
So we can assume that gcd(x,y) = 1.
Let x = sy  Fz (2).
Replacing in the original equation:
A(sy  Fz)^{2} + B(sy  Fz)y + Cy^{2} + F = 0
As^{2}y^{2}  2AFsyz + AF^{2}z^{2} + Bsy^{2}  BFyz + Cy^{2} = F
(As^{2} + Bs + C) y^{2} + (2As  B)Fyz + AF^{2}z^{2} =  F
Dividing by F:
(As^{2} + Bs + C) y^{2} / F + (2As + B)yz  AFz^{2} = 1 (3)
Now we must determine the values of s between 0 and F  1 such that As^{2} + Bs + C 0 (mod F). Once the values of y and z are found using continued fraction expansions of the roots of (As^{2} + Bs + C) t^{2} / F + (2As + B)t  AF = 0, the value of x is found by (2). If no solutions are found amongst the convergents, there will be no solutions to (1).
If the original equation has solutions, there should be a solution to the previous congruence, except when gcd(A,B,F) > 1. In this case, if gcd(B,C,F) = 1 we should make the substitution y = sx  Fz (4), so replacing in the original equation:
Ax^{2} + Bx(sx  Fz) + C(sx  Fz)^{2} + F = 0
Ax^{2} + Bsx^{2}  BFxz + Cs^{2}x^{2}  2CFsxz + CF^{2}z^{2} = F
(Cs^{2} + Bs + A) x^{2} + (2Cs  B)Fxz + CF^{2}z^{2} =  F
Dividing by F:
(Cs^{2} + Bs + A) x^{2} / F + (2Cs + B)xz  CFz^{2} = 1 (5).
Now we must determine the values of s between 0 and F  1 such that Cs^{2} + Bs + A 0 (mod F). Once the values of x and z are found using continued fraction expansions of the roots of (Cs^{2} + Bs + A) t^{2} / F + (2Cs + B)t  CF = 0, the value of y is found by (4). If no solutions are found amongst the convergents, there will be no solutions to (1).
The equations (4) and (5) have no solutions when both gcd(A,B,F) and gcd(B,C,F) are greater than 1. In this case we will use the following approach:
Let i, j, m and n be four integers such that in  jm = 1 (6).
If x = iX + jY and y = mX + nY (7) we obtain X = nx  jy and Y = mx + iy (8).
Since the transformation is reversible, we can convert any (x,y) to (X,Y) and vice versa. So we will work with (X,Y) and with these solutions will can compute the values of (x,y) that satisfies the original equation.
Ax^{2} + Bxy + Cy^{2} =where:
= A(iX+jY)^{2} + B(iX+jY)(mX+nY) + C(mX+nY)^{2} =
= aX^{2} + bXY + cY^{2}
a = Ai^{2} + Bim + Cm^{2} (9)So we have to find the values of i and m such that a = Ai^{2} + Bim + Cm^{2} is relatively prime to F.
b = 2Aij + Bin + Bjm + 2Cmn (10)
c = Aj^{2} + Bjn + Cn^{2} (11)
Since gcd(C, F) > 1 we have gcd(Ai^{2} + Bim + Cm^{2}, C) = 1, so gcd(i, C) = 1 and gcd(Ai+Bm, C) = 1.
Since gcd(A, F) > 1 we have gcd(Ai^{2} + Bim + Cm^{2}, A) = 1, so gcd(m, A) = 1 and gcd(Bi+Cm, A) = 1.
From (6), gcd(i, m) = 1.
If F 0 (mod p) (p prime):
A  B  C  i, m  Examples 

A 0  B 0  C 0  Not applicable (gcd(A, B, C) = 1)  
A 0  B 0  C 0  m 0  i 0, m 1 
A 0  B 0  C 0  i 0, m 0  i 1, m 1 
A 0  B 0  C 0  m 0, i Cm/B  i 1C, m B 
A 0  B 0  C 0  i 0  i 1, m 0 
A 0  B 0  C 0  i 0 or m 0  i 1, m 1 
A 0  B 0  C 0  i 0, m Ai/B  i B, m 1A 
A 0  B 0  C 0  i 0 or m 0  i 1, m 1 
While it is possible to generate the values of i and m from their values modulo different primes, it is very tedious and it is not necessary, because from the table above, almost all values of i and m can be used. So it is better to use the following pseudocode in order to find both values:
for i=0 to F1 for m=0 to i+1 if gcd(i, m) = 1 and gcd(Ai^{2} + Bim + Cm^{2}, F) = 1, end. next m next i
With the values of i and m just found, we can compute the values of j and n from (6) using the methods for the linear equation. Then we have to compute a, b and c using (9), (10) and (11), from which the set of solutions (X,Y) can be found. With the formula (7) we can find the set of solutions (x,y).
Credits: This method was emailed to me by Iain Davidson.
Example 5: Find some solutions for 18 x^{2} + 41 xy + 19 y^{2}  24 = 0
First of all we must determine the gcd of all coefficients but the constant term, that is: gcd(18, 41, 19) = 1.
Dividing the equation by the greatest common divisor we obtain:
18 x^{2} + 41 xy + 19 y^{2}  24 = 0
Let x = sy  fz, so [(as^{2} + bs + c)/f]y^{2} + (2sa + b)yz  afz^{2} = 1.
So 18 s^{2} + 41 s + 19 should be multiple of 24.
This holds for s = 19.
We have to find the continued fraction expansion of the roots of 304 t^{2} + 725 t + 432 = 0, that is, (sqrt(313)  725) /608
The continued fraction expansion is:
2 + //1, 5, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2, 17, 2, 1//
where the periodic part is marked in bold (the period has 19 coefficients).
The following table shows how the values of Y_{0} and Z_{0} are found (the third column are the values for 304 y^{2} + 725 yz + 432 z^{2}):
c_{n}  y_{n}  z_{n}  

1  0  
2  2  1  198 
1  1  1  11 
5  7  6  2 
8  57  49  3 
5  292  251  12 
1  349  300  4 
3  1339  1151  9 
1  1688  1451  8 
1  3027  2602  6 
2  7742  6655  6 
2  18511  15912  8 
1  26253  22567  9 
1  44764  38479  4 
3  160545  138004  12 
1  205309  176483  3 
5  1 187090  1 020419  2 
8  9 702029  8 339835  11 
1  10 889119  9 360254  6 
2  31 480267  27 060343  1 
17  546 053658  469 386085  6 
2  1123 587583  965 832513  11 
1  1669 641241  1435 218598  2 
8  14480 717511  12447 581297  3 
5  74073 228796  63673 125083  12 
1  88553 946307  76120 706380  4 
3  339735 067717  292035 244223  9 
1  428289 014024  368155 950603  8 
1  768024 081741  660191 194826  6 
2  1 964337 177506  1 688538 340255  6 
2  4 696698 436753  4 037267 875336  8 
1  6 661035 614259  5 725806 215591  9 
1  11 357734 051012  9 763074 090927  4 
3  40 734237 767295  35 015028 488372  12 
1  52 091971 818307  44 778102 579299  3 
5  301 194096 858830  258 905541 384867  2 
8  2461 644746 688947  2116 022433 658235  11 
1  2762 838843 547777  2374 927975 043102  6 
2  7987 322433 784501  6865 878383 744439  1 
17  138547 320217 884294  119094 860498 698565  6 
2  285081 962869 553089  245055 599381 141569  11 
1  423629 283087 437383  364150 459879 840134  2 
Notice that y_{n} = c_{n} y_{n1} + y_{n2} and z_{n} = c_{n} z_{n1} + z_{n2}.
The signs in the third column are alternated, so the numbers will repeat after an even number of convergents. Therefore two entire periods should be considered if the period length is odd. If it is even, only one period should be considered. With these solutions and the recurrence relation to be developed in the next section we can find all solutions of the homogeneous equation.
Y_{0} = 7987 322433 784501 (16 digits)
Z_{0} = 6865 878383 744439 (16 digits)
Since X_{0} = 19 Y_{0} + 24 Z_{0}:
X_{0} = 13021 954967 961017 (17 digits) Y_{0} = 7987 322433 784501 (16 digits) 

Since the equation Ax^{2} + Bxy + Cy^{2} + F = 0 does not change when x is replaced by x and y is replaced by y simultaneously we have another solution:
X_{0} = 13021 954967 961017 (17 digits) Y_{0} = 7987 322433 784501 (16 digits) 

Now we must consider the continued fraction of the other root: (sqrt(313)  725) /608
The expansion is 2 + //1, 3, 1, 1, 17, 2, 1, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2// (the period has 19 coefficients).
The following table shows how the values of Y_{0} and Z_{0} are found (the third column are the values for 304 y^{2} + 725 yz + 432 z^{2}):
c_{n}  y_{n}  z_{n}  

1  0  
2  2  1  198 
1  1  1  11 
3  5  4  12 
1  6  5  6 
1  11  9  1 
17  193  158  6 
2  397  325  11 
1  590  483  2 
8  5117  4189  3 
5  26175  21428  12 
This table is not complete, but there are no solutions in the section not shown.
Y_{0} = 11
Z_{0} = 9
Since X_{0} = 19 Y_{0} + 24 Z_{0}:
X_{0} = 7 Y_{0} = 11 

X_{0} = 7 Y_{0} = 11 

We have to find the continued fraction expansion of the roots of 18 t^{2} + 41 t + 19 = 0, that is, (sqrt(313)  41) / 38
The continued fraction expansion is:
1 + //2, 1, 5, 8, 1, 2, 17, 2, 1, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3//
where the periodic part is marked in bold (the period has 19 coefficients).
The following table shows how the values of U_{0} and V_{0} are found (the third column are the values for 18 u^{2} + 41 uv + 19 v^{2}):
c_{n}  u_{n}  v_{n}  

1  0  
1  1  1  4 
2  1  2  12 
1  2  3  3 
5  11  17  2 
8  90  139  11 
1  101  156  6 
2  292  451  1 
17  5065  7823  6 
2  10422  16097  11 
1  15487  23920  2 
8  134318  207457  3 
5  687077  1 061205  12 
1  821395  1 268662  4 
3  3 151262  4 867191  9 
1  3 972657  6 135853  8 
1  7 123919  11 003044  6 
2  18 220495  28 141941  6 
2  43 564909  67 286926  8 
1  61 785404  95 428867  9 
1  105 350313  162 715793  4 
3  377 836343  583 576246  12 
1  483 186656  746 292039  3 
5  2793 769623  4315 036441  2 
8  22833 343640  35266 583567  11 
1  25627 113263  39581 620008  6 
2  74087 570166  114429 823583  1 
17  1 285115 806085  1 984888 620919  6 
2  2 644319 182336  4 084207 065421  11 
1  3 929434 988421  6 069095 686340  2 
8  34 079799 089704  52 636972 556141  3 
5  174 328430 436941  269 253958 467045  12 
1  208 408229 526645  321 890931 023186  4 
3  799 553119 016876  1234 926751 536603  9 
1  1007 961348 543521  1556 817682 559789  8 
1  1807 514467 560397  2791 744434 096392  6 
2  4622 990283 664315  7140 306550 752573  6 
2  11053 495034 889027  17072 357535 601538  8 
1  15676 485318 553342  24212 664086 354111  9 
1  26729 980353 442369  41285 021621 955649  4 
3  95866 426378 880449  148067 728952 221058  12 
As explained above, x = 2u and y = 2v, so:
X_{0} = 202 Y_{0} = 312 

X_{0} = 202 Y_{0} = 312 

X_{0} = 10130 Y_{0} = 15646 

X_{0} = 10130 Y_{0} = 15646 

X_{0} = 14 247838 (8 digits) Y_{0} = 22 006088 (8 digits) 

X_{0} = 14 247838 (8 digits) Y_{0} = 22 006088 (8 digits) 

X_{0} = 9245 980567 328630 (16 digits) Y_{0} = 14280 613101 505146 (17 digits) 

X_{0} = 9245 980567 328630 (16 digits) Y_{0} = 14280 613101 505146 (17 digits) 

The other root of the equation 18 t^{2} + 41 t + 19 = 0 is (sqrt(313)  41) / 38
Its continued fraction expansion is:
2 + //2, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2, 17, 2, 1, 8, 5, 1, 3, 1//
where the periodic part is marked in bold (the period has 19 coefficients).
The following table shows how the values of U_{0} and V_{0} are found (the third column are the values for 18 u^{2} + 41 uv + 19 v^{2}):
c_{n}  u_{n}  v_{n}  

1  0  
2  2  1  9 
2  3  2  8 
1  5  3  6 
2  13  8  6 
2  31  19  8 
1  44  27  9 
1  75  46  4 
3  269  165  12 
1  344  211  3 
5  1989  1220  2 
8  16256  9971  11 
1  18245  11191  6 
2  52746  32353  1 
17  914927  561192  6 
2  1 882600  1 154737  11 
1  2 797527  1 715929  2 
8  24 262816  14 882169  3 
5  124 111607  76 126774  12 
1  148 374423  91 008943  4 
3  569 234876  349 153603  9 
1  717 609299  440 162546  8 
1  1286 844175  789 316149  6 
2  3291 297649  2018 794844  6 
2  7869 439473  4826 905837  8 
1  11160 737122  6845 700681  9 
1  19030 176595  11672 606518  4 
3  68251 266907  41863 520235  12 
1  87281 443502  53536 126753  3 
5  504658 484417  309544 154000  2 
8  4 124549 318838  2 529889 358753  11 
1  4 629207 803255  2 839433 512753  6 
2  13 382964 925348  8 208756 384259  1 
17  232 139611 534171  142 388292 045156  6 
2  477 662187 993690  292 985340 474571  11 
1  709 801799 527861  435 373632 519727  2 
8  6156 076584 216578  3775 974400 632387  3 
5  31490 184720 610751  19315 245635 681662  12 
1  37646 261304 827329  23091 220036 314049  4 
3  144428 968635 092738  88588 905744 623809  9 
1  182075 229939 920067  111680 125780 937858  8 
As explained above, x = 2u and y = 2v, so:
X_{0} = 10 Y_{0} = 6 

X_{0} = 10 Y_{0} = 6 

X_{0} = 6582 595298 (10 digits) Y_{0} = 4037 589688 (10 digits) 

X_{0} = 6582 595298 (10 digits) Y_{0} = 4037 589688 (10 digits) 

X_{0} = 9 258415 606510 (13 digits) Y_{0} = 5 678867 025506 (13 digits) 

X_{0} = 9 258415 606510 (13 digits) Y_{0} = 5 678867 025506 (13 digits) 

X_{0} = 464 279223 068342 (15 digits) Y_{0} = 284 776584 090312 (15 digits) 

X_{0} = 464 279223 068342 (15 digits) Y_{0} = 284 776584 090312 (15 digits) 

X_{n+1} = P X_{n} + Q Y_{n}
Y_{n+1} = R X_{n} + S Y_{n}
where P, Q, R and S should be determined.
Let M(x, y) = Ax^{2} + Bxy + Cy^{2} = M and N(u, v) = u^{2} + Buv + ACv^{2} = N.
M(p, q) = Ap^{2} + Bpq + Cq^{2}
M(p, q)/A = p^{2} + (B/A)pq + (C/A)q^{2}
M(p, q)/A = p^{2} + Dpq + Eq^{2}
M(p/q, 1)/A = (p/q)^{2} + D(p/q) + E (12)
The roots of M(p/q, 1)/A = (p/q  J) (p/q  J') = 0 (13) are:
J =  B + sqrt(B^{2}  4AC) 2A  and J' =  B  sqrt(B^{2}  4AC) 2A 

It can be easily shown by equating (12) and (13) that:
J^{2} = DJ  E (14)
J'^{2} = DJ'  E (15)
J + J' = D (16)
JJ' = E (17)
The roots of N(p/q, 1) = (p/q  K) (p/q  K') = 0 are:
K =  B + sqrt(B^{2}  4AC) 2  and K' =  B  sqrt(B^{2}  4AC) 2 

so K = AJ, K' = AJ' (18)
M(p, q)/A = (p  Jq)(p  J'q) = M (19)
N(r, s) = (r  Ks)(r  K's) = N (20)
From (18) we obtain:
(p  Jq)(r  Ks) = (p  Jq)(r  AJs) = (pr  AJps  Jqr + AJ^{2}qs)
From (14) we obtain:
[pr  AJps  Jqr + A(DJ  E)qs] = (pr  AEqs)  (Aps + qr + AEqs)J = (pr  Cqs)  (Aps + qr + Bqs)J (21)
From (18) we obtain:
(p  J'q)(r  K's) = (p  J'q)(r  AJ's) = (pr  AJ'ps  J'qr + AJ'^{2}qs)
From (15) we obtain:
[pr  AJ'ps  J'qr + A(DJ'  E)qs] = (pr  AEqs)  (Aps + qr + AEqs)J' = (pr  Cqs)  (Aps + qr + Bqs)J' (22)
Let X = pr  Cqs and Y = Aps + qr + Bqs (23).
Multiplying (21) by (22) we obtain:
(M(p, q)/A) N(r, s) = (X  YJ) (X  YJ') = X^{2}  (J + J')XY + JJ'Y^{2}
Multiplying equations (16) and (17) we obtain: (M(p, q)/A) N(r, s) = X^{2} + DY + EY^{2}
Multiplying by A we get (from (19) and (20)):
AX^{2} + BXY + CY^{2} = MN
Letting M = F and N = 1 we can see that X and Y are also solutions of the original equation.
Let r and s be a solution to N(r, s) = r^{2} + Brs + ACs^{2} = 1,
X_{n} = p, Y_{n} = q, X_{n+1} = X and Y_{n+1} = Y (since the last two pairs of numbers are solutions to the original equation).
From (23) we obtain:
X_{n+1} = rX_{n}  CsY_{n+1}
Y_{n+1} = AsX_{n} + rY_{n+1} + BsY_{n+1}
This means that:
X_{n+1} = P X_{n} + Q Y_{n} Y_{n+1} = R X_{n} + S Y_{n}  P = r (24) Q = Cs (25) R = As (26) S = r + Bs (27) where r^{2} + Brs + ACs^{2} = 1 (28) 

Credits: This method was emailed to me by Iain Davidson. I've made some modifications.
Multiplying the equation by 4A:
4A^{2}x^{2} + 4ABxy + 4ACy^{2} + 4ADx + 4AEy + 4AF = 0
(2Ax + By + D)^{2}  (By + D)^{2} + 4ACy^{2} + 4AEy + 4AF = 0
(2Ax + By + D)^{2} + (4AC  B^{2})y^{2} + (4AE  2BD)y + (4AF  D^{2}) = 0
Let x_{1} = 2Ax + By + D
and g = gcd(4AC  B^{2}, 2AE  BD).
Multiplying by (4AC  B^{2})/g:
4AC  B^{2} g  x_{1}^{2} +  (4AC  B^{2})^{2} g  y^{2} +  2(4AC  B^{2}) (2AE  BD) g  y +  (4AC  B^{2}) (4AF  D^{2}) g  = 0 

4AC  B^{2} g  x_{1}^{2} + g y_{1}^{2} +  (4AC  B^{2}) (4AF  D^{2})  (2AE  BD)^{2} g  = 0 

4AC  B^{2} g  x_{1}^{2} + g y_{1}^{2} +  4A(4ACF  AE^{2}  B^{2}F + BDE  CD^{2}) g  = 0 

where:
y_{1} =  4AC  B^{2} g  y +  2AE  BD g 

X_{n+1} = P X_{n} + Q Y_{n} + K
Y_{n+1} = R X_{n} + S Y_{n} + L
Replacing in the original equation x by Px + Qy + K and y by Rx + Sy + L:
A(Px + Qy + K)^{2} + B(Px + Qy + K) (Rx + Sy + L) + C(Rx + Sy + L)^{2} + D(Px + Qy + K) + E(Rx + Sy + L) + F = 0
(AP^{2} + BPR + CR^{2})x^{2} + (2APQ + B(PS+QR) + 2CRS)xy + (AQ^{2} + BQS + CS^{2})y^{2} + (2ABP + B(KR+LP) + 2CLR + DP + ER)x + (2AKQ + B(KS+LQ) + 2CLS + DQ + ES)y + (AK^{2} + BKL + CL^{2} + DK + EL + F) = 0 (29)
Now we will investigate the values inside the parentheses.
AP^{2} + BPR + CR^{2} = Ar^{2} + BrAs + CA^{2}s^{2} = A(r^{2} + Brs + ACs^{2})
From equation (28) we obtain:
AP^{2} + BPR + CR^{2} = A (30)
2APQ + B(PS+QR) + 2CRS = 2Ar(Cs) + B[r(r+Bs)+(Cs)As] + 2CAs(r+Bs) = 2ACrs + B(r^{2}+BsACs^{2}) + 2ACrs + 2ABCs^{2} = B(r^{2}+Bs+ACs^{2})
From equation (28) we obtain:
2APQ + B(PS+QR) + 2CRS = B (31)
AQ^{2} + BQS + CS^{2} = AC^{2}s^{2} + B(Cs)(r+Bs) + C(r+Bs)^{2} = AC^{2}s^{2}  BCrs  B^{2}Cs^{2} + Cr^{2} + 2BCrs + B^{2}Cs^{2} = AC^{2} + BCrs + Cr^{2} = C(r^{2} + Brs + ACs^{2})
From equation (28) we obtain:
AQ^{2} + BQS + CS^{2} = C (32)
These two equations are equivalent to:
(2AP+BR)K + (BP+2CR)L = D(P1)  ER
and (2AQ+BS)K + (BQ+2CS)L = DQ  E(S1)
Solving the equation system for K and L:
K =  D[BQ  2C(PSQRS)] + E[B(PSRQP)  2CR] 4AC (PS  QR) + B^{2} (QR  PS) 

L =  D[B(PSRQS)  2AQ] + E[BR  2A(PSRQP)] 4AC (PS  QR) + B^{2} (QR  PS) 

Since PS  QR = r(r+Bs)  (Cs)As = r^{2} + Brs + ACs^{2} = 1, these equations can be simplified to:


Now we must show that the expression inside the right parentheses of (29) equals F. This means that we have to prove that the values of K and L just found verify the equation Z = AK^{2} + BKL + CL^{2} + DK + EL = 0 (33).
The expansion is very complicated and will not be reproduced here, but fortunately it is a multiple of 4ACB^{2}, so it cancels the square in the denominator, since it is (4ACB^{2})^{2}.
This means that Z(4ACB^{2}) is an integer number and it is equal to:
AD^{2}Q^{2}  2ADEPQ + AE^{2}P^{2}  AE^{2} + BD^{2}QS  BDEPS  BDEQR + BDE + BE^{2}PR + CD^{2}S^{2}  CD^{2}  2CDERS + CE^{2}R^{2}
Reordering terms:
AD^{2}Q^{2} + BD^{2}QS + CD^{2}S^{2}  CD^{2}  2ADEPQ  BDEPS  BDEQR  2CDERS + BDE + AE^{2}P^{2} + BE^{2}PR + CE^{2}R^{2}  AE^{2}
D^{2}(AQ^{2} + BQS + CS^{2})  CD^{2}  DE(2APQ + BPS + BQR + 2CRS) + BDE + E^{2}(AP^{2} + BPR + CR^{2})  AE^{2}
From (30), (31) and (32):
Z(4ACB^{2}) = CD^{2}  CD^{2}  BDE + BDE + AE^{2}  AE^{2} = 0
This means that Z = 0, so (33) holds, then (29) holds too.
Let K =  K_{D}D + K_{E}E 4AC  B^{2}  and L =  L_{D}D + L_{E}E 4AC  B^{2}  (34) 

To continue simplifying the expressions we should note the following:
K_{D} = BQ  2C(1  S)
K_{D} = B(Cs)  2C(1  r  Bs)
K_{D} = BCs  2C + 2Cr + 2BCs
K_{D} = C(2 + 2r + Bs) (35)
K_{D} = C(P + S  2)
L_{E} = BR  2A(1  P)
L_{E} = ABs  2A + 2Ar
L_{E} = A(2 + 2r + Bs)
L_{E} = A(P + S  2)
K_{E} = B(1  P)  2CR
K_{E} = B(1  r)  2ACs
K_{E} = B  Br  2ACs (36)
L_{D} = B(1  S)  2AQ
L_{D} = B(1  r  Bs) + 2ACs
L_{D} = B  Br  B^{2}s + 2ACs
L_{D}  K_{E} = (4AC  B^{2})s (37)
So:

Generally the numerators will not be multiple of 4AC  B^{2}, so using this formula we cannot find a recurrence for all values of D and E.
For some values of D and E there will be solutions, as shown below. Using equations (24)  (27):
K_{D}L_{E}  K_{E}L_{D} = 4ACr^{2} + 4ABCrs + 4A^{2}C^{2}s^{2}  B^{2}r^{2}  B^{3}rs  AB^{2}Cs^{2}  4ABCs  B^{3}s + 4AC  B^{2}  8ACr + 2B^{2}r =
= (4AC  B^{2}) (r^{2} + Brs + ACs^{2})  (4AC  B^{2})Bs + (4AC  B^{2})  (4AC  B^{2})2r =
= (4AC  B^{2}) (2  2r  Bs)
The equal signs shown below mean congruence mod 4AC  B^{2}.
K_{D}L_{E}  K_{E}L_{D} = 0 => K_{D}/K_{E} = L_{D}/L_{E} (38)
Since K and L must be integers they should be (from (34)):
K_{D}D + K_{E}E = 0 => E = (K_{D}/K_{E})D (39)
L_{D}D + L_{E}E = 0 => E = (L_{D}/L_{E})D
These equations are consistent because of equation (38).
In some cases (see example 6) we can find a recurrence by using the solutions r and s since (r)^{2} + B(r)(s) + AC(s)^{2} = r^{2} + Brs + ACs^{2} = 1.
If no solutions were found (as in example 7), we should use the next pair of solutions (r_{1}, s_{1}) of r^{2} + Brs + ACs^{2} = 1 because there will always be solutions as shown below.
First we should find r_{1} and s_{1} from r and s. To do that we use the formulas (24)  (28).
r_{1} = r r + (ACs)s = r^{2}  ACs^{2}
s_{1} = s r + (r + Bs)s = 2rs + Bs^{2}
Now the values of r and s should be replaced by r_{1} and s_{1}.
From (24): P_{1} = r_{1} = r^{2}  ACs^{2}
From (25): Q_{1} = Cs_{1} = C(2rs + Bs^{2})
From (26): R_{1} = As_{1} = A(2rs + Bs^{2})
From (27): S_{1} = r_{1} + Bs_{1} = r^{2} + 2Brs + (B^{2}  AC)s^{2}
From (35):
K_{1D} = C(2 + 2r_{1} + Bs_{1})
K_{1D} = C[2 + 2(r^{2}  ACs^{2}) + B(2rs + Bs^{2})]
K_{1D} = C[2 + 2(r^{2} + Brs + ACs^{2})  4ACs^{2} + B^{2}s^{2}]
K_{1D} = C[2 + 2 + (B^{2}  4AC)s^{2}]
K_{1D} = (B^{2}  4AC)Cs^{2}
From (36):
K_{1E} = B  Br_{1}  2ACs_{1}
K_{1E} = B  Br^{2} + ABCr^{2}  4ACrs  2ABCr^{2}
K_{1E} = B  B(r^{2} + ACs^{2})  4ACrs
K_{1E} = B  B(r^{2} + ACs^{2} + Brs  Brs)  4ACrs
K_{1E} = B  B(1  Brs)  4ACrs
K_{1E} = (B^{2}  4AC)rs
From (37):
L_{1D}  K_{1E} = (4AC  B^{2})s_{1}
L_{1D} = (B^{2}  4AC) (rs  2rs  Bs^{2})
L_{1D} = (B^{2}  4AC) (rs  Bs^{2})
Therefore:
K_{1} =  K_{1D}D + K_{1E}E 4AC  B^{2}  = CDs^{2}  Ers 

L_{1} =  L_{1D}D + L_{1E}E 4AC  B^{2}  = Ds(r + Bs)  AEs^{2} 

So, finally:
X_{n+1} = (r^{2}  ACs^{2})X_{n}  Cs(2r+Bs)Y_{n}  CDs^{2}  Ers (40) Y_{n+1} = As(2r+Bs)X_{n} + [r^{2} + 2Brs + (B^{2}AC)s^{2}]Y_{n} + Ds(r+Bs)  AEs^{2} (41) 
Notice that in this case, in order to find the solutions using the continued fraction method, we will need to compute two entire periods if the period length is even and four if it is odd.
P = r = 8351
Q = Cs = 32625
R = As = 19575
S = r + Bs = 76474
K =  CD(P+S2) + E(BBr2ACs) 4ACB^{2}  =  340605 109  D +  87174 109  E 

L =  D(BBr2ACs) + AE(P+S2) 4AC  B^{2}  + Ds =  798399 109  D   204363 109  E 

The numerator of K (or L) is not a multiple of the denominator (4AC  B^{2} = 109), so there is no recurrence with the values of P, Q, R, S shown above, except for special cases (according to (39), when E 93 D (mod 109)).
Using the solution r = 8351, s = 6525 we get:
P = r = 8351
Q = Cs = 32625
R = As = 19575
S = r + Bs = 76474
K =  CD(P+S2) + E(BBr2ACs) 4ACB^{2}  =  3125 D  800 E 

L =  D(BBr2ACs) + AE(P+S2) 4AC  B^{2}  + Ds = 7325 D + 1875 E 

So, the recursive relation between solutions is:
X_{n+1} = 8351 X_{n}  32625 Y_{n} + (3125 D  800 E) Y_{n+1} = 19575 X_{n}  76474 Y_{n} + (7325 D + 1875 E) 

Check: Knowing that x = 2, y = 3 is a solution of 3x^{2} + 13xy + 5y^{2}  11x  7y  92 = 0, find other two solutions.
Replacing D = 11 and E = 7 in the previous equations:
X_{n+1} = 8351 X_{n}  32625 Y_{n}  28775
Y_{n+1} = 19575 X_{n} + 76474 Y_{n} + 67450
So, replacing here X_{0} = 2 and Y_{0} = 3, we find X_{1} = 85802 and Y_{1} = 201122.
and replacing X_{1} = 85802 and Y_{1} = 201122, we find X_{2} = 5845 101523 and Y_{2} = 13701 097128.
Replacing these values in the original equation we can check that these values are correct.
P = r = 391
Q = Cs = 1638
R = As = 819
S = r + Bs = 3431
K =  CD(P+S2) + E(BBr2ACs) 4ACB^{2}  =  147 D + 35 E 

L =  D(BBr2ACs) + AE(P+S2) 4AC  B^{2}  + Ds = 308 D   147 2  E 

The numerator of L is not a multiple of the denominator (4AC  B^{2} = 124), so there is no recurrence with the values of P, Q, R, S shown above, except for special cases (when E is even).
Using the solution r = 391, s = 273 we get:
P = r = 391
Q = Cs = 1638
R = As = 819
S = r + Bs = 3431
K =  CD(P+S2) + E(BBr2ACs) 4ACB^{2}  =  4563 31  D   1092 31  E 

L =  D(BBr2ACs) + AE(P+S2) 4AC  B^{2}  + Ds =  4563 62  D   9555 31  E 

The numerator of K (or L) is not a multiple of the denominator (4AC  B^{2} = 124), so there is no recurrence with the values of P, Q, R, S shown above, except for special cases.
Using (40) and (41):
P_{1} = r^{2}  ACs^{2} = 1188641
Q_{1} = Cs(2r+Bs) = 4979520
R_{1} = As(2r+Bs) = 2489760
S_{1} = r^{2} + 2Brs + (B^{2}AC)s^{2} = 10430239
K_{1} = CDs^{2}  Ers = 106743 D + 447174 E
L_{1} = Ds(r+Bs)  AEs^{2} = 936663 D  223587 E
So, the recursive relation between solutions is:
X_{n+1} = 1188641 X_{n}  4979520 Y_{n} + (106743 D  447174 E) Y_{n+1} = 2489760 X_{n} + 10430239 Y_{n} + (936663 D  223587 E) 

Check: Knowing that x = 4, y = 7 is a solution of 3x^{2} + 14xy + 6y^{2}  17x  23y  505 = 0, find other two solutions.
Replacing D = 17 and E = 23 in the previous equations:
X_{n+1} = 1188641 X_{n}  4979520 Y_{n} + 5146869
Y_{n+1} = 2489760 X_{n} + 10430239 Y_{n}  10780770
So, replacing here X_{0} = 4 and Y_{0} = 7, we find X_{1} = 34 464335 and Y_{1} = 72 189943.
and replacing X_{1} = 34 464335 and Y_{1} = 72 189943, we find X_{2} = 318 505538 201756 and Y_{2} = 667 150425 396007.
Replacing these values in the original equation we can check that these values are correct.
Press here to run a JavaScript program compatible with Internet Explorer 3.0 or Netscape Navigator 3.0 or later (not finished yet).
Press here to run a Java program compatible with Internet Explorer 4.0 or Netscape Navigator 4.0 or later. It is about 50 times faster than the previous program.
Both programs run in two modes: solution only and stepbystep. In the last mode, the program explains how the results are found.
Click here to send me your comments.