extend Cholesky decomposition method for finite fields

src/sage/matrix/matrix2.pyx

@@ -12943,7 +12943,179 @@ cdef class Matrix(Matrix1):
         else:
             return subspace
 
-    def cholesky(self):
+    def _cholesky_extended_ff(self):
+        r"""
+        Performs the extended Cholesky decomposition of a Hermitian matrix over a finite field of square order.
+
+        INPUT:
+
+        - ``self`` -- a square matrix with entries from a finite field of square order
+
+        OUTPUT:
+
+        A square matrix over the same finite such that multiplying itself by its conjugate-transpose yields the input matrix
+
+        ALGORITHM:
+
+        First, we ensure the matrix is square and defined over a finite field of square order. Then we can perform the conjugate-symmetric
+        version of Gaussian elimination, but the resulting decomposition matrix `L` might not be lower triangular.
+
+        This is a translation of ``BaseChangeToCanonical`` from the GAP ``forms`` package (for a Hermitian form).
+
+        EXAMPLES:
+
+        Here we use the extended decomposition, where the result may not be a lower triangular matrix::
+
+            sage: U = matrix(GF(17**2),[[0,1],[1,0]])
+            sage: B = U._cholesky_extended_ff(); B
+            [13*z2 + 6 3*z2 + 16]
+            [13*z2 + 6 14*z2 + 1]
+            sage: U == B * B.H
+            True
+            sage: U = matrix(GF(13**2),[[1,4,7],[4,1,4],[7,4,1]])
+            sage: B = U._cholesky_extended_ff(); B
+            [        1         0         0]
+            [        4  7*z2 + 3         0]
+            [        7 6*z2 + 10 12*z2 + 6]
+            sage: U == B * B.H
+            True
+            sage: U = matrix(GF(7**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
+            sage: B = U._cholesky_extended_ff(); B
+            [4*z2 + 2     6*z2        0]
+            [4*z2 + 2       z2        0]
+            [5*z2 + 6       z2       z2]
+            sage: U == B * B.H
+            True
+
+        TESTS:
+
+        If the matrix is not full rank, we compute the rank and throw an exception.
+
+            sage: U = matrix(GF(3**2),[[1,4,7],[4,1,4],[7,4,1]])
+            sage: U._cholesky_extended_ff()
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix is not full rank
+            sage: U = matrix(GF(3**2),[[0,4,7],[4,1,4],[7,4,1]])
+            sage: U._cholesky_extended_ff()
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix is not full rank
+        """
+        from sage.matrix.constructor import identity_matrix
+
+        if not self.is_hermitian():
+            raise ValueError("matrix is not Hermitian")
+
+        F = self._base_ring
+        n = self.nrows()
+        if self.fetch("rank") not in [n, None]:
+            raise ValueError("matrix is not full rank")
+        if not (F.is_finite() and F.order().is_square()):
+            raise ValueError("the base ring must be a finite field of square order")
+
+        q = F.order().sqrt(extend=False)
+
+        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 = -1
+
+        # Diagonalize A
+        while True:
+            row += 1
+
+            # Look for a non-zero element on the main diagonal, starting from `row`
+            i = row
+            while i < n and A[i, i].is_zero():
+                i += 1
+
+            if i == row:
+                # Do nothing since A[row, row] != 0
+                pass
+            elif i < n:
+                # Swap to ensure A[row, row] != 0
+                A.swap_rows(row, i)
+                A.swap_columns(row, i)
+                D.swap_rows(row, i)
+            else:
+                # All entries on the main diagonal are zero; look for an off-diagonal element
+                i = row
+                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 + 1
+                    break
+
+                # Fetch the non-zero element and place it at A[row, row + 1]
+                if i != row:
+                    A.swap_rows(row, i)
+                    A.swap_columns(row, i)
+                    D.swap_rows(row, i)
+
+                A.swap_rows(row + 1, k)
+                A.swap_columns(row + 1, k)
+                D.swap_rows(row + 1, k)
+
+                b = ~A[row + 1, row]
+                A.add_multiple_of_column(row, row + 1, b**q)
+                A.add_multiple_of_row(row, row + 1, b)
+                D.add_multiple_of_row(row, row + 1, b)
+
+            # Eliminate below-diagonal entries in the current column
+            a = ~(-A[row, row])
+            for i in range(row + 1, n):
+                b = A[i, row] * a
+                if not b.is_zero():
+                    A.add_multiple_of_column(i, row, b**q)
+                    A.add_multiple_of_row(i, row, b)
+                    D.add_multiple_of_row(i, row, 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
+
+        if r < n:
+            self.cache('rank', r)
+            raise ValueError("matrix is not full rank")
+
+        # 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, extended=False):
         r"""
         Return the Cholesky decomposition of a Hermitian matrix.
 
@@ -12989,6 +13161,15 @@ cdef class Matrix(Matrix1):
         closure or the algebraic reals, depending on whether or not
         imaginary numbers are required.
 
+        Over finite fields, the Cholesky decomposition might
+        not exist, but when the field has square order (i.e.,
+        `\GF{q^2}`), then we can perform the conjugate-symmetric
+        version of Gaussian elimination, but the resulting
+        decomposition matrix `L` might not be lower triangular.
+
+        This is a translation of ``BaseChangeToCanonical`` from
+        the GAP ``forms`` package (for a Hermitian form).
+
         EXAMPLES:
 
         This simple example has a result with entries that remain
@@ -13159,6 +13340,30 @@ cdef class Matrix(Matrix1):
             ...
             ValueError: matrix is not positive definite
 
+        Here we use the extended decomposition, where the result
+        may not be a lower triangular matrix::
+
+            sage: U = matrix(GF(5**2),[[0,1],[1,0]])
+            sage: B = U.cholesky(extended=True); B
+            [3*z2 4*z2]
+            [3*z2   z2]
+            sage: U == B * B.H
+            True
+            sage: U = matrix(GF(11**2),[[1,4,7],[4,1,4],[7,4,1]])
+            sage: B = U.cholesky(extended=True); B
+            [        1         0         0]
+            [        4  9*z2 + 2         0]
+            [        7 10*z2 + 1  3*z2 + 3]
+            sage: U == B * B.H
+            True
+            sage: U = matrix(GF(3**2), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
+            sage: B = U.cholesky(extended=True); B
+            [2*z2    2    0]
+            [2*z2    1    0]
+            [   0    1   z2]
+            sage: U == B * B.H
+            True
+
         TESTS:
 
         This verifies that :issue:`11274` is resolved::
@@ -13206,6 +13411,9 @@ cdef class Matrix(Matrix1):
             ....:      for R in (RR,CC,RDF,CDF,ZZ,QQ,AA,QQbar) )
             True
         """
+        if extended:
+            return self._cholesky_extended_ff()
+
         cdef Matrix C  # output matrix
         C = self.fetch('cholesky')
         if C is not None: