# 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 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.

#### 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: 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.

#### References

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