In my project proposal for a Diophantine equation for SymPy, I mentioned the following five deliverables.
1. Linear Diophantine equation :
I implemented solutions for linear diophantine equations, you can access this functionality through `diop_linear()`.
2. Simplified Pell equation, :
Not only I implemented solutions for simplified Pell equation, I completely solved the general binary quadratic equation .
3. The equation, :
I implemented solutions for more general ternary quadratic equation .
4. Extended Pythagorean equation, :
I implemented solutions for slightly more general equation .
5. General sum of squares, :
This is a computationally hard problem and method I implemented finds only one solution. It’s quick and work for large but not complete. I also implemented a brute force version which finds all the solutions but it doesn’t work for larger .
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 , then return
2. Write as
3. if then return the hard coded representation of . Here, . Representations for these numbers can be found here.
4. If is a perfect square then return
5. if , find an odd number such that is a prime. Set and . Find two numbers such that . (You can use Cornacchia’s algorithm for this. Return .
6. If or then find an odd number such that is a prime. Else find an even number with the above requirement. Set and . Find such that . Return .
Note that above algorithm can not be used if for some integer . That is if is in the form .
Algorithm for representing a positive number n as a sum of four squares
Every non-negative integer can be represented as a sum of four squares.
1. Write as .
2. If then set and . If or then set and . Otherwise set .
3. Use the algorithm described earlier to represent as a sum of three squares. Say, .
4. Return .
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 needs to be represented as a sum of squares, we can divide variables into segments of four and divide into that same number of segments and represent each segment by four squares using the algorithm given above.
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 . The solutions for the equation can be given by using parameters like below,
I implemented solutions for slightly more general equation . Solutions for this equation can be constructed from the solutions of the former equation, multiplying each solution of it by .
General sum of n squares
I also implemented solutions for where is an integer. It is obvious that equation is solvable only when .
Lagrange’s four square theorem states 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 can be expressed as a sum of three squares where are non-negative integers. Also any integer which doesn’t contain, in it’s canonical representation, odd powers of a prime of the form can be expressed as a sum of two squares. For example, can’t be expressed as a sum of two squares since it contains and odd power of , which is a prime of the form . But we can express 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,
So, if we can represent each prime divisor of a number as a sum of four squares, we can then use this identity to construct such a representation for. 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 as a sum of squares. Most obvious thing to do is, represent as a sum of four squares and set other variables to zero. We can also partition into approximately 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. 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.
 Lagrange’s four square theorem, [online], Available: http://en.wikipedia.org/wiki/Lagrange%27s_four-square_theorem
 Integer as a sum of three squares, [online], Available: http://www.proofwiki.org/wiki/Integer_as_Sum_of_Three_Squares
 Sum of Squares, [online], Available: http://mathworld.wolfram.com/SumofSquaresFunction.html
 Euler’s four square identity, [online], Available: http://en.wikipedia.org/wiki/Euler%27s_four-square_identity
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 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 can be converted to the form only in the case and is not a perfect square. Here . It occurred to me only when I read  more carefully that any binary quadratic with and 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 is not a perfect square and , we transform the equation to the form . 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 with and . In our form under the case , and , so these two conditions are automatically fulfilled. Below is a rough outline of the Cornacchia’s algorithm. Refer  and  for more details.
Solving by Cornacchia’s method
1. Let be such that .
2. Construct a set such that . We can ommit solutions greater than .
3. For each set .
4. Define the sequences, and for .
5. Find the terms until the condition is met.
6. Set .
7. If and is a perfect square we have the solution and .
8. Add to a set .
9. When we are done with all contains all the solutions to the equation .
Cornacchia’s method finds solutions such that . So it won’t find solutions to the equation since the solutions of this equation is and . So we have to extract out square factors in and find solutions to the equation with this new 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 , i.e for the equation 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 , 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 and solve and separately. After that, we can construct the general solution by adding parameters. The general solution will be constructed after solving and will be added to the solution set for . 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 `test_diophantine.py` 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 `test_diophantine.py` 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.
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 . 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: https://code.google.com/p/sympy/issues/detail?id=3968
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.
 Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003. http://www.jpr2718.org/ax2p.pdf
 Solving ax^2 + by^2 = m by Cornacchia’s method, Number theory web, http://www.numbertheory.org/php/cornacchia.html
 A. Nitaj, L’algorithme de Cornacchia, Expositiones Mathematicae 13 (1995), 358-365
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 is concerned with reducing solutions of the quadratic ternary nonrmal equation . The Holzer’s theorem says that if the above equation is solvable there is always a solution such that . The algorithm for this is explained in . Below is a rough sketch. Before applying the algorithm we have to make square free and we assume that positive and is negative(We can always select in such a way that this is the case). Suppose there is a solution that is not Holzer reduced.
- If is even set and let be any solution to the equation .
- If is odd, c is odd let and let be any solution to the equation .
- Let be the nearest integer to .
- Then the following expressions for gives a solution such that .
Continuing this manner one would find a solution such that .
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 , it will solve and 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.
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 . 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 can be used to improve the speed. Current algorithms for solving binary quadratic equation employs a general solving method so it miss out some speed enhancements for the special case . My plan is if can be converted to the form , then solve it using methods for that, otherwise apply the general method. I found a good reference  for this. A huge thank should go to Pernici for pointing out the improvements.
 Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
 Binary Quadratic Forms: An Algorithmic Approach, J. Buchmann and U. Vollmer, Springer-Verlag, Berlin, Heidelberg, 2007.
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 . Here is a sketch of the algorithm.
1. If then swap and , solve the resulting equation, then swap and in the solution obtained.
2. If then set and stop.
3. If then set and stop.
4. If there is no solution (since must be ).
5. If then set and stop.
6. If then let be a solution of , set and stop.
7. Let be a solution to . with , and set , so that .
8. Use lattice reduction to find a new nontrivial solution to the congruence with as small as possible.
9. Set and write with square-free.
10. Let be a solution to then
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  which uses a method specific to this case.
 .. Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
 .. Gaussian lattice Reduction [online]. Available: http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=404
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, (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 and . In other words, we want to find non-trivial integer solutions to the equation when it has integer coefficients (clearly 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 to a normal form by two rational transformations of the variables. Lets assume that (We can always switch variables if ). With an easy simplification followed by clearing of the denominators , one can easily see that the transformation given by puts our original equation in the form where . Finally, the transformation puts the equation in the desired form, . Since finding integer solutions is equivalent to finding rational solutions and the transformations we used were rational, If an integer solution exist for we can conclude that the new form 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 to have a non-trivial integer solution. If has a non-trivial integer solution, we can use the inverse transformation to recover the solutions for . It turns out that the conditions for solubility of are,
1. All three variables, should not have the same sign.
2. The three quadratic congruence equations, 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 to the form (We can easily do this by multiplying entire equation by ). If then is a solution to our equation. If then is a solution to our equation. Otherwise we continue reducing and 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, which are co-prime to each other. Suppose we have found a non-trivial solution for the equation . We can always assume that by a switching of the variables. Then assuming and substituting in the original equation, one can find an expression for using and . So we can find parametrized solutions for by substituting the value of .
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.
 The algorithmic resolution of Diophantine equations, Nigel P. Smart, London Mathematical Society Student Texts 41, Cambridge University Press, Cambridge, 1998.
 Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.
 L.J. Mordell, Diophantine Equations, Academic Press, 1969.
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.
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  and . 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.
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 . 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 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 .
 Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003. http://www.jpr2718.org/ax2p.pdf
 Solving the generalized Pell equation x**2 – D*y**2 = N, John P. Robertson, July 31, 2004. http://www.jpr2718.org/pell.pdf
 H.W Lenstra, Jr. http://www.ams.org/notices/200202/fea-lenstra.pdf
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.