32.16 Dense matrices over the integer ring

Module: sage.matrix.matrix_integer_dense

Dense matrices over the integer ring.

Author Log:

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: loads(dumps(a)) == a
True

Module-level Functions

clear_mpz_globals( )

gmp_randrange( )

init_mpz_globals( )

tune_multiplication( )

Compare various multiplication algorithms.

Input:

k
- integer; affects numbers of trials
nmin
- integer; smallest matrix to use
nmax
- integer; largest matrix to use
bitmin
- integer; smallest bitsize
bitmax
- integer; largest bitsize

Output:
prints what doing then who wins
- multimodular or classical

sage: from sage.matrix.matrix_integer_dense import tune_multiplication
sage: tune_multiplication(2, nmin=10, nmax=60, bitmin=2,bitmax=8)
10 2 0.2
...

Class: Matrix_integer_dense

class Matrix_integer_dense
Matrix over the integers.

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.

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

Functions: charpoly,$ \,$ decomposition,$ \,$ determinant,$ \,$ echelon_form,$ \,$ elementary_divisors,$ \,$ frobenius,$ \,$ gcd,$ \,$ height,$ \,$ hermite_form,$ \,$ index_in_saturation,$ \,$ insert_row,$ \,$ kernel,$ \,$ kernel_matrix,$ \,$ LLL,$ \,$ LLL_gram,$ \,$ minpoly,$ \,$ pivots,$ \,$ prod_of_row_sums,$ \,$ randomize,$ \,$ rank,$ \,$ rational_reconstruction,$ \,$ saturation,$ \,$ smith_form,$ \,$ stack,$ \,$ symplectic_form

charpoly( )

Input:

var
- a variable name
algorithm
- 'linbox' (default) 'generic'

NOTE: Linbox charpoly disabled on 64-bit machines, since it hangs in many cases.

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: n=20; A = Mat(ZZ,n)(range(n^2))
sage: A.charpoly()
x^20 - 3990*x^19 - 266000*x^18
sage: A.minpoly()                # optional -- os x only right now
x^3 - 3990*x^2 - 266000*x

decomposition( )

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.

sage: t = ModularSymbols(11,sign=1).hecke_matrix(2)
sage: w = t.change_ring(ZZ)
sage: w.list()
[3, -1, 0, -2]

determinant( )

Return the determinant of this matrix.

Input: algorithm -

'padic'
- (the default) uses a p-adic / multimodular algorithm that relies on code in IML and linbox
'linbox'
- calls linbox det
'ntl'
- calls NTL's det function
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$ .

TIMINGS: This is perhaps the fastest implementation of determinants in the world. E.g., for a 500x500 random matrix with 32-bit entries on a core2 duo 2.6Ghz running OS X, Sage takes 4.12 seconds, whereas Magma takes 62.87 seconds (both with proof False). With proof=True on the same problem Sage takes 5.73 seconds. For another example, a 200x200 random matrix with 1-digit entries takes 4.18 seconds in pari, 0.18 in Sage with proof True, 0.11 in Sage with proof False, and 0.21 seconds in Magma with proof True and 0.18 in Magma with proof False.

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

echelon_form( )

Return the echelon form of this matrix over the integers also known as the hermit normal form (HNF).

Input: algorithm -

'padic'
- (default) a new fast p-adic modular algorithm,
'pari'
- use PARI with flag 1
'ntl'
- use NTL (only works for square matrices of full rank!)
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 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:
matrix
- the Hermite normal form (=echelon form over ZZ) of self.

NOTE: The result is not cached.

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]

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]

elementary_divisors( )

Return the elementary divisors of self, in order.

IMPLEMENTATION: Uses linbox, except sometimes linbox doesn't work (errors about pre-conditioning), in which case PARI is used.

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 this matrix. They are ordered in reverse by divisibility.

Input:

self
- matrix
algorithm
- (default: 'pari') 'pari': works robustless, but is slower.
'linbox'
- use linbox (currently off, broken)

Output: list of int's

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]

SEE ALSO: smith_form

frobenius( )

Return the Frobenius form (rational canonical form) of this matrix.

Input: flag -an integer:

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} F B$ .
var
- a string (default: 'x')

Input:

flag
- 0 (default), 1 or 2 as described above

ALGORITHM: uses pari's matfrobenius()

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: A.frobenius(2)
([ 0  0  0]
[ 1  0 18]
[ 0  1 12],
[    -1      2     -1]
[     0  23/15 -14/15]
[     0  -2/15   1/15])
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.

Author: 2006-04-02: Martin Albrecht

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.

gcd( )

Return the gcd of all entries of self; very fast.

sage: a = matrix(ZZ,2, [6,15,-6,150])
sage: a.gcd()
3

height( )

Return the height of this matrix, i.e., the max absolute value of the entries of the matrix.

Output: A nonnegative integer.

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

hermite_form( )

Return the Hermite normal form of self.

This is a synonym for self.echelon_form(...). See the documentation for self.echelon_form for more details.

sage: A = matrix(ZZ, 3, 5, [-1, -1, -2, 2, -2, -4, -19, -17, 1, 2, -3, 1, 1, -4, 1])
sage: E, U = A.hermite_form(transformation=True)
sage: E
[   1    0   52 -133  109]
[   0    1   19  -47   38]
[   0    0   69 -178  145]
sage: U
[-46   3  11]
[-16   1   4]
[-61   4  15]
sage: U*A
[   1    0   52 -133  109]
[   0    1   19  -47   38]
[   0    0   69 -178  145]
sage: A.hermite_form()
[   1    0   52 -133  109]
[   0    1   19  -47   38]
[   0    0   69 -178  145]

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]

index_in_saturation( )

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.

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]

insert_row( )

Create a new matrix from self with.

Input:

index
- integer
row
- a vector

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]

kernel( )

Return the left kernel of this matrix, as a module over the integers. This is the saturated ZZ-module spanned by all the row vectors v such that v*self = 0.

Input:

algorithm
- 'padic': a new p-adic based algorithm 'pari': use PARI
LLL
- bool (default: False); if True the basis is an LLL reduced basis; otherwise, it is an echelon basis.
proof
- None (default: proof.linear_algebra()); if False, impacts how determinants are computed.

By convention if self has 0 rows, the kernel is of dimension 0, whereas the kernel is the whole domain if self has 0 columns.

sage: M = MatrixSpace(ZZ,4,2)(range(8))
sage: M.kernel()
Free module of degree 4 and rank 2 over Integer Ring
Echelon basis matrix:
[ 1  0 -3  2]
[ 0  1 -2  1]

kernel_matrix( )

The options are exactly like self.kernel(...), but returns a matrix A whose rows form a basis for the left kernel, i.e., so that A*self = 0.

This is mainly useful to avoid all overhead associated with creating a free module.

sage: A = matrix(ZZ, 3, 3, [1..9])
sage: A.kernel_matrix()
[-1  2 -1]

Note that the basis matrix returned above is not in Hermite form.

sage: A.kernel()
Free module of degree 3 and rank 1 over Integer Ring
Echelon basis matrix:
[ 1 -2  1]

We compute another kernel:

sage: A = matrix(ZZ, 4, 2, [2, -1, 1, 1, -18, -1, -1, -5])
sage: K = A.kernel_matrix(); K
[-17 -20  -3   0]
[  7   3   1  -1]

K is a basis for the left kernel:

sage: K*A
[0 0]
[0 0]

We illustrate the LLL flag:

sage: L = A.kernel_matrix(LLL=True); L
[  7   3   1  -1]
[  4 -11   0  -3]
sage: K.hermite_form()
[ 1 64  3 12]
[ 0 89  4 17]
sage: L.hermite_form()
[ 1 64  3 12]
[ 0 89  4 17]

LLL( )

Returns 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:

(a) For any $ i>j$ , we have $ \vert mu_{i, j}\vert <= \eta$ , (b) For any $ i<d$ , we have $ \delta \vert b_i^*\vert^2 <= \vert b_{i + 1}^* + mu_{i + 1, i} b_{i + 1}^* \vert^2$ ,

where $ mu_{i,j} = <b_i, b_j^*>/<b_j^*,b_j^*>$ 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 <= 1.0$ and $ 0.5 <= \eta <
sqrt(\delta)$ . Polynomial time complexity is only guaranteed for $ \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
- parameter as described above (default: 3/4)
eta
- parameter as described above (default: 0.501), ignored by NTL
algorithm
- string (default: "fpLLL:wrapper") one of the algorithms mentioned below
fp
- None - NTL's exact reduction or fpLLL's wrapper
'fp'
- double precision: NTL's FP or fpLLL's double
'qd'
- quad doubles: NTL's QP
'xd'
- extended exponent: NTL's XD or fpLLL's dpe
'rr'
- arbitrary precision: NTL'RR or fpLLL's MPFR
prec
- precision, ignored by NTL (default: auto choose)
early_red
- perform early reduction, ignored by NTL (default: False)
use_givens
- use Givens orthogonalization (default: False) only applicable to approximate reductions and NTL. This is more stable but slower.

Also, if the verbose level is >= 2, some more verbose output is printed during the calculation if NTL is used.

AVAILABLE ALGORITHMS:

NTL:LLL
- NTL's LLL + fp
fpLLL:heuristic
- fpLLL's heuristic + fp
fpLLL:fast
- fpLLL's fast
fpLLL:wrapper
- fpLLL's automatic choice (default)

Output: a matrix over the integers

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 xrange(n):
...       X[i,i + 1] = 1
sage: for i in xrange(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

ALGORITHM: Uses the NTL library by Victor Shoup or fpLLL library by Damien Stehle depending on the chosen algorithm.

REFERENCES: ntl.mat_ZZ or sage.libs.fplll.fplll for details on the used algorithms.

LLL_gram( )

LLL reduction of the lattice whose gram matrix is self.

Input:

M
- gram matrix of a definite quadratic form

Output:
U
- unimodular transformation matrix such that

U.transpose() * M * U

is LLL-reduced

ALGORITHM: Use PARI

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]

Semidefinite and indefinite forms raise a ValueError:

sage: Matrix(ZZ,2,2,[2,6,6,3]).LLL_gram()
Traceback (most recent call last):
...
ValueError: not a definite matrix
sage: Matrix(ZZ,2,2,[1,0,0,-1]).LLL_gram()
Traceback (most recent call last):
...
ValueError: not a definite matrix

BUGS: should work for semidefinite forms (PARI is ok)

minpoly( )

Input:

var
- a variable name
algorithm
- 'linbox' (default) 'generic'

NOTE: Linbox charpoly disabled on 64-bit machines, since it hangs in many cases.

sage: A = matrix(ZZ,6, range(36))
sage: A.minpoly()           # optional -- os x only right now
x^3 - 105*x^2 - 630*x
sage: n=6; A = Mat(ZZ,n)([k^2 for k in range(n^2)])
sage: A.minpoly()           # optional -- os x only right now
x^4 - 2695*x^3 - 257964*x^2 + 1693440*x

pivots( )

Return the pivot column positions of this matrix as a list of Python integers.

This returns a list, of the position of the first nonzero entry in each row of the echelon form.

Output:

list
- a list of Python ints

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]

prod_of_row_sums( )

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

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

randomize( )

Randomize density proportion of the entries of this matrix, leaving the rest unchanged.

The parameters are the same as the integer ring's random_element function.

If x and y are given, randomized entries of this matrix 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 are passed to 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

Output:
- modifies this matrix in place

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]

rank( )

Return the rank of this matrix.

Output:

nonnegative integer
- the rank

NOTE: The rank is cached.

ALGORITHM: First check if the matrix has maxim posible rank by working modulo one random prime. If not call Linbox's rank function.

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's a bigger example - the rank is of course still 2:

sage: a = matrix(ZZ,100,[1..100^2]); a.rank()
2

rational_reconstruction( )

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

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 of 4 (mod 11) 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

saturation( )

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 $ \mathbf{Z}$ module $ M$ embedded in $ \mathbf{Z}^n$ is the a module $ S$ that contains $ M$ with finite index such that $ \mathbf{Z}^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 $ \mathbf{Z}$ -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.

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)

smith_form( )

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.

sage: A = MatrixSpace(IntegerRing(), 3)(range(9))
sage: D, U, V = A.smith_form()
sage: D
[0 0 0]
[0 3 0]
[0 0 1]
sage: U
[-1  2 -1]
[ 0 -1  1]
[ 0  1  0]
sage: V
[ 1  4 -1]
[-2 -3  1]
[ 1  0  0]
sage: U*A*V
[0 0 0]
[0 3 0]
[0 0 1]

It also makes sense for nonsquare matrices:

sage: A = Matrix(ZZ,3,2,range(6))
sage: D, U, V = A.smith_form()
sage: D
[0 0]
[2 0]
[0 1]
sage: U
[-1  2 -1]
[ 0 -1  1]
[ 0  1  0]
sage: V
[ 3 -1]
[-2  1]
sage: U * A * V
[0 0]
[2 0]
[0 1]

SEE ALSO: elementary_divisors

stack( )

Return the matrix self on top of other: [ self ] [ other ]

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]

symplectic_form( )

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 $ z_1, \ldots, z_i, e_1,
\ldots, e_j, f_1, \ldots f_j$ 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 = d_i$ for all $ i$ , and $ d_{i+1} \vert d_{i}$ for all $ i$ ; * $ e_i M {f_j}^t = 0$ for all $ i$ not equal $ j$ .

The ordering for the factors $ d_{i+1} \vert d_{i}$ 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.

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  0  0  0]
[ 0  0  0  2  0]
[ 0  0  0  0  1]
[ 0 -2  0  0  0]
[ 0  0 -1  0  0]
sage: F == C * E * C.transpose()
True
sage: E.smith_form()[0]
[0 0 0 0 0]
[0 2 0 0 0]
[0 0 2 0 0]
[0 0 0 1 0]
[0 0 0 0 1]

Special Functions: __copy__,$ \,$ __eq__,$ \,$ __ge__,$ \,$ __gt__,$ \,$ __init__,$ \,$ __le__,$ \,$ __lt__,$ \,$ __ne__,$ \,$ _add_row_and_maintain_echelon_form,$ \,$ _adjoint,$ \,$ _change_ring,$ \,$ _charpoly_linbox,$ \,$ _delete_zero_columns,$ \,$ _det_linbox,$ \,$ _det_ntl,$ \,$ _echelon_in_place_classical,$ \,$ _echelon_strassen,$ \,$ _elementary_divisors_linbox,$ \,$ _factor_out_common_factors_from_each_row,$ \,$ _hnf_mod,$ \,$ _insert_zero_columns,$ \,$ _invert_iml,$ \,$ _kernel_gens_using_pari,$ \,$ _kernel_matrix_using_padic_algorithm,$ \,$ _linbox_dense,$ \,$ _linbox_modn,$ \,$ _linbox_modn_det,$ \,$ _linbox_sparse,$ \,$ _minpoly_linbox,$ \,$ _mod_int,$ \,$ _multiply_classical,$ \,$ _multiply_linbox,$ \,$ _multiply_multi_modular,$ \,$ _ntl_,$ \,$ _pickle,$ \,$ _poly_linbox,$ \,$ _rank_linbox,$ \,$ _rank_modp,$ \,$ _rational_echelon_via_solve,$ \,$ _rational_kernel_iml,$ \,$ _reduce,$ \,$ _solve_iml,$ \,$ _solve_right_nonsingular_square,$ \,$ _unpickle

__copy__( )

Returns a new copy of this matrix.

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

_add_row_and_maintain_echelon_form( )

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.

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])

_adjoint( )

Return the adjoint of this matrix.

Assumes self is a square matrix (checked in adjoint).

_change_ring( )

Return the matrix obtained by coercing the entries of this matrix into the given ring.

sage: a = matrix(ZZ,2,[1,-7,3,5])
sage: a._change_ring(RDF)
[ 1.0 -7.0]
[ 3.0  5.0]

_charpoly_linbox( )

_delete_zero_columns( )

Return matrix obtained from self by deleting all zero columns along with the positions of those columns.

Output: matrix list of integers

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])

_det_linbox( )

Compute the determinant of this matrix using Linbox.

_det_ntl( )

Compute the determinant of this matrix using NTL.

_echelon_in_place_classical( )

_echelon_strassen( )

_elementary_divisors_linbox( )

_factor_out_common_factors_from_each_row( )

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.

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]

_hnf_mod( )

Input:

D
- a small integer that is assumed to be a multiple of 2*det(self)
Output:
matrix
- the Hermite normal form of self.

_insert_zero_columns( )

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

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]

_invert_iml( )

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.

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]

Author: William Stein

_kernel_gens_using_pari( )

Compute an LLL reduced list of independent generators that span the kernel of self.

ALGORITHM: Call pari's matkerint function.

_kernel_matrix_using_padic_algorithm( )

Compute a list of independent generators that span the right kernel of self.

ALGORITHM: Use IML to compute the kernel over QQ, clear denominators, then saturate.

_linbox_dense( )

_linbox_modn( )

Return modn linbox object associated to this integer matrix.

sage: a = matrix(ZZ, 3, [1,2,5,-7,8,10,192,5,18])
sage: b = a._linbox_modn(19); b
<sage.libs.linbox.linbox.Linbox_modn_dense object at ...>
sage: b.charpoly()
[2L, 10L, 11L, 1L]

_linbox_modn_det( )

Input:

n
- a prime (at most 67108879)

sage: a = matrix(ZZ, 3, [1,2,5,-7,8,10,192,5,18])
sage: a.det()
-3669
sage: a._linbox_modn_det(5077)
1408
sage: a._linbox_modn_det(3)
0
sage: a._linbox_modn_det(2)
1
sage: a.det()%5077
1408
sage: a.det()%2
1
sage: a.det()%3
0

_linbox_sparse( )

_minpoly_linbox( )

_mod_int( )

_multiply_classical( )

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]

_multiply_linbox( )

Multiply matrices over ZZ using linbox.

WARNING: This is very slow right now, i.e., linbox is very slow.

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]

_multiply_multi_modular( )

_ntl_( )

ntl.mat_ZZ representation of self.

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.

_pickle( )

sage: S = ModularSymbols(250,4,sign=1).cuspidal_submodule().new_subspace().decomposition() # long
sage: S == loads(dumps(S)) # long
True

_poly_linbox( )

Input:

var
- 'x'
typ
- 'minpoly' or 'charpoly'

_rank_linbox( )

Compute the rank of this matrix using Linbox.

_rank_modp( )

_rational_echelon_via_solve( )

Computes information that gives the reduced row echelon form (over QQ!) of a matrix with integer entries.

Input:

self
- a matrix over the integers.

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.

Author: 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.

1
Compute r = rank(A).

2
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).

3
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.

4
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\vert 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.

5
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.

6
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.

_rational_kernel_iml( )

IML: Return the rational kernel of this matrix (acting from the left), considered as a matrix over QQ. I.e., returns a matrix K such that self*K = 0, and the number of columns of K equals the nullity of self.

Author: William Stein

_reduce( )

_solve_iml( )

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

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

ALGORITHM: Uses IML.

Author: Martin Albrecht

_solve_right_nonsingular_square( )

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.
Output: a matrix or vector over $ \mathbf{Q}$

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 $ \mathbf{Q}$ .

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: 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

_unpickle( )

See About this document... for information on suggesting changes.