GCD of polyomials over polynomial rings

[Bruno Grenet]

In this ticket, we add a mechanism for computing GCD of polynomials over polynomial ring: For instance, to compute the gcd of polynomials in x over QQ['y'], we use the implementation of GCD in QQ['x,y']. This fasten a lot the computations.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​df736ce	#20220: use multivariate gcd for univariate polys over polyrings
​269d10a	#20220: use multivariate gcd for multivariate polys over polyrings
​e83510a	#20220: Correct indentation (unrelated)

[Bruno Grenet]

Note: third commit simply correct an indentation error in the file src/sage/rings/polynomial/multi_polynomial.pyx: The content of the method inverse_mod used 3 spaces rather than 4 for indentation.

[Alyson Deines]

The following doctest now fail:

aly@aly-laptop:~/Sage/sage$ ./sage -t src/sage/schemes/elliptic_curves/ell_curve_isogeny.py
Running doctests with ID 2016-03-20-12-44-23-b7ac0a73.
Git branch: ticket/20220
Using --optional=ccache,mpir,python2,sage
Doctesting 1 file.
sage -t --warn-long 17.3 src/sage/schemes/elliptic_curves/ell_curve_isogeny.py
**********************************************************************
File "src/sage/schemes/elliptic_curves/ell_curve_isogeny.py", line 847, in sage.schemes.elliptic_curves.ell_curve_isogeny.EllipticCurveIsogeny
Failed example:
    isogs[0].rational_maps()
Expected:
    ((x^2 - t^2)/x, (x^3*y + t^2*x*y)/x^3)
Got:
    ((x^2 - t^2)/x, (x^2*y + t^2*y)/x^2)
**********************************************************************
File "src/sage/schemes/elliptic_curves/ell_curve_isogeny.py", line 852, in 
sage.schemes.elliptic_curves.ell_curve_isogeny.EllipticCurveIsogeny
Failed example:
    duals[0].rational_maps()
Expected:
    ((1/4*x^2 + t^2)/x, (1/8*x^3*y + (-1/2*t^2)*x*y)/x^3)
Got:
    ((1/4*x^2 + t^2)/x, (1/8*x^2*y + (-1/2*t^2)*y)/x^2)
**********************************************************************
1 item had failures:
   2 of 139 in sage.schemes.elliptic_curves.ell_curve_isogeny.EllipticCurveIsogeny
    [979 tests, 2 failures, 19.24 s]
----------------------------------------------------------------------
sage -t --warn-long 17.3 src/sage/schemes/elliptic_curves/ell_curve_isogeny.py  # 2 doctests failed
----------------------------------------------------------------------
Total time for all tests: 19.4 seconds
    cpu time: 19.2 seconds
    cumulative wall time: 19.2 seconds

and

aly@aly-laptop:~/Sage/sage$ ./sage -t src/sage/schemes/affine/affine_morphism.py 
Running doctests with ID 2016-03-20-12-44-52-88984cb0.
Git branch: ticket/20220
Using --optional=ccache,mpir,python2,sage
Doctesting 1 file.
sage -t --warn-long 17.3 src/sage/schemes/affine/affine_morphism.py
**********************************************************************
File "src/sage/schemes/affine/affine_morphism.py", line 418, in sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.homogenize
Failed example:
    f.homogenize((2, 0))
Exception raised:
    Traceback (most recent call last):
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.homogenize[7]>", line 1, in <module>
        f.homogenize((Integer(2), Integer(0)))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/schemes/affine/affine_morphism.py", line 526, in homogenize
        g = gcd(F)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1614, in gcd
        return __GCD_sequence(seq, **kwargs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1656, in __GCD_sequence
        g = vi.gcd(g, **kwargs)
      File "sage/rings/polynomial/multi_polynomial.pyx", line 1758, in sage.rings.polynomial.multi_polynomial.MPolynomial.gcd (build/cythonized/sage/rings/polynomial/multi_polynomial.c:17781)
        return self._parent(unibase._gcd_univariate_polynomial(uniself, other.polynomial(x)))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/categories/unique_factorization_domains.py", line 181, in _gcd_univariate_polynomial
        d = a.gcd(b)
      File "sage/rings/polynomial/multi_polynomial.pyx", line 1758, in sage.rings.polynomial.multi_polynomial.MPolynomial.gcd (build/cythonized/sage/rings/polynomial/multi_polynomial.c:17781)
        return self._parent(unibase._gcd_univariate_polynomial(uniself, other.polynomial(x)))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/categories/unique_factorization_domains.py", line 182, in _gcd_univariate_polynomial
        A = parent(A/a)
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 439, in _element_constructor_
        return C(self, x, check, is_gen, construct=construct, **kwds)
      File "sage/rings/polynomial/polynomial_element.pyx", line 7974, in sage.rings.polynomial.polynomial_element.Polynomial_generic_dense.__init__ (build/cythonized/sage/rings/polynomial/polynomial_element.c:68167)
        self.__coeffs = [R(a, **kwds) for a in x.list()]
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 426, in _element_constructor_
        raise TypeError("denominator must be a unit")
    TypeError: denominator must be a unit
**********************************************************************
File "src/sage/schemes/affine/affine_morphism.py", line 583, in sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.dynatomic_polynomial
Failed example:
    f.dynatomic_polynomial(3)
Exception raised:
    Traceback (most recent call last):
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.dynatomic_polynomial[11]>", line 1, in <module>
        f.dynatomic_polynomial(Integer(3))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/schemes/affine/affine_morphism.py", line 621, in dynatomic_polynomial
        F = self.homogenize(1).dynatomic_polynomial(period)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/schemes/affine/affine_morphism.py", line 526, in homogenize
        g = gcd(F)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1614, in gcd
        return __GCD_sequence(seq, **kwargs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1656, in __GCD_sequence
        g = vi.gcd(g, **kwargs)
      File "sage/rings/polynomial/multi_polynomial.pyx", line 1758, in sage.rings.polynomial.multi_polynomial.MPolynomial.gcd (build/cythonized/sage/rings/polynomial/multi_polynomial.c:17781)
        return self._parent(unibase._gcd_univariate_polynomial(uniself, other.polynomial(x)))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/categories/unique_factorization_domains.py", line 182, in _gcd_univariate_polynomial
        A = parent(A/a)
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 439, in _element_constructor_
        return C(self, x, check, is_gen, construct=construct, **kwds)
      File "sage/rings/polynomial/polynomial_element.pyx", line 7974, in sage.rings.polynomial.polynomial_element.Polynomial_generic_dense.__init__ (build/cythonized/sage/rings/polynomial/polynomial_element.c:68167)
        self.__coeffs = [R(a, **kwds) for a in x.list()]
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 426, in _element_constructor_
        raise TypeError("denominator must be a unit")
    TypeError: denominator must be a unit
**********************************************************************
File "src/sage/schemes/affine/affine_morphism.py", line 600, in sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.dynatomic_polynomial
Failed example:
    f.dynatomic_polynomial(2)
Exception raised:
    Traceback (most recent call last):
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 496, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 858, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.dynatomic_polynomial[19]>", line 1, in <module>
        f.dynatomic_polynomial(Integer(2))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/schemes/affine/affine_morphism.py", line 621, in dynatomic_polynomial
        F = self.homogenize(1).dynatomic_polynomial(period)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/schemes/affine/affine_morphism.py", line 526, in homogenize
        g = gcd(F)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1614, in gcd
        return __GCD_sequence(seq, **kwargs)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/arith/misc.py", line 1656, in __GCD_sequence
        g = vi.gcd(g, **kwargs)
      File "sage/rings/polynomial/multi_polynomial.pyx", line 1758, in sage.rings.polynomial.multi_polynomial.MPolynomial.gcd (build/cythonized/sage/rings/polynomial/multi_polynomial.c:17781)
        return self._parent(unibase._gcd_univariate_polynomial(uniself, other.polynomial(x)))
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/categories/unique_factorization_domains.py", line 183, in _gcd_univariate_polynomial
        B = parent(B/b)
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 439, in _element_constructor_
        return C(self, x, check, is_gen, construct=construct, **kwds)
      File "sage/rings/polynomial/polynomial_element.pyx", line 7974, in sage.rings.polynomial.polynomial_element.Polynomial_generic_dense.__init__ (build/cythonized/sage/rings/polynomial/polynomial_element.c:68167)
        self.__coeffs = [R(a, **kwds) for a in x.list()]
      File "sage/structure/parent.pyx", line 1111, in sage.structure.parent.Parent.__call__ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/parent.c:9828)
        return mor._call_(x)
      File "sage/structure/coerce_maps.pyx", line 109, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4542)
        raise
      File "sage/structure/coerce_maps.pyx", line 104, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (/home/aly/Sage/sage/src/build/cythonized/sage/structure/coerce_maps.c:4435)
        return C._element_constructor(x)
      File "/home/aly/Sage/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_ring.py", line 426, in _element_constructor_
        raise TypeError("denominator must be a unit")
    TypeError: denominator must be a unit
**********************************************************************
2 items had failures:
   2 of  25 in sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.dynatomic_polynomial
   1 of  31 in sage.schemes.affine.affine_morphism.SchemeMorphism_polynomial_affine_space.homogenize
    [268 tests, 3 failures, 0.81 s]
----------------------------------------------------------------------
sage -t --warn-long 17.3 src/sage/schemes/affine/affine_morphism.py  # 3 doctests failed
----------------------------------------------------------------------
Total time for all tests: 0.9 seconds
    cpu time: 0.8 seconds
    cumulative wall time: 0.8 seconds
aly@aly-laptop:~/Sage/sage$ 

[Bruno Grenet]

Thanks for noticing this, and my apologies for not having run the tests myself!

The first failed doctest is not a real failure: I would say that my changes improve the situation since the rational functions in this test were not reduced, and they are now.

The second failed doctest is a real failure. It requires a change in the method _gcd_univariate_polynomial of the category UniqueFactorizationDomains (replace truediv by floordiv). This change breaks other doctests. tl;dr: I am working on it!

[edit] There is another change to make: Right now, for polynomial rings such as R[x][y][z] where R is some "uncommon" ring (not ZZ, QQ, but CC or Frac(QQ['z']) for instance), the algorithm first converts the input to R[x,y,z], and then in the multivariate algorithm converts back to R[x][y][z]. This works, but is of course silly. I am also trying to find the best solution to avoid the problem!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​e943a17	#20220: Replace truediv by floordiv in _gcd_univariate_polynomial
​1670b32	#20220: Correct failings doctest in schemes/elliptic_curves/ell_curve_isogeny.py

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​2fb1836

#20220: use singular for multivariate gcd

[Bruno Grenet]

My three commits should correct the problems noticed above. The problem mentioned in the "[edit]" part of comment:5 is corrected by the latest commit: Instead of trying to use the generic method gcd for multivariate polynomials, it directly tests whether Singular can perform the computation.

My main concern now is performance: At some point, I run tests that were very slow, but I am not able to reproduce the behavior anymore. That's why I put Needs Review: Please perform tests to see whether my changes imply slowdown in some parts. It should'nt...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​fec268c

#20220: Better use of Singular

[Bruno Grenet]

Replying to bruno:

My main concern now is performance: At some point, I run tests that were very slow, but I am not able to reproduce the behavior anymore. That's why I put Needs Review: Please perform tests to see whether my changes imply slowdown in some parts. It should'nt...

This is not valid (anymore)!

[Alyson Deines]

Everything looks good. I'm not seeing any slow doctests (besides the ones we already know about!)