Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/doc/en/reference/matrices/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
8 changes: 4 additions & 4 deletions src/sage/matrix/matrix2.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -14883,21 +14883,21 @@ 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()
A = A.transpose()
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):
Expand Down
2 changes: 1 addition & 1 deletion src/sage/matrix/matrix_cyclo_dense.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/sage/matrix/matrix_integer_dense.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
12 changes: 7 additions & 5 deletions src/sage/matrix/matrix_integer_sparse.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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):
...
Expand Down
209 changes: 23 additions & 186 deletions src/sage/matrix/misc.pyx
Original file line number Diff line number Diff line change
@@ -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]
Expand Down Expand Up @@ -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)
Expand All @@ -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()
Loading