Faster determinant for matrices over gf2e (M4RIE)

src/sage/libs/m4rie.pxd

@@ -5,7 +5,7 @@
 #                  http://www.gnu.org/licenses/
 ##############################################################################
 
-from sage.libs.m4ri cimport mzd_t, m4ri_word
+from sage.libs.m4ri cimport mzd_t, m4ri_word, mzp_t
 
 
 cdef extern from "m4rie/m4rie.h":
@@ -109,6 +109,8 @@ cdef extern from "m4rie/m4rie.h":
 
     size_t mzed_echelonize_ple(mzed_t *, size_t)
 
+    int mzed_ple(mzed_t *A, mzp_t *P, mzp_t *Q)
+
 #cdef extern from "m4rie/strassen.h":
     mzed_t *mzed_mul_strassen(mzed_t *, mzed_t *, mzed_t *, size_t cutoff)
 

src/sage/matrix/matrix_gf2e_dense.pyx

@@ -95,7 +95,7 @@
 from sage.matrix.matrix_mod2_dense cimport Matrix_mod2_dense
 from sage.matrix.args cimport SparseEntry, MatrixArgs_init
 
-from sage.libs.m4ri cimport m4ri_word, mzd_copy
+from sage.libs.m4ri cimport m4ri_word, mzd_copy, mzp_t, mzp_init, mzp_free
 from sage.libs.m4rie cimport *
 from sage.libs.m4rie cimport mzed_t
 
@@ -1010,7 +1010,7 @@
 
             sage: K.<a> = GF(2^3)
             sage: A = random_matrix(K, 3, 3)
             sage: while A.rank() != 3:
             ....:     A = random_matrix(K, 3, 3)
             sage: B = ~A
             sage: A*B
@@ -1586,6 +1586,88 @@
         mzed_cling(self._entries, v)
         mzd_slice_free(v)
 
+    def determinant(self):
+        """
+        Return the determinant of this matrix.
+
+        Relies directly on M4RIE's PLE decomposition, and incidentally caches
+        the rank of ``self``.
+
+        EXAMPLES::
+
+            sage: gf4.<z> = GF(4)
+            sage: mat = matrix(gf4, 2, 2, [[z + 1, z + 1], [z, 1]])
+            sage: mat
+            [z + 1 z + 1]
+            [    z     1]
+            sage: mat.determinant()
+            z
+            sage: gf256.<t> = GF(2**8)
+            sage: mat = matrix(gf256, 3, 3, [[1, t, t**2],
+            ....:                            [t**2, 1, t],
+            ....:                            [t, t**2, 1]])
+            sage: mat.determinant()
+            t^6 + 1
+            sage: mat = matrix(gf256, 3, 3, [[1, t, t**2],
+            ....:                            [t**2, 1, t],
+            ....:                            [t**2 + 1, t + 1, t**2 + t]])
+            sage: mat.determinant()
+            0
+
+        Non-square matrices and the `0 \times 0` matrix are taken care of::
+
+            sage: matrix(gf4, 0, 0).determinant()
+            1
+            sage: matrix(gf4, 3, 2).determinant()
+            Traceback (most recent call last):
+            ...
+            ValueError: self must be a square matrix
+        """
+        cdef size_t m = self._nrows 
+
+        if m != self._ncols:
+            raise ValueError("self must be a square matrix")
+        if m == 0:
+            return self._one
+
+        x = self.fetch('det')
+        if x is not None:
+            return x
+
+        cdef mzed_t * A = mzed_copy(NULL, self._entries)
+        cdef mzp_t * P = mzp_init(m)
+        cdef mzp_t * Q = mzp_init(m)
+
+        sig_on()
+        cdef int r = mzed_ple(A, P, Q)
+        sig_off()
+
+        self.cache('rank', r)
+
+        if r < m:
+            mzp_free(P)
+            mzp_free(Q)
+            mzed_free(A)
+            self.cache('det', self._zero)
+            return self._zero
+
+        cdef Cache_base cache = <Cache_base> self._base_ring._cache
+
+        # characteristic 2, so det(P) == det(Q) == 1
+        cdef Py_ssize_t i
+        cdef int elt
+        cdef det = self._one
+        for i from 0 <= i < m:
+            elt = mzed_read_elem(A, i, i)
+            det = det * cache.fetch_int(elt)
+
+        mzp_free(P)
+        mzp_free(Q)
+        mzed_free(A)
+
+        self.cache('det', det)
+        return det
+
 
 def unpickle_matrix_gf2e_dense_v0(Matrix_mod2_dense a, base_ring, nrows, ncols):
     r"""