Trac #30892: Characteristic polynomials over complete discrete valuation rings/fields

We implement Hessenberg algorithm (with choice of pivot) for computing
the characteristic polynomial of a matrix with coefficients in a
complete discrete valuation ring/field.

URL: https://trac.sagemath.org/30892
Reported by: caruso
Ticket author(s): Xavier Caruso, Raphaël Pagès
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/categories/discrete_valuation.py

@@ -64,6 +64,48 @@ def residue_field(self):
                 Rational Field
             """
 
+        def _matrix_charpoly(self, M, var):
+            r"""
+            Return the characteristic polynomial of `M`.
+
+            EXAMPLES::
+
+                sage: R.<t> = PowerSeriesRing(GF(5))
+                sage: M = matrix(4, 4, [ (t^(i+j)).add_bigoh(10)
+                ....:                    for i in range(4) for j in range(4) ])
+                sage: M
+                [  1 + O(t^10)   t + O(t^10) t^2 + O(t^10) t^3 + O(t^10)]
+                [  t + O(t^10) t^2 + O(t^10) t^3 + O(t^10) t^4 + O(t^10)]
+                [t^2 + O(t^10) t^3 + O(t^10) t^4 + O(t^10) t^5 + O(t^10)]
+                [t^3 + O(t^10) t^4 + O(t^10) t^5 + O(t^10) t^6 + O(t^10)]
+                sage: M.charpoly()   # indirect doctest
+                x^4 + (4 + 4*t^2 + 4*t^4 + 4*t^6 + O(t^10))*x^3
+
+            Note that this function uses a Hessenberg-like algorithm
+            that performs divisions. Hence, truncations may show up
+            even if the input matrix is exact::
+
+                sage: M = matrix(3, 3, [ 1, t, t^2, 1+t, t^2, t^3, t^2, t^3, t^4 ])
+                sage: M
+                [    1     t   t^2]
+                [1 + t   t^2   t^3]
+                [  t^2   t^3   t^4]
+                sage: M.charpoly()
+                x^3 + (4 + 4*t^2 + 4*t^4 + O(t^25))*x^2 + (4*t + O(t^24))*x
+
+            Another example over the p-adics::
+
+                sage: R = Zp(5, print_mode="digits", prec=5)
+                sage: M = matrix(R, 3, 3, range(9))
+                sage: M
+                [        0  ...00001  ...00002]
+                [ ...00003  ...00004 ...000010]
+                [ ...00011  ...00012  ...00013]
+                sage: M.charpoly()
+                ...00001*x^3 + ...44423*x^2 + ...44412*x + ...00000
+            """
+            return M._charpoly_hessenberg(var)
+
     class ElementMethods:
         @abstract_method
         def valuation(self):
@@ -200,7 +242,7 @@ def uniformizer(self):
         @abstract_method
         def residue_field(self):
             """
-            Return the residue field of the ring of integers of 
+            Return the residue field of the ring of integers of
             this discrete valuation field.
 
             EXAMPLES::
@@ -213,6 +255,45 @@ def residue_field(self):
                 Rational Field
             """
 
+        def _matrix_hessenbergize(self, H):
+            r"""
+            Replace `H` with an Hessenberg form of it.
+
+            EXAMPLES::
+
+                sage: R.<t> = PowerSeriesRing(GF(5))
+                sage: K = R.fraction_field()
+                sage: H = matrix(K, 4, 4, [ (t^(i+j)).add_bigoh(10)
+                ....:                       for i in range(4) for j in range(4) ])
+                sage: H
+                [  1 + O(t^10)   t + O(t^10) t^2 + O(t^10) t^3 + O(t^10)]
+                [  t + O(t^10) t^2 + O(t^10) t^3 + O(t^10) t^4 + O(t^10)]
+                [t^2 + O(t^10) t^3 + O(t^10) t^4 + O(t^10) t^5 + O(t^10)]
+                [t^3 + O(t^10) t^4 + O(t^10) t^5 + O(t^10) t^6 + O(t^10)]
+                sage: H.hessenbergize()
+                sage: H
+                [              1 + O(t^10)   t + t^3 + t^5 + O(t^10)             t^2 + O(t^10)             t^3 + O(t^10)]
+                [              t + O(t^10) t^2 + t^4 + t^6 + O(t^10)             t^3 + O(t^10)             t^4 + O(t^10)]
+                [                  O(t^10)                   O(t^10)                   O(t^10)                   O(t^10)]
+                [                  O(t^10)                   O(t^10)                   O(t^10)                   O(t^10)]
+
+            Another example over the p-adics::
+
+                sage: K = Qp(5, print_mode="digits", prec=5)
+                sage: H = matrix(K, 3, 3, range(9))
+                sage: H
+                [        0  ...00001  ...00002]
+                [ ...00003  ...00004 ...000010]
+                [ ...00011  ...00012  ...00013]
+                sage: H.hessenbergize()
+                sage: H
+                [        0  ...00010  ...00002]
+                [ ...00003  ...00024 ...000010]
+                [ ...00000  ...44440  ...44443]
+            """
+            from sage.matrix.matrix_cdv import hessenbergize_cdvf
+            hessenbergize_cdvf(H)
+
     class ElementMethods:
         @abstract_method
         def valuation(self):

src/sage/matrix/matrix2.pyx

@@ -2700,6 +2700,8 @@ cdef class Matrix(Matrix1):
 
         ALGORITHM:
 
+        If the base ring has a method `_matrix_charpoly`, we use it.
+
         In the generic case of matrices over a ring (commutative and with
         unity), there is a division-free algorithm, which can be accessed
         using ``"df"``, with complexity `O(n^4)`.  Alternatively, by
@@ -2862,18 +2864,16 @@ cdef class Matrix(Matrix1):
         if f is not None:
             return f.change_variable_name(var)
 
-        if algorithm is None:
-            from sage.rings.finite_rings.integer_mod_ring import is_IntegerModRing
+        R = self._base_ring
 
-            R = self._base_ring
-            if is_NumberField(R):
-                f = self._charpoly_over_number_field(var)
-            elif is_IntegerModRing(R):
-                f = self.lift().charpoly(var).change_ring(R)
-            elif R in _Fields and R.is_exact():
-                f = self._charpoly_hessenberg(var)
-            else:
-                f = self._charpoly_df(var)
+        if algorithm is None:
+            if hasattr(R, '_matrix_charpoly'):
+                f = R._matrix_charpoly(self, var)
+            if f is None:
+                if R in _Fields and R.is_exact():
+                    f = self._charpoly_hessenberg(var)
+                else:
+                    f = self._charpoly_df(var)
         else:
             if algorithm == "hessenberg":
                 f = self._charpoly_hessenberg(var)
@@ -3049,53 +3049,6 @@ cdef class Matrix(Matrix1):
 
         return f
 
-    def _charpoly_over_number_field(self, var='x'):
-        r"""
-        Use PARI to compute the characteristic polynomial of self as a
-        polynomial over the base ring.
-
-        EXAMPLES::
-
-            sage: x = QQ['x'].gen()
-            sage: K.<a> = NumberField(x^2 - 2)
-            sage: m = matrix(K, [[a-1, 2], [a, a+1]])
-            sage: m._charpoly_over_number_field('Z')
-            Z^2 - 2*a*Z - 2*a + 1
-            sage: m._charpoly_over_number_field('a')(m) == 0
-            True
-            sage: m = matrix(K, [[0, a, 0], [-a, 0, 0], [0, 0, 0]])
-            sage: m._charpoly_over_number_field('Z')
-            Z^3 + 2*Z
-
-        The remaining tests are indirect::
-
-            sage: L.<b> = K.extension(x^3 - a)
-            sage: m = matrix(L, [[b+a, 1], [a, b^2-2]])
-            sage: m.charpoly('Z')
-            Z^2 + (-b^2 - b - a + 2)*Z + a*b^2 - 2*b - 2*a
-            sage: m.charpoly('a')
-            a^2 + (-b^2 - b - a + 2)*a + a*b^2 - 2*b - 2*a
-            sage: m.charpoly('a')(m) == 0
-            True
-
-        ::
-
-            sage: M.<c> = L.extension(x^2 - a*x + b)
-            sage: m = matrix(M, [[a+b+c, 0, b], [0, c, 1], [a-1, b^2+1, 2]])
-            sage: f = m.charpoly('Z'); f
-            Z^3 + (-2*c - b - a - 2)*Z^2 + ((b + 2*a + 4)*c - b^2 + (-a + 2)*b + 2*a - 1)*Z + (b^2 + (a - 3)*b - 4*a + 1)*c + a*b^2 + 3*b + 2*a
-            sage: f(m) == 0
-            True
-            sage: f.base_ring() is M
-            True
-        """
-        K = self.base_ring()
-        if not is_NumberField(K):
-            raise ValueError("_charpoly_over_number_field called with base ring (%s) not a number field" % K)
-
-        paripoly = self.__pari__().charpoly()
-        return K[var](paripoly)
-
     def fcp(self, var='x'):
         """
         Return the factorization of the characteristic polynomial of self.
@@ -3373,11 +3326,16 @@ cdef class Matrix(Matrix1):
         if not self.is_square():
             raise TypeError("self must be square")
 
+        self.check_mutability()
+
+        base = self._base_ring
+        if hasattr(base, '_matrix_hessenbergize'):
+            base._matrix_hessenbergize(self)
+            return
+
         if self._base_ring not in _Fields:
             raise TypeError("Hessenbergize only possible for matrices over a field")
 
-        self.check_mutability()
-
         zero = self._base_ring(0)
         one = self._base_ring(1)
         for m from 1 <= m < n-1:

src/sage/matrix/matrix_cdv.pxd

@@ -0,0 +1,3 @@
+from sage.matrix.matrix_generic_dense cimport Matrix_generic_dense
+
+cpdef hessenbergize_cdvf(Matrix_generic_dense)

src/sage/matrix/matrix_cdv.pyx

@@ -0,0 +1,83 @@
+r"""
+Special methods for matrices over discrete valuation rings/fields.
+"""
+
+# ****************************************************************************
+#       Copyright (C) 2020 Xavier Caruso <[email protected]>
+#                          Raphaël Pagès <[email protected]>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License 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.matrix_generic_dense cimport Matrix_generic_dense
+from sage.structure.element cimport RingElement
+from sage.rings.infinity import Infinity
+
+
+# We assume that H is square
+cpdef hessenbergize_cdvf(Matrix_generic_dense H):
+    r"""
+    Replace `H` with an Hessenberg form of it.
+
+    .. NOTE::
+
+        This function assumes that H is a matrix over
+        a complete discrete valuation field.
+
+        The pivot on each column is always chosen
+        with maximal relative precision, which ensures 
+        the numerical stability of the algorithm.
+
+    TESTS::
+
+        sage: K = Qp(5, print_mode="digits", prec=5)
+        sage: H = matrix(K, 3, 3, range(9))
+        sage: H
+        [        0  ...00001  ...00002]
+        [ ...00003  ...00004 ...000010]
+        [ ...00011  ...00012  ...00013]
+        sage: H.hessenbergize()
+        sage: H
+        [        0  ...00010  ...00002]
+        [ ...00003  ...00024 ...000010]
+        [ ...00000  ...44440  ...44443]
+
+    ::
+
+        sage: M = random_matrix(K, 6, 6)
+        sage: M.charpoly()[0] == M.determinant()
+        True
+    """
+    cdef Py_ssize_t n, m, i, j, k
+    cdef Matrix_generic_dense c
+    cdef RingElement pivot, inv, scalar
+
+    n = H.nrows()
+    for j in range(n-1):
+        k = j + 1
+        maxi = H.get_unsafe(k, j).precision_relative()
+        i = j + 2
+        while maxi is not Infinity and i < n:
+            entry = H.get_unsafe(i, j)
+            if entry:
+                m = entry.precision_relative()
+                if m > maxi:
+                    maxi = m
+                    k = i
+            i += 1
+
+        if k != j + 1:
+            H.swap_rows_c(j+1, k)
+            H.swap_columns_c(j+1, k)
+        pivot = H.get_unsafe(j+1, j)
+        if pivot:
+            inv = ~pivot
+            for i in range(j+2, n):
+                scalar = inv * H.get_unsafe(i, j)
+                H.add_multiple_of_row_c(i, j+1, -scalar, j)
+                H.add_multiple_of_column_c(j+1, i, scalar, 0)

src/sage/matrix/matrix_cyclo_dense.pyx

@@ -1285,7 +1285,8 @@ cdef class Matrix_cyclo_dense(Matrix_dense):
         if algorithm == 'multimodular':
             f = self._charpoly_multimodular(var, proof=proof)
         elif algorithm == 'pari':
-            f = self._charpoly_over_number_field(var)
+            paripoly = self.__pari__().charpoly()
+            f = self.base_ring()[var](paripoly)
         elif algorithm == 'hessenberg':
             f = self._charpoly_hessenberg(var)
         else:

src/sage/rings/laurent_series_ring_element.pyx

@@ -155,15 +155,15 @@ cdef class LaurentSeries(AlgebraElement):
 
         # self is that t^n * u:
         if not f:
-            if n == infinity:
+            if n is infinity:
                 self.__n = 0
                 self.__u = parent._power_series_ring.zero()
             else:
                 self.__n = n
                 self.__u = f
         else:
             val = f.valuation()
-            if val == infinity:
+            if val is infinity:
                 self.__n = 0
                 self.__u = f
             elif val == 0:
@@ -328,7 +328,7 @@ cdef class LaurentSeries(AlgebraElement):
             '2 + 2/3*t^3'
         """
         if self.is_zero():
-            if self.prec() == infinity:
+            if self.prec() is infinity:
                 return "0"
             else:
                 return "O(%s^%s)"%(self._parent.variable_name(),self.prec())
@@ -451,7 +451,7 @@ cdef class LaurentSeries(AlgebraElement):
             \left(a + b\right)x
         """
         if self.is_zero():
-            if self.prec() == infinity:
+            if self.prec() is infinity:
                 return "0"
             else:
                 return "0 + \\cdots"
@@ -835,7 +835,7 @@ cdef class LaurentSeries(AlgebraElement):
             sage: (t^(-2)).add_bigoh(-3)
             O(t^-3)
         """
-        if prec == infinity or prec >= self.prec():
+        if prec is infinity or prec >= self.prec():
             return self
         P = self._parent
         if not self or prec < self.__n:
@@ -1691,7 +1691,7 @@ cdef class LaurentSeries(AlgebraElement):
         """
         if prec is None:
             prec = self.prec()
-            if prec == infinity:
+            if prec is infinity:
                 prec = self.parent().default_prec()
         else:
             prec = min(self.prec(), prec)