Repository: sympy/sympy

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

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 55b874e6f8..86b8e0da17 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1550,6 +1550,8 @@ 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 <123136644+zylipku@users.noreply.github.com>
+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 3df4fe9b0d..d02d0689d8 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 a == 1:
+            # y = 0
+            s, _exact = iroot(m // a, 2)
+            if _exact:
+                sols.add((int(s), 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((0, int(s)))
+        return sols
+
+    # the original cornacchia
     for t in sqrt_mod_iter(-b*invert(a, m), m):
         if t < m // 2:
             continue


Test patch:
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index 094770b7bb..a500903fc0 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)}