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

Advertisements

2 Comments

  1. I don’t remember proposing that…

  2. Karen

    Wow I looked up the name Cornacchia because that is my last name and found you. I find it amazing that someone can understand math like Giuseppe Cornacchia . The mind is an amazing thing. Sympy…. I will have to look this word up in the dictionary.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: