Status of the Diophantine module

Hi All,

In my project proposal for a Diophantine equation for SymPy, I mentioned the following five deliverables.

1. Linear Diophantine equation a_1x_1 + a_2x_2 + . . . + a_nx_n = b:
I implemented solutions for linear diophantine equations, you can access this functionality through `diop_linear()`.

2. Simplified Pell equation, x^2 - Dy^2 = 1:
Not only I implemented solutions for simplified Pell equation, I completely solved the general binary quadratic equation ax^2 + bxy + cy^2 + dx + ey + f = 0.

3. The equation, x^2 + axy + y^2 = z^2:
I implemented solutions for more general ternary quadratic equation ax^2 + by^2 + cz^2 + dxy + eyz + fxz = 0.

4. Extended Pythagorean equation,  x_1^2 + x_2^2 + . . . + x_n^2 = x_{n+1}^2:
I implemented solutions for slightly more general equation a_1^2x_1^2 + a_2^2x_2^2 + . . . + a_n^2x_n^2 = a_{n+1}^2x_{n+1}^2.

5. General sum of squares, x_1^2 + x_2^2 + . . . + x_k^2 = n:
This is a computationally hard problem and method I implemented finds only one solution. It’s quick and work for large n but not complete. I also implemented a brute force version which finds all the solutions but it doesn’t work for larger n.


Change of plans for general sum of squares

Hi All,

In my last post I described about Euler’s four square identity to for solving general sum of squares equation. My idea was to factorize the number, represent each prime as a sum of four squares, then use Euler’s four square identity to construct the sum of four squares representation for the original Number. But I found it slower than the algorithm I found here. So I adapted the latter. The idea is to reduce the problem of representing a number as a sum of four squares to representing a number as a sum of three squares.

Algorithm for representing a positive number n as a sum of three squares

1. If n == 0, then return (0, 0, 0)
2. Write n as 4^vn_{1}
3. if n_{1} \in S then return the hard coded representation of n_{1}. Here, S = {2, 3, 10, 34, 58, 85, 130, 214, 226, 370, 526, 706, 730, 1414, 1906, 2986, 9634}. Representations for these numbers can be found here.
4. If n_{1} is a perfect square then return (2^v\sqrt{n_{1}}, 0, 0)
5. if n_{1} = 3 (mod 8), find an odd number i, i < \sqrt{n_{1}} such that \frac{n_{1} - i^2}{2} is a prime. Set x = i and p = \frac{n_{1} - i^2}{2}. Find two numbers y, z such that y^2 + z^2 = p. (You can use Cornacchia’s  algorithm for this. Return (2^vx, 2^v(y + z), 2^v|y - z|).
6. If n_{1} = 2 (mod 8) or n_{1} = 6 (mod 8) then find an odd number i, i < \sqrt{n_{1}} such that n_{1} - i^2 is a prime. Else find an even number i, i < \sqrt{n_{1}} with the above requirement. Set x = i and p = n_{1} - i^2. Find y, z such that y^2 + z^2 = p . Return (2^vx, 2^vy, 2^vz).

Note that above algorithm can not be used if n_{1} = 8k + 7 for some integer k \in Z. That is if n is in the form 4^v(8k + 7).

Algorithm for representing a positive number n as a sum of four squares

Every non-negative integer n can be represented as a sum of four squares.

1. Write n as 4^vn_{1}.
2. If n_{1} = 7 (mod 8) then set d = 2 and n_{1} = n_{1} - 4. If n_{1} = 6 (mod 8) or n_{1} = 2 (mod 8) then set d = 1 and n_{1} = n_{1} - 1. Otherwise set d = 0.
3. Use the algorithm described earlier to represent n_{1} as a sum of three squares. Say, x^2 + y^2 + z^2 = n.
4. Return (2^vd, 2^vx, 2^vy, 2^vz).

By using these two algorithms one can represent any non-negative integer as a sum of four squares. In the general sum of squares equation, where a given integer n needs to be represented as a sum of k squares, we can divide k  variables into segments of four and divide n into that same number of segments and represent each segment by four squares using the algorithm given above.

Status of the Diophantine Module

Hi All,

I am really pleased to say that my work (proposed) with Diophantine equation is coming to an end. I have implemented all the deliverables in my project report. When the current PR gets merged into the master, I can start pushing the rest of my work. I also hope to improve the documentation during the weekend and get them in quickly. Ondrej is currently reviewing the PR. I invite all of you to go through the work and give your feedback. I am really thankful to Aaron, Pernici, Stephen and Julien for the support given so far. First of all let me give you a rough idea about the work currently not in Github.

General Pythagorean equation

A general Pythagorean equation is an equation of the form x_{1}^2 + x_{2}^2 + . . . + x_{k}^2 = x_{k+1}^2. The solutions for the equation can be given by using k parameters like below,

x_{1} = m_{1}^2 + m_{2}^2 + . . . + m_{k-1}^2 - m_{k}^2
x_{2} = 2m_{1}m_{k}
x_{k} = 2m_{k-1}m_{k}x_{k +1} = m_{1}^2 + m_{2}^2 + . . . + m_{k-1}^2 + m_{k}^2

I implemented solutions for slightly more  general equation a_{1}^2x_{1}^2 + a_{2}^2x_{2}^2 + . . . + a_{k}^2x_{k}^2 = a_{k+1}^2x_{k+1}^2. Solutions for this equation can be constructed from the solutions of the former equation, multiplying each solution x_{i} of it by \frac{lcm(a_{1}, a_{2}, . . . a_{n})}{a_{i}}.

General sum of n squares

I also implemented solutions for x_{i}^2 + x_{2}^2 + . . . + x_{n}^2 = k where k is an integer. It is obvious that equation is solvable only when k \geq 0.

Lagrange’s four square theorem states\textsuperscript{[1]} that every non-negative number can be expressed as a sum of four squares. It is also known that every integer not in the the form 4^k(8m + 7) can be expressed as a sum of three squares\textsuperscript{[2]} where m, n are non-negative integers.  Also any integer which doesn’t contain, in it’s canonical representation, odd powers of a prime of the form 4m + 3 can be expressed as a sum of two squares\textsuperscript{[3]}. For example, N = 2^2.5^3.11^2.3^3 can’t be expressed as a sum of two squares since it contains and odd power of 3, which is a prime of the form 4m+3. But we can express N = 2^2.5^3. 11^2.3^4 as a sum of two squares.

Also, there is an interesting identity, found by Euler, known as Euler’s four square identity which can be stated as,

(a_1^2+a_2^2+a_3^2+a_4^2)(b_1^2+b_2^2+b_3^2+b_4^2) = (a_1 b_1 - a_2 b_2 - a_3 b_3 - a_4 b_4)^2 + (a_1 b_2 + a_2 b_1 + a_3 b_4 - a_4 b_3)^2 + (a_1 b_3 - a_2 b_4 + a_3 b_1 + a_4 b_2)^2 + (a_1 b_4 + a_2 b_3 - a_3 b_2 + a_4 b_1)^2.

So, if we can represent each prime divisor of a number N as a sum of four squares, we can then use this identity to construct such a representation forn. This is the idea behind my implementation of representing number as a sum of four squares. When we know how to do this, we have several approaches for representing a given non-negative integer k as a sum of n squares. Most obvious thing to do is, represent k as a sum of four squares and set other n -4 variables to zero. We can also partition k into approximately n/4 groups and represent each as a sum of four squares and later combine those results. I still haven’t decided on how to do this, I would like to know the ideas of the community.  I used a slightly modified version of  the algorithm found in\textsuperscript{[5]}.  I’ll describe the algorithm in my next blog post.

Apart from that, I did a lot of bug fixing and reviewed Pernici’s Pull request. My tests are now 2x faster after using his functions. A huge thank should go to Pernici for doing a great Job with solving quadratic congruences.


[1] Lagrange’s four square theorem, [online], Available:
[2] Integer as a sum of three squares, [online], Available:
[3] Sum of Squares, [online], Available:
[4] Euler’s four square identity, [online], Available:

Improving the solving process using Cornacchia

Hi All, It has been sometime since I last wrote about my project.  I spent most of my time improving the solving methods used in binary quadratic forms. I also did a fair bit of refactoring of the earlier work and bug fixing. The implementation of ternary quadratic forms is now almost done and once the pull request by Pernici gets merged, we can use his function for finding the solutions for the quadratic congruence x^2 \equiv a \ (mod \ p) and the time taken to find the solutions for ternary / binary quadratic forms will be reduced considerably.

Tweaks to the algorithms used in solving quadratic binary forms

I had a misunderstanding earlier that the binary quadratic form ax^2 + bxy + cy^2 + dx + ey + f = 0 can be converted to the form x^2 - Dy^2 = N only in the case \Delta > 0 and \Delta is not a perfect square. Here \Delta = b^2 - 4ac. It occurred to me only when I read [1] more carefully that any binary quadratic with a \neq 0 and c \neq 0 can be transformed to the latter form. With this new discovery, our treatment of the binary quadratic becomes more general and efficient.  Now in both the cases \Delta > 0, \Delta is not a perfect square and \Delta < 0, we transform the equation to the form x^2 -Dy^2 = N. The former case was earlier addressed under solving generalized Pell equation as the equation reduces to a generalized Pell equation in that case. My earlier approach with the latter case was to use brute force. I managed to implement the Cornacchia’s algorithm for this case as suggested by Pernici. Cornacchia’s algorithm can be used to solve the equation ax^2 + by^2 = m with gcd(a, b) = 1 = gcd(a, m) and a > 0, b > 0. In our form under the case \Delta < 0, a = 1 and b = -D > 0, so these two conditions are automatically fulfilled. Below is a rough outline of the Cornacchia’s algorithm. Refer [2] and [3] for more details.

Solving ax^2 + by^2 = m by Cornacchia’s method

1. Let p be such that ap + mq = gcd(a, m).
2. Construct a set v such that v = \{ x | x^2 \equiv -bp \ (mod \ m) \}. We can ommit solutions greater than m/2.
3. For each t \in v set u_{1} = t, r_{1} = m.
4. Define the sequences, u_{i+1} = r_{i} and r_{i+1} = u_{i} \ mod \ r_{i} for i > 1.
5. Find the terms u_{i}, r_{i} until the condition ar_{i}^2 < m is met.
6. Set m_{1} = m - ar_{i}^2.
7. If b | m_{1} and m_{2} = m_{1}/b is a perfect square we have the solution x = r_{i} and y = \sqrt{m_{2}}.
8. Add (x, y) to a set S.
9. When we are done with all t, S contains all the solutions to the equation ax^2 + by^2 = m.

Cornacchia’s method finds solutions such that gcd(x, y) = 1. So it won’t find solutions to the equation x^2 + y^2 = 20 since the solutions of this equation is (x, y) = (4, 2) and gcd(x, y) > 1. So we have to extract out square factors in m and find solutions to the equation with this new m value and reconstruct the solutions. For example in the above example, 2 is a square factor of 20 (i.e square of two divides 20) so, we find solutions to m = 20 /2^2 = 5, i.e for the equation x^2 + y^2 = 5 and then reconstruct the solutions for the original equation.

diophantine() and check_solutions()

I replaced the main routine of the module by a new method `diophantine()`. What the old method,`diop_solve()` did was, when given with an input equation, it tried to find  the type of the equation using `classify_diop()` and called the appropriate solving function accordingly. So when given an equation like x^2 + 3xy + 5x = 0 , diop_solve will identify this as a binary quadratic and try to solve this using methods for binary quadratics. But the best thing to do in this case is to factorize the above equation as x(x + 3y + 5) = 0 and solve x = 0 and x+3y+5 = 0 separately. After that, we can construct the general solution by adding parameters. The general solution (x, y, z) = (0, n_{1}, n_{2}) will be constructed after solving x = 0 and will be added to the solution set for x^2 + 3xy + 5x = 0 . This was done by introducing the new routine `diophantine()`. This is now the main routine for the module and when we call `diophantine(eq)`, first `diophantine()` tries to factor `eq` and then solve each factor using `diop_solve()`. Then it creates the general solution and that will be added to the solution set for `eq`. This is a very good approach and it is also a very common method used when we try to solve Diophantine equations manually.

I also unified various kinds of  functions in `` which were used to find whether the solutions returned by a function satisfies the input equation. I wrote a new method `check_solutions()`. Aaron should be thanked for proposing such a method and now the `` looks nicer and whoever trying to hack it doesn’t need to study various kinds of solution checking functions. To be honest, I was running out of names for these functions because a new one had to be created for each equation type.

Bug Fixing

I spent a considerable time fixing the bugs in the module, especially towards the end of the last week. I had done some mistakes when transforming general ternary quadratic equation into the form ax^2 + by^2 + cz^2 = 0. Also, in `diophantine()` method, I hadn’t check for the case when None is returned to indicate the insolubility of the equation. I also found a bug in the `factor_list()` function under Python3. Here is the issue I created for it:

Hmm, That is pretty much everything I did during the past weeks. I hope in the upcoming week I can start on solving extended Pythagorean equation.


[1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003.
[2] Solving ax^2 + by^2 = m by Cornacchia’s method, Number theory web,
[3] A. Nitaj, L’algorithme de Cornacchia, Expositiones Mathematicae 13 (1995), 358-365

Holzer reduction and some modifications

Hi All, This week I managed to implement Holzer’s reduction and change the behaviour of DE module a little bit.  First let’s take a look at Holzer’s reduction.

Holzer Reduction

Holzer reduction is concerned with reducing solutions of the quadratic ternary nonrmal equation ax^2 + by^2 +cz^2 = 0. The Holzer’s theorem says that if the above equation is solvable there is always a solution (x, y, z) such that x \leq \sqrt{|bc|}, y \leq \sqrt{|ac|}, z \leq \sqrt{|ab|}. The algorithm for this is explained in [1]. Below is a rough sketch. Before applying the algorithm we have to make abc square free and we assume that a, b positive and c is negative(We can always select a, b, c in such a way that this is the case). Suppose there is a solution (x_{0}, y_{0}, z_{0}) that is not Holzer reduced.

  • If c is even set k = c /2 and let u_{0}, v_{0} be any solution to the equation k = uy_{0} - vx_{0}.
  • If c is odd, c is odd let k = 2c and  let u_{0}, v_{0} be any solution to the equation c = uy_{0} - vx_{0}.
  • Let w be the nearest integer to -(aux_{0} + bvy_{0})/(cz_{0}).
  • Then the following expressions for x, y, z gives a solution such that z < z_{0}.
    x = (x_{0}(au_{0}^2 + bv_{0}^2 + cw^2) - 2u_{0}(au_{0}x_{0} + bv_{0}y_{0} + cwz_{0})) / k
    y = (y_{0}(au_{0}^2 + bv_{0}^2 + cw^2) - 2v_{0}(au_{0}x_{0} + bv_{0}y_{0} + cwz_{0})) / k
    z = (z_{0}(au_{0}^2 + bv_{0}^2 + cw^2) - 2w(au_{0}x_{0} + bv_{0}y_{0} + cwz_{0}))) / k

Continuing this manner one would find a solution such that z_{0} \leq \sqrt{|ab|}.

Some changes to DE module

I made some changes to the structure of DE module too. Now every type of equation returns a tuple ordered according to the sorted order of variable names used in the equation. Quadratic ternary forms and linear equations returns only one tuple, the parametric solution,  and quadratic binary equation returns a set of tuples. This format may be changed in future. The most important change I did to the DE module is that, now before solving any given equation, DE module checks whether that can be factored and If so, solves those factors separately. For example, given the equation y^2 - 7xy + 4yz = 0, it will solve y = 0 and y-7x+4z = 0 separately and combine those results. Also I managed to implement a general method for testing that would check whether the solutions returned by various methods satisfy the original equation. So I hope to get rid of the other redundant methods like check_ternary_quadratic(), solution_ok_quadratic() in the test file. A thank should go to Aaron for proposing such a methodology.

Future work

Pernici has pointed out that some of the solution methodologies in quadratic binary form can be improved so that they take less time. He also has done a great job in implementing the solutions for quadratic congruence x^2 \equiv a (mod \ b). I hope to use those results when his PR gets merged into master. He proposed a method due to cornacchia to improve the speed of the solutions but after a little surveying, I found that implementing solutions for ax^2 + bxy + cy^2 = N can be used to improve the speed. Current algorithms for solving binary quadratic equation ax^2 + bxy + cy^2 + dx + ey + f = 0 employs a general solving method so it miss out some speed enhancements for the special case ax^2 + bxy + cy^2 = N . My plan is if ax^2 + bxy + cy^2 + dx + ey + f = 0 can be converted to the form ax^2 + bxy + cy^2 = N, then solve it using methods for that, otherwise apply the general method. I found a good reference [2] for this. A huge thank should go to Pernici for pointing out the improvements.


[1] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
[2] Binary Quadratic Forms: An Algorithmic Approach, J. Buchmann and U. Vollmer, Springer-Verlag, Berlin, Heidelberg, 2007.

Descent method with gaussian reduction.

Sorry for a very late blog post, I had a really busy week and finally I managed to find sometime to write about the progress I made in the last week. I implemented an improved version of the descent method which uses gaussian reduction. I found this algorithm in [1]. Here is a sketch of the algorithm.

descent(a, b)

1. If |a| > |b| then swap a and b, solve the resulting equation, then swap y and z in the solution obtained.
2. If b = 1 then set (x, y, z) = (1, 1, 0) and stop.
3. If a = 1  then set (x, y, z) = (1, 0, 1) and stop.
4. If b = -1 there is no solution (since a must be -1).
5. If b = -a then set (x, y, z) = (0, 1, 1) and stop.
6. If b = a then let (x_{1}, y_{1}, z_{1}) be a solution of X_{1}^2 + Z_{1}^2 = aY_{1}^2 , set (x, y, z) = (ay_{1}, x_{1}, z_{1}) and stop.
7. Let w be a solution to x^2 \equiv a (mod \ b). with w \leq |b|/2, and set (x_{0}, z_{0}) = (w, 1), so that x_{0}^2 - az_{0}^2 \equiv 0 (mod \ b).
8. Use lattice reduction to find a new nontrivial solution (x_{0}, z_{0}) to the congruence X^2 - aZ^2 \equiv 0 (mod \ b) with x_{0}^2 + |a|z_{0}^2 as small as possible.
9. Set t = (x_{0}^2 - az_{0}^2)/b and write t = t_{1}t_{2}^2 with t_{1} square-free.
10. Let (x_{1}, y_{1}, z_{1}) be a solution to X^2 - aZ^2 = t_{1}Y^2 then
(x, y, z) = (x_{0 }x_{1} + az_{0}z_{1}, t_{1}t_{2}y_{1}, z_{0}x_{1} + x_{0}z_{1})

I spent a lot of time finding an algorithm for lattice reduction. There were generalized algorithms but I looked for something which is faster because we are concerned with vectors with two bases. I found  gaussian reduction algorithm in [2] which uses a method specific to this case.


[1] .. Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
[2] .. Gaussian lattice Reduction [online]. Available:

Solving Ternary Quadratic forms by descent method

Now SymPy can solve Ternary Quadratic forms !!!! I thought of adopting a new style for giving more details on this amazing news, So Let’s try a Q&A approach and see how it goes.

Q1. What is a “Ternary Quadratic form”?

A ternary quadratic form is an equation of the form, Q(x, y, z) \equiv Ax^2 + Bxy + Cxz + Dy^2 + Eyz + Fz^2 = 0 (Yes, It’s a little bit odd representation with square terms scattered but you will see later why this is a better representation). One important thing worth observing is that the degree of each term is precisely two. This makes it a lot easier to solve the equation and moreover, it makes finding integer solutions to the equation equivalent to finding the rational solutions. We are particularly interested in the case where A, B, C, D, E, F, x, y, z \in \mathbb{Z} and  (x, y, z) \neq (0, 0, 0). In other words, we want to find non-trivial integer solutions to the equation when it has integer coefficients (clearly (x, y, z) = (0, 0, 0) is always a solution to this equation so we call it the `trivial` solution).

Q2. Can we make the equation more simpler to deal with?

Yes, It turns out that we can. We can convert Q to a normal form Q' \equiv ax^2 + by^2 + cz^2 = 0 by two rational transformations of the variables. Lets assume that A \neq 0(We can always switch variables if A = 0).  With an easy simplification followed by clearing of the denominators , one can easily see that the transformation given by x \rightarrow x - (By + Cz)/2A puts our original equation in the form A'x ^2 + D'y^ 2 + E'yz + F'z^ 2= 0 where A', D', E' , F' \in \mathbb{Z}. Finally,  the transformation y \rightarrow y - E'z/2D' puts the equation in the desired form, ax^2 + by^2 + cz^2 = 0. Since finding integer solutions is equivalent to finding rational solutions and the transformations we used were rational, If an integer solution exist for Q we can conclude that the new form Q' should also have an integer solution.

Q3. What are the required conditions to have a non-trivial solution(s)?

It’s sufficient to consider the conditions on  Q' to have a non-trivial integer solution. If Q' has a non-trivial integer solution, we can use the inverse transformation to recover the solutions for Q. It turns out that the conditions for solubility of Q' are,
1. All three variables, a, b, c should not have the same sign.
2. The three quadratic congruence equations, x^2 + bc \equiv 0 \ (mod \ a), x^2 + ca \equiv 0 \ (mod \ b), x^2 + ab \equiv 0 \ (mod \ c) should be solvable.

You can find the proofs for  the sufficiency of these conditions in the references. It is also a proven fact that If one non-trivial solution exist for one of the equations, there will be infinitely many non-trivial solutions.

Q4. What are the algorithms that can be used ?

There are quite a few algorithms which can be used to tackle down a non-trivial solution. There is a method which uses sieving and another one based on quadratic number fields. But the descent method due to Lagrange beats the above methods when it comes to performance. The idea is like this. Suppose there is a solution, then under some conditions we can deduce that an even smaller solution should exist. We continue making the solution smaller and smaller until we find an equation which has a solution we can easily spot.  First we convert ax^2 + by^2 + cz^2 = 0 to the form w^2 = Ax^2 + By^2 (We can easily do this by multiplying entire equation by -c). If A = 1 then (x, y, w) = (1, 0, 1) is a solution to our equation. If B = 1 then (x, y, z) = (0, 1, 1) is a solution to our equation. Otherwise we continue reducing A and B until we arrive at a situation like above.

Q5. How can one find all the non-trivial solutions?

Well, All the non-trivial solutions can be represented using two parameters, p, q which are co-prime to each other. Suppose we have found a non-trivial solution (x_{0}, y_{0}, z_{0}) for the equation Q(x, y, z) = 0. We can always assume that x_{0} \neq 0 by a switching of the variables. Then assuming (x, y, z) = (x_{0}r, y_{0}r + p, z_{0}r + q) and substituting in the original equation, one can find an expression for r using p and q. So we can find parametrized solutions for (x, y, z) by substituting the value of r.

Q6. Any Improvements to current methodology?

Yes, An efficient algorithm for solving the quadratic congruence should be implemented. Also, there is an algorithm in MAGMA which currently solves this kind of equations. I plan to implement that also so that I can verify my results and if it is faster than the descent method, replace the descent method with it.

And a special thank should go to Mario for helping me out with the bugs in my code.


[1] The algorithmic resolution of Diophantine equations, Nigel P. Smart, London Mathematical Society Student Texts 41, Cambridge University Press, Cambridge, 1998.
[2] Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
[3] L.J. Mordell, Diophantine Equations, Academic Press, 1969.

From Euclid To Gauss via Pell

While the fourth week of GSoC is coming to an end, I managed to finish two of the five deliverables of my project proposal. So the fourth week is more of a transition week for me as I had to correct few bugs associated with current implementation of QDEs (Quadratic Diophantine Equation) and LDEs (Linear Diophantine Equation) and find new resources for the future work, solving Ternary Quadratic Forms.

In QDE I had a little problem dealing with the case where B**2 - 4*A*C is a perfect square. General solution procedure for this case is to use a transformation which converts this case to the Simple Hyperbolic case which I had solved previously. Then the original solutions can be recovered from the solutions of the transformed equation using two divisible criteria which should be satisfied to honour the transformation. Things went wrong when the solutions of the transformed equation involved parameters because I had to check finitely many equivalence classes with respect to modulo 4*A*sqrt(B**2 - 4*A*C in that case and find whether there are any classes satisfying the criteria and find a nice representation for the solutions if at least one such class exist. I found this hard going but finally managed to implement it. Work related to the QDEs are almost finished and there are a few optimization steps to be carried out. For example I use brute force to solve the congruence equation x**2 = D (mod a) and there is an algorithm which is more efficient. I hope to implement it since this equation will come up in my future work too. Several extensions can also be carried out like finding the rational solutions satisfying a given QDE and implementing other algorithms like the cyclic method which also solves QDEs.

I also corrected an error I had done in linear Diophantine equations due to which the solutions returned by the solver were only a subset of the exact solutions. I corrected this and checked the results with the Wolfram Alpha. Now, the solutions returned by the solver can be made identical to that returned by the Wolfram Alpha by a single shift or by inverting the parameter. Since the parameters returned in solutions for linear Diophantine equations has no boundary conditions and can be an any integer, this does not make the two solution sets different.

During the week, I also studied the ternary quadratic forms, on which I plan to work in the weeks to come. A ternary quadratic form is a homogeneous equation in three variables and having a degree two, i.e. a equation of the form Ax**2 + By**2 + Cz**2 + Dyz + Exz + Fxy = 0. I found a good resource on this, “Algorithmic resolution of Diophantine equations” by Nigel P. Smart.

Solving the generalized Pell equation

Continuing the work from last week, I was aiming to solve quadratic Diophantine equations with delta = B**2 – 4*A*C > 0. I found two better references on this case which can be found in [1] and [2]. These two papers describe all the necessary algorithms for solving the generalized Pell equation, to which the quadratic Diophantine equation reduces in the case delta > 0. So let’s take a look at the generalized Pell equation and methods I used to solve the equation.

The generalized Pell equation is an equation of the form x**2 – D * y**2 = N. The main case of the equation is D > 0, D is not a perfect square and abs(N) > 0. Other cases are fairly easy and solutions are straightforward.

Case D < 0

Since D < 0,  equation becomes x**2 + |D| * y**2 = N, which is always positive or zero. So there are no solutions when  N < 0 and only the trivial solution x = 0 and y = 0 when N = 0. If N > 0, we can note that |D| * y**2 <= N which implies y <= sqrt(N/|D|), so a brute force search on this interval would find us all the solutions. That’s what I did in this case. But there are more better ways like using methods of binary quadratic forms for positive definite forms. I hope to implement it in the near future.

Case D = 0

Equation reduces to x**2 = N, solution procedure is straight forward and no solutions exist if N < 0 or N is not a perfect square.

Case D > 0

If D is a square, then let D = r**2 for some integer r. Then the LHS is x**2 – r**2 * y**2 = N –> (x – r*y) (x + r*y) = N. Now if N = 0 either x – r*y = 0 or x+ r*y = 0 and the parametric solutions for the above two cases are (r*t, t) and (-r*t, t) respectively. If N does not equal zero, then both x – r*y and x + r*y should be factors of N. We should find an upper bound M for y so that we can search through the solutions in the range 0 <= y <= M (We don’t have to consider negative solutions, reason is explained later in the article).  Let x – r*y = a and x + r*y = b. Solving for y we get y = (b – a)/(2*r). Now N = (x – r*y) * (x + r*y) = a * b so b <= N and a >= 1 since a is a positive factor o f N. which implies b <= N and -a <= -1. Now adding up the two inequalities gives b – a <= N – 1 which gives y <= (N – 1)/ (2*r) which is the required upper bound.

Now suppose D is not a perfect square. If  N = 0, x**2 – D*y**2 = 0 –> x**2 = D * y**2.  Suppose y is non zero. Then y should divide x (Otherwise y**2 can’t divide x**2 which is the case here). Putting x = k*y on the LHS one can see that D = k**2 which is in contradiction to the fact that D is not a perfect square. So y = 0 and consequently x = 0. Hence if N = 0 we will have only the trivial solution x = 0 and y = 0. Now if N is non zero comes our main case.

PQa Algorithm

This algorithm computes the continued fraction expansion of the quadratic irrational (P + sqrt(D))/Q and some auxiliary variables, i.e of numbers like (2 + sqrt(3))/5 and (-34 + sqrt(19))/24. Note that D should not be perfect square otherwise the expression (P + sqrt(D))/Q will not be a quadratic irrational. For more details on algorithm refer [2]. Algorithm takes three variables P, Q, D as input and output six quantities P_i, Q_i, a_i, A_i, B_i and G_i for each time it’s called for a given set of initial values (This is implemented as a generator, please see python documentation). a_i is the ith term of continued fraction expansion of (P + sqrt(D))/Q and others are required in the LMM algorithm described below. Here i >= 0.

It is also important to know that continued fraction representation of a quadratic irrationality is periodic (most of the time) after a few terms. In the solution we need to know the sum of the lengths of the aperiodic part and the periodic part. I implemented a separate function, length() to find this.

LMM Algorithm

LMM algorithm finds primitive solutions of the each equivalence class of the solution set of x**2 – D * y**2 = N. We have to construct the general solution of each class afterwards. Here is a rough description how the algorithm works.

1. Make a list of f > 0 such that f**2 divides N
2. Set m = N / f**2 and find all z such that z**2 ≡ D( mod |m|) and |m|/2 < z <= |m|/2
3. Apply PQa algorithm with P = z, Q = |m| and D = D
4. Continue till there is an i >= 1 such that |Q_i| = 1 or reach the end of first period of continued fraction expansion without such a i, i.e continue till i exceeds the value returned by length()
5. If you reached an i with |Q_i| = 1 then set r = G_i-1 and s = B_i-1. if r**2 – D * s**2 = m, then add x = f*r and y = f*s to the list of solutions. Otherwise r**2 – D * s**2 = -m. Now if x**2 – D * y**2 = -1 has solutions, find it’s minimal positive solution (u, v) and add x = f*(r*u+ s*v*d) and y = f*(r*v + s*u) to the list of solutions. If x**2 – D * y**2 = -1 does not have solution, continue with next z
6.Repeat this for every z for corresponding f, and then for every f corresponding to N

We can solve the x**2 – D * y**2 = -1 by reapplying LMM algorithm but in the Diophantine module I separately solved the equation for |N| = 1 and used the solutions returned by them in the LMM algorithm.

If you are interested, a lot of cool information about solving process and the history of the Pell equation can be found in [3].


[1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003.
[2] Solving the generalized Pell equation x**2 – D*y**2 = N, John P. Robertson, July 31, 2004.
[3] H.W Lenstra, Jr.

Quadratic Diophantine equation – I

Quadratic Diophantine equation is an equation of the form Ax**2 + Bxy + Cy**2 + Dx + Ey + F = 0 where A, B, C, D, E, and F are integer constants and x and y being integer variables. Study of this equation has always been an interesting area among number theorists. The famous pell equation is a special case of the above with delta = B**2-4AC > 0 and delta not being a perfect square. Normally, this equation is broken down into five cases for analytical purposes.

1) A = B = C = 0 (Linear case): Reduces to a linear Diophantine equation of two variables.

2) A = C = 0 and B != 0  (Simple hyperbolic case): Equation reduces to (Bx + E) (By + D) = DE – BF, which can be solved by considering the factors of DE – BF.

3) B**2 – 4AC < 0 (Elliptical case): In this case, values of the x should lie between the roots of the equation (B**2-4AC)x**2 + 2(BE – 2CD)x + E**2 – 4CF = 0. Values for x should be selected so that y is an integer.

4) B**2 – 4AC = 0 (Parabolic case): Solution procedure is rather complex in this case.

I will describe these cases in detail in the future posts. I had almost completed the above cases at the start of Week 2.

5) delta = B**2 – 4AC > 0: This is split into several subcases.

Case delta = B**2 – 4AC > 0:

Subcase 1: D = E = 0:

This is the homogeneous case and again considered under two cases F = 0 and F != 0. If F = 0, then x = 0 and y = 0 are solutions. More solutions may exist if B**2 – 4AC is a perfect square. Otherwise x = 0 and y = 0 is the only solution. I implemented this case in the module. If F != 0 the solution procedure is rather complex and involves continued fractions. I am currently working on this.