diff --git a/src/doc/en/reference/matrices/index.rst b/src/doc/en/reference/matrices/index.rst index fb06992c1e1..89453635472 100644 --- a/src/doc/en/reference/matrices/index.rst +++ b/src/doc/en/reference/matrices/index.rst @@ -94,6 +94,8 @@ objects like operation tables (e.g. the multiplication table of a group). sage/matrix/matrix_misc sage/matrix/matrix_window sage/matrix/misc + sage/matrix/misc_mpfr + sage/matrix/misc_flint sage/matrix/symplectic_basis sage/matrix/compute_J_ideal diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index aa865ee08c3..d5402d5c3b0 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -14883,7 +14883,6 @@ cdef class Matrix(Matrix1): 12215 """ from sage.rings.real_double import RDF - from sage.rings.real_mpfr import RealField try: A = self.change_ring(RDF) m1 = A._hadamard_row_bound() @@ -14891,13 +14890,14 @@ cdef class Matrix(Matrix1): m2 = A._hadamard_row_bound() return min(m1, m2) except (OverflowError, TypeError): + from sage.rings.real_mpfr import RealField # Try using MPFR, which handles large numbers much better, but is slower. - from . import misc + from .misc_mpfr import hadamard_row_bound_mpfr R = RealField(53, rnd='RNDU') A = self.change_ring(R) - m1 = misc.hadamard_row_bound_mpfr(A) + m1 = hadamard_row_bound_mpfr(A) A = A.transpose() - m2 = misc.hadamard_row_bound_mpfr(A) + m2 = hadamard_row_bound_mpfr(A) return min(m1, m2) def find(self,f, indices=False): diff --git a/src/sage/matrix/matrix_cyclo_dense.pyx b/src/sage/matrix/matrix_cyclo_dense.pyx index 363fe637e93..0b78bb7fd78 100644 --- a/src/sage/matrix/matrix_cyclo_dense.pyx +++ b/src/sage/matrix/matrix_cyclo_dense.pyx @@ -63,7 +63,7 @@ from .matrix cimport Matrix from . import matrix_dense from .matrix_integer_dense cimport _lift_crt from sage.structure.element cimport Matrix as baseMatrix -from .misc import matrix_integer_dense_rational_reconstruction +from .misc_flint import matrix_integer_dense_rational_reconstruction from sage.arith.misc import binomial, previous_prime from sage.rings.rational_field import QQ diff --git a/src/sage/matrix/matrix_integer_dense.pyx b/src/sage/matrix/matrix_integer_dense.pyx index a2c8f46247f..89bd3075db8 100644 --- a/src/sage/matrix/matrix_integer_dense.pyx +++ b/src/sage/matrix/matrix_integer_dense.pyx @@ -3410,7 +3410,7 @@ cdef class Matrix_integer_dense(Matrix_dense): ... ZeroDivisionError: The modulus cannot be zero """ - from .misc import matrix_integer_dense_rational_reconstruction + from .misc_flint import matrix_integer_dense_rational_reconstruction return matrix_integer_dense_rational_reconstruction(self, N) def randomize(self, density=1, x=None, y=None, distribution=None, diff --git a/src/sage/matrix/matrix_integer_sparse.pyx b/src/sage/matrix/matrix_integer_sparse.pyx index e8f374d0e71..49ed3403fca 100644 --- a/src/sage/matrix/matrix_integer_sparse.pyx +++ b/src/sage/matrix/matrix_integer_sparse.pyx @@ -399,13 +399,15 @@ cdef class Matrix_integer_sparse(Matrix_sparse): def rational_reconstruction(self, N): """ - Use rational reconstruction to lift self to a matrix over the - rational numbers (if possible), where we view self as a matrix - modulo N. + Use rational reconstruction to lift ``self`` to a matrix over the + rational numbers (if possible), where we view ``self`` as a matrix + modulo `N`. EXAMPLES:: - sage: A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, 7, 2, 2, 3, 4, 3, 4, (5/7)%500], sparse=True) + sage: A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, + ....: 7, 2, 2, 3, + ....: 4, 3, 4, (5/7)%500], sparse=True) sage: A.rational_reconstruction(500) [1/3 2 3 -4] [ 7 2 2 3] @@ -415,7 +417,7 @@ cdef class Matrix_integer_sparse(Matrix_sparse): Check that :trac:`9345` is fixed:: - sage: A = random_matrix(ZZ, 3, 3, sparse = True) + sage: A = random_matrix(ZZ, 3, 3, sparse=True) sage: A.rational_reconstruction(0) Traceback (most recent call last): ... diff --git a/src/sage/matrix/misc.pyx b/src/sage/matrix/misc.pyx index 8d5bd2c9b4c..1819d45591e 100644 --- a/src/sage/matrix/misc.pyx +++ b/src/sage/matrix/misc.pyx @@ -1,143 +1,46 @@ """ Misc matrix algorithms - -Code goes here mainly when it needs access to the internal structure -of several classes, and we want to avoid circular cimports. - -NOTE: The whole problem of avoiding circular imports -- the reason for -existence of this file -- is now a non-issue, since some bugs in -Cython were fixed. Probably all this code should be moved into the -relevant classes and this file deleted. """ from cysignals.signals cimport sig_check -cimport sage.rings.abc - from sage.arith.misc import CRT_basis, previous_prime from sage.arith.rational_reconstruction cimport mpq_rational_reconstruction from sage.data_structures.binary_search cimport * from sage.ext.mod_int cimport * -from sage.libs.flint.fmpq cimport fmpq_set_mpq, fmpq_canonicalise -from sage.libs.flint.fmpq_mat cimport fmpq_mat_entry_num, fmpq_mat_entry_den, fmpq_mat_entry -from sage.libs.flint.fmpz cimport fmpz_set_mpz, fmpz_one from sage.libs.gmp.mpq cimport * from sage.libs.gmp.mpz cimport * -from sage.libs.mpfr cimport * +from sage.misc.lazy_import import LazyImport from sage.misc.verbose import verbose from sage.modules.vector_integer_sparse cimport * from sage.modules.vector_modn_sparse cimport * from sage.modules.vector_rational_sparse cimport * from sage.rings.integer cimport Integer from sage.rings.rational_field import QQ -from sage.rings.real_mpfr cimport RealNumber from .matrix0 cimport Matrix -from .matrix_integer_dense cimport Matrix_integer_dense from .matrix_integer_sparse cimport Matrix_integer_sparse -from .matrix_rational_dense cimport Matrix_rational_dense from .matrix_rational_sparse cimport Matrix_rational_sparse - -def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer N): - """ - Given a matrix over the integers and an integer modulus, do - rational reconstruction on all entries of the matrix, viewed as - numbers mod N. This is done efficiently by assuming there is a - large common factor dividing the denominators. - - INPUT: - - A -- matrix - N -- an integer - - EXAMPLES:: - - sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ) - sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(B, 500) - [ 1/3 2/3 1 -4/3] - [ 7/3 2/3 6 1] - [ 4/3 1 4/3 5/3] - - TESTS: - - Check that :trac:`9345` is fixed:: - - sage: A = random_matrix(ZZ, 3) - sage: sage.matrix.misc.matrix_integer_dense_rational_reconstruction(A, 0) - Traceback (most recent call last): - ... - ZeroDivisionError: The modulus cannot be zero - """ - if not N: - raise ZeroDivisionError("The modulus cannot be zero") - cdef Matrix_rational_dense R - R = Matrix_rational_dense.__new__(Matrix_rational_dense, - A.parent().change_ring(QQ), 0,0,0) - - cdef mpz_t a, bnd, other_bnd, denom, tmp - cdef mpq_t qtmp - cdef Integer _bnd - cdef Py_ssize_t i, j - cdef int do_it - - mpz_init_set_si(denom, 1) - mpz_init(a) - mpz_init(tmp) - mpz_init(other_bnd) - mpq_init(qtmp) - - _bnd = (N//2).isqrt() - mpz_init_set(bnd, _bnd.value) - mpz_sub(other_bnd, N.value, bnd) - - for i in range(A._nrows): - for j in range(A._ncols): - sig_check() - A.get_unsafe_mpz(i, j, a) - if mpz_cmp_ui(denom, 1) != 0: - mpz_mul(a, a, denom) - mpz_fdiv_r(a, a, N.value) - do_it = 0 - if mpz_cmp(a, bnd) <= 0: - do_it = 1 - elif mpz_cmp(a, other_bnd) >= 0: - mpz_sub(a, a, N.value) - do_it = 1 - if do_it: - fmpz_set_mpz(fmpq_mat_entry_num(R._matrix, i, j), a) - if mpz_cmp_ui(denom, 1) != 0: - fmpz_set_mpz(fmpq_mat_entry_den(R._matrix, i, j), denom) - fmpq_canonicalise(fmpq_mat_entry(R._matrix, i, j)) - else: - fmpz_one(fmpq_mat_entry_den(R._matrix, i, j)) - else: - # Otherwise have to do it the hard way - A.get_unsafe_mpz(i, j, tmp) - mpq_rational_reconstruction(qtmp, tmp, N.value) - mpz_lcm(denom, denom, mpq_denref(qtmp)) - fmpq_set_mpq(fmpq_mat_entry(R._matrix, i, j), qtmp) - - mpz_clear(denom) - mpz_clear(a) - mpz_clear(tmp) - mpz_clear(other_bnd) - mpz_clear(bnd) - mpq_clear(qtmp) - - return R +matrix_integer_dense_rational_reconstruction = \ + LazyImport('sage.matrix.misc_flint', 'matrix_integer_dense_rational_reconstruction', + deprecation=35758) +hadamard_row_bound_mpfr = \ + LazyImport('sage.matrix.misc_mpfr', 'hadamard_row_bound_mpfr', + deprecation=35758) def matrix_integer_sparse_rational_reconstruction(Matrix_integer_sparse A, Integer N): - """ + r""" Given a sparse matrix over the integers and an integer modulus, do rational reconstruction on all entries of the matrix, viewed as - numbers mod N. + numbers mod `N`. EXAMPLES:: sage: A = matrix(ZZ, 3, 4, [(1/3)%500, 2, 3, (-4)%500, 7, 2, 2, 3, 4, 3, 4, (5/7)%500], sparse=True) - sage: sage.matrix.misc.matrix_integer_sparse_rational_reconstruction(A, 500) + sage: from sage.matrix.misc import matrix_integer_sparse_rational_reconstruction + sage: matrix_integer_sparse_rational_reconstruction(A, 500) [1/3 2 3 -4] [ 7 2 2 3] [ 4 3 4 5/7] @@ -424,32 +327,33 @@ def matrix_rational_echelon_form_multimodular(Matrix self, height_guess=None, pr def cmp_pivots(x, y): - """ + r""" Compare two sequences of pivot columns. - If x is shorter than y, return -1, i.e., x < y, "not as good". - If x is longer than y, then x > y, so "better" and return +1. - If the length is the same, then x is better, i.e., x > y - if the entries of x are correspondingly <= those of y with + If `x` is shorter than `y`, return `-1`, i.e., `x < y`, "not as good". + If `x` is longer than `y`, then `x > y`, so "better" and return `+1`. + If the length is the same, then `x` is better, i.e., `x > y` + if the entries of `x` are correspondingly `\leq` those of `y` with one being strictly less. INPUT: - - x, y -- lists or tuples of integers + - ``x``, ``y`` -- lists or tuples of integers EXAMPLES: We illustrate each of the above comparisons. :: - sage: sage.matrix.misc.cmp_pivots([1,2,3], [4,5,6,7]) + sage: from sage.matrix.misc import cmp_pivots + sage: cmp_pivots([1,2,3], [4,5,6,7]) -1 - sage: sage.matrix.misc.cmp_pivots([1,2,3,5], [4,5,6]) + sage: cmp_pivots([1,2,3,5], [4,5,6]) 1 - sage: sage.matrix.misc.cmp_pivots([1,2,4], [1,2,3]) + sage: cmp_pivots([1,2,4], [1,2,3]) -1 - sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,3]) + sage: cmp_pivots([1,2,3], [1,2,3]) 0 - sage: sage.matrix.misc.cmp_pivots([1,2,3], [1,2,4]) + sage: cmp_pivots([1,2,3], [1,2,4]) 1 """ x = tuple(x) @@ -464,70 +368,3 @@ def cmp_pivots(x, y): return 0 else: return -1 - - -def hadamard_row_bound_mpfr(Matrix A): - """ - Given a matrix A with entries that coerce to RR, compute the row - Hadamard bound on the determinant. - - INPUT: - - A -- a matrix over RR - - OUTPUT: - - integer -- an integer n such that the absolute value of the - determinant of this matrix is at most $10^n$. - - EXAMPLES: - - We create a very large matrix, compute the row Hadamard bound, - and also compute the row Hadamard bound of the transpose, which - happens to be sharp. :: - - sage: a = matrix(ZZ, 2, [2^10000,3^10000,2^50,3^19292]) - sage: import sage.matrix.misc - sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.change_ring(RR)) - 13976 - sage: len(str(a.det())) - 12215 - sage: sage.matrix.misc.hadamard_row_bound_mpfr(a.transpose().change_ring(RR)) - 12215 - - Note that in the above example using RDF would overflow:: - - sage: b = a.change_ring(RDF) - sage: b._hadamard_row_bound() - Traceback (most recent call last): - ... - OverflowError: cannot convert float infinity to integer - """ - if not isinstance(A.base_ring(), sage.rings.abc.RealField): - raise TypeError("A must have base field an mpfr real field.") - - cdef RealNumber a, b - cdef mpfr_t s, d, pr - cdef Py_ssize_t i, j - - mpfr_init(s) - mpfr_init(d) - mpfr_init(pr) - mpfr_set_si(d, 0, MPFR_RNDU) - - for i in range(A._nrows): - mpfr_set_si(s, 0, MPFR_RNDU) - for j in range(A._ncols): - sig_check() - a = A.get_unsafe(i, j) - mpfr_mul(pr, a.value, a.value, MPFR_RNDU) - mpfr_add(s, s, pr, MPFR_RNDU) - mpfr_log10(s, s, MPFR_RNDU) - mpfr_add(d, d, s, MPFR_RNDU) - b = a._new() - mpfr_set(b.value, d, MPFR_RNDU) - b /= 2 - mpfr_clear(s) - mpfr_clear(d) - mpfr_clear(pr) - return b.ceil() diff --git a/src/sage/matrix/misc_flint.pyx b/src/sage/matrix/misc_flint.pyx new file mode 100644 index 00000000000..37d326c9a74 --- /dev/null +++ b/src/sage/matrix/misc_flint.pyx @@ -0,0 +1,107 @@ +r""" +Misc matrix algorithms using FLINT +""" + +from cysignals.signals cimport sig_check + +from sage.arith.rational_reconstruction cimport mpq_rational_reconstruction +from sage.libs.gmp.mpq cimport * +from sage.libs.gmp.mpz cimport * +from sage.libs.flint.fmpq cimport fmpq_set_mpq, fmpq_canonicalise +from sage.libs.flint.fmpq_mat cimport fmpq_mat_entry_num, fmpq_mat_entry_den, fmpq_mat_entry +from sage.libs.flint.fmpz cimport fmpz_set_mpz, fmpz_one +from sage.rings.integer cimport Integer +from sage.rings.rational_field import QQ + +from .matrix_integer_dense cimport Matrix_integer_dense +from .matrix_rational_dense cimport Matrix_rational_dense + + +def matrix_integer_dense_rational_reconstruction(Matrix_integer_dense A, Integer N): + r""" + Given a matrix over the integers and an integer modulus, do + rational reconstruction on all entries of the matrix, viewed as + numbers mod `N`. This is done efficiently by assuming there is a + large common factor dividing the denominators. + + INPUT: + + - ``A`` -- matrix + - ``N`` -- an integer + + EXAMPLES:: + + sage: B = ((matrix(ZZ, 3,4, [1,2,3,-4,7,2,18,3,4,3,4,5])/3)%500).change_ring(ZZ) + sage: from sage.matrix.misc_flint import matrix_integer_dense_rational_reconstruction + sage: matrix_integer_dense_rational_reconstruction(B, 500) + [ 1/3 2/3 1 -4/3] + [ 7/3 2/3 6 1] + [ 4/3 1 4/3 5/3] + + TESTS: + + Check that :trac:`9345` is fixed:: + + sage: A = random_matrix(ZZ, 3) + sage: matrix_integer_dense_rational_reconstruction(A, 0) + Traceback (most recent call last): + ... + ZeroDivisionError: The modulus cannot be zero + """ + if not N: + raise ZeroDivisionError("The modulus cannot be zero") + cdef Matrix_rational_dense R + R = Matrix_rational_dense.__new__(Matrix_rational_dense, + A.parent().change_ring(QQ), 0,0,0) + + cdef mpz_t a, bnd, other_bnd, denom, tmp + cdef mpq_t qtmp + cdef Integer _bnd + cdef Py_ssize_t i, j + cdef int do_it + + mpz_init_set_si(denom, 1) + mpz_init(a) + mpz_init(tmp) + mpz_init(other_bnd) + mpq_init(qtmp) + + _bnd = (N//2).isqrt() + mpz_init_set(bnd, _bnd.value) + mpz_sub(other_bnd, N.value, bnd) + + for i in range(A._nrows): + for j in range(A._ncols): + sig_check() + A.get_unsafe_mpz(i, j, a) + if mpz_cmp_ui(denom, 1) != 0: + mpz_mul(a, a, denom) + mpz_fdiv_r(a, a, N.value) + do_it = 0 + if mpz_cmp(a, bnd) <= 0: + do_it = 1 + elif mpz_cmp(a, other_bnd) >= 0: + mpz_sub(a, a, N.value) + do_it = 1 + if do_it: + fmpz_set_mpz(fmpq_mat_entry_num(R._matrix, i, j), a) + if mpz_cmp_ui(denom, 1) != 0: + fmpz_set_mpz(fmpq_mat_entry_den(R._matrix, i, j), denom) + fmpq_canonicalise(fmpq_mat_entry(R._matrix, i, j)) + else: + fmpz_one(fmpq_mat_entry_den(R._matrix, i, j)) + else: + # Otherwise have to do it the hard way + A.get_unsafe_mpz(i, j, tmp) + mpq_rational_reconstruction(qtmp, tmp, N.value) + mpz_lcm(denom, denom, mpq_denref(qtmp)) + fmpq_set_mpq(fmpq_mat_entry(R._matrix, i, j), qtmp) + + mpz_clear(denom) + mpz_clear(a) + mpz_clear(tmp) + mpz_clear(other_bnd) + mpz_clear(bnd) + mpq_clear(qtmp) + + return R diff --git a/src/sage/matrix/misc_mpfr.pyx b/src/sage/matrix/misc_mpfr.pyx new file mode 100644 index 00000000000..48b6ade6379 --- /dev/null +++ b/src/sage/matrix/misc_mpfr.pyx @@ -0,0 +1,79 @@ +r""" +Misc matrix algorithms using MPFR +""" + +from cysignals.signals cimport sig_check + +cimport sage.rings.abc + +from sage.libs.mpfr cimport * +from sage.rings.real_mpfr cimport RealNumber + +from .matrix0 cimport Matrix + + +def hadamard_row_bound_mpfr(Matrix A): + r""" + Given a matrix `A` with entries that coerce to ``RR``, compute the row + Hadamard bound on the determinant. + + INPUT: + + - ``A`` -- a matrix over ``RR`` + + OUTPUT: + + integer -- an integer n such that the absolute value of the + determinant of this matrix is at most `10^n`. + + EXAMPLES: + + We create a very large matrix, compute the row Hadamard bound, + and also compute the row Hadamard bound of the transpose, which + happens to be sharp. :: + + sage: a = matrix(ZZ, 2, [2^10000, 3^10000, 2^50, 3^19292]) + sage: from sage.matrix.misc_mpfr import hadamard_row_bound_mpfr + sage: hadamard_row_bound_mpfr(a.change_ring(RR)) + 13976 + sage: len(str(a.det())) + 12215 + sage: hadamard_row_bound_mpfr(a.transpose().change_ring(RR)) + 12215 + + Note that in the above example using RDF would overflow:: + + sage: b = a.change_ring(RDF) + sage: b._hadamard_row_bound() + Traceback (most recent call last): + ... + OverflowError: cannot convert float infinity to integer + """ + if not isinstance(A.base_ring(), sage.rings.abc.RealField): + raise TypeError("A must have base field an mpfr real field.") + + cdef RealNumber a, b + cdef mpfr_t s, d, pr + cdef Py_ssize_t i, j + + mpfr_init(s) + mpfr_init(d) + mpfr_init(pr) + mpfr_set_si(d, 0, MPFR_RNDU) + + for i in range(A._nrows): + mpfr_set_si(s, 0, MPFR_RNDU) + for j in range(A._ncols): + sig_check() + a = A.get_unsafe(i, j) + mpfr_mul(pr, a.value, a.value, MPFR_RNDU) + mpfr_add(s, s, pr, MPFR_RNDU) + mpfr_log10(s, s, MPFR_RNDU) + mpfr_add(d, d, s, MPFR_RNDU) + b = a._new() + mpfr_set(b.value, d, MPFR_RNDU) + b /= 2 + mpfr_clear(s) + mpfr_clear(d) + mpfr_clear(pr) + return b.ceil()