diff --git a/src/sage/matrix/matrix1.pxd b/src/sage/matrix/matrix1.pxd index a74cae197ef..8d0545f47f7 100644 --- a/src/sage/matrix/matrix1.pxd +++ b/src/sage/matrix/matrix1.pxd @@ -1,4 +1,4 @@ cimport matrix0 cdef class Matrix(matrix0.Matrix): - pass + cdef _stack_impl(self, bottom) diff --git a/src/sage/matrix/matrix1.pyx b/src/sage/matrix/matrix1.pyx index eb42a18cc9e..1ae2d24c103 100644 --- a/src/sage/matrix/matrix1.pyx +++ b/src/sage/matrix/matrix1.pyx @@ -23,6 +23,7 @@ include "sage/ext/python.pxi" import sage.modules.free_module + cdef class Matrix(matrix0.Matrix): ################################################### # Coercion to Various Systems @@ -1116,8 +1117,10 @@ cdef class Matrix(matrix0.Matrix): ########################################################################### def stack(self, bottom, subdivide=False): r""" - Returns a new matrix formed by appending the matrix - (or vector) ``bottom`` beneath ``self``. + Return the matrix ``self`` on top of ``bottom``:: + + [ self ] + [ bottom ] INPUT: @@ -1142,6 +1145,7 @@ cdef class Matrix(matrix0.Matrix): in the result. .. warning:: + If ``subdivide`` is ``True`` then unequal column subdivisions will be discarded, since it would be ambiguous how to interpret them. If the subdivision behavior is not what you need, @@ -1187,13 +1191,13 @@ cdef class Matrix(matrix0.Matrix): sage: A.stack(B) Traceback (most recent call last): ... - TypeError: number of columns must be the same, 2 != 3 + TypeError: number of columns must be the same, not 2 and 3 sage: v = vector(RR, [100, 200, 300]) sage: A.stack(v) Traceback (most recent call last): ... - TypeError: number of columns must be the same, 2 != 3 + TypeError: number of columns must be the same, not 2 and 3 Setting ``subdivide`` to ``True`` will, in its simplest form, add a subdivision between ``self`` and ``bottom``. :: @@ -1251,16 +1255,17 @@ cdef class Matrix(matrix0.Matrix): [ 5 6 7 8 9] [10 11 12 13 14] - The result retains the base ring of ``self`` by coercing - the elements of ``bottom`` into the base ring of ``self``. :: + The base ring of the result is the common parent for the base + rings of ``self`` and ``bottom``. In particular, the parent for + ``A.stack(B)`` and ``B.stack(A)`` should be equal:: sage: A = matrix(QQ, 1, 2, [1,2]) sage: B = matrix(RR, 1, 2, [sin(1.1), sin(2.2)]) sage: C = A.stack(B); C - [ 1 2] - [183017397/205358938 106580492/131825561] + [ 1.00000000000000 2.00000000000000] + [0.891207360061435 0.808496403819590] sage: C.parent() - Full MatrixSpace of 2 by 2 dense matrices over Rational Field + Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision sage: D = B.stack(A); D [0.891207360061435 0.808496403819590] @@ -1268,32 +1273,52 @@ cdef class Matrix(matrix0.Matrix): sage: D.parent() Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision - Sometimes it is not possible to coerce into the base ring - of ``self``. A solution is to change the base ring of ``self`` - to a more expansive ring. Here we mix the rationals with - a ring of polynomials with rational coefficients. :: + :: - sage: R = PolynomialRing(QQ, 'y') - sage: A = matrix(QQ, 1, 2, [1,2]) - sage: B = matrix(R, 1, 2, ['y', 'y^2']) + sage: R. = PolynomialRing(ZZ) + sage: A = matrix(QQ, 1, 2, [1, 2/3]) + sage: B = matrix(R, 1, 2, [y, y^2]) - sage: C = B.stack(A); C + sage: C = A.stack(B); C + [ 1 2/3] [ y y^2] - [ 1 2] sage: C.parent() Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field - sage: D = A.stack(B) - Traceback (most recent call last): - ... - TypeError: not a constant polynomial + Stacking a dense matrix atop a sparse one returns a sparse + matrix:: - sage: E = A.change_ring(R) - sage: F = E.stack(B); F - [ 1 2] - [ y y^2] - sage: F.parent() - Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in y over Rational Field + sage: M = Matrix(ZZ, 2, 3, range(6), sparse=False) + sage: N = diagonal_matrix([10,11,12], sparse=True) + sage: P = M.stack(N); P + [ 0 1 2] + [ 3 4 5] + [10 0 0] + [ 0 11 0] + [ 0 0 12] + sage: P.is_sparse() + True + sage: P = N.stack(M); P + [10 0 0] + [ 0 11 0] + [ 0 0 12] + [ 0 1 2] + [ 3 4 5] + sage: P.is_sparse() + True + + One can stack matrices over different rings (:trac:`16399`). :: + + sage: M = Matrix(ZZ, 2, 3, range(6)) + sage: N = Matrix(QQ, 1, 3, [10,11,12]) + sage: M.stack(N) + [ 0 1 2] + [ 3 4 5] + [10 11 12] + sage: N.stack(M) + [10 11 12] + [ 0 1 2] + [ 3 4 5] TESTS: @@ -1306,43 +1331,80 @@ cdef class Matrix(matrix0.Matrix): [ 3 4 5] [10 11 12] - AUTHOR: + Non-matrices fail gracefully:: - - Rob Beezer (2011-03-19) - rewritten to mirror code for :meth:`augment` - """ - from sage.matrix.constructor import matrix + sage: M.stack(polygen(QQ)) + Traceback (most recent call last): + ... + TypeError: a matrix must be stacked with another matrix or a vector + + AUTHORS: - if hasattr(bottom, '_vector_'): - bottom = bottom.row() - if not isinstance(bottom, sage.matrix.matrix1.Matrix): - raise TypeError('a matrix must be stacked with another matrix, or ' - 'a vector') + - Rob Beezer (2011-03-19): rewritten to mirror code for :meth:`augment` + - Jeroen Demeyer (2015-01-06): refactor, see :trac:`16399`. + Put all boilerplate in one place (here) and put the actual + type-dependent implementation in ``_stack_impl``. + """ cdef Matrix other - other = bottom + if isinstance(bottom, Matrix): + other = bottom + else: + if hasattr(bottom, '_vector_'): + bottom = bottom.row() + else: + raise TypeError('a matrix must be stacked with ' + 'another matrix or a vector') + other = bottom if self._ncols != other._ncols: - raise TypeError('number of columns must be the same, ' - '{0} != {1}'.format(self._ncols, other._ncols)) - if not (self._base_ring is other.base_ring()): - other = other.change_ring(self._base_ring) + raise TypeError("number of columns must be the same, not %s and %s" % + (self.ncols(), bottom.ncols()) ) + + top_ring = self._base_ring + bottom_ring = other._base_ring + if top_ring is not bottom_ring: + from sage.structure.element import get_coercion_model + coercion_model = get_coercion_model() + R = coercion_model.common_parent(top_ring, bottom_ring) + if top_ring is not R: + self = self.change_ring(R) + if bottom_ring is not R: + other = other.change_ring(R) + + if type(self) is not type(other): + # If one of the matrices is sparse, return a sparse matrix + if self.is_sparse_c() and not other.is_sparse_c(): + other = other.sparse_matrix() + elif other.is_sparse_c() and not self.is_sparse_c(): + self = self.sparse_matrix() + + Z = self._stack_impl(other) + if subdivide: + Z._subdivide_on_stack(self, other) + return Z + cdef _stack_impl(self, bottom): + """ + Implementation of :meth:`stack`. + + Assume that ``self`` and ``other`` are compatible in the sense + that they have the same base ring and that both are either + dense or sparse. + """ + cdef Matrix other = bottom cdef Matrix Z - Z = self.new_matrix(nrows = self._nrows + other._nrows) + Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols) cdef Py_ssize_t r, c - for r from 0 <= r < self._nrows: - for c from 0 <= c < self._ncols: - Z.set_unsafe(r,c, self.get_unsafe(r,c)) - nr = self.nrows() - - for r from 0 <= r < other._nrows: - for c from 0 <= c < other._ncols: + cdef Py_ssize_t nr = self._nrows + for r in range(self._nrows): + for c in range(self._ncols): + Z.set_unsafe(r, c, self.get_unsafe(r,c)) + for r in range(other._nrows): + for c in range(other._ncols): Z.set_unsafe(r+nr, c, other.get_unsafe(r,c)) - if subdivide: - Z._subdivide_on_stack(self, other) - return Z def augment(self, right, subdivide=False): diff --git a/src/sage/matrix/matrix_integer_dense.pyx b/src/sage/matrix/matrix_integer_dense.pyx index a352e96638d..c4d6104f2d0 100644 --- a/src/sage/matrix/matrix_integer_dense.pyx +++ b/src/sage/matrix/matrix_integer_dense.pyx @@ -95,7 +95,7 @@ from sage.rings.integer_ring cimport IntegerRing_class from sage.rings.finite_rings.integer_mod_ring import IntegerModRing from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.polynomial.polynomial_integer_dense_flint cimport Polynomial_integer_dense_flint -from sage.structure.element cimport ModuleElement, RingElement, Element, Vector, CoercionModel +from sage.structure.element cimport ModuleElement, RingElement, Element, Vector from sage.structure.element import is_Vector from sage.structure.sequence import Sequence @@ -4679,9 +4679,9 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse ################################################################# # operations with matrices ################################################################# - def stack(self, bottom, subdivide=False): + cdef _stack_impl(self, bottom): r""" - Return the matrix ``self`` on top of ``bottom``: + Return the matrix ``self`` on top of ``bottom``:: [ self ] [ bottom ] @@ -4718,73 +4718,21 @@ cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse [-----------] [ 0 1 2 3] [ 4 5 6 7] + """ + cdef Matrix_integer_dense other = bottom + cdef Matrix_integer_dense Z + Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols) - TESTS: - - Stacking a dense matrix atop a sparse one should work:: - - sage: M = Matrix(ZZ, 2, 3, range(6)) - sage: M.is_sparse() - False - sage: N = diagonal_matrix([10,11,12], sparse=True) - sage: N.is_sparse() - True - sage: P = M.stack(N); P - [ 0 1 2] - [ 3 4 5] - [10 0 0] - [ 0 11 0] - [ 0 0 12] - sage: P.is_sparse() - False - - One can stack matrices over different rings (:trac:`16399`). :: - - sage: M = Matrix(ZZ, 2, 3, range(6)) - sage: N = Matrix(QQ, 1, 3, [10,11,12]) - sage: M.stack(N) - [0 1 2 3 4] - [5 6 7 8 9] - [0 1 2 3 4] - sage: N.stack(M) - ? - sage: M2 = Matrix(ZZ['x'], 2, 3, range(6)) - sage: N.stack(M2) - ? - """ - if hasattr(bottom, '_vector_'): - bottom = bottom.row() - if self._ncols != bottom.ncols(): - raise TypeError("number of columns must be the same") - cdef CoercionModel coercion_model - - top_ring = self._base_ring - bottom_ring = bottom.base_ring() - if not (top_ring is bottom_ring): - if top_ring.has_coerce_map_from(bottom_ring): - bottom = bottom.change_ring(top_ring) - elif bottom_ring.has_coerce_map_from(top_ring): - new_top = self.change_ring(bottom_ring) - return new_top.stack(bottom, subdivide=subdivide): - else: - from sage.structure.element import get_coercion_model - coercion_model = get_coercion_model() - com_ring = coercion_model.common_parent(top_ring, bottom_ring) - self = self.change_ring(com_ring) # allowed ? - bottom = other.change_ring(com_ring) + cdef Py_ssize_t r, c + cdef Py_ssize_t nr = self._nrows + for r in range(self._nrows): + for c in range(self._ncols): + fmpz_set(fmpz_mat_entry(Z._matrix, r, c),fmpz_mat_entry(self._matrix, r, c)) + for r in range(other._nrows): + for c in range(other._ncols): + fmpz_set(fmpz_mat_entry(Z._matrix, r+nr, c),fmpz_mat_entry(other._matrix, r, c)) - cdef Matrix_integer_dense other = bottom.dense_matrix() - cdef Matrix_integer_dense M - M = self.new_matrix(nrows = self._nrows + other._nrows, ncols = self.ncols()) - cdef Py_ssize_t i, j, k - for j from 0 <= j < self._ncols: - for i from 0 <= i < self._nrows: - fmpz_set(fmpz_mat_entry(M._matrix,i,j),fmpz_mat_entry(self._matrix,i,j)) - for i from 0 <= i < other._nrows: - fmpz_set(fmpz_mat_entry(M._matrix,self._nrows + i,j),fmpz_mat_entry(other._matrix,i,j)) - if subdivide: - M._subdivide_on_stack(self, other) - return M + return Z def augment(self, right, subdivide=False): r""" diff --git a/src/sage/matrix/matrix_mod2_dense.pyx b/src/sage/matrix/matrix_mod2_dense.pyx index 445427394ce..b7934b51f9b 100644 --- a/src/sage/matrix/matrix_mod2_dense.pyx +++ b/src/sage/matrix/matrix_mod2_dense.pyx @@ -1575,7 +1575,7 @@ cdef class Matrix_mod2_dense(matrix_dense.Matrix_dense): # dense or sparse Z._subdivide_on_augment(self, other) return Z - def stack(self, bottom, subdivide=False): + cdef _stack_impl(self, bottom): r""" Stack ``self`` on top of ``bottom``. @@ -1641,25 +1641,11 @@ cdef class Matrix_mod2_dense(matrix_dense.Matrix_dense): # dense or sparse sage: M.stack(N) [] """ - if hasattr(bottom, '_vector_'): - bottom = bottom.row() - cdef Matrix_mod2_dense other = bottom - - if self._ncols != other._ncols: - raise TypeError("Both numbers of columns must match.") - - if self._nrows == 0: - return other.__copy__() - if other._nrows == 0: - return self.__copy__() - + cdef Matrix_mod2_dense other = bottom cdef Matrix_mod2_dense Z - Z = self.new_matrix(nrows = self._nrows + other._nrows) - if self._ncols == 0: - return Z - Z._entries = mzd_stack(Z._entries, self._entries, other._entries) - if subdivide: - Z._subdivide_on_stack(self, other) + Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols) + if self._ncols > 0: + Z._entries = mzd_stack(Z._entries, self._entries, other._entries) return Z def submatrix(self, lowr, lowc, nrows , ncols): @@ -1763,7 +1749,8 @@ cdef class Matrix_mod2_dense(matrix_dense.Matrix_dense): # dense or sparse '1 0 1 1 1 0' """ cdef Py_ssize_t i, j, k, n - cdef char *s, *t + cdef char *s + cdef char *t if self._nrows == 0 or self._ncols == 0: data = '' @@ -1845,7 +1832,8 @@ cdef class Matrix_mod2_dense(matrix_dense.Matrix_dense): # dense or sparse if self._nrows == 0 or self._ncols == 0: return 0 cdef mzd_t *A = mzd_copy(NULL, self._entries) - cdef mzp_t *P, *Q + cdef mzp_t *P + cdef mzp_t *Q if algorithm == 'ple': P = mzp_init(self._entries.nrows) diff --git a/src/sage/matrix/matrix_sparse.pyx b/src/sage/matrix/matrix_sparse.pyx index ed5ee0747a6..8356904c614 100644 --- a/src/sage/matrix/matrix_sparse.pyx +++ b/src/sage/matrix/matrix_sparse.pyx @@ -4,11 +4,10 @@ Base class for sparse matrices cimport matrix cimport matrix0 -from sage.structure.element cimport (Element, RingElement, ModuleElement, Vector, CoercionModel) +from sage.structure.element cimport Element, RingElement, ModuleElement, Vector from sage.rings.ring import is_Ring from sage.misc.misc import verbose -include 'sage/ext/cdefs.pxi' include 'sage/ext/stdsage.pxi' include 'sage/ext/python.pxi' include 'sage/ext/interrupt.pxi' @@ -18,9 +17,9 @@ cdef extern from "Python.h": PyObject* PyList_GET_ITEM0 "PyList_GET_ITEM" (PyObject* list, Py_ssize_t i) Py_ssize_t PyNumber_AsSsize_t(PyObject* o, PyObject* exc) - import sage.matrix.matrix_space + cdef class Matrix_sparse(matrix.Matrix): cdef bint is_sparse_c(self): @@ -855,9 +854,9 @@ cdef class Matrix_sparse(matrix.Matrix): A.set_unsafe(new_row, new_col, entry) return A - def stack(self, bottom, subdivide=False): + cdef _stack_impl(self, bottom): r""" - Stack ``self`` on top of ``bottom``. + Stack ``self`` on top of ``bottom``:: [ self ] [ bottom ] @@ -902,50 +901,28 @@ cdef class Matrix_sparse(matrix.Matrix): sage: M = Matrix(ZZ, 2, 3, range(6), sparse=True) sage: N = Matrix(QQ, 1, 3, [10,11,12], sparse=True) sage: M.stack(N) - [0 1 2 3 4] - [5 6 7 8 9] - [0 1 2 3 4] + [ 0 1 2] + [ 3 4 5] + [10 11 12] sage: N.stack(M) - ? + [10 11 12] + [ 0 1 2] + [ 3 4 5] sage: M2 = Matrix(ZZ['x'], 2, 3, range(6), sparse=True) sage: N.stack(M2) - ? + [10 11 12] + [ 0 1 2] + [ 3 4 5] """ - if hasattr(bottom, '_vector_'): - bottom = bottom.row() - if not isinstance(bottom, matrix.Matrix): - raise TypeError("other must be a matrix") - - cdef Matrix_sparse other = bottom.sparse_matrix() - cdef CoercionModel coercion_model - - if self._ncols != other._ncols: - raise TypeError("number of columns must be the same") - - top_ring = self._base_ring - bottom_ring = other.base_ring() - if not (top_ring is bottom_ring): - if top_ring.has_coerce_map_from(bottom_ring): - other = other.change_ring(top_ring) - elif bottom_ring.has_coerce_map_from(top_ring): - self = self.change_ring(bottom_ring) # allowed ? - else: - from sage.structure.element import get_coercion_model - coercion_model = get_coercion_model() - com_ring = coercion_model.common_parent(top_ring, bottom_ring) - self = self.change_ring(com_ring) # allowed ? - other = other.change_ring(com_ring) - - + cdef Matrix_sparse other = bottom cdef Matrix_sparse Z - Z = self.new_matrix(nrows = self._nrows + other.nrows()) + Z = self.new_matrix(nrows=self._nrows + other._nrows, ncols=self._ncols) for i, j in self.nonzero_positions(copy=False): Z.set_unsafe(i, j, self.get_unsafe(i,j)) for i, j in other.nonzero_positions(copy=False): Z.set_unsafe(i + self._nrows, j, other.get_unsafe(i,j)) - if subdivide: - Z._subdivide_on_stack(self, other) + return Z def augment(self, right, subdivide=False): diff --git a/src/sage/modules/diamond_cutting.py b/src/sage/modules/diamond_cutting.py index c6098fa5838..ea98105c2c6 100644 --- a/src/sage/modules/diamond_cutting.py +++ b/src/sage/modules/diamond_cutting.py @@ -259,7 +259,7 @@ def calculate_voronoi_cell(basis, radius=None, verbose=False): artificial_length = None if dim[0] < dim[1]: # introduce "artificial" basis points (representing infinity) - artificial_length = ceil(max(abs(v) for v in basis)) * 2 + artificial_length = max(abs(v) for v in basis).ceil() * 2 additional_vectors = identity_matrix(dim[1]) * artificial_length basis = basis.stack(additional_vectors) # LLL-reduce to get quadratic matrix diff --git a/src/sage/modules/fg_pid/fgp_morphism.py b/src/sage/modules/fg_pid/fgp_morphism.py index 80fe02ba823..86bbe370955 100644 --- a/src/sage/modules/fg_pid/fgp_morphism.py +++ b/src/sage/modules/fg_pid/fgp_morphism.py @@ -423,6 +423,7 @@ def lift(self, x): # Stack it on top of the basis for W'. Wp = CD.V().coordinate_module(CD.W()).basis_matrix() + Wp = Wp.change_ring(A.base_ring()) B = A.stack(Wp) # Compute Hermite form of C with transformation