gh-36368: Laurent polynomials, Fitting ideals and characteristic vari…

…eties

    
<!-- Describe your changes here in detail -->
 Recently in #36128 (already in develop) characteristic varieties of
finitely presented fundamental groups were introduced. Its computation
is based on Fitting ideals of Laurent polynomial matrices. In #36299,
Fitting ideals were implemented for generic rings with some improvements
for PID and polynomial rings.

There are two original goals in this PR: to improve the output of
characteristic varieties and to use the cited implementation.  In order
to make computations faster, the implementation of Fitting ideals should
apply to Laurent polynomial rings and for this goal, several changes
should be applied to Laurent polynomials in `Sagemath`.

I am not sure if a deeper change should be made, since I applied only
the changes I needed for the above goal:

- src/sage/groups/finitely_presented.py:
    - Changes in `fitting_ideals`.
- src/sage/matrix/matrix2.py: Style changes.
- src/sage/matrix/matrix_laurent_mpolynomial_dense.pyx: This is a new
file to create the class `Matrix_laurent_mpolynomial_dense`.
    - A method `laurent_matrix_reduction` to obtain a matrix of
polynomials where the variables are non common factors for neither the
rows nor the columns.
    - A methord `_fitting_ideal` to use the same method for matrices of
polynomials.
- src/sage/matrix/matrix__mpolynomial_dense.pyx: Style changes.

The main changes are for Laurent polynomials to avoid errors in the
above implementations.

- src/sage/rings/polynomials/laurent_polynomial.pyx:
    - Style changes.
    - Create `xgcd` needed for `inverse_mod`.
    - Create `inverse_mod`.
    - Create `divides`, I copied the code for polynomials with minor
changes.
- src/sage/rings/polynomials/laurent_polynomial_ideal.py:
    - Style changes.
    - Changes in hint keyword in  `__init__`, the previous code create
issues, e.g., impossible to sum ideals of univariate Laurent polynomial
rings. They involve changes in doctests for `hint`
    - Changes in `__contains__` since `__reduce__` is different for
univariate and multivaraite case.
    - Create `gens_reduced`.
    - Changes in `polynomial_ideal` to deal differently if uni- and
multi-variate.
- src/sage/rings/polynomials/laurent_polynomial_mpair.py:
    - Style changes.
    - Create `divides`, I copied the code for polynomials with minor
changes.
- src/sage/rings/polynomials/laurent_polynomial_ring.py:
    - Style changes.
    - Some `TestSuite`'s applied to domains as base_rings; the
corresponding `TestSuite`'s for polynomials also failed if applied to
polynomial rings.
- src/sage/rings/polynomials/laurent_polynomial_ring_base.py:
    - Style changes.
    - Implement `is_noetherian`.
- src/sage/rings/polynomials/polynomial_element.pyx: Style changes.

<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #36368
Reported by: Enrique Manuel Artal Bartolo
Reviewer(s): Enrique Manuel Artal Bartolo, kedlaya, miguelmarco, Travis Scrimshaw

src/sage/groups/finitely_presented.py

@@ -145,6 +145,7 @@
 from sage.sets.set import Set
 from sage.structure.unique_representation import UniqueRepresentation
 
+
 class GroupMorphismWithGensImages(SetMorphism):
     r"""
     Class used for morphisms from finitely presented groups to
@@ -1702,12 +1703,18 @@ def abelian_alexander_matrix(self, ring=QQ, simplified=True):
             sage: g = FreeGroup(1) / []
             sage: g.abelian_alexander_matrix()
             ([], [])
+            sage: g.abelian_alexander_matrix()[0].base_ring()
+            Univariate Laurent Polynomial Ring in f1 over Rational Field
             sage: g = FreeGroup(0) / []
-            sage: g.abelian_alexander_matrix()
-            ([], [])
+            sage: A, ideal = g.abelian_alexander_matrix(); A
+            []
+            sage: A.base_ring()
+            Rational Field
         """
         ab, R, ideal, images = self.abelianization_to_algebra(ring=ring)
         A = self.alexander_matrix(im_gens=images)
+        if A.base_ring() != R:
+            A = A.change_ring(R)
         if simplified:
             n, m = A.dimensions()
             if n == 0 or m == 0:
@@ -1761,66 +1768,111 @@ def characteristic_varieties(self, ring=QQ, matrix_ideal=None, groebner=False):
 
         OUTPUT:
 
-        If ``groebner`` is ``False`` a list of ideals defining the characteristic varieties.
-        If it is ``True``, a list of lists for Gröbner bases for each ideal.
+        A dictionary with keys the indices of the varieties. If ``groebner`` is ``False``
+        the values are the ideals defining the characteristic varieties.
+        If it is ``True``, lists for Gröbner bases for the ideal of each irreducible
+        component, stopping when the first time a characteristic variety is empty.
 
         EXAMPLES::
 
             sage: L = [2*(i, j) + 2* (-i, -j) for i, j in ((1, 2), (2, 3), (3, 1))]
             sage: G = FreeGroup(3) / L
             sage: G.characteristic_varieties(groebner=True)
-            [[(f1 - 1, f2 - 1, f3 - 1),
-             (f1 + 1, f2 - 1, f3 - 1),
-             (f1 - 1, f2 - 1, f3 + 1),
-             (f3^2 + 1, f1 - f3, f2 - f3),
-             (f1 - 1, f2 + 1, f3 - 1)],
-             [(f1 - 1, f2 - 1, f3 - 1),
-             (f1*f3 + 1, f2 - 1),
-             (f1*f2 + 1, f3 - 1),
-             (f2*f3 + 1, f1 - 1),
-             (f2*f3 + 1, f1 - f2),
-             (f2*f3 + 1, f1 - f3),
-             (f1*f3 + 1, f2 - f3)]]
+            {0: [(0,)],
+             1: [(f1 - 1, f2 - 1, f3 - 1), (f1*f3 + 1, f2 - 1), (f1*f2 + 1, f3 - 1), (f2*f3 + 1, f1 - 1),
+                 (f2*f3 + 1, f1 - f2), (f2*f3 + 1, f1 - f3), (f1*f3 + 1, f2 - f3)],
+             2: [(f1 - 1, f2 - 1, f3 - 1), (f1 + 1, f2 - 1, f3 - 1), (f1 - 1, f2 - 1, f3 + 1),
+                 (f3^2 + 1, f1 - f3, f2 - f3), (f1 - 1, f2 + 1, f3 - 1)],
+             3: [(f1 - 1, f2 - 1, f3 - 1)],
+             4: []}
             sage: G = FreeGroup(2)/[2*(1,2,-1,-2)]
             sage: G.characteristic_varieties()
-            [Ideal (-2*f2 + 2, 2*f1 - 2) of Multivariate Laurent Polynomial Ring in f1, f2 over Rational Field]
+            {0: Ideal (0) of Multivariate Laurent Polynomial Ring in f1, f2 over Rational Field,
+             1: Ideal (f2 - 1, f1 - 1) of Multivariate Laurent Polynomial Ring in f1, f2 over Rational Field,
+             2: Ideal (f2 - 1, f1 - 1) of Multivariate Laurent Polynomial Ring in f1, f2 over Rational Field,
+             3: Ideal (1) of Multivariate Laurent Polynomial Ring in f1, f2 over Rational Field}
             sage: G.characteristic_varieties(ring=ZZ)
-            [Ideal (-2*f2 + 2, 2*f1 - 2) of Multivariate Laurent Polynomial Ring in f1, f2 over Integer Ring]
+            {0: Ideal (0) of Multivariate Laurent Polynomial Ring in f1, f2 over Integer Ring,
+             1: Ideal (2*f2 - 2, 2*f1 - 2) of Multivariate Laurent Polynomial Ring in f1, f2 over Integer Ring,
+             2: Ideal (f2 - 1, f1 - 1) of Multivariate Laurent Polynomial Ring in f1, f2 over Integer Ring,
+             3: Ideal (1) of Multivariate Laurent Polynomial Ring in f1, f2 over Integer Ring}
             sage: G = FreeGroup(2)/[(1,2,1,-2,-1,-2)]
             sage: G.characteristic_varieties()
-            [Ideal (1 - f2 + f2^2, -1 + f2 - f2^2) of Univariate Laurent Polynomial Ring in f2 over Rational Field]
+            {0: Ideal (0) of Univariate Laurent Polynomial Ring in f2 over Rational Field,
+             1: Ideal (-1 + 2*f2 - 2*f2^2 + f2^3) of Univariate Laurent Polynomial Ring in f2 over Rational Field,
+             2: Ideal (1) of Univariate Laurent Polynomial Ring in f2 over Rational Field}
+            sage: G.characteristic_varieties(groebner=True)
+            {0: [0], 1: [-1 + f2, 1 - f2 + f2^2], 2: []}
+            sage: G = FreeGroup(2)/[3 * (1, ), 2 * (2, )]
             sage: G.characteristic_varieties(groebner=True)
-            [[1 - f2 + f2^2]]
+            {0: [-1 + F1, 1 + F1, 1 - F1 + F1^2, 1 + F1 + F1^2], 1: [1 - F1 + F1^2],  2: []}
+            sage: G = FreeGroup(2)/[2 * (2, )]
+            sage: G.characteristic_varieties(groebner=True)
+            {0: [(f1 + 1,), (f1 - 1,)], 1: [(f1 + 1,), (f1 - 1, f2 - 1)], 2: []}
+            sage: G = (FreeGroup(0) / [])
+            sage: G.characteristic_varieties()
+            {0: Principal ideal (0) of Rational Field,
+             1: Principal ideal (1) of Rational Field}
+            sage: G.characteristic_varieties(groebner=True)
+            {0: [(0,)], 1: [(1,)]}
         """
-        A, ideal = self.abelian_alexander_matrix(ring=ring, simplified=True)
+        if self.ngens() == 0:
+            if groebner:
+                return {j: [(ring(j),)] for j in (0, 1)}
+            return {j: ring.ideal(j) for j in (0, 1)}
+        A, rels = self.abelian_alexander_matrix(ring=ring, simplified=True)
         R = A.base_ring()
-        res = []
+        eval_1 = {x: ring(1) for x in R.gens()}
+        A_scalar = A.apply_map(lambda p: p.subs(eval_1))
+        n = A.ncols()
+        n1 = n - A_scalar.rank()
+        ideal_1 = R.ideal([x - 1 for x in R.gens()])
         S = R.polynomial_ring()
-        ideal = [S(elt) for elt in ideal]
-        for j in range(1, A.ncols()):
-            L = [p.monomial_reduction()[0] for p in A.minors(j)]
-            J = R.ideal(L + ideal)
-            res.append(J)
-        if not groebner or not R.base_ring().is_field():
+        K = R.base_ring()
+        id_rels = R.ideal(rels)
+        res = dict()
+        bound = n + 1
+        for j in range(bound + 1):
+            J = id_rels + A.fitting_ideal(j)
+            # J = R.ideal(id_rels.gens() + A.fitting_ideal(j).gens())
+            if j <= n1:
+                J1 = K.ideal([K(p.subs(eval_1)) for p in J.gens()])
+                if J1:
+                    J *= ideal_1
+            res[j] = R.ideal(J.gens_reduced())
+            if R(1) in res[j].gens():
+                bound = j
+                break
+        if not groebner or not ring.is_field():
             return res
         if R.ngens() == 1:
-            res0 = [gcd(S(p) for p in J.gens()) for J in res]
-            res1 = []
-            for p in res0:
-                if p == 0:
-                    res1.append([R(0)])
+            res = {j: gcd(S(p) for p in res[j].gens()) for j in range(bound + 1)}
+            char_var = dict()
+            strict = True
+            j = 0
+            while strict and j <= bound:
+                if res[j] == 0:
+                    char_var[j] = [R(0)]
                 else:
-                    fct = [q[0] for q in R(p).factor()]
+                    fct = [q[0] for q in R(res[j]).factor()]
                     if fct:
-                        res1.append(fct)
-            return res1
-        res1 = []
-        for J in res:
-            LJ = J.minimal_associated_primes()
+                        char_var[j] = fct
+                    else:
+                        char_var[j] = []
+                        strict = False
+                j += 1
+            return char_var
+        char_var = dict()
+        strict = True
+        j = 0
+        while strict and j <= bound:
+            LJ = res[j].minimal_associated_primes()
             fct = [id.groebner_basis() for id in LJ]
-            if fct != [(S.one(),)]:
-                res1.append(fct)
-        return res1
+            char_var[j] = fct
+            if not fct:
+                strict = False
+            j += 1
+        return char_var
 
     def rewriting_system(self):
         """

src/sage/matrix/matrix_laurent_mpolynomial_dense.pxd

@@ -0,0 +1,4 @@
+from sage.matrix.matrix_generic_dense cimport Matrix_generic_dense
+
+cdef class Matrix_Laurent_mpolynomial_dense(Matrix_generic_dense):
+    pass

src/sage/matrix/matrix_laurent_mpolynomial_dense.pyx

@@ -0,0 +1,114 @@
+"""
+Dense matrices over multivariate polynomials over fields.
+
+AUTHOR:
+
+- Enrique Artal (2023-??): initial version
+"""
+
+# *****************************************************************************
+#       Copyright (C) 2023 Enrique Artal <artal@unizar.es>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#  as published by the Free Software Foundation; either version 2 of
+#  the License, or (at your option) any later version.
+#                  https://www.gnu.org/licenses/
+# *****************************************************************************
+from sage.matrix.constructor import identity_matrix
+from sage.rings.polynomial.laurent_polynomial_ring_base import LaurentPolynomialRing_generic
+
+cdef class Matrix_laurent_mpolynomial_dense(Matrix_generic_dense):
+    """
+    Dense matrix over a Laurent multivariate polynomial ring over a field.
+    """
+    def laurent_matrix_reduction(self):
+        """
+        From a matrix `self` of Laurent polynomials, apply elementary operations
+        to obtain a matrix `P` of polynomials such that the variables do not divide
+        no column and no row.
+
+        OUTPUT:
+
+        Three matrices `L`, `P`, `R` such that ``self` equals `L P R`, where `L` and
+        `R` are diagonal with monomial entries.
+
+        EXAMPLES:
+
+            sage: R.<x, y> = LaurentPolynomialRing(QQ)
+            sage: L = [1/3*x^-1*y - 6*x^-2*y^2 - 1/2*x^-2*y, 1/5*x + 1/2*y + 1/6]
+            sage: L += [1/2 - 5*x^-1*y - 2*x^-1, -1/3*y^-2 - 4*x^-1*y^-1 + 11*x^-1*y^-2]
+            sage: A = matrix(R, 2, L)
+            sage: lf, P, rg = A.laurent_matrix_reduction()
+            sage: lf
+            [     x^-2         0]
+            [        0 x^-1*y^-2]
+            sage: P
+            [            1/3*x - 6*y - 1/2 1/5*x^3 + 1/2*x^2*y + 1/6*x^2]
+            [        1/2*x*y - 5*y^2 - 2*y             -1/3*x - 4*y + 11]
+            sage: rg
+            [y 0]
+            [0 1]
+        """
+        R = self.base_ring()
+        n_rows, n_cols = self.dimensions()
+        mat_l = identity_matrix(R, n_rows)
+        mat_r = identity_matrix(R, n_cols)
+        res = self.__copy__()
+        for j, rw in enumerate(res.rows()):
+            for t in R.gens():
+                n = min(mon.degree(t) for a in rw for cf, mon in a)
+                res.rescale_row(j, t ** -n)
+                mat_l.rescale_col(j, t ** n)
+        for j, cl in enumerate(res.columns()):
+            for t in R.gens():
+                n = min(mon.degree(t) for a in cl for cf, mon in a)
+                res.rescale_col(j, t ** -n)
+                mat_r.rescale_row(j, t ** n)
+        res = res.change_ring(R.polynomial_ring())
+        return mat_l, res, mat_r
+
+    def _fitting_ideal(self, i):
+        r"""
+        Return the `i`-th Fitting ideal of the matrix. This is the ideal generated
+        by the `n - i` minors, where `n` is the number of columns.
+
+        INPUT:
+
+        ``i`` -- an integer
+
+        OUTPUT:
+
+        An ideal on the base ring.
+
+        EXAMPLES::
+
+            sage: R.<x,y,z> = LaurentPolynomialRing(QQ)
+            sage: M = matrix(R, [[2*x^-1-z, 0, y-z^-2, 0], [0, z - y^-1, z - x, 0],[z - y, x^-2 - y, 0, z]])
+            sage: M
+            [-z + 2*x^-1           0    y - z^-2           0]
+            [          0    z - y^-1      -x + z           0]
+            [     -y + z   -y + x^-2           0           z]
+            sage: M.fitting_ideal(0)
+            Ideal (0) of Multivariate Laurent Polynomial Ring in x, y, z over Rational Field
+            sage: M.fitting_ideal(1) == M._fitting_ideal(1)
+            True
+            sage: M.fitting_ideal(1).groebner_basis()
+            (x^4 - 2*x^3*y - x*z^3 - 4*x^2*y + 8*x*y^2 + 4*x*y*z + 2*z^2 - 8*y,
+            x*y*z^2 - x*z - 2*y*z + 2,
+            x^2*z - x*z^2 - 2*x + 2*z,
+            y^2*z + 1/4*x^2 - 1/2*x*y - 1/4*x*z - y + 1/2)
+            sage: M.fitting_ideal(2).groebner_basis()
+            (1,)
+            sage: M.fitting_ideal(3).groebner_basis()
+            (1,)
+            sage: M.fitting_ideal(4).groebner_basis()
+            (1,)
+            sage: [R.ideal(M.minors(i)) == M._fitting_ideal(4 - i) for i in range(5)]
+            [True, True, True, True, True]
+
+        """
+        R = self.base_ring()
+        S = R.polynomial_ring()
+        A = self.laurent_matrix_reduction()[1].change_ring(S)
+        J = A._fitting_ideal(i)
+        return J.change_ring(R)

src/sage/matrix/matrix_space.py

@@ -306,6 +306,13 @@ def get_matrix_class(R, nrows, ncols, sparse, implementation):
                         pass
                     else:
                         return matrix_mpolynomial_dense.Matrix_mpolynomial_dense
+                elif isinstance(R, sage.rings.polynomial.laurent_polynomial_ring.LaurentPolynomialRing_mpair) and R.base_ring() in _Fields:
+                    try:
+                        from . import matrix_laurent_mpolynomial_dense
+                    except ImportError:
+                        pass
+                    else:
+                        return matrix_laurent_mpolynomial_dense.Matrix_laurent_mpolynomial_dense
 
             # The fallback
             from sage.matrix.matrix_generic_dense import Matrix_generic_dense

src/sage/rings/polynomial/ideal.py

@@ -85,3 +85,16 @@ def groebner_basis(self, algorithm=None):
         gb = self.gens_reduced()
         from sage.rings.polynomial.multi_polynomial_sequence import PolynomialSequence_generic
         return PolynomialSequence_generic([gb], self.ring(), immutable=True)
+
+    def change_ring(self, R):
+        """
+        Coerce an ideal into a new ring.
+
+        EXAMPLES::
+
+            sage: R.<q> = QQ[]
+            sage: I = R.ideal([q^2+q-1])
+            sage: I.change_ring(RR['q'])
+            Principal ideal (q^2 + q - 1.00000000000000) of Univariate Polynomial Ring in q over Real Field with 53 bits of precision
+        """
+        return R.ideal(self.gens())