Context Navigation

source: sage/matrix/matrix_integer_dense.pyx @ 2830:bf1870252c66

Visit:

Revision 2830:bf1870252c66, 45.0 KB checked in by William Stein <wstein@…>, 6 years ago (diff)
linear algebra fixes for linux for now

Line
1	"""
2	Dense matrices over the integer ring.
3	"""
4
5	######################################################################
6	# Copyright (C) 2006,2007 William Stein
7	#
8	# Distributed under the terms of the GNU General Public License (GPL)
9	#
10	# The full text of the GPL is available at:
11	# http://www.gnu.org/licenses/
12	######################################################################
13
14	from sage.misc.misc import verbose, get_verbose, UNAME
15
16	include "../ext/interrupt.pxi"
17	include "../ext/stdsage.pxi"
18	include "../ext/gmp.pxi"
19
20	cdef extern from "matrix_integer_dense_linbox.h":
21	void linbox_integer_dense_minpoly_hacked(mpz_t* minpoly, size_t degree,
22	size_t n, mpz_t** matrix, int do_minpoly)
23	void linbox_integer_dense_minpoly(mpz_t* minpoly, size_t degree,
24	size_t n, mpz_t** matrix)
25	void linbox_integer_dense_charpoly(mpz_t* charpoly, size_t degree,
26	size_t n, mpz_t** matrix)
27	void linbox_integer_dense_delete_array(mpz_t* f)
28	int linbox_integer_dense_matrix_matrix_multiply(mpz_t ans, mpz_t A, mpz_t **B,
29	size_t A_nr, size_t A_nc, size_t B_nr, size_t B_nc)
30
31	ctypedef unsigned int uint
32
33	from sage.ext.multi_modular cimport MultiModularBasis
34	cdef MultiModularBasis mm
35	mm = MultiModularBasis()
36
37	from sage.rings.integer cimport Integer
38	from sage.rings.rational_field import QQ
39	from sage.rings.integer_ring import ZZ
40	from sage.rings.integer_mod_ring import IntegerModRing
41	from sage.rings.polynomial_ring import PolynomialRing
42	from sage.structure.element cimport ModuleElement, RingElement
43
44	from matrix_modn_dense import Matrix_modn_dense
45	from matrix_modn_dense cimport Matrix_modn_dense
46
47	import sage.modules.free_module
48
49	from matrix cimport Matrix
50
51	import matrix_space
52
53
54	cdef class Matrix_integer_dense(matrix_dense.Matrix_dense): # dense or sparse
55	r"""
56	Matrix over the integers.
57
58	On a 32-bit machine, they can have at most $2^{32}-1$ rows or
59	columns. On a 64-bit machine, matrices can have at most
60	$2^{64}-1$ rows or columns.
61
62	EXAMPLES:
63	sage: a = MatrixSpace(ZZ,3)(2); a
64	[2 0 0]
65	[0 2 0]
66	[0 0 2]
67	sage: a = matrix(ZZ,1,3, [1,2,-3]); a
68	[ 1 2 -3]
69	sage: a = MatrixSpace(ZZ,2,4)(2); a
70	Traceback (most recent call last):
71	...
72	TypeError: nonzero scalar matrix must be square
73	"""
74	########################################################################
75	# LEVEL 1 functionality
76	# x * __new__
77	# x * __dealloc__
78	# x * __init__
79	# x * set_unsafe
80	# x * get_unsafe
81	# x * def _pickle
82	# x * def _unpickle
83	########################################################################
84
85	def __new__(self, parent, entries, coerce, copy):
86	"""
87	Create and allocate memory for the matrix. Does not actually initialize
88	any of the memory.
89
90	INPUT:
91	parent, entries, coerce, copy -- as for __init__.
92
93	EXAMPLES:
94	sage: from sage.matrix.matrix_integer_dense import Matrix_integer_dense
95	sage: a = Matrix_integer_dense.__new__(Matrix_integer_dense, Mat(ZZ,3), 0,0,0)
96	sage: type(a)
97	<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
98
99	WARNING: This is for internal use only, or if you really know what you're doing.
100	"""
101	matrix_dense.Matrix_dense.__init__(self, parent)
102	self._initialized = 0
103	self._nrows = parent.nrows()
104	self._ncols = parent.ncols()
105	self._pivots = None
106
107	# Allocate an array where all the entries of the matrix are stored.
108	_sig_on
109	self._entries = <mpz_t >sage_malloc(sizeof(mpz_t) (self._nrows * self._ncols))
110	_sig_off
111	if self._entries == NULL:
112	raise MemoryError, "out of memory allocating a matrix"
113
114	# Allocate an array of pointers to the rows, which is useful for
115	# certain algorithms.
116	self._matrix = <mpz_t *> sage_malloc(sizeof(mpz_t)*self._nrows)
117	if self._matrix == NULL:
118	sage_free(self._entries)
119	self._entries = NULL
120	raise MemoryError, "out of memory allocating a matrix"
121
122	# Set each of the pointers in the array self._matrix to point
123	# at the memory for the corresponding row.
124	cdef Py_ssize_t i, k
125	k = 0
126	for i from 0 <= i < self._nrows:
127	self._matrix[i] = self._entries + k
128	k = k + self._ncols
129
130	def __copy__(self):
131	cdef Matrix_integer_dense A
132	A = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent,
133	0, 0, 0)
134	cdef Py_ssize_t i
135	_sig_on
136	for i from 0 <= i < self._nrows * self._ncols:
137	mpz_init_set(A._entries[i], self._entries[i])
138	_sig_off
139	return A
140
141	def __dealloc__(self):
142	"""
143	Frees all the memory allocated for this matrix.
144	EXAMPLE:
145	sage: a = Matrix(ZZ,2,[1,2,3,4])
146	sage: del a
147	"""
148	if self._entries == NULL: return
149	cdef Py_ssize_t i
150	if self._initialized:
151	for i from 0 <= i < (self._nrows * self._ncols):
152	mpz_clear(self._entries[i])
153	sage_free(self._entries)
154	sage_free(self._matrix)
155
156	def __init__(self, parent, entries, copy, coerce):
157	r"""
158	Initialize a dense matrix over the integers.
159
160	INPUT:
161	parent -- a matrix space
162	entries -- list - create the matrix with those entries along the rows.
163	other -- a scalar; entries is coerced to an integer and the diagonal
164	entries of this matrix are set to that integer.
165	coerce -- whether need to coerce entries to the integers (program may crash
166	if you get this wrong)
167	copy -- ignored (since integers are immutable)
168
169	EXAMPLES:
170
171	The __init__ function is called implicitly in each of the
172	examples below to actually fill in the values of the matrix.
173
174	We create a $2 \times 2$ and a $1\times 4$ matrix:
175	sage: matrix(ZZ,2,2,range(4))
176	[0 1]
177	[2 3]
178	sage: Matrix(ZZ,1,4,range(4))
179	[0 1 2 3]
180
181	If the number of columns isn't given, it is determined from the number of
182	elements in the list.
183	sage: matrix(ZZ,2,range(4))
184	[0 1]
185	[2 3]
186	sage: matrix(ZZ,2,range(6))
187	[0 1 2]
188	[3 4 5]
189
190	Another way to make a matrix is to create the space of
191	matrices and coerce lists into it.
192	sage: A = Mat(ZZ,2); A
193	Full MatrixSpace of 2 by 2 dense matrices over Integer Ring
194	sage: A(range(4))
195	[0 1]
196	[2 3]
197
198	Actually it is only necessary that the input can be coerced to a list, so
199	the following also works:
200	sage: v = reversed(range(4)); type(v)
201	<type 'listreverseiterator'>
202	sage: A(v)
203	[3 2]
204	[1 0]
205
206	Matrices can have many rows or columns (in fact, on a 64-bit machine they could
207	have up to $2^64-1$ rows or columns):
208	sage: v = matrix(ZZ,1,10^5, range(10^5))
209	sage: v.parent()
210	Full MatrixSpace of 1 by 100000 dense matrices over Integer Ring
211	"""
212	cdef Py_ssize_t i, j
213	cdef int is_list
214	cdef Integer x
215
216	if not isinstance(entries, list): # todo -- change to PyObject_TypeCheck???
217	try:
218	entries = list(entries)
219	is_list = 1
220	except TypeError:
221	try:
222	# Try to coerce entries to a scalar (an integer)
223	x = ZZ(entries)
224	is_list = 0
225	except TypeError:
226	raise TypeError, "entries must be coercible to a list or integer"
227	else:
228	is_list = 1
229
230	if is_list:
231
232	# Create the matrix whose entries are in the given entry list.
233	if len(entries) != self._nrows * self._ncols:
234	raise TypeError, "entries has the wrong length"
235	if coerce:
236	for i from 0 <= i < self._nrows * self._ncols:
237	# TODO: Should use an unsafe un-bounds-checked array access here.
238	x = ZZ(entries[i])
239	# todo -- see integer.pyx and the TODO there; perhaps this could be
240	# sped up by creating a mpz_init_set_sage function.
241	mpz_init_set(self._entries[i], x.value)
242	else:
243	for i from 0 <= i < self._nrows * self._ncols:
244	# TODO: Should use an unsafe un-bounds-checked array access here.
245	mpz_init_set(self._entries[i], (<Integer> entries[i]).value)
246
247	else:
248
249	# If x is zero, make the zero matrix and be done.
250	if mpz_cmp_si(x.value, 0) == 0:
251	self._zero_out_matrix()
252	return
253
254	# the matrix must be square:
255	if self._nrows != self._ncols:
256	sage_free(self._entries)
257	sage_free(self._matrix)
258	self._entries = NULL
259	raise TypeError, "nonzero scalar matrix must be square"
260
261	# Now we set all the diagonal entries to x and all other entries to 0.
262	self._zero_out_matrix()
263	j = 0
264	for i from 0 <= i < self._nrows:
265	mpz_init_set(self._entries[j], x.value)
266	j = j + self._nrows + 1
267	self._initialized = 1
268
269
270	cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
271	"""
272	Set position i,j of this matrix to x.
273
274	(VERY UNSAFE -- value must be of type Integer).
275
276	INPUT:
277	ij -- tuple (i,j), where i is the row and j the column
278	Alternatively, ij can be an integer, and the ij-th row is set.
279
280	EXAMPLES:
281	sage: a = matrix(ZZ,2,3, range(6)); a
282	[0 1 2]
283	[3 4 5]
284	sage: a[0,0] = 10
285	sage: a
286	[10 1 2]
287	[ 3 4 5]
288	"""
289	#cdef Integer Z
290	#Z = value
291	#mpz_set(self._matrix[i][j], Z.value)
292	mpz_set(self._matrix[i][j], (<Integer>value).value)
293
294	cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
295	"""
296	EXAMPLES:
297	sage: a = MatrixSpace(ZZ,3)(range(9)); a
298	[0 1 2]
299	[3 4 5]
300	[6 7 8]
301	sage: a[1,2]
302	5
303	sage: a[0]
304	(0, 1, 2)
305	sage: a[4,7]
306	Traceback (most recent call last):
307	...
308	IndexError: matrix index out of range
309	sage: a[-1,0]
310	Traceback (most recent call last):
311	...
312	IndexError: matrix index out of range
313	"""
314	cdef Integer z
315	z = PY_NEW(Integer)
316	mpz_set(z.value, self._matrix[i][j])
317	return z
318
319	def _pickle(self):
320	return self._pickle_version0(), 0
321
322	cdef _pickle_version0(self):
323	# TODO: maybe redo this to use mpz_import and mpz_export
324	# from sec 5.14 of the GMP manual. ??
325	cdef int i, j, len_so_far, m, n
326	cdef char *a
327	cdef char s, t, *tmp
328
329	if self._nrows == 0 or self._ncols == 0:
330	data = ''
331	else:
332	n = self._nrowsself._ncols10
333	s = <char> sage_malloc(n sizeof(char))
334	t = s
335	len_so_far = 0
336
337	_sig_on
338	for i from 0 <= i < self._nrows * self._ncols:
339	m = mpz_sizeinbase (self._entries[i], 32)
340	if len_so_far + m + 1 >= n:
341	# copy to new string with double the size
342	n = 2*n + m + 1
343	tmp = <char> sage_malloc(n sizeof(char))
344	strcpy(tmp, s)
345	sage_free(s)
346	s = tmp
347	t = s + len_so_far
348	#endif
349	mpz_get_str(t, 32, self._entries[i])
350	m = strlen(t)
351	len_so_far = len_so_far + m + 1
352	t = t + m
353	t[0] = <char>32
354	t[1] = <char>0
355	t = t + 1
356	_sig_off
357	data = str(s)[:-1]
358	free(s)
359	return data
360
361	def _unpickle(self, data, int version):
362	if version == 0:
363	self._unpickle_version0(data)
364	else:
365	raise RuntimeError, "unknown matrix version (=%s)"%version
366
367	cdef _unpickle_version0(self, data):
368	cdef Py_ssize_t i, n
369	data = data.split()
370	n = self._nrows * self._ncols
371	if len(data) != n:
372	raise RuntimeError, "invalid pickle data."
373	for i from 0 <= i < n:
374	s = data[i]
375	if mpz_init_set_str(self._entries[i], s, 32):
376	raise RuntimeError, "invalid pickle data"
377
378
379	def __richcmp__(Matrix self, right, int op): # always need for mysterious reasons.
380	return self._richcmp(right, op)
381
382	def __hash__(self):
383	return self._hash()
384
385	########################################################################
386	# LEVEL 1 helpers:
387	# These function support the implementation of the level 1 functionality.
388	########################################################################
389	cdef _zero_out_matrix(self):
390	"""
391	Set this matrix to be the zero matrix.
392	This is only for internal use.
393	"""
394	# TODO: This is about 6-10 slower than MAGMA doing what seems to be the same thing.
395	# Moreover, NTL can also do this quickly. Why? I think both have specialized
396	# small integer classes.
397	_sig_on
398	cdef Py_ssize_t i
399	for i from 0 <= i < self._nrows * self._ncols:
400	mpz_init(self._entries[i])
401	_sig_off
402	self._initialized = 1
403
404	cdef _new_unitialized_matrix(self, Py_ssize_t nrows, Py_ssize_t ncols):
405	"""
406	Return a new matrix over the integers with the given number of rows and columns.
407	All memory is allocated for this matrix, but its entries have not yet been
408	filled in.
409	"""
410	P = self._parent.matrix_space(nrows, ncols)
411	return Matrix_integer_dense.__new__(Matrix_integer_dense, P, None, None, None)
412
413
414	########################################################################
415	# LEVEL 2 functionality
416	# x * cdef _add_c_impl
417	# x * cdef _sub_c_impl
418	# x * cdef _mul_c_impl
419	# * cdef _cmp_c_impl
420	# * __neg__
421	# * __invert__
422	# * __copy__
423	# x * _multiply_classical
424	# * _list -- list of underlying elements (need not be a copy)
425	# * _dict -- sparse dictionary of underlying elements (need not be a copy)
426	########################################################################
427
428	# cdef _mul_c_impl(self, Matrix right):
429	# cdef int _cmp_c_impl(self, Matrix right) except -2:
430	# def __neg__(self):
431	# def __invert__(self):
432	# def __copy__(self):
433	# def _multiply_classical(left, matrix.Matrix _right):
434	# def _list(self):
435	# def _dict(self):
436
437	def _multiply_linbox(self, Matrix right):
438	"""
439	Multiply matrices over ZZ using linbox.
440
441	WARNING: This is very slow right now, i.e., linbox is very slow.
442
443	EXAMPLES:
444	sage: A = matrix(ZZ,2,3,range(6))
445	sage: A*A.transpose()
446	[ 5 14]
447	[14 50]
448	sage: A._multiply_linbox(A.transpose())
449	[ 5 14]
450	[14 50]
451	"""
452	cdef int e
453	cdef Matrix_integer_dense ans, B
454	ans = self.new_matrix(nrows = self.nrows(), ncols = right.ncols())
455	B = right
456	_sig_on
457	e = linbox_integer_dense_matrix_matrix_multiply(ans._matrix, self._matrix, B._matrix,
458	self._nrows, self._ncols,
459	right._nrows, right._ncols)
460	_sig_off
461	if e:
462	raise RuntimeError
463	return ans
464
465
466	def _multiply_classical(self, Matrix right):
467	"""
468	EXAMPLE:
469	sage: n = 3
470	sage: a = MatrixSpace(ZZ,n,n)(range(n^2))
471	sage: a*a
472	[ 15 18 21]
473	[ 42 54 66]
474	[ 69 90 111]
475	"""
476	if self._ncols != right._nrows:
477	raise IndexError, "Number of columns of self must equal number of rows of right."
478
479	cdef Py_ssize_t i, j, k, l, nr, nc, snc
480	cdef mpz_t *v
481
482	nr = self._nrows
483	nc = right._ncols
484	snc = self._ncols
485
486	parent = self.matrix_space(nr, nc)
487
488	cdef Matrix_integer_dense M, _right
489	_right = right
490
491	M = Matrix_integer_dense.__new__(Matrix_integer_dense, parent, None, None, None)
492	Matrix.__init__(M, parent)
493
494	# M has memory allocated but entries are not initialized
495	cdef mpz_t *entries
496	entries = M._entries
497
498	cdef mpz_t s, z
499	mpz_init(s)
500	mpz_init(z)
501
502	_sig_on
503	l = 0
504	for i from 0 <= i < nr:
505	for j from 0 <= j < nc:
506	mpz_set_si(s,0) # set s = 0
507	v = self._matrix[i]
508	for k from 0 <= k < snc:
509	mpz_mul(z, v[k], _right._matrix[k][j])
510	mpz_add(s, s, z)
511	mpz_init_set(entries[l], s)
512	l += 1
513	_sig_off
514	mpz_clear(s)
515	mpz_clear(z)
516	return M
517
518	cdef ModuleElement _lmul_c_impl(self, RingElement right):
519	"""
520	EXAMPLES:
521	sage: a = matrix(QQ,2,range(6))
522	sage: (3/4) * a
523	[ 0 3/4 3/2]
524	[ 9/4 3 15/4]
525	"""
526	cdef Py_ssize_t i
527	cdef Integer _x
528	_x = Integer(right)
529	cdef Matrix_integer_dense M
530	M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None)
531	for i from 0 <= i < self._nrows * self._ncols:
532	mpz_init(M._entries[i])
533	mpz_mul(M._entries[i], self._entries[i], _x.value)
534	M._initialized = 1
535	return M
536
537	cdef ModuleElement _add_c_impl(self, ModuleElement right):
538	"""
539	Add two dense matrices over ZZ.
540
541	EXAMPLES:
542	sage: a = MatrixSpace(ZZ,3)(range(9))
543	sage: a+a
544	[ 0 2 4]
545	[ 6 8 10]
546	[12 14 16]
547	sage: b = MatrixSpace(ZZ,3)(range(9))
548	sage: b.swap_rows(1,2)
549	sage: a+b
550	[ 0 2 4]
551	[ 9 11 13]
552	[ 9 11 13]
553	"""
554	cdef Py_ssize_t i
555
556	cdef Matrix_integer_dense M
557	M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None)
558	Matrix.__init__(M, self._parent)
559
560	_sig_on
561
562	cdef mpz_t *entries
563	entries = M._entries
564	for i from 0 <= i < self._ncols * self._nrows:
565	mpz_init(entries[i])
566	mpz_add(entries[i], self._entries[i], (<Matrix_integer_dense> right)._entries[i])
567	_sig_off
568	return M
569
570	cdef ModuleElement _sub_c_impl(self, ModuleElement right):
571	"""
572	Subtract two dense matrices over ZZ.
573
574	EXAMPLES:
575	sage: M = Mat(ZZ,3)
576	sage: a = M(range(9)); b = M(reversed(range(9)))
577	sage: a - b
578	[-8 -6 -4]
579	[-2 0 2]
580	[ 4 6 8]
581
582	"""
583	cdef Py_ssize_t i
584
585	cdef Matrix_integer_dense M
586	M = Matrix_integer_dense.__new__(Matrix_integer_dense, self._parent, None, None, None)
587	Matrix.__init__(M, self._parent)
588
589	_sig_on
590
591	cdef mpz_t *entries
592	entries = M._entries
593	for i from 0 <= i < self._ncols * self._nrows:
594	mpz_init(entries[i])
595	mpz_sub(entries[i], self._entries[i], (<Matrix_integer_dense> right)._entries[i])
596
597	_sig_off
598	return M
599
600
601
602	########################################################################
603	# LEVEL 3 functionality (Optional)
604	# * cdef _sub_c_impl
605	# * __deepcopy__
606	# * __invert__
607	# * Matrix windows -- only if you need strassen for that base
608	# * Other functions (list them here):
609	# * Specialized echelon form
610	########################################################################
611
612	def charpoly(self, var='x', algorithm='linbox'):
613	"""
614	INPUT:
615	var -- a variable name
616	algorithm -- 'linbox' (default)
617	'generic'
618
619	NOTE: Linbox only used on Darwin OS X right now.
620
621	EXAMPLES:
622	sage: A = matrix(ZZ,6, range(36))
623	sage: f = A.charpoly(); f
624	x^6 - 105x^5 - 630x^4
625	sage: f(A) == 0
626	True
627	sage: n=20; A = Mat(ZZ,n)(range(n^2))
628	sage: A.charpoly()
629	x^20 - 3990x^19 - 266000x^18
630	sage: A.minpoly() # optional -- os x only right now
631	x^3 - 3990x^2 - 266000x
632	"""
633	key = 'charpoly_%s_%s'%(algorithm, var)
634	x = self.fetch(key)
635	if x: return x
636
637	if UNAME != "Darwin" and algorithm == "linbox":
638	algorithm = 'generic'
639
640	if algorithm == 'linbox':
641	g = self._charpoly_linbox(var)
642	elif algorithm == 'generic':
643	g = matrix_dense.Matrix_dense.charpoly(self, var)
644	else:
645	raise ValueError, "no algorithm '%s'"%algorithm
646
647	self.cache(key, g)
648	return g
649
650	def minpoly(self, var='x', algorithm='linbox'):
651	"""
652	INPUT:
653	var -- a variable name
654	algorithm -- 'linbox' (default)
655	'generic'
656
657	NOTE: Linbox only used on Darwin OS X right now.
658
659	EXAMPLES:
660	sage: A = matrix(ZZ,6, range(36))
661	sage: A.minpoly() # optional -- os x only right now
662	x^3 - 105x^2 - 630x
663	sage: n=6; A = Mat(ZZ,n)([k^2 for k in range(n^2)])
664	sage: A.minpoly() # optional -- os x only right now
665	x^4 - 2695x^3 - 257964x^2 + 1693440*x
666	"""
667	key = 'minpoly_%s_%s'%(algorithm, var)
668	x = self.fetch(key)
669	if x: return x
670
671	if UNAME != "Darwin" and algorithm == "linbox":
672	algorithm = 'generic'
673
674	if algorithm == 'linbox':
675	g = self._minpoly_linbox(var)
676	elif algorithm == 'generic':
677	g = matrix_dense.Matrix_dense.minpoly(self, var)
678	else:
679	raise ValueError, "no algorithm '%s'"%algorithm
680
681	self.cache(key, g)
682	return g
683
684	def _minpoly_linbox(self, var='x'):
685	return self._poly_linbox(var=var, typ='minpoly')
686
687	def _charpoly_linbox(self, var='x'):
688	return self._poly_linbox(var=var, typ='charpoly')
689
690	def _poly_linbox(self, var='x', typ='minpoly'):
691	"""
692	INPUT:
693	var -- 'x'
694	typ -- 'minpoly' or 'charpoly'
695	"""
696	time = verbose('computing %s of %s x %s matrix using linbox'%(typ, self._nrows, self._ncols))
697	if self._nrows != self._ncols:
698	raise ValueError, "matrix must be square"
699	if self._nrows <= 1:
700	return matrix_dense.Matrix_dense.charpoly(self, var)
701	cdef mpz_t* poly
702	cdef size_t degree
703	if self._nrows % 4 == 0 and UNAME == "Darwin":
704	verbose("using hack to get around bug in linbox on OS X since n is divisible by 4")
705	if typ == 'minpoly':
706	_sig_on
707	linbox_integer_dense_minpoly_hacked(&poly, &degree, self._nrows, self._matrix,1)
708	_sig_off
709	else:
710	_sig_on
711	linbox_integer_dense_minpoly_hacked(&poly, &degree, self._nrows, self._matrix,0)
712	_sig_off
713	else:
714	if typ == 'minpoly':
715	_sig_on
716	linbox_integer_dense_minpoly(&poly, &degree, self._nrows, self._matrix)
717	_sig_off
718	else:
719	_sig_on
720	linbox_integer_dense_charpoly(&poly, &degree, self._nrows, self._matrix)
721	_sig_off
722
723	v = []
724	cdef Integer k
725	cdef size_t n
726	for n from 0 <= n <= degree:
727	k = PY_NEW(Integer)
728	mpz_set(k.value, poly[n])
729	mpz_clear(poly[n])
730	v.append(k)
731	linbox_integer_dense_delete_array(poly)
732	R = self._base_ring[var]
733	verbose('finished computing %s'%typ, time)
734	return R(v)
735
736
737	def height(self):
738	"""
739	Return the height of this matrix, i.e., the max absolute value
740	of the entries of the matrix.
741
742	OUTPUT:
743	A nonnegative integer.
744
745	EXAMPLE:
746	sage: a = Mat(ZZ,3)(range(9))
747	sage: a.height()
748	8
749	sage: a = Mat(ZZ,2,3)([-17,3,-389,15,-1,0]); a
750	[ -17 3 -389]
751	[ 15 -1 0]
752	sage: a.height()
753	389
754	"""
755	cdef mpz_t h
756	cdef Integer x
757
758	self.mpz_height(h)
759	x = Integer()
760	x.set_from_mpz(h)
761	mpz_clear(h)
762
763	return x
764
765	cdef int mpz_height(self, mpz_t height) except -1:
766	"""
767	Used to compute the height of this matrix.
768
769	INPUT:
770	height -- a GMP mpz_t (that has not been initialized!)
771	OUTPUT:
772	sets the value of height to the height of this matrix, i.e., the max absolute
773	value of the entries of the matrix.
774	"""
775	cdef mpz_t x, h
776	cdef Py_ssize_t i
777
778	mpz_init_set_si(h, 0)
779	mpz_init(x)
780
781	_sig_on
782
783	for i from 0 <= i < self._nrows * self._ncols:
784	mpz_abs(x, self._entries[i])
785	if mpz_cmp(h, x) < 0:
786	mpz_set(h, x)
787
788	_sig_off
789
790	mpz_init_set(height, h)
791	mpz_clear(h)
792	mpz_clear(x)
793
794	return 0 # no error occured.
795
796	def _multiply_multi_modular(left, Matrix_integer_dense right):
797
798	cdef Integer h
799	cdef mod_int *moduli
800	cdef int i, n
801
802	h = left.height() * right.height()
803	n = mm.moduli_list_c(&moduli, h.value)
804	res = []
805	for i from 0 <= i < n:
806	res.append(left._mod_int_c(moduli[i]) * right._mod_int_c(moduli[i]))
807	sage_free(moduli)
808	return left._lift_crt(res)
809
810	def _mod_int(self, modulus):
811	return self._mod_int_c(modulus)
812
813	cdef _mod_int_c(self, mod_int p):
814	cdef Py_ssize_t i, j
815	cdef Matrix_modn_dense res
816	cdef mpz_t* self_row
817	cdef mod_int* res_row
818	res = Matrix_modn_dense.__new__(Matrix_modn_dense, matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
819	for i from 0 <= i < self._nrows:
820	self_row = self._matrix[i]
821	res_row = res._matrix[i]
822	for j from 0 <= j < self._ncols:
823	res_row[j] = mpz_fdiv_ui(self_row[j], p)
824	return res
825
826	def _lift_crt(self, residues):
827
828	cdef size_t n, i, j, k
829	cdef Py_ssize_t nr, nc
830
831	n = len(residues)
832	nr = residues[0].nrows()
833	nc = residues[0].ncols()
834
835	for b in residues:
836	if not PY_TYPE_CHECK(b, Matrix_modn_dense):
837	raise TypeError, "Can only perform CRT on list of type Matrix_modn_dense."
838	cdef PyObject** res
839	res = FAST_SEQ_UNSAFE(residues)
840
841	cdef mod_int **row_list
842	row_list = <mod_int*>sage_malloc(sizeof(mod_int) * n)
843	if row_list == NULL:
844	raise MemoryError, "out of memory allocating multi-modular coefficent list"
845
846	cdef Matrix_integer_dense M
847	M = Matrix_integer_dense.__new__(Matrix_integer_dense, self.matrix_space(nr, nc), None, None, None)
848
849	_sig_on
850	for i from 0 <= i < nr:
851	for k from 0 <= k < n:
852	row_list[k] = (<Matrix_modn_dense>res[k])._matrix[i]
853	mm.mpz_crt_vec(M._matrix[i], row_list, n, nc)
854	_sig_off
855
856	sage_free(row_list)
857	return M
858
859
860	def _echelon_in_place_classical(self):
861	cdef Matrix_integer_dense E
862	E = self.echelon_form()
863
864	cdef int i
865	for i from 0 <= i < self._ncols * self._nrows:
866	mpz_set(self._entries[i], E._entries[i])
867
868	self.clear_cache()
869
870	def _echelon_strassen(self):
871	raise NotImplementedError
872
873	def echelon_form(self, algorithm="default", cutoff=0, include_zero_rows=True):
874	r"""
875	Return the echelon form of this matrix over the integers.
876
877	INPUT:
878	algorithm, cutoff -- ignored currently
879	include_zero_rows -- (default: True) if False, don't include zero rows.
880
881	EXAMPLES:
882	sage: A = MatrixSpace(ZZ,2)([1,2,3,4])
883	sage: A.echelon_form()
884	[1 0]
885	[0 2]
886
887	sage: A = MatrixSpace(ZZ,5)(range(25))
888	sage: A.echelon_form()
889	[ 5 0 -5 -10 -15]
890	[ 0 1 2 3 4]
891	[ 0 0 0 0 0]
892	[ 0 0 0 0 0]
893	[ 0 0 0 0 0]
894	"""
895	x = self.fetch('echelon_form')
896	if not x is None:
897	return x
898
899	if self._nrows == 0 or self._ncols == 0:
900	self.cache('echelon_form', self)
901	self.cache('pivots', [])
902	self.cache('rank', 0)
903	return self
904
905	cdef Py_ssize_t nr, nc, n, i
906	nr = self._nrows
907	nc = self._ncols
908
909	# The following complicated sequence of column reversals
910	# and transposes is needed since PARI's Hermite Normal Form
911	# does column operations instead of row operations.
912	n = nc
913	r = []
914	for i from 0 <= i < n:
915	r.append(n-i)
916	v = self._pari_()
917	v = v.vecextract(r) # this reverses the order of columns
918	v = v.mattranspose()
919	w = v.mathnf(1)
920
921	cdef Matrix_integer_dense H_m
922	H = convert_parimatrix(w[0])
923	if nc == 1:
924	H = [H]
925
926	# We do a 'fast' change of the above into a list of ints,
927	# since we know all entries are ints:
928	num_missing_rows = (nr*nc - len(H)) / nc
929	rank = nr - num_missing_rows
930
931	if include_zero_rows:
932	H = H + ['0'](num_missing_rowsnc)
933	H_m = self.new_matrix(nrows=nr, ncols=nc, entries=H, coerce=True)
934	else:
935	H_m = self.new_matrix(nrows=rank, ncols=nc, entries=H, coerce=True)
936
937	H_m.set_immutable()
938	H_m.cache('rank', rank)
939	self.cache('rank',rank)
940	H_m.cache('echelon_form',H_m)
941	self.cache('echelon_form',H_m)
942	return H_m
943
944	def pivots(self):
945	"""
946	Return the pivot column positions of this matrix as a list of Python integers.
947
948	This returns a list, of the position of the first nonzero entry in each row
949	of the echelon form.
950
951	OUTPUT:
952	list -- a list of Python ints
953
954	EXAMPLES:
955	sage: n = 3; A = matrix(ZZ,n,range(n^2)); A
956	[0 1 2]
957	[3 4 5]
958	[6 7 8]
959	sage: A.pivots()
960	[0, 1]
961	sage: A.echelon_form()
962	[ 3 0 -3]
963	[ 0 1 2]
964	[ 0 0 0]
965	"""
966	p = self.fetch('pivots')
967	if not p is None: return p
968
969	cdef Matrix_integer_dense E
970	E = self.echelon_form()
971
972	# Now we determine the pivots from the matrix E as quickly as we can.
973	# For each row, we find the first nonzero position in that row -- it is the pivot.
974	cdef Py_ssize_t i, j, k
975	cdef mpz_t *row
976	p = []
977	k = 0
978	for i from 0 <= i < E._nrows:
979	row = E._matrix[i] # pointer to ith row
980	for j from k <= j < E._ncols:
981	if mpz_cmp_si(row[j], 0) != 0: # nonzero position
982	p.append(j)
983	k = j+1 # so start at next position next time
984	break
985	self.cache('pivots', p)
986	return p
987
988	def elementary_divisors(self):
989	"""
990	Return the elementary divisors of self, in order.
991
992	The elementary divisors are the invariants of the finite
993	abelian group that is the cokernel of this matrix. They are
994	ordered in reverse by divisibility.
995
996	INPUT:
997	matrix
998	OUTPUT:
999	list of int's
1000
1001	EXAMPLES:
1002	sage: A = MatrixSpace(IntegerRing(), 3)(range(9))
1003	sage: A.elementary_divisors()
1004	[0, 3, 1]
1005	sage: C = MatrixSpace(ZZ,4)([3,4,5,6,7,3,8,10,14,5,6,7,2,2,10,9])
1006	sage: C.elementary_divisors()
1007	[687, 1, 1, 1]
1008
1009	SEE ALSO: smith_form
1010	"""
1011	d = self.fetch('elementary_divisors')
1012	if not d is None:
1013	return d
1014	if self._nrows == 0 or self._ncols == 0:
1015	d = []
1016	else:
1017	d = self._pari_().matsnf(0).python()
1018	self.cache('elementary_divisors', d)
1019	return d
1020
1021	def smith_form(self):
1022	"""
1023	Returns matrices S, U, and V such that S = UselfV, and S
1024	is in Smith normal form. Thus S is diagonal with diagonal
1025	entries the ordered elementary divisors of S.
1026
1027	The elementary divisors are the invariants of the finite
1028	abelian group that is the cokernel of this matrix. They are
1029	ordered in reverse by divisibility.
1030
1031	EXAMPLES:
1032	sage: A = MatrixSpace(IntegerRing(), 3)(range(9))
1033	sage: D, U, V = A.smith_form()
1034	sage: D
1035	[0 0 0]
1036	[0 3 0]
1037	[0 0 1]
1038	sage: U
1039	[-1 2 -1]
1040	[ 0 -1 1]
1041	[ 0 1 0]
1042	sage: V
1043	[ 1 4 -1]
1044	[-2 -3 1]
1045	[ 1 0 0]
1046	sage: UAV
1047	[0 0 0]
1048	[0 3 0]
1049	[0 0 1]
1050
1051	It also makes sense for nonsquare matrices:
1052
1053	sage: A = Matrix(ZZ,3,2,range(6))
1054	sage: D, U, V = A.smith_form()
1055	sage: D
1056	[0 0]
1057	[2 0]
1058	[0 1]
1059	sage: U
1060	[-1 2 -1]
1061	[ 0 -1 1]
1062	[ 0 1 0]
1063	sage: V
1064	[ 3 -1]
1065	[-2 1]
1066	sage: U * A * V
1067	[0 0]
1068	[2 0]
1069	[0 1]
1070
1071	SEE ALSO: elementary_divisors
1072	"""
1073	v = self._pari_().matsnf(1).python()
1074	D = self.matrix_space()(v[2])
1075	U = self.matrix_space(ncols = self._nrows)(v[0])
1076	V = self.matrix_space(nrows = self._ncols)(v[1])
1077	return D, U, V
1078
1079	def frobenius(self,flag=0, var='x'):
1080	"""
1081	Return the Frobenius form (rational canonical form) of this matrix.
1082
1083	INPUT:
1084	flag --an integer:
1085	0 -- (default) return the Frobenius form of this matrix.
1086	1 -- return only the elementary divisor polynomials, as
1087	polynomials in var.
1088	2 -- return a two-components vector [F,B] where F is the
1089	Frobenius form and B is the basis change so that $M=B^{-1}FB$.
1090	var -- a string (default: 'x')
1091
1092	INPUT:
1093	flag -- 0 (default), 1 or 2 as described above
1094
1095	ALGORITHM: uses pari's matfrobenius()
1096
1097	EXAMPLE:
1098	sage: A = MatrixSpace(ZZ, 3)(range(9))
1099	sage: A.frobenius(0)
1100	[ 0 0 0]
1101	[ 1 0 18]
1102	[ 0 1 12]
1103	sage: A.frobenius(1)
1104	[x^3 - 12x^2 - 18x]
1105	sage: A.frobenius(1, var='y')
1106	[y^3 - 12y^2 - 18y]
1107	sage: A.frobenius(2)
1108	([ 0 0 0]
1109	[ 1 0 18]
1110	[ 0 1 12],
1111	[ -1 2 -1]
1112	[ 0 23/15 -14/15]
1113	[ 0 -2/15 1/15])
1114
1115	AUTHOR:
1116	-- 2006-04-02: Martin Albrecht
1117
1118	TODO:
1119	-- move this to work for more general matrices than just over Z.
1120	This will require fixing how PARI polynomials are coerced
1121	to SAGE polynomials.
1122	"""
1123	if not self.is_square():
1124	raise ArithmeticError, "frobenius matrix of non-square matrix not defined."
1125
1126	v = self._pari_().matfrobenius(flag)
1127	if flag==0:
1128	return self.matrix_space()(v.python())
1129	elif flag==1:
1130	r = PolynomialRing(self.base_ring(), names=var)
1131	retr = []
1132	for f in v:
1133	retr.append(eval(str(f).replace("^","**"), {'x':r.gen()}, r.gens_dict()))
1134	return retr
1135	elif flag==2:
1136	F = matrix_space.MatrixSpace(QQ, self.nrows())(v[0].python())
1137	B = matrix_space.MatrixSpace(QQ, self.nrows())(v[1].python())
1138	return F, B
1139
1140	def kernel(self, LLL=False):
1141	r"""
1142	Return the kernel of this matrix, as a module over the integers.
1143
1144	INPUT:
1145	LLL -- bool (default: False); if True the basis is an LLL
1146	reduced basis; otherwise, it is an echelon basis.
1147
1148	EXAMPLES:
1149	sage: M = MatrixSpace(ZZ,4,2)(range(8))
1150	sage: M.kernel()
1151	Free module of degree 4 and rank 2 over Integer Ring
1152	Echelon basis matrix:
1153	[ 1 0 -3 2]
1154	[ 0 1 -2 1]
1155	"""
1156	if self._nrows == 0: # from a 0 space
1157	M = sage.modules.free_module.FreeModule(ZZ, self._nrows)
1158	return M.zero_subspace()
1159
1160	elif self._ncols == 0: # to a 0 space
1161	return sage.modules.free_module.FreeModule(ZZ, self._nrows)
1162
1163	A = self._pari_().mattranspose()
1164	B = A.matkerint()
1165	n = self._nrows
1166	M = sage.modules.free_module.FreeModule(ZZ, n)
1167
1168	if B.ncols() == 0:
1169	return M.zero_submodule()
1170
1171	# Now B is a basis for the LLL-reduced integer kernel as a
1172	# PARI object. The basis vectors or B[0], ..., B[n-1],
1173	# where n is the dimension of the kernel.
1174	X = []
1175	for b in B:
1176	tmp = []
1177	for x in b:
1178	tmp.append(ZZ(x))
1179	X.append(M(tmp))
1180
1181	if LLL:
1182	return M.span_of_basis(X)
1183	else:
1184	return M.span(X)
1185
1186	def _adjoint(self):
1187	"""
1188	Return the adjoint of this matrix.
1189
1190	Assumes self is a square matrix (checked in adjoint).
1191	"""
1192	return self.parent()(self._pari_().matadjoint().python())
1193
1194	def _ntl_(self):
1195	r"""
1196	ntl.mat_ZZ representation of self.
1197
1198	\note{NTL only knows dense matrices, so if you provide a
1199	sparse matrix NTL will allocate memory for every zero entry.}
1200	"""
1201	return mat_ZZ(self._nrows,self._ncols, self.list())
1202
1203
1204	####################################################################################
1205	# LLL
1206	####################################################################################
1207	def lllgram(self):
1208	"""
1209	LLL reduction of the lattice whose gram matrix is self.
1210
1211	INPUT:
1212	M -- gram matrix of a definite quadratic form
1213
1214	OUTPUT:
1215	U -- unimodular transformation matrix such that
1216
1217	U.transpose() * M * U
1218
1219	is LLL-reduced
1220
1221	ALGORITHM:
1222	Use PARI
1223
1224	EXAMPLES:
1225	sage: M = Matrix(ZZ, 2, 2, [5,3,3,2]) ; M
1226	[5 3]
1227	[3 2]
1228	sage: U = M.lllgram(); U
1229	[-1 1]
1230	[ 1 -2]
1231	sage: U.transpose() * M * U
1232	[1 0]
1233	[0 1]
1234
1235	Semidefinite and indefinite forms raise a ValueError:
1236
1237	sage: Matrix(ZZ,2,2,[2,6,6,3]).lllgram()
1238	Traceback (most recent call last):
1239	...
1240	ValueError: not a definite matrix
1241	sage: Matrix(ZZ,2,2,[1,0,0,-1]).lllgram()
1242	Traceback (most recent call last):
1243	...
1244	ValueError: not a definite matrix
1245
1246	BUGS:
1247	should work for semidefinite forms (PARI is ok)
1248	"""
1249	if self._nrows != self._ncols:
1250	raise ArithmeticError, "matrix must be square"
1251
1252	n = self.nrows()
1253	# maybe should be /unimodular/ matrices ?
1254	P = self._pari_()
1255	try:
1256	U = P.lllgramint()
1257	except (RuntimeError, ArithmeticError), msg:
1258	raise ValueError, "not a definite matrix"
1259	MS = matrix_space.MatrixSpace(ZZ,n)
1260	U = MS(U.python())
1261	# Fix last column so that det = +1
1262	if U.det() == -1:
1263	for i in range(n):
1264	U[i,n-1] = - U[i,n-1]
1265	return U
1266
1267	def prod_of_row_sums(self, cols):
1268	cdef Py_ssize_t c, row
1269	cdef mpz_t s, pr
1270	mpz_init(s)
1271	mpz_init(pr)
1272
1273	mpz_set_si(pr, 1)
1274	for row from 0 <= row < self._nrows:
1275	tmp = []
1276	mpz_set_si(s, 0)
1277	for c in cols:
1278	if c<0 or c >= self._ncols:
1279	raise IndexError, "matrix column index out of range"
1280	mpz_add(s, s, self._matrix[row][c])
1281	mpz_mul(pr, pr, s)
1282	cdef Integer z
1283	z = PY_NEW(Integer)
1284	mpz_set(z.value, pr)
1285	mpz_clear(s)
1286	mpz_clear(pr)
1287	return z
1288
1289	def _linbox_sparse(self):
1290	cdef Py_ssize_t i, j
1291	s = '%s %s M\n'%(self._nrows, self._ncols)
1292	for i from 0 <= i < self._nrows:
1293	for j from 0 <= j < self._ncols:
1294	if mpz_cmp_si(self._matrix[i][j], 0):
1295	s += '%s %s %s\n'%(i+1,j+1,self.get_unsafe(i,j))
1296	s += '0 0 0\n'
1297	return s
1298
1299	def _linbox_dense(self):
1300	cdef Py_ssize_t i, j
1301	s = '%s %s x'%(self._nrows, self._ncols)
1302	for i from 0 <= i < self._nrows:
1303	for j from 0 <= j < self._ncols:
1304	s += ' %s'%self.get_unsafe(i,j)
1305	return s
1306
1307	def rational_reconstruction(self, N):
1308	"""
1309	Use rational reconstruction to lift self to a matrix over the
1310	rational numbers (if possible), where we view self as a matrix
1311	modulo N.
1312	"""
1313	import misc
1314	return misc.matrix_integer_dense_rational_reconstruction(self, N)
1315
1316	def randomize(self, density=1, x=None, y=None):
1317	"""
1318	Randomize density proportion of the entries of this matrix,
1319	leaving the rest unchanged.
1320
1321	The randomized entries of this matrix to be between x and y
1322	and have density 1.
1323	"""
1324	self.check_mutability()
1325	self.clear_cache()
1326
1327	cdef int _min, _max
1328	if y is None:
1329	if x is None:
1330	min = -2
1331	max = 3
1332	else:
1333	min = 0
1334	max = x
1335	else:
1336	min = x
1337	max = y
1338
1339	density = float(density)
1340
1341	cdef int min_is_zero
1342	min_is_nonzero = (min != 0)
1343
1344	cdef Integer n_max, n_min, n_width
1345	n_max = Integer(max)
1346	n_min = Integer(min)
1347	n_width = n_max - n_min
1348
1349	cdef Py_ssize_t i, j, k, nc, num_per_row
1350	global state
1351
1352	_sig_on
1353	if density == 1:
1354	for i from 0 <= i < self._nrows*self._ncols:
1355	mpz_urandomm(self._entries[i], state, n_width.value)
1356	if min_is_nonzero:
1357	mpz_add(self._entries[i], self._entries[i], n_min.value)
1358	else:
1359	nc = self._ncols
1360	num_per_row = int(density * nc)
1361	for i from 0 <= i < self._nrows:
1362	for j from 0 <= j < num_per_row:
1363	k = random()%nc
1364	mpz_urandomm(self._matrix[i][k], state, n_width.value)
1365	if min_is_nonzero:
1366	mpz_add(self._matrix[i][k], self._matrix[i][j], n_min.value)
1367	_sig_off
1368
1369
1370	###############################################################
1371
1372	###########################################
1373	# Helper code for Echelon form algorithm.
1374	###########################################
1375	def _parimatrix_to_strlist(A):
1376	s = str(A)
1377	s = s.replace('Mat(','').replace(')','')
1378	s = s.replace(';',',').replace(' ','')
1379	s = s.replace(",", "','")
1380	s = s.replace("[", "['")
1381	s = s.replace("]", "']")
1382	return eval(s)
1383
1384	def _parimatrix_to_reversed_strlist(A):
1385	s = str(A)
1386	if s.find('Mat') != -1:
1387	return _parimatrix_to_strlist(A)
1388	s = s.replace('[','').replace(']','').replace(' ','')
1389	v = s.split(';')
1390	v.reverse()
1391	s = "['" + (','.join(v)) + "']"
1392	s = s.replace(",", "','")
1393	return eval(s)
1394
1395	def convert_parimatrix(z):
1396	n = z.ncols();
1397	r = []
1398	for i from 0 <= i < n:
1399	r.append(n-i)
1400	z = z.vecextract(r)
1401	z = z.mattranspose()
1402	n = z.ncols();
1403	r = []
1404	for i from 0 <= i < n:
1405	r.append(n-i)
1406	z = z.vecextract(r)
1407	return _parimatrix_to_strlist(z)
1408
1409
1410	##########################################################
1411	# Setup the c-library and GMP random number generators.
1412	# seed it when module is loaded.
1413	from random import randrange
1414	cdef extern from "stdlib.h":
1415	long random()
1416	void srandom(unsigned int seed)
1417	k = randrange(0,2**32)
1418	srandom(k)
1419
1420	cdef gmp_randstate_t state
1421	gmp_randinit_mt(state)
1422	gmp_randseed_ui(state,k)
1423

Note: See TracBrowser for help on using the repository browser.

Download in other formats: