implement doubly_lexical_ordering in 01-matrix

Author: kappybar

src/doc/en/reference/references/index.rst

@@ -269,6 +269,9 @@ REFERENCES:
              finite Drinfeld modules.* manuscripta mathematica 93, 1 (01 Aug 1997),
              369–379. https://doi.org/10.1007/BF02677478
 
+.. [Anna1987] Lubiw, Anna. Doubly Lexical Orderings of Matrices.
+              SIAM Journal on Computing 16.5 (1987): 854-879.
+
 .. [ANR2023] Robert Angarone, Anastasia Nathanson, and Victor Reiner. *Chow rings of
              matroids as permutation representations*, 2023. :arxiv:`2309.14312`.
 
@@ -3421,6 +3424,10 @@ REFERENCES:
 .. [Hoc] Winfried Hochstaettler, "About the Tic-Tac-Toe Matroid",
          preprint.
 
+.. [Hoffman1985] Hoffman, Alan J., Anthonius Wilhelmus Johannes Kolen, and Michel Sakarovitch.
+                 Totally-balanced and greedy matrices.
+                 SIAM Journal on Algebraic Discrete Methods 6.4 (1985): 721-730.
+
 .. [HJ18] Thorsten Holm and  Peter Jorgensen
     *A p-angulated generalisation of Conway and Coxeter's theorem on frieze patterns*,
     International Mathematics Research Notices (2018)

src/sage/matrix/matrix_mod2_dense.pyx

@@ -2132,6 +2132,169 @@ cdef class Matrix_mod2_dense(matrix_dense.Matrix_dense):   # dense or sparse
         verbose("done computing right kernel matrix over integers mod 2 for %sx%s matrix" % (self.nrows(), self.ncols()),level=1, t=tm)
         return 'computed-pluq', M
 
+    def doubly_lexical_ordering(self, inplace=False):
+        r"""
+        Return a doubly lexical ordering of the matrix.
+
+        A doubly lexical ordering of a matrix is an ordering of the rows
+        and of the columns of the matrix so that both the rows and the
+        columns, as vectors, are lexically increasing. See [Anna1987]_.
+        A lexical ordering of vectors is the standard dictionary ordering,
+        except that vectors will be read from highest to lowest coordinate.
+        Thus row vectors will be compared from right to left, and column
+        vectors from bottom to top.
+
+        INPUT:
+
+        - ``inplace`` -- boolean (default: ``False``); using ``inplace=True``
+          will permute the rows and columns of the current matrix
+          according to a doubly lexical ordering. This will modify the matrix.
+
+        OUTPUT:
+
+        Returns a pair (``row_ordering``, ``col_ordering``). Each item
+        is a ``PermutationGroupElement`` that represents a doubly lexical
+        ordering of the rows or columns.
+
+        .. SEEALSO::
+
+            - :meth:`~sage.matrix.matrix2.Matrix.permutation_normal_form` --
+              a similar matrix normal form
+
+        ALGORITHM:
+
+        The algorithm is adapted from section 3 of [Hoffman1985]_. The time
+        complexity of this algorithm is `O(n \cdot m^2)` for a `n \times m`
+        matrix.
+
+        EXAMPLES:
+
+            sage: A = Matrix(GF(2), [
+            ....:                    [0, 1],
+            ....:                    [1, 0]])
+            sage: r, c = A.doubly_lexical_ordering()
+            sage: r
+            (1,2)
+            sage: c
+            ()
+            sage: A.permute_rows_and_columns(r, c); A
+            [1 0]
+            [0 1]
+
+        ::
+
+            sage: A = Matrix(GF(2), [
+            ....:                    [0, 1],
+            ....:                    [1, 0]])
+            sage: r, c = A.doubly_lexical_ordering(inplace=True); A
+            [1 0]
+            [0 1]
+
+        TESTS:
+
+            sage: A = Matrix(GF(2), [
+            ....:                    [1, 1, 0, 0, 0, 0, 0],
+            ....:                    [1, 1, 0, 0, 0, 0, 0],
+            ....:                    [1, 1, 0, 1, 0, 0, 0],
+            ....:                    [0, 0, 1, 1, 0, 0, 0],
+            ....:                    [0, 1, 1, 1, 1, 0, 0],
+            ....:                    [0, 0, 0, 0, 0, 1, 1],
+            ....:                    [0, 0, 0, 0, 0, 1, 1],
+            ....:                    [0, 0, 0, 0, 1, 1, 1],
+            ....:                    [0, 0, 0, 1, 1, 1, 0]])
+            sage: r, c = A.doubly_lexical_ordering()
+            sage: B = A.with_permuted_rows_and_columns(r, c)
+            sage: flag = True
+            sage: for i in range(B.ncols()):
+            ....:     for j in range(i):
+            ....:         for k in reversed(range(B.nrows())):
+            ....:             if B[k][j] > B[k][i]:
+            ....:                 flag = False
+            ....:                 break
+            ....:             if B[k][j] < B[k][i]:
+            ....:                 break
+            ....:
+            sage: for i in range(B.nrows()):
+            ....:     for j in range(i):
+            ....:         for k in reversed(range(B.ncols())):
+            ....:             if B[j][k] > B[i][k]:
+            ....:                 flag = False
+            ....:                 break
+            ....:             if B[j][k] < B[i][k]:
+            ....:                 break
+            ....:
+            sage: flag
+            True
+            sage: r, c = A.doubly_lexical_ordering(inplace=True)
+            sage: A == B
+            True
+
+        """
+
+        partition_rows = [False for _ in range(self._nrows - 1)]
+        partition_num = 1
+        row_swapped = list(range(1, self._nrows + 1))
+        col_swapped = list(range(1, self._ncols + 1))
+
+        cdef Matrix_mod2_dense A = self if inplace else self.__copy__()
+
+        for i in reversed(range(1, A._ncols + 1)):
+
+            # count 1 for each partition and column
+            count1 = [[0 for _ in range(partition_num)] for _ in range(i)]
+            for col in range(i):
+                parition_i = 0
+                for row in reversed(range(A._nrows)):
+                    count1[col][parition_i] += 1 if mzd_read_bit(A._entries, row, col) else 0
+                    if row > 0 and partition_rows[row - 1]:
+                        parition_i += 1
+
+            # calculate largest_col = col s.t. count1[col] is lexicographically largest (0 <= col < i)
+            largest_col = 0
+            largest_count1 = count1[0]
+            for col in range(1, i):
+                if count1[col] >= largest_count1:
+                    largest_col = col
+                    largest_count1 = count1[col]
+
+            # We refine each partition of rows according to the value of A[:][largest_col].
+            # and also move down rows that satisfy A[row][largest_col] = 1 in each partition.
+            partition_start = 0
+            for _ in range(partition_num):
+                partition_end = partition_start
+                while partition_end < A._nrows - 1 and not partition_rows[partition_end]:
+                    partition_end += 1
+                row_start = partition_start
+                row_end = partition_end
+                while row_start < row_end:
+                    while row_start < row_end and not mzd_read_bit(A._entries, row_start, largest_col):
+                        row_start += 1
+                    while row_start < row_end and mzd_read_bit(A._entries, row_end, largest_col):
+                        row_end -= 1
+                    if row_start < row_end: # swap row
+                        A.swap_rows_c(row_start, row_end)
+                        row_swapped[row_start], row_swapped[row_end] = row_swapped[row_end], row_swapped[row_start]
+                partition_start = partition_end + 1 # for next partition
+
+            for row in range(A._nrows - 1):
+                if mzd_read_bit(A._entries, row, largest_col) != mzd_read_bit(A._entries, row + 1, largest_col):
+                    if not partition_rows[row]:
+                        partition_rows[row] = True
+                        partition_num += 1
+
+            # swap column
+            A.swap_columns_c(largest_col, i - 1)
+            col_swapped[largest_col], col_swapped[i - 1] = col_swapped[i - 1], col_swapped[largest_col]
+
+        from sage.groups.perm_gps.permgroup_named import SymmetricGroup
+        from sage.groups.perm_gps.permgroup_element import make_permgroup_element_v2
+        symmetric_group_nrows = SymmetricGroup(self._nrows)
+        symmetric_group_ncols = SymmetricGroup(self._ncols)
+        row_ordering = make_permgroup_element_v2(symmetric_group_nrows, row_swapped, symmetric_group_nrows.domain())
+        col_ordering = make_permgroup_element_v2(symmetric_group_ncols, col_swapped, symmetric_group_ncols.domain())
+
+        return row_ordering, col_ordering
+
 # Used for hashing
 cdef int i, k
 cdef unsigned long parity_table[256]