diff --git a/.gitignore b/.gitignore index 0f4ae48e..213e156f 100644 --- a/.gitignore +++ b/.gitignore @@ -17,3 +17,4 @@ MANIFEST *.swp .python-version *.DS_Store +.venv diff --git a/doc/source/fmpz_mod_mat.rst b/doc/source/fmpz_mod_mat.rst new file mode 100644 index 00000000..73e9f012 --- /dev/null +++ b/doc/source/fmpz_mod_mat.rst @@ -0,0 +1,8 @@ +**fmpz_mod_mat** -- matrices over integers mod n for arbitrary n +=============================================================================== + +.. autoclass :: flint.fmpz_mod_mat + :members: + :inherited-members: + :undoc-members: + diff --git a/doc/source/index.rst b/doc/source/index.rst index 2df7aa3f..0a1acb33 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -58,6 +58,7 @@ Matrix types fmpz_mat.rst fmpq_mat.rst nmod_mat.rst + fmpz_mod_mat.rst arb_mat.rst acb_mat.rst diff --git a/setup.py b/setup.py index 60b2d04b..fa727ab7 100644 --- a/setup.py +++ b/setup.py @@ -76,18 +76,30 @@ ext_files = [ ("flint.pyflint", ["src/flint/pyflint.pyx"]), + + ("flint.flint_base.flint_base", ["src/flint/flint_base/flint_base.pyx"]), + ("flint.flint_base.flint_context", ["src/flint/flint_base/flint_context.pyx"]), + ("flint.types.fmpz", ["src/flint/types/fmpz.pyx"]), ("flint.types.fmpz_poly", ["src/flint/types/fmpz_poly.pyx"]), + ("flint.types.fmpz_mpoly", ["src/flint/types/fmpz_mpoly.pyx"]), ("flint.types.fmpz_mat", ["src/flint/types/fmpz_mat.pyx"]), ("flint.types.fmpz_series", ["src/flint/types/fmpz_series.pyx"]), + ("flint.types.fmpq", ["src/flint/types/fmpq.pyx"]), ("flint.types.fmpq_poly", ["src/flint/types/fmpq_poly.pyx"]), ("flint.types.fmpq_mat", ["src/flint/types/fmpq_mat.pyx"]), ("flint.types.fmpq_series", ["src/flint/types/fmpq_series.pyx"]), + ("flint.types.nmod", ["src/flint/types/nmod.pyx"]), ("flint.types.nmod_poly", ["src/flint/types/nmod_poly.pyx"]), ("flint.types.nmod_mat", ["src/flint/types/nmod_mat.pyx"]), ("flint.types.nmod_series", ["src/flint/types/nmod_series.pyx"]), + + ("flint.types.fmpz_mod", ["src/flint/types/fmpz_mod.pyx"]), + ("flint.types.fmpz_mod_poly", ["src/flint/types/fmpz_mod_poly.pyx"]), + ("flint.types.fmpz_mod_mat", ["src/flint/types/fmpz_mod_mat.pyx"]), + ("flint.types.arf", ["src/flint/types/arf.pyx"]), ("flint.types.arb", ["src/flint/types/arb.pyx"]), ("flint.types.arb_poly", ["src/flint/types/arb_poly.pyx"]), @@ -97,15 +109,10 @@ ("flint.types.acb_poly", ["src/flint/types/acb_poly.pyx"]), ("flint.types.acb_mat", ["src/flint/types/acb_mat.pyx"]), ("flint.types.acb_series", ["src/flint/types/acb_series.pyx"]), - ("flint.types.fmpz_mpoly", ["src/flint/types/fmpz_mpoly.pyx"]), - ("flint.types.fmpz_mod", ["src/flint/types/fmpz_mod.pyx"]), - ("flint.types.fmpz_mod_poly", ["src/flint/types/fmpz_mod_poly.pyx"]), + ("flint.types.dirichlet", ["src/flint/types/dirichlet.pyx"]), - ("flint.flint_base.flint_base", ["src/flint/flint_base/flint_base.pyx"]), - ("flint.flint_base.flint_context", ["src/flint/flint_base/flint_context.pyx"]), - # Helper for unittests - ("flint.functions.showgood", ["src/flint/functions/showgood.pyx"]), + ("flint.functions.showgood", ["src/flint/functions/showgood.pyx"]), ] ext_options = { diff --git a/src/flint/__init__.py b/src/flint/__init__.py index bda516f9..d13ba26e 100644 --- a/src/flint/__init__.py +++ b/src/flint/__init__.py @@ -1,16 +1,25 @@ from .pyflint import * + from .types.fmpz import * from .types.fmpz_poly import * from .types.fmpz_mat import * from .types.fmpz_series import * + from .types.fmpq import * from .types.fmpq_poly import * from .types.fmpq_mat import * from .types.fmpq_series import * + from .types.nmod import * from .types.nmod_poly import * from .types.nmod_mat import * from .types.nmod_series import * + +from .types.fmpz_mpoly import * +from .types.fmpz_mod import * +from .types.fmpz_mod_poly import * +from .types.fmpz_mod_mat import fmpz_mod_mat + from .types.arf import * from .types.arb import * from .types.arb_poly import * @@ -20,9 +29,7 @@ from .types.acb_poly import * from .types.acb_mat import * from .types.acb_series import * -from .types.fmpz_mpoly import * -from .types.fmpz_mod import * -from .types.fmpz_mod_poly import * + from .types.dirichlet import * from .functions.showgood import good, showgood diff --git a/src/flint/flint_base/flint_base.pyx b/src/flint/flint_base/flint_base.pyx index fbb61478..7d42d397 100644 --- a/src/flint/flint_base/flint_base.pyx +++ b/src/flint/flint_base/flint_base.pyx @@ -1,5 +1,7 @@ +from flint.flintlib.flint cimport FLINT_BITS as _FLINT_BITS from flint.flint_base.flint_context cimport thectx + cdef class flint_elem: def __repr__(self): if thectx.pretty: @@ -134,8 +136,6 @@ cdef class flint_mat(flint_elem): """ def repr(self): - if thectx.pretty: - return str(self) # XXX return "%s(%i, %i, [%s])" % (type(self).__name__, self.nrows(), self.ncols(), (", ".join(map(str, self.entries())))) diff --git a/src/flint/flintlib/fmpz_mod_mat.pxd b/src/flint/flintlib/fmpz_mod_mat.pxd new file mode 100644 index 00000000..c1b08f48 --- /dev/null +++ b/src/flint/flintlib/fmpz_mod_mat.pxd @@ -0,0 +1,75 @@ +from flint.flintlib.flint cimport ulong, slong, fmpz_struct, flint_rand_t +from flint.flintlib.fmpz cimport fmpz_t +from flint.flintlib.fmpz_mod cimport fmpz_mod_ctx_t +from flint.flintlib.fmpz_mat cimport fmpz_mat_t +from flint.flintlib.fmpz_mod_poly cimport fmpz_mod_poly_t + + +cdef extern from "flint/fmpz_mod_mat.h": + ctypedef struct fmpz_mod_mat_struct: + fmpz_mat_t mat + fmpz_t mod + ctypedef fmpz_mod_mat_struct fmpz_mod_mat_t[1] + + +cdef extern from "flint/fmpz_mod_mat.h": + # This is not exposed in the docs: + int fmpz_mod_mat_equal(const fmpz_mod_mat_t mat1, const fmpz_mod_mat_t mat2); + + +cdef extern from "flint/fmpz_mod_mat.h": + fmpz_struct * fmpz_mod_mat_entry(const fmpz_mod_mat_t mat, slong i, slong j) + void fmpz_mod_mat_set_entry(fmpz_mod_mat_t mat, slong i, slong j, const fmpz_t val) + void fmpz_mod_mat_init(fmpz_mod_mat_t mat, slong rows, slong cols, const fmpz_t n) + void fmpz_mod_mat_init_set(fmpz_mod_mat_t mat, const fmpz_mod_mat_t src) + void fmpz_mod_mat_clear(fmpz_mod_mat_t mat) + slong fmpz_mod_mat_nrows(const fmpz_mod_mat_t mat) + slong fmpz_mod_mat_ncols(const fmpz_mod_mat_t mat) + void _fmpz_mod_mat_set_mod(fmpz_mod_mat_t mat, const fmpz_t n) + void fmpz_mod_mat_one(fmpz_mod_mat_t mat) + void fmpz_mod_mat_zero(fmpz_mod_mat_t mat) + void fmpz_mod_mat_swap(fmpz_mod_mat_t mat1, fmpz_mod_mat_t mat2) + void fmpz_mod_mat_swap_entrywise(fmpz_mod_mat_t mat1, fmpz_mod_mat_t mat2) + int fmpz_mod_mat_is_empty(const fmpz_mod_mat_t mat) + int fmpz_mod_mat_is_square(const fmpz_mod_mat_t mat) + void _fmpz_mod_mat_reduce(fmpz_mod_mat_t mat) + void fmpz_mod_mat_randtest(fmpz_mod_mat_t mat, flint_rand_t state) + void fmpz_mod_mat_window_init(fmpz_mod_mat_t window, const fmpz_mod_mat_t mat, slong r1, slong c1, slong r2, slong c2) + void fmpz_mod_mat_window_clear(fmpz_mod_mat_t window) + void fmpz_mod_mat_concat_horizontal(fmpz_mod_mat_t res, const fmpz_mod_mat_t mat1, const fmpz_mod_mat_t mat2) + void fmpz_mod_mat_concat_vertical(fmpz_mod_mat_t res, const fmpz_mod_mat_t mat1, const fmpz_mod_mat_t mat2) + void fmpz_mod_mat_print_pretty(const fmpz_mod_mat_t mat) + int fmpz_mod_mat_is_zero(const fmpz_mod_mat_t mat) + void fmpz_mod_mat_set(fmpz_mod_mat_t B, const fmpz_mod_mat_t A) + void fmpz_mod_mat_transpose(fmpz_mod_mat_t B, const fmpz_mod_mat_t A) + void fmpz_mod_mat_set_fmpz_mat(fmpz_mod_mat_t A, const fmpz_mat_t B) + void fmpz_mod_mat_get_fmpz_mat(fmpz_mat_t A, const fmpz_mod_mat_t B) + void fmpz_mod_mat_add(fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + void fmpz_mod_mat_sub(fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + void fmpz_mod_mat_neg(fmpz_mod_mat_t B, const fmpz_mod_mat_t A) + void fmpz_mod_mat_scalar_mul_si(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, slong c) + void fmpz_mod_mat_scalar_mul_ui(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, ulong c) + void fmpz_mod_mat_scalar_mul_fmpz(fmpz_mod_mat_t B, const fmpz_mod_mat_t A, fmpz_t c) + void fmpz_mod_mat_mul(fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + # unimported types {'thread_pool_handle'} + # void _fmpz_mod_mat_mul_classical_threaded_pool_op(fmpz_mod_mat_t D, const fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B, int op, thread_pool_handle * threads, slong num_threads) + void _fmpz_mod_mat_mul_classical_threaded_op(fmpz_mod_mat_t D, const fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B, int op) + void fmpz_mod_mat_mul_classical_threaded(fmpz_mod_mat_t C, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + void fmpz_mod_mat_sqr(fmpz_mod_mat_t B, const fmpz_mod_mat_t A) + void fmpz_mod_mat_mul_fmpz_vec(fmpz_struct * c, const fmpz_mod_mat_t A, const fmpz_struct * b, slong blen) + void fmpz_mod_mat_mul_fmpz_vec_ptr(fmpz_struct * const * c, const fmpz_mod_mat_t A, const fmpz_struct * const * b, slong blen) + void fmpz_mod_mat_fmpz_vec_mul(fmpz_struct * c, const fmpz_struct * a, slong alen, const fmpz_mod_mat_t B) + void fmpz_mod_mat_fmpz_vec_mul_ptr(fmpz_struct * const * c, const fmpz_struct * const * a, slong alen, const fmpz_mod_mat_t B) + void fmpz_mod_mat_trace(fmpz_t trace, const fmpz_mod_mat_t mat) + slong fmpz_mod_mat_rref(slong * perm, fmpz_mod_mat_t mat) + void fmpz_mod_mat_strong_echelon_form(fmpz_mod_mat_t mat) + slong fmpz_mod_mat_howell_form(fmpz_mod_mat_t mat) + int fmpz_mod_mat_inv(fmpz_mod_mat_t B, fmpz_mod_mat_t A) + slong fmpz_mod_mat_lu(slong * P, fmpz_mod_mat_t A, int rank_check) + void fmpz_mod_mat_solve_tril(fmpz_mod_mat_t X, const fmpz_mod_mat_t L, const fmpz_mod_mat_t B, int unit) + void fmpz_mod_mat_solve_triu(fmpz_mod_mat_t X, const fmpz_mod_mat_t U, const fmpz_mod_mat_t B, int unit) + int fmpz_mod_mat_solve(fmpz_mod_mat_t X, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + int fmpz_mod_mat_can_solve(fmpz_mod_mat_t X, const fmpz_mod_mat_t A, const fmpz_mod_mat_t B) + void fmpz_mod_mat_similarity(fmpz_mod_mat_t M, slong r, fmpz_t d) + void fmpz_mod_mat_charpoly(fmpz_mod_poly_t p, const fmpz_mod_mat_t M, const fmpz_mod_ctx_t ctx) + void fmpz_mod_mat_minpoly(fmpz_mod_poly_t p, const fmpz_mod_mat_t M, const fmpz_mod_ctx_t ctx) diff --git a/src/flint/test/__main__.py b/src/flint/test/__main__.py index c22e5f48..637fdb46 100644 --- a/src/flint/test/__main__.py +++ b/src/flint/test/__main__.py @@ -60,6 +60,7 @@ def run_doctests(verbose=None): flint.types.fmpz_series, flint.types.fmpz_mod, flint.types.fmpz_mod_poly, + flint.types.fmpz_mod_mat, flint.types.fmpq, flint.types.fmpq_poly, flint.types.fmpq_mat, diff --git a/src/flint/test/test.py b/src/flint/test/test.py index aab785ee..13a3481f 100644 --- a/src/flint/test/test.py +++ b/src/flint/test/test.py @@ -557,7 +557,7 @@ def test_fmpz_mat(): assert raises(lambda: M([[1],[2,3]]), ValueError) assert raises(lambda: M(None), TypeError) assert raises(lambda: M(2,2,[1,2,3]), ValueError) - assert raises(lambda: M(2,2,2,2), ValueError) + assert raises(lambda: M(2,2,2,2), TypeError) assert M([[1,2,3],[4,5,6]]) == M(2,3,[1,2,3,4,5,6]) assert raises(lambda: M([[1]]) < M([[2]]), TypeError) assert (M([[1]]) == 1) is False @@ -571,15 +571,15 @@ def test_fmpz_mat(): def set_bad(i,j): D[i,j] = -1 # XXX: Should be IndexError - raises(lambda: set_bad(2,0), ValueError) - raises(lambda: set_bad(0,2), ValueError) - raises(lambda: D[0,2], ValueError) - raises(lambda: D[0,2], ValueError) + raises(lambda: set_bad(2,0), IndexError) + raises(lambda: set_bad(0,2), IndexError) + raises(lambda: D[0,2], IndexError) + raises(lambda: D[0,2], IndexError) # XXX: Negative indices? - raises(lambda: set_bad(-1,0), ValueError) - raises(lambda: set_bad(0,-1), ValueError) - raises(lambda: D[-1,0], ValueError) - raises(lambda: D[0,-1], ValueError) + raises(lambda: set_bad(-1,0), IndexError) + raises(lambda: set_bad(0,-1), IndexError) + raises(lambda: D[-1,0], IndexError) + raises(lambda: D[0,-1], IndexError) assert M.hadamard(2) == M([[1,1],[1,-1]]) assert raises(lambda: M.hadamard(3), ValueError) assert M.hadamard(2).is_hadamard() is True @@ -1067,7 +1067,7 @@ def test_fmpq_mat(): assert raises(lambda: Q(2,3,[1,2,3,4,5]), ValueError) assert raises(lambda: Q([[1,2,3],[4,[],6]]), TypeError) assert raises(lambda: Q(2,3,[1,2,3,4,[],6]), TypeError) - assert raises(lambda: Q(2,3,[1,2],[3,4]), ValueError) + assert raises(lambda: Q(2,3,[1,2],[3,4]), TypeError) assert bool(Q([[1]])) is True assert bool(Q([[0]])) is False assert raises(lambda: Q([[1]]) < Q([[0]]), TypeError) @@ -1078,13 +1078,12 @@ def test_fmpq_mat(): # XXX: Negative indices should probably be allowed def set_bad(i): M[i,0] = -1 - raises(lambda: M[-1,0], ValueError) - raises(lambda: M[0,-1], ValueError) - raises(lambda: set_bad(-1), ValueError) - # XXX: Should be IndexError - raises(lambda: M[2,0], ValueError) - raises(lambda: M[0,2], ValueError) - raises(lambda: set_bad(2), ValueError) + raises(lambda: M[-1,0], IndexError) + raises(lambda: M[0,-1], IndexError) + raises(lambda: set_bad(-1), IndexError) + raises(lambda: M[2,0], IndexError) + raises(lambda: M[0,2], IndexError) + raises(lambda: set_bad(2), IndexError) assert Q([[1,2,3],[4,5,6]]).transpose() == Q([[1,4],[2,5],[3,6]]) raises(lambda: M + [], TypeError) raises(lambda: M - [], TypeError) @@ -1488,9 +1487,12 @@ def test_nmod_mat(): assert A*(B*C) == (A*B)*C assert bool(M(2,2,[0,0,0,0],17)) == False assert bool(M(2,2,[0,0,0,1],17)) == True - ctx.pretty = False - assert repr(M(2,2,[1,2,3,4],17)) == 'nmod_mat(2, 2, [1, 2, 3, 4], 17)' - ctx.pretty = True + pretty = ctx.pretty + try: + ctx.pretty = False + assert repr(M(2,2,[1,2,3,4],17)) == 'nmod_mat(2, 2, [1, 2, 3, 4], 17)' + finally: + ctx.pretty = pretty assert str(M(2,2,[1,2,3,4],17)) == '[1, 2]\n[3, 4]' assert repr(M(2,2,[1,2,3,4],17)) == '[1, 2]\n[3, 4]' assert M(1,2,[3,4],17) / 3 == M(1,2,[3,4],17) * (~G(3,17)) @@ -1504,7 +1506,7 @@ def test_nmod_mat(): assert raises(lambda: M(None,17), TypeError) assert M(2,3,17) == M(2,3,[0,0,0,0,0,0],17) assert raises(lambda: M(2,3,[0,0,0,0,0],17), ValueError) - assert raises(lambda: M(2,3,[0,1],[1,2],17), ValueError) + assert raises(lambda: M(2,3,[0,1],[1,2],17), TypeError) assert M([[1,2,3],[4,5,6]], 5) == M(2,3,[1,2,3,4,5,6], 5) assert raises(lambda: M([[0]],13) < M([[1]],13), TypeError) assert (M([[1]],17) == M([[1]],13)) is False @@ -1522,18 +1524,17 @@ def test_nmod_mat(): def set_bad(i,j): M3[i,j] = 2 # XXX: negative indices should be allowed - assert raises(lambda: M3[-1,0], ValueError) - assert raises(lambda: M3[0,-1], ValueError) - assert raises(lambda: set_bad(-1,0), ValueError) - assert raises(lambda: set_bad(0,-1), ValueError) - # XXX: Should be IndexError - assert raises(lambda: M3[2,0], ValueError) - assert raises(lambda: M3[0,2], ValueError) - assert raises(lambda: set_bad(2,0), ValueError) - assert raises(lambda: set_bad(0,2), ValueError) + assert raises(lambda: M3[-1,0], IndexError) + assert raises(lambda: M3[0,-1], IndexError) + assert raises(lambda: set_bad(-1,0), IndexError) + assert raises(lambda: set_bad(0,-1), IndexError) + assert raises(lambda: M3[2,0], IndexError) + assert raises(lambda: M3[0,2], IndexError) + assert raises(lambda: set_bad(2,0), IndexError) + assert raises(lambda: set_bad(0,2), IndexError) def set_bad2(): M3[0,0] = 1.5 - assert raises(set_bad2, ValueError) + assert raises(set_bad2, TypeError) assert raises(lambda: M3 + [], TypeError) assert raises(lambda: M3 - [], TypeError) assert raises(lambda: M3 * [], TypeError) @@ -2217,6 +2218,64 @@ def test_fmpz_mod_poly(): assert [f(x) for x in l] == f.multipoint_evaluate(l) +def test_fmpz_mod_mat(): + c11 = flint.fmpz_mod_ctx(11) + c13 = flint.fmpz_mod_ctx(13) + + M11_1 = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c11) + M11_2 = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c11) + M13 = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c13) + assert (M11_1 == M11_2) is True + assert (M11_1 != M11_2) is False + assert (M11_1 == M13) is False + assert (M11_1 != M13) is True + + Mn_11 = flint.nmod_mat([[1,2,3],[4,5,6],[7,8,9]], 11) + Mn_13 = flint.nmod_mat([[1,2,3],[4,5,6],[7,8,9]], 13) + assert (M11_1 == Mn_11) is True + assert (M11_1 != Mn_11) is False + assert (M11_1 == Mn_13) is False + assert (M11_1 != Mn_13) is True + assert (Mn_11 == M11_1) is True + assert (Mn_11 != M11_1) is False + assert (Mn_13 == M11_1) is False + assert (Mn_13 != M11_1) is True + + M = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c11) + assert M.nrows() == 3 + assert M.ncols() == 3 + assert M.entries() == [flint.fmpz_mod(i, c11) for i in [1,2,3,4,5,6,7,8,9]] + assert type(M[0,0]) == flint.fmpz_mod + assert M[0,0] == 1 + assert raises(lambda: flint.fmpz_mod_mat([[1]]), TypeError) + assert raises(lambda: flint.fmpz_mod_mat([[1,2],[3,4]], c11, c11), TypeError) + assert raises(lambda: flint.fmpz_mod_mat([[1,2],[3,4]], None), TypeError) + + M11 = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c11) + M13 = flint.fmpz_mod_mat([[1,2,3],[4,5,6],[7,8,9]], c13) + Mz = flint.fmpz_mat([[1,2,3],[4,5,6],[7,8,9]]) + Mn_11 = flint.nmod_mat([[1,2,3],[4,5,6],[7,8,9]], 11) + Mn_13 = flint.nmod_mat([[1,2,3],[4,5,6],[7,8,9]], 13) + assert raises(lambda: flint.fmpz_mod_mat(Mz), TypeError) + assert flint.fmpz_mod_mat(Mz, c11) == M11 + assert flint.fmpz_mod_mat(Mz, c13) == M13 + assert flint.fmpz_mod_mat(Mn_11) == M11 + assert flint.fmpz_mod_mat(Mn_13) == M13 + assert flint.fmpz_mod_mat(Mn_11, c11) == M11 + assert flint.fmpz_mod_mat(Mn_13, c13) == M13 + assert raises(lambda: flint.fmpz_mod_mat(Mn_11, c13), TypeError) + assert raises(lambda: flint.fmpz_mod_mat(Mn_13, c11), TypeError) + assert flint.fmpz_mod_mat(M11) == M11 + assert flint.fmpz_mod_mat(M13) == M13 + assert flint.fmpz_mod_mat(M11, c11) == M11 + assert flint.fmpz_mod_mat(M13, c13) == M13 + assert raises(lambda: flint.fmpz_mod_mat(M11, c13), TypeError) + assert raises(lambda: flint.fmpz_mod_mat(M13, c11), TypeError) + + A = flint.arb_mat([[1,2,3],[4,5,6],[7,8,9]]) + assert raises(lambda: flint.fmpz_mod_mat(A, c11), TypeError) + + def _all_polys(): return [ # (poly_type, scalar_type, is_field) @@ -2425,14 +2484,521 @@ def setbad(obj, i, val): if is_field: assert P([1, 2, 1]).integral() == P([0, 1, 1, S(1)/3]) + +def _all_matrices(): + """Return a list of matrix types and scalar types.""" + R163 = flint.fmpz_mod_ctx(163) + R127 = flint.fmpz_mod_ctx(2**127 - 1) + R255 = flint.fmpz_mod_ctx(2**255 - 19) + return [ + # (matrix_type, scalar_type, is_field) + (flint.fmpz_mat, flint.fmpz, False), + (flint.fmpq_mat, flint.fmpq, True), + (lambda *a: flint.nmod_mat(*a, 17), lambda x: flint.nmod(x, 17), True), + (lambda *a: flint.fmpz_mod_mat(*a, R163), lambda x: flint.fmpz_mod(x, R163), True), + (lambda *a: flint.fmpz_mod_mat(*a, R127), lambda x: flint.fmpz_mod(x, R127), True), + (lambda *a: flint.fmpz_mod_mat(*a, R255), lambda x: flint.fmpz_mod(x, R255), True), + ] + + +def _coercible_matrix_types(mat_type): + """Return a list of matrix types that can be coerced to `mat_type`.""" + M = mat_type([[1, 2], [3, 4]]) + + if type(M) is flint.fmpz_mat: + return [ + flint.fmpz_mat + ] + elif type(M) is flint.fmpq_mat: + return [ + flint.fmpz_mat, + flint.fmpq_mat + ] + elif type(M) is flint.nmod_mat: + m = M.modulus() + return [ + flint.fmpz_mat, + lambda *a: flint.nmod_mat(*a, m) + ] + elif type(M) is flint.fmpz_mod_mat: + m = M.modulus() + c = flint.fmpz_mod_ctx(m) + mat_types = [ + flint.fmpz_mat, + lambda *a: flint.fmpz_mod_mat(*a, c) + ] + if m < 2**64 - 1: + mat_types.append(lambda *a: flint.nmod_mat(*a, m)) + return mat_types + else: + assert False + + +def _compatible_matrix_types(mat_type): + """Return a list of matrix types that can be added to `mat_type`.""" + M = mat_type([[1, 2], [3, 4]]) + + if type(M) is flint.fmpz_mat: + return [ + flint.fmpz_mat + ] + elif type(M) is flint.fmpq_mat: + return [ + # flint.fmpz_mat, + flint.fmpq_mat + ] + elif type(M) is flint.nmod_mat: + m = M.modulus() + return [ + flint.fmpz_mat, + lambda *a: flint.nmod_mat(*a, m) + ] + elif type(M) is flint.fmpz_mod_mat: + m = M.modulus() + c = flint.fmpz_mod_ctx(m) + mat_types = [ + flint.fmpz_mat, + lambda *a: flint.fmpz_mod_mat(*a, c) + ] + if m < 2**64 - 1: + mat_types.append(lambda *a: flint.nmod_mat(*a, m)) + return mat_types + else: + assert False + + +def _incompatible_matrix_types(mat_type): + M = mat_type([[1, 2], [3, 4]]) + + if type(M) is flint.fmpz_mat: + return [] + elif type(M) is flint.fmpq_mat: + return [ + lambda *a: flint.nmod_mat(*a, 17), + lambda *a: flint.fmpz_mod_mat(*a, flint.fmpz_mod_ctx(163)), + ] + elif type(M) is flint.nmod_mat: + m = M.modulus() + m2 = m - 1 + return [ + flint.fmpq_mat, + lambda *a: flint.nmod_mat(*a, m2), + ] + elif type(M) is flint.fmpz_mod_mat: + m = M.modulus() + m2 = m - 1 + mat_types = [ + flint.fmpq_mat, + lambda *a: flint.fmpz_mod_mat(*a, flint.fmpz_mod_ctx(m2)), + ] + if m < 2**64 - 1: + mat_types.append(lambda *a: flint.nmod_mat(*a, m2)) + return mat_types + else: + assert False + + +def _poly_type_from_matrix_type(mat_type): + M = mat_type([[1, 2], [3, 4]]) + + if type(M) is flint.fmpz_mat: + return flint.fmpz_poly + elif type(M) is flint.fmpq_mat: + return flint.fmpq_poly + elif type(M) is flint.nmod_mat: + m = M.modulus() + return lambda *a: flint.nmod_poly(*a, m) + elif type(M) is flint.fmpz_mod_mat: + m = M.modulus() + ctx = flint.fmpz_mod_ctx(m) + pctx = flint.fmpz_mod_poly_ctx(ctx) + return lambda *a: flint.fmpz_mod_poly(*a, pctx) + else: + assert False + + +def test_matrices_eq(): + for M, S, is_field in _all_matrices(): + A1 = M([[1, 2], [3, 4]]) + A2 = M([[1, 2], [3, 4]]) + B = M([[5, 6], [7, 8]]) + assert (A1 == A2) is True + assert (A1 != A2) is False + assert (A1 == B) is False + assert (A1 != B) is True + assert (A1 == None) is False + assert (A1 != None) is True + assert (None == A1) is False + assert (None != A1) is True + assert raises(lambda: A1 < A2, TypeError) + assert raises(lambda: A1 <= A2, TypeError) + assert raises(lambda: A1 > A2, TypeError) + assert raises(lambda: A1 >= A2, TypeError) + + for M2 in _incompatible_matrix_types(M): + A1 = M([[1, 2], [3, 4]]) + A2 = M2([[1, 2], [3, 4]]) + assert (A1 == A2) is False + assert (A1 != A2) is True + + +def test_matrices_constructor(): + for M, S, is_field in _all_matrices(): + assert raises(lambda: M(), TypeError) + + # Empty matrices + assert M([]).nrows() == 0 + assert M([]).ncols() == 0 + assert M([]) != M([[]]) + assert M([[]]).nrows() == 1 + assert M([[]]).ncols() == 0 + assert M([]) == M(0, 0, []) == M(0, 0) + assert M([[]]) == M(1, 0, []) == M(1, 0) + assert M([[], []]) == M(2, 0, []) == M(2, 0) + assert M(0, 2, []) == M(0, 2) != M(2, 0) + + # Basic construction + M1234 = M([[1, 2], [3, 4]]) + assert M1234.nrows() == 2 + assert M1234.ncols() == 2 + assert M1234 == M([[1, 2], [3, 4]]) + assert M1234 != M([[1, 2], [3, 5]]) + assert M1234 == M(2, 2, [1, 2, 3, 4]) + assert M1234.entries() == [S(1), S(2), S(3), S(4)] + assert M1234.table() == [[S(1), S(2)], [S(3), S(4)]] + assert M1234[0, 0] == S(1) + assert M1234[0, 1] == S(2) + assert M1234[1, 0] == S(3) + assert M1234[1, 1] == S(4) + assert type(M1234[0, 0]) is type(S(1)) + assert M1234 == M([[S(1), S(2)], [S(3), S(4)]]) + assert M1234 == M(2, 2, [S(1), S(2), S(3), S(4)]) + + # Non-square matrices + M6 = M([[1, 2, 3], [4, 5, 6]]) + assert M6.nrows() == 2 + assert M6.ncols() == 3 + assert M6 == M(2, 3, [1, 2, 3, 4, 5, 6]) + assert M6.entries() == [S(1), S(2), S(3), S(4), S(5), S(6)] + assert M6 != M1234 + + # Bad list inputs + assert raises(lambda: M(None), TypeError) + assert raises(lambda: M([None]), TypeError) + assert raises(lambda: M(2, 3, [1, 2, 3, 4, 5]), ValueError) + assert raises(lambda: M([[1, 2], [3, 4], [5, 6, 7]]), ValueError) + assert raises(lambda: M([[1, 2], [3, 4], [5, 6, 7, 8]]), ValueError) + assert raises(lambda: M([[1, 2], [3, None]]), TypeError) + + # Construct from a coercible matrix + M1234 = M([[1, 2], [3, 4]]) + for M2 in _coercible_matrix_types(M): + if type(M2([[1]])) is flint.nmod_mat: + continue + M1234_2 = M(M2([[1, 2], [3, 4]])) + assert M1234_2 == M1234 + assert type(M1234_2) is type(M1234) + + +def _matrix_repr(M): + item = M[0, 0] + if type(item) is flint.fmpz: + return f"fmpz_mat({M.nrows()}, {M.ncols()}, {M.entries()!r})" + elif type(item) is flint.fmpq: + return f"fmpq_mat({M.nrows()}, {M.ncols()}, {M.entries()!r})" + elif type(item) is flint.nmod: + return f"nmod_mat({M.nrows()}, {M.ncols()}, {M.entries()!r}, {M.modulus()})" + elif type(item) is flint.fmpz_mod: + return f"fmpz_mod_mat({M.nrows()}, {M.ncols()}, {M.entries()!r}, {M.modulus()})" + else: + assert False + + +def test_matrices_strrepr(): + for M, S, is_field in _all_matrices(): + A = M([[1, 2], [3, 4]]) + A_str = "[1, 2]\n[3, 4]" + A_repr = _matrix_repr(A) + + assert A.str() == A_str, type(A).__name__ + assert A.repr() == A_repr, type(A).__name__ + + # str always returns a pretty result + assert str(A) == A_str, type(A).__name__ + + # repr depends on ctx.pretty + pretty = ctx.pretty + try: + ctx.pretty = True + assert repr(A) == A_str, type(A).__name__ + ctx.pretty = False + assert repr(A) == A_repr, type(A).__name__ + finally: + ctx.pretty = pretty + + +def test_matrices_getitem(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + assert M1234[0, 0] == S(1) + assert M1234[0, 1] == S(2) + assert M1234[1, 0] == S(3) + assert M1234[1, 1] == S(4) + assert raises(lambda: M1234[0, 2], IndexError) + assert raises(lambda: M1234[2, 0], IndexError) + assert raises(lambda: M1234[2, 2], IndexError) + # XXX: Should negative indices be allowed? + assert raises(lambda: M1234[-1, 0], IndexError) + assert raises(lambda: M1234[0, -1], IndexError) + assert raises(lambda: M1234[-1, -1], IndexError) + + +def test_matrices_setitem(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + + assert M1234[0, 0] == S(1) + M1234[0, 0] = S(5) + assert M1234[0, 0] == S(5) + assert M1234.table() == [[S(5), S(2)], [S(3), S(4)]] + + def setbad(obj, key, val): + obj[key] = val + + assert raises(lambda: setbad(M1234, (0,0), None), TypeError) + assert raises(lambda: setbad(M1234, (0,None), 1), TypeError) + assert raises(lambda: setbad(M1234, (None,0), 1), TypeError) + assert raises(lambda: setbad(M1234, None, 1), TypeError) + + assert raises(lambda: setbad(M1234, (0,2), 1), IndexError) + assert raises(lambda: setbad(M1234, (2,0), 1), IndexError) + assert raises(lambda: setbad(M1234, (2,2), 1), IndexError) + # XXX: Should negative indices be allowed? + assert raises(lambda: setbad(M1234, (-1,0), 1), IndexError) + assert raises(lambda: setbad(M1234, (0,-1), 1), IndexError) + assert raises(lambda: setbad(M1234, (-1,-1), 1), IndexError) + + +def test_matrices_bool(): + for M, S, is_field in _all_matrices(): + assert bool(M([])) is False + assert bool(M([[0]])) is False + assert bool(M([[1]])) is True + assert bool(M([[0, 0], [0, 0]])) is False + assert bool(M([[1, 0], [0, 0]])) is True + assert bool(M([[0, 0], [0, 1]])) is True + assert bool(M([[1, 0], [0, 1]])) is True + + +def test_matrices_pos_neg(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + assert +M1234 == M1234 + assert -M1234 == M([[-1, -2], [-3, -4]]) + + +def test_matrices_add(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + M5678 = M([[5, 6], [7, 8]]) + assert M1234 + M5678 == M([[6, 8], [10, 12]]) + assert raises(lambda: M1234 + 1, TypeError) + assert raises(lambda: 1 + M1234, TypeError) + assert raises(lambda: M1234 + None, TypeError) + assert raises(lambda: None + M1234, TypeError) + assert raises(lambda: M1234 + M([[1, 2, 3], [4, 5, 6]]), ValueError) + for M2 in _compatible_matrix_types(M): + assert M1234 + M2([[1, 2], [3, 4]]) == M([[2, 4], [6, 8]]) + assert M2([[1, 2], [3, 4]]) + M1234 == M([[2, 4], [6, 8]]) + assert raises(lambda: M1234 + M2([[1, 2, 3], [4, 5, 6]]), ValueError) + assert raises(lambda: M2([[1, 2, 3], [4, 5, 6]]) + M1234, ValueError) + for M2 in _incompatible_matrix_types(M): + assert raises(lambda: M1234 + M2([[1, 2], [3, 4]]), (TypeError, ValueError)) + assert raises(lambda: M2([[1, 2], [3, 4]]) + M1234, (TypeError, ValueError)) + + +def test_matrices_sub(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + M5678 = M([[5, 6], [7, 8]]) + assert M1234 - M5678 == M([[-4, -4], [-4, -4]]) + assert raises(lambda: M1234 - 1, TypeError) + assert raises(lambda: 1 - M1234, TypeError) + assert raises(lambda: M1234 - None, TypeError) + assert raises(lambda: None - M1234, TypeError) + assert raises(lambda: M1234 - M([[1, 2, 3], [4, 5, 6]]), ValueError) + for M2 in _compatible_matrix_types(M): + assert M1234 - M2([[1, 2], [3, 4]]) == M([[0, 0], [0, 0]]) + assert M2([[1, 2], [3, 4]]) - M1234 == M([[0, 0], [0, 0]]) + assert raises(lambda: M1234 - M2([[1, 2, 3], [4, 5, 6]]), ValueError) + assert raises(lambda: M2([[1, 2, 3], [4, 5, 6]]) - M1234, ValueError) + for M2 in _incompatible_matrix_types(M): + assert raises(lambda: M1234 - M2([[1, 2], [3, 4]]), (TypeError, ValueError)) + assert raises(lambda: M2([[1, 2], [3, 4]]) - M1234, (TypeError, ValueError)) + + +def test_matrices_mul(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + M5678 = M([[5, 6], [7, 8]]) + assert M1234 * M5678 == M([[19, 22], [43, 50]]) + M21 = M([[2], [1]]) + assert M1234 * M21 == M([[4], [10]]) + + assert 2 * M1234 == M([[2, 4], [6, 8]]) + assert M1234 * 2 == M([[2, 4], [6, 8]]) + assert S(2) * M1234 == M([[2, 4], [6, 8]]) + assert M1234 * S(2) == M([[2, 4], [6, 8]]) + assert raises(lambda: M1234 * None, TypeError) + assert raises(lambda: None * M1234, TypeError) + assert raises(lambda: M1234 * M([[1, 2], [3, 4], [5, 6]]), ValueError) + assert raises(lambda: M([[1, 2, 3], [4, 5, 6]]) * M1234, ValueError) + + for M2 in _compatible_matrix_types(M): + assert M1234 * M2([[1, 2], [3, 4]]) == M([[7, 10], [15, 22]]) + assert M2([[1, 2], [3, 4]]) * M1234 == M([[7, 10], [15, 22]]) + + for M2 in _incompatible_matrix_types(M): + assert raises(lambda: M1234 * M2([[1, 2], [3, 4]]), (TypeError, ValueError)) + assert raises(lambda: M2([[1, 2], [3, 4]]) * M1234, (TypeError, ValueError)) + + +def test_matrices_pow(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + assert M1234**0 == M([[1, 0], [0, 1]]) + assert M1234**1 == M1234 + assert M1234**2 == M([[7, 10], [15, 22]]) + assert M1234**3 == M([[37, 54], [81, 118]]) + if is_field: + assert M1234**-1 == M([[-4, 2], [3, -1]]) / 2 + assert M1234**-2 == M([[22, -10], [-15, 7]]) / 4 + assert M1234**-3 == M([[-118, 54], [81, -37]]) / 8 + Ms = M([[1, 2], [3, 6]]) + assert raises(lambda: Ms**-1, ZeroDivisionError) + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert raises(lambda: Mr**0, ValueError) + assert raises(lambda: Mr**1, ValueError) + assert raises(lambda: Mr**2, ValueError) + assert raises(lambda: M1234**None, TypeError) + assert raises(lambda: None**M1234, TypeError) + + +def test_matrices_div(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + if is_field: + assert M1234 / 2 == M([[S(1)/2, S(1)], [S(3)/2, 2]]) + assert M1234 / S(2) == M([[S(1)/2, S(1)], [S(3)/2, 2]]) + assert raises(lambda: M1234 / 0, ZeroDivisionError) + assert raises(lambda: M1234 / S(0), ZeroDivisionError) + raises(lambda: M1234 / None, TypeError) + raises(lambda: None / M1234, TypeError) + + +def test_matrices_inv(): + for M, S, is_field in _all_matrices(): + if is_field: + M1234 = M([[1, 2], [3, 4]]) + assert M1234.inv() == M([[-2, 1], [S(3)/2, -S(1)/2]]) + M1236 = M([[1, 2], [3, 6]]) + assert raises(lambda: M1236.inv(), ZeroDivisionError) + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert raises(lambda: Mr.inv(), ValueError) + # XXX: Test non-field matrices. unimodular? + + +def test_matrices_det(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + assert M1234.det() == S(-2) + M9 = M([[1, 2, 3], [4, 5, 6], [7, 8, 10]]) + assert M9.det() == S(-3) + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert raises(lambda: Mr.det(), ValueError) + + +def test_matrices_charpoly(): + for M, S, is_field in _all_matrices(): + P = _poly_type_from_matrix_type(M) + M1234 = M([[1, 2], [3, 4]]) + assert M1234.charpoly() == P([-2, -5, 1]) + M9 = M([[1, 2, 3], [4, 5, 6], [7, 8, 10]]) + assert M9.charpoly() == P([3, -12, -16, 1]) + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert raises(lambda: Mr.charpoly(), ValueError) + + +def test_matrices_minpoly(): + for M, S, is_field in _all_matrices(): + P = _poly_type_from_matrix_type(M) + M1234 = M([[1, 2], [3, 4]]) + assert M1234.minpoly() == P([-2, -5, 1]) + M9 = M([[2, 1, 0], [0, 2, 0], [0, 0, 2]]) + assert M9.minpoly() == P([4, -4, 1]) + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert raises(lambda: Mr.minpoly(), ValueError) + + +def test_matrices_rank(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2], [3, 4]]) + assert M1234.rank() == 2 + Mr = M([[1, 2, 3], [4, 5, 6]]) + assert Mr.rank() == 2 + Ms = M([[1, 2], [2, 4]]) + assert Ms.rank() == 1 + Mz = M([[0, 0], [0, 0]]) + assert Mz.rank() == 0 + + +def test_matrices_rref(): + for M, S, is_field in _all_matrices(): + if is_field: + Mr = M([[1, 2, 3], [4, 5, 6]]) + Mr_rref = M([[1, 0, -1], [0, 1, 2]]) + assert Mr.rref() == (Mr_rref, 2) + assert Mr == M([[1, 2, 3], [4, 5, 6]]) + assert Mr.rref(inplace=True) == (Mr_rref, 2) + assert Mr == Mr_rref + + +def test_matrices_solve(): + for M, S, is_field in _all_matrices(): + if is_field: + A = M([[1, 2], [3, 4]]) + x = M([[1], [2]]) + b = M([[5], [11]]) + assert A*x == b + assert A.solve(b) == x + A22 = M([[1, 2], [3, 4]]) + A23 = M([[1, 2, 3], [4, 5, 6]]) + b2 = M([[5], [11]]) + b3 = M([[5], [11], [17]]) + assert raises(lambda: A22.solve(b3), ValueError) + assert raises(lambda: A23.solve(b2), ValueError) + assert raises(lambda: A.solve(None), TypeError) + A = M([[1, 2], [2, 4]]) + assert raises(lambda: A.solve(b), ZeroDivisionError) + + +def test_matrices_transpose(): + for M, S, is_field in _all_matrices(): + M1234 = M([[1, 2, 3], [4, 5, 6]]) + assert M1234.transpose() == M([[1, 4], [2, 5], [3, 6]]) + + def test_all_tests(): test_funcs = {f for name, f in globals().items() if name.startswith("test_")} untested = test_funcs - set(all_tests) assert not untested, f"Untested functions: {untested}" + all_tests = [ + test_pyflint, test_showgood, + test_fmpz, test_fmpz_factor, test_fmpz_functions, @@ -2441,19 +3007,48 @@ def test_all_tests(): test_fmpz_poly_functions, test_fmpz_mat, test_fmpz_series, + test_fmpq, test_fmpq_poly, test_fmpq_mat, test_fmpq_series, + test_nmod, test_nmod_poly, test_nmod_mat, test_nmod_series, - test_arb, + test_fmpz_mod, test_fmpz_mod_dlog, test_fmpz_mod_poly, + test_fmpz_mod_mat, + + test_arb, + test_polys, + + test_matrices_eq, + test_matrices_constructor, + test_matrices_strrepr, + test_matrices_getitem, + test_matrices_setitem, + test_matrices_bool, + test_matrices_transpose, + test_matrices_pos_neg, + test_matrices_add, + test_matrices_sub, + test_matrices_mul, + test_matrices_pow, + test_matrices_div, + test_matrices_inv, + test_matrices_det, + test_matrices_charpoly, + test_matrices_minpoly, + test_matrices_rank, + test_matrices_rref, + test_matrices_solve, + test_pickling, test_all_tests, + ] diff --git a/src/flint/types/fmpq_mat.pyx b/src/flint/types/fmpq_mat.pyx index bb322729..35d32ff3 100644 --- a/src/flint/types/fmpq_mat.pyx +++ b/src/flint/types/fmpq_mat.pyx @@ -90,7 +90,7 @@ cdef class fmpq_mat(flint_mat): x = fmpq(entries[i*n + j]) fmpq_set(fmpq_mat_entry(self.val, i, j), (x).val) else: - raise ValueError("fmpq_mat: expected 1-3 arguments") + raise TypeError("fmpq_mat: expected 1-3 arguments") def __nonzero__(self): return not fmpq_mat_is_zero(self.val) @@ -121,7 +121,7 @@ cdef class fmpq_mat(flint_mat): cdef fmpq x i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) x = fmpq.__new__(fmpq) fmpq_set(x.val, fmpq_mat_entry(self.val, i, j)) return x @@ -130,7 +130,7 @@ cdef class fmpq_mat(flint_mat): cdef long i, j i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) c = fmpq(value) # XXX fmpq_set(fmpq_mat_entry(self.val, i, j), (c).val) @@ -297,7 +297,7 @@ cdef class fmpq_mat(flint_mat): def transpose(self): """ Returns the transpose of *self*. - + >>> fmpq_mat(2,3,range(6)).transpose() [0, 3] [1, 4] @@ -393,6 +393,17 @@ cdef class fmpq_mat(flint_mat): rank = fmpq_mat_rref((res).val, self.val) return res, rank + def rank(self): + """ + Returns the rank of *self*. + + >>> from flint import fmpq_mat + >>> M = fmpq_mat([[1,2,3],[4,5,6],[7,8,9]]) + >>> M.rank() + 2 + """ + return self.rref()[1] + @classmethod def hilbert(cls, long n, long m): """ @@ -429,14 +440,46 @@ cdef class fmpq_mat(flint_mat): return num, den def charpoly(self): + """Returns the characteristic polynomial of *self* as an *fmpq_poly*. + + >>> from flint import fmpq_mat, fmpq + >>> A = fmpq_mat(3,3,[1,3,5,2,4,6,fmpq(2,3),2,4]) + >>> A + [ 1, 3, 5] + [ 2, 4, 6] + [2/3, 2, 4] + >>> A.charpoly() + x^3 + (-9)*x^2 + 8/3*x + 4/3 + """ cdef fmpq_poly u + + if self.nrows() != self.ncols(): + raise ValueError("matrix must be square") + u = fmpq_poly.__new__(fmpq_poly) fmpq_poly_init(u.val) fmpq_mat_charpoly(u.val, self.val) return u def minpoly(self): + """Returns the minimal polynomial of *self* as an *fmpq_poly*. + + >>> from flint import fmpq_mat, fmpq + >>> A = fmpq_mat([[2, 1, 0], [0, 2, 0], [0, 0, 2]]) + >>> A + [2, 1, 0] + [0, 2, 0] + [0, 0, 2] + >>> A.charpoly() + x^3 + (-6)*x^2 + 12*x + (-8) + >>> A.minpoly() + x^2 + (-4)*x + 4 + """ cdef fmpq_poly u + + if self.nrows() != self.ncols(): + raise ValueError("matrix must be square") + u = fmpq_poly.__new__(fmpq_poly) fmpq_poly_init(u.val) fmpq_mat_minpoly(u.val, self.val) @@ -444,12 +487,16 @@ cdef class fmpq_mat(flint_mat): def __pow__(self, n, z): cdef fmpq_mat v - assert z is None + + if self.nrows() != self.ncols(): + raise ValueError("matrix must be square") + if z is not None: + raise TypeError("fmpq_mat does not support modular exponentiation") + n = int(n) if n == 0: - r, c = self.nrows(), self.ncols() - assert r == c - v = fmpq_mat(r, c) + r = self.nrows() + v = fmpq_mat(r, r) fmpq_mat_one(v.val) return v if n == 1: @@ -458,6 +505,7 @@ cdef class fmpq_mat(flint_mat): return self * self if n < 0: return self.inv() ** (-n) + v = self ** (n // 2) v = v * v if n % 2: diff --git a/src/flint/types/fmpz_mat.pyx b/src/flint/types/fmpz_mat.pyx index 121f2ce8..a2609c2a 100644 --- a/src/flint/types/fmpz_mat.pyx +++ b/src/flint/types/fmpz_mat.pyx @@ -126,7 +126,7 @@ cdef class fmpz_mat(flint_mat): x = fmpz(entries[i*n + j]) fmpz_set(fmpz_mat_entry(self.val, i, j), (x).val) else: - raise ValueError("fmpz_mat: expected 1-3 arguments") + raise TypeError("fmpz_mat: expected 1-3 arguments") def __nonzero__(self): return not fmpz_mat_is_zero(self.val) @@ -163,7 +163,7 @@ cdef class fmpz_mat(flint_mat): cdef fmpz x i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) x = fmpz.__new__(fmpz) fmpz_set(x.val, fmpz_mat_entry(self.val, i, j)) return x @@ -172,7 +172,7 @@ cdef class fmpz_mat(flint_mat): cdef long i, j i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) c = fmpz(value) # XXX fmpz_set(fmpz_mat_entry(self.val, i, j), (c).val) @@ -716,16 +716,47 @@ cdef class fmpz_mat(flint_mat): return bool(fmpz_mat_is_in_snf(self.val)) def charpoly(self): + """Returns the characteristic polynomial of *self* as an *fmpz_poly*. + + >>> from flint import fmpz_mat + >>> A = fmpz_mat(3, 3, range(9)) + >>> A + [0, 1, 2] + [3, 4, 5] + [6, 7, 8] + >>> A.charpoly() + x^3 + (-12)*x^2 + (-18)*x + """ cdef fmpz_poly u + + if not fmpz_mat_is_square(self.val): + raise ValueError("matrix must be square") + u = fmpz_poly.__new__(fmpz_poly) fmpz_poly_init(u.val) fmpz_mat_charpoly(u.val, self.val) return u def minpoly(self): + """Returns the minimal polynomial of *self* as an *fmpz_poly*. + + >>> from flint import fmpz_mat + >>> A = fmpz_mat([[2, 1, 0], [0, 2, 0], [0, 0, 2]]) + >>> A + [2, 1, 0] + [0, 2, 0] + [0, 0, 2] + >>> A.charpoly() + x^3 + (-6)*x^2 + 12*x + (-8) + >>> A.minpoly() + x^2 + (-4)*x + 4 + """ cdef fmpz_poly u + + if not fmpz_mat_is_square(self.val): + raise ValueError("matrix must be square") + u = fmpz_poly.__new__(fmpz_poly) fmpz_poly_init(u.val) fmpz_mat_minpoly(u.val, self.val) return u - diff --git a/src/flint/types/fmpz_mod.pyx b/src/flint/types/fmpz_mod.pyx index 9a15807d..63c4c1e3 100644 --- a/src/flint/types/fmpz_mod.pyx +++ b/src/flint/types/fmpz_mod.pyx @@ -409,21 +409,27 @@ cdef class fmpz_mod(flint_scalar): def __bool__(self): return not self.is_zero() - def __repr__(self): + def repr(self): return "fmpz_mod({}, {})".format( fmpz_get_intlong(self.val), self.ctx.modulus() ) + def __repr__(self): + return self.repr() + def __hash__(self): return hash((int(self))) def __int__(self): return fmpz_get_intlong(self.val) - def __str__(self): + def str(self): return str(int(self)) + def __str__(self): + return self.str() + # ---------------- # # Arithmetic # # ---------------- # diff --git a/src/flint/types/fmpz_mod_mat.pxd b/src/flint/types/fmpz_mod_mat.pxd new file mode 100644 index 00000000..085a523d --- /dev/null +++ b/src/flint/types/fmpz_mod_mat.pxd @@ -0,0 +1,27 @@ +from flint.flintlib.flint cimport slong +from flint.flintlib.fmpz cimport fmpz_t +from flint.flintlib.fmpz_mod_mat cimport fmpz_mod_mat_t + +from flint.flint_base.flint_base cimport flint_mat +from flint.types.fmpz_mod cimport fmpz_mod_ctx, fmpz_mod + + +cdef class fmpz_mod_mat(flint_mat): + cdef fmpz_mod_mat_t val + cdef bint _initialized + cdef fmpz_mod_ctx ctx + + cdef void _init_empty(self, slong m, slong n, list args) + cdef void _init_empty_ctx(self, slong m, slong n, fmpz_mod_ctx ctx) + cdef void _init_from_list(self, slong m, slong n, list entries, list args) + #cdef void _init_from_matrix(self, flint_mat M, list args) + cdef void _init_from_matrix(self, M, list args) + cdef fmpz_mod_ctx _parse_args(self, list args) + cdef fmpz_mod_mat _new(self, slong m, slong n, fmpz_mod_ctx ctx) + cdef fmpz_mod_mat _newlike(self) + cdef fmpz_mod_mat _copy(self) + + cpdef slong nrows(self) + cpdef slong ncols(self) + cdef fmpz_mod _getitem(self, slong i, slong j) + cdef void _setitem(self, slong i, slong j, fmpz_t e) diff --git a/src/flint/types/fmpz_mod_mat.pyx b/src/flint/types/fmpz_mod_mat.pyx new file mode 100644 index 00000000..2635f3ac --- /dev/null +++ b/src/flint/types/fmpz_mod_mat.pyx @@ -0,0 +1,656 @@ +from flint.utils.typecheck cimport ( + typecheck, +) +from flint.flintlib.fmpz cimport ( + fmpz_struct, + fmpz_t, + fmpz_init_set, + fmpz_set, +) +from flint.flintlib.fmpz_mod_mat cimport ( + fmpz_mod_mat_init, + fmpz_mod_mat_init_set, + fmpz_mod_mat_clear, + fmpz_mod_mat_set, + fmpz_mod_mat_nrows, + fmpz_mod_mat_ncols, + fmpz_mod_mat_entry, + fmpz_mod_mat_set_entry, + fmpz_mod_mat_one, + fmpz_mod_mat_equal, + fmpz_mod_mat_is_zero, + fmpz_mod_mat_neg, + fmpz_mod_mat_add, + fmpz_mod_mat_sub, + fmpz_mod_mat_mul, + fmpz_mod_mat_scalar_mul_fmpz, + fmpz_mod_mat_inv, + fmpz_mod_mat_transpose, + fmpz_mod_mat_solve, + fmpz_mod_mat_rref, + fmpz_mod_mat_charpoly, + fmpz_mod_mat_minpoly, +) + +from flint.flint_base.flint_base cimport ( + flint_mat, +) +from flint.types.fmpz cimport ( + fmpz, +) +from flint.types.nmod cimport ( + nmod, +) +from flint.types.fmpz_mod cimport ( + fmpz_mod_ctx, + fmpz_mod, +) +from flint.types.fmpz_mod_poly cimport ( + fmpz_mod_poly, + fmpz_mod_poly_ctx, +) +from flint.types.fmpz_mat cimport ( + fmpz_mat, +) +from flint.types.nmod_mat cimport ( + nmod_mat, +) + + +cdef any_as_fmpz_mod_mat(x): + if typecheck(x, fmpz_mod_mat): + return x + elif typecheck(x, nmod_mat): + return fmpz_mod_mat(x) + else: + return NotImplemented + + +cdef any_as_fmpz_mod_mat_ctx(x, fmpz_mod_ctx ctx): + y = any_as_fmpz_mod_mat(x) + if y is not NotImplemented: + return y + if typecheck(x, fmpz_mat): + return fmpz_mod_mat(x, ctx) + return NotImplemented + + +cdef class fmpz_mod_mat(flint_mat): + """ + The ``fmpz_mod_mat`` type represents dense matrices over ``Z/nZ`` for + arbitrary ``n`` (not necessarily word-sized unlike ``nmod_mat``). Some + operations may assume that n is a prime. + """ + def __dealloc__(self): + if self._initialized: + fmpz_mod_mat_clear(self.val) + + def __init__(self, *args): + """Construct an ``fmpz_mod_mat`` matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx, fmpz_mat + >>> ctx = fmpz_mod_ctx(7) + >>> fmpz_mod_mat(2, 2, [1, 2, 3, 4], ctx) + [1, 2] + [3, 4] + >>> fmpz_mod_mat([[1, 2], [3, 4]], ctx) + [1, 2] + [3, 4] + >>> fmpz_mod_mat(2, 2, ctx) + [0, 0] + [0, 0] + >>> zmat = fmpz_mat(2, 2, [1, 2, 3, 4]) + >>> fmpz_mod_mat(zmat, ctx) + [1, 2] + [3, 4] + """ + # XXX: This logic should be moved to flint_mat. + if len(args) >= 2 and isinstance(args[0], int) and isinstance(args[1], int): + if len(args) >= 3 and isinstance(args[2], (list, tuple)): + m, n, entries, *arg = args + if len(entries) != m * n: + raise ValueError("fmpz_mod_mat: invalid number of entries") + self._init_from_list(m, n, entries, arg) + else: + m, n, *arg = args + self._init_empty(m, n, arg) + + elif (len(args) >= 1 and isinstance(args[0], (list, tuple)) + and all(isinstance(arg, (list, tuple)) for arg in args[0])): + lol, *arg = args + m = len(lol) + if m == 0: + n = 0 + entries = [] + else: + n = len(lol[0]) + if any(len(row) != n for row in lol): + raise ValueError("fmpz_mod_mat: inconsistent row lengths") + entries = [x for row in lol for x in row] + self._init_from_list(m, n, entries, arg) + + # XXX: nmod_mat is not a subclass of flint_mat (it should be) + elif len(args) >= 1 and isinstance(args[0], (flint_mat, nmod_mat)): + M, *arg = args + self._init_from_matrix(M, arg) + + else: + raise TypeError("fmpz_mod_mat: invalid arguments") + + cdef fmpz_mod_ctx _parse_args(self, list args): + """Parse the context argument.""" + if len(args) != 1: + if len(args) == 0: + raise TypeError("fmpz_mod_mat: missing modulus argument") + else: + raise TypeError("fmpz_mod_mat: too many arguments") + arg = args[0] + if not typecheck(arg, fmpz_mod_ctx): + raise TypeError("fmpz_mod_mat: invalid modulus argument") + return arg + + cdef void _init_empty_ctx(self, slong m, slong n, fmpz_mod_ctx ctx): + """Initialize an empty matrix with a given modulus context.""" + self.ctx = ctx + fmpz_mod_mat_init(self.val, m, n, ctx.val.n) + self._initialized = True + + cdef void _init_empty(self, slong m, slong n, list args): + """Initialize an empty matrix using the constructor arguments.""" + cdef fmpz_mod_ctx ctx + ctx = self._parse_args(args) + self._init_empty_ctx(m, n, ctx) + + cdef void _init_from_list(self, slong m, slong n, list entries, list args): + """Initialize a matrix from a list of entries.""" + cdef fmpz_mod_ctx ctx + cdef fmpz_mod e + ctx = self._parse_args(args) + entries = [ctx.any_as_fmpz_mod(x) for x in entries] + + self._init_empty_ctx(m, n, ctx) + for i in range(m): + for j in range(n): + val = ctx.any_as_fmpz_mod(entries[i*n + j]) + if val is NotImplemented: + raise TypeError("fmpz_mod_mat: incompatible entries") + e = val + fmpz_mod_mat_set_entry(self.val, i, j, e.val) + + # XXX: Should be possible to type this as flint_mat but nmod_mat is not a subclass + # cdef void _init_from_matrix(self, flint_mat M, list args): + cdef void _init_from_matrix(self, M, list args): + """Initialize from another matrix.""" + cdef fmpz_mod_ctx ctx + cdef fmpz_mod_mat N1 + + if typecheck(M, fmpz_mod_mat): + N1 = M + if args: + ctx = self._parse_args(args) + else: + ctx = N1.ctx + if N1.ctx != ctx: + raise TypeError("fmpz_mod_mat: incompatible moduli") + fmpz_mod_mat_init_set(self.val, N1.val) + self.ctx = ctx + self._initialized = True + elif typecheck(M, fmpz_mat): + # XXX: This is inefficient. + self._init_from_list(M.nrows(), M.ncols(), M.entries(), args) + elif typecheck(M, nmod_mat): + m = M.modulus() + if args: + ctx = self._parse_args(args) + if m != ctx.modulus(): + raise TypeError("fmpz_mod_mat: incompatible moduli") + else: + ctx = fmpz_mod_ctx(m) + # XXX: This is inefficient. + entries = [int(x) for x in M.entries()] + self._init_from_list(M.nrows(), M.ncols(), entries, [ctx]) + else: + raise TypeError("fmpz_mod_mat: unrocognized matrix type") + + cdef fmpz_mod_mat _new(self, slong m, slong n, fmpz_mod_ctx ctx): + """Create an initialized matrix of given shape and context.""" + cdef fmpz_mod_mat res + res = fmpz_mod_mat.__new__(fmpz_mod_mat) + fmpz_mod_mat_init(res.val, m, n, ctx.val.n) + res.ctx = ctx + res._initialized = True + return res + + cdef fmpz_mod_mat _newlike(self): + """Create an uninitialized matrix of the same shape and context.""" + return self._new(self.nrows(), self.ncols(), self.ctx) + + cdef fmpz_mod_mat _copy(self): + """Create a copy of the matrix.""" + cdef fmpz_mod_mat res + res = self._newlike() + fmpz_mod_mat_set(res.val, self.val) + return res + + cpdef slong nrows(self): + """Return the number of rows.""" + return fmpz_mod_mat_nrows(self.val) + + cpdef slong ncols(self): + """Return the number of columns.""" + return fmpz_mod_mat_ncols(self.val) + + def modulus(self): + """Return the modulus.""" + cdef fmpz mod + mod = fmpz.__new__(fmpz) + fmpz_init_set(mod.val, self.val.mod) + return mod + + cdef fmpz_mod _getitem(self, slong i, slong j): + """Retrieve an element as an ``fmpz_mod``.""" + cdef fmpz_struct * p_e + cdef fmpz_mod e + p_e = fmpz_mod_mat_entry(self.val, i, j) + e = fmpz_mod.__new__(fmpz_mod) + fmpz_set(e.val, p_e) + e.ctx = self.ctx + return e + + cdef void _setitem(self, slong i, slong j, fmpz_t e): + """Set an element from a raw ``fmpz_t``.""" + fmpz_mod_mat_set_entry(self.val, i, j, e) + + def repr(self): + """Return a representation string.""" + # XXX: Generic repr does not include modulus argument. + m = self.nrows() + n = self.ncols() + mod = self.modulus() + entries = (", ".join([x.repr() for x in self.entries()])) + return f"fmpz_mod_mat({m}, {n}, [{entries}], {mod})" + + def entries(self): + """Return a list of entries.""" + # XXX: Ideally move this to flint_mat + cdef slong i, j, m, n + m = self.nrows() + n = self.ncols() + L = [None] * (m * n) + for i from 0 <= i < m: + for j from 0 <= j < n: + L[i*n + j] = self._getitem(i, j) + return L + + def __getitem__(self, index): + """Retrieve an element.""" + # XXX: Ideally move this to flint_mat + cdef slong i, j + i, j = index + if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) + return self._getitem(i, j) + + def __setitem__(self, index, value): + """Set an element.""" + cdef slong i, j + cdef fmpz_mod e + if not (isinstance(index, tuple) and len(index) == 2): + raise TypeError("fmpz_mod_mat indices must be 2-tuples") + i, j = index + if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) + e = self.ctx.any_as_fmpz_mod(value) + self._setitem(i, j, e.val) + + def __nonzero__(self): + """Return ``True`` if the matrix has any nonzero entries.""" + cdef bint zero + zero = fmpz_mod_mat_is_zero(self.val) + return not zero + + def __richcmp__(self, other, int op): + """Compare two matrices.""" + cdef bint res + + if op != 2 and op != 3: + raise TypeError("matrices cannot be ordered") + + other = any_as_fmpz_mod_mat_ctx(other, self.ctx) + if other is NotImplemented: + return other + + res = fmpz_mod_mat_equal((self).val, (other).val) + + if op == 2: + return res + if op == 3: + return not res + + def __pos__(self): + """``+M``: Return self (not a copy).""" + return self + + def __neg__(self): + """``-M``""" + res = self._newlike() + fmpz_mod_mat_neg(( res).val, self.val) + return res + + def _as_fmpz_mod_mat(self, other, same_shape=True): + """Convert to ``fmpz_mod_mat`` but raise if shape or modulus do not match.""" + cdef fmpz_mod_mat o + other = any_as_fmpz_mod_mat_ctx(other, self.ctx) + if other is NotImplemented: + return NotImplemented + o = other + if same_shape and not (self.nrows() == o.nrows() and self.ncols() == o.ncols()): + raise ValueError("Shape mismatch: cannot add matrices") + if self.ctx != o.ctx: + raise ValueError("fmpz_mod_mat: incompatible moduli") + return o + + def _add(self, fmpz_mod_mat other): + """Add two ``fmpz_mod_mat`` matrices.""" + res = self._newlike() + fmpz_mod_mat_add(res.val, self.val, other.val) + return res + + def _sub(self, fmpz_mod_mat other): + """Subtract two ``fmpz_mod_mat`` matrices.""" + res = self._newlike() + fmpz_mod_mat_sub(res.val, self.val, other.val) + return res + + def _matmul(self, fmpz_mod_mat other): + """Multiply two ``fmpz_mod_mat`` matrices.""" + if self.ncols() != other.nrows(): + raise ValueError("Shape mismatch: cannot multiply matrices") + res = self._new(self.nrows(), other.ncols(), self.ctx) + fmpz_mod_mat_mul(res.val, self.val, other.val) + return res + + def _scalarmul(self, fmpz_mod other): + """Multiply an ``fmpz_mod_mat`` matrix by an ``fmpz_mod`` scalar.""" + res = self._newlike() + fmpz_mod_mat_scalar_mul_fmpz(res.val, self.val, other.val) + return res + + def _pow(self, slong other): + """Raise an ``fmpz_mod_mat`` matrix to an integer power.""" + cdef fmpz_mod_mat res, tmp + + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat pow: matrix must be square") + if other < 0: + self = self.inv() + other = -other + + res = self._newlike() + fmpz_mod_mat_one(res.val) + + tmp = self._copy() + + while other > 0: + if other % 2 == 1: + fmpz_mod_mat_mul(res.val, res.val, tmp.val) + fmpz_mod_mat_mul(tmp.val, tmp.val, tmp.val) + other //= 2 + + return res + + def _div(self, fmpz_mod other): + """Divide an ``fmpz_mod_mat`` matrix by an ``fmpz_mod`` scalar.""" + return self._scalarmul(other.inverse()) + + def __add__(self, other): + """``M + N``: Add two matrices.""" + cdef fmpz_mod_mat o + other = self._as_fmpz_mod_mat(other) + if other is NotImplemented: + return NotImplemented + o = other + return self._add(other) + + def __radd__(self, other): + """``N + M``: Add two matrices.""" + cdef fmpz_mod_mat o + other = self._as_fmpz_mod_mat(other) + if other is NotImplemented: + return NotImplemented + o = other + return other._add(self) + + def __sub__(self, other): + """``M - N``: Subtract two matrices.""" + cdef fmpz_mod_mat o + other = self._as_fmpz_mod_mat(other) + if other is NotImplemented: + return NotImplemented + o = other + return self._sub(other) + + def __rsub__(self, other): + """``N - M``: Subtract two matrices.""" + cdef fmpz_mod_mat o + other = self._as_fmpz_mod_mat(other) + if other is NotImplemented: + return NotImplemented + o = other + return other._sub(self) + + def __mul__(self, other): + """``M * N``: Multiply two matrices.""" + cdef fmpz_mod_mat o + cdef fmpz_mod e + other_mat = self._as_fmpz_mod_mat(other, same_shape=False) + if other_mat is not NotImplemented: + o = other_mat + return self._matmul(other_mat) + other_scalar = self.ctx.any_as_fmpz_mod(other) + if other_scalar is not NotImplemented: + e = other_scalar + return self._scalarmul(e) + return NotImplemented + + def __rmul__(self, other): + """``N * M``: Multiply two matrices.""" + cdef fmpz_mod e + other_mat = self._as_fmpz_mod_mat(other) + if other_mat is not NotImplemented: + o = other_mat + return other_mat._matmul(self) + other_scalar = self.ctx.any_as_fmpz_mod(other) + if other_scalar is not NotImplemented: + e = other_scalar + return self._scalarmul(e) + return NotImplemented + + def __pow__(self, other): + """``M ** n``: Raise a matrix to an integer power.""" + if not isinstance(other, int): + return NotImplemented + return self._pow(other) + + def __truediv__(self, other): + """``M / n``: Divide a matrix by a scalar.""" + cdef fmpz_mod e + other_scalar = self.ctx.any_as_fmpz_mod(other) + if other_scalar is not NotImplemented: + e = other_scalar + return self._div(e) + return NotImplemented + + def inv(self): + """Return the inverse of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2], [3, 4]], ctx) + >>> M.inv() + [5, 1] + [5, 3] + + Assumes that the modulus is prime. + """ + cdef fmpz_mod_mat res + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat inv: matrix must be square") + res = self._newlike() + r = fmpz_mod_mat_inv(res.val, self.val) + if r == 0: + raise ZeroDivisionError("fmpz_mod_mat inv: matrix is not invertible") + return res + + def det(self): + """Return the determinant of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2], [3, 4]], ctx) + >>> M.det() + fmpz_mod(5, 7) + + """ + # XXX: No fmpz_mod_mat_det function... + p = self.charpoly() + p0 = p[0] + if self.nrows() % 2: + p0 = -p0 + return p0 + + def charpoly(self): + """Return the characteristic polynomial of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2], [3, 4]], ctx) + >>> M.charpoly() + x^2 + 2*x + 5 + + """ + cdef fmpz_mod_poly_ctx pctx + cdef fmpz_mod_poly res + + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat charpoly: matrix must be square") + + pctx = fmpz_mod_poly_ctx(self.ctx) + res = fmpz_mod_poly(0, pctx) + fmpz_mod_mat_charpoly(res.val, self.val, self.ctx.val) + return res + + def minpoly(self): + """Return the minimal polynomial of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[2, 1, 0], [0, 2, 0], [0, 0, 2]], ctx) + >>> M.charpoly() + x^3 + x^2 + 5*x + 6 + >>> M.minpoly() + x^2 + 3*x + 4 + + """ + cdef fmpz_mod_poly_ctx pctx + cdef fmpz_mod_poly res + + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat minpoly: matrix must be square") + + pctx = fmpz_mod_poly_ctx(self.ctx) + res = fmpz_mod_poly(0, pctx) + fmpz_mod_mat_minpoly(res.val, self.val, self.ctx.val) + return res + + def transpose(self): + """Return the transpose of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2], [3, 4]], ctx) + >>> M + [1, 2] + [3, 4] + >>> M.transpose() + [1, 3] + [2, 4] + """ + cdef fmpz_mod_mat res + res = self._new(self.ncols(), self.nrows(), self.ctx) + fmpz_mod_mat_transpose(res.val, self.val) + return res + + def solve(self, rhs): + """Solve a linear system. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2], [3, 4]], ctx) + >>> rhs = fmpz_mod_mat([[5], [6]], ctx) + >>> M.solve(rhs) + [3] + [1] + + The matrix must be square and invertible. + + Assumes that the modulus is prime. + """ + cdef bint success + cdef fmpz_mod_mat res + + rhs = any_as_fmpz_mod_mat(rhs) + + if rhs is NotImplemented: + raise TypeError("fmpz_mod_mat solve: invalid rhs") + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat solve: matrix must be square") + if self.nrows() != rhs.nrows(): + raise ValueError("fmpz_mod_mat solve: shape mismatch") + + res = self._new(rhs.nrows(), rhs.ncols(), self.ctx) + success = fmpz_mod_mat_solve(res.val, self.val, ( rhs).val) + + if not success: + raise ZeroDivisionError("Singular matrix in solve") + + return res + + def rank(self): + """Return the rank of a matrix. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(11) + >>> M = fmpz_mod_mat([[1, 2, 3], [4, 5, 6], [7, 8, 9]], ctx) + >>> M.rank() + 2 + + Assumes that the modulus is prime. + """ + return self.rref()[1] + + def rref(self, inplace=False): + """Return the reduced row echelon form of a matrix and the rank. + + >>> from flint import fmpz_mod_mat, fmpz_mod_ctx + >>> ctx = fmpz_mod_ctx(7) + >>> M = fmpz_mod_mat([[1, 2, 3], [4, 5, 6]], ctx) + >>> Mr, r = M.rref() + >>> Mr + [1, 0, 6] + [0, 1, 2] + >>> r + 2 + + If ``inplace`` is ``True``, the matrix is modified in place. + + Assumes that the modulus is prime. + """ + cdef fmpz_mod_mat res + cdef slong r + if inplace: + res = self + else: + res = self._copy() + r = fmpz_mod_mat_rref(NULL, res.val) + return (res, r) diff --git a/src/flint/types/nmod_mat.pxd b/src/flint/types/nmod_mat.pxd index fc930383..0b8336d1 100644 --- a/src/flint/types/nmod_mat.pxd +++ b/src/flint/types/nmod_mat.pxd @@ -3,7 +3,7 @@ from flint.flint_base.flint_base cimport flint_mat from flint.flintlib.nmod_mat cimport nmod_mat_t from flint.flintlib.flint cimport mp_limb_t -cdef class nmod_mat: +cdef class nmod_mat(flint_mat): cdef nmod_mat_t val cpdef long nrows(self) cpdef long ncols(self) diff --git a/src/flint/types/nmod_mat.pyx b/src/flint/types/nmod_mat.pyx index 350bf73d..35e016cb 100644 --- a/src/flint/types/nmod_mat.pyx +++ b/src/flint/types/nmod_mat.pyx @@ -1,20 +1,60 @@ -from flint.utils.conversion cimport matrix_to_str +cimport cython + +from flint.flintlib.flint cimport ulong, mp_limb_t +from flint.flintlib.nmod cimport nmod_t + +from flint.flintlib.nmod_poly cimport ( + nmod_poly_init, +) + +from flint.flintlib.fmpz_mat cimport fmpz_mat_nrows, fmpz_mat_ncols +from flint.flintlib.fmpz_mat cimport fmpz_mat_get_nmod_mat + +from flint.flintlib.nmod_mat cimport ( + nmod_mat_struct, + nmod_mat_init, + nmod_mat_init_set, + nmod_mat_clear, + nmod_mat_nrows, + nmod_mat_ncols, + nmod_mat_is_square, + nmod_mat_entry, + nmod_mat_set_entry, + nmod_mat_equal, + nmod_mat_is_zero, + nmod_mat_transpose, + nmod_mat_scalar_mul, + nmod_mat_neg, + nmod_mat_add, + nmod_mat_sub, + nmod_mat_mul, + nmod_mat_pow, + nmod_mat_inv, + nmod_mat_solve, + nmod_mat_nullspace, + nmod_mat_det, + nmod_mat_rref, + nmod_mat_rank, + nmod_mat_charpoly, + nmod_mat_minpoly, + nmod_mat_randtest, +) + from flint.utils.typecheck cimport typecheck from flint.types.fmpz_mat cimport any_as_fmpz_mat from flint.types.fmpz_mat cimport fmpz_mat from flint.types.nmod cimport nmod from flint.types.nmod cimport any_as_nmod +from flint.types.nmod_poly cimport nmod_poly from flint.pyflint cimport global_random_state from flint.flint_base.flint_context cimport thectx -cimport cython +from flint.flint_base.flint_base cimport flint_mat + -from flint.flintlib.flint cimport ulong -from flint.flintlib.nmod_mat cimport * -from flint.flintlib.fmpz_mat cimport fmpz_mat_nrows, fmpz_mat_ncols -from flint.flintlib.fmpz_mat cimport fmpz_mat_get_nmod_mat ctx = thectx + cdef any_as_nmod_mat(obj, nmod_t mod): cdef nmod_mat r cdef mp_limb_t v @@ -29,10 +69,13 @@ cdef any_as_nmod_mat(obj, nmod_t mod): return r return NotImplemented -cdef class nmod_mat: + +cdef class nmod_mat(flint_mat): """ - The nmod_poly type represents dense matrices over Z/nZ for - word-size n. Some operations may assume that n is a prime. + The nmod_mat type represents dense matrices over Z/nZ for word-size n (see + fmpz_mod_mat for larger moduli). + + Some operations may assume that n is a prime. """ # cdef nmod_mat_t val @@ -91,7 +134,7 @@ cdef class nmod_mat: x = nmod(entries[i*n + j], mod) # XXX: slow self.val.rows[i][j] = (x).val else: - raise ValueError("nmod_mat: expected 1-3 arguments plus modulus") + raise TypeError("nmod_mat: expected 1-3 arguments plus modulus") def __nonzero__(self): return not nmod_mat_is_zero(self.val) @@ -131,14 +174,11 @@ cdef class nmod_mat: nmod_mat_randtest(mat.val, global_random_state) return mat - def __repr__(self): - if ctx.pretty: - return str(self) - return "nmod_mat(%i, %i, %s, %i)" % (self.nrows(), self.ncols(), - [int(c) for c in self.entries()], self.modulus()) - - def __str__(self): - return matrix_to_str(self.table()) + def repr(self): + m = self.nrows() + n = self.ncols() + entries = ', '.join(map(str, self.entries())) + return f"nmod_mat({m}, {n}, [{entries}], {self.modulus()})" def entries(self): cdef long i, j, m, n @@ -153,19 +193,12 @@ cdef class nmod_mat: L[i*n + j] = t return L - def table(self): - cdef long i, m, n - m = self.nrows() - n = self.ncols() - L = self.entries() - return [L[i*n:(i+1)*n] for i in range(m)] - def __getitem__(self, index): cdef long i, j cdef nmod x i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) x = nmod(nmod_mat_entry(self.val, i, j), self.modulus()) # XXX: slow return x @@ -174,26 +207,11 @@ cdef class nmod_mat: cdef mp_limb_t v i, j = index if i < 0 or i >= self.nrows() or j < 0 or j >= self.ncols(): - raise ValueError("index %i,%i exceeds matrix dimensions" % (i, j)) + raise IndexError("index %i,%i exceeds matrix dimensions" % (i, j)) if any_as_nmod(&v, value, self.val.mod): - self.val.rows[i][j] = v + nmod_mat_set_entry(self.val, i, j, v) else: - raise ValueError("cannot set item of type %s" % type(value)) - - def det(self): - """ - Returns the determinant of self as an nmod. - - >>> nmod_mat(2,2,[1,2,3,4],17).det() - 15 - - """ - if not nmod_mat_is_square(self.val): - raise ValueError("matrix must be square") - return nmod(nmod_mat_det(self.val), self.modulus()) - - def rank(self): - return nmod_mat_rank(self.val) + raise TypeError("cannot set item of type %s" % type(value)) def __pos__(self): return self @@ -313,6 +331,21 @@ cdef class nmod_mat: return u return u * s + def __pow__(self, e, m): + cdef nmod_mat t + cdef ulong ee + if not self.nrows() == self.ncols(): + raise ValueError("matrix must be square") + if m is not None: + raise NotImplementedError("modular matrix exponentiation") + if e < 0: + self = self.inv() + e = -e + ee = e + t = nmod_mat(self) # XXX + nmod_mat_pow(t.val, t.val, ee) + return t + @staticmethod def _div_(nmod_mat s, t): cdef mp_limb_t v @@ -327,7 +360,29 @@ cdef class nmod_mat: def __div__(s, t): return nmod_mat._div_(s, t) + def det(self): + """ + Returns the determinant of self as an nmod. + + >>> nmod_mat(2,2,[1,2,3,4],17).det() + 15 + + """ + if not nmod_mat_is_square(self.val): + raise ValueError("matrix must be square") + return nmod(nmod_mat_det(self.val), self.modulus()) + def inv(self): + """ + Returns the inverse of self. + + >>> from flint import nmod_mat + >>> A = nmod_mat(2,2,[1,2,3,4],17) + >>> A.inv() + [15, 1] + [10, 8] + + """ cdef nmod_mat u if not nmod_mat_is_square(self.val): raise ValueError("matrix must be square") @@ -422,6 +477,16 @@ cdef class nmod_mat: rank = nmod_mat_rref((res).val) return res, rank + def rank(self): + """Return the rank of a matrix. + + >>> from flint import nmod_mat + >>> M = nmod_mat([[1, 2], [3, 4]], 11) + >>> M.rank() + 2 + """ + return nmod_mat_rank(self.val) + def nullspace(self): """ Computes a basis for the nullspace of self. Returns (X, nullity) @@ -451,3 +516,42 @@ cdef class nmod_mat: nullity = nmod_mat_nullspace(res.val, self.val) return res, nullity + def charpoly(self): + """Return the characteristic polynomial of a matrix. + + >>> from flint import nmod_mat + >>> M = nmod_mat([[1, 2], [3, 4]], 11) + >>> M.charpoly() + x^2 + 6*x + 9 + + """ + cdef nmod_poly res + + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat charpoly: matrix must be square") + + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init(res.val, self.val.mod.n) + nmod_mat_charpoly(res.val, self.val) + return res + + def minpoly(self): + """Return the minimal polynomial of a matrix. + + >>> from flint import nmod_mat + >>> M = nmod_mat([[2, 1, 0], [0, 2, 0], [0, 0, 2]], 7) + >>> M.charpoly() + x^3 + x^2 + 5*x + 6 + >>> M.minpoly() + x^2 + 3*x + 4 + + """ + cdef nmod_poly res + + if self.nrows() != self.ncols(): + raise ValueError("fmpz_mod_mat minpoly: matrix must be square") + + res = nmod_poly.__new__(nmod_poly) + nmod_poly_init(res.val, self.val.mod.n) + nmod_mat_minpoly(res.val, self.val) + return res diff --git a/src/flint/utils/conversion.pxd b/src/flint/utils/conversion.pxd index 9641969b..2ae46001 100644 --- a/src/flint/utils/conversion.pxd +++ b/src/flint/utils/conversion.pxd @@ -1,5 +1,3 @@ -from cpython.version cimport PY_MAJOR_VERSION - cdef inline long prec_to_dps(n): return max(1, int(round(int(n)/3.3219280948873626)-1)) @@ -7,29 +5,10 @@ cdef inline long dps_to_prec(n): return max(1, int(round((int(n)+1)*3.3219280948873626))) cdef inline chars_from_str(s): - if PY_MAJOR_VERSION < 3: - return s - else: - return s.encode('ascii') + return s.encode('ascii') cdef inline str_from_chars(s): - if PY_MAJOR_VERSION < 3: - return str(s) - else: - return bytes(s).decode('ascii') - -cdef inline matrix_to_str(tab): - if len(tab) == 0 or len(tab[0]) == 0: - return "[]" - tab = [[str(c) for c in row] for row in tab] - widths = [] - for i in xrange(len(tab[0])): - w = max([len(row[i]) for row in tab]) - widths.append(w) - for i in xrange(len(tab)): - tab[i] = [s.rjust(widths[j]) for j, s in enumerate(tab[i])] - tab[i] = "[" + (", ".join(tab[i])) + "]" - return "\n".join(tab) + return bytes(s).decode('ascii') cdef inline _str_trunc(s, trunc=0): if trunc > 0 and len(s) > 3 * trunc: @@ -37,4 +16,3 @@ cdef inline _str_trunc(s, trunc=0): omitted = len(s) - left - right return s[:left] + ("{...%s digits...}" % omitted) + s[-right:] return s -