Merge #34519

src/sage/rings/polynomial/msolve.py

@@ -28,14 +28,93 @@
 from sage.misc.converting_dict import KeyConvertingDict
 from sage.misc.sage_eval import sage_eval
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+from sage.rings.finite_rings.finite_field_base import FiniteField
 from sage.rings.rational_field import QQ
 from sage.rings.real_arb import RealBallField
 from sage.rings.real_double import RealDoubleField_class
 from sage.rings.real_mpfr import RealField_class
 from sage.rings.real_mpfi import RealIntervalField_class, RealIntervalField
+from sage.structure.sequence import Sequence
 
+def _run_msolve(ideal, options):
+    r"""
+    Internal utility function
+    """
+
+    base = ideal.base_ring()
+    if not (base is QQ or isinstance(base, FiniteField) and
+            base.is_prime_field() and base.characteristic() < 2**31):
+        raise NotImplementedError(f"unsupported base field: {base}")
 
-def _variety(ideal, ring, proof):
+    # Run msolve
+
+    msolve().require()
+
+    drlpolring = ideal.ring().change_ring(order='degrevlex')
+    polys = ideal.change_ring(drlpolring).gens()
+    msolve_in = tempfile.NamedTemporaryFile(mode='w',
+                                            encoding='ascii', delete=False)
+    command = ["msolve", "-f", msolve_in.name] + options
+    try:
+        print(",".join(drlpolring.variable_names()), file=msolve_in)
+        print(base.characteristic(), file=msolve_in)
+        print(*(pol._repr_().replace(" ", "") for pol in polys),
+                sep=',\n', file=msolve_in)
+        msolve_in.close()
+        msolve_out = subprocess.run(command, capture_output=True, text=True)
+    finally:
+        os.unlink(msolve_in.name)
+    msolve_out.check_returncode()
+
+    return msolve_out.stdout
+
+def groebner_basis_degrevlex(ideal):
+    r"""
+    Compute a degrevlex Gröbner basis using msolve
+
+    EXAMPLES::
+
+        sage: from sage.rings.polynomial.msolve import groebner_basis_degrevlex
+
+        sage: R.<a,b,c> = PolynomialRing(GF(101), 3, order='lex')
+        sage: I = sage.rings.ideal.Katsura(R,3)
+        sage: gb = groebner_basis_degrevlex(I); gb # optional - msolve
+        [a + 2*b + 2*c - 1, b*c - 19*c^2 + 10*b + 40*c,
+        b^2 - 41*c^2 + 20*b - 20*c, c^3 + 28*c^2 - 37*b + 13*c]
+        sage: gb.universe() is R
+        False
+        sage: gb.universe().term_order()
+        Degree reverse lexicographic term order
+        sage: ideal(gb).transformed_basis(other_ring=R)
+        [c^4 + 38*c^3 - 6*c^2 - 6*c, 30*c^3 + 32*c^2 + b - 14*c,
+        a + 2*b + 2*c - 1]
+
+    TESTS::
+
+        sage: R.<foo, bar> = PolynomialRing(GF(536870909), 2)
+        sage: I = Ideal([ foo^2 - 1, bar^2 - 1 ])
+        sage: I.groebner_basis(algorithm='msolve') # optional - msolve
+        [bar^2 - 1, foo^2 - 1]
+
+        sage: R.<x, y> = PolynomialRing(QQ, 2)
+        sage: I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
+        sage: I.groebner_basis(algorithm='msolve') # optional - msolve
+        Traceback (most recent call last):
+        ...
+        NotImplementedError: unsupported base field: Rational Field
+    """
+
+    base = ideal.base_ring()
+    if not (isinstance(base, FiniteField) and base.is_prime_field() and
+            base.characteristic() < 2**31):
+        raise NotImplementedError(f"unsupported base field: {base}")
+
+    drlpolring = ideal.ring().change_ring(order='degrevlex')
+    msolve_out = _run_msolve(ideal, ["-g", "2"])
+    gbasis = sage_eval(msolve_out[:-2], locals=drlpolring.gens_dict())
+    return Sequence(gbasis)
+
+def variety(ideal, ring, *, proof=True):
     r"""
     Compute the variety of a zero-dimensional ideal using msolve.
 
@@ -44,7 +123,27 @@ def _variety(ideal, ring, proof):
 
     TESTS::
 
-        sage: K.<x, y> = PolynomialRing(QQ, 2, order='lex')
+        sage: p = 536870909
+        sage: R.<x, y> = PolynomialRing(GF(p), 2, order='lex')
+        sage: I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
+
+        sage: sorted(I.variety(algorithm="msolve", proof=False), key=str) # optional - msolve
+        [{x: 1, y: 1}, {x: 267525699, y: 473946006}]
+
+        sage: K.<a> = GF(p^2)
+        sage: sorted(I.variety(K, algorithm="msolve", proof=False), key=str) # optional - msolve
+        [{x: 1, y: 1},
+         {x: 118750849*a + 194048031, y: 510295713*a + 18174854},
+         {x: 267525699, y: 473946006},
+         {x: 418120060*a + 75297182, y: 26575196*a + 44750050}]
+
+        sage: R.<x, y> = PolynomialRing(GF(2147483659), 2, order='lex')
+        sage: ideal([x, y]).variety(algorithm="msolve", proof=False)
+        Traceback (most recent call last):
+        ...
+        NotImplementedError: unsupported base field: Finite Field of size 2147483659
+
+        sage: R.<x, y> = PolynomialRing(QQ, 2, order='lex')
         sage: I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
 
         sage: I.variety(algorithm='msolve', proof=False) # optional - msolve
@@ -98,67 +197,46 @@ def _variety(ideal, ring, proof):
         ...
         ValueError: positive-dimensional ideal
 
-        sage: K.<x, y> = PolynomialRing(RR, 2, order='lex')
+        sage: R.<x, y> = PolynomialRing(RR, 2, order='lex')
         sage: Ideal(x, y).variety(algorithm='msolve', proof=False)
         Traceback (most recent call last):
         ...
         NotImplementedError: unsupported base field: Real Field with 53 bits of precision
 
-        sage: K.<x, y> = PolynomialRing(QQ, 2, order='lex')
+        sage: R.<x, y> = PolynomialRing(QQ, 2, order='lex')
         sage: Ideal(x, y).variety(ZZ, algorithm='msolve', proof=False)
         Traceback (most recent call last):
         ...
         ValueError: no coercion from base field Rational Field to output ring Integer Ring
     """
 
-    # Normalize and check input
+    proof = sage.structure.proof.proof.get_flag(proof, "polynomial")
+    if proof:
+        raise ValueError("msolve relies on heuristics; please use proof=False")
 
     base = ideal.base_ring()
     if ring is None:
         ring = base
-    proof = sage.structure.proof.proof.get_flag(proof, "polynomial")
-    if proof:
-        raise ValueError("msolve relies on heuristics; please use proof=False")
-    # As of msolve 0.2.4, prime fields seem to be supported, by I cannot
-    # make sense of msolve's output in the positive characteristic case.
-    # if not (base is QQ or isinstance(base, FiniteField) and
-    #         base.is_prime_field() and base.characteristic() < 2**31):
-    if base is not QQ:
-        raise NotImplementedError(f"unsupported base field: {base}")
     if not ring.has_coerce_map_from(base):
         raise ValueError(
             f"no coercion from base field {base} to output ring {ring}")
 
-    # Run msolve
-
-    msolve().require()
-
-    drlpolring = ideal.ring().change_ring(order='degrevlex')
-    polys = ideal.change_ring(drlpolring).gens()
-    msolve_in = tempfile.NamedTemporaryFile(mode='w',
-                                            encoding='ascii', delete=False)
-    command = ["msolve", "-f", msolve_in.name]
     if isinstance(ring, (RealIntervalField_class, RealBallField,
                          RealField_class, RealDoubleField_class)):
         parameterization = False
-        command += ["-p", str(ring.precision())]
+        options = ["-p", str(ring.precision())]
     else:
         parameterization = True
-        command += ["-P", "1"]
-    try:
-        print(",".join(drlpolring.variable_names()), file=msolve_in)
-        print(base.characteristic(), file=msolve_in)
-        print(*(pol._repr_().replace(" ", "") for pol in polys),
-                sep=',\n', file=msolve_in)
-        msolve_in.close()
-        msolve_out = subprocess.run(command, capture_output=True, text=True)
-    finally:
-        os.unlink(msolve_in.name)
-    msolve_out.check_returncode()
+        options = ["-P", "1"]
+
+    msolve_out = _run_msolve(ideal, options)
 
     # Interpret output
 
-    data = sage_eval(msolve_out.stdout[:-2])
+    try:
+        data = sage_eval(msolve_out[:-2])
+    except SyntaxError:
+        raise NotImplementedError(f"unsupported msolve output format: {data}")
 
     dim = data[0]
     if dim == -1:
@@ -172,23 +250,32 @@ def _variety(ideal, ring, proof):
 
     if parameterization:
 
-        def to_poly(p, upol=PolynomialRing(base, 't')):
-            assert len(p[1]) == p[0] + 1
-            return upol(p[1])
+        def to_poly(p, d=1, *, upol=PolynomialRing(base, 't')):
+            assert len(p[1]) == p[0] + 1 or p == [-1, [0]]
+            return upol(p[1])/d
 
-        if len(data) != 3:
+        try:
+            [char, nvars, deg, vars, _, [one, [elim, den, param]]] = data[1]
+        except (IndexError, ValueError):
             raise NotImplementedError(
                 f"unsupported msolve output format: {data}")
-        [dim1, nvars, _, vars, _, [one, elim, den, param]] = data[1]
-        assert dim1.is_zero()
+        assert char == ideal.base_ring().characteristic()
         assert one.is_one()
         assert len(vars) == nvars
         ringvars = out_ring.variable_names()
         assert sorted(vars[:len(ringvars)]) == sorted(ringvars)
         vars = [out_ring(name) for name in vars[:len(ringvars)]]
         elim = to_poly(elim)
+        # Criterion suggested by Mohab Safey El Din to avoid cases where there
+        # is no rational parameterization or where the one returned by msolve
+        # has a significant probability of being incorrect.
+        if deg >= char > 0 or 0 < char <= 2**17 and deg != elim.degree():
+            raise NotImplementedError(f"characteristic {char} too small")
         den = to_poly(den)
-        param = [to_poly(f)/d for [f, d] in param]
+        # As of msolve 0.4.4, param is of the form [pol, denom] in char 0, but
+        # [pol] in char p > 0. My understanding is that both cases will
+        # eventually use the same format, so let's not be too picky.
+        param = [to_poly(*f) for f in param]
         elim_roots = elim.roots(ring, multiplicities=False)
         variety = []
         for rt in elim_roots:
@@ -201,10 +288,9 @@ def to_poly(p, upol=PolynomialRing(base, 't')):
 
     else:
 
-        if len(data) != 2 or data[1][0] != 1:
+        if len(data[1]) < 2 or len(data[1]) != data[1][0] + 1:
             raise NotImplementedError(
                 f"unsupported msolve output format: {data}")
-        _, [_, variety] = data
         if isinstance(ring, (RealIntervalField_class, RealBallField)):
             to_out_ring = ring
         else:
@@ -213,7 +299,8 @@ def to_poly(p, upol=PolynomialRing(base, 't')):
             to_out_ring = lambda iv: ring.coerce(myRIF(iv).center())
         vars = out_ring.gens()
         variety = [[to_out_ring(iv) for iv in point]
-                   for point in variety]
+                   for l in data[1][1:]
+                   for point in l]
 
     return [KeyConvertingDict(out_ring, zip(vars, point)) for point in variety]
 

src/sage/rings/polynomial/multi_polynomial_ideal.py

@@ -2373,9 +2373,6 @@ def variety(self, ring=None, *, algorithm="triangular_decomposition", proof=True
             sage: I.variety(ring=AA)
             [{y: 1, x: 1},
              {y: 0.3611030805286474?, x: 2.769292354238632?}]
-            sage: I.variety(RBF, algorithm='msolve', proof=False) # optional - msolve
-            [{x: [2.76929235423863 +/- 2.08e-15], y: [0.361103080528647 +/- 4.53e-16]},
-             {x: 1.000000000000000, y: 1.000000000000000}]
 
         and a total of four intersections::
 
@@ -2394,6 +2391,13 @@ def variety(self, ring=None, *, algorithm="triangular_decomposition", proof=True
              {y: 0.3611030805286474?, x: 2.769292354238632?},
              {y: 1, x: 1}]
 
+        We can also use the external program msolve to compute the variety.
+        See :mod:`~sage.rings.polynomial.msolve` for more information.
+
+            sage: I.variety(RBF, algorithm='msolve', proof=False) # optional - msolve
+            [{x: [2.76929235423863 +/- 2.08e-15], y: [0.361103080528647 +/- 4.53e-16]},
+             {x: 1.000000000000000, y: 1.000000000000000}]
+
         Computation over floating point numbers may compute only a partial solution,
         or even none at all. Notice that x values are missing from the following variety::
 
@@ -2437,6 +2441,26 @@ def variety(self, ring=None, *, algorithm="triangular_decomposition", proof=True
             sage: v["y"]
             -7.464101615137755?
 
+        msolve also works over finite fields::
+
+            sage: R.<x, y> = PolynomialRing(GF(536870909), 2, order='lex')
+            sage: I = Ideal([ x^2 - 1, y^2 - 1 ])
+            sage: sorted(I.variety(algorithm='msolve', proof=False), key=str) # optional - msolve
+            [{x: 1, y: 1},
+             {x: 1, y: 536870908},
+             {x: 536870908, y: 1},
+             {x: 536870908, y: 536870908}]
+
+        but may fail in small characteristic, especially with ideals of high
+        degree with respect to the characteristic::
+
+            sage: R.<x, y> = PolynomialRing(GF(3), 2, order='lex')
+            sage: I = Ideal([ x^2 - 1, y^2 - 1 ])
+            sage: I.variety(algorithm='msolve', proof=False) # optional - msolve
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: characteristic 3 too small
+
         ALGORITHM:
 
         - With ``algorithm`` = ``"triangular_decomposition"`` (default),
@@ -2453,7 +2477,7 @@ def variety(self, ring=None, *, algorithm="triangular_decomposition", proof=True
             return self._variety_triangular_decomposition(ring)
         elif algorithm == "msolve":
             from . import msolve
-            return msolve._variety(self, ring, proof)
+            return msolve.variety(self, ring, proof=proof)
         else:
             raise ValueError(f"unknown algorithm {algorithm!r}")
 
@@ -3988,6 +4012,10 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
         'macaulay2:mgb'
             Macaulay2's ``GroebnerBasis`` command with the strategy "MGB" (if available)
 
+        'msolve'
+            `msolve <https://msolve.lip6.fr/>`_ (if available, degrevlex order,
+            prime fields)
+
         'magma:GroebnerBasis'
             Magma's ``Groebnerbasis`` command (if available)
 
@@ -4111,6 +4139,15 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
             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]
 
+        Over prime fields of small characteristic, we can also use
+        `msolve <https://msolve.lip6.fr/>`_ (if available in the system program
+        search path)::
+
+            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('msolve')  # optional - msolve
+            [a + 2*b + 2*c - 1, b*c - 19*c^2 + 10*b + 40*c, b^2 - 41*c^2 + 20*b - 20*c, c^3 + 28*c^2 - 37*b + 13*c]
+
         ::
 
             sage: I = sage.rings.ideal.Katsura(P,3) # regenerate to prevent caching
@@ -4329,6 +4366,15 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
             sage: I = sage.rings.ideal.Katsura(P,3) # regenerate to prevent caching
             sage: I.groebner_basis('magma:GroebnerBasis') # optional - magma
             [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]
+
+        msolve currently supports the degrevlex order only::
+
+            sage: R.<a,b,c> = PolynomialRing(GF(101), 3, order='lex')
+            sage: I = sage.rings.ideal.Katsura(R,3) # regenerate to prevent caching
+            sage: I.groebner_basis('msolve')  # optional - msolve
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: msolve only supports the degrevlex order (use transformed_basis())
         """
         from sage.rings.polynomial.multi_polynomial_sequence import PolynomialSequence
         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
@@ -4414,7 +4460,14 @@ def groebner_basis(self, algorithm='', deg_bound=None, mult_bound=None, prot=Fal
         elif algorithm == 'giac:gbasis':
             from sage.libs.giac import groebner_basis as groebner_basis_libgiac
             gb = groebner_basis_libgiac(self, prot=prot, *args, **kwds)
-
+        elif algorithm == 'msolve':
+            if self.ring().term_order() != 'degrevlex':
+                raise NotImplementedError("msolve only supports the degrevlex order "
+                                          "(use transformed_basis())")
+            if not (deg_bound is mult_bound is None) or prot:
+                raise NotImplementedError("unsupported options for msolve")
+            from . import msolve
+            return msolve.groebner_basis_degrevlex(self, *args, **kwds)
         else:
             raise NameError("Algorithm '%s' unknown."%algorithm)