Trac #33734: variety() for polynomial systems over ℚ using msolve

Add the option of using the external package [https://msolve.lip6.fr/
msolve] for computing the variety of a polynomial system. Only systems
over ℚ are supported at the moment, because I cannot make sense of
msolve's output in the positive characteristic case.

Note that this ticket does not add msolve as a dependency—see #31664 for
that.

URL: https://trac.sagemath.org/33734
Reported by: mmezzarobba
Ticket author(s): Marc Mezzarobba
Reviewer(s): Sébastien Labbé, Dima Pasechnik

src/sage/features/msolve.py

@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+r"""
+Feature for testing the presence of msolve
+
+`msolve <https://msolve.lip6.fr/>`_ is a multivariate polynomial system solver
+developed mainly by Jérémy Berthomieu (Sorbonne University), Christian Eder
+(TU Kaiserslautern), and Mohab Safey El Din (Sorbonne University).
+"""
+
+import subprocess
+from . import Executable
+from . import FeatureTestResult
+
+class msolve(Executable):
+    r"""
+    A :class:`~sage.features.Feature` describing the presence of msolve
+
+    EXAMPLES::
+
+        sage: from sage.features.msolve import msolve
+        sage: msolve().is_present() # optional - msolve
+        FeatureTestResult('msolve', True)
+    """
+    def __init__(self):
+        r"""
+        TESTS::
+
+            sage: from sage.features.msolve import msolve
+            sage: isinstance(msolve(), msolve)
+            True
+        """
+        Executable.__init__(self, "msolve", executable="msolve",
+                            url="https://msolve.lip6.fr/")
+
+    def is_functional(self):
+        r"""
+        Test if our installation of msolve is working
+
+        EXAMPLES::
+
+            sage: from sage.features.msolve import msolve
+            sage: msolve().is_functional() # optional - msolve
+            FeatureTestResult('msolve', True)
+        """
+        msolve_out = subprocess.run(["msolve", "-h"], capture_output=True)
+
+        if msolve_out.returncode != 0:
+            return FeatureTestResult(self, False, reason="msolve -h returned "
+                                f"non-zero exit status {msolve_out.returncode}")
+        elif (msolve_out.stdout[:46] !=
+              b'\nmsolve library for polynomial system solving\n'):
+            return FeatureTestResult(self, False,
+                                     reason="output of msolve -h not recognized")
+        return FeatureTestResult(self, True)
+
+def all_features():
+    return [msolve()]

src/sage/rings/polynomial/msolve.py

@@ -0,0 +1,220 @@
+# coding: utf-8
+r"""
+Solution of polynomial systems using msolve
+
+`msolve <https://msolve.lip6.fr/>`_ is a multivariate polynomial system solver
+developed mainly by Jérémy Berthomieu (Sorbonne University), Christian Eder
+(TU Kaiserslautern), and Mohab Safey El Din (Sorbonne University).
+
+This module provide implementations of some operations on polynomial ideals
+based on msolve. Currently the only supported operation is the computation of
+the variety of zero-dimensional ideal over the rationals.
+
+Note that msolve must be installed separately.
+
+.. SEEALSO::
+
+- :mod:`sage.features.msolve`
+- :mod:`sage.rings.polynomial.multi_polynomial_ideal`
+"""
+
+import os
+import tempfile
+import subprocess
+
+import sage.structure.proof.proof
+
+from sage.features.msolve import msolve
+from sage.misc.all import SAGE_TMP
+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.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
+
+
+def _variety(ideal, ring, proof):
+    r"""
+    Compute the variety of a zero-dimensional ideal using msolve.
+
+    Part of the initial implementation was loosely based on the example
+    interfaces available as part of msolve, with the authors' permission.
+
+    TESTS::
+
+        sage: K.<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
+        [{x: 1, y: 1}]
+        sage: I.variety(RealField(100), algorithm='msolve', proof=False) # optional - msolve
+        [{x: 2.7692923542386314152404094643, y: 0.36110308052864737763464656216},
+         {x: 1.0000000000000000000000000000, y: 1.0000000000000000000000000000}]
+        sage: I.variety(RealIntervalField(100), algorithm='msolve', proof=False) # optional - msolve
+        [{x: 2.76929235423863141524040946434?, y: 0.361103080528647377634646562159?},
+         {x: 1, y: 1}]
+        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}]
+        sage: I.variety(RDF, algorithm='msolve', proof=False) # optional - msolve
+        [{x: 2.7692923542386314, y: 0.36110308052864737}, {x: 1.0, y: 1.0}]
+        sage: I.variety(AA, algorithm='msolve', proof=False) # optional - msolve
+        [{x: 2.769292354238632?, y: 0.3611030805286474?},
+         {x: 1.000000000000000?, y: 1.000000000000000?}]
+        sage: I.variety(QQbar, algorithm='msolve', proof=False) # optional - msolve
+        [{x: 2.769292354238632?, y: 0.3611030805286474?},
+         {x: 1, y: 1},
+         {x: 0.11535382288068429? + 0.5897428050222055?*I, y: 0.3194484597356763? - 1.633170240915238?*I},
+         {x: 0.11535382288068429? - 0.5897428050222055?*I, y: 0.3194484597356763? + 1.633170240915238?*I}]
+        sage: I.variety(ComplexField(100))
+        [{y: 1.0000000000000000000000000000, x: 1.0000000000000000000000000000},
+         {y: 0.36110308052864737763464656216, x: 2.7692923542386314152404094643},
+         {y: 0.31944845973567631118267671892 - 1.6331702409152376561188467320*I, x: 0.11535382288068429237979526783 + 0.58974280502220550164728074602*I},
+         {y: 0.31944845973567631118267671892 + 1.6331702409152376561188467320*I, x: 0.11535382288068429237979526783 - 0.58974280502220550164728074602*I}]
+
+        sage: Ideal(x^2 + y^2 - 1, x - y).variety(RBF, algorithm='msolve', proof=False) # optional - msolve
+        [{x: [-0.707106781186547 +/- 6.29e-16], y: [-0.707106781186547 +/- 6.29e-16]},
+         {x: [0.707106781186547 +/- 6.29e-16], y: [0.707106781186547 +/- 6.29e-16]}]
+        sage: sorted(Ideal(x^2 - 1, y^2 - 1).variety(QQ, algorithm='msolve', proof=False), key=str) # optional - msolve
+        [{x: -1, y: -1}, {x: -1, y: 1}, {x: 1, y: -1}, {x: 1, y: 1}]
+        sage: Ideal(x^2-1, y^2-2).variety(CC, algorithm='msolve', proof=False) # optional - msolve
+        [{x: 1.00000000000000, y: 1.41421356237310},
+         {x: -1.00000000000000, y: 1.41421356237309},
+         {x: 1.00000000000000, y: -1.41421356237309},
+         {x: -1.00000000000000, y: -1.41421356237310}]
+
+        sage: Ideal([x, y, x + y]).variety(algorithm='msolve', proof=False) # optional - msolve
+        [{x: 0, y: 0}]
+
+        sage: Ideal([x, y, x + y - 1]).variety(algorithm='msolve', proof=False) # optional - msolve
+        []
+        sage: Ideal([x, y, x + y - 1]).variety(RR, algorithm='msolve', proof=False) # optional - msolve
+        []
+
+        sage: Ideal([x*y - 1]).variety(QQbar, algorithm='msolve', proof=False) # optional - msolve
+        Traceback (most recent call last):
+        ...
+        ValueError: positive-dimensional ideal
+
+        sage: K.<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: 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
+
+    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(dir=SAGE_TMP, 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())]
+    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()
+
+    # Interpret output
+
+    data = sage_eval(msolve_out.stdout[:-2])
+
+    dim = data[0]
+    if dim == -1:
+        return []
+    elif dim > 0:
+        raise ValueError("positive-dimensional ideal")
+    else:
+        assert dim.is_zero()
+
+    out_ring = ideal.ring().change_ring(ring)
+
+    if parameterization:
+
+        def to_poly(p, upol=PolynomialRing(base, 't')):
+            assert len(p[1]) == p[0] + 1
+            return upol(p[1])
+
+        if len(data) != 3:
+            raise NotImplementedError(
+                f"unsupported msolve output format: {data}")
+        [dim1, nvars, _, vars, _, [one, elim, den, param]] = data[1]
+        assert dim1.is_zero()
+        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)
+        den = to_poly(den)
+        param = [to_poly(f)/d for [f, d] in param]
+        elim_roots = elim.roots(ring, multiplicities=False)
+        variety = []
+        for rt in elim_roots:
+            den_of_rt = den(rt)
+            point = [-p(rt)/den_of_rt for p in param]
+            if len(param) != len(vars):
+                point.append(rt)
+            assert len(point) == len(vars)
+            variety.append(point)
+
+    else:
+
+        if len(data) != 2 or 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:
+            assert isinstance(ring, (RealField_class, RealDoubleField_class))
+            myRIF = RealIntervalField(ring.precision())
+            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]
+
+    return [KeyConvertingDict(out_ring, zip(vars, point)) for point in variety]
+

src/sage/rings/polynomial/multi_polynomial_ideal.py

@@ -232,6 +232,7 @@
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
 
+
 from sage.interfaces.singular import singular as singular_default
 from sage.interfaces.magma import magma as magma_default
 
@@ -2283,7 +2284,7 @@ def saturation(self, other):
         return (R.ideal(ideal), ZZ(expo))
 
     @require_field
-    def variety(self, ring=None):
+    def variety(self, ring=None, *, algorithm="triangular_decomposition", proof=True):
         r"""
         Return the variety of this ideal.
 
@@ -2317,6 +2318,8 @@ def variety(self, ring=None):
 
         - ``ring`` - return roots in the ``ring`` instead of the base
           ring of this ideal (default: ``None``)
+        - ``algorithm`` - algorithm or implementation to use; see below for
+          supported values
         - ``proof`` - return a provably correct result (default: ``True``)
 
         EXAMPLES::
@@ -2370,7 +2373,9 @@ def variety(self, ring=None):
             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::
 
@@ -2432,6 +2437,34 @@ def variety(self, ring=None):
             sage: v["y"]
             -7.464101615137755?
 
+        ALGORITHM:
+
+        - With ``algorithm`` = ``"triangular_decomposition"`` (default),
+          uses triangular decomposition, via Singular if possible, falling back
+          on a toy implementation otherwise.
+
+        - With ``algorithm`` = ``"msolve"``, calls the external program
+          `msolve <https://msolve.lip6.fr/>`_ (if available in the system
+          program search path). Note that msolve uses heuristics and therefore
+          requires setting the ``proof`` flag to ``False``. See
+          :mod:`~sage.rings.polynomial.msolve` for more information.
+        """
+        if algorithm == "triangular_decomposition":
+            return self._variety_triangular_decomposition(ring)
+        elif algorithm == "msolve":
+            from . import msolve
+            return msolve._variety(self, ring, proof)
+        else:
+            raise ValueError(f"unknown algorithm {algorithm!r}")
+
+    def _variety_triangular_decomposition(self, ring):
+        r"""
+        Compute the variety of this ideal by triangular decomposition.
+
+        The triangular decomposition is computed using Singular when conversion
+        of the ideal to Singular is supported, falling back to a toy Python
+        implementation otherwise.
+
         TESTS::
 
             sage: K.<w> = GF(27)
@@ -2531,10 +2564,8 @@ def variety(self, ring=None):
             sage: len(I.variety())
             4
 
-        ALGORITHM:
-
-        Uses triangular decomposition.
         """
+
         def _variety(T, V, v=None):
             """
             Return variety ``V`` for one triangular set of
@@ -2571,7 +2602,7 @@ def _variety(T, V, v=None):
             return []
 
         if isinstance(self.base_ring(), sage.rings.abc.ComplexField):
-            verbose("Warning: computations in the complex field are inexact; variety may be computed partially or incorrectly.", level=0)
+            verbose("Warning: computations in the complex field are inexact; variety may be computed partially or incorrectly.", level=0, caller_name="variety")
         P = self.ring()
         if ring is not None:
             P = P.change_ring(ring)
@@ -2580,7 +2611,7 @@ def _variety(T, V, v=None):
             T = [list(each.gens()) for each in TI]
         except TypeError:  # conversion to Singular not supported
             if self.ring().term_order().is_global():
-                verbose("Warning: falling back to very slow toy implementation.", level=0)
+                verbose("Warning: falling back to very slow toy implementation.", level=0, caller_name="variety")
                 T = toy_variety.triangular_factorization(self.groebner_basis())
             else:
                 raise TypeError("Local/unknown orderings not supported by 'toy_buchberger' implementation.")