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
# -*- coding: utf-8 -*- """ Dense matrices over the integer ring
AUTHORS:
- William Stein
- Robert Bradshaw
- Marc Masdeu (August 2014). Implemented using FLINT, see :trac:`16803`.
- Jeroen Demeyer (October 2014): lots of fixes, see :trac:`17090` and :trac:`17094`.
- Vincent Delecroix (February 2015): make it faster, see :trac:`17822`.
- Vincent Delecroix (May 2017): removed duplication of entries and cleaner linbox interface
EXAMPLES::
sage: a = matrix(ZZ, 3,3, range(9)); a [0 1 2] [3 4 5] [6 7 8] sage: a.det() 0 sage: a[0,0] = 10; a.det() -30 sage: a.charpoly() x^3 - 22*x^2 + 102*x + 30 sage: b = -3*a sage: a == b False sage: b < a True
TESTS::
sage: a = matrix(ZZ,2,range(4), sparse=False) sage: TestSuite(a).run() sage: Matrix(ZZ,0,0).inverse() []
"""
#***************************************************************************** # Copyright (C) 2006,2007 William Stein # Copyright (C) 2014 Marc Masdeu # Copyright (C) 2014 Jeroen Demeyer # Copyright (C) 2015,2016,2017 Vincent Delecroix # # 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.stdint cimport int64_t from libc.string cimport strcpy, strlen
from sage.cpython.string cimport char_to_str, str_to_bytes from sage.ext.stdsage cimport PY_NEW from cysignals.signals cimport sig_check, sig_on, sig_str, sig_off from cysignals.memory cimport sig_malloc, sig_free, check_allocarray
from sage.libs.gmp.mpz cimport *
from sage.modules.vector_integer_dense cimport Vector_integer_dense
from sage.misc.misc import verbose, get_verbose, cputime
from sage.arith.all import previous_prime from sage.arith.power cimport generic_power from sage.structure.element cimport Element from sage.structure.proof.proof import get_flag as get_proof_flag from sage.misc.randstate cimport randstate, current_randstate
from sage.matrix.matrix_rational_dense cimport Matrix_rational_dense
######################################################### # PARI C library from cypari2.gen cimport Gen from sage.libs.pari.convert_gmp cimport INT_to_mpz from sage.libs.pari.convert_flint cimport (_new_GEN_from_fmpz_mat_t, _new_GEN_from_fmpz_mat_t_rotate90, integer_matrix) from cypari2.stack cimport clear_stack from cypari2.paridecl cimport * #########################################################
from sage.arith.multi_modular cimport MultiModularBasis
from sage.libs.flint.fmpz cimport * from sage.libs.flint.fmpz_mat cimport *
from sage.rings.integer cimport Integer from sage.rings.rational_field import QQ from sage.rings.real_double import RDF from sage.rings.integer_ring import ZZ, IntegerRing_class 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 from sage.structure.element import is_Vector from sage.structure.sequence import Sequence
from .matrix_modn_dense_float cimport Matrix_modn_dense_template from .matrix_modn_dense_float cimport Matrix_modn_dense_float from .matrix_modn_dense_double cimport Matrix_modn_dense_double
from .matrix_mod2_dense import Matrix_mod2_dense from .matrix_mod2_dense cimport Matrix_mod2_dense
from .matrix2 import decomp_seq
from .matrix cimport Matrix
cimport sage.structure.element
import sage.matrix.matrix_space as matrix_space
################ # Used for modular HNF from sage.rings.fast_arith cimport arith_int cdef arith_int ai = arith_int()
######### linbox interface ########## from sage.libs.linbox.linbox_flint_interface cimport *
########## iml -- integer matrix library ########### from sage.libs.iml cimport *
fplll_fp_map = {None: None, 'fp': 'double', 'ld': 'long double', 'xd': 'dpe', 'rr': 'mpfr'}
cdef inline mpz_t * fmpz_mat_to_mpz_array(fmpz_mat_t m) except? NULL: cdef size_t i, j
cdef inline void mpz_array_clear(mpz_t * a, size_t length): cdef size_t i
cdef class Matrix_integer_dense(Matrix_dense): r""" Matrix over the integers, implemented using FLINT.
On a 32-bit machine, they can have at most `2^{32}-1` rows or columns. On a 64-bit machine, matrices can have at most `2^{64}-1` rows or columns.
EXAMPLES::
sage: a = MatrixSpace(ZZ,3)(2); a [2 0 0] [0 2 0] [0 0 2] sage: a = matrix(ZZ,1,3, [1,2,-3]); a [ 1 2 -3] sage: a = MatrixSpace(ZZ,2,4)(2); a Traceback (most recent call last): ... TypeError: nonzero scalar matrix must be square
TESTS:
Test hashing::
sage: a = Matrix(ZZ, 2, [1,2,3,4]) sage: hash(a) Traceback (most recent call last): ... TypeError: mutable matrices are unhashable sage: a.set_immutable() sage: hash(a) 1846857684291126914 # 64-bit 1591707266 # 32-bit """ def __cinit__(self, parent, entries, coerce, copy): """ Create and allocate memory for the matrix. Does not actually initialize any of the memory.
INPUT:
- ``parent, entries, coerce, copy`` - as for __init__.
EXAMPLES::
sage: from sage.matrix.matrix_integer_dense import Matrix_integer_dense sage: a = Matrix_integer_dense.__new__(Matrix_integer_dense, Mat(ZZ,3), 0,0,0) sage: type(a) <type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
TESTS::
sage: Matrix(ZZ, sys.maxsize, sys.maxsize) Traceback (most recent call last): ... RuntimeError: FLINT exception """
def __dealloc__(self): """ Frees all the memory allocated for this matrix.
EXAMPLES::
sage: a = Matrix(ZZ,2,[1,2,3,4]) sage: del a """
def __init__(self, parent, entries, copy, coerce): r""" Initialize a dense matrix over the integers.
INPUT:
- ``parent`` - a matrix space
- ``entries`` - list - create the matrix with those entries along the rows.
- ``other`` - a scalar; entries is coerced to an integer and the diagonal entries of this matrix are set to that integer.
- ``coerce`` - whether need to coerce entries to the integers (program may crash if you get this wrong)
- ``copy`` - ignored (since integers are immutable)
EXAMPLES:
The __init__ function is called implicitly in each of the examples below to actually fill in the values of the matrix.
We create a `2 \times 2` and a `1\times 4` matrix::
sage: matrix(ZZ,2,2,range(4)) [0 1] [2 3] sage: Matrix(ZZ,1,4,range(4)) [0 1 2 3]
If the number of columns isn't given, it is determined from the number of elements in the list.
::
sage: matrix(ZZ,2,range(4)) [0 1] [2 3] sage: matrix(ZZ,2,range(6)) [0 1 2] [3 4 5]
Another way to make a matrix is to create the space of matrices and coerce lists into it.
::
sage: A = Mat(ZZ,2); A Full MatrixSpace of 2 by 2 dense matrices over Integer Ring sage: A(range(4)) [0 1] [2 3]
Actually it is only necessary that the input can be coerced to a list, so the following also works::
sage: v = reversed(range(4)) sage: A(v) [3 2] [1 0]
Matrices can have many rows or columns (in fact, on a 64-bit machine they could have up to `2^64-1` rows or columns)::
sage: v = matrix(ZZ,1,10^5, range(10^5)) sage: v.parent() Full MatrixSpace of 1 by 100000 dense matrices over Integer Ring """ cdef Py_ssize_t i, j, k cdef bint is_list cdef Integer x cdef list entries_list
else: # Create the matrix whose entries are in the given entry list. # todo -- see integer.pyx and the TODO there; perhaps this could be # sped up by creating a mpz_init_set_sage function. else: else: # If x is zero, make the zero matrix and be done.
# the matrix must be square:
# Now we set all the diagonal entries to x and all other entries to 0.
cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, object x): """ Set position i,j of this matrix to ``x``.
The object ``x`` must be of type ``Integer``.
INPUT:
- ``i`` -- row
- ``j`` -- column
- ``x`` -- must be Integer! The value to set self[i,j] to.
EXAMPLES::
sage: a = matrix(ZZ,2,3, range(6)); a [0 1 2] [3 4 5] sage: a[0,0] = 10 sage: a [10 1 2] [ 3 4 5] """
cdef void set_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, const mpz_t value): """ Set position i,j of this matrix to ``value``.
INPUT:
- ``i`` -- row
- ``j`` -- column
- ``value`` -- The value to set self[i,j] to. This will make a copy of ``value``.
EXAMPLES::
sage: a = matrix(ZZ,2,3, range(6)); a [0 1 2] [3 4 5] sage: a[0,0] = 10 sage: a [10 1 2] [ 3 4 5] """
cdef void set_unsafe_si(self, Py_ssize_t i, Py_ssize_t j, long value): """ Set position i,j of this matrix to ``value``. """
cdef void set_unsafe_double(self, Py_ssize_t i, Py_ssize_t j, double value): """ Set position i,j of this matrix to ``value``. """
cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): """ Returns (i, j) entry of self as a new Integer.
.. warning::
This is very unsafe; it assumes i and j are in the right range.
EXAMPLES::
sage: a = MatrixSpace(ZZ,3)(range(9)); a [0 1 2] [3 4 5] [6 7 8] sage: a[1,2] 5 sage: a[4,7] Traceback (most recent call last): ... IndexError: matrix index out of range sage: a[-1,0] 6 """
cdef inline void get_unsafe_mpz(self, Py_ssize_t i, Py_ssize_t j, mpz_t value): """ Copy entry i,j of the matrix ``self`` to ``value``.
.. warning::
This is very unsafe; it assumes i and j are in the right range.
EXAMPLES::
sage: a = MatrixSpace(ZZ,3)(range(9)); a [0 1 2] [3 4 5] [6 7 8] sage: a[1,2] 5 sage: a[4,7] Traceback (most recent call last): ... IndexError: matrix index out of range sage: a[-1,0] 6 """
cdef inline double get_unsafe_double(self, Py_ssize_t i, Py_ssize_t j): """ Returns (j, i) entry of self as a new Integer.
.. warning::
This is very unsafe; it assumes i and j are in the right range.
EXAMPLES::
sage: a = MatrixSpace(ZZ,3)(range(9)); a [0 1 2] [3 4 5] [6 7 8] sage: a[1,2] 5 sage: a[4,7] Traceback (most recent call last): ... IndexError: matrix index out of range sage: a[-1,0] 6 """
def _pickle(self): """ EXAMPLES::
sage: a = matrix(ZZ,2,3,[1,193,15,-2,3,0]) sage: a._pickle() ('1 61 f -2 3 0', 0)
sage: S = ModularSymbols(250,4,sign=1).cuspidal_submodule().new_subspace().decomposition() # long time sage: S == loads(dumps(S)) # long time True """
cdef _pickle_version0(self): """ EXAMPLES::
sage: matrix(ZZ,1,3,[1,193,15])._pickle() # indirect doctest ('1 61 f', 0)
"""
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 integer = 36; (default: 10)
EXAMPLES::
sage: m = matrix(ZZ,2,3,[1,2,-3,1,-2,-45]) sage: m._export_as_string(10) '1 2 -3 1 -2 -45' sage: m._export_as_string(16) '1 2 -3 1 -2 -2d' """ # TODO: *maybe* redo this to use mpz_import and mpz_export # from sec 5.14 of the GMP manual. ?? cdef int i, j, len_so_far, m, n cdef char *a cdef char *s cdef char *t cdef char *tmp
else:
# mat_entry = fmpz_mat_entry(self._matrix,i,j) # copy to new string with double the size n = 2*n + m + 1 tmp = <char*> sig_malloc(n * sizeof(char)) strcpy(tmp, s) sig_free(s) s = tmp t = s + len_so_far #endif
def _unpickle(self, data, int version): elif isinstance(data, list): self._unpickle_matrix_2x2_version0(data) else: raise RuntimeError("invalid pickle data") else: raise RuntimeError("unknown matrix version (=%s)"%version)
cdef _unpickle_version0(self, data): cdef Py_ssize_t i, j, n, k raise RuntimeError("invalid pickle data") raise RuntimeError("invalid pickle data")
def _unpickle_matrix_2x2_version0(self, data): if len(data) != 4 or self._nrows != 2 or self._ncols != 2: raise RuntimeError("invalid pickle data") self.set_unsafe(0, 0, data[0]) self.set_unsafe(0, 1, data[1]) self.set_unsafe(1, 0, data[2]) self.set_unsafe(1, 1, data[3])
######################################################################## # LEVEL 1 helpers: # These function support the implementation of the level 1 functionality. ######################################################################## cdef Matrix_integer_dense _new(self, Py_ssize_t nrows, Py_ssize_t ncols): """ Return a new matrix over the integers from given parent All memory is allocated for this matrix, but its entries have not yet been filled in. """ else:
######################################################################## # LEVEL 2 functionality # x * cdef _add_ # x * cdef _sub_ # x * cdef _mul_ # x * cpdef _cmp_ # x * __neg__ # x * __invert__ -> SEE LEVEL 3 FUNCTIONALITIES # 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) ########################################################################
# cdef _mul_(self, Matrix right): # def _multiply_classical(left, matrix.Matrix _right): # def _list(self): # def _dict(self):
def __copy__(self): r""" Returns a new copy of this matrix.
EXAMPLES::
sage: a = matrix(ZZ,1,3, [1,2,-3]); a [ 1 2 -3] sage: b = a.__copy__(); b [ 1 2 -3] sage: b is a False sage: b == a True
sage: M = MatrixSpace(ZZ,2,3) sage: m = M([1,2,3,3,2,1]) sage: mc = m.__copy__() sage: mc == m and mc is not m True """ cdef Matrix_integer_dense A
def __nonzero__(self): r""" Tests whether self is not the zero matrix.
EXAMPLES::
sage: a = MatrixSpace(ZZ, 2, 3)(range(6)); a [0 1 2] [3 4 5] sage: a.__nonzero__() True sage: (a - a).__nonzero__() False
::
sage: a = MatrixSpace(ZZ, 0, 3)() sage: a.__nonzero__() False sage: a = MatrixSpace(ZZ, 3, 0)() sage: a.__nonzero__() False sage: a = MatrixSpace(ZZ, 0, 0)() sage: a.__nonzero__() False """
def is_one(self): r""" Tests whether self is the identity matrix.
EXAMPLES::
sage: matrix(2, [1,0,0,1]).is_one() True sage: matrix(2, [1,1,0,1]).is_one() False sage: matrix(2, 3, [1,0,0,0,1,0]).is_one() False """
def _multiply_linbox(self, Matrix_integer_dense right): """ Multiply matrices over ZZ using linbox.
.. warning::
This is very slow right now, i.e., linbox is very slow.
EXAMPLES::
sage: A = matrix(ZZ, 2, 3, range(6)) sage: A * A.transpose() [ 5 14] [14 50] sage: A._multiply_linbox(A.transpose()) [ 5 14] [14 50]
TESTS:
This fixes a bug found in :trac:`17094`::
sage: A = identity_matrix(ZZ, 3) sage: A._multiply_linbox(A) [1 0 0] [0 1 0] [0 0 1] """ cdef Matrix_integer_dense ans
def _multiply_classical(self, Matrix_integer_dense right): """ EXAMPLES::
sage: n = 3 sage: a = MatrixSpace(ZZ,n,n)(range(n^2)) sage: b = MatrixSpace(ZZ,n,n)(range(1, n^2 + 1)) sage: a._multiply_classical(b) [ 18 21 24] [ 54 66 78] [ 90 111 132]
TESTS::
sage: for _ in range(100): ....: nrows = randint(0, 10) ....: nmid = randint(0, 10) ....: ncols = randint(0, 10) ....: m1 = random_matrix(ZZ, nrows, nmid) ....: m2 = random_matrix(ZZ, nmid, ncols) ....: ans_flint = m1 * m2 ....: ans_classical = m1._multiply_classical(m2) ....: ans_linbox = m1._multiply_linbox(m2) ....: if ans_flint != ans_classical: ....: raise RuntimeError("ERROR\nm1=\n{}\nm2=\n{}\nans_flint=\n{}\nans_classical=\n{}".format( ....: m1.str(), m2.str(), ans_flint.str(), ans_classical.str())) ....: if ans_flint != ans_linbox: ....: raise RuntimeError("ERROR\nm1=\n{}\nm2=\n{}\nans_flint=\n{}\nans_linbox=\n{}".format( ....: m1.str(), m2.str(), ans_flint.str(), ans_linbox.str())) """ raise IndexError("Number of columns of self must equal number of rows of right.")
cdef Py_ssize_t i, j, k, nr, nc, snc cdef object parent
# self acts on the space of right # right acts on the space of self else:
cdef Matrix_integer_dense M, _right
cdef fmpz_t s
cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right): cdef Matrix_integer_dense M
raise IndexError("Number of columns of self must equal number of rows of right.")
cpdef _lmul_(self, Element right): """ EXAMPLES::
sage: a = matrix(ZZ, 2, range(6)) sage: 5 * a [ 0 5 10] [15 20 25] """ cdef fmpz_t z
cpdef _add_(self, right): """ Add two dense matrices over ZZ.
EXAMPLES::
sage: a = MatrixSpace(ZZ,3)(range(9)) sage: a+a [ 0 2 4] [ 6 8 10] [12 14 16] sage: b = MatrixSpace(ZZ,3)(range(9)) sage: b.swap_rows(1,2) sage: a+b [ 0 2 4] [ 9 11 13] [ 9 11 13] """
cpdef _sub_(self, right): """ Subtract two dense matrices over ZZ.
EXAMPLES::
sage: M = Mat(ZZ,3) sage: a = M(range(9)); b = M(reversed(range(9))) sage: a - b [-8 -6 -4] [-2 0 2] [ 4 6 8] """
def __pow__(sself, n, dummy): r""" Return the ``n``-th power of this matrix.
EXAMPLES::
sage: M = MatrixSpace(ZZ,3) sage: m = M([1, 1, 1, 2, 1, 1, -3, -2, -1]) sage: m ** 3 [-3 -2 -1] [-3 -2 0] [ 2 1 -3] sage: m ** -2 [ 2 -3 -1] [-4 4 1] [ 1 0 0] sage: M(range(9)) ** -1 Traceback (most recent call last): ... ZeroDivisionError: Matrix is singular
TESTS::
sage: m ** 3 == m ** 3r == (~m) ** (-3) == (~m) ** (-3r) True
The following exponents do not fit in an unsigned long and the multiplication method fall back to the generic power implementation in :mod:`sage.structure.element`::
sage: m = M.identity_matrix() sage: m ** (2**256) [1 0 0] [0 1 0] [0 0 1] sage: m ** (2r**256r) [1 0 0] [0 1 0] [0 0 1]
In this case, the second argument to ``__pow__`` is a matrix, which should raise the correct error::
sage: M = Matrix(2, 2, range(4)) sage: None^M Traceback (most recent call last): ... TypeError: Cannot convert NoneType to sage.matrix.matrix_integer_dense.Matrix_integer_dense sage: M^M Traceback (most recent call last): ... NotImplementedError: the given exponent is not supported """
raise ValueError raise ArithmeticError("self must be a square matrix")
cdef unsigned long e
else: else:
else: # it is very likely that the following will never finish except # if self has only eigenvalues 0, 1 or -1.
def __neg__(self): r""" Return the negative of this matrix.
TESTS::
sage: a = matrix(ZZ,2,range(4)) sage: a.__neg__() [ 0 -1] [-2 -3] sage: -a [ 0 -1] [-2 -3] """
cpdef int _cmp_(self, right) except -2: r""" Compare ``self`` with ``right``, examining entries in lexicographic (row major) ordering.
EXAMPLES::
sage: Matrix(ZZ, [[0, 10], [20, 30]]) == (Matrix(ZZ, [[0, 10], [20, 30]])) True sage: Matrix(ZZ, [[0, 10], [20, 30]]) < (Matrix(ZZ, [[0, 15], [20, 30]])) True sage: Matrix(ZZ, [[5, 10], [20, 30]]) > (Matrix(ZZ, [[0, 15], [20, 30]])) True sage: Matrix(ZZ, [[5, 10], [20, 30]]) <= (Matrix(ZZ, [[0, 10], [25, 30]])) False """ cdef Py_ssize_t i, j cdef int k
sig_on() for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: k = fmpz_cmp(fmpz_mat_entry(self._matrix,i,j),fmpz_mat_entry((<Matrix_integer_dense>right)._matrix,i,j)) if k: sig_off() if k < 0: return -1 else: return 1 sig_off() return 0
# TODO: Implement better 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(ZZ,2, [1,2,3,4]) sage: V = ZZ^2 sage: w = V([-1,5]) sage: w*B (14, 18) """ cdef Vector_integer_dense w, ans cdef Py_ssize_t i, j cdef fmpz_t x cdef fmpz_t z
######################################################################## # LEVEL 3 functionality (Optional) # * __deepcopy__ # x * __invert__ # * Matrix windows -- only if you need strassen for that base # * Other functions (list them here): # * Specialized echelon form ########################################################################
def _clear_denom(self): """ INPUT:
- ``self`` - a matrix
OUTPUT: self, 1
EXAMPLES::
sage: a = matrix(ZZ,2,[1,2,3,4]) sage: a._clear_denom() ( [1 2] [3 4], 1 ) """
def charpoly(self, var='x', algorithm=None): """ .. NOTE::
The characteristic polynomial is defined as `\det(xI-A)`.
INPUT:
- ``var`` - a variable name
- ``algorithm`` - (optional) either 'generic', 'flint' or 'linbox'. Default is set to 'linbox'.
EXAMPLES::
sage: A = matrix(ZZ,6, range(36)) sage: f = A.charpoly(); f x^6 - 105*x^5 - 630*x^4 sage: f(A) == 0 True sage: g = A.charpoly(algorithm='flint') sage: f == g True sage: n=20; A = Mat(ZZ,n)(range(n^2)) sage: A.charpoly() x^20 - 3990*x^19 - 266000*x^18 sage: A.minpoly() x^3 - 3990*x^2 - 266000*x
On non square matrices, this method raises an ArithmeticError::
sage: matrix(ZZ, 2, 1).charpoly() Traceback (most recent call last): ... ArithmeticError: only valid for square matrix
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(ZZ, 2) sage: A = M(range(0, 2^2)) sage: type(A) <type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'> sage: A.charpoly('x') x^2 - 3*x - 2 sage: A.charpoly('y') y^2 - 3*y - 2 sage: A._cache['charpoly_linbox'] x^2 - 3*x - 2
Test corner cases::
sage: matrix([5]).charpoly('z', 'flint') z - 5 sage: matrix([5]).charpoly('z', 'linbox') z - 5 sage: matrix([5]).charpoly('z', 'generic') z - 5
sage: matrix([]).charpoly('y', 'flint') 1 sage: matrix([]).charpoly('y', 'linbox') 1 sage: matrix([]).charpoly('y', 'generic') 1
sage: matrix([0]).charpoly('x', 'flint') x sage: matrix([0]).charpoly('x', 'linbox') x sage: matrix([0]).charpoly('x', 'generic') x
Consistency on random inputs::
sage: for _ in range(100): ....: dim = randint(1, 20) ....: m = random_matrix(ZZ, dim) ....: m._clear_cache(); ans_flint = m.charpoly(algorithm='flint') ....: m._clear_cache(); ans_linbox = m.charpoly(algorithm='linbox') ....: m._clear_cache(); ans_generic = m.charpoly(algorithm='generic') ....: if ans_flint != ans_linbox or ans_flint != ans_generic: ....: raise RuntimeError("ans_flint = {}, ans_linbox = {} and ans_generic = {} for\n{}".format( ....: ans_flint, ans_linbox, ans_generic, m.str())) """
cdef Polynomial_integer_dense_flint g
else: raise ValueError("no algorithm '%s'"%algorithm)
def minpoly(self, var='x', algorithm=None): """ INPUT:
- ``var`` - a variable name
- ``algorithm`` - (optional) either 'linbox' (default) or 'generic'
EXAMPLES::
sage: A = matrix(ZZ, 6, range(36)) sage: A.minpoly() x^3 - 105*x^2 - 630*x
sage: A = Mat(ZZ, 6)([k^2 for k in range(36)]) sage: A.minpoly(algorithm='linbox') x^4 - 2695*x^3 - 257964*x^2 + 1693440*x sage: A.minpoly(algorithm='generic') x^4 - 2695*x^3 - 257964*x^2 + 1693440*x
On non square matrices, this method raises an ArithmeticError::
sage: matrix(ZZ, 2, 1).minpoly() Traceback (most recent call last): ... ArithmeticError: only valid for square matrix
TESTS:
Corner cases::
sage: matrix([5]).minpoly('z', 'linbox') z - 5 sage: matrix([5]).minpoly('z', 'generic') z - 5
sage: matrix([]).minpoly('y', 'linbox') 1 sage: matrix([]).minpoly('y', 'generic') 1
sage: matrix(ZZ, 2).minpoly('x', 'linbox') x sage: matrix(ZZ, 2).minpoly('x', 'generic') x
Consistency on random inputs::
sage: for _ in range(100): ....: dim = randint(1, 20) ....: m = random_matrix(ZZ, dim) ....: m._clear_cache(); ans_generic = m.minpoly(algorithm='generic') ....: m._clear_cache(); ans_linbox = m.minpoly(algorithm='linbox') ....: if ans_generic != ans_linbox: ....: raise RuntimeError("ans_generic = {} and ans_linbox = {} for\n{}".format( ....: ans_generic, ans_linbox, m.str())) """
cdef Polynomial_integer_dense_flint g
else: raise ValueError("no algorithm '%s'"%algorithm)
def height(self): """ Return the height of this matrix, i.e., the max absolute value of the entries of the matrix.
OUTPUT: A nonnegative integer.
EXAMPLES::
sage: a = Mat(ZZ,3)(range(9)) sage: a.height() 8 sage: a = Mat(ZZ,2,3)([-17,3,-389,15,-1,0]); a [ -17 3 -389] [ 15 -1 0] sage: a.height() 389 """
cdef int mpz_height(self, mpz_t height) except -1: """ Used to compute the height of this matrix.
INPUT:
- ``height`` -- a GMP mpz_t which has been initialized
OUTPUT: sets the value of height to the height of this matrix, i.e., the max absolute value of the entries of the matrix. """ cdef fmpz_t x,h cdef Py_ssize_t i,j
def _multiply_multi_modular(self, Matrix_integer_dense right): """ Multiply this matrix by ``left`` using a multi modular algorithm.
EXAMPLES::
sage: M = Matrix(ZZ, 2, 3, range(5,11)) sage: N = Matrix(ZZ, 3, 2, range(15,21)) sage: M._multiply_multi_modular(N) [310 328] [463 490] sage: M._multiply_multi_modular(-N) [-310 -328] [-463 -490] """ cdef Integer h cdef mod_int *moduli cdef int i, n, k cdef object parent
cdef Matrix_integer_dense result
def _mod_int(self, modulus): """ Reduce the integer matrix modulo a positive integer.
EXAMPLES::
sage: M = Matrix(ZZ, 2, [1,2,-2,3]) sage: M._mod_int(2) [1 0] [0 1] sage: M._mod_int(1000000) [ 1 2] [999998 3] """ raise OverflowError else:
cdef _mod_two(self): cdef Matrix_mod2_dense res
cdef _mod_int_c(self, mod_int p):
cdef Py_ssize_t i, j cdef mpz_t* self_row
cdef float* res_row_f cdef Matrix_modn_dense_float res_f
cdef double* res_row_d cdef Matrix_modn_dense_double res_d
else: raise ValueError("p to big.")
def _reduce(self, moduli):
moduli = MultiModularBasis(moduli)
cdef MultiModularBasis mm
res.append( Matrix_modn_dense_float.__new__(Matrix_modn_dense_float, matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None) ) None, None, None) ) else: raise ValueError("p=%d too big."%p)
cdef size_t i, k, n cdef Py_ssize_t nr, nc cdef mpz_t tmp
cdef mod_int *entry_list raise MemoryError("out of memory allocating multi-modular coefficient list")
(<Matrix_modn_dense_float>res[k])._matrix[i][j] = (<float>entry_list[k])%(<Matrix_modn_dense_float>res[k]).p else:
def _echelon_in_place_classical(self): cdef Matrix_integer_dense E E = self.echelon_form() sig_on() fmpz_mat_set(self._matrix,E._matrix) sig_off() self.clear_cache()
def _echelon_strassen(self): raise NotImplementedError
def _magma_init_(self, magma): """ EXAMPLES::
sage: m = matrix(ZZ,2,3,[1,2,-3,1,-2,-45]) sage: m._magma_init_(magma) 'Matrix(IntegerRing(),2,3,StringToIntegerSequence("1 2 -3 1 -2 -45"))' sage: magma(m) # optional - magma [ 1 2 -3] [ 1 -2 -45] """
def symplectic_form(self): r""" Find a symplectic basis for self if self is an anti-symmetric, alternating matrix.
Returns a pair (F, C) such that the rows of C form a symplectic basis for self and F = C \* self \* C.transpose().
Raises a ValueError if self is not anti-symmetric, or self is not alternating.
Anti-symmetric means that `M = -M^t`. Alternating means that the diagonal of `M` is identically zero.
A symplectic basis is a basis of the form `e_1, \ldots, e_j, f_1, \ldots f_j, z_1, \dots, z_k` such that
- `z_i M v^t` = 0 for all vectors `v`
- `e_i M {e_j}^t = 0` for all `i, j`
- `f_i M {f_j}^t = 0` for all `i, j`
- `e_i M {f_i}^t = 1` for all `i`
- `e_i M {f_j}^t = 0` for all `i` not equal `j`.
The ordering for the factors `d_{i} | d_{i+1}` and for the placement of zeroes was chosen to agree with the output of ``smith_form``.
See the example for a pictorial description of such a basis.
EXAMPLES::
sage: E = matrix(ZZ, 5, 5, [0, 14, 0, -8, -2, -14, 0, -3, -11, 4, 0, 3, 0, 0, 0, 8, 11, 0, 0, 8, 2, -4, 0, -8, 0]); E [ 0 14 0 -8 -2] [-14 0 -3 -11 4] [ 0 3 0 0 0] [ 8 11 0 0 8] [ 2 -4 0 -8 0] sage: F, C = E.symplectic_form() sage: F [ 0 0 1 0 0] [ 0 0 0 2 0] [-1 0 0 0 0] [ 0 -2 0 0 0] [ 0 0 0 0 0] sage: F == C * E * C.transpose() True sage: E.smith_form()[0] [1 0 0 0 0] [0 1 0 0 0] [0 0 2 0 0] [0 0 0 2 0] [0 0 0 0 0] """
hermite_form = echelon_form
def echelon_form(self, algorithm="default", proof=None, include_zero_rows=True, transformation=False, D=None): r""" Return the echelon form of this matrix over the integers, also known as the hermite normal form (HNF).
INPUT:
- ``algorithm`` -- String. The algorithm to use. Valid options are:
- ``'default'`` -- Let Sage pick an algorithm (default). Up to 75 rows or columns with no transformation matrix, use pari with flag 0; otherwise, use flint.
- ``'flint'`` - use flint
- ``'ntl'`` - use NTL (only works for square matrices of full rank!)
- ``'padic'`` - an asymptotically fast p-adic modular algorithm, If your matrix has large coefficients and is small, you may also want to try this.
- ``'pari'`` - use PARI with flag 1
- ``'pari0'`` - use PARI with flag 0
- ``'pari1'`` - use PARI with flag 1
- ``'pari4'`` - use PARI with flag 4 (use heuristic LLL)
- ``proof`` - (default: True); if proof=False certain determinants are computed using a randomized hybrid p-adic multimodular strategy until it stabilizes twice (instead of up to the Hadamard bound). It is *incredibly* unlikely that one would ever get an incorrect result with proof=False.
- ``include_zero_rows`` - (default: True) if False, don't include zero rows
- ``transformation`` - if given, also compute transformation matrix; only valid for flint and padic algorithm
- ``D`` - (default: None) if given and the algorithm is 'ntl', then D must be a multiple of the determinant and this function will use that fact.
OUTPUT:
The Hermite normal form (=echelon form over `\ZZ`) of self as an immutable matrix.
EXAMPLES::
sage: A = MatrixSpace(ZZ,2)([1,2,3,4]) sage: A.echelon_form() [1 0] [0 2] sage: A = MatrixSpace(ZZ,5)(range(25)) sage: A.echelon_form() [ 5 0 -5 -10 -15] [ 0 1 2 3 4] [ 0 0 0 0 0] [ 0 0 0 0 0] [ 0 0 0 0 0]
Getting a transformation matrix in the nonsquare case::
sage: A = matrix(ZZ,5,3,[1..15]) sage: H, U = A.hermite_form(transformation=True, include_zero_rows=False) sage: H [1 2 3] [0 3 6] sage: U [ 0 0 0 4 -3] [ 0 0 0 13 -10] sage: U*A == H True
TESTS:
Make sure the zero matrices are handled correctly::
sage: m = matrix(ZZ,3,3,[0]*9) sage: m.echelon_form() [0 0 0] [0 0 0] [0 0 0] sage: m = matrix(ZZ,3,1,[0]*3) sage: m.echelon_form() [0] [0] [0] sage: m = matrix(ZZ,1,3,[0]*3) sage: m.echelon_form() [0 0 0]
The ultimate border case!
::
sage: m = matrix(ZZ,0,0,[]) sage: m.echelon_form() []
.. NOTE::
If 'ntl' is chosen for a non square matrix this function raises a ValueError.
Special cases: 0 or 1 rows::
sage: a = matrix(ZZ, 1,2,[0,-1]) sage: a.hermite_form() [0 1] sage: a.pivots() (1,) sage: a = matrix(ZZ, 1,2,[0,0]) sage: a.hermite_form() [0 0] sage: a.pivots() () sage: a = matrix(ZZ,1,3); a [0 0 0] sage: a.echelon_form(include_zero_rows=False) [] sage: a.echelon_form(include_zero_rows=True) [0 0 0]
Illustrate using various algorithms.::
sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='pari') [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='pari0') [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='pari4') [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='padic') [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='default') [1 2 3] [0 3 6] [0 0 0]
The 'ntl' algorithm doesn't work on matrices that do not have full rank.::
sage: matrix(ZZ,3,[1..9]).hermite_form(algorithm='ntl') Traceback (most recent call last): ... ValueError: ntl only computes HNF for square matrices of full rank. sage: matrix(ZZ,3,[0] +[2..9]).hermite_form(algorithm='ntl') [1 0 0] [0 1 0] [0 0 3]
TESTS:
This example illustrated :trac:`2398`::
sage: a = matrix([(0, 0, 3), (0, -2, 2), (0, 1, 2), (0, -2, 5)]) sage: a.hermite_form() [0 1 2] [0 0 3] [0 0 0] [0 0 0]
Check that :trac:`12280` is fixed::
sage: m = matrix([(-2, 1, 9, 2, -8, 1, -3, -1, -4, -1), ....: (5, -2, 0, 1, 0, 4, -1, 1, -2, 0), ....: (-11, 3, 1, 0, -3, -2, -1, -11, 2, -2), ....: (-1, 1, -1, -2, 1, -1, -1, -1, -1, 7), ....: (-2, -1, -1, 1, 1, -2, 1, 0, 2, -4)]).stack( ....: 200 * identity_matrix(ZZ, 10)) sage: matrix(ZZ,m).hermite_form(algorithm='pari', include_zero_rows=False) [ 1 0 2 0 13 5 1 166 72 69] [ 0 1 1 0 20 4 15 195 65 190] [ 0 0 4 0 24 5 23 22 51 123] [ 0 0 0 1 23 7 20 105 60 151] [ 0 0 0 0 40 4 0 80 36 68] [ 0 0 0 0 0 10 0 100 190 170] [ 0 0 0 0 0 0 25 0 100 150] [ 0 0 0 0 0 0 0 200 0 0] [ 0 0 0 0 0 0 0 0 200 0] [ 0 0 0 0 0 0 0 0 0 200] sage: matrix(ZZ,m).hermite_form(algorithm='padic', include_zero_rows=False) [ 1 0 2 0 13 5 1 166 72 69] [ 0 1 1 0 20 4 15 195 65 190] [ 0 0 4 0 24 5 23 22 51 123] [ 0 0 0 1 23 7 20 105 60 151] [ 0 0 0 0 40 4 0 80 36 68] [ 0 0 0 0 0 10 0 100 190 170] [ 0 0 0 0 0 0 25 0 100 150] [ 0 0 0 0 0 0 0 200 0 0] [ 0 0 0 0 0 0 0 0 200 0] [ 0 0 0 0 0 0 0 0 0 200]
Check that the output is correct in corner cases, see :trac:`18613`::
sage: m = matrix(2, 0) sage: m.parent() Full MatrixSpace of 2 by 0 dense matrices over Integer Ring sage: H, U = m.echelon_form(transformation=True) sage: H.parent() Full MatrixSpace of 2 by 0 dense matrices over Integer Ring sage: H.is_immutable() True sage: U [1 0] [0 1] sage: H == U * m True sage: H, U = m.echelon_form(transformation=True, ....: include_zero_rows=False) sage: H.parent() Full MatrixSpace of 0 by 0 dense matrices over Integer Ring sage: U.parent() Full MatrixSpace of 0 by 2 dense matrices over Integer Ring sage: H == U * m True sage: m = random_matrix(ZZ, 100, 100, x=-1000, y=1000, density=.1) sage: m.parent() Full MatrixSpace of 100 by 100 dense matrices over Integer Ring sage: H, U = m.hermite_form(algorithm="flint", transformation=True) sage: H == U*m True """
cdef Matrix_integer_dense H_m,w,U cdef Py_ssize_t nr, nc, n, i, j else:
else:
else: H_m, U = matrix_integer_dense_hnf.hnf_with_transformation(self, proof=proof) if not include_zero_rows: r = H_m.rank() H_m = H_m[:r] U = U[:r] else: raise ValueError("transformation matrix only available with p-adic algorithm") else: raise ValueError("ntl only computes HNF for square matrices of full rank.")
else: H_m = self.new_matrix(nrows=w1.nrows())
else: raise ValueError("algorithm %r not understood" % algorithm)
else:
def saturation(self, p=0, proof=None, max_dets=5): r""" Return a saturation matrix of self, which is a matrix whose rows span the saturation of the row span of self. This is not unique.
The saturation of a `\ZZ` module `M` embedded in `\ZZ^n` is the a module `S` that contains `M` with finite index such that `\ZZ^n/S` is torsion free. This function takes the row span `M` of self, and finds another matrix of full rank with row span the saturation of `M`.
INPUT:
- ``p`` - (default: 0); if nonzero given, saturate only at the prime `p`, i.e., return a matrix whose row span is a `\ZZ`-module `S` that contains self and such that the index of `S` in its saturation is coprime to `p`. If `p` is None, return full saturation of self.
- ``proof`` - (default: use proof.linear_algebra()); if False, the determinant calculations are done with proof=False.
- ``max_dets`` - (default: 5); technical parameter - max number of determinant to compute when bounding prime divisor of self in its saturation.
OUTPUT:
- ``matrix`` - a matrix over ZZ
.. NOTE::
The result is *not* cached.
ALGORITHM: 1. Replace input by a matrix of full rank got from a subset of the rows. 2. Divide out any common factors from rows. 3. Check max_dets random dets of submatrices to see if their GCD (with p) is 1 - if so matrix is saturated and we're done. 4. Finally, use that if A is a matrix of full rank, then `hnf(transpose(A))^{-1}*A` is a saturation of A.
EXAMPLES::
sage: A = matrix(ZZ, 3, 5, [-51, -1509, -71, -109, -593, -19, -341, 4, 86, 98, 0, -246, -11, 65, 217]) sage: A.echelon_form() [ 1 5 2262 20364 56576] [ 0 6 35653 320873 891313] [ 0 0 42993 386937 1074825] sage: S = A.saturation(); S [ -51 -1509 -71 -109 -593] [ -19 -341 4 86 98] [ 35 994 43 51 347]
Notice that the saturation spans a different module than A.
::
sage: S.echelon_form() [ 1 2 0 8 32] [ 0 3 0 -2 -6] [ 0 0 1 9 25] sage: V = A.row_space(); W = S.row_space() sage: V.is_submodule(W) True sage: V.index_in(W) 85986 sage: V.index_in_saturation() 85986
We illustrate each option::
sage: S = A.saturation(p=2) sage: S = A.saturation(proof=False) sage: S = A.saturation(max_dets=2) """
def index_in_saturation(self, proof=None): """ Return the index of self in its saturation.
INPUT:
- ``proof`` - (default: use proof.linear_algebra()); if False, the determinant calculations are done with proof=False.
OUTPUT:
- ``positive integer`` - the index of the row span of this matrix in its saturation
ALGORITHM: Use Hermite normal form twice to find an invertible matrix whose inverse transforms a matrix with the same row span as self to its saturation, then compute the determinant of that matrix.
EXAMPLES::
sage: A = matrix(ZZ, 2,3, [1..6]); A [1 2 3] [4 5 6] sage: A.index_in_saturation() 3 sage: A.saturation() [1 2 3] [1 1 1] """
def pivots(self): """ Return the pivot column positions of this matrix.
OUTPUT: a tuple of Python integers: the position of the first nonzero entry in each row of the echelon form.
EXAMPLES::
sage: n = 3; A = matrix(ZZ,n,range(n^2)); A [0 1 2] [3 4 5] [6 7 8] sage: A.pivots() (0, 1) sage: A.echelon_form() [ 3 0 -3] [ 0 1 2] [ 0 0 0] """
cdef Matrix_integer_dense E
# Now we determine the pivots from the matrix E as quickly as we can. # For each row, we find the first nonzero position in that row -- it is the pivot. cdef Py_ssize_t i, j, k
#### Elementary divisors
def elementary_divisors(self, algorithm='pari'): """ Return the elementary divisors of self, in order.
.. warning::
This is MUCH faster than the smith_form function.
The elementary divisors are the invariants of the finite abelian group that is the cokernel of *left* multiplication of this matrix. They are ordered in reverse by divisibility.
INPUT:
- ``self`` - matrix
- ``algorithm`` - (default: 'pari')
- ``'pari'``: works robustly, but is slower.
- ``'linbox'`` - use linbox (currently off, broken)
OUTPUT: list of integers
.. NOTE::
These are the invariants of the cokernel of *left* multiplication::
sage: M = Matrix([[3,0,1],[0,1,0]]) sage: M [3 0 1] [0 1 0] sage: M.elementary_divisors() [1, 1] sage: M.transpose().elementary_divisors() [1, 1, 0]
EXAMPLES::
sage: matrix(3, range(9)).elementary_divisors() [1, 3, 0] sage: matrix(3, range(9)).elementary_divisors(algorithm='pari') [1, 3, 0] sage: C = MatrixSpace(ZZ,4)([3,4,5,6,7,3,8,10,14,5,6,7,2,2,10,9]) sage: C.elementary_divisors() [1, 1, 1, 687]
::
sage: M = matrix(ZZ, 3, [1,5,7, 3,6,9, 0,1,2]) sage: M.elementary_divisors() [1, 1, 6]
This returns a copy, which is safe to change::
sage: edivs = M.elementary_divisors() sage: edivs.pop() 6 sage: M.elementary_divisors() [1, 1, 6]
.. SEEALSO::
:meth:`smith_form` """ raise ValueError("linbox too broken -- currently Linbox SNF is disabled.") else: raise ValueError("algorithm (='%s') unknown"%algorithm)
def smith_form(self): r""" Returns matrices S, U, and V such that S = U\*self\*V, and S is in Smith normal form. Thus S is diagonal with diagonal entries the ordered elementary divisors of S.
.. warning::
The elementary_divisors function, which returns the diagonal entries of S, is VASTLY faster than this function.
The elementary divisors are the invariants of the finite abelian group that is the cokernel of this matrix. They are ordered in reverse by divisibility.
EXAMPLES::
sage: A = MatrixSpace(IntegerRing(), 3)(range(9)) sage: D, U, V = A.smith_form() sage: D [1 0 0] [0 3 0] [0 0 0] sage: U [ 0 1 0] [ 0 -1 1] [-1 2 -1] sage: V [-1 4 1] [ 1 -3 -2] [ 0 0 1] sage: U*A*V [1 0 0] [0 3 0] [0 0 0]
It also makes sense for nonsquare matrices::
sage: A = Matrix(ZZ,3,2,range(6)) sage: D, U, V = A.smith_form() sage: D [1 0] [0 2] [0 0] sage: U [ 0 1 0] [ 0 -1 1] [-1 2 -1] sage: V [-1 3] [ 1 -2] sage: U * A * V [1 0] [0 2] [0 0]
Empty matrices are handled sensibly (see :trac:`3068`)::
sage: m = MatrixSpace(ZZ, 2,0)(0); d,u,v = m.smith_form(); u*m*v == d True sage: m = MatrixSpace(ZZ, 0,2)(0); d,u,v = m.smith_form(); u*m*v == d True sage: m = MatrixSpace(ZZ, 0,0)(0); d,u,v = m.smith_form(); u*m*v == d True
.. SEEALSO::
:meth:`elementary_divisors` """ # need to reverse order of rows of U, columns of V, and both of D.
# silly special cases for matrices with 0 columns (PARI has a unique empty matrix) else:
# silly special cases for matrices with 0 rows (PARI has a unique empty matrix) else:
def frobenius(self, flag=0, var='x'): """ Return the Frobenius form (rational canonical form) of this matrix.
INPUT:
- ``flag`` -- 0 (default), 1 or 2 as follows:
- ``0`` -- (default) return the Frobenius form of this matrix.
- ``1`` -- return only the elementary divisor polynomials, as polynomials in var.
- ``2`` -- return a two-components vector [F,B] where F is the Frobenius form and B is the basis change so that `M=B^{-1}FB`.
- ``var`` -- a string (default: 'x')
ALGORITHM: uses PARI's matfrobenius()
EXAMPLES::
sage: A = MatrixSpace(ZZ, 3)(range(9)) sage: A.frobenius(0) [ 0 0 0] [ 1 0 18] [ 0 1 12] sage: A.frobenius(1) [x^3 - 12*x^2 - 18*x] sage: A.frobenius(1, var='y') [y^3 - 12*y^2 - 18*y] sage: F, B = A.frobenius(2) sage: A == B^(-1)*F*B True sage: a=matrix([]) sage: a.frobenius(2) ([], []) sage: a.frobenius(0) [] sage: a.frobenius(1) [] sage: B = random_matrix(ZZ,2,3) sage: B.frobenius() Traceback (most recent call last): ... ArithmeticError: frobenius matrix of non-square matrix not defined.
AUTHORS:
- Martin Albrect (2006-04-02)
TODO: - move this to work for more general matrices than just over Z. This will require fixing how PARI polynomials are coerced to Sage polynomials. """
def _right_kernel_matrix(self, **kwds): r""" Returns a pair that includes a matrix of basis vectors for the right kernel of ``self``.
INPUT:
- ``algorithm`` - determines which algorithm to use, options are:
- 'flint' - use the algorithm from the FLINT library - 'pari' - use the ``matkerint()`` function from the PARI library - 'padic' - use the p-adic algorithm from the IML library - 'default' - use a heuristic to decide which of the three above routines is fastest. This is the default value.
- ``proof`` - this is passed to the p-adic IML algorithm. If not specified, the global flag for linear algebra will be used.
OUTPUT:
Returns a pair. First item is the string is either 'computed-flint-int', 'computed-pari-int', 'computed-flint-int', which identifies the nature of the basis vectors.
Second item is a matrix whose rows are a basis for the right kernel, over the integers, as computed by either the FLINT, IML or PARI libraries.
EXAMPLES::
sage: A = matrix(ZZ, [[4, 7, 9, 7, 5, 0], ....: [1, 0, 5, 8, 9, 1], ....: [0, 1, 0, 1, 9, 7], ....: [4, 7, 6, 5, 1, 4]])
sage: result = A._right_kernel_matrix(algorithm='pari') sage: result[0] 'computed-pari-int' sage: X = result[1]; X [ 26 -31 30 -21 -2 10] [ 47 13 -48 14 11 -18] sage: A*X.transpose() == zero_matrix(ZZ, 4, 2) True
sage: result = A._right_kernel_matrix(algorithm='padic') sage: result[0] 'computed-iml-int' sage: X = result[1]; X [-469 214 -30 119 -37 0] [ 370 -165 18 -91 30 -2] sage: A*X.transpose() == zero_matrix(ZZ, 4, 2) True
sage: result = A._right_kernel_matrix(algorithm='default') sage: result[0] 'computed-flint-int' sage: result[1] [ 469 -214 30 -119 37 0] [-370 165 -18 91 -30 2]
sage: result = A._right_kernel_matrix(algorithm='flint') sage: result[0] 'computed-flint-int' sage: result[1] [ 469 -214 30 -119 37 0] [-370 165 -18 91 -30 2]
With the 'default' given as the algorithm, several heuristics are used to determine if FLINT, PARI or IML ('padic') is used. The code has exact details, but roughly speaking, relevant factors are: the absolute size of the matrix, or the relative dimensions, or the magnitude of the entries. ::
sage: A = random_matrix(ZZ, 18, 11) sage: A._right_kernel_matrix(algorithm='default')[0] 'computed-pari-int' sage: A = random_matrix(ZZ, 18, 11, x = 10^200) sage: A._right_kernel_matrix(algorithm='default')[0] 'computed-iml-int' sage: A = random_matrix(ZZ, 60, 60) sage: A._right_kernel_matrix(algorithm='default')[0] 'computed-iml-int' sage: A = random_matrix(ZZ, 60, 55) sage: A._right_kernel_matrix(algorithm='default')[0] 'computed-pari-int'
TESTS:
We test three trivial cases. FLINT is used for small matrices, but we let the heuristic decide that. ::
sage: A = matrix(ZZ, 0, 2) sage: A._right_kernel_matrix()[1] [] sage: A = matrix(ZZ, 2, 0) sage: A._right_kernel_matrix()[1].parent() Full MatrixSpace of 0 by 0 dense matrices over Integer Ring sage: A = zero_matrix(ZZ, 4, 3) sage: A._right_kernel_matrix()[1] [1 0 0] [0 1 0] [0 0 1] """
# The heuristic here could be auto-tuned, stored for # different architecture, etc. What I've done below here # I just got by playing around with examples. This is # *dramatically* better than doing absolutely nothing # (i.e., always choosing 'padic'), but is of course # far from optimal. -- William Stein
# I sometimes favor FLINT over PARI, but this should be better tuned. -- Marc Masdeu # Use FLINT for very small matrices, as long as entries aren't huge. # when entries are huge, padic relatively good. else: # the padic algorithm is much better for bigger # matrices if there are nearly more columns than rows # (that is its forte) else:
else: raise ValueError('unknown algorithm: %s'%algorithm)
# TODO: implement using flint function def _adjoint(self): """ Return the adjoint of this matrix.
Assumes self is a square matrix (checked in adjoint).
EXAMPLES::
sage: m = matrix(ZZ,3,[1..9]) sage: m.adjoint() [ -3 6 -3] [ 6 -12 6] [ -3 6 -3] """
def _ntl_(self): r""" ntl.mat_ZZ representation of self.
EXAMPLES::
sage: a = MatrixSpace(ZZ,200).random_element(x=-2, y=2) # -2 to 2 sage: A = a._ntl_()
.. NOTE::
NTL only knows dense matrices, so if you provide a sparse matrix NTL will allocate memory for every zero entry. """
#################################################################################### # LLL #################################################################################### def LLL_gram(self, flag = 0): """ LLL reduction of the lattice whose gram matrix is ``self``, assuming that ``self`` is positive definite.
.. WARNING::
The algorithm does not check if ``self`` is positive definite.
INPUT:
- ``self`` -- a gram matrix of a positive definite quadratic form
- ``flag`` -- an optional flag passed to ``qflllgram``. According to :pari:`qflllgram`'s documentation the options are:
- ``0`` -- (default), assume that ``self`` has either exact (integral or rational) or real floating point entries. The matrix is rescaled, converted to integers and the behavior is then as in ``flag = 1``.
- ``1`` -- assume that G is integral. Computations involving Gram-Schmidt vectors are approximate, with precision varying as needed.
OUTPUT:
A dense matrix ``U`` over the integers that represents a unimodular transformation matrix such that ``U.T * M * U`` is LLL-reduced.
ALGORITHM:
Calls PARI's :pari:`qflllgram`.
EXAMPLES::
sage: M = Matrix(ZZ, 2, 2, [5,3,3,2]) ; M [5 3] [3 2] sage: U = M.LLL_gram(); U [-1 1] [ 1 -2] sage: U.transpose() * M * U [1 0] [0 1]
The algorithm might work for some semidefinite and indefinite forms::
sage: Matrix(ZZ,2,2,[2,6,6,3]).LLL_gram() [-3 -1] [ 1 0] sage: Matrix(ZZ,2,2,[1,0,0,-1]).LLL_gram() [ 0 -1] [ 1 0]
However, it might fail for others by raising a ``ValueError``::
sage: Matrix(ZZ, 1,1,[0]).LLL_gram() Traceback (most recent call last): ... ValueError: qflllgram did not return a square matrix, perhaps the matrix is not positive definite
sage: Matrix(ZZ, 2, 2, [0,1,1,0]).LLL_gram() Traceback (most recent call last): ... ValueError: qflllgram did not return a square matrix, perhaps the matrix is not positive definite
or by running forever::
sage: Matrix(ZZ, 2, 2, [-5, -1, -1, -5]).LLL_gram() # not tested Traceback (most recent call last): ... RuntimeError: infinite loop while calling qflllgram
""" raise ArithmeticError("self must be a square matrix")
# maybe should be /unimodular/ matrices ? except (RuntimeError, ArithmeticError) as msg: raise ValueError("qflllgram failed, " "perhaps the matrix is not positive definite") "perhaps the matrix is not positive definite"); # Fix last column so that det = +1
def BKZ(self, delta=None, algorithm="fpLLL", fp=None, block_size=10, prune=0, use_givens=False, precision=0, proof=None, **kwds): """ Block Korkin-Zolotarev reduction.
INPUT:
- ``delta`` -- (default: ``0.99``) LLL parameter
- ``algorithm`` -- (default: ``"fpLLL"``) ``"fpLLL"`` or ``"NTL"``
- ``fp`` -- floating point number implementation
- ``None`` -- NTL's exact reduction or fpLLL's wrapper (default)
- ``'fp'`` -- double precision: NTL's FP or fpLLL's double
- ``'ld'`` -- long doubles (fpLLL only)
- ``'qd'`` -- NTL's QP
-``'qd1'`` -- quad doubles: Uses ``quad_float`` precision to compute Gram-Schmidt, but uses double precision in the search phase of the block reduction algorithm. This seems adequate for most purposes, and is faster than ``'qd'``, which uses quad_float precision uniformly throughout (NTL only).
- ``'xd'`` -- extended exponent: NTL's XD or fpLLL's dpe
- ``'rr'`` -- arbitrary precision: NTL'RR or fpLLL's MPFR
- ``block_size`` -- (default: ``10``) Specifies the size of the blocks in the reduction. High values yield shorter vectors, but the running time increases double exponentially with ``block_size``. ``block_size`` should be between 2 and the number of rows of ``self``.
- ``proof`` -- (default: same as ``proof.linear_algebra()``) Insist on full BKZ reduction. If disabled and fplll is called, reduction is much faster but the result is not fully BKZ reduced.
NLT SPECIFIC INPUT:
- ``prune`` -- (default: ``0``) The optional parameter ``prune`` can be set to any positive number to invoke the Volume Heuristic from [SH1995]_. This can significantly reduce the running time, and hence allow much bigger block size, but the quality of the reduction is of course not as good in general. Higher values of ``prune`` mean better quality, and slower running time. When ``prune`` is ``0``, pruning is disabled. Recommended usage: for ``block_size==30``, set ``10 <= prune <=15``.
- ``use_givens`` -- Use Given's orthogonalization. This is a bit slower, but generally much more stable, and is really the preferred orthogonalization strategy. For a nice description of this, see Chapter 5 of [GL1996]_.
fpLLL SPECIFIC INPUT:
- ``precision`` -- (default: ``0`` for automatic choice) bit precision to use if ``fp='rr'`` is set
- ``**kwds`` -- kwds are passed through to fpylll. See `fpylll.fplll.BKZ.Param` for details.
Also, if the verbose level is at least `2`, some output is printed during the computation.
EXAMPLES::
sage: A = Matrix(ZZ,3,3,range(1,10)) sage: A.BKZ() [ 0 0 0] [ 2 1 0] [-1 1 3]
sage: A = Matrix(ZZ,3,3,range(1,10)) sage: A.BKZ(use_givens=True) [ 0 0 0] [ 2 1 0] [-1 1 3]
sage: A = Matrix(ZZ,3,3,range(1,10)) sage: A.BKZ(fp="fp") [ 0 0 0] [ 2 1 0] [-1 1 3]
ALGORITHM:
Calls either NTL or fpLLL. """ elif delta <= 0.25: raise TypeError("delta must be > 0.25") elif delta > 1: raise TypeError("delta must be <= 1")
raise TypeError("prune must be >= 0")
if fp is None: fp = "rr"
if fp == "fp": algorithm = "BKZ_FP" elif fp == "qd": algorithm = "BKZ_QP" elif fp == "qd1": algorithm = "BKZ_QP1" elif fp == "xd": algorithm = "BKZ_XD" elif fp == "rr": algorithm = "BKZ_RR" else: raise TypeError("fp parameter not understood.")
A = self._ntl_()
if algorithm == "BKZ_FP": if not use_givens: r = A.BKZ_FP(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose) else: r = A.G_BKZ_FP(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose)
elif algorithm == "BKZ_QP": if not use_givens: r = A.BKZ_QP(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose) else: r = A.G_BKZ_QP(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose)
elif algorithm == "BKZ_QP1": if not use_givens: r = A.BKZ_QP1(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose) else: r = A.G_BKZ_QP1(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose)
elif algorithm == "BKZ_XD": if not use_givens: r = A.BKZ_XD(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose) else: r = A.G_BKZ_XD(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose)
elif algorithm == "BKZ_RR": if not use_givens: r = A.BKZ_RR(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose) else: r = A.G_BKZ_RR(U=None, delta=delta, BlockSize=block_size, prune=prune, verbose=verbose)
self.cache("rank",ZZ(r)) R = <Matrix_integer_dense>self.new_matrix(entries=map(ZZ,A.list()))
kwds["flags"] = kwds.get("flags", BKZ.DEFAULT) | BKZ.VERBOSE
# enable performance improvements unless # 1. provable results are requested or # 2. the user has specified the relevant parameters already kwds["strategies"] = load_strategies_json(BKZ.DEFAULT_STRATEGY)
kwds["flags"] = kwds.get("flags", BKZ.DEFAULT) | BKZ.AUTO_ABORT kwds["auto_abort"] = True
def LLL(self, delta=None, eta=None, algorithm="fpLLL:wrapper", fp=None, prec=0, early_red=False, use_givens=False, use_siegel=False, **kwds): r""" Return LLL reduced or approximated LLL reduced lattice `R` for this matrix interpreted as a lattice.
A lattice `(b_1, b_2, ..., b_d)` is `(\delta, \eta)`-LLL-reduced if the two following conditions hold:
- For any `i > j`, we have `\lvert \mu_{i,j} \rvert \leq \eta`.
- For any `i < d`, we have `\delta \lvert b_i^* \rvert^2 \leq \lvert b_{i + 1}^* + \mu_{i+1, i} b_i^* \rvert^2`,
where `μ_{i,j} = \langle b_i, b_j^* \rangle / \langle b_j^*, b_j^* \rangle` and `b_i^*` is the `i`-th vector of the Gram-Schmidt orthogonalisation of `(b_1, b_2, ..., b_d)`.
The default reduction parameters are `\delta = 3/4` and `\eta = 0.501`. The parameters `\delta` and `\eta` must satisfy: `0.25 < \delta \leq 1.0` and `0.5 \leq \eta < \sqrt{\delta}`. Polynomial time complexity is only guaranteed for `\delta < 1`. Not every algorithm admits the case `\delta = 1`.
The lattice is returned as a matrix. Also the rank (and the determinant) of ``self`` are cached if those are computed during the reduction. Note that in general this only happens when ``self.rank() == self.ncols()`` and the exact algorithm is used.
INPUT:
- ``delta`` -- (default: ``0.99``) `\delta` parameter as described above
- ``eta`` -- (default: ``0.501``) `\eta` parameter as described above, ignored by NTL
- ``algorithm`` -- string one of the algorithms listed below (default: ``"fpLLL:wrapper"``).
- ``fp`` -- floating point number implementation:
- ``None`` -- NTL's exact reduction or fpLLL's wrapper - ``'fp'`` -- double precision: NTL's FP or fpLLL's double - ``'ld'`` -- long doubles (fpLLL only) - ``'qd'`` -- NTL's QP - ``'xd'`` -- extended exponent: NTL's XD or fpLLL's dpe - ``'rr'`` -- arbitrary precision: NTL's RR or fpLLL's MPFR
- ``prec`` -- (default: auto choose) precision, ignored by NTL
- ``early_red`` -- (default: ``False``) perform early reduction, ignored by NTL
- ``use_givens`` -- (default: ``False``) use Givens orthogonalization only applicable to approximate reductions and NTL; this is more stable but slower
- ``use_siegel`` -- (default: ``False``) use Siegel's condition instead of Lovasz's condition, ignored by NTL
- ``**kwds`` -- kwds are passed through to fpylll. See `fpylll.fplll.LLL.reduction` for details.
Also, if the verbose level is at least `2`, some output is printed during the computation.
AVAILABLE ALGORITHMS:
- ``NTL:LLL`` - NTL's LLL + choice of ``fp``.
- ``fpLLL:heuristic`` - fpLLL's heuristic + choice of ``fp``.
- ``fpLLL:fast`` - fpLLL's fast + choice of ``fp``.
- ``fpLLL:proved`` - fpLLL's proved + choice of ``fp``.
- ``fpLLL:wrapper`` - fpLLL's automatic choice (default).
OUTPUT:
A matrix over the integers.
EXAMPLES::
sage: A = Matrix(ZZ,3,3,range(1,10)) sage: A.LLL() [ 0 0 0] [ 2 1 0] [-1 1 3]
We compute the extended GCD of a list of integers using LLL, this example is from the Magma handbook::
sage: Q = [ 67015143, 248934363018, 109210, 25590011055, 74631449, ....: 10230248, 709487, 68965012139, 972065, 864972271 ] sage: n = len(Q) sage: S = 100 sage: X = Matrix(ZZ, n, n + 1) sage: for i in range(n): ....: X[i, i + 1] = 1 sage: for i in range(n): ....: X[i, 0] = S * Q[i] sage: L = X.LLL() sage: M = L.row(n-1).list()[1:] sage: M [-3, -1, 13, -1, -4, 2, 3, 4, 5, -1] sage: add(Q[i]*M[i] for i in range(n)) -1
The case `\delta = 1` is not always supported::
sage: L = X.LLL(delta=2) Traceback (most recent call last): ... TypeError: delta must be <= 1 sage: L = X.LLL(delta=1) # not tested, will eat lots of ram Traceback (most recent call last): ... RuntimeError: infinite loop in LLL sage: L = X.LLL(delta=1, algorithm='NTL:LLL') sage: L[-1] (-100, -3, -1, 13, -1, -4, 2, 3, 4, 5, -1)
TESTS::
sage: matrix(ZZ, 0, 0).LLL() [] sage: matrix(ZZ, 3, 0).LLL() [] sage: matrix(ZZ, 0, 3).LLL() []
sage: M = matrix(ZZ, [[1,2,3],[31,41,51],[101,201,301]]) sage: A = M.LLL() sage: A [ 0 0 0] [-1 0 1] [ 1 1 1] sage: B = M.LLL(algorithm='NTL:LLL') sage: C = M.LLL(algorithm='NTL:LLL', fp=None) sage: D = M.LLL(algorithm='NTL:LLL', fp='fp') sage: F = M.LLL(algorithm='NTL:LLL', fp='xd') sage: G = M.LLL(algorithm='NTL:LLL', fp='rr') sage: A == B == C == D == F == G True sage: H = M.LLL(algorithm='NTL:LLL', fp='qd') Traceback (most recent call last): ... TypeError: algorithm NTL:LLL_QD not supported
.. NOTE::
See ``ntl.mat_ZZ`` or ``fpylll.fplll.lll`` for details on the used algorithms.
Albeit LLL is a deterministic algorithm, the output for different implementations and on CPUs (32-bit vs. 64-bit) may vary, while still being correct. """
raise TypeError("precision prec must be >= 0")
raise TypeError("delta must be > 1/4") raise TypeError("delta must be <= 1")
else: raise TypeError("delta must be > 0.25")
elif eta < 0.5: raise TypeError("eta must be >= 0.5")
self.cache("det", det) pass r = A.G_LLL_FP(delta, verbose=verb) else: if use_givens: r = A.G_LLL_QP(delta, verbose=verb) else: r = A.LLL_QP(delta, verbose=verb) r = A.G_LLL_XD(delta, verbose=verb) else: r = A.G_LLL_RR(delta, verbose=verb) else: else:
kwds["flags"] = kwds.get("flags", LLL.DEFAULT) | LLL.SIEGEL kwds["flags"] = kwds.get("flags", LLL.DEFAULT) | LLL.EARLY_RED
else: raise TypeError("algorithm %s not supported"%algorithm)
def is_LLL_reduced(self, delta=None, eta=None): r""" Return ``True`` if this lattice is `(\delta, \eta)`-LLL reduced. See ``self.LLL`` for a definition of LLL reduction.
INPUT:
- ``delta`` -- (default: `0.99`) parameter `\delta` as described above
- ``eta`` -- (default: `0.501`) parameter `\eta` as described above
EXAMPLES::
sage: A = random_matrix(ZZ, 10, 10) sage: L = A.LLL() sage: A.is_LLL_reduced() False sage: L.is_LLL_reduced() True """
raise TypeError("delta must be > 1/4") raise TypeError("delta must be <= 1")
raise TypeError("eta must be >= 0.5")
# this is pretty slow #For any $i>j$, we have $|mu_{i, j}| <= \eta$
#For any $i<d$, we have $\delta |b_i^*|^2 <= |b_{i+1}^* + mu_{i+1, i} b_i^* |^2$ return False
def prod_of_row_sums(self, cols): """ Return the product of the sums of the entries in the submatrix of ``self`` with given columns.
INPUT:
- ``cols`` -- a list (or set) of integers representing columns of ``self``
OUTPUT: an integer
EXAMPLES::
sage: a = matrix(ZZ,2,3,[1..6]); a [1 2 3] [4 5 6] sage: a.prod_of_row_sums([0,2]) 40 sage: (1+3)*(4+6) 40 sage: a.prod_of_row_sums(set([0,2])) 40
""" cdef Py_ssize_t c, row cdef fmpz_t s,pr
fmpz_clear(s) fmpz_clear(pr) raise IndexError("matrix column index out of range")
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.
INPUT:
- ``N`` - an integer
OUTPUT:
- ``matrix`` - over QQ or raise a ValueError
EXAMPLES: We create a random 4x4 matrix over ZZ.
::
sage: A = matrix(ZZ, 4, [4, -4, 7, 1, -1, 1, -1, -12, -1, -1, 1, -1, -3, 1, 5, -1])
There isn't a unique rational reconstruction of it::
sage: A.rational_reconstruction(11) Traceback (most recent call last): ... ValueError: rational reconstruction does not exist
We throw in a denominator and reduce the matrix modulo 389 - it does rationally reconstruct.
::
sage: B = (A/3 % 389).change_ring(ZZ) sage: B.rational_reconstruction(389) == A/3 True
TESTS:
Check that :trac:`9345` is fixed::
sage: A = random_matrix(ZZ, 3, 3) sage: A.rational_reconstruction(0) Traceback (most recent call last): ... ZeroDivisionError: The modulus cannot be zero """
def randomize(self, density=1, x=None, y=None, distribution=None, \ nonzero=False): """ Randomize ``density`` proportion of the entries of this matrix, leaving the rest unchanged.
The parameters are the same as the ones for the integer ring's ``random_element`` function.
If ``x`` and ``y`` are given, randomized entries of this matrix have to be between ``x`` and ``y`` and have density 1.
INPUT:
- ``self`` - a mutable matrix over ZZ
- ``density`` - a float between 0 and 1
- ``x, y`` - if not ``None``, these are passed to the ``ZZ.random_element`` function as the upper and lower endpoints in the uniform distribution
- ``distribution`` - would also be passed into ``ZZ.random_element`` if given
- ``nonzero`` - bool (default: ``False``); whether the new entries are guaranteed to be zero
OUTPUT:
- None, the matrix is modified in-place
EXAMPLES::
sage: A = matrix(ZZ, 2,3, [1..6]); A [1 2 3] [4 5 6] sage: A.randomize() sage: A [-8 2 0] [ 0 1 -1] sage: A.randomize(x=-30,y=30) sage: A [ 5 -19 24] [ 24 23 -9] """ return density = float(1)
cdef mpz_t tmp cdef Py_ssize_t i, j, k, nc, num_per_row global state, ZZ
# Original code, before adding the ``nonzero`` option. distribution) else: nc = self._ncols num_per_row = int(density * nc) for i from 0 <= i < self._nrows: for j from 0 <= j < num_per_row: k = rstate.c_random()%nc the_integer_ring._randomize_mpz(tmp, \ x, y, distribution) self.set_unsafe_mpz(i,k,tmp) else: # New code, to implement the ``nonzero`` option. Note that this # code is almost the same as above, the only difference being that # each entry is set until it's non-zero. for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: while fmpz_sgn(fmpz_mat_entry(self._matrix,i,j)) == 0: the_integer_ring._randomize_mpz(tmp, \ x, y, distribution) self.set_unsafe_mpz(i,j,tmp) else: x, y, distribution)
#### Rank
def rank(self, algorithm='modp'): """ Return the rank of this matrix.
INPUT:
- ``algorithm`` -- either ``'modp'`` (default) or ``'flint'`` or ``'linbox'``
OUTPUT:
- a nonnegative integer -- the rank
.. NOTE::
The rank is cached.
ALGORITHM:
If set to ``'modp'``, first check if the matrix has maximum possible rank by working modulo one random prime. If not, call LinBox's rank function.
EXAMPLES::
sage: a = matrix(ZZ,2,3,[1..6]); a [1 2 3] [4 5 6] sage: a.rank() 2 sage: a = matrix(ZZ,3,3,[1..9]); a [1 2 3] [4 5 6] [7 8 9] sage: a.rank() 2
Here is a bigger example - the rank is of course still 2::
sage: a = matrix(ZZ,100,[1..100^2]); a.rank() 2
TESTS::
sage: a.rank(algorithm='funky') Traceback (most recent call last): ... ValueError: algorithm must be one of 'modp', 'flint' or 'linbox' """ "or 'linbox'")
# Can very quickly detect full rank by working modulo p. # Algorithm is 'linbox' or detecting full rank didn't work -- # use LinBox's general algorithm.
def _rank_linbox(self): """ Compute the rank of this matrix using Linbox.
TESTS::
sage: matrix(ZZ, 4, 6, 0)._rank_linbox() 0 sage: matrix(ZZ, 3, 4, range(12))._rank_linbox() 2 sage: matrix(ZZ, 5, 10, [1+i+i^2 for i in range(50)])._rank_linbox() 3 """
def _rank_modp(self, p=46337):
#### Determinant
def determinant(self, algorithm='default', proof=None, stabilize=2): r""" Return the determinant of this matrix.
INPUT:
- ``algorithm``
- ``'default'`` -- use ``flint``
- ``'flint'`` -- let flint do the determinant
- ``'padic'`` - uses a p-adic / multimodular algorithm that relies on code in IML and linbox
- ``'linbox'`` - calls linbox det (you *must* set proof=False to use this!)
- ``'ntl'`` - calls NTL's det function
- ``'pari'`` - uses PARI
- ``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.
- ``stabilize`` - if proof is False, require det to be the same for this many CRT primes in a row. Ignored if proof is True.
ALGORITHM: The p-adic algorithm works by first finding a random vector v, then solving A\*x = v and taking the denominator `d`. This gives a divisor of the determinant. Then we compute `\det(A)/d` using a multimodular algorithm and the Hadamard bound, skipping primes that divide `d`.
EXAMPLES::
sage: A = matrix(ZZ,8,8,[3..66]) sage: A.determinant() 0
::
sage: A = random_matrix(ZZ,20,20) sage: D1 = A.determinant() sage: A._clear_cache() sage: D2 = A.determinant(algorithm='ntl') sage: D1 == D2 True
We have a special-case algorithm for 4 x 4 determinants::
sage: A = matrix(ZZ,4,[1,2,3,4,4,3,2,1,0,5,0,1,9,1,2,3]) sage: A.determinant() 270
Next we try the Linbox det. Note that we must have proof=False.
::
sage: A = matrix(ZZ,5,[1,2,3,4,5,4,6,3,2,1,7,9,7,5,2,1,4,6,7,8,3,2,4,6,7]) sage: A.determinant(algorithm='linbox') Traceback (most recent call last): ... RuntimeError: you must pass the proof=False option to the determinant command to use LinBox's det algorithm sage: A.determinant(algorithm='linbox', proof=False) -21 sage: A._clear_cache() sage: A.determinant() -21
Try the other algorithms on the same example::
sage: A._clear_cache(); A.determinant(algorithm='padic') -21 sage: A._clear_cache(); A.determinant(algorithm='pari') -21 sage: A._clear_cache(); A.determinant(algorithm='ntl') -21 sage: A._clear_cache(); A.determinant(algorithm='padic') -21
A bigger example::
sage: A = random_matrix(ZZ,30) sage: d = A.determinant() sage: A._clear_cache() sage: A.determinant(algorithm='linbox',proof=False) == d True
TESTS:
This shows that we can compute determinants for all sizes up to 80. The check that the determinant of a squared matrix is a square is a sanity check that the result is probably correct::
sage: for s in [1..80]: # long time ....: M = random_matrix(ZZ, s) ....: d = (M*M).determinant() ....: assert d.is_square()
Check consistency::
sage: all(matrix(ZZ, 0).det(algorithm=algo).is_one() for algo in ['flint', 'padic', 'pari', 'ntl']) True sage: for _ in range(100): ....: dim = randint(1, 10) ....: m = random_matrix(ZZ, dim) ....: det_flint = m.__copy__().det(algorithm='flint') ....: det_padic = m.__copy__().det(algorithm='padic') ....: det_pari = m.__copy__().det(algorithm='pari') ....: det_ntl = m.__copy__().det(algorithm='ntl') ....: if type(det_flint) is not Integer: ....: raise RuntimeError("type(det_flint) = {}".format(type(det_flint))) ....: if type(det_padic) is not Integer: ....: raise RuntimeError("type(det_padic) = {}".format(type(det_padic))) ....: if type(det_pari) is not Integer: ....: raise RuntimeError("type(det_pari) = {}".format(type(det_pari))) ....: if type(det_ntl) is not Integer: ....: raise RuntimeError("type(det_ntl) = {}".format(type(det_ntl))) ....: if det_flint != det_padic or det_flint != det_pari or det_flint != det_ntl: ....: raise RuntimeError("ERROR\ndet_flint = {}\ndet_padic={}\ndet_pari={}\ndet_ntl={}\n{}".format( ....: det_flint, det_padic, det_pari, det_ntl, self.str())) """
cdef fmpz_t e
else: raise TypeError("algorithm '%s' not known"%(algorithm))
def _det_linbox(self): """ Compute the determinant of this matrix using Linbox.
TESTS::
sage: matrix(ZZ, 0)._det_linbox() 1 """ raise ArithmeticError("self must be a square matrix")
cdef fmpz_t tmp
def _det_pari(self, int flag=0): """ Determinant of this matrix using Gauss-Bareiss. If (optional) flag is set to 1, use classical Gaussian elimination.
For efficiency purposes, this det is computed entirely on the PARI stack then the PARI stack is cleared. This function is most useful for very small matrices.
EXAMPLES::
sage: matrix(ZZ, 0)._det_pari() 1 sage: matrix(ZZ, 0)._det_pari(1) 1 sage: matrix(ZZ, 3, [1..9])._det_pari() 0 sage: matrix(ZZ, 3, [1..9])._det_pari(1) 0 """ # now convert d to a Sage integer e
def _det_ntl(self): """ Compute the determinant of this matrix using NTL.
EXAMPLES::
sage: matrix(ZZ, 0)._det_ntl() 1 sage: matrix(ZZ, 3, [1..9])._det_ntl() 0 sage: matrix(ZZ, 3, [1,3,6,2,7,8,2,1,0])._det_ntl() -32 """
def _rational_kernel_iml(self): """ Return the rational (left) kernel of this matrix
OUTPUT:
A matrix ``K`` such that ``self * K = 0``, and the number of columns of K equals the nullity of self.
EXAMPLES::
sage: m = matrix(ZZ, 5, 5, [1+i+i^2 for i in range(25)]) sage: m._rational_kernel_iml() [ 1 3] [-3 -8] [ 3 6] [-1 0] [ 0 -1]
sage: V1 = m._rational_kernel_iml().column_space().change_ring(QQ) sage: V2 = m._rational_kernel_flint().column_space().change_ring(QQ) sage: assert V1 == V2 """
cdef long dim cdef unsigned long i,j,k cdef mpz_t *mp_N # Now read the answer as a matrix. cdef Matrix_integer_dense M
def _rational_kernel_flint(self): """ Return the rational (left) kernel of this matrix
OUTPUT:
A matrix ``K`` such that ``self * K = 0``, and the number of columns of K equals the nullity of self.
EXAMPLES::
sage: m = matrix(ZZ, 4, 6, [i^2-2*i for i in range(24)]) sage: m._rational_kernel_flint() [ 1728 5184 10368] [ -5184 -13824 -25920] [ 5184 10368 17280] [ -1728 0 0] [ 0 -1728 0] [ 0 0 -1728]
sage: V1 = m._rational_kernel_iml().column_space().change_ring(QQ) sage: V2 = m._rational_kernel_flint().column_space().change_ring(QQ) sage: assert V1 == V2 """
cdef long dim cdef fmpz_mat_t M0 # Now read the answer as a matrix. cdef size_t i,j
def _invert_iml(self, use_nullspace=False, check_invertible=True): """ Invert this matrix using IML. The output matrix is an integer matrix and a denominator.
INPUT:
- ``self`` - an invertible matrix
- ``use_nullspace`` - (default: False): whether to use nullspace algorithm, which is slower, but doesn't require checking that the matrix is invertible as a precondition.
- ``check_invertible`` - (default: True) whether to check that the matrix is invertible.
OUTPUT: A, d such that A\*self = d
- ``A`` - a matrix over ZZ
- ``d`` - an integer
ALGORITHM: Uses IML's p-adic nullspace function.
EXAMPLES::
sage: a = matrix(ZZ,3,[1,2,5, 3,7,8, 2,2,1]) sage: b, d = a._invert_iml(); b,d ( [ 9 -8 19] [-13 9 -7] [ 8 -2 -1], 23 ) sage: a*b [23 0 0] [ 0 23 0] [ 0 0 23]
AUTHORS:
- William Stein """ raise TypeError("self must be a square matrix.")
A = self.augment(P.identity_matrix()) K = A._rational_kernel_iml() d = -K[self._nrows,0] if K.ncols() != self._ncols or d == 0: raise ZeroDivisionError("input matrix must be nonsingular") B = K[:self._nrows] verbose("finished computing inverse using IML", time) return B, d else:
def _invert_flint(self): """ Invert this matrix using FLINT. The output matrix is an integer matrix and a denominator.
INPUT:
- ``self`` - an invertible matrix
OUTPUT: A, d such that A\*self = d
- ``A`` - a matrix over ZZ
- ``d`` - an integer
EXAMPLES::
sage: a = matrix(ZZ,3,[1,2,5, 3,7,8, 2,2,1]) sage: b, d = a._invert_flint(); b,d ( [ 9 -8 19] [-13 9 -7] [ 8 -2 -1], 23 ) sage: a*b [23 0 0] [ 0 23 0] [ 0 0 23]
AUTHORS:
- William Stein
- Marc Masdeu -- (08/2014) Use FLINT """
cdef Matrix_integer_dense M cdef int res cdef fmpz_t fden else:
def __invert__(self): r""" Return the inverse of self.
EXAMPLES::
sage: M = MatrixSpace(ZZ,3) sage: m = M([1,2,3,3,4,5,1,2,-3]) sage: ~m [-11/6 1 -1/6] [ 7/6 -1/2 1/3] [ 1/6 0 -1/6] sage: ~m * m == m * ~m == M.identity_matrix() True
Note that inverse of determinant one integer matrices do not belong to the same parent::
sage: (~M.identity_matrix()).parent() Full MatrixSpace of 3 by 3 dense matrices over Rational Field
This is consistent with::
sage: (~1).parent() Rational Field
TESTS::
sage: ~M.zero_matrix() Traceback (most recent call last): ... ZeroDivisionError: Matrix is singular """
def _invert_unit(self): r""" If self is a matrix with determinant `1` or `-1` return the inverse of ``self`` as a matrix over `ZZ`.
EXAMPLES::
sage: a = matrix(2, [1,2,1,1]) sage: a^(-1) [-1 2] [ 1 -1] sage: m = a._invert_unit(); m [-1 2] [ 1 -1] sage: m.parent() Full MatrixSpace of 2 by 2 dense matrices over Integer Ring sage: matrix(2, [2,1,0,1])._invert_unit() Traceback (most recent call last): ... ZeroDivisionError: matrix is not invertible over Integer Ring """
def _solve_right_nonsingular_square(self, B, check_rank=True, algorithm = 'iml'): r""" If self is a matrix `A` of full rank, then this function returns a vector or matrix `X` such that `A X = B`. If `B` is a vector then `X` is a vector and if `B` is a matrix, then `X` is a matrix. The base ring of `X` is the integers unless a denominator is needed in which case the base ring is the rational numbers.
.. NOTE::
In Sage one can also write ``A B`` for ``A.solve_right(B)``, i.e., Sage implements the "the MATLAB/Octave backslash operator".
.. NOTE::
This is currently only implemented when A is square.
INPUT:
- ``B`` - a matrix or vector
- ``check_rank`` - bool (default: True); if True verify that in fact the rank is full.
- ``algorithm`` - ``'iml'`` (default) or ``'flint'``
OUTPUT: a matrix or vector over `\QQ`
EXAMPLES::
sage: a = matrix(ZZ, 2, [0, -1, 1, 0]) sage: v = vector(ZZ, [2, 3]) sage: a \ v (3, -2)
Note that the output vector or matrix is always over `\QQ`.
::
sage: parent(a\v) Vector space of dimension 2 over Rational Field
We solve a bigger system where the answer is over the rationals.
::
sage: a = matrix(ZZ, 3, 3, [1,2,3,4, 5, 6, 8, -2, 3]) sage: v = vector(ZZ, [1,2,3]) sage: w = a \ v; w (2/15, -4/15, 7/15) sage: parent(w) Vector space of dimension 3 over Rational Field sage: a * w (1, 2, 3)
We solve a system where the right hand matrix has multiple columns.
::
sage: a = matrix(ZZ, 3, 3, [1,2,3,4, 5, 6, 8, -2, 3]) sage: b = matrix(ZZ, 3, 2, [1,5, 2, -3, 3, 0]) sage: w = a \ b; w [ 2/15 -19/5] [-4/15 -27/5] [ 7/15 98/15] sage: a * w [ 1 5] [ 2 -3] [ 3 0]
TESTS:
We create a random 100x100 matrix and solve the corresponding system, then verify that the result is correct. (Note that this test is very risky without having a seeded random number generator!)
::
sage: n = 100 sage: a = random_matrix(ZZ,n) sage: v = vector(ZZ,n,range(n)) sage: x = a \ v sage: a * x == v True
sage: n = 100 sage: a = random_matrix(ZZ,n) sage: v = vector(ZZ,n,range(n)) sage: x = a._solve_right_nonsingular_square(v,algorithm = 'flint') sage: a * x == v True
"""
# It would probably be much better to rewrite linbox so it # throws an error instead of ** going into an infinite loop ** # in the non-full rank case. In any case, we do this for now, # since rank is very fast and infinite loops are evil. raise ValueError("self must be of full rank.")
raise NotImplementedError("the input matrix must be square.")
raise ValueError("number of rows of self must equal degree of B.") raise ValueError("number of rows of self must equal number of rows of B.")
else: raise NotImplementedError
# likely much better to just invert then multiply X = self**(-1)*C verbose('finished solve_right (via inverse)', t) return X
else: raise ValueError("Unknown algorithm '%s'"%algorithm) # Convert back to a vector
def _solve_iml(self, Matrix_integer_dense B, right=True): """ Let A equal self be a square matrix. Given B return an integer matrix C and an integer d such that self C\*A == d\*B if right is False or A\*C == d\*B if right is True.
OUTPUT:
- ``C`` - integer matrix
- ``d`` - integer denominator
EXAMPLES::
sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) sage: B = matrix(ZZ,3,4, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) sage: C,d = A._solve_iml(B,right=False); C [ 6 -18 -15 27] [ 0 24 24 -36] [ 4 -12 -6 -2]
::
sage: d 12
::
sage: C*A == d*B True
::
sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) sage: B = matrix(ZZ,4,3, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) sage: C,d = A._solve_iml(B) sage: C [ 12 40 28] [-12 -4 -4] [ -6 -25 -16] [ 12 34 16]
::
sage: d 12
::
sage: A*C == d*B True
Test wrong dimensions::
sage: A = random_matrix(ZZ, 4, 4) sage: B = random_matrix(ZZ, 2, 3) sage: B._solve_iml(A) Traceback (most recent call last): ... ValueError: self must be a square matrix sage: A._solve_iml(B, right=False) Traceback (most recent call last): ... ArithmeticError: B's number of columns must match self's number of rows sage: A._solve_iml(B, right=True) Traceback (most recent call last): ... ArithmeticError: B's number of rows must match self's number of columns
Check that this can be interrupted properly (:trac:`15453`)::
sage: A = random_matrix(ZZ, 2000, 2000) sage: B = random_matrix(ZZ, 2000, 2000) sage: t0 = walltime() sage: alarm(2); A._solve_iml(B) # long time Traceback (most recent call last): ... AlarmInterrupt sage: t = walltime(t0) sage: t < 10 or t True
ALGORITHM: Uses IML.
AUTHORS:
- Martin Albrecht """ cdef unsigned long i, j, k cdef mpz_t *mp_N cdef mpz_t mp_D cdef Matrix_integer_dense M cdef Integer D
# This is *required* by the IML function we call below.
cdef SOLU_POS solu_pos
return P.zero_matrix(), Integer(1)
else: # left
return P.zero_matrix(), Integer(1)
return self.new_matrix(nrows = m, ncols = n), Integer(1)
finally:
def _solve_flint(self, Matrix_integer_dense B, right=True): """ Let A equal self be a square matrix. Given B return an integer matrix C and an integer d such that self C\*A == d\*B if right is False or A\*C == d\*B if right is True.
OUTPUT:
- ``C`` - integer matrix
- ``d`` - integer denominator
EXAMPLES::
sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) sage: B = matrix(ZZ,3,4, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) sage: C,d = A._solve_flint(B,right=False); C [ 6 -18 -15 27] [ 0 24 24 -36] [ 4 -12 -6 -2]
::
sage: d 12
::
sage: C*A == d*B True
::
sage: A = matrix(ZZ,4,4,[0, 1, -2, -1, -1, 1, 0, 2, 2, 2, 2, -1, 0, 2, 2, 1]) sage: B = matrix(ZZ,4,3, [-1, 1, 1, 0, 2, 0, -2, -1, 0, -2, -2, -2]) sage: C,d = A._solve_flint(B) sage: C [ 12 40 28] [-12 -4 -4] [ -6 -25 -16] [ 12 34 16]
::
sage: d 12
::
sage: A*C == d*B True
Test wrong dimensions::
sage: A = random_matrix(ZZ, 4, 4) sage: B = random_matrix(ZZ, 2, 3) sage: B._solve_flint(A) Traceback (most recent call last): ... ValueError: self must be a square matrix sage: A._solve_flint(B, right=False) Traceback (most recent call last): ... ArithmeticError: B's number of columns must match self's number of rows sage: A._solve_flint(B, right=True) Traceback (most recent call last): ... ArithmeticError: B's number of rows must match self's number of columns
Check that this can be interrupted properly (:trac:`15453`)::
sage: A = random_matrix(ZZ, 2000, 2000) sage: B = random_matrix(ZZ, 2000, 2000) sage: t0 = walltime() sage: alarm(2); A._solve_flint(B) # long time Traceback (most recent call last): ... AlarmInterrupt sage: t = walltime(t0) sage: t < 10 or t True
AUTHORS:
- Marc Masdeu (08/2014) following _solve_iml implementation of Martin Albrecht """ cdef Matrix_integer_dense M cdef fmpz_t tmp # This is *required* by the FLINT function we call below.
return B, self[0,0]
return self.new_matrix(nrows = n, ncols = m), Integer(1) mpz_neg(den.value,den.value) fmpz_mat_neg(M._matrix,M._matrix)
else: # left
return self.new_matrix(nrows = n, ncols = m), Integer(1)
def _rational_echelon_via_solve(self, solver = 'iml'): r""" Computes information that gives the reduced row echelon form (over QQ!) of a matrix with integer entries.
INPUT:
- ``self`` - a matrix over the integers.
- ``solver`` - either ``'iml'`` (default) or ``'flint'``
OUTPUT:
- ``pivots`` - ordered list of integers that give the pivot column positions
- ``nonpivots`` - ordered list of the nonpivot column positions
- ``X`` - matrix with integer entries
- ``d`` - integer
If you put standard basis vectors in order at the pivot columns, and put the matrix (1/d)\*X everywhere else, then you get the reduced row echelon form of self, without zero rows at the bottom.
.. NOTE::
IML is the actual underlying `p`-adic solver that we use.
AUTHORS:
- William Stein
ALGORITHM: I came up with this algorithm from scratch. As far as I know it is new. It's pretty simple, but it is ... (fast?!).
Let A be the input matrix.
#. Compute r = rank(A).
#. Compute the pivot columns of the transpose `A^t` of `A`. One way to do this is to choose a random prime `p` and compute the row echelon form of `A^t` modulo `p` (if the number of pivot columns is not equal to `r`, pick another prime).
#. Let `B` be the submatrix of `A` consisting of the rows corresponding to the pivot columns found in the previous step. Note that, aside from zero rows at the bottom, `B` and `A` have the same reduced row echelon form.
#. Compute the pivot columns of `B`, again possibly by choosing a random prime `p` as in [2] and computing the Echelon form modulo `p`. If the number of pivot columns is not `r`, repeat with a different prime. Note - in this step it is possible that we mistakenly choose a bad prime `p` such that there are the right number of pivot columns modulo `p`, but they are at the wrong positions - e.g., imagine the augmented matrix `[pI|I]` - modulo `p` we would miss all the pivot columns. This is OK, since in the next step we would detect this, as the matrix we would obtain would not be in echelon form.
#. Let `C` be the submatrix of `B` of pivot columns. Let `D` be the complementary submatrix of `B` of all all non-pivot columns. Use a `p`-adic solver to find the matrix `X` and integer `d` such that `C (1/d) X=D`. I.e., solve a bunch of linear systems of the form `Cx = v`, where the columns of `X` are the solutions.
#. Verify that we had chosen the correct pivot columns. Inspect the matrix `X` obtained in step 5. If when used to construct the echelon form of `B`, `X` indeed gives a matrix in reduced row echelon form, then we are done - output the pivot columns, `X`, and `d`. To tell if `X` is correct, do the following: For each column of `X` (corresponding to non-pivot column `i` of `B`), read up the column of `X` until finding the first nonzero position `j`; then verify that `j` is strictly less than the number of pivot columns of `B` that are strictly less than `i`. Otherwise, we got the pivots of `B` wrong - try again starting at step 4, but with a different random prime. """
# Step 1: Compute the rank
# The input matrix already has full rank. else: # Steps 2 and 3: Extract out a submatrix of full rank. else: i += 1 if i == 50: raise RuntimeError("Bug in _rational_echelon_via_solve in finding linearly independent rows.")
raise RuntimeError("Bug in _rational_echelon_via_solve -- pivot columns keep being wrong.")
# Step 4: Now we instead worry about computing the reduced row echelon form of B. else: i += 1 if i == 50: raise RuntimeError("Bug in _rational_echelon_via_solve in finding pivot columns.")
# Step 5: Apply p-adic solver else: X, d = C._solve_flint(D, right=True)
# Step 6: Verify that we had chosen the correct pivot columns. break else: pivots_are_right = False break
#end while
# Finally, return the answer.
def decomposition(self, **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, and are saturated as ZZ modules.
INPUT:
- ``self`` - a matrix over the integers
- ``**kwds`` - these are passed onto to the decomposition over QQ command.
EXAMPLES::
sage: t = ModularSymbols(11,sign=1).hecke_matrix(2) sage: w = t.change_ring(ZZ) sage: w.list() [3, -1, 0, -2] """
D, E = X D = [(W.intersection(V), t) for W, t in D] E = [(W.intersection(V), t) for W, t in E] return decomp_seq(D), decomp_seq(E) else:
def _add_row_and_maintain_echelon_form(self, row, pivots): """ Assuming self is a full rank n x m matrix in reduced row Echelon form over ZZ and row is a vector of degree m, this function creates a new matrix that is the echelon form of self with row appended to the bottom.
.. warning::
It is assumed that self is in echelon form.
INPUT:
- ``row`` - a vector of degree m over ZZ
- ``pivots`` - a list of integers that are the pivot columns of self.
OUTPUT:
- ``matrix`` - a matrix of in reduced row echelon form over ZZ
- ``pivots`` - list of integers
ALGORITHM: For each pivot column of self, we use the extended Euclidean algorithm to clear the column. The result is a new matrix B whose row span is the same as self.stack(row), and whose last row is 0 if and only if row is in the QQ-span of the rows of self. If row is not in the QQ-span of the rows of self, then row is nonzero and suitable to be inserted into the top n rows of A to form a new matrix that is in reduced row echelon form. We then clear that corresponding new pivot column.
EXAMPLES::
sage: a = matrix(ZZ, 3, [1, 0, 110, 0, 3, 112, 0, 0, 221]); a [ 1 0 110] [ 0 3 112] [ 0 0 221] sage: a._add_row_and_maintain_echelon_form(vector(ZZ,[1,2,3]),[0,1,2]) ( [1 0 0] [0 1 0] [0 0 1], [0, 1, 2] ) sage: a._add_row_and_maintain_echelon_form(vector(ZZ,[0,0,0]),[0,1,2]) ( [ 1 0 110] [ 0 3 112] [ 0 0 221], [0, 1, 2] ) sage: a = matrix(ZZ, 2, [1, 0, 110, 0, 3, 112]) sage: a._add_row_and_maintain_echelon_form(vector(ZZ,[1,2,3]),[0,1]) ( [ 1 0 110] [ 0 1 219] [ 0 0 545], [0, 1, 2] ) """
# 0. Base case else: pivots = []
# 1. Create a new matrix that has row as the last row. # 2. Working from the left, clear each column to put # the resulting matrix back in echelon form. # p is the i-th pivot
# (a). Take xgcd of pivot positions in last row and in ith # row.
# TODO (optimize) -- change to use direct call to gmp and # no bounds checking! raise ZeroDivisionError("claimed pivot is not a pivot") # (b) Subtract a multiple of row i from row n. else: # (b). More elaborate. # Replace the ith row by s*A[i] + t*A[n], which will # have g in the i,p position, and replace the last row by # (b//g)*A[i] - (a//g)*A[n], which will have 0 in the i,p # position. raise ZeroDivisionError("claimed pivot is not a pivot (got a 0 gcd)")
# OK -- now we have to make sure the top part of the matrix # but with row i replaced by # r = s*row_i[j] + t*row_n[j] # is put in rref. We do this by recursively calling this # function with the top part of A (all but last row) and the # row r.
# 3. If it turns out that the last row is nonzero, # insert last row in A sliding other rows down.
# If pivot entry is negative negate this row.
# Determine where the last row should be inserted. # The new row should go *before* row j, so it becomes row j except RuntimeError: raise ZeroDivisionError("mistake in claimed pivots")
##################################################################################### # Hermite form modulo D # This code below is by E. Burcin. Thanks! ##################################################################################### def _hnf_mod(self, D): """ INPUT:
- ``D`` - a small integer that is assumed to be a multiple of 2\*det(self)
OUTPUT:
- ``matrix`` - the Hermite normal form of self. """
cdef int _hnf_modn(Matrix_integer_dense self, Matrix_integer_dense res, unsigned int det) except -1: """ Puts self into HNF form modulo det. Changes self in place. """ cdef int* res_l cdef Py_ssize_t i,j,k
cdef int* _hnf_modn_impl(Matrix_integer_dense self, unsigned int det, Py_ssize_t nrows, Py_ssize_t ncols) except NULL: # NOTE: det should be at most 2^31-1, such that anything modulo # det fits in a 32-bit signed integer. To avoid overflow, we # need 64-bit arithmetic for some intermediate computations. cdef int* res cdef int* T_ent cdef int* *res_rows cdef int** T_rows cdef int* B cdef Py_ssize_t i, j, k cdef int T_i_i, T_j_i, c1, c2, q, t cdef int u, v, d
# allocate memory for result matrix raise MemoryError("out of memory allocating a matrix") sig_free(res) raise MemoryError("out of memory allocating a matrix")
# allocate memory for temporary matrix sig_free(res) sig_free(res_rows) raise MemoryError("out of memory allocating a matrix") sig_free(res) sig_free(res_rows) sig_free(T_ent) raise MemoryError("out of memory allocating a matrix")
# allocate memory for temporary row vector sig_free(res) sig_free(res_rows) sig_free(T_ent) sig_free(T_rows) raise MemoryError("out of memory allocating a matrix")
# initialize the row pointers
cdef mpz_t tmp # copy entries from self to temporary matrix
# initialize variables
res_rows[i][i] = R
################################################################# # operations with matrices ################################################################# cdef _stack_impl(self, bottom): r""" Return the matrix ``self`` on top of ``bottom``::
[ self ] [ bottom ]
EXAMPLES::
sage: M = Matrix(ZZ, 2, 3, range(6)) sage: N = Matrix(ZZ, 1, 3, [10,11,12]) sage: M.stack(N) [ 0 1 2] [ 3 4 5] [10 11 12]
A vector may be stacked below a matrix. ::
sage: A = matrix(ZZ, 2, 4, range(8)) sage: v = vector(ZZ, 4, range(4)) sage: A.stack(v) [0 1 2 3] [4 5 6 7] [0 1 2 3]
The ``subdivide`` option will add a natural subdivision between ``self`` and ``bottom``. For more details about how subdivisions are managed when stacking, see :meth:`sage.matrix.matrix1.Matrix.stack`. ::
sage: A = matrix(ZZ, 3, 4, range(12)) sage: B = matrix(ZZ, 2, 4, range(8)) sage: A.stack(B, subdivide=True) [ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [-----------] [ 0 1 2 3] [ 4 5 6 7] """ cdef Matrix_integer_dense Z
cdef Py_ssize_t r, c
def augment(self, right, subdivide=False): r""" Returns a new matrix formed by appending the matrix (or vector) ``right`` on the right side of ``self``.
INPUT:
- ``right`` - a matrix, vector or free module element, whose dimensions are compatible with ``self``.
- ``subdivide`` - default: ``False`` - request the resulting matrix to have a new subdivision, separating ``self`` from ``right``.
OUTPUT:
A new matrix formed by appending ``right`` onto the right side of ``self``. If ``right`` is a vector (or free module element) then in this context it is appropriate to consider it as a column vector. (The code first converts a vector to a 1-column matrix.)
EXAMPLES::
sage: A = matrix(ZZ, 4, 5, range(20)) sage: B = matrix(ZZ, 4, 3, range(12)) sage: A.augment(B) [ 0 1 2 3 4 0 1 2] [ 5 6 7 8 9 3 4 5] [10 11 12 13 14 6 7 8] [15 16 17 18 19 9 10 11]
A vector may be augmented to a matrix. ::
sage: A = matrix(ZZ, 3, 5, range(15)) sage: v = vector(ZZ, 3, range(3)) sage: A.augment(v) [ 0 1 2 3 4 0] [ 5 6 7 8 9 1] [10 11 12 13 14 2]
The ``subdivide`` option will add a natural subdivision between ``self`` and ``right``. For more details about how subdivisions are managed when augmenting, see :meth:`sage.matrix.matrix1.Matrix.augment`. ::
sage: A = matrix(ZZ, 3, 5, range(15)) sage: B = matrix(ZZ, 3, 3, range(9)) sage: A.augment(B, subdivide=True) [ 0 1 2 3 4| 0 1 2] [ 5 6 7 8 9| 3 4 5] [10 11 12 13 14| 6 7 8]
Errors are raised if the sizes are incompatible. ::
sage: A = matrix(ZZ, [[1, 2],[3, 4]]) sage: B = matrix(ZZ, [[10, 20], [30, 40], [50, 60]]) sage: A.augment(B) Traceback (most recent call last): ... TypeError: number of rows must be the same, not 2 != 3 """ right = right.change_ring(self._base_ring)
cdef Matrix_integer_dense Z cdef Py_ssize_t i, j, p, qs, qa
def insert_row(self, Py_ssize_t index, row): """ Create a new matrix from self with.
INPUT:
- ``index`` - integer
- ``row`` - a vector
EXAMPLES::
sage: X = matrix(ZZ,3,range(9)); X [0 1 2] [3 4 5] [6 7 8] sage: X.insert_row(1, [1,5,-10]) [ 0 1 2] [ 1 5 -10] [ 3 4 5] [ 6 7 8] sage: X.insert_row(0, [1,5,-10]) [ 1 5 -10] [ 0 1 2] [ 3 4 5] [ 6 7 8] sage: X.insert_row(3, [1,5,-10]) [ 0 1 2] [ 3 4 5] [ 6 7 8] [ 1 5 -10] """ cdef Py_ssize_t j cdef Integer z cdef fmpz_t zflint raise ValueError("index must be nonnegative") raise ValueError("index must be less than number of rows")
def _delete_zero_columns(self): """ Return matrix obtained from self by deleting all zero columns along with the positions of those columns.
OUTPUT: matrix list of integers
EXAMPLES::
sage: a = matrix(ZZ, 2,3, [1,0,3,-1,0,5]); a [ 1 0 3] [-1 0 5] sage: a._delete_zero_columns() ( [ 1 3] [-1 5], [1] ) """
def _insert_zero_columns(self, cols): """ Return matrix obtained by self by inserting zero columns so that the columns with positions specified in cols are all 0.
INPUT:
- ``cols`` - list of nonnegative integers
OUTPUT: matrix
EXAMPLES::
sage: a = matrix(ZZ, 2,3, [1,0,3,-1,0,5]); a [ 1 0 3] [-1 0 5] sage: b, cols = a._delete_zero_columns() sage: b [ 1 3] [-1 5] sage: cols [1] sage: b._insert_zero_columns(cols) [ 1 0 3] [-1 0 5] """ # Now fill in the entries of A that come from self. # The following does this quickly: A[r,c] = self[r,i]
def _factor_out_common_factors_from_each_row(self): """ Very very quickly modifies self so that the gcd of the entries in each row is 1 by dividing each row by the common gcd.
EXAMPLES::
sage: a = matrix(ZZ, 3, [-9, 3, -3, -36, 18, -5, -40, -5, -5, -20, -45, 15, 30, -15, 180]) sage: a [ -9 3 -3 -36 18] [ -5 -40 -5 -5 -20] [-45 15 30 -15 180] sage: a._factor_out_common_factors_from_each_row() sage: a [ -3 1 -1 -12 6] [ -1 -8 -1 -1 -4] [ -3 1 2 -1 12] """
cdef fmpz_t g cdef fmpz_t tmp cdef long i, j
# divide through row
def gcd(self): """ Return the gcd of all entries of self; very fast.
EXAMPLES::
sage: a = matrix(ZZ,2, [6,15,-6,150]) sage: a.gcd() 3 """ cdef Py_ssize_t i, j cdef fmpz_t tmpgcd
def _change_ring(self, ring): """ Return the matrix obtained by coercing the entries of this matrix into the given ring.
EXAMPLES::
sage: a = matrix(ZZ,2,[1,-7,3,5]) sage: a._change_ring(RDF) [ 1.0 -7.0] [ 3.0 5.0] """ else:
def _singular_(self, singular=None): r""" Return Singular representation of this integer matrix.
INPUT:
- ``singular`` - Singular interface instance (default: None)
EXAMPLES::
sage: A = random_matrix(ZZ,3,3) sage: As = singular(A); As -8 2 0 0 1 -1 2 1 -95 sage: As.type() 'intmat' """ from sage.interfaces.singular import singular as singular_default singular = singular_default
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(ZZ, 2, 3, range(6)) sage: type(A) <type 'sage.matrix.matrix_integer_dense.Matrix_integer_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: 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_integer_dense A
def antitranspose(self): """ Returns the antitranspose of self, without changing self.
EXAMPLES::
sage: A = matrix(2,3,range(6)) sage: type(A) <type 'sage.matrix.matrix_integer_dense.Matrix_integer_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] """
cdef Matrix_integer_dense A cdef Py_ssize_t i,j cdef Py_ssize_t ri,rj # reversed i and j
def __pari__(self): """ Return PARI C-library version of this matrix.
EXAMPLES::
sage: a = matrix(ZZ,2,2,[1,2,3,4]) sage: a.__pari__() [1, 2; 3, 4] sage: pari(a) [1, 2; 3, 4] sage: type(pari(a)) <type 'cypari2.gen.Gen'> """
def _rank_pari(self): """ Rank of this matrix, computed using PARI. The computation is done entirely on the PARI stack, then the PARI stack is cleared. This function is most useful for very small matrices.
EXAMPLES:: sage: matrix(ZZ,3,[1..9])._rank_pari() 2 """
def _hnf_pari(self, int flag=0, bint include_zero_rows=True): """ Hermite form of this matrix, computed using PARI. The computation is done entirely on the PARI stack, then the PARI stack is cleared. This function is only useful for small matrices, and can crash on large matrices (e.g., if the PARI stack overflows).
INPUT:
- ``flag`` -- 0 (default), 1, 3 or 4 (see docstring for gp.mathnf).
- ``include_zero_rows`` -- boolean. if False, do not include any of the zero rows at the bottom of the matrix in the output.
.. NOTE::
In no cases is the transformation matrix returned by this function.
EXAMPLES::
sage: matrix(ZZ,3,[1..9])._hnf_pari() [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9])._hnf_pari(1) [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9])._hnf_pari(3) [1 2 3] [0 3 6] [0 0 0] sage: matrix(ZZ,3,[1..9])._hnf_pari(4) [1 2 3] [0 3 6] [0 0 0]
Check that ``include_zero_rows=False`` works correctly::
sage: matrix(ZZ,3,[1..9])._hnf_pari(0, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari(1, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari(3, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari(4, include_zero_rows=False) [1 2 3] [0 3 6]
Check that :trac:`12346` is fixed::
sage: pari('mathnf(Mat([0,1]), 4)') [Mat(1), [1, 0; 0, 1]] """ cdef GEN A
def _hnf_pari_big(self, int flag=0, bint include_zero_rows=True): """ Hermite form of this matrix, computed using PARI.
INPUT:
- ``flag`` -- 0 (default), 1, 3 or 4 (see docstring for gp.mathnf).
- ``include_zero_rows`` -- boolean. if False, do not include any of the zero rows at the bottom of the matrix in the output.
.. NOTE::
In no cases is the transformation matrix returned by this function.
EXAMPLES::
sage: a = matrix(ZZ,3,3,[1..9]) sage: a._hnf_pari_big(flag=0, include_zero_rows=True) [1 2 3] [0 3 6] [0 0 0] sage: a._hnf_pari_big(flag=1, include_zero_rows=True) [1 2 3] [0 3 6] [0 0 0] sage: a._hnf_pari_big(flag=3, include_zero_rows=True) [1 2 3] [0 3 6] [0 0 0] sage: a._hnf_pari_big(flag=4, include_zero_rows=True) [1 2 3] [0 3 6] [0 0 0]
Check that ``include_zero_rows=False`` works correctly::
sage: matrix(ZZ,3,[1..9])._hnf_pari_big(0, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari_big(1, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari_big(3, include_zero_rows=False) [1 2 3] [0 3 6] sage: matrix(ZZ,3,[1..9])._hnf_pari_big(4, include_zero_rows=False) [1 2 3] [0 3 6] """
cdef extract_hnf_from_pari_matrix(self, GEN H, int flag, bint include_zero_rows): # Throw away the transformation matrix (yes, we should later # code this to keep track of it). cdef mpz_t tmp
# Figure out how many columns we got back. # Now get the resulting Hermite form matrix back to Sage, suitably re-arranged. cdef Matrix_integer_dense B else:
def p_minimal_polynomials(self, p, s_max=None): r""" Compute `(p^s)`-minimal polynomials `\nu_s` of this matrix.
For `s \ge 0`, a `(p^s)`-minimal polynomial of a matrix `B` is a monic polynomial `f \in \ZZ[X]` of minimal degree such that all entries of `f(B)` are divisible by `p^s`.
Compute a finite subset `\mathcal{S}` of the positive integers and `(p^s)`-minimal polynomials `\nu_s` for `s \in \mathcal{S}`.
For `0 < t \le \max \mathcal{S}`, a `(p^t)`-minimal polynomial is given by `\nu_s` where `s = \min\{ r \in \mathcal{S} \mid r\ge t \}`. For `t > \max\mathcal{S}`, the minimal polynomial of `B` is also a `(p^t)`-minimal polynomial.
INPUT:
- ``p`` -- a prime in `\ZZ`
- ``s_max`` -- a positive integer (default: ``None``); if set, only `(p^s)`-minimal polynomials for ``s <= s_max`` are computed (see below for details)
OUTPUT:
A dictionary. Keys are the finite set `\mathcal{S}`, the values are the associated `(p^s)`-minimal polynomials `\nu_s`, `s\in\mathcal{S}`.
Setting ``s_max`` only affects the output if ``s_max`` is at most `\max\mathcal{S}` where `\mathcal{S}` denotes the full set. In that case, only those `\nu_s` with ``s <= s_max`` are returned where ``s_max`` is always included even if it is not included in the full set `\mathcal{S}`.
EXAMPLES::
sage: B = matrix(ZZ, [[1, 0, 1], [1, -2, -1], [10, 0, 0]]) sage: B.p_minimal_polynomials(2) {2: x^2 + 3*x + 2}
.. SEEALSO::
:mod:`~sage.matrix.compute_J_ideal`, :meth:`~sage.matrix.compute_J_ideal.ComputeMinimalPolynomials.p_minimal_polynomials` """
def null_ideal(self, b=0): r""" Return the `(b)`-ideal of this matrix.
Let `B` be a `n \times n` matrix. The *null ideal* modulo `b`, or `(b)`-ideal, is
.. MATH::
N_{(b)}(B) = \{f \in \ZZ[X] \mid f(B) \in M_n(b\ZZ)\}.
INPUT:
- ``b`` -- an element of `\ZZ` (default: 0)
OUTPUT:
An ideal in `\ZZ[X]`.
EXAMPLES::
sage: B = matrix(ZZ, [[1, 0, 1], [1, -2, -1], [10, 0, 0]]) sage: B.null_ideal() Principal ideal (x^3 + x^2 - 12*x - 20) of Univariate Polynomial Ring in x over Integer Ring sage: B.null_ideal(8) Ideal (8, x^3 + x^2 - 12*x - 20, 2*x^2 + 6*x + 4) of Univariate Polynomial Ring in x over Integer Ring sage: B.null_ideal(6) Ideal (6, 2*x^3 + 2*x^2 - 24*x - 40, 3*x^2 + 3*x) of Univariate Polynomial Ring in x over Integer Ring
.. SEEALSO::
:mod:`~sage.matrix.compute_J_ideal`, :meth:`~sage.matrix.compute_J_ideal.ComputeMinimalPolynomials.null_ideal` """
def integer_valued_polynomials_generators(self): r""" Determine the generators of the ring of integer valued polynomials on this matrix.
OUTPUT:
A pair ``(mu_B, P)`` where ``P`` is a list of polynomials in `\QQ[X]` such that
.. MATH::
\{f \in \QQ[X] \mid f(B) \in M_n(\ZZ)\} = \mu_B \QQ[X] + \sum_{g\in P} g \ZZ[X]
where `B` is this matrix.
EXAMPLES::
sage: B = matrix(ZZ, [[1, 0, 1], [1, -2, -1], [10, 0, 0]]) sage: B.integer_valued_polynomials_generators() (x^3 + x^2 - 12*x - 20, [1, 1/4*x^2 + 3/4*x + 1/2])
.. SEEALSO::
:mod:`~sage.matrix.compute_J_ideal`, :meth:`~sage.matrix.compute_J_ideal.ComputeMinimalPolynomials.integer_valued_polynomials_generators` """
cdef inline GEN pari_GEN(Matrix_integer_dense B): r""" Create the PARI GEN object on the stack defined by the integer matrix B. This is used internally by the function for conversion of matrices to PARI.
For internal use only; this directly uses the PARI stack. One should call ``sig_on()`` before and ``sig_off()`` after. """
#####################################################################################
cdef _clear_columns(Matrix_integer_dense A, pivots, Py_ssize_t n): # Clear all columns cdef fmpz_t c,t # subtract off c*v from row k; resulting A[k,i] entry will be < b, hence in Echelon form.
###############################################################
cpdef _lift_crt(Matrix_integer_dense M, residues, moduli=None): """ INPUT:
- ``M`` -- A ``Matrix_integer_dense``. Will be modified to hold the output.
- ``residues`` -- a list of ``Matrix_modn_dense_template``. The matrix to reconstruct modulo primes.
OUTPUT:
The matrix whose reductions modulo primes are the input ``residues``.
TESTS::
sage: from sage.matrix.matrix_integer_dense import _lift_crt sage: T1 = Matrix(Zmod(5), 4, 4, [1, 4, 4, 0, 2, 0, 1, 4, 2, 0, 4, 1, 1, 4, 0, 3]) sage: T2 = Matrix(Zmod(7), 4, 4, [1, 4, 6, 0, 2, 0, 1, 2, 4, 0, 6, 6, 1, 6, 0, 5]) sage: T3 = Matrix(Zmod(11), 4, 4, [1, 4, 10, 0, 2, 0, 1, 9, 8, 0, 10, 6, 1, 10, 0, 9]) sage: _lift_crt(Matrix(ZZ, 4, 4), [T1, T2, T3]) [ 1 4 -1 0] [ 2 0 1 9] [-3 0 -1 6] [ 1 -1 0 -2]
sage: from sage.arith.multi_modular import MultiModularBasis sage: mm = MultiModularBasis([5,7,11]) sage: _lift_crt(Matrix(ZZ, 4, 4), [T1, T2, T3], mm) [ 1 4 -1 0] [ 2 0 1 9] [-3 0 -1 6] [ 1 -1 0 -2]
The modulus must be smaller than the maximum for the multi-modular reconstruction (using ``mod_int``) and also smaller than the limit for ``Matrix_modn_dense_double`` to be able to represent the ``residues`` ::
sage: from sage.arith.multi_modular import MAX_MODULUS as MAX_multi_modular sage: from sage.matrix.matrix_modn_dense_double import MAX_MODULUS as MAX_modn_dense_double sage: MAX_MODULUS = min(MAX_multi_modular, MAX_modn_dense_double) sage: p0 = previous_prime(MAX_MODULUS) sage: p1 = previous_prime(p0) sage: mmod = [matrix(GF(p0), [[-1, 0, 1, 0, 0, 1, 1, 0, 0, 0, p0-1, p0-2]]), ....: matrix(GF(p1), [[-1, 0, 1, 0, 0, 1, 1, 0, 0, 0, p1-1, p1-2]])] sage: _lift_crt(Matrix(ZZ, 1, 12), mmod) [-1 0 1 0 0 1 1 0 0 0 -1 -2] """
cdef size_t i, j, k cdef Py_ssize_t nr, n return M
else: raise IndexError("Number of residues (%s) does not match number of moduli (%s)"%(len(residues), len(moduli)))
cdef MultiModularBasis mm
raise TypeError("Can only perform CRT on list of matrices mod n.")
cdef mod_int **row_list raise MemoryError("out of memory allocating multi-modular coefficient list")
raise MemoryError("out of memory allocating multi-modular coefficient list")
|