gh-pr-34968: compute matrix kernels modulo composites
    
Fixes #34862.
    
URL: #34968
Reported by: Lorenz Panny
Reviewer(s): Lorenz Panny, Travis Scrimshaw

Author: Release Manager

src/sage/geometry/polyhedron/ppl_lattice_polytope.py

@@ -855,8 +855,8 @@ def base_projection_matrix(self, fiber):
             sage: poly = LatticePolytope_PPL((-9,-6,-1,-1),(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0))
             sage: fiber = next(poly.fibration_generator(2))
             sage: poly.base_projection_matrix(fiber)
-            [0 0 1 0]
-            [0 0 0 1]
+            [ 0  0 -1  0]
+            [ 0  0  0 -1]
 
         Note that the basis choice in :meth:`base_projection` for the
         quotient is usually different::
@@ -866,7 +866,7 @@ def base_projection_matrix(self, fiber):
             sage: [ proj(p) for p in poly.integral_points() ]
             [(-1, -1), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (1, 0)]
             sage: [ proj_matrix*p for p in poly.integral_points() ]
-            [(-1, -1), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (1, 0)]
+            [(1, 1), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, 0), (0, -1), (-1, 0)]
         """
         return matrix(ZZ, fiber.vertices()).right_kernel_matrix()
 

src/sage/matrix/matrix2.pyx

@@ -3670,7 +3670,7 @@ cdef class Matrix(Matrix1):
 
     def _right_kernel_matrix_over_number_field(self):
         r"""
-        Returns a pair that includes a matrix of basis vectors
+        Return a pair that includes a matrix of basis vectors
         for the right kernel of ``self``.
 
         OUTPUT:
@@ -3724,7 +3724,7 @@ cdef class Matrix(Matrix1):
 
     def _right_kernel_matrix_over_field(self, *args, **kwds):
         r"""
-        Returns a pair that includes a matrix of basis vectors
+        Return a pair that includes a matrix of basis vectors
         for the right kernel of ``self``.
 
         OUTPUT:
@@ -3868,6 +3868,38 @@ cdef class Matrix(Matrix1):
                 % (self.nrows(), self.ncols()), level=1, t=tm)
         return 'computed-smith-form', self.new_matrix(nrows=len(basis), ncols=self._ncols, entries=basis)
 
+    def _right_kernel_matrix_over_integer_mod_ring(self):
+        r"""
+        Return a pair that includes a matrix of basis vectors
+        for the right kernel of ``self``.
+
+        OUTPUT:
+
+        Returns a pair. First item is the string 'computed-pari-matkermod'
+        that identifies the nature of the basis vectors.
+
+        Second item is a matrix whose rows are a basis for the right kernel,
+        over the integer modular ring, as computed by :pari:`matkermod`.
+
+        EXAMPLES::
+
+            sage: A = matrix(Zmod(24480), [[1,2,3,4,5],[7,7,7,7,7]])
+            sage: result = A._right_kernel_matrix_over_integer_mod_ring()
+            sage: result[0]
+            'computed-pari-matkermod'
+            sage: P = result[1]; P
+            [    1 24478     1     0     0]
+            [    2 24477     0     1     0]
+            [    3 24476     0     0     1]
+            sage: A*P.transpose() == 0
+            True
+        """
+        R = self.base_ring()
+        pariM = self.change_ring(ZZ).__pari__()
+        basis = list(pariM.matkermod(R.characteristic()))
+        from sage.matrix.matrix_space import MatrixSpace
+        return 'computed-pari-matkermod', MatrixSpace(R, len(basis), ncols=self._ncols)(basis)
+
     def right_kernel_matrix(self, *args, **kwds):
         r"""
         Returns a matrix whose rows form a basis
@@ -3888,10 +3920,11 @@ cdef class Matrix(Matrix1):
             over the rationals and integers
           - 'pluq' - PLUQ matrix factorization for matrices mod 2
 
-        - ``basis`` - default: 'echelon' - a keyword that describes
+        - ``basis`` - default: 'default' - a keyword that describes
           the format of the basis returned.  Allowable values are:
 
-          - 'echelon': the basis matrix is returned in echelon form
+          - 'default': uses 'echelon' over fields; 'computed' otherwise.
+          - 'echelon': the basis matrix is returned in echelon form.
           - 'pivot' : each basis vector is computed from the reduced
             row-echelon form of ``self`` by placing a single one in a
             non-pivot column and zeros in the remaining non-pivot columns.
@@ -4372,10 +4405,6 @@ cdef class Matrix(Matrix1):
             Traceback (most recent call last):
             ...
             ValueError: 'pari' matrix kernel algorithm only available over non-trivial number fields and the integers, not over Rational Field
-            sage: matrix(Integers(6), 2, 2).right_kernel_matrix(algorithm='generic')
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: Echelon form not implemented over 'Ring of integers modulo 6'.
             sage: matrix(QQ, 2, 2).right_kernel_matrix(algorithm='pluq')
             Traceback (most recent call last):
             ...
@@ -4430,13 +4459,15 @@ cdef class Matrix(Matrix1):
         # Determine the basis format of independent spanning set to return
         basis = kwds.pop('basis', None)
         if basis is None:
-            basis = 'echelon'
-        elif basis not in ['computed', 'echelon', 'pivot', 'LLL']:
+            basis = 'default'
+        elif basis not in ['default', 'computed', 'echelon', 'pivot', 'LLL']:
             raise ValueError("matrix kernel basis format '%s' not recognized" % basis )
         elif basis == 'pivot' and R not in _Fields:
             raise ValueError('pivot basis only available over a field, not over %s' % R)
         elif basis == 'LLL' and not is_IntegerRing(R):
             raise ValueError('LLL-reduced basis only available over the integers, not over %s' % R)
+        if basis == 'default':
+            basis = 'echelon' if R in _Fields else 'computed'
 
         # Determine proof keyword for integer matrices
         proof = kwds.pop('proof', None)
@@ -4484,6 +4515,9 @@ cdef class Matrix(Matrix1):
         if M is None and R.is_integral_domain():
             format, M = self._right_kernel_matrix_over_domain()
 
+        if M is None and isinstance(R, sage.rings.abc.IntegerModRing):
+            format, M = self._right_kernel_matrix_over_integer_mod_ring()
+
         if M is None:
             raise NotImplementedError("Cannot compute a matrix kernel over %s" % R)
 
@@ -4523,6 +4557,24 @@ cdef class Matrix(Matrix1):
             else:
                 return M
 
+    def left_kernel_matrix(self, *args, **kwds):
+        r"""
+        Returns a matrix whose rows form a basis for the left kernel
+        of ``self``.
+
+        This method is a thin wrapper around :meth:`right_kernel_matrix`.
+        For supported parameters and input/output formats, see there.
+
+        EXAMPLES::
+
+            sage: M = matrix([[1,2],[3,4],[5,6]])
+            sage: K = M.left_kernel_matrix(); K
+            [ 1 -2  1]
+            sage: K * M
+            [0 0]
+        """
+        return self.transpose().right_kernel_matrix(*args, **kwds)
+
     def right_kernel(self, *args, **kwds):
         r"""
         Returns the right kernel of this matrix, as a vector space or

src/sage/matrix/matrix_modn_dense_template.pxi

@@ -1840,7 +1840,11 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
     def right_kernel_matrix(self, algorithm='linbox', basis='echelon'):
         r"""
         Returns a matrix whose rows form a basis for the right kernel
-        of ``self``, where ``self`` is a matrix over a (small) finite field.
+        of ``self``.
+
+        If the base ring is the ring of integers modulo a composite,
+        the keyword arguments are ignored and the computation is
+        delegated to :meth:`Matrix_dense.right_kernel_matrix`.
 
         INPUT:
 
@@ -1894,7 +1898,10 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
             [0 0 0 0]
         """
         if self.fetch('in_echelon_form') is None:
-            self = self.echelon_form(algorithm=algorithm)
+            try:
+                self = self.echelon_form(algorithm=algorithm)
+            except NotImplementedError:  # composite modulus
+                return Matrix_dense.right_kernel_matrix(self)
 
         cdef Py_ssize_t r = self.rank()
         cdef Py_ssize_t nrows = self._nrows