Merge pull request #261 from oscarbenjamin/pr_fflu

fmpz_mat: add fraction-free LU decomposition

src/flint/flintlib/types/flint.pxd

@@ -55,6 +55,9 @@ cdef extern from *:
     """
 
 cdef extern from "flint/flint.h":
+    # These defines are needed to work around a Cython bug.
+    # Otherwise sizeof(ulong) will give the wrong size on 64 bit Windows.
+    # https://github.com/cython/cython/issues/6339
     """
     #define SIZEOF_ULONG sizeof(ulong)
     #define SIZEOF_SLONG sizeof(slong)

src/flint/test/test_all.py

@@ -1,9 +1,8 @@
-import sys
 import math
 import operator
 import pickle
-import doctest
 import platform
+import random
 
 from flint.utils.flint_exceptions import DomainError, IncompatibleContextError
 
@@ -1871,7 +1870,6 @@ def test_fmpz_mod_dlog():
     p = 2**e2 * 3**e3 + 1
     F = fmpz_mod_ctx(p)
 
-    import random
     for _ in range(10):
         g = F(random.randint(0,p))
         for _ in range(10):
@@ -4084,6 +4082,50 @@ def test_matrices_div():
         raises(lambda: None / M1234, TypeError)
 
 
+def test_matrices_properties():
+    for M, S, is_field in _all_matrices():
+        # XXX: Add these properties to all matrix types
+        if M is not flint.fmpz_mat:
+            continue
+
+        assert M([[1, 2], [3, 4]]).is_square() is True
+        assert M([[1, 2]]).is_square() is False
+        assert M(0, 2, []).is_square() is False
+        assert M(2, 0, []).is_square() is False
+
+        assert M([[1, 2], [3, 4]]).is_empty() is False
+        assert M(0, 2, []).is_empty() is True
+        assert M(2, 0, []).is_empty() is True
+
+        assert M([[1, 2], [3, 4]]).is_zero() is False
+        assert M([[0, 0], [0, 0]]).is_zero() is True
+
+        assert M([[1, 0], [0, 1]]).is_one() is True
+        assert M([[1, 2], [3, 4]]).is_one() is False
+        assert M([[1, 0], [0, 1], [0, 0]]).is_one() is True # ??
+        assert M(0, 0, []).is_one() is True
+
+        assert M([[-1, 0], [0, -1]]).is_neg_one() is True
+        assert M([[-1, 0], [0, 1]]).is_neg_one() is False
+        assert M([[-1, -1], [-1, -1]]).is_neg_one() is False
+        assert M([[-1, 0], [0, -1], [0, 0]]).is_neg_one() is False # ??
+        assert M(0, 0, []).is_neg_one() is True
+
+        assert M([[2, 0], [0, 2]]).is_scalar() is True
+        assert M([[2, 0], [0, 3]]).is_scalar() is False
+        assert M([[1, 0], [0, 1], [0, 0]]).is_scalar() is False
+
+        assert M([[1, 0], [0, 2]]).is_diagonal() is True
+        assert M([[1, 0], [1, 2]]).is_diagonal() is False
+        assert M([[1, 0], [0, 1], [0, 0]]).is_diagonal() is True
+
+        assert M([[1, 1, 1], [0, 2, 2]]).is_upper_triangular() is True
+        assert M([[1, 1, 1], [1, 2, 2]]).is_upper_triangular() is False
+
+        assert M([[1, 0, 0], [1, 2, 0]]).is_lower_triangular() is True
+        assert M([[1, 1, 0], [1, 2, 0]]).is_lower_triangular() is False
+
+
 def test_matrices_inv():
     for M, S, is_field in _all_matrices():
         if is_field:
@@ -4151,6 +4193,63 @@ def test_matrices_rref():
             assert Mr == Mr_rref
 
 
+def test_matrices_fflu():
+
+    QQ = flint.fmpq_mat
+    shape = lambda A: (A.nrows(), A.ncols())
+
+    def is_permutation(P):
+        if not P.is_square():
+            return False
+        n = P.nrows()
+        for i, row in enumerate(sorted(P.tolist(), reverse=True)):
+            if row != [int(i == j) for j in range(n)]:
+                return False
+        return True
+
+    def check_fflu(A):
+        m, n = shape(A)
+        P, L, D, U = A.fflu()
+        Dq = QQ(D)
+        assert P*A == L*Dq.inv()*U
+        assert shape(P) == shape(L) == shape(D) == (m, m)
+        assert shape(A) == shape(U) == (m, n)
+        assert is_permutation(P)
+        assert L.is_lower_triangular()
+        assert U.is_upper_triangular()
+        assert D.is_diagonal()
+
+    for M, S, is_field in _all_matrices():
+        # XXX: Add this to more matrix types...
+        if M is not flint.fmpz_mat:
+            continue
+
+        A = M([[1, 2], [3, 4]])
+        P, L, D, U = A.fflu()
+        assert P == M([[1, 0], [0, 1]])
+        assert L == M([[1, 0], [3, -2]])
+        assert D == M([[1, 0], [0, -2]])
+        assert U == M([[1, 2], [0, -2]])
+
+        check_fflu(M(0, 0, []))
+        check_fflu(M(2, 0, []))
+        check_fflu(M(0, 2, []))
+        check_fflu(M([[1]]))
+
+        check_fflu(M([[1, 2], [3, 4]]))
+        check_fflu(M([[1, 2, 3], [4, 5, 6]]))
+        check_fflu(M([[1, 2], [3, 4], [5, 6]]))
+        check_fflu(M([[1, 2], [2, 4]]))
+        check_fflu(M([[0, 0], [0, 0]]))
+        check_fflu(M([[1, 1, 1], [1, 1, 1]]))
+
+        for _ in range(10):
+            for m in range(1, 5):
+                for n in range(1, 5):
+                    A = M.randbits(m, n, 10)
+                    check_fflu(A)
+
+
 def test_matrices_solve():
     for M, S, is_field in _all_matrices():
         if is_field:
@@ -4619,13 +4718,15 @@ def test_all_tests():
     test_matrices_mul,
     test_matrices_pow,
     test_matrices_div,
+    test_matrices_properties,
     test_matrices_inv,
     test_matrices_det,
     test_matrices_charpoly,
     test_matrices_minpoly,
     test_matrices_rank,
     test_matrices_rref,
     test_matrices_solve,
+    test_matrices_fflu,
 
     test_fq_default,
     test_fq_default_poly,

src/flint/types/fmpz_mat.pyx

@@ -7,8 +7,19 @@ from flint.types.fmpz cimport any_as_fmpz
 from flint.pyflint cimport global_random_state
 from flint.types.fmpq cimport any_as_fmpq
 cimport cython
-
-from flint.flintlib.functions.fmpz cimport fmpz_set, fmpz_init, fmpz_clear
+cimport libc.stdlib
+
+from flint.flintlib.types.flint cimport SIZEOF_SLONG
+
+from flint.flintlib.functions.fmpz cimport (
+    fmpz_set,
+    fmpz_init,
+    fmpz_clear,
+    fmpz_set_si,
+    fmpz_mul,
+    fmpz_equal,
+    fmpz_equal_si,
+)
 from flint.flintlib.functions.fmpz cimport fmpz_is_zero, fmpz_is_pm1
 from flint.flintlib.types.fmpz cimport (
     fmpz_mat_struct,
@@ -316,6 +327,62 @@ cdef class fmpz_mat(flint_mat):
         fmpz_mat_pow(t.val, t.val, ee)
         return t
 
+    def is_square(self):
+        """Return whether *self* is a square *NxN* matrix."""
+        return bool(fmpz_mat_is_square(self.val))
+
+    def is_empty(self):
+        """Return whether *self* is an empty *0xN* or *Nx0* matrix."""
+        return bool(fmpz_mat_is_empty(self.val))
+
+    def is_zero(self):
+        """Return whether *self* is a zero matrix."""
+        return bool(fmpz_mat_is_zero(self.val))
+
+    def is_one(self):
+        """Return whether *self* is the identity matrix."""
+        return bool(fmpz_mat_is_one(self.val))
+
+    def is_neg_one(self):
+        """Return whether *self* is the negative identity matrix."""
+        if not self.is_square():
+            return False
+        elif not self.is_scalar():
+            return False
+        elif self.is_empty():
+            return True
+        else:
+            return bool(fmpz_equal_si(fmpz_mat_entry(self.val, 0, 0), -1))
+
+    def is_upper_triangular(self):
+        """Return whether *self* is an upper triangular matrix."""
+        for i in range(1, fmpz_mat_nrows(self.val)):
+            for j in range(min(i, fmpz_mat_ncols(self.val))):
+                if not fmpz_is_zero(fmpz_mat_entry(self.val, i, j)):
+                    return False
+        return True
+
+    def is_lower_triangular(self):
+        """Return whether *self* is a lower triangular matrix."""
+        for i in range(fmpz_mat_nrows(self.val)):
+            for j in range(i + 1, fmpz_mat_ncols(self.val)):
+                if not fmpz_is_zero(fmpz_mat_entry(self.val, i, j)):
+                    return False
+        return True
+
+    def is_diagonal(self):
+        """Return whether *self* is a diagonal matrix."""
+        return self.is_upper_triangular() and self.is_lower_triangular()
+
+    def is_scalar(self):
+        """Return whether *self* is a scalar multiple of the identity matrix."""
+        if not (self.is_square() and self.is_diagonal()):
+            return False
+        for i in range(fmpz_mat_nrows(self.val)):
+            if not fmpz_equal(fmpz_mat_entry(self.val, i, i), fmpz_mat_entry(self.val, 0, 0)):
+                return False
+        return True
+
     @classmethod
     def hadamard(cls, ulong n):
         """
@@ -563,6 +630,128 @@ cdef class fmpz_mat(flint_mat):
                 raise ZeroDivisionError("singular matrix in solve()")
             return u
 
+    def _fflu(self):
+        """
+        Fraction-free LU decomposition of *self*.
+
+            >>> A = fmpz_mat([[1, 2], [3, 4]])
+            >>> LU, d, perm, rank = A._fflu()
+            >>> LU
+            [1,  2]
+            [3, -2]
+            >>> d
+            -2
+            >>> perm
+            [0, 1]
+            >>> rank
+            2
+
+        The matrix *LU* is the LU contains both the lower and upper parts of
+        the decomposition. The integer *d* is the divisor and is up to a sign
+        the determinant when *self* is square. The list *perm* is the
+        permutation of the rows of *self*.
+
+        This is the raw output from the underlying FLINT function fmpz_mat_fflu.
+        The method :meth:`fflu` provides a more understandable representation
+        of the decomposition.
+
+        """
+        cdef fmpz d
+        cdef slong* perm
+        cdef slong r, c, rank
+        cdef fmpz_mat LU
+        r = fmpz_mat_nrows(self.val)
+        c = fmpz_mat_ncols(self.val)
+        perm = <slong*>libc.stdlib.malloc(r * SIZEOF_SLONG)
+        if perm is NULL:
+            raise MemoryError("malloc failed")
+        try:
+            for i in range(r):
+                perm[i] = i
+            LU = fmpz_mat.__new__(fmpz_mat)
+            fmpz_mat_init((<fmpz_mat>LU).val, r, c)
+            d = fmpz.__new__(fmpz)
+            rank = fmpz_mat_fflu(LU.val, d.val, perm, self.val, 0)
+            perm_int = []
+            for i in range(r):
+                perm_int.append(perm[i])
+        finally:
+            libc.stdlib.free(perm)
+
+        return LU, d, perm_int, rank
+
+    def fflu(self):
+        """
+        Fraction-free LU decomposition of *self*.
+
+        Returns a tuple (*P*, *L*, *D*, *U*) representing the the fraction-free
+        LU decomposition of a matrix *A* as
+
+            P*A = L*inv(D)*U
+
+        where *P* is a permutation matrix, *L* is lower triangular, *D* is
+        diagonal and *U* is upper triangular.
+
+            >>> A = fmpz_mat([[1, 2], [3, 4]])
+            >>> P, L, D, U = A.fflu()
+            >>> P
+            [1, 0]
+            [0, 1]
+            >>> L
+            [1,  0]
+            [3, -2]
+            >>> D
+            [1,  0]
+            [0, -2]
+            >>> U
+            [1,  2]
+            [0, -2]
+            >>> P*A == L*D.inv()*U
+            True
+
+        This method works for matrices of any shape and rank.
+
+        """
+        cdef slong r, c
+        cdef slong i, j, k, l
+        cdef fmpz di
+        cdef fmpz_mat P, L, U, D
+        r = fmpz_mat_nrows(self.val)
+        c = fmpz_mat_ncols(self.val)
+
+        U, _d, perm, _rank = self._fflu()
+
+        P = fmpz_mat(r, r)
+        for i, pi in enumerate(perm):
+            fmpz_set_si(fmpz_mat_entry(P.val, i, pi), 1)
+
+        L = fmpz_mat(r, r)
+
+        i = j = k = 0
+        while i < r and j < c:
+            if not fmpz_is_zero(fmpz_mat_entry(U.val, i, j)):
+                fmpz_set(fmpz_mat_entry(L.val, i, k), fmpz_mat_entry(U.val, i, j))
+                for l in range(i + 1, r):
+                    fmpz_set(fmpz_mat_entry(L.val, l, k), fmpz_mat_entry(U.val, l, j))
+                    fmpz_set_si(fmpz_mat_entry(U.val, l, j), 0)
+                i += 1
+                k += 1
+            j += 1
+
+        for k in range(k, r):
+            fmpz_set_si(fmpz_mat_entry(L.val, k, k), 1)
+
+        D = fmpz_mat(r, r)
+
+        if r >= 1:
+            fmpz_set(fmpz_mat_entry(D.val, 0, 0), fmpz_mat_entry(L.val, 0, 0))
+        di = fmpz(1)
+        for i in range(1, r):
+            fmpz_mul(di.val, fmpz_mat_entry(L.val, i-1, i-1), fmpz_mat_entry(L.val, i, i))
+            fmpz_set(fmpz_mat_entry(D.val, i, i), di.val)
+
+        return P, L, D, U
+
     def rref(self, inplace=False):
         """
         Computes the reduced row echelon form (rref) of *self*,