Coverage for local/lib/python2.7/site-packages/sage/matrix/matrix_rational_dense.pyx : 78%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" Dense matrices over the rational field
EXAMPLES:
We create a 3x3 matrix with rational entries and do some operations with it.
::
sage: a = matrix(QQ, 3,3, [1,2/3, -4/5, 1,1,1, 8,2, -3/19]); a [ 1 2/3 -4/5] [ 1 1 1] [ 8 2 -3/19] sage: a.det() 2303/285 sage: a.charpoly() x^3 - 35/19*x^2 + 1259/285*x - 2303/285 sage: b = a^(-1); b [ -615/2303 -426/2303 418/2303] [ 2325/2303 1779/2303 -513/2303] [-1710/2303 950/2303 95/2303] sage: b.det() 285/2303 sage: a == b False sage: a < b False sage: b < a True sage: a > b True sage: a*b [1 0 0] [0 1 0] [0 0 1]
TESTS::
sage: a = matrix(QQ, 2, range(4), sparse=False) sage: TestSuite(a).run()
Test hashing::
sage: m = matrix(QQ, 2, [1/2, -1, 2, 3]) sage: hash(m) Traceback (most recent call last): ... TypeError: mutable matrices are unhashable sage: m.set_immutable() sage: hash(m) 2212268000387745777 # 64-bit 1997752305 # 32-bit """
#***************************************************************************** # Copyright (C) 2004,2005,2006 William Stein <wstein@gmail.com> # 2017 Vincent Delecroix <20100.delecroix@gmail.com> # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 2 of the License, or # (at your option) any later version. # http://www.gnu.org/licenses/ #*****************************************************************************
from __future__ import absolute_import, print_function
from libc.string cimport strcpy, strlen
from sage.cpython.string cimport char_to_str, str_to_bytes
from sage.modules.vector_rational_dense cimport Vector_rational_dense from sage.ext.stdsage cimport PY_NEW from sage.misc.randstate cimport randstate, current_randstate
from sage.modules.vector_rational_dense cimport Vector_rational_dense
from cysignals.signals cimport sig_on, sig_off, sig_check from cysignals.memory cimport sig_malloc, sig_free
from sage.arith.rational_reconstruction cimport mpq_rational_reconstruction
from sage.libs.gmp.types cimport mpz_t, mpq_t from sage.libs.gmp.mpz cimport mpz_init, mpz_clear, mpz_cmp_si from sage.libs.gmp.mpq cimport mpq_init, mpq_clear, mpq_set_si, mpq_mul, mpq_add, mpq_set from sage.libs.gmp.randomize cimport (mpq_randomize_entry, mpq_randomize_entry_as_int, mpq_randomize_entry_recip_uniform, mpq_randomize_entry_nonzero, mpq_randomize_entry_as_int_nonzero, mpq_randomize_entry_recip_uniform_nonzero)
from sage.libs.flint.fmpz cimport * from sage.libs.flint.fmpq cimport * from sage.libs.flint.fmpz_mat cimport * from sage.libs.flint.fmpq_mat cimport *
cimport sage.structure.element
from sage.structure.sequence import Sequence from sage.rings.rational cimport Rational from .matrix cimport Matrix from .matrix_integer_dense cimport Matrix_integer_dense, _lift_crt from sage.structure.element cimport ModuleElement, RingElement, Element, Vector from sage.rings.integer cimport Integer from sage.rings.ring import is_Ring from sage.rings.integer_ring import ZZ, is_IntegerRing from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF from sage.rings.finite_rings.integer_mod_ring import is_IntegerModRing from sage.rings.rational_field import QQ from sage.arith.all import gcd
from .matrix2 import decomp_seq from .matrix0 import Matrix as Matrix_base
from sage.misc.all import verbose, get_verbose, prod
######################################################### # PARI C library from cypari2.gen cimport Gen from sage.libs.pari.all import PariError from sage.libs.pari.convert_gmp cimport INTFRAC_to_mpq from sage.libs.pari.convert_flint cimport rational_matrix, _new_GEN_from_fmpq_mat_t from cypari2.stack cimport clear_stack from cypari2.paridecl cimport * #########################################################
cdef class Matrix_rational_dense(Matrix_dense): def __cinit__(self, parent, entries, copy, coerce): """ Create and allocate memory for the matrix.
EXAMPLES::
sage: from sage.matrix.matrix_rational_dense import Matrix_rational_dense sage: a = Matrix_rational_dense.__new__(Matrix_rational_dense, Mat(ZZ,3), 0,0,0) sage: type(a) <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'>
.. warning::
This is for internal use only, or if you really know what you're doing. """
cdef inline Matrix_rational_dense _new_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols): else:
def __dealloc__(self):
def __init__(self, parent, entries=None, coerce=True, copy=True): r""" TESTS::
sage: matrix(QQ, 2, 2, 1/4) [1/4 0] [ 0 1/4] sage: matrix(QQ, 3, 1, [1/2, -3/4, 0]) [ 1/2] [-3/4] [ 0] """ cdef Py_ssize_t i, j, k cdef Rational z
entries = list(entries) raise TypeError("entries has the wrong length")
# TODO: Should use an unsafe un-bounds-checked array access here. else: # TODO: Should use an unsafe un-bounds-checked array access here.
else: # is it a scalar? # Try to coerce entries to a scalar (an integer)
def matrix_from_columns(self, columns): """ Return the matrix constructed from self using columns with indices in the columns list.
EXAMPLES::
sage: A = matrix(QQ, 3, range(9)) sage: A [0 1 2] [3 4 5] [6 7 8] sage: A.matrix_from_columns([2,1]) [2 1] [5 4] [8 7] sage: A.matrix_from_columns((2,1,0,2)) [2 1 0 2] [5 4 3 5] [8 7 6 8] """ cdef Matrix_rational_dense A cdef Py_ssize_t ncols, k, r, col
raise IndexError("column out of range")
def add_to_entry(self, Py_ssize_t i, Py_ssize_t j, elt): r""" Add ``elt`` to the entry at position ``(i,j)``
EXAMPLES::
sage: m = matrix(QQ, 2, 2) sage: m.add_to_entry(0, 0, -1/3) sage: m [-1/3 0] [ 0 0] """ elt = Rational(elt) i += self._nrows raise IndexError("row index out of range") j += self._ncols raise IndexError("column index out of range") cdef fmpq_t tmp fmpq_mat_entry(self._matrix, i, j), tmp)
cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): cdef Rational x
cdef _add_ui_unsafe_assuming_int(self, Py_ssize_t i, Py_ssize_t j, unsigned long int n): # doesn't check immutability # doesn't do bounds checks. # assumes that self[i,j] is an integer.
cdef _sub_ui_unsafe_assuming_int(self, Py_ssize_t i, Py_ssize_t j, unsigned long int n): # doesn't check immutability # doesn't do bounds checks. # assumes that self[i,j] is an integer.
def _pickle(self):
def _unpickle(self, data, int version): else: raise RuntimeError("unknown matrix version (=%s)"%version)
cdef _pickle_version0(self):
cpdef _export_as_string(self, int base=10): """ Return space separated string of the entries in this matrix, in the given base. This is optimized for speed.
INPUT:
- ``base`` - an optional integer (default is ``10``)
EXAMPLES::
sage: m = matrix(QQ,2,3,[1,2/3,-3/4,1,-2/3,-45/17]) sage: m._export_as_string(10) '1 2/3 -3/4 1 -2/3 -45/17' sage: m._export_as_string(16) '1 2/3 -3/4 1 -2/3 -2d/11' """ cdef Py_ssize_t i, j, len_so_far, m, n cdef char *a cdef char *s cdef char *t cdef char *tmp
else:
m = fmpz_sizeinbase (fmpq_mat_entry_num(self._matrix, i, j), base) + \ # copy to new string with double the size
cdef _unpickle_version0(self, data): r""" TESTS::
sage: a = random_matrix(QQ, 4, 3, num_bound=2**500, den_bound=2**500) sage: loads(dumps(a)) == a # indirect doctest True """ cdef Py_ssize_t i, j, k raise RuntimeError("invalid pickle data") raise RuntimeError("invalid pickle data") else: raise RuntimeError("invalid pickle data")
######################################################################## # LEVEL 2 functionality # x * cdef _add_ # x * cdef _mul_ # x * cdef _vector_times_matrix_ # x * cpdef _cmp_ # x * __neg__ # * __invert__ # x * __copy__ # x * _multiply_classical # * _list -- list of underlying elements (need not be a copy) # * _dict -- sparse dictionary of underlying elements (need not be a copy) ########################################################################
cpdef _lmul_(self, Element right): """ EXAMPLES::
sage: a = matrix(QQ, 2, range(6)) sage: (3/4) * a [ 0 3/4 3/2] [ 9/4 3 15/4] """ cdef Matrix_rational_dense M cdef fmpq_t x
cpdef _add_(self, right): """ Add two dense matrices over QQ.
EXAMPLES::
sage: a = MatrixSpace(QQ,3)(range(9)) sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)]) sage: a+b [ 1 3/2 7/3] [13/4 21/5 31/6] [43/7 57/8 73/9] sage: b.swap_rows(1,2) sage: #a+b """ cdef Matrix_rational_dense ans
cpdef _sub_(self, right): """ Subtract two dense matrices over QQ.
EXAMPLES::
sage: a = MatrixSpace(QQ,3)(range(9)) sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)]) sage: a-b [ -1 1/2 5/3] [11/4 19/5 29/6] [41/7 55/8 71/9] """ cdef Matrix_rational_dense ans
cpdef int _cmp_(self, right) except -2: r""" TESTS::
sage: M = MatrixSpace(QQ, 1) sage: M(1) < M(2) True sage: M(1/3) >= M(5/2) False sage: M(2) == M(2) True sage: M(3/4) != M(2) True
sage: matrix(QQ, 2, 3) == matrix(QQ, 2, 3) True sage: matrix(QQ, 2, 2) == matrix(QQ, 2, 3) False sage: matrix(QQ, 2, 2) == matrix(QQ, 3, 2) False sage: matrix(QQ, 2, 3) == matrix(QQ, 3, 2) False
sage: mats = [matrix(QQ, 2, 2, 1), matrix(QQ, 2, 2, -1), matrix(QQ, 2, 2, 0)] sage: mats.sort() sage: mats == [-1, 0, 1] True """ cdef Py_ssize_t i, j cdef int k for i in range(self._nrows): for j in range(self._ncols): k = fmpq_cmp(fmpq_mat_entry(self._matrix, i, j), fmpq_mat_entry((<Matrix_rational_dense> right)._matrix, i, j)) if k: return (k > 0) - (k < 0) return 0
cdef _vector_times_matrix_(self, Vector v): """ Returns the vector times matrix product.
INPUT:
- ``v`` - a free module element.
OUTPUT: The vector times matrix product v\*A.
EXAMPLES::
sage: B = matrix(QQ,2, [1,2,3,4]) sage: V = QQ^2 sage: w = V([-1,5/2]) sage: w * B (13/2, 8) """ cdef Vector_rational_dense w, ans cdef Py_ssize_t i, j cdef mpq_t x, y, z
def __neg__(self): """ Negate a matrix over QQ.
EXAMPLES::
sage: a = matrix(QQ, 3, [1/n for n in range(1,10)]) sage: -a [ -1 -1/2 -1/3] [-1/4 -1/5 -1/6] [-1/7 -1/8 -1/9] """ cdef Matrix_rational_dense ans
def __copy__(self): """ Copy a matrix over QQ.
TESTS::
sage: a = matrix(QQ, 3, [1/n for n in range(1,10)]) sage: b = a.__copy__() sage: a == b True sage: a is b False sage: b[0,0] = 5 sage: a == b False
sage: a.subdivide(2, 1) sage: b = a.__copy__() sage: b.subdivisions() ([2], [1]) sage: a.subdivide(2, 2) sage: b.subdivisions() ([2], [1]) """ cdef Matrix_rational_dense ans
######################################################################## # LEVEL 3 functionality (Optional) # x * cdef _sub_ # * __deepcopy__ # * __invert__ # * Matrix windows -- only if you need strassen for that base # * Other functions (list them here): # x * denom(self): # x * mpz_denom(self, mpz_t d): # x * _clear_denom(self): # o * echelon_modular(self, height_guess=None): ######################################################################## def __invert__(self): """ EXAMPLES::
sage: a = matrix(QQ,3,range(9)) sage: a.inverse() Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular sage: a = matrix(QQ, 2, [1, 5, 17, 3]) sage: a.inverse() [-3/82 5/82] [17/82 -1/82] sage: ~matrix(QQ, 2, 3) Traceback (most recent call last): ... ArithmeticError: self must be a square matrix """
def _invert_flint(self): r""" TESTS::
sage: matrix(QQ, 2, [1,2,3,4])._invert_flint() [ -2 1] [ 3/2 -1/2] sage: matrix(QQ, 1)._invert_flint() Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular """ cdef int ret cdef Matrix_rational_dense ans
def inverse(self, algorithm=None, check_invertible=True): """ Return the inverse of this matrix
INPUT:
- ``algorithm`` -- an optional specification of an algorithm. It can be one of
- ``None``: (default) uses flint
- ``'flint'``: uses flint library
- ``'pari'``: uses PARI library
- ``'iml'``: uses IML library
- ``check_invertible`` - only used when ``algorithm=iml``. Whether to check that matrix is invertible
EXAMPLES::
sage: a = matrix(QQ,3,[1,2,5,3,2,1,1,1,1,]) sage: a.inverse() [1/2 3/2 -4] [ -1 -2 7] [1/2 1/2 -2]
sage: a = matrix(QQ, 2, [1, 5, 17, 3]) sage: a.inverse(algorithm="flint") [-3/82 5/82] [17/82 -1/82] sage: a.inverse(algorithm="flint") * a [1 0] [0 1]
sage: a = matrix(QQ, 2, [-1, 5, 12, -3]) sage: a.inverse(algorithm="iml") [1/19 5/57] [4/19 1/57] sage: a.inverse(algorithm="iml") * a [1 0] [0 1]
sage: a = matrix(QQ, 4, primes_first_n(16)) sage: a.inverse(algorithm="pari") [ 3/11 -12/55 -1/5 2/11] [ -5/11 -2/55 3/10 -3/22] [ -13/22 307/440 -1/10 -9/88] [ 15/22 -37/88 0 7/88]
On singular matrices this method raises a ``ZeroDivisionError``::
sage: a = matrix(QQ, 2) sage: a.inverse(algorithm="flint") Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular sage: a.inverse(algorithm="iml") Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular sage: a.inverse(algorithm="pari") Traceback (most recent call last): ... ZeroDivisionError: input matrix must be nonsingular
TESTS::
sage: a = matrix(QQ, 2) sage: a.inverse(algorithm="IAmNotAnAlgorithm") Traceback (most recent call last): ... ValueError: unknown algorithm 'IAmNotAnAlgorithm'
sage: for _ in range(30): ....: dim = randint(1, 20) ....: a = random_matrix(QQ, dim, num_bound=10, den_bound=10) ....: while a.rank() != dim: a = random_matrix(QQ, dim) ....: inv_flint = a.inverse(algorithm='flint') ....: inv_pari = a.inverse(algorithm='pari') ....: inv_iml = a.inverse(algorithm='iml') ....: assert inv_flint == inv_pari == inv_iml """
else:
def determinant(self, algorithm=None, proof=None): """ Return the determinant of this matrix.
INPUT:
- ``algorithm`` -- an optional specification of an algorithm. It can be one of
- ``None``: (default) uses flint
- ``'flint'``: uses flint library
- ``'pari'``: uses PARI library
- ``'integer'``: removes denominator and call determinant on the corresponding integer matrix
- ``'generic'``: calls the generic Sage implementation
- ``proof`` - bool or None; if None use proof.linear_algebra(); only relevant for the padic algorithm.
.. NOTE::
It would be *VERY VERY* hard for det to fail even with proof=False.
EXAMPLES::
sage: m = matrix(QQ,3,[1,2/3,4/5, 2,2,2, 5,3,2/5]) sage: m.determinant() -34/15 sage: m.charpoly() x^3 - 17/5*x^2 - 122/15*x + 34/15
sage: m = matrix(QQ, 3, [(1/i)**j for i in range(2,5) for j in range(3)]) sage: m.determinant(algorithm="flint") -1/288
sage: m = matrix(QQ, 4, [(-1)**n/n for n in range(1,17)]) sage: m.determinant(algorithm="pari") 2/70945875
sage: m = matrix(QQ, 5, [1/(i+j+1) for i in range(5) for j in range(5)]) sage: m.determinant(algorithm="integer") 1/266716800000
On non-square matrices, the method raises a ``ValueError``::
sage: matrix(QQ, 2, 3).determinant(algorithm='flint') Traceback (most recent call last): ... ValueError: non square matrix sage: matrix(QQ, 2, 3).determinant(algorithm='pari') Traceback (most recent call last): ... ValueError: non square matrix sage: matrix(QQ, 2, 3).determinant(algorithm='integer') Traceback (most recent call last): ... ValueError: non square matrix sage: matrix(QQ, 2, 3).determinant(algorithm='generic') Traceback (most recent call last): ... ValueError: non square matrix
TESTS:
Check that the four algorithms agree::
sage: for _ in range(20): ....: dim = randint(0, 30) ....: m = random_matrix(QQ, dim, num_bound=10, den_bound=10) ....: det_flint = m.determinant("flint"); m._clear_cache() ....: det_pari = m.determinant("pari"); m._clear_cache() ....: det_int = m.determinant("integer"); m._clear_cache() ....: det_gen = m.determinant("generic") ....: assert det_flint == det_pari == det_int == det_gen """
else: raise ValueError("unknown algorithm '%s'"%algorithm)
def _det_flint(self): r""" Return the determinant (computed using flint)
EXAMPLES::
sage: matrix(QQ, 2, [1/3, 2/5, 3/4, 7/8])._det_flint() -1/120 sage: matrix(QQ, 0)._det_flint() 1 sage: matrix(QQ, 1, [0])._det_flint() 0 """ cdef fmpq_t e
def denominator(self): """ Return the denominator of this matrix.
OUTPUT: a Sage Integer
EXAMPLES::
sage: b = matrix(QQ,2,range(6)); b[0,0]=-5007/293; b [-5007/293 1 2] [ 3 4 5] sage: b.denominator() 293
sage: matrix(QQ, 2, [1/2, 1/3, 1/4, 1/5]).denominator() 60 """ cdef fmpz_t tmp
cdef int fmpz_denom(self, fmpz_t d) except -1: cdef Py_ssize_t i, j
def _clear_denom(self): """ INPUT:
- ``self`` - a matrix
OUTPUT: D\*self, D
The product is a matrix over ZZ
EXAMPLES::
sage: a = matrix(QQ,2,[-1/6,-7,3,5/4]); a [-1/6 -7] [ 3 5/4] sage: b, d = a._clear_denom() sage: b [ -2 -84] [ 36 15] sage: d 12 sage: b == d * a True """
cdef Py_ssize_t i, j cdef Matrix_integer_dense A cdef fmpz * entry cdef fmpz_t denom
def charpoly(self, var='x', algorithm=None): """ Return the characteristic polynomial of this matrix.
.. NOTE::
The characteristic polynomial is defined as `\det(xI-A)`.
INPUT:
- ``var`` - (optional) name of the variable as a string
- ``algorithm`` -- an optional specification of an algorithm. It can be one of:
- ``None``: (default) will use flint for small dimensions and linbox otherwise
- ``'flint'``: uses flint library
- ``'linbox'``: uses linbox library
- ``'generic'``: uses Sage generic implementation
OUTPUT: a polynomial over the rational numbers.
EXAMPLES::
sage: a = matrix(QQ, 3, [4/3, 2/5, 1/5, 4, -3/2, 0, 0, -2/3, 3/4]) sage: f = a.charpoly(); f x^3 - 7/12*x^2 - 149/40*x + 97/30 sage: f(a) [0 0 0] [0 0 0] [0 0 0]
TESTS:
The cached polynomial should be independent of the ``var`` argument (:trac:`12292`). We check (indirectly) that the second call uses the cached value by noting that its result is not cached::
sage: M = MatrixSpace(QQ, 2) sage: A = M(range(0, 2^2)) sage: type(A) <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'> sage: A.charpoly('x') x^2 - 3*x - 2 sage: A.charpoly('y') y^2 - 3*y - 2 sage: A._cache['charpoly'] x^2 - 3*x - 2
Check consistency::
sage: for _ in range(100): ....: dim = randint(0, 10) ....: m = random_matrix(QQ, dim, num_bound=8, den_bound=8) ....: p_flint = m.charpoly(algorithm='flint'); m._clear_cache() ....: p_linbox = m.charpoly(algorithm='linbox'); m._clear_cache() ....: p_generic = m.charpoly(algorithm='generic') ....: assert p_flint == p_linbox == p_generic """
else:
def minpoly(self, var='x', algorithm=None): """ Return the minimal polynomial of this matrix
INPUT:
- ``var`` - (optional) the variable name as a string (default is 'x')
- ``algorithm`` - an optional specification of an algorithm. It can be one of
- ``None``: (default) will use linbox
- ``'linbox'``: uses the linbox library
- ``'generic'``: uses the generic Sage implementation
OUTPUT: a polynomial over the rationals
EXAMPLES::
sage: a = matrix(QQ, 3, [4/3, 2/5, 1/5, 4, -3/2, 0, 0, -2/3, 3/4]) sage: f = a.minpoly(); f x^3 - 7/12*x^2 - 149/40*x + 97/30 sage: a = Mat(ZZ,4)(range(16)) sage: f = a.minpoly(); f.factor() x * (x^2 - 30*x - 80) sage: f(a) == 0 True
::
sage: a = matrix(QQ, 4, [1..4^2]) sage: factor(a.minpoly()) x * (x^2 - 34*x - 80) sage: factor(a.minpoly('y')) y * (y^2 - 34*y - 80) sage: factor(a.charpoly()) x^2 * (x^2 - 34*x - 80) sage: b = matrix(QQ, 4, [-1, 2, 2, 0, 0, 4, 2, 2, 0, 0, -1, -2, 0, -4, 0, 4]) sage: a = matrix(QQ, 4, [1, 1, 0,0, 0,1,0,0, 0,0,5,0, 0,0,0,5]) sage: c = b^(-1)*a*b sage: factor(c.minpoly()) (x - 5) * (x - 1)^2 sage: factor(c.charpoly()) (x - 5)^2 * (x - 1)^2
Check consistency::
sage: for _ in range(100): ....: dim = randint(0, 10) ....: m = random_matrix(QQ, dim, num_bound=8, den_bound=8) ....: p_linbox = m.charpoly(algorithm='linbox'); m._clear_cache() ....: p_generic = m.charpoly(algorithm='generic') ....: assert p_linbox == p_generic """
elif algorithm == 'generic': g = Matrix_dense.minpoly(self, var) else: raise ValueError("no algorithm '%s'"%algorithm)
cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right): """ EXAMPLES::
sage: a = matrix(QQ, 3, range(9))/3 sage: b = matrix(QQ, 3, range(1, 10))/5 sage: a * b # indirect doctest [ 6/5 7/5 8/5] [18/5 22/5 26/5] [ 6 37/5 44/5]
sage: matrix(QQ, 2, 3) * matrix(QQ, 4, 5) Traceback (most recent call last): ... TypeError: unsupported operand parent(s) for *: 'Full MatrixSpace of 2 by 3 dense matrices over Rational Field' and 'Full MatrixSpace of 4 by 5 dense matrices over Rational Field' """
def _multiply_flint(self, Matrix_rational_dense right): r""" Multiply this matrix by ``right`` using the flint library.
EXAMPLES::
sage: n = 3 sage: a = matrix(QQ,n,range(n^2))/3 sage: b = matrix(QQ,n,range(1, n^2 + 1))/5 sage: a._multiply_flint(b) [ 6/5 7/5 8/5] [18/5 22/5 26/5] [ 6 37/5 44/5] """ # self acts on the space of right # right acts on the space of self else:
cdef Matrix_rational_dense ans
def _multiply_over_integers(self, Matrix_rational_dense right, algorithm='default'): """ Multiply this matrix by right using a multimodular algorithm and return the result.
INPUT:
- ``self`` - matrix over QQ
- ``right`` - matrix over QQ
- ``algorithm``
- 'default': use whatever is the default for A\*B when A, B are over ZZ.
- 'multimodular': use a multimodular algorithm
EXAMPLES::
sage: a = MatrixSpace(QQ,10,5)(range(50)) sage: b = MatrixSpace(QQ,5,12)([1/n for n in range(1,61)]) sage: a._multiply_over_integers(b) == a._multiply_over_integers(b, algorithm='multimodular') True
::
sage: a = MatrixSpace(QQ,3)(range(9)) sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)]) sage: c = a._multiply_over_integers(b, algorithm = 'multimodular') sage: c [ 15/28 9/20 7/18] [ 33/7 117/40 20/9] [249/28 27/5 73/18] sage: c == a._multiply_flint(b) True """ cdef Matrix_integer_dense A, B, AB cdef Matrix_rational_dense res cdef Integer D else: sig_off() raise ValueError("unknown algorithm '%s'"%algorithm) # self acts on the space of right # right acts on the space of self else:
def height(self): """ Return the height of this matrix, which is the maximum of the absolute values of all numerators and denominators of entries in this matrix.
OUTPUT: an Integer
EXAMPLES::
sage: b = matrix(QQ,2,range(6)); b[0,0]=-5007/293; b [-5007/293 1 2] [ 3 4 5] sage: b.height() 5007 """ cdef Integer z cdef fmpz_t tmp
cdef int fmpz_height(self, fmpz_t h) except -1: cdef fmpz_t x cdef int i, j
def _adjoint(self): """ Return the adjoint of this matrix.
Assumes self is a square matrix (checked in adjoint).
EXAMPLES::
sage: m = matrix(QQ,3,[1..9])/9; m [1/9 2/9 1/3] [4/9 5/9 2/3] [7/9 8/9 1] sage: m.adjoint() [-1/27 2/27 -1/27] [ 2/27 -4/27 2/27] [-1/27 2/27 -1/27] """
def _magma_init_(self, magma): """ EXAMPLES::
sage: m = matrix(QQ,2,3,[1,2/3,-3/4,1,-2/3,-45/17]) sage: m._magma_init_(magma) 'Matrix(RationalField(),2,3,StringToIntegerSequence("204 136 -153 204 -136 -540"))/204' sage: magma(m) # optional - magma [ 1 2/3 -3/4] [ 1 -2/3 -45/17] """
def prod_of_row_sums(self, cols): cdef Py_ssize_t i, c cdef fmpq_t s, pr
raise IndexError("matrix column index out of range") cdef Rational ans
def _right_kernel_matrix(self, **kwds): r""" Returns a pair that includes a matrix of basis vectors for the right kernel of ``self``.
INPUT:
- ``kwds`` - these are provided for consistency with other versions of this method. Here they are ignored as there is no optional behavior available.
OUTPUT:
Returns a pair. First item is the string 'computed-iml-rational' that identifies the nature of the basis vectors.
Second item is a matrix whose rows are a basis for the right kernel, over the rationals, as computed by the IML library. Notice that the IML library returns a matrix that is in the 'pivot' format, once the whole matrix is multiplied by -1. So the 'computed' format is very close to the 'pivot' format.
EXAMPLES::
sage: A = matrix(QQ, [ ....: [1, 0, 1, -3, 1], ....: [-5, 1, 0, 7, -3], ....: [0, -1, -4, 6, -2], ....: [4, -1, 0, -6, 2]]) sage: result = A._right_kernel_matrix() sage: result[0] 'computed-iml-rational' sage: result[1] [-1 2 -2 -1 0] [ 1 2 0 0 -1] sage: X = result[1].transpose() sage: A*X == zero_matrix(QQ, 4, 2) True
Computed result is the negative of the pivot basis, which is just slightly more efficient to compute. ::
sage: A.right_kernel_matrix(basis='pivot') == -A.right_kernel_matrix(basis='computed') True
TESTS:
We test three trivial cases. ::
sage: A = matrix(QQ, 0, 2) sage: A._right_kernel_matrix()[1] [1 0] [0 1] sage: A = matrix(QQ, 2, 0) sage: A._right_kernel_matrix()[1].parent() Full MatrixSpace of 0 by 0 dense matrices over Rational Field sage: A = zero_matrix(QQ, 4, 3) sage: A._right_kernel_matrix()[1] [1 0 0] [0 1 0] [0 0 1] """ # _rational_kernel_flint() gets the zero-row case wrong, fix it there else:
################################################ # Change ring ################################################ def change_ring(self, R): """ Create the matrix over R with entries the entries of self coerced into R.
EXAMPLES::
sage: a = matrix(QQ,2,[1/2,-1,2,3]) sage: a.change_ring(GF(3)) [2 2] [2 0] sage: a.change_ring(ZZ) Traceback (most recent call last): ... TypeError: matrix has denominators so can't change to ZZ. sage: b = a.change_ring(QQ['x']); b [1/2 -1] [ 2 3] sage: b.parent() Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in x over Rational Field
TESTS:
Make sure that subdivisions are preserved when changing rings::
sage: a = matrix(QQ, 3, range(9)) sage: a.subdivide(2,1); a [0|1 2] [3|4 5] [-+---] [6|7 8] sage: a.change_ring(ZZ).change_ring(QQ) [0|1 2] [3|4 5] [-+---] [6|7 8] sage: a.change_ring(GF(3)) [0|1 2] [0|1 2] [-+---] [0|1 2] """ raise TypeError("R must be a ring")
# fallback to the generic version
################################################ # Echelon form ################################################ def echelonize(self, algorithm=None, height_guess=None, proof=None, **kwds): """ Transform the matrix ``self`` into reduced row echelon form in place.
INPUT:
- ``algorithm`` -- an optional specification of an algorithm. One of
- ``None``: (default) uses flint for small dimension and multimodular otherwise
- ``'flint'``: use the flint library,
- ``'padic'``: an algorithm based on the IML p-adic solver,
- ``'multimodular'``: uses a multimodular algorithm the uses linbox modulo many primes (likely to be faster when coefficients are huge),
- ``'classical'``: just clear each column using Gauss elimination.
- ``height_guess``, ``**kwds`` - all passed to the multimodular algorithm; ignored by other algorithms.
- ``proof`` - bool or None (default: None, see proof.linear_algebra or sage.structure.proof). Passed to the multimodular algorithm. Note that the Sage global default is ``proof=True``.
EXAMPLES::
sage: a = matrix(QQ, 4, range(16)); a[0,0] = 1/19; a[0,1] = 1/5; a [1/19 1/5 2 3] [ 4 5 6 7] [ 8 9 10 11] [ 12 13 14 15] sage: a.echelonize() sage: a [ 1 0 0 -76/157] [ 0 1 0 -5/157] [ 0 0 1 238/157] [ 0 0 0 0]
::
sage: a = matrix(QQ, 4, range(16)); a[0,0] = 1/19; a[0,1] = 1/5 sage: a.echelonize(algorithm='multimodular') sage: a [ 1 0 0 -76/157] [ 0 1 0 -5/157] [ 0 0 1 238/157] [ 0 0 0 0]
TESTS:
Echelonizing a matrix in place throws away the cache of the old matrix (:trac:`14506`)::
sage: for algo in ["flint", "padic", "multimodular", "classical"]: ....: a = Matrix(QQ, [[1,2],[3,4]]) ....: _ = a.det() # fills the cache ....: _ = a._clear_denom() # fills the cache ....: a.echelonize(algorithm=algo) ....: assert sorted(a._cache.keys()) == ['echelon_form', 'in_echelon_form', 'pivots', 'rank'], (algo, a._cache.keys()) """
else:
else: raise ValueError("no algorithm '%s'"%algorithm)
raise RuntimeError("BUG: pivots must get set as a tuple. Got {} for algo {} with {}x{} matrix.".format( type(pivots), algorithm, self._nrows, self._ncols))
def echelon_form(self, algorithm=None, height_guess=None, proof=None, **kwds): r""" Return the echelon form of this matrix.
The (row) echelon form of a matrix, see :wikipedia:`Row_echelon_form`, is the matrix obtained by performing Gauss elimination on the rows of the matrix.
INPUT: See :meth:`echelonize` for the options.
EXAMPLES::
sage: a = matrix(QQ, 4, range(16)); a[0,0] = 1/19; a[0,1] = 1/5; a [1/19 1/5 2 3] [ 4 5 6 7] [ 8 9 10 11] [ 12 13 14 15] sage: a.echelon_form() [ 1 0 0 -76/157] [ 0 1 0 -5/157] [ 0 0 1 238/157] [ 0 0 0 0] sage: a.echelon_form(algorithm='multimodular') [ 1 0 0 -76/157] [ 0 1 0 -5/157] [ 0 0 1 238/157] [ 0 0 0 0]
The result is an immutable matrix, so if you want to modify the result then you need to make a copy. This checks that :trac:`10543` is fixed.::
sage: A = matrix(QQ, 2, range(6)) sage: E = A.echelon_form() sage: E.is_mutable() False sage: F = copy(E) sage: F[0,0] = 50 sage: F [50 0 -1] [ 0 1 2]
TESTS:
Check consistency::
sage: for _ in range(100): ....: nrows = randint(0, 30) ....: ncols = randint(0, 30) ....: m = random_matrix(QQ, nrows, ncols, num_bound=10, den_bound=10) ....: ech_flint = m.echelon_form('flint'); m._clear_cache() ....: ech_padic = m.echelon_form('padic'); m._clear_cache() ....: ech_multi = m.echelon_form('multimodular'); m._clear_cache() ....: ech_class = m.echelon_form('classical') ....: assert ech_flint == ech_padic == ech_multi == ech_class """ raise RuntimeError('in_echelon_form set but not echelon_form')
def _echelonize_flint(self): r""" EXAMPLES::
sage: m = matrix(QQ, 4, range(16)) sage: m._echelonize_flint() (0, 1) sage: m [ 1 0 -1 -2] [ 0 1 2 3] [ 0 0 0 0] [ 0 0 0 0] sage: m = matrix(QQ, 4, 6, [-1,0,0,-2,-1,-2,-1,0,0,-2,-1,0,3,3,-2,0,0,3,-2,-3,1,1,-2,3]) sage: m._echelonize_flint() (0, 1, 2, 5) sage: m [ 1 0 0 2 1 0] [ 0 1 0 -4/3 1 0] [ 0 0 1 1 3 0] [ 0 0 0 0 0 1] """ cdef long r
# compute pivots cdef long i, j, k else: break
def _echelonize_padic(self): """ Echelonize self using a p-adic nullspace algorithm.
EXAMPLES::
sage: m = matrix(QQ, 4, range(16)) sage: m._echelonize_padic() (0, 1) sage: m [ 1 0 -1 -2] [ 0 1 2 3] [ 0 0 0 0] [ 0 0 0 0]
sage: m = matrix(QQ, 4, 6, [-1,0,0,-2,-1,-2,-1,0,0,-2,-1,0,3,3,-2,0,0,3,-2,-3,1,1,-2,3]) sage: m._echelonize_padic() (0, 1, 2, 5) sage: m [ 1 0 0 2 1 0] [ 0 1 0 -4/3 1 0] [ 0 0 1 1 3 0] [ 0 0 0 0 0 1] """ cdef Matrix_integer_dense X cdef Integer d cdef fmpq * entry
# FIXME: we should always have X.nrows() == len(pivots)
cdef Py_ssize_t i,j # 1 at pivot
# nonzero part
# zeros on the left of the pivot
# zeros on top of the other pivots
# Fill in the 0-rows at the bottom.
# FIXME: pivots should already be a tuple in all cases
def _echelonize_multimodular(self, height_guess=None, proof=None): """ Echelonize ``self`` using multimodular recomposition.
REFERENCE:
- Chapter 7 of Stein's "Explicitly Computing Modular Forms".
INPUT:
- ``height_guess`` - integer or None
- ``proof`` - boolean (default: None, see proof.linear_algebra or sage.structure.proof) Note that the Sage global default is proof=True.
EXAMPLES::
sage: m = matrix(QQ, 4, range(16)) sage: m._echelonize_multimodular() (0, 1) sage: m [ 1 0 -1 -2] [ 0 1 2 3] [ 0 0 0 0] [ 0 0 0 0] sage: m = matrix(QQ, 4, 6, [-1,0,0,-2,-1,-2,-1,0,0,-2,-1,0,3,3,-2,0,0,3,-2,-3,1,1,-2,3]) sage: m._echelonize_multimodular() (0, 1, 2, 5) sage: m [ 1 0 0 2 1 0] [ 0 1 0 -4/3 1 0] [ 0 0 1 1 3 0] [ 0 0 0 0 0 1] """
cdef swap_rows_c(self, Py_ssize_t r1, Py_ssize_t r2): """ EXAMPLES::
sage: a = matrix(QQ,2,[1..6]) sage: a.swap_rows(0,1) # indirect doctest sage: a [4 5 6] [1 2 3] """ # no bounds checking! cdef Py_ssize_t c fmpq_mat_entry(self._matrix, r2, c))
cdef swap_columns_c(self, Py_ssize_t c1, Py_ssize_t c2): """ EXAMPLES::
sage: a = matrix(QQ,2,[1..6]) sage: a.swap_columns(0,1) # indirect doctest sage: a [2 1 3] [5 4 6] """ # no bounds checking! fmpq_mat_entry(self._matrix, r, c2))
def decomposition(self, is_diagonalizable=False, dual=False, algorithm=None, height_guess=None, proof=None): """ Returns the decomposition of the free module on which this matrix A acts from the right (i.e., the action is x goes to x A), along with whether this matrix acts irreducibly on each factor. The factors are guaranteed to be sorted in the same way as the corresponding factors of the characteristic polynomial.
Let A be the matrix acting from the on the vector space V of column vectors. Assume that A is square. This function computes maximal subspaces W_1, ..., W_n corresponding to Galois conjugacy classes of eigenvalues of A. More precisely, let f(X) be the characteristic polynomial of A. This function computes the subspace `W_i = ker(g_(A)^n)`, where g_i(X) is an irreducible factor of f(X) and g_i(X) exactly divides f(X). If the optional parameter is_diagonalizable is True, then we let W_i = ker(g(A)), since then we know that ker(g(A)) = `ker(g(A)^n)`.
If dual is True, also returns the corresponding decomposition of V under the action of the transpose of A. The factors are guaranteed to correspond.
INPUT:
- ``is_diagonalizable`` - ignored
- ``dual`` - whether to also return decompositions for the dual
- ``algorithm`` - an optional specification of an algorithm
- ``None`` - (default) use default algorithm for computing Echelon forms
- 'multimodular': much better if the answers factors have small height
- ``height_guess`` - positive integer; only used by the multimodular algorithm
- ``proof`` - bool or None (default: None, see proof.linear_algebra or sage.structure.proof); only used by the multimodular algorithm. Note that the Sage global default is proof=True.
.. NOTE::
IMPORTANT: If you expect that the subspaces in the answer are spanned by vectors with small height coordinates, use algorithm='multimodular' and height_guess=1; this is potentially much faster than the default. If you know for a fact the answer will be very small, use algorithm='multimodular', height_guess=bound on height, proof=False.
You can get very very fast decomposition with proof=False.
EXAMPLES::
sage: a = matrix(QQ,3,[1..9]) sage: a.decomposition() [ (Vector space of degree 3 and dimension 1 over Rational Field Basis matrix: [ 1 -2 1], True), (Vector space of degree 3 and dimension 2 over Rational Field Basis matrix: [ 1 0 -1] [ 0 1 2], True) ]
""" Y = self.transpose()._decomposition_rational(is_diagonalizable=is_diagonalizable, echelon_algorithm = algorithm, height_guess = height_guess, proof=proof) return X, Y
def _decomposition_rational(self, is_diagonalizable = False, echelon_algorithm=None, kernel_algorithm='default', **kwds): """ Returns the decomposition of the free module on which this matrix A acts from the right (i.e., the action is x goes to x A), along with whether this matrix acts irreducibly on each factor. The factors are guaranteed to be sorted in the same way as the corresponding factors of the characteristic polynomial.
INPUT:
- ``self`` - a square matrix over the rational numbers
- ``echelon_algorithm`` - an optional algorithm to be passed to the method ``echelon_form``
- ``'multimodular'`` - use this if the answers have small height
- ``**kwds`` - passed on to echelon function.
.. NOTE::
IMPORTANT: If you expect that the subspaces in the answer are spanned by vectors with small height coordinates, use algorithm='multimodular' and height_guess=1; this is potentially much faster than the default. If you know for a fact the answer will be very small, use algorithm='multimodular', height_guess=bound on height, proof=False
OUTPUT:
- ``Sequence`` - list of tuples (V,t), where V is a vector spaces and t is True if and only if the charpoly of self on V is irreducible. The tuples are in order corresponding to the elements of the sorted list self.charpoly().factor(). """ cdef Py_ssize_t k
raise ArithmeticError("self must be a square matrix")
# Just use kernel -- much easier.
# General case, i.e., deg(g) > 1:
# Compute the complementary factor of the charpoly.
# Compute one element of the kernel of g(A)**m. caller_name='rational decomp')
# Get the rest of the kernel. else:
else: verbose('we have not yet generated all the kernel (rank so far=%s, target rank=%s)'%( W.rank(), m*g.degree()), level=2, caller_name='rational decomp') tries += 1 if tries > 5*m: raise RuntimeError("likely bug in decomposition") # end if #end while #end for
## def simple_decomposition(self, echelon_algorithm='default', **kwds): ## """ ## Returns the decomposition of the free module on which this ## matrix A acts from the right (i.e., the action is x goes to x ## A), as a direct sum of simple modules.
## NOTE: self *must* be diagonalizable.
## INPUT: ## self -- a square matrix that is assumed to be diagonalizable ## echelon_algorithm -- 'default' ## 'multimodular' -- use this if the answers ## have small height ## **kwds -- passed on to echelon function.
## IMPORTANT NOTE: ## If you expect that the subspaces in the answer are spanned by vectors ## with small height coordinates, use algorithm='multimodular' and ## height_guess=1; this is potentially much faster than the default. ## If you know for a fact the answer will be very small, use ## algorithm='multimodular', height_guess=bound on height, proof=False
## OUTPUT: ## Sequence -- list of tuples (V,g), where V is a subspace ## and an irreducible polynomial g, which is the ## charpoly (=minpoly) of self acting on V. ## """ ## cdef Py_ssize_t k
## if not self.is_square(): ## raise ArithmeticError("self must be a square matrix")
## if self.nrows() == 0: ## return decomp_seq([])
## A, _ = self._clear_denom()
## f = A.charpoly('x') ## E = decomp_seq([])
## t = verbose('factoring the characteristic polynomial', level=2, caller_name='simple decomp') ## F = f.factor() ## G = [g for g, _ in F] ## minpoly = prod(G) ## squarefree_degree = sum([g.degree() for g in G]) ## verbose('done factoring', t=t, level=2, caller_name='simple decomp')
## V = ZZ**self.nrows() ## v = V.random_element() ## num_iterates = max([squarefree_degree - g.degree() for g in G]) + 1
## S = [ ]
## F.sort() ## for i in range(len(F)): ## g, m = F[i]
## if g.degree() == 1: ## # Just use kernel -- much easier. ## B = A.__copy__() ## for k from 0 <= k < A.nrows(): ## B[k,k] += g[0] ## if m > 1 and not is_diagonalizable: ## B = B**m ## W = B.change_ring(QQ).kernel() ## for b in W.basis(): ## E.append((W.span(b), g)) ## continue
## # General case, i.e., deg(g) > 1: ## W = None ## while True:
## # Compute the complementary factor of the charpoly. ## h = minpoly // g ## v = h.list()
## while len(S) < m: ## t = verbose('%s-spinning %s-th random vector'%(num_iterates, len(S)), ## level=2, caller_name='simple decomp') ## S.append(A.iterates(V.random_element(x=-10,y=10), num_iterates)) ## verbose('done spinning', level=2, t=t, caller_name='simple decomp')
## for j in range(len(S)): ## # Compute one element of the kernel of g(A). ## t = verbose('compute element of kernel of g(A), for g of degree %s'%g.degree(),level=2, ## caller_name='simple decomp') ## w = S[j].linear_combination_of_rows(h.list()) ## t = verbose('done computing element of kernel of g(A)', t=t,level=2, caller_name='simple decomp')
## # Get the rest of the kernel. ## t = verbose('fill out rest of kernel',level=2, caller_name='simple decomp') ## if W is None: ## W = A.iterates(w, g.degree()) ## else: ## W = W.stack(A.iterates(w, g.degree())) ## t = verbose('finished filling out more of kernel',level=2, t=t, caller_name='simple decomp')
## if W.rank() == m * g.degree(): ## W = W.change_ring(QQ) ## t = verbose('now computing row space', level=2, caller_name='simple decomp') ## W.echelonize(algorithm = echelon_algorithm, **kwds) ## E.append((W.row_space(), m==1)) ## verbose('computed row space', level=2,t=t, caller_name='simple decomp') ## break ## else: ## verbose('we have not yet generated all the kernel (rank so far=%s, target rank=%s)'%( ## W.rank(), m*g.degree()), level=2, caller_name='simple decomp') ## j += 1 ## if j > 3*m: ## raise RuntimeError("likely bug in decomposition") ## # end if ## #end while ## #end for ## return E
def _lift_crt_rr(self, res, mm): cdef Integer m cdef Matrix_integer_dense ZA cdef Matrix_rational_dense QA cdef Py_ssize_t i, j cdef mpz_t* Z_row cdef mpq_t* Q_row cdef mpz_t tmp cdef mpq_t tmp2 mpz_init(tmp) mpq_init(tmp2) ZA = _lift_crt(res, mm) QA = Matrix_rational_dense.__new__(Matrix_rational_dense, self.parent(), None, None, None) m = mm.prod() for i in range(ZA._nrows): for j in range(ZA._ncols): fmpz_get_mpz(tmp, fmpz_mat_entry(ZA._matrix,i,j)) mpq_rational_reconstruction(tmp2, tmp, m.value) fmpq_set_mpq(fmpq_mat_entry(QA._matrix, i, j), tmp2) mpz_clear(tmp) mpq_clear(tmp2) return QA
def randomize(self, density=1, num_bound=2, den_bound=2, \ distribution=None, nonzero=False): """ Randomize ``density`` proportion of the entries of this matrix, leaving the rest unchanged.
If ``x`` and ``y`` are given, randomized entries of this matrix have numerators and denominators bounded by ``x`` and ``y`` and have density 1.
INPUT:
- ``density`` - number between 0 and 1 (default: 1)
- ``num_bound`` - numerator bound (default: 2)
- ``den_bound`` - denominator bound (default: 2)
- ``distribution`` - ``None`` or '1/n' (default: ``None``); if '1/n' then ``num_bound``, ``den_bound`` are ignored and numbers are chosen using the GMP function ``mpq_randomize_entry_recip_uniform``
OUTPUT:
- None, the matrix is modified in-space
EXAMPLES::
sage: a = matrix(QQ,2,4); a.randomize(); a [ 0 -1 2 -2] [ 1 -1 2 1] sage: a = matrix(QQ,2,4); a.randomize(density=0.5); a [ -1 -2 0 0] [ 0 0 1/2 0] sage: a = matrix(QQ,2,4); a.randomize(num_bound=100, den_bound=100); a [ 14/27 21/25 43/42 -48/67] [-19/55 64/67 -11/51 76] sage: a = matrix(QQ,2,4); a.randomize(distribution='1/n'); a [ 3 1/9 1/2 1/4] [ 1 1/39 2 -1955/2]
TESTS:
Check that the option ``nonzero`` is meaningful (:trac:`22970`)::
sage: a = matrix(QQ, 10, 10, 1) sage: b = a.__copy__() sage: b.randomize(nonzero=True) sage: a == b False sage: any(b[i,j].is_zero() for i in range(10) for j in range(10)) False """ return
cdef Integer B, C cdef Py_ssize_t i, j, k, num_per_row cdef randstate rstate cdef mpq_t tmp
else: sig_on() for i in range(self._nrows): for j in range(self._ncols): mpq_randomize_entry_as_int(tmp, B.value) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, j), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(num_per_row): k = rstate.c_random() % self._ncols mpq_randomize_entry_recip_uniform(tmp) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, k), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(num_per_row): k = rstate.c_random() % self._ncols mpq_randomize_entry_as_int(tmp, B.value) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, k), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(self._ncols): mpq_randomize_entry_recip_uniform_nonzero(tmp) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, j), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(self._ncols): mpq_randomize_entry_as_int_nonzero(tmp, B.value) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, j), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(num_per_row): k = rstate.c_random() % self._ncols mpq_randomize_entry_recip_uniform_nonzero(tmp) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, k), tmp) sig_off() else: sig_on() for i in range(self._nrows): for j in range(num_per_row): k = rstate.c_random() % self._ncols mpq_randomize_entry_as_int_nonzero(tmp, B.value) fmpq_set_mpq(fmpq_mat_entry(self._matrix, i, k), tmp) sig_off()
def rank(self, algorithm=None): """ Return the rank of this matrix.
INPUT:
- ``algorithm`` - an optional specification of an algorithm. One of
- ``None``: (default) will use flint
- ``'flint'``: uses the flint library
- ``'pari'``: uses the PARI library
- ``'integer'``: eliminate denominators and calls the rank function on the corresponding integer matrix
EXAMPLES::
sage: matrix(QQ,3,[1..9]).rank() 2 sage: matrix(QQ,100,[1..100^2]).rank() 2
TESTS::
sage: for _ in range(100): ....: dim = randint(0, 30) ....: m = random_matrix(QQ, dim, num_bound=2, density=0.5) ....: r_pari = m.rank('pari'); m._clear_cache() ....: r_flint = m.rank('flint'); m._clear_cache() ....: r_int = m.rank('integer'); m._clear_cache() ....: assert r_pari == r_flint == r_int """
else: raise ValueError("unknown algorithm %s" % algorithm)
def transpose(self): """ Returns the transpose of self, without changing self.
EXAMPLES:
We create a matrix, compute its transpose, and note that the original matrix is not changed.
::
sage: A = matrix(QQ, 2, 3, range(6)) sage: type(A) <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'> sage: B = A.transpose() sage: print(B) [0 3] [1 4] [2 5] sage: print(A) [0 1 2] [3 4 5]
``.T`` is a convenient shortcut for the transpose::
sage: print(A.T) [0 3] [1 4] [2 5]
::
sage: A.subdivide(None, 1); A [0|1 2] [3|4 5] sage: A.transpose() [0 3] [---] [1 4] [2 5] """ cdef Matrix_rational_dense ans else:
def antitranspose(self): """ Returns the antitranspose of self, without changing self.
EXAMPLES::
sage: A = matrix(QQ,2,3,range(6)) sage: type(A) <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'> sage: A.antitranspose() [5 2] [4 1] [3 0] sage: A [0 1 2] [3 4 5]
sage: A.subdivide(1,2); A [0 1|2] [---+-] [3 4|5] sage: A.antitranspose() [5|2] [-+-] [4|1] [3|0] """ parent = self._parent else:
cdef Matrix_rational_dense ans
cdef Py_ssize_t i,j cdef Py_ssize_t ri,rj # reversed i and j fmpq_mat_entry(self._matrix, i, j))
def set_row_to_multiple_of_row(self, Py_ssize_t i, Py_ssize_t j, s): """ Set row i equal to s times row j.
EXAMPLES::
sage: a = matrix(QQ,2,3,range(6)); a [0 1 2] [3 4 5] sage: a.set_row_to_multiple_of_row(1,0,-3) sage: a [ 0 1 2] [ 0 -3 -6] """ cdef Py_ssize_t k cdef fmpq_t ss fmpq_mat_entry(self._matrix, j, k), ss)
def _set_row_to_negative_of_row_of_A_using_subset_of_columns(self, Py_ssize_t i, Matrix A, Py_ssize_t r, cols, cols_index=None): """ Set row i of self to -(row r of A), but where we only take the given column positions in that row of A. We do not zero out the other entries of self's row i either.
.. NOTE::
This function exists just because it is useful for modular symbols presentations.
INPUT:
- ``i`` - integer, index into the rows of self
- ``A`` - a matrix
- ``r`` - integer, index into rows of A
- ``cols`` - a *sorted* list of integers.
EXAMPLES::
sage: a = matrix(QQ,2,3,range(6)); a [0 1 2] [3 4 5] sage: a._set_row_to_negative_of_row_of_A_using_subset_of_columns(0,a,1,[1,2]) sage: a [-4 -5 2] [ 3 4 5] """ cdef Matrix_rational_dense _A raise IndexError("invalid row")
A = A.change_ring(QQ) A = A.dense_matrix()
def _add_col_j_of_A_to_col_i_of_self(self, Py_ssize_t i, Matrix_rational_dense A, Py_ssize_t j): """ Unsafe technical function that very quickly adds the j-th column of A to the i-th column of self.
Does not check mutability. """ raise TypeError("nrows of self and A must be the same") cdef Py_ssize_t r fmpq_mat_entry(self._matrix, r, i), fmpq_mat_entry(A._matrix, r, j))
################################################# # Methods using PARI library # #################################################
def __pari__(self): """ Return pari version of this matrix.
EXAMPLES::
sage: matrix(QQ,2,[1/5,-2/3,3/4,4/9]).__pari__() [1/5, -2/3; 3/4, 4/9] """
def _det_pari(self, int flag=0): """ Return the determinant of this matrix computed using pari.
EXAMPLES:: sage: matrix(QQ,3,[1..9])._det_pari() 0 sage: matrix(QQ,3,[1..9])._det_pari(1) 0 sage: matrix(QQ,3,[0]+[2..9])._det_pari() 3 """ # now convert d to a Sage rational
def _rank_pari(self): """ Return the rank of this matrix computed using pari.
EXAMPLES::
sage: matrix(QQ,3,[1..9])._rank_pari() 2 sage: matrix(QQ, 0, 0)._rank_pari() 0 """
def _multiply_pari(self, Matrix_rational_dense right): """ Return the product of self and right, computed using PARI.
EXAMPLES::
sage: matrix(QQ,2,[1/5,-2/3,3/4,4/9])._multiply_pari(matrix(QQ,2,[1,2,3,4])) [ -9/5 -34/15] [ 25/12 59/18]
We verify that 0 rows or columns works::
sage: x = matrix(QQ,2,0); y= matrix(QQ,0,2); x*y [0 0] [0 0] sage: matrix(ZZ, 0, 0) * matrix(QQ, 0, 5) [] """ raise ArithmeticError("self must be a square matrix") # pari doesn't work in case of 0 rows or columns # This case is easy, since the answer must be the 0 matrix. return self.matrix_space(self._nrows, right._ncols).zero_matrix().__copy__() _new_GEN_from_fmpq_mat_t(right._matrix))
def _invert_pari(self): """ Return the inverse of this matrix computed using PARI.
EXAMPLES::
sage: matrix(QQ,2,[1,2,3,4])._invert_pari() [ -2 1] [ 3/2 -1/2] sage: matrix(QQ,2,[1,2,2,4])._invert_pari() Traceback (most recent call last): ... PariError: impossible inverse in ginv: [1, 2; 2, 4] """ raise ValueError("self must be a square matrix") cdef GEN M, d
# Convert matrix back to Sage.
def row(self, Py_ssize_t i, from_list=False): """ Return the i-th row of this matrix as a dense vector.
INPUT:
- ``i`` - integer
- ``from_list`` - ignored
EXAMPLES::
sage: m = matrix(QQ, 2, [1/5, -2/3, 3/4, 4/9]) sage: m.row(0) (1/5, -2/3) sage: m.row(1) (3/4, 4/9) sage: m.row(1, from_list=True) (3/4, 4/9) sage: m.row(-2) (1/5, -2/3)
sage: m.row(2) Traceback (most recent call last): ... IndexError: row index out of range sage: m.row(-3) Traceback (most recent call last): ... IndexError: row index out of range """
cdef Py_ssize_t j
def column(self, Py_ssize_t i, from_list=False): """ Return the i-th column of this matrix as a dense vector.
INPUT:
- ``i`` - integer
- ``from_list`` - ignored
EXAMPLES::
sage: m = matrix(QQ, 3, 2, [1/5,-2/3,3/4,4/9,-1,0]) sage: m.column(1) (-2/3, 4/9, 0) sage: m.column(1,from_list=True) (-2/3, 4/9, 0) sage: m.column(-1) (-2/3, 4/9, 0) sage: m.column(-2) (1/5, 3/4, -1)
sage: m.column(2) Traceback (most recent call last): ... IndexError: column index out of range sage: m.column(-3) Traceback (most recent call last): ... IndexError: column index out of range """
cdef Py_ssize_t j
################################################ # LLL ################################################
def LLL(self, *args, **kwargs): """ Return an LLL reduced or approximated LLL reduced lattice for ``self`` interpreted as a lattice.
For details on input parameters, see :meth:`sage.matrix.matrix_integer_dense.Matrix_integer_dense.LLL`.
EXAMPLES::
sage: A = Matrix(QQ, 3, 3, [1/n for n in range(1, 10)]) sage: A.LLL() [ 1/28 -1/40 -1/18] [ 1/28 -1/40 1/18] [ 0 -3/40 0] """
cdef new_matrix_from_pari_GEN(parent, GEN d): """ Given a PARI GEN with ``t_INT`` or ``t_FRAC entries, create a :class:`Matrix_rational_dense` from it.
EXAMPLES::
sage: matrix(QQ,2,[1..4])._multiply_pari(matrix(QQ,2,[2..5])) # indirect doctest [10 13] [22 29] """ cdef Py_ssize_t i, j Matrix_rational_dense, parent, None, None, None) cdef mpq_t tmp |