Bug in polynomial GCD computation

[bluescarni]

This is incorrect:

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

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

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

[smichr]

Thanks for the valuable feedback! It looks like the algorithm for heuristic gcd is describeed on pages 320 ff. I can't see pages 327-330, however. Someone else will have to see if they can get the reference or some other way address the issue by finding the problem in the existing SymPy code..

[bluescarni]

@smichr No problem. I have an electronic copy of the book (legally obtained). Do you think I can extract the subsection describing the algorithm and send it to you? Not sure if it falls under any fair use clause.

[jksuom]

It seems that the version of heuristic gcd implemented in SymPy makes use of following upper bound for the roots r of a polynomial f:
|r| < f_norm / |f.LC| + 1,
where f_norm is the maximum absolute value of the coefficients and f.LC is the leading coefficient. However, this is not exactly how it is implemented. Instead of the rational quotient f_norm / |f.LC| the truncated integer quotient f_norm // |f.LC| is used. Hence it is necessary to add 1 to preserve the validity of the strict inequality:
|r| < f_norm // |f.LC| + 2.

In the code the test values have to be chosen strictly greater than twice the absolute values of the roots. Hence the choice should be modified:

--- a/sympy/polys/heuristicgcd.py
+++ b/sympy/polys/heuristicgcd.py
@@ -69,7 +69,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)