Trac #28567: allow user to specify strategy for Groebner basis computations with Macaulay2

The {{{groebner_basis}}} method for an ideal in a polynomial ring can be
computed using the default Groebner basis command in Macaulay2. Here we
extend the functionality by allowing the user to specify the strategy to
be used by Macaulay2.

URL: https://trac.sagemath.org/28567
Reported by: saliola
Ticket author(s): Dima Pasechnik, Franco Saliola
Reviewer(s): Markus Wageringel, Dima Pasechnik

Author: Release Manager

src/sage/misc/cachefunc.pyx

@@ -925,7 +925,7 @@ cdef class CachedFunction(object):
             sage: P.<x,y> = QQ[]
             sage: I = P*[x,y]
             sage: from sage.misc.sageinspect import sage_getsourcelines
-            sage: l = "        elif algorithm == 'macaulay2:gb':\n"
+            sage: l = '        elif algorithm.startswith("macaulay2:"):\n'
             sage: l in sage_getsourcelines(I.groebner_basis)[0] # indirect doctest
             True
 

src/sage/rings/polynomial/multi_polynomial_ideal.py

@@ -2997,16 +2997,25 @@ def _macaulay2_(self, macaulay2=None):
         self.__macaulay2[macaulay2] = z
         return z
 
-    def _groebner_basis_macaulay2(self):
+    def _groebner_basis_macaulay2(self, strategy=None):
         r"""
         Return the Groebner basis for this ideal, computed using
         Macaulay2.
 
         ALGORITHM:
 
-        Computed using Macaulay2. A big advantage of Macaulay2 is that
-        it can compute the Groebner basis of ideals in polynomial
-        rings over the integers.
+        Compute the Groebner basis using the specified strategy in Macaulay2.
+        With no strategy option, the Macaulay2 ``gb`` command is used; other
+        possible strategies are "f4" and "mgb", which correspond to the "F4"
+        and "MGB" strategies in Macaulay2.
+
+        A big advantage of Macaulay2 is that it can compute the Groebner basis
+        of ideals in polynomial rings over the integers.
+
+        INPUT:
+
+        - ``strategy`` -- (default: ``'gb'``) argument specifying the strategy
+          to be used by Macaulay2; possiblities: ``'f4'``, ``'gb'``, ``'mgb'``.
 
         EXAMPLES::
 
@@ -3027,11 +3036,41 @@ def _groebner_basis_macaulay2(self):
             sage: I = ideal([y^2*z - x^3 - 19^2*x*z, y^2, 19^2])
             sage: I.groebner_basis('macaulay2')               # optional - macaulay2
             [x^3, y^2, 361]
+
+        Over finite fields, Macaulay2 supports different algorithms to compute
+        Gröbner bases::
+
+            sage: R = PolynomialRing(GF(101), 'x', 4)
+            sage: I = sage.rings.ideal.Cyclic(R)
+            sage: gb1 = I.groebner_basis('macaulay2:gb')  # optional - macaulay2
+            sage: I = sage.rings.ideal.Cyclic(R)
+            sage: gb2 = I.groebner_basis('macaulay2:mgb')  # optional - macaulay2
+            sage: I = sage.rings.ideal.Cyclic(R)
+            sage: gb3 = I.groebner_basis('macaulay2:f4')  # optional - macaulay2
+            sage: gb1 == gb2 == gb3  # optional - macaulay2
+            True
+
+        TESTS::
+
+            sage: R.<x,y,z> = ZZ[]
+            sage: I = ideal([y^2*z - x^3 - 19*x*z, y^2, 19^2])
+            sage: I.groebner_basis('macaulay2:gibberish')     # optional - macaulay2
+            Traceback (most recent call last):
+            ...
+            ValueError: unsupported Macaulay2 strategy
         """
         from sage.rings.polynomial.multi_polynomial_sequence import PolynomialSequence
 
         I = self._macaulay2_()
-        G = str(I.gb().generators().external_string()).replace('\n','')
+        if strategy == "gb" or strategy is None:
+            m2G = I.gb().generators()
+        elif strategy == 'f4':
+            m2G = I.groebnerBasis('Strategy=>"F4"')
+        elif strategy == 'mgb':
+            m2G = I.groebnerBasis('Strategy=>"MGB"')
+        else:
+            raise ValueError("unsupported Macaulay2 strategy")
+        G = str(m2G.external_string()).replace('\n','')
         i = G.rfind('{{')
         j = G.rfind('}}')
         G = G[i+2:j].split(',')
@@ -3817,6 +3856,12 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
         'macaulay2:gb'
             Macaulay2's ``gb`` command (if available)
 
+        'macaulay2:f4'
+            Macaulay2's ``GroebnerBasis`` command with the strategy "F4" (if available)
+
+        'macaulay2:mgb'
+            Macaulay2's ``GroebnerBasis`` command with the strategy "MGB" (if available)
+
         'magma:GroebnerBasis'
             Magma's ``Groebnerbasis`` command (if available)
 
@@ -3916,19 +3961,29 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
             sage: gb == gb.reduced()
             False
 
-        but that ``toy:buchberger2`` does.::
+        but that ``toy:buchberger2`` does. ::
 
             sage: I = sage.rings.ideal.Katsura(P,3) # regenerate to prevent caching
             sage: gb = I.groebner_basis('toy:buchberger2'); gb
             [a - 60*c^3 + 158/7*c^2 + 8/7*c - 1, b + 30*c^3 - 79/7*c^2 + 3/7*c, c^4 - 10/21*c^3 + 1/84*c^2 + 1/84*c]
             sage: gb == gb.reduced()
             True
 
-        ::
+        Here we use Macaulay2 with three different strategies over a finite
+        field. ::
 
-            sage: I = sage.rings.ideal.Katsura(P,3) # regenerate to prevent caching
-            sage: I.groebner_basis('macaulay2:gb') # optional - macaulay2
-            [a - 60*c^3 + 158/7*c^2 + 8/7*c - 1, b + 30*c^3 - 79/7*c^2 + 3/7*c, c^4 - 10/21*c^3 + 1/84*c^2 + 1/84*c]
+            sage: R.<a,b,c> = PolynomialRing(GF(101), 3)
+            sage: I = sage.rings.ideal.Katsura(R,3) # regenerate to prevent caching
+            sage: I.groebner_basis('macaulay2:gb')  # optional - macaulay2
+            [c^3 + 28*c^2 - 37*b + 13*c, b^2 - 41*c^2 + 20*b - 20*c, b*c - 19*c^2 + 10*b + 40*c, a + 2*b + 2*c - 1]
+
+            sage: I = sage.rings.ideal.Katsura(R,3) # regenerate to prevent caching
+            sage: I.groebner_basis('macaulay2:f4')  # optional - macaulay2
+            [c^3 + 28*c^2 - 37*b + 13*c, b^2 - 41*c^2 + 20*b - 20*c, b*c - 19*c^2 + 10*b + 40*c, a + 2*b + 2*c - 1]
+
+            sage: I = sage.rings.ideal.Katsura(R,3) # regenerate to prevent caching
+            sage: I.groebner_basis('macaulay2:mgb') # optional - macaulay2
+            [c^3 + 28*c^2 - 37*b + 13*c, b^2 - 41*c^2 + 20*b - 20*c, b*c - 19*c^2 + 10*b + 40*c, a + 2*b + 2*c - 1]
 
         ::
 
@@ -3937,7 +3992,7 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
             [a - 60*c^3 + 158/7*c^2 + 8/7*c - 1, b + 30*c^3 - 79/7*c^2 + 3/7*c, c^4 - 10/21*c^3 + 1/84*c^2 + 1/84*c]
 
         Singular and libSingular can compute Groebner basis with degree
-        restrictions.::
+        restrictions. ::
 
             sage: R.<x,y> = QQ[]
             sage: I = R*[x^3+y^2,x^2*y+1]
@@ -3980,7 +4035,7 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
         The list of available options is provided at
         :class:`~sage.libs.singular.option.LibSingularOptions`.
 
-        Note that Groebner bases over `\ZZ` can also be computed.::
+        Note that Groebner bases over `\ZZ` can also be computed. ::
 
             sage: P.<a,b,c> = PolynomialRing(ZZ,3)
             sage: I = P * (a + 2*b + 2*c - 1, a^2 - a + 2*b^2 + 2*c^2, 2*a*b + 2*b*c - b)
@@ -4212,8 +4267,8 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
             if prot == "sage":
                 warn("The libsingular interface does not support prot='sage', reverting to 'prot=True'.")
             gb = self._groebner_basis_libsingular(algorithm[len('libsingular:'):], deg_bound=deg_bound, mult_bound=mult_bound, prot=prot, *args, **kwds)
-        elif algorithm == 'macaulay2:gb':
-            gb = self._groebner_basis_macaulay2(*args, **kwds)
+        elif algorithm.startswith("macaulay2:"):
+            gb = self._groebner_basis_macaulay2(strategy=algorithm.split(":")[1], *args, **kwds)
         elif algorithm == 'magma:GroebnerBasis':
             gb = self._groebner_basis_magma(prot=prot, deg_bound=deg_bound, *args, **kwds)
         elif algorithm == 'toy:buchberger':
@@ -4745,7 +4800,7 @@ def random_element(self, degree, compute_gb=False, *args, **kwds):
 
         EXAMPLES:
 
-        We compute a uniformly random element up to the provided degree.::
+        We compute a uniformly random element up to the provided degree. ::
 
             sage: P.<x,y,z> = GF(127)[]
             sage: I = sage.rings.ideal.Katsura(P)