From f0754e974d7df6a0b3d0bdf76281f9d2a3cc38ed Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Sun, 9 Feb 2025 20:12:54 +0000 Subject: [PATCH 1/2] fmpz_mat: add fraction-free LU decomposition --- src/flint/test/test_all.py | 107 +++++++++++++++++++- src/flint/types/fmpz_mat.pyx | 191 ++++++++++++++++++++++++++++++++++- 2 files changed, 293 insertions(+), 5 deletions(-) diff --git a/src/flint/test/test_all.py b/src/flint/test/test_all.py index 37022cdc..eac44e83 100644 --- a/src/flint/test/test_all.py +++ b/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,6 +4718,7 @@ 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, @@ -4626,6 +4726,7 @@ def test_all_tests(): test_matrices_rank, test_matrices_rref, test_matrices_solve, + test_matrices_fflu, test_fq_default, test_fq_default_poly, diff --git a/src/flint/types/fmpz_mat.pyx b/src/flint/types/fmpz_mat.pyx index 097a90ae..683501b0 100644 --- a/src/flint/types/fmpz_mat.pyx +++ b/src/flint/types/fmpz_mat.pyx @@ -7,8 +7,17 @@ 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.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 +325,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 +628,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 = 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((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*, From 09850de78c81204090dc04874adf81eada9f7423 Mon Sep 17 00:00:00 2001 From: Oscar Benjamin Date: Sun, 9 Feb 2025 22:48:01 +0000 Subject: [PATCH 2/2] Use SIZEOF_SLONG macro to fix Windows bug --- src/flint/flintlib/types/flint.pxd | 3 +++ src/flint/types/fmpz_mat.pyx | 4 +++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/flint/flintlib/types/flint.pxd b/src/flint/flintlib/types/flint.pxd index 40fc127c..319771e7 100644 --- a/src/flint/flintlib/types/flint.pxd +++ b/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) diff --git a/src/flint/types/fmpz_mat.pyx b/src/flint/types/fmpz_mat.pyx index 683501b0..d236675c 100644 --- a/src/flint/types/fmpz_mat.pyx +++ b/src/flint/types/fmpz_mat.pyx @@ -9,6 +9,8 @@ from flint.types.fmpq cimport any_as_fmpq cimport cython cimport libc.stdlib +from flint.flintlib.types.flint cimport SIZEOF_SLONG + from flint.flintlib.functions.fmpz cimport ( fmpz_set, fmpz_init, @@ -660,7 +662,7 @@ cdef class fmpz_mat(flint_mat): cdef fmpz_mat LU r = fmpz_mat_nrows(self.val) c = fmpz_mat_ncols(self.val) - perm = libc.stdlib.malloc(r * sizeof(slong)) + perm = libc.stdlib.malloc(r * SIZEOF_SLONG) if perm is NULL: raise MemoryError("malloc failed") try: