add hermitian_decomposition method

src/sage/matrix/matrix2.pyx

@@ -12952,6 +12952,156 @@ cdef class Matrix(Matrix1):
         else:
             return subspace
 
+    def hermitian_decomposition(self):
+        r"""
+        Return a conjugate-transpose factorization of a Hermitian matrix.
+
+        .. NOTE::
+
+            The base ring needs to be a finite field `\GF{q^2}`.
+
+        OUTPUT:
+
+        For a Hermitian matrix `A` the routine returns a matrix `B` such that,
+
+        .. MATH::
+
+            A = B^* B,
+
+        where `B^\ast` is the conjugate-transpose.
+
+        ALGORITHM:
+
+        A conjugate-symmetric version of Gaussian elimination.
+        This is a translation of ``BaseChangeToCanonical`` from
+        the GAP ``forms`` for a Hermitian form.
+
+        TESTS:
+
+        Verify we can decompose a Hermitian matrix::
+
+            sage: A = matrix(GF(11**2),[[1,4,7],[4,1,4],[7,4,1]])
+            sage: B = A.hermitian_decomposition()
+            sage: A == B * B.H
+            True
+
+        Verify we can decompose a Hermitian matrix for which
+        there is no triangular decomposition::
+
+            sage: A = matrix(GF(3**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
+            sage: B = A.hermitian_decomposition()
+            sage: A == B * B.H
+            True
+        """
+        from sage.matrix.constructor import identity_matrix
+        from sage.misc.functional import sqrt
+
+        if not self.is_hermitian():
+            raise ValueError("matrix is not Hermitian")
+
+        F = self._base_ring
+        n = self.nrows()
+
+        if not (F.is_finite() and F.order().is_square()):
+            raise ValueError("the base ring must be a finite field of square order")
+
+        q = sqrt(F.order())
+
+        def conj_square_root(u):
+            if u == 0:
+                return 0
+            z = F.multiplicative_generator()
+            k = u.log(z)
+            if k % (q+1) != 0:
+                raise ValueError(f"unable to factor: {u} is not in base field GF({q})")
+            return z ** (k//(q+1))
+
+        if self.nrows() == 1 and self.ncols() == 1:
+            return self.__class__(F,[conj_square_root(self[0][0])])
+
+        A = copy(self)
+        D = identity_matrix(F, n)
+        row = 0
+
+        # Diagonalize A
+        while True:
+            row += 1
+
+            # Look for a non-zero element on the main diagonal, starting from `row`
+            i = row - 1  # Adjust for zero-based indexing in Sage
+            while i < n and A[i, i].is_zero():
+                i += 1
+
+            if i == row - 1:
+                # Do nothing since A[row, row] != 0
+                pass
+            elif i < n:
+                # Swap to ensure A[row, row] != 0
+                A.swap_rows(row - 1, i)
+                A.swap_columns(row - 1, i)
+                D.swap_rows(row - 1, i)
+            else:
+                # All entries on the main diagonal are zero; look for an off-diagonal element
+                i = row - 1
+                while i < n - 1:
+                    k = i + 1
+                    while k < n and A[i, k].is_zero():
+                        k += 1
+                    if k == n:
+                        i += 1
+                    else:
+                        break
+
+                if i == n - 1:
+                    # All elements are zero; terminate
+                    row -= 1
+                    r = row
+                    break
+
+                # Fetch the non-zero element and place it at A[row, row + 1]
+                if i != row - 1:
+                    A.swap_rows(row - 1, i)
+                    A.swap_columns(row - 1, i)
+                    D.swap_rows(row - 1, i)
+
+                A.swap_rows(row, k)
+                A.swap_columns(row, k)
+                D.swap_rows(row, k)
+
+                b = ~A[row, row-1]
+                A.add_multiple_of_column(row - 1, row, b**q)
+                A.add_multiple_of_row(row - 1, row, b)
+                D.add_multiple_of_row(row - 1, row, b)
+
+            # Eliminate below-diagonal entries in the current column
+            a = ~(-A[row-1, row-1])
+            for i in range(row, n):
+                b = A[i, row - 1] * a
+                if not b.is_zero():
+                    A.add_multiple_of_column(i, row-1, b**q)
+                    A.add_multiple_of_row(i, row - 1, b)
+                    D.add_multiple_of_row(i, row - 1, b)
+
+            if row == n - 1:
+                break
+
+        # Count how many variables have been used
+        if row == n - 1:
+            if A[n-1, n-1]:  # nonzero entry
+                r = n
+            else:
+                r = n - 1
+
+        # Normalize diagonal elements to 1
+        for i in range(r):
+            a = A[i, i]
+            if not a.is_one():
+                # Find an element `b` such that `a = b*b^q = b^(q+1)`
+                b = conj_square_root(a)
+                D.rescale_row(i, 1 / b)
+
+        return D.inverse()
+
     def cholesky(self):
         r"""
         Return the Cholesky decomposition of a Hermitian matrix.