Repository: sympy/sympy

Issue:
Bug in polynomial GCD computation
This is incorrect:

``` python
>>> from sympy import Symbol as S, gcd
>>> x,y,z = S('x'),S('y'),S('z')
>>> num = 12*x**6*y**7*z**3 - 3*x**4*y**9*z**3 + 12*x**3*y**5*z**4 
>>> den = -48*x**7*y**8*z**3 + 12*x**5*y**10*z**3 - 48*x**5*y**7*z**2 + 36*x**4*y**7*z - 48*x**4*y**6*z**4 + 12*x**3*y**9*z**2 - 48*x**3*y**4 - 9*x**2*y**9*z - 48*x**2*y**5*z**3 + 12*x*y**6 + 36*x*y**5*z**2 - 48*y**2*z
>>> gcd(num,den)
3*y**2
```

The correct result is `-12*y**2*z-12*x**3*y**4+3*x*y**6`. This boils down to an error in the heuristic GCD. Indeed, the modular GCD works:

``` python
>>> from sympy.polys.modulargcd import modgcd_multivariate
>>> from sympy.polys.rings import ring
>>> from sympy.polys.domains import ZZ
>>> R, x, y, z = ring("x,y,z", ZZ)
>>> num = 12*x**6*y**7*z**3 - 3*x**4*y**9*z**3 + 12*x**3*y**5*z**4
>>> den = -48*x**7*y**8*z**3 + 12*x**5*y**10*z**3 - 48*x**5*y**7*z**2 + 36*x**4*y**7*z - 48*x**4*y**6*z**4 + 12*x**3*y**9*z**2 - 48*x**3*y**4 - 9*x**2*y**9*z - 48*x**2*y**5*z**3 + 12*x*y**6 + 36*x*y**5*z**2 - 48*y**2*z
>>> modgcd_multivariate(num,den) 
(12*x**3*y**4 - 3*x*y**6 + 12*y**2*z,
 x**3*y**3*z**3,
 -4*x**4*y**4*z**3 - 4*x**2*y**3*z**2 + 3*x*y**3*z - 4)
```

The heuristic GCD does not:

``` python
>>> from sympy.polys.heuristicgcd import heugcd
>>> heugcd(num,den)
(3*y**2,
 4*x**6*y**5*z**3 - x**4*y**7*z**3 + 4*x**3*y**3*z**4,
 -16*x**7*y**6*z**3 + 4*x**5*y**8*z**3 - 16*x**5*y**5*z**2 + 12*x**4*y**5*z - 16*x**4*y**4*z**4 + 4*x**3*y**7*z**2 - 16*x**3*y**2 - 3*x**2*y**7*z - 16*x**2*y**3*z**3 + 4*x*y**4 + 12*x*y**3*z**2 - 16*z)
```

I encountered the bug while implementing heuristic GCD in piranha (https://github.com/bluescarni/piranha). The test case above was generated randomly by fuzzy testing and it is the only one that fails. When I used sympy as a verification tool, sympy gave the same (wrong) result. It turns out that I used the same paper as sympy for the implementation of the heuristic GCD:

http://dl.acm.org/citation.cfm?id=220376

So at the very least it looks like in the paper the algorithm has been described in a way which leads to incorrect implementations. In piranha I ended up using the implementation of heuristic GCD as described in the book by Geddes:

https://books.google.de/books/about/Algorithms_for_Computer_Algebra.html?id=TXfY891RiiEC&redir_esc=y

This works so far in all test cases.


Discussion:

Patch:
diff --git a/sympy/polys/heuristicgcd.py b/sympy/polys/heuristicgcd.py
index e2cb12a013..6a25fb215b 100644
--- a/sympy/polys/heuristicgcd.py
+++ b/sympy/polys/heuristicgcd.py
@@ -68,7 +68,7 @@ def heugcd(f, g):
 
     x = max(min(B, 99*domain.sqrt(B)),
             2*min(f_norm // abs(f.LC),
-                  g_norm // abs(g.LC)) + 2)
+                  g_norm // abs(g.LC)) + 4)
 
     for i in range(0, HEU_GCD_MAX):
         ff = f.evaluate(x0, x)


Test patch:
diff --git a/sympy/polys/tests/test_heuristicgcd.py b/sympy/polys/tests/test_heuristicgcd.py
index cc2d8438bb..4f2d6c8f55 100644
--- a/sympy/polys/tests/test_heuristicgcd.py
+++ b/sympy/polys/tests/test_heuristicgcd.py
@@ -126,3 +126,16 @@ def test_heugcd_multivariate_integers():
     H, cff, cfg = heugcd(f, g)
 
     assert H == h and H*cff == f and H*cfg == g
+
+def test_issue_10996():
+    R, x, y, z = ring("x,y,z", ZZ)
+
+    f = 12*x**6*y**7*z**3 - 3*x**4*y**9*z**3 + 12*x**3*y**5*z**4
+    g = -48*x**7*y**8*z**3 + 12*x**5*y**10*z**3 - 48*x**5*y**7*z**2 + \
+    36*x**4*y**7*z - 48*x**4*y**6*z**4 + 12*x**3*y**9*z**2 - 48*x**3*y**4 \
+    - 9*x**2*y**9*z - 48*x**2*y**5*z**3 + 12*x*y**6 + 36*x*y**5*z**2 - 48*y**2*z
+
+    H, cff, cfg = heugcd(f, g)
+
+    assert H == 12*x**3*y**4 - 3*x*y**6 + 12*y**2*z
+    assert H*cff == f and H*cfg == g