Repository: sympy/sympy
Pull request: Fix a bug with diophantine() when cornacchia() is not called properly
Occasionally diophantine() and diop_quadratic() result in a DN form $X^2 - DY^2 = N$ which cannot be handled by the current cornacchia() properly. For example, previously the diophantine solver will miss some solutions:
```
>>> eq1 = x**2 - 15*x + y**2 - 8*y # also mentioned in #18628
>>> diophantine(eq1)
{(0, 0), (15, 0), (15, 8), (0, 8)}
>>> eq1.subs({x: 16, y: 4})
0
```
```
>>> eq2 = 2*x**2 - 9*x + 4*y**2 - 8*y + 14
>>> diophantine(eq2)
set()
>>> eq2.subs({x: 2, y: 1})
0
```
Both the two examples lead to a function call of cornacchia(1, X, 1) with X > 1, which incorrrectly turns out to have no solution by using cornacchia(). The solver can only be applied when $a+b\le m$. See reference from docs:
http://www.numbertheory.org/php/cornacchia.html 
An additional check has been added before the algorithm based on the fact that $ax^2 + by^2 = m$ admits a solution for $a + b > m$ only if $x = 0$ or $y = 0$. A test function including both the two examples has been appended as well.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #18628 

#### Brief description of what is fixed or changed
An additional check has been added before the cornacchia algorithm body based on the fact that $ax^2 + by^2 = m$ admits a solution for $a + b > m$ only if $x = 0$ or $y = 0$.

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * Fix a bug with cornacchia(). Formerly diophantine() and diop_quadratic() may miss some solutions for quadratic equations.
<!-- END RELEASE NOTES -->


Problem statement:
diophantine(), diop_quadratic() missing solutions
With sympy version 1.5.1:

```
from sympy import *
from sympy.solvers.diophantine import diop_quadratic

a, b, x, y = symbols('a b x y', integer=True, positive=True)
f = -a*x + x**2 - b*y + y**2
g = f.subs([(a, 15), (b, 8)])
print('diophantine   (', g, ') solutions:', diophantine(g))
print('diop_quadratic(', g, ') solutions:', diop_quadratic(g))
h = f.subs([(a, 8), (b, 15)])
print('diophantine   (', h, ') solutions:', diophantine(h))
print('diop_quadratic(', h, ') solutions:', diop_quadratic(h))
print(g.subs([(x, 16), (y, 4)]))

```
Output:

```
diophantine   ( x**2 - 15*x + y**2 - 8*y ) solutions: {(15, 8)}
diop_quadratic( x**2 - 15*x + y**2 - 8*y ) solutions: {(15, 0), (0, 8), (15, 8), (0, 0)}
diophantine   ( x**2 - 8*x + y**2 - 15*y ) solutions: {(8, 15), (4, 16)}
diop_quadratic( x**2 - 8*x + y**2 - 15*y ) solutions: {(0, 0), (0, 15), (8, 0), (4, 16), (4, -1), (8, 15)}
0

```
First, I expect diophantine() to return the same solutions as diop_quadratic(), but it doesn't.

Second, the function f() is symmetric in x and y, so the solutions to g() and h() should be symmetric, but they aren't.  For example, (16, 4) is a solution to g(), but is missing from the solution set.

Hints from issue discussion:
`diophantine` is doing additional filtering on the results to match the assumptions on `x` and `y` as positive integers. `diop_quadratic` has no such checks.

The boundary between public and private methods in the diophantine solver is extremely ill-defined.

I agree that there are solutions missing from your example.
I see, thanks.  Zero isn't considered positive, so I should have:

`a, b, x, y = symbols('a b x y', integer=True, negative=False)
`

Then the output is:

```
diophantine   ( x**2 - 15*x + y**2 - 8*y ) solutions: {(15, 0), (0, 8), (15, 8), (0, 0)}
diop_quadratic( x**2 - 15*x + y**2 - 8*y ) solutions: {(15, 0), (0, 8), (15, 8), (0, 0)}
diophantine   ( x**2 - 8*x + y**2 - 15*y ) solutions: {(0, 0), (0, 15), (8, 0), (4, 16), (8, 15)}
diop_quadratic( x**2 - 8*x + y**2 - 15*y ) solutions: {(0, 0), (0, 15), (8, 0), (4, 16), (4, -1), (8, 15)}
0

```
So the issue is with diop_quadratic(g), which is missing the solution (16, 4).
I've investigated this further.  The problem appears to actually be in `transformation_to_DN` and `find_DN`:

```
from sympy import *
from sympy.solvers.diophantine import transformation_to_DN, find_DN, diop_DN

d, e, x, y = symbols('d e x y', integer=True)

f = x**2 + y**2 + d*x + e*y         
g = f.subs([(d, -15), (e, -8)])
h = f.subs([(d, -8), (e, -15)])

print(f'transformation_to_DN({g}) =')
A, B = transformation_to_DN(g)
print('A:', A)
print('B:', B)
D, N = find_DN(g)
print(f'find_DN({g}) = {D}, {N}')
print(f'diop_DN({D}, {N}) = {diop_DN(D, N)}')
print('')

print(f'transformation_to_DN({h}) =')
A, B = transformation_to_DN(h)
print('A:', A)
print('B:', B)
D, N = find_DN(h)
print(f'find_DN({h}) = {D}, {N}')
print(f'diop_DN({D}, {N}) = {diop_DN(D, N)}')
```

Output:

```
transformation_to_DN(x**2 - 15*x + y**2 - 8*y) =
A: Matrix([[-1/2, 0], [0, -1]])
B: Matrix([[15/2], [4]])
find_DN(x**2 - 15*x + y**2 - 8*y) = -4, 289
diop_DN(-4, 289) = [(15, 4)]

transformation_to_DN(x**2 - 8*x + y**2 - 15*y) =
A: Matrix([[-1/4, 0], [0, -1/2]])
B: Matrix([[4], [15/2]])
find_DN(x**2 - 8*x + y**2 - 15*y) = -4, 1156
diop_DN(-4, 1156) = [(16, 15), (30, 8), (0, 17)]

```
My math isn't strong, but shouldn't `transformation_to_DN(x**2 - 15*x + y**2 - 8*y)` return:

`A: Matrix([[-1/2, 0], [0, -1/4]])
`

and `find_DN(x**2 - 15*x + y**2 - 8*y)` also return `-4, 1156`?

I don't understand the method used by these functions, and the link to the paper _Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003_  ([http://www.jpr2718.org/ax2p.pdf](url)) is broken.

[Here's a link that works: [https://vdocuments.mx/solving-the-equation-ax2-bxy-cy2-dx-ey-f-0.html](url)]
The transformation is not unique. It seems that the one defined by `Matrix([[-1/2, 0], [0, -1]])` would also be valid. However, `diop_DN(-4, 289)` does not return the non-primitive solution `(17, 0)` of `x**2 + 4*y**2 = 289`because `cornacchia` does not find the trivial solution `(1, 0)` of the equation `x**2 + 4*y**2 = 1`:
```
>>> from sympy.solvers.diophantine.diophantine import cornacchia
>>> cornacchia(1, 4, 1)
set()
```
That looks like a bug in `cornacchia`.
From [http://www.numbertheory.org/php/cornacchia.html](url), it appears that `cornacchia` isn't applicable because a + b > m.
From the documentation for `cornacchia`:

Solving the diophantine equation ax**2 + by**2 = m by Cornacchia’s method, [online], Available: [http://www.numbertheory.org/php/cornacchia.html](http://www.numbertheory.org/php/cornacchia.html)

The page at that link states: _Here a > 0, b > 0,m ≥ a+b, gcd(a,m)=1=gcd(a,b)_, so Cornacchia’s method fails because `1 + 4 > 1`, making me think that `cornacchia `isn't applicable in this case.
> From the documentation for `cornacchia`:
> 
> Solving the diophantine equation ax**2 + by**2 = m by Cornacchia’s method, [online], Available: http://www.numbertheory.org/php/cornacchia.html
> 
> The page at that link states: _Here a > 0, b > 0,m ≥ a+b, gcd(a,m)=1=gcd(a,b)_, so Cornacchia’s method fails because `1 + 4 > 1`, making me think that `cornacchia `isn't applicable in this case.

I've redacted my statements above since I now realize that I did not fully understand the constraints on the coefficients in the equation <code>a*x^2 + b*y^2 = m</code> when solving it using Cornacchia’s method. Thank you for clarifying, I now understand the problem much better. 

> The transformation is not unique. It seems that the one defined by `Matrix([[-1/2, 0], [0, -1]])` would also be valid. However, `diop_DN(-4, 289)` does not return the non-primitive solution `(17, 0)` of `x**2 + 4*y**2 = 289`because `cornacchia` does not find the trivial solution `(1, 0)` of the equation `x**2 + 4*y**2 = 1`:
> 
> ```
> >>> from sympy.solvers.diophantine.diophantine import cornacchia
> >>> cornacchia(1, 4, 1)
> set()
> ```
> 
> That looks like a bug in `cornacchia`.

There is not a bug in the <code>cornacchia</code> method (at least not for this reason). @kgorlen is correct, we can not use Cornacchia’s method in this instance since a + b > m.

Commit message: fix typos in the comments of cornacchia()

Patch to review:
diff --git a/.mailmap b/.mailmap
index 55b874e6f897..d67a21ac36df 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1550,6 +1550,7 @@ Zeel Shah <kshah215@gmail.com>
 Zhenxu Zhu <xzdlj@outlook.com> xzdlj <xzdlj@outlook.com>
 Zhi-Qiang Zhou <zzq_890709@hotmail.com> zhouzq-thu <zzq_890709@hotmail.com>
 Zhongshi <zj495@nyu.edu>
+Zhuoyuan Li <zy.li@stu.pku.edu.cn> zylipku <zy.li@stu.pku.edu.cn>
 Zlatan Vasović <zlatanvasovic@gmail.com>
 Zoufiné Lauer-Baré <raszoufine@gmail.com> Zoufiné Lauer-Baré <82505312+zolabar@users.noreply.github.com>
 Zoufiné Lauer-Baré <raszoufine@gmail.com> zolabar <raszoufine@gmail.com>
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index 3df4fe9b0df1..fc38c7f19d35 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -2174,6 +2174,25 @@ def cornacchia(a:int, b:int, m:int) -> set[tuple[int, int]]:
     """
     # Assume gcd(a, b) = gcd(a, m) = 1 and a, b > 0 but no error checking
     sols = set()
+
+    if a + b > m:
+        # xy = 0 must hold if there exists a solution
+        if m % a == 0:
+            # y = 0
+            s, _exact = iroot(m // a, 2)
+            if _exact:
+                sols.add((int(s), int(0)))
+            if a == b:
+            # only keep one solution
+                return sols
+        if m % b == 0:
+            # x = 0
+            s, _exact = iroot(m // b, 2)
+            if _exact:
+                sols.add((int(0), int(s)))
+        return sols
+
+    # the original cornacchia
     for t in sqrt_mod_iter(-b*invert(a, m), m):
         if t < m // 2:
             continue
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index 094770b7bba7..a500903fc06f 100644
--- a/sympy/solvers/diophantine/tests/test_diophantine.py
+++ b/sympy/solvers/diophantine/tests/test_diophantine.py
@@ -1049,3 +1049,11 @@ def test_quadratic_parameter_passing():
     # test that parameters are passed all the way to the final solution
     assert solution == {(t, 11*t), (t, -22*t)}
     assert solution(0, 0) == {(0, 0)}
+
+def test_issue_18628():
+    eq1 = x**2 - 15*x + y**2 - 8*y
+    sol = diophantine(eq1)
+    assert sol == {(15, 0), (15, 8), (-1, 4), (0, 0), (0, 8), (16, 4)}
+    eq2 = 2*x**2 - 9*x + 4*y**2 - 8*y + 14
+    sol = diophantine(eq2)
+    assert sol == {(2, 1)}


Code review comments:
diff_hunk: @@ -2174,6 +2174,25 @@ def cornacchia(a:int, b:int, m:int) -> set[tuple[int, int]]:
     """
     # Assume gcd(a, b) = gcd(a, m) = 1 and a, b > 0 but no error checking
     sols = set()
+
+    if a + b > m:
+        # xy = 0 must hold if there exists a solution
+        if m % a == 0:
+            # y = 0
+            s, _exact = iroot(m // a, 2)
+            if _exact:
+                sols.add((int(s), int(0)))
+            if a == b:
+            # only keep one solution
+                return sols
+        if m % b == 0:
+            # x = 0
+            s, _exact = iroot(m // b, 2)
+            if _exact:
+                sols.add((int(0), int(s)))
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
```suggestion
                sols.add((0, int(s)))
```

diff_hunk: @@ -2174,6 +2174,25 @@ def cornacchia(a:int, b:int, m:int) -> set[tuple[int, int]]:
     """
     # Assume gcd(a, b) = gcd(a, m) = 1 and a, b > 0 but no error checking
     sols = set()
+
+    if a + b > m:
+        # xy = 0 must hold if there exists a solution
+        if m % a == 0:
+            # y = 0
+            s, _exact = iroot(m // a, 2)
+            if _exact:
+                sols.add((int(s), int(0)))
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
```suggestion
                sols.add((int(s), 0))
```

diff_hunk: @@ -2174,6 +2174,25 @@ def cornacchia(a:int, b:int, m:int) -> set[tuple[int, int]]:
     """
     # Assume gcd(a, b) = gcd(a, m) = 1 and a, b > 0 but no error checking
     sols = set()
+
+    if a + b > m:
+        # xy = 0 must hold if there exists a solution
+        if m % a == 0:
+            # y = 0
+            s, _exact = iroot(m // a, 2)
+            if _exact:
+                sols.add((int(s), int(0)))
+            if a == b:
+            # only keep one solution
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
```suggestion
                # only keep one solution
```

diff_hunk: @@ -2174,6 +2174,25 @@ def cornacchia(a:int, b:int, m:int) -> set[tuple[int, int]]:
     """
     # Assume gcd(a, b) = gcd(a, m) = 1 and a, b > 0 but no error checking
     sols = set()
+
+    if a + b > m:
+        # xy = 0 must hold if there exists a solution
+        if m % a == 0:
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
Given the constraint `gcd(a, m) = 1`, I think the condition `if a == 1` should be sufficient.