Laurent polynomials, Fitting ideals and characteristic varieties

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 @@
                         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())