Merge pull request sympy#26780 from zylipku/fix-cornacchia-miss

Fix a bug with diophantine() when cornacchia() is not called properly

Author: oscarbenjamin

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

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

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)}