Context Navigation

matrix_rational_sparse.pyx @ 3092:227119d39b73

Visit:

Revision 3092:227119d39b73, 47.4 KB checked in by William Stein <wstein@…>, 6 years ago (diff)
matrix coding sprint -- snapshot.

Rev	Line
[1775]	1	r"""nodoctest
[1687]	2	Sparse matrices over $\Q$
[6]	3
[3092]	4	This is an optimzed implementation of sparse matrices over
	5	the rational numbers.
[771]	6
[6]	7	AUTHOR:
	8	-- William Stein (2004): first version
[11]	9	-- William Stein (2006-02-12): added set_row_to_multiple_of_row
[1687]	10	-- William Stein (2006-03-04): added multimodular echelon, __reduce__, etc.
[3092]	11	-- William Stein (2007-02-20): update to new SAGE matrix class structure.
[1687]	12
	13	EXAMPLES:
	14	sage: M = MatrixSpace(QQ, 2, 3, sparse=True)
	15	sage: A = M([1,2,3, 1,1,1])
	16	sage: A
	17	[1 2 3]
	18	[1 1 1]
	19	sage: A.echelon_form()
	20	[ 1 0 -1]
	21	[ 0 1 2]
	22
	23	sage: M = MatrixSpace(QQ, 1000,1000, sparse=True)
	24	sage: A = M(0)
	25	sage: A[1,1] = 5
	26
	27
	28	sage: from sage.matrix.sparse_matrix import SparseMatrix
	29	sage: x = SparseMatrix(QQ, 5,10)
	30	sage: x.randomize(5)
	31	sage: x.echelon_form() # random output
	32	[
	33	1, 0, 0, 0, 0, 0, 0, 0, 0, 10/29,
	34	0, 1, 0, 0, 0, 0, 0, 0, 0, -4/29,
	35	0, 0, 1, 0, 0, 0, 0, 0, 0, -12/29,
	36	0, 0, 0, 0, 0, 0, 1, 0, 0, 24/29,
	37	0, 0, 0, 0, 0, 0, 0, 0, 1, 4/29
	38	]
	39
[0]	40	"""
	41
[10]	42	#############################################################################
[3092]	43	# Copyright (C) 2004, 2007 William Stein <wstein@gmail.com>
[0]	44	# Distributed under the terms of the GNU General Public License (GPL)
	45	# The full text of the GPL is available at:
	46	# http://www.gnu.org/licenses/
[10]	47	#############################################################################
[0]	48
	49	import random
	50
[11]	51	import sage.misc.all
[1593]	52	import sage.matrix.matrix_rational_sparse
[11]	53	import sage.rings.integer_ring
[823]	54	import sage.rings.finite_field
[11]	55	import sage.rings.arith
[0]	56
[1319]	57	cimport sage.rings.rational
	58	import sage.rings.rational
[11]	59
[1319]	60	cimport sage.rings.integer
	61	import sage.rings.integer
[11]	62
[1319]	63	cimport sage.ext.arith
	64	import sage.ext.arith
	65	cdef sage.ext.arith.arith_int ai
	66	ai = sage.ext.arith.arith_int()
[0]	67
[1664]	68	cimport matrix_field_sparse
	69	import matrix_field_sparse
[0]	70
[1593]	71	cimport matrix_modn_sparse
	72	import matrix_modn_sparse
	73
	74	cimport matrix_rational_dense
	75	import matrix_rational_dense
	76
[1319]	77	include "../ext/gmp.pxi"
[3092]	78	include '../ext/interrupt.pxi'
[0]	79
[10]	80	START_PRIME = 20011 # used for multi-modular algorithms
[0]	81
	82	cdef class Vector_mpq
	83
	84	cdef void Vector_mpq_rescale(Vector_mpq w, mpq_t x):
	85	scale_mpq_vector(&w.v, x)
	86
	87	cdef class Vector_mpq:
	88	"""
	89	Vector_mpq -- a sparse vector of GMP rationals. This is a Python
	90	extension type that wraps the C implementation of sparse vectors
	91	modulo a small prime.
	92	"""
	93	cdef mpq_vector v
	94
	95	def __init__(self, int degree, int num_nonzero=0, entries=[], sort=True):
	96	cdef int i
	97	init_mpq_vector(&self.v, degree, num_nonzero)
	98	if entries != []:
	99	if len(entries) != num_nonzero:
	100	raise ValueError, "length of entries (=%s) must equal num_nonzero (=%s)"%(len(entries), num_nonzero)
	101	if sort:
	102	entries = list(entries) # copy so as not to modify passed in list
	103	entries.sort()
	104	for i from 0 <= i < num_nonzero:
	105	s = str(entries[i][1])
	106	mpq_set_str(self.v.entries[i], s, 0)
	107	self.v.positions[i] = entries[i][0]
	108
	109	def __dealloc__(self):
	110	clear_mpq_vector(&self.v)
	111
	112	def __getitem__(self, int n):
	113	cdef mpq_t x
[1319]	114	cdef sage.rings.rational.Rational a
[0]	115	mpq_init(x)
	116	mpq_vector_get_entry(&x, &self.v, n)
[1319]	117	a = sage.rings.rational.Rational()
[0]	118	a.set_from_mpq(x)
	119	mpq_clear(x)
	120	return a
	121
	122	def cmp(self, Vector_mpq other):
	123	return mpq_vector_cmp(&self.v, &other.v)
	124
	125	def __richcmp__(Vector_mpq self, x, int op):
	126	if not isinstance(x, Vector_mpq):
	127	return -1
	128	cdef int n
	129	n = self.cmp(x)
	130	if op == 0:
	131	return bool(n < 0)
	132	elif op == 1:
	133	return bool(n <= 0)
	134	elif op == 2:
	135	return bool(n == 0)
	136	elif op == 3:
	137	return bool(n != 0)
	138	elif op == 4:
	139	return bool(n > 0)
	140	elif op == 5:
	141	return bool(n >= 0)
	142
	143	def __setitem__(self, int n, x):
	144	cdef object s
	145	s = str(x)
	146	mpq_vector_set_entry_str(&self.v, n, s)
	147
	148	def __repr__(self):
	149	return str(list(self))
	150
	151	def degree(self):
	152	return self.v.degree
	153
	154	def num_nonzero(self):
	155	return self.v.num_nonzero
	156
	157	def list(self):
	158	return mpq_vector_to_list(&self.v)
	159
	160	cdef void rescale(self, mpq_t x):
	161	scale_mpq_vector(&self.v, x)
	162
	163	def __add__(Vector_mpq self, Vector_mpq other):
	164	cdef mpq_vector z1, *z2
	165	cdef Vector_mpq w
	166	cdef mpq_t ONE
	167	mpq_init(ONE)
	168	mpq_set_si(ONE,1,1)
	169
	170	add_mpq_vector_init(&z1, &self.v, &other.v, ONE)
	171	mpq_clear(ONE)
	172	w = Vector_mpq(self.v.degree)
	173	z2 = &(w.v)
	174	clear_mpq_vector(z2) # free memory wasted on allocated w
	175	z2.entries = z1.entries
	176	z2.positions = z1.positions
	177	z2.num_nonzero = z1.num_nonzero
	178	# at this point we do not free z1, since it is referenced by w.
	179	return w
	180
	181	def __sub__(Vector_mpq self, Vector_mpq other):
	182	return self + other*(-1)
	183
	184	def copy(self):
	185	cdef int i
	186	cdef Vector_mpq w
	187	w = Vector_mpq(self.v.degree, self.v.num_nonzero)
	188	for i from 0 <= i < self.v.num_nonzero:
	189	mpq_set(w.v.entries[i], self.v.entries[i])
	190	w.v.positions[i] = self.v.positions[i]
	191	return w
	192
	193	def __mul__(x, y):
	194	if isinstance(x, Vector_mpq):
	195	self = x
	196	other = y
	197	elif isinstance(y, Vector_mpq):
	198	self = y
	199	other = x
	200	else:
	201	raise TypeError, "Invalid types."
	202	cdef object s, z
	203	cdef mpq_t t
	204	z = self.copy()
	205	mpq_init(t)
	206	s = str(other)
	207	mpq_set_str(t, s, 0)
	208	Vector_mpq_rescale(z, t)
	209	mpq_clear(t)
	210	return z
	211
	212	def randomize(self, int sparcity, bound=3):
	213	"""
	214	randomize(self, int sparcity, exact=False):
	215
	216	The sparcity is a bound on the number of nonzeros per row.
	217	"""
	218	cdef int i
	219	for i from 0 <= i < sparcity:
	220	self[random.randrange(self.v.degree)] = random.randrange(1,bound)
	221
	222
	223	#############################################################
	224	#
	225	# Sparse Matrix over mpq_t (the GMP rationals)
	226	#
	227	#############################################################
[2732]	228	cdef class Matrix_rational_sparse(matrix_sparse.Matrix_sparse):
[0]	229
[11]	230	def __new__(self, int nrows, int ncols, object entries=[], init=True, coerce=False):
[0]	231	# allocate memory
	232	cdef int i
[1793]	233	self.rows = <mpq_vector> sage_malloc(nrowssizeof(mpq_vector))
[11]	234	self.is_init = init
	235	if self.is_init:
	236	self.is_init = True
	237	for i from 0 <= i < nrows:
	238	init_mpq_vector(&self.rows[i], ncols, 0)
[0]	239
	240	def __dealloc__(self):
[11]	241	if not self.is_init:
	242	return
[0]	243	cdef int i
	244	for i from 0 <= i < self.nr:
	245	clear_mpq_vector(&self.rows[i])
	246
[11]	247	def __init__(self, int nrows, int ncols, object entries=[], init=True, coerce=False):
[0]	248	"""
[10]	249	INPUT:
	250	nrows -- number of rows
	251	ncols -- number of columns
	252	entries -- list of triples (i,j,x), where 0 <= i < nrows,
	253	0 <= j < ncols, and x is a rational number.
	254	Then the i,j entry of the matrix is set to x.
	255	It is OK for some x to be zero.
	256	or, list of all entries of the matrix (dense representation).
[11]	257	init -- bool (default: True); if False, don't allocate anything (for
	258	internal use only!)
	259	coerce -- bool (default: False), if True, entries might not be of type Rational.
[0]	260	"""
	261	cdef object s
[11]	262	cdef int ii, jj, k, n
[1319]	263	cdef sage.rings.rational.Rational z
[10]	264
[0]	265	self.nr = nrows
	266	self.nc = ncols
[11]	267	self.__pivots = None
	268
	269	if isinstance(entries, str):
	270	v = entries.split('\n')
[10]	271	for ii from 0 <= ii < nrows:
[11]	272	w = v[ii].split()
	273	n = int(w[0])
	274	init_mpq_vector(&self.rows[ii], ncols, n)
	275	for jj from 0 <= jj < n:
	276	self.rows[ii].positions[jj] = int(w[jj+1])
	277	t = w[jj+n+1]
	278	mpq_set_str(self.rows[ii].entries[jj], t, 32)
	279	self.is_init = True
	280	return
	281
	282	if not init:
	283	return
	284
	285	if not coerce:
	286	if isinstance(entries, list):
	287	if len(entries) == 0:
	288	return
	289	if not isinstance(entries[0], tuple):
	290	# dense input representation
	291	k = 0
	292	for ii from 0 <= ii < nrows:
	293	_sig_check
	294	for jj from 0 <= jj < ncols:
	295	z = entries[k]
	296	if mpq_sgn(z.value): # if z is nonzero
	297	mpq_vector_set_entry(&self.rows[ii], jj, z.value)
	298	k = k + 1
	299	else:
	300	# sparse input rep
	301	for i, j, x in entries:
	302	z = x
	303	if mpq_sgn(z.value): # if z is nonzero
	304	mpq_vector_set_entry(&self.rows[i], j, z.value)
	305	return
	306
	307	if isinstance(entries, dict):
	308	for (i,j), x in entries.iteritems():
	309	z = x
	310	mpq_vector_set_entry(&self.rows[i], j, z.value)
	311	return
	312
[10]	313	else:
[11]	314	# copying the above code is very ugly. but this
	315	# function must be very fast...
	316	if isinstance(entries, list):
	317	if len(entries) == 0:
	318	return
	319	if not isinstance(entries[0], tuple):
	320	# dense input representation
	321	k = 0
	322	for ii from 0 <= ii < nrows:
	323	_sig_check
	324	for jj from 0 <= jj < ncols:
[1319]	325	z = sage.rings.rational.Rational(entries[k])
[11]	326	if mpq_sgn(z.value): # if z is nonzero
	327	mpq_vector_set_entry(&self.rows[ii], jj, z.value)
	328	k = k + 1
	329	else:
	330	# one sparse rep
	331	for i, j, x in entries:
[1319]	332	z = sage.rings.rational.Rational(x)
[11]	333	if mpq_sgn(z.value): # if z is nonzero
	334	mpq_vector_set_entry(&self.rows[i], j, z.value)
	335	return
	336
	337	if isinstance(entries, dict):
	338	for (i,j), x in entries.iteritems():
[1319]	339	z = sage.rings.rational.Rational(x)
[10]	340	mpq_vector_set_entry(&self.rows[i], j, z.value)
[11]	341	return
	342
	343	# Now assume entries is a rationl number and matrix should be scalar
	344	if nrows == ncols:
[1319]	345	x = sage.rings.rational.Rational(entries)
[11]	346	for ii from 0 <= ii < nrows:
	347	self[ii,ii] = x
	348	return
	349
	350	raise TypeError, "no way to make matrix from entries (=%s)"%entries
	351
	352
[1593]	353	def __cmp__(self, Matrix_rational_sparse other):
[11]	354	if self.nr != other.nr:
	355	return -1
	356	if self.nc != other.nc:
	357	return -1
	358	cdef int i, j, c
	359	for i from 0 <= i < self.nr:
	360	if self.rows[i].num_nonzero != other.rows[i].num_nonzero:
	361	return 1
	362	for j from 0 <= j < self.rows[i].num_nonzero:
	363	c = mpq_cmp(self.rows[i].entries[j], other.rows[i].entries[j])
	364	if c:
	365	return c
	366	return 0
	367
	368	def __reduce__(self):
	369	"""
	370	EXAMPLES:
[1663]	371	sage: A = MatrixSpace(QQ, 2, sparse=True)([1,2,-1/5,19/397])
[11]	372	sage: loads(dumps(A)) == A
	373	True
	374	"""
	375	# Format of reduced serialized representation of a sparse matrix
	376	# number_nonzer_entries_in_row nonzero_positions nonzero_entries_in_base32[newline]
	377	#
	378	# 3 1 4 7 3/4 7/83 39/45
	379	# 2 1 193 -17 38
	380	#
	381
	382	cdef char s, t, *tmp
	383	cdef int m, n, ln, i, j, z, len_so_far
	384
	385	n = self.nr * 200 + 30
[1793]	386	s = <char> sage_malloc(n sizeof(char))
[11]	387	len_so_far = 0
	388	t = s
	389	s[0] = <char>0 # make s a null-terminated string
	390	for i from 0 <= i < self.nr:
	391	ln = self.rows[i].num_nonzero
	392	if len_so_far + 20*ln >= n:
	393	# copy to new string with double the size
	394	n = 2n + 20ln
[1793]	395	tmp = <char> sage_malloc(n sizeof(char))
[11]	396	strcpy(tmp, s)
[1793]	397	sage_free(s)
[11]	398	s = tmp
	399	t = s + len_so_far
	400	#endif
	401	z = sprintf(t, '%d ', ln)
	402	t = t + z
	403	len_so_far = len_so_far + z
	404	for j from 0 <= j < ln:
	405	z = sprintf(t, '%d ', self.rows[i].positions[j])
	406	t = t + z
	407	len_so_far = len_so_far + z
	408
	409	for j from 0 <= j < ln:
	410	m = mpz_sizeinbase(mpq_numref(self.rows[i].entries[j]), 32) + \
	411	mpz_sizeinbase(mpq_denref(self.rows[i].entries[j]), 32) + 3
	412	if len_so_far + m >= n:
	413	# copy to new string with double the size
	414	n = 2*n + m + 1
[1793]	415	tmp = <char*> sage_malloc(n)
[11]	416	strcpy(tmp, s)
[1793]	417	sage_free(s)
[11]	418	s = tmp
	419	t = s + len_so_far
	420	mpq_get_str(t, 32, self.rows[i].entries[j])
	421	m = strlen(t)
	422	len_so_far = len_so_far + m + 1
	423	t = t + m
	424	if j <= ln-1:
	425	t[0] = <char>32
	426	t[1] = <char>0
	427	t = t + 1
	428
	429	z = sprintf(t, '\n')
	430	t = t + z
	431	len_so_far = len_so_far + z
	432	# end for
	433
	434	entries = str(s)[:-1]
[1793]	435	sage_free(s)
[1663]	436	return make_sparse_rational_matrix, (self.nr, self.nc, entries)
[11]	437
[10]	438
	439	def copy(self):
[11]	440	"""
	441	Return a copy of this matrix.
[0]	442
[11]	443	EXAMPLES:
	444	sage: A = MatrixSpace(QQ, 2, sparse=True)([1,2,1/3,-4/39])._sparse_matrix_mpq_()
	445	sage: A.copy()
	446	[
	447	1, 2,
	448	1/3, -4/39
	449	]
	450	"""
[1593]	451	cdef Matrix_rational_sparse A
	452	A = Matrix_rational_sparse(self.nr, self.nc, init=False)
[11]	453	cdef int i, j, k
	454	for i from 0 <= i < self.nr:
	455	# _sig_check # not worth it since whole function is so fast?
	456	k = self.rows[i].num_nonzero
	457	if init_mpq_vector(&A.rows[i], self.nc, k) == -1:
	458	raise MemoryError, "Error allocating memory"
	459	for j from 0 <= j < k:
	460	mpq_set(A.rows[i].entries[j], self.rows[i].entries[j])
	461	A.rows[i].positions[j] = self.rows[i].positions[j]
	462	A.is_init = True
	463	return A
	464
[689]	465	def matrix_from_rows(self, rows):
[11]	466	"""
[689]	467	Return the matrix obtained from the rows of self given by the
	468	input list of integers.
	469
[11]	470	INPUT:
	471	rows -- list of integers
	472
[689]	473	EXAMPLES:
	474	sage: M = MatrixSpace(QQ,2,2,sparse=True)(range(1,5))
	475	sage: m = M._Matrix_sparse_rational__matrix
	476	sage: s = m.matrix_from_rows([1,1,1])
	477	sage: s
	478	[
	479	3, 4,
	480	3, 4,
	481	3, 4
	482	]
[11]	483	"""
	484	cdef int i, j, k, nr
	485	nr = len(rows)
	486
[1593]	487	cdef Matrix_rational_sparse A
	488	A = Matrix_rational_sparse(nr, self.nc, init=False)
[11]	489
	490	for i from 0 <= i < nr:
	491	j = rows[i]
	492	init_mpq_vector(&A.rows[i], self.nc, self.rows[j].num_nonzero)
	493	for k from 0 <= k < self.rows[j].num_nonzero:
	494	mpq_set(A.rows[i].entries[k], self.rows[j].entries[k])
[689]	495	A.rows[i].positions[k] = self.rows[j].positions[k]
[11]	496
	497	A.is_init = True
	498
	499	return A
	500
	501	def dense_matrix(self):
	502	"""
	503	Return corresponding dense matrix.
	504	"""
[1593]	505	cdef matrix_rational_dense.Matrix_rational_dense A
	506	A = matrix_rational_dense.Matrix_rational_dense(self.nr, self.nc)
[11]	507	# now A is the zero matrix. We fill in the entries.
	508	cdef int i, j, k
	509	for i from 0 <= i < self.nr:
	510	k = self.rows[i].num_nonzero
	511	for j from 0 <= j < k:
	512	# now set A[i,self.rows[i].positions[j]] equal to self.rows[i].entries[j]
[1593]	513	mpq_set(A._matrix[i][self.rows[i].positions[j]], self.rows[i].entries[j])
[11]	514	return A
	515
[0]	516	def linear_combination_of_rows(self, Vector_mpq v):
	517	if self.nr != v.degree():
	518	raise ArithmeticError, "Incompatible vector * matrix multiply."
	519	cdef mpq_vector w, sum, sum2
	520	cdef int i, r, nr
	521	cdef Vector_mpq ans
	522	nr = self.nr
	523	w = v.v
	524	init_mpq_vector(&sum, self.nc, 0)
	525	_sig_on
	526	for i from 0 <= i < w.num_nonzero:
	527	r = w.positions[i]
	528	add_mpq_vector_init(&sum2, &sum, &self.rows[r], w.entries[i])
	529	# Now sum2 is initialized and equals sum + w[i]*self.rows[i]
	530	# We want sum to equal this.
	531	clear_mpq_vector(&sum)
	532	sum = sum2
	533	_sig_off
	534	# Now sum is a sparse C-vector that gives the linear combination of rows.
	535	# Convert to a Vector_mpq and return.
	536	ans = Vector_mpq(nr)
	537	clear_mpq_vector(&ans.v)
	538	ans.v = sum
	539	return ans
[6]	540
[1319]	541	def set_row_to_multiple_of_row(self, int row_to, int row_from,
	542	sage.rings.rational.Rational multiple):
[6]	543	"""
	544	Set row row_to equal to multiple times row row_from.
	545
	546	EXAMPLES:
[1593]	547	sage: m = Matrix_rational_sparse(3,3,entries=[(1,1,10/3)])
[6]	548	sage: m
	549	[
	550	0, 0, 0,
	551	0, 10/3, 0,
	552	0, 0, 0
	553	]
	554	sage: m.set_row_to_multiple_of_row(0, 1, 6/1) # third argument must be a rational!
	555	sage: m
	556	[
	557	0, 20, 0,
	558	0, 10/3, 0,
	559	0, 0, 0
	560	]
	561	sage: m.set_row_to_multiple_of_row(2,1,-10/1)
	562	sage: m
	563	[
	564	0, 20, 0,
	565	0, 10/3, 0,
	566	0, -100/3, 0
	567	]
[1593]	568	sage: m = Matrix_rational_sparse(3,3,entries=[(1,1,10/3)])
[11]	569	sage: m.set_row_to_multiple_of_row(1,1,-10/1); m
	570	[
	571	0, 0, 0,
	572	0, -100/3, 0,
	573	0, 0, 0
	574	]
[6]	575	sage: m.set_row_to_multiple_of_row(-1, 1, 6/1)
	576	Traceback (most recent call last):
	577	...
	578	IndexError: row_to is -1 but must be >= 0 and < 3
	579	sage: m.set_row_to_multiple_of_row(0, 3, 6/1)
	580	Traceback (most recent call last):
	581	...
	582	IndexError: row_from is 3 but must be >= 0 and < 3
	583	"""
	584	# A sparse matrix is an array of pointers to mpq_vector's.
	585	# 1. Delete the vector in position row_to
	586	# 2. Initialize a new one in its place.
	587	# 3. Fill in the entries with appropriate multiples of the entries in row_from.
[11]	588	cdef int i
[6]	589
	590	if row_from < 0 or row_from >= self.nr:
	591	raise IndexError, "row_from is %s but must be >= 0 and < %s"%(row_from, self.nr)
	592	if row_to < 0 or row_to >= self.nr:
	593	raise IndexError, "row_to is %s but must be >= 0 and < %s"%(row_to, self.nr)
[11]	594
	595	if row_from == row_to:
	596	scale_mpq_vector(&self.rows[row_from], multiple.value)
	597	return
[6]	598
	599	clear_mpq_vector(&self.rows[row_to])
[11]	600	init_mpq_vector(&self.rows[row_to], self.nc, self.rows[row_from].num_nonzero)
	601
[6]	602	for i from 0 <= i < self.rows[row_from].num_nonzero:
[11]	603	mpq_mul(self.rows[row_to].entries[i], multiple.value, self.rows[row_from].entries[i])
	604	self.rows[row_to].positions[i] = self.rows[row_from].positions[i]
	605
	606
[1593]	607	def set_row_to_negative_of_row_of_A_using_subset_of_columns(self, int i, Matrix_rational_sparse A, int r, cols):
[11]	608	# the ints cols are assumed sorted.
	609	# this function exists just because it is useful for modular symbols presentations.
	610	cdef int l
	611	cdef mpq_t x
	612	mpq_init(x)
	613	l = 0
	614	for k in cols:
	615	mpq_vector_get_entry(&x, &A.rows[r], k)
	616	if mpz_cmp_si(mpq_numref(x), 0): # x is nonzero
	617	mpz_mul_si(mpq_numref(x), mpq_numref(x), -1)
	618	mpq_vector_set_entry(&self.rows[i], l, x)
	619	l = l + 1
	620	mpq_clear(x)
[0]	621
	622	def parent(self):
	623	import sage.matrix.matrix_space
	624	import sage.rings.rings
	625	return sage.matrix.matrix_space.MatrixSpace(
	626	sage.rings.rings.RationalField(),self.nr,self.nc)
	627
	628	def pivots(self):
[11]	629	if self.__pivots is None:
	630	raise NotImplementedError
[0]	631	return self.__pivots
	632
	633	def row_to_dict(self, int i):
	634	"""
	635	Return an associative arrow of pairs
	636	n:x
	637	where the keys n run through the nonzero positions of the row,
	638	and the x are nonzero and of type Integer.
	639	"""
	640	cdef int j, n
[1319]	641	cdef sage.rings.rational.Rational x
[0]	642	cdef object entries
	643	if i < 0 or i >= self.nr: raise IndexError
	644	X = {}
	645	for j from 0 <= j < self.rows[i].num_nonzero:
	646	n = self.rows[i].positions[j]
[1319]	647	x = sage.rings.rational.Rational()
[0]	648	x.set_from_mpq(self.rows[i].entries[j])
	649	X[n] = x
	650	return X
	651
	652	def dict(self):
	653	"""
	654	Return an associative arrow of pairs
	655	(i,j):x
	656	where the keys (i,j) run through the nonzero positions of the matrix
	657	and the x are nonzero and of type Integer.
	658	"""
	659	cdef int i, j, n
[1319]	660	cdef sage.rings.rational.Rational x
[0]	661	cdef mpq_t t
	662	cdef object entries
	663	X = {}
	664	for i from 0 <= i < self.nr:
	665	if PyErr_CheckSignals(): raise KeyboardInterrupt
	666	for j from 0 <= j < self.rows[i].num_nonzero:
	667	n = self.rows[i].positions[j]
[1319]	668	x = sage.rings.rational.Rational()
[0]	669	x.set_from_mpq(self.rows[i].entries[j])
	670	X[(i,n)] = x
	671	return X
	672
	673	def randomize(self, int sparcity, bound=2):
	674	"""
	675	randomize(self, int sparcity):
	676
	677	The sparcity is a bound on the number of nonzeros per row.
	678	"""
	679	cdef int i, j, k, r
	680	for i from 0 <= i < self.nr:
	681	if PyErr_CheckSignals(): raise KeyboardInterrupt
	682	for j from 0 <= j <= sparcity:
	683	self[i, random.randrange(0,self.nc)] = random.randrange(-bound,bound)
	684
	685	def __repr__(self):
	686	cdef int i, j
	687	cdef mpq_t x
	688	cdef char *buf
	689
	690	mpq_init(x)
	691	s = "[\n"
	692	for i from 0 <= i < self.nr:
	693	if PyErr_CheckSignals(): raise KeyboardInterrupt
	694	for j from 0 <= j < self.nc:
	695	mpq_vector_get_entry(&x, &self.rows[i], j)
	696	buf = mpq_get_str(NULL, 10, x)
	697	s = s + str(buf) + ", "
[1793]	698	sage_free(buf) # use c's malloc/free
[0]	699	s = s + "\n"
	700	s = s[:-3] + "\n]"
	701	mpq_clear(x)
	702	return s
	703
	704	def list(self):
	705	cdef int i
	706	X = []
	707	for i from 0 <= i < self.nr:
	708	for j, x in mpq_vector_to_list(&self.rows[i]):
	709	X.append((i,j,x))
	710	return X
	711
	712	def __getitem__(self, t):
	713	if not isinstance(t, tuple) or len(t) != 2:
[11]	714	raise IndexError, "Index (=%s) of matrix access must be a row and a column."%t
[1319]	715	cdef sage.rings.rational.Rational y
[0]	716	cdef int i, j
	717	i, j = t
	718	cdef mpq_t x
	719	mpq_init(x)
	720	mpq_vector_get_entry(&x, &self.rows[i], j)
[1319]	721	y = sage.rings.rational.Rational()
[0]	722	y.set_from_mpq(x)
	723	mpq_clear(x)
	724	return y
	725
	726	def __setitem__(self, t, x):
	727	if not isinstance(t, tuple) or len(t) != 2:
[11]	728	raise IndexError, "index (=%s) for setting matrix item must be a 2-tuple."%t
[0]	729	cdef int i, j
	730	i, j = t
	731	if i<0 or i >= self.nr or j<0 or j >= self.nc:
	732	raise IndexError, "Array index out of bounds."
	733	cdef object s
	734	s = str(x)
	735	mpq_vector_set_entry_str(&self.rows[i], j, s)
[11]	736
[1593]	737	def matrix_multiply(self, Matrix_rational_sparse B):
[11]	738	"""
	739	Return the matrix product of self and B.
	740	"""
	741	cdef int i, j, k
	742
	743	cdef mpq_t x, t
	744	mpq_init(x)
	745	mpq_init(t)
	746
[1593]	747	cdef Matrix_rational_sparse A
	748	A = Matrix_rational_sparse(self.nr, B.nc)
[11]	749
	750	for i from 0 <= i < self.nr:
	751	if sage.misc.all.get_verbose()>=3 and i%50 == 0:
	752	sage.misc.all.verbose('row %s of %s'%(i, self.nr), level=3)
	753	for j from 0 <= j < B.nc:
	754	# dot of i-th row with j-th column
	755	mpq_set_si(x, 0, 1)
	756	for k from 0 <= k < self.rows[i].num_nonzero:
	757	mpq_vector_get_entry(&t, &B.rows[self.rows[i].positions[k]], j)
	758	if mpz_cmp_si(mpq_numref(t), 0) != 0: # is nonzero
	759	mpq_mul(t, t, self.rows[i].entries[k])
	760	mpq_add(x, x, t)
	761	if mpz_cmp_si(mpq_numref(x), 0) != 0:
	762	mpq_vector_set_entry(&A.rows[i], j, x)
	763
	764	mpq_clear(x)
	765	mpq_clear(t)
	766	return A
	767
[0]	768
	769	def nrows(self):
	770	return self.nr
	771
	772	def ncols(self):
	773	return self.nc
	774
[11]	775	def matrix_modint(self, int n, denoms=True):
[10]	776	"""
	777	Return reduction of this matrix modulo the integer $n$.
[11]	778
	779	INPUT:
	780	n -- int
	781	denoms -- bool (default: True) if True reduce denominators;
	782	if False assume all denoms are 1.
[10]	783	"""
[11]	784	cdef int i, j, d
[1593]	785	cdef matrix_modn_sparse.Matrix_modn_sparse A
	786	#cdef matrix_sparse.Matrix_sparse A
[10]	787	cdef unsigned int num, den
	788	cdef mpq_vector* v
[11]	789	d = denoms
[10]	790
[1593]	791	A = matrix_modn_sparse.Matrix_modn_sparse(MatrixSpace(GF(n), self.nr, self.nc), n, self.nr, self.nc)
[10]	792	for i from 0 <= i < self.nr:
	793	v = &self.rows[i]
	794	for j from 0 <= j < v.num_nonzero:
	795	if mpz_cmp_si(mpq_denref(v.entries[j]), 1) == 0:
[1593]	796	#matrix_rational_sparse.pyx:1480:30: Cannot assign type 'c_vector_modint ()' to 'c_vector_modint ()' ?!
	797	set_entry(<c_vector_modint *>&A.rows[i], v.positions[j],
[10]	798	mpz_fdiv_ui(mpq_numref(v.entries[j]), n))
	799	else:
	800	num = mpz_fdiv_ui(mpq_numref(v.entries[j]), n)
[11]	801	if denoms:
	802	den = mpz_fdiv_ui(mpq_denref(v.entries[j]), n)
[1593]	803	set_entry(<c_vector_modint *>&A.rows[i], v.positions[j],
[11]	804	int((num * ai.inverse_mod_int(den, n)) % n))
	805	else:
[1593]	806	set_entry(<c_vector_modint *>&A.rows[i], v.positions[j], num)
[11]	807
[10]	808	return A
	809
[0]	810	def swap_rows(self, int n1, int n2):
	811	"""
	812	Swap the rows in positions n1 and n2
	813	"""
	814	if n1 < 0 or n1 >= self.nr or n2 < 0 or n2 >= self.nr:
[11]	815	raise IndexError, "Invalid row number (n1=%s, n2=%s)"%(n1,n2)
[0]	816	if n1 == n2:
	817	return
	818	cdef mpq_vector tmp
	819	tmp = self.rows[n1]
	820	self.rows[n1] = self.rows[n2]
	821	self.rows[n2] = tmp
[11]	822
[0]	823
[11]	824
	825	def height(self, scale=1):
	826	"""
	827	Returns the height of scale*self, which is the maximum of the
	828	absolute values of all numerators and denominators of the
	829	entries of scale*self.
	830
	831	OUTPUT:
	832	-- Integer
	833
	834	NOTE: Since 0 = 0/1 has denominator 1, the height is at least 1.
	835
	836	EXAMPLES:
	837	sage: A = MatrixSpace(QQ, 3,3, sparse=True)([1,2,3,4,5/13,6,-7/17,8,9])._sparse_matrix_mpq_()
	838	sage: A.height()
	839	17
	840	sage: A = MatrixSpace(QQ, 2,2, sparse=True)([1,2,-197/13,4])._sparse_matrix_mpq_()
	841	sage: A.height()
	842	197
	843	sage: A.height(26)
	844	394
	845	"""
	846	cdef mpz_t h, a
	847	mpz_init(h)
	848	mpz_init(a)
	849	mpz_set_si(h, 1)
	850
	851	cdef mpq_t s, v2
[1319]	852	cdef sage.rings.rational.Rational tmp
	853	tmp = sage.rings.rational.Rational(abs(scale))
[11]	854	mpq_init(s)
	855	mpq_set(s,tmp.value)
	856
	857	cdef int i, j
	858	cdef mpq_vector* v
	859
	860	if scale == 1:
	861	for i from 0 <= i < self.nr:
	862	v = &self.rows[i]
	863	for j from 0 <= j < v.num_nonzero:
	864	mpz_abs(a, mpq_numref(v.entries[j]))
	865	if mpz_cmp(a, h) > 0:
	866	mpz_set(h, a)
	867	if mpz_cmp(mpq_denref(v.entries[j]), h) > 0:
	868	mpz_set(h, mpq_denref(v.entries[j]))
	869	else:
	870	mpq_init(v2)
	871	for i from 0 <= i < self.nr:
	872	v = &self.rows[i]
	873	for j from 0 <= j < v.num_nonzero:
	874	mpq_mul(v2, v.entries[j], s)
	875	mpz_abs(a, mpq_numref(v2))
	876	if mpz_cmp(a, h) > 0:
	877	mpz_set(h, a)
	878	if mpz_cmp(mpq_denref(v2), h) > 0:
	879	mpz_set(h, mpq_denref(v2))
	880	mpq_clear(v2)
	881
	882	#endif
	883	mpz_clear(a)
	884	mpq_clear(s)
	885
[1319]	886	cdef sage.rings.integer.Integer r
	887	r = sage.rings.integer.Integer()
[11]	888	r.set_from_mpz(h)
	889	mpz_clear(h)
	890	return r
	891
	892	def denom(self):
	893	"""
	894	Returns the denominator of self, which is the least common
	895	multiple of the denominators of the entries of self.
	896
	897	OUTPUT:
	898	-- Integer
	899
	900	NOTE: Since 0 = 0/1 has denominator 1, the height is at least 1.
	901
	902	EXAMPLES:
	903	sage: A = MatrixSpace(QQ, 3,3, sparse=True)([1/17,2,3,4,5/13,6,-7/17,8,9])._sparse_matrix_mpq_()
	904	sage: A.denom()
	905	221
	906	sage: A = MatrixSpace(QQ, 2,2, sparse=True)([1,2,-197/13,4])._sparse_matrix_mpq_()
	907	sage: A.denom()
	908	13
	909	"""
	910	cdef mpz_t d
	911	mpz_init(d)
	912	mpz_set_si(d, 1)
	913
	914	cdef int i, j
	915	cdef mpq_vector* v
	916
	917	_sig_on
	918	for i from 0 <= i < self.nr:
	919	v = &self.rows[i]
	920	for j from 0 <= j < v.num_nonzero:
	921	mpz_lcm(d, d, mpq_denref(v.entries[j]))
	922	_sig_off
	923
[1319]	924	cdef sage.rings.integer.Integer _d
	925	_d = sage.rings.integer.Integer()
[11]	926	_d.set_from_mpz(d)
	927	mpz_clear(d)
	928	return _d
	929
	930	def clear_denom(self):
	931	"""
	932	Replace self by self times the denominator of self.
	933
	934	EXAMPLES:
	935	sage: A = MatrixSpace(QQ,2, sparse=True)([1/3,2,3/2,4])._sparse_matrix_mpq_()
	936	sage: A.clear_denom ()
	937	6
	938	sage: A
	939	[
	940	2, 12,
	941	9, 24
	942	]
	943	"""
[1319]	944	cdef sage.rings.integer.Integer _d
[11]	945	_d = self.denom()
	946
	947	cdef mpq_t d
	948	mpq_init(d)
	949	mpz_set(mpq_numref(d), _d.value)
	950
	951	cdef int i, j
	952	cdef mpq_vector* v
	953
	954	_sig_on
	955	for i from 0 <= i < self.nr:
	956	v = &self.rows[i]
	957	for j from 0 <= j < v.num_nonzero:
	958	mpq_mul(v.entries[j], v.entries[j], d)
	959	_sig_off
	960
	961	mpq_clear(d)
	962	return _d
	963
	964
[1319]	965	def divide_by(self, sage.rings.integer.Integer d):
[11]	966	"""
	967	Replace self by self/d.
	968
	969	EXAMPLES:
	970	sage: A = MatrixSpace(QQ,2, sparse=True)([1,2,3,4])._sparse_matrix_mpq_()
	971	sage: A.divide_by(5)
	972	sage: A
	973	[
	974	1/5, 2/5,
	975	3/5, 4/5
	976	]
	977	"""
	978	cdef mpq_t dd
	979	mpq_init(dd)
	980	mpq_set_si(dd, 1, 1)
	981	mpz_set(mpq_denref(dd), d.value)
	982
	983	cdef int i, j
	984	cdef mpq_vector* v
	985
	986	_sig_on
	987	for i from 0 <= i < self.nr:
	988	v = &self.rows[i]
	989	for j from 0 <= j < v.num_nonzero:
	990	mpq_mul(v.entries[j], v.entries[j], dd)
	991	_sig_off
	992	mpq_clear(dd)
	993
	994	def echelon_multimodular(self, height_guess=None, proof=True):
	995	"""
	996	Returns reduced row-echelon form using a multi-modular
	997	algorithm. Does not change self.
	998
	999	INPUT:
	1000	height_guess -- integer or None
	1001	proof -- boolean (default: True)
	1002
	1003	EXAMPLES:
	1004	sage: A = MatrixSpace(QQ,2, sparse=True)([1/3,2,3/2,4])._sparse_matrix_mpq_(); A
	1005	[
	1006	1/3, 2,
	1007	3/2, 4
	1008	]
	1009	sage: B = A.echelon_multimodular(); B
	1010	[
	1011	1, 0,
	1012	0, 1
	1013	]
	1014
	1015	Note that A is unchanged.
	1016	sage: A
	1017	[
	1018	1/3, 2,
	1019	3/2, 4
	1020	]
	1021
	1022
	1023	ALGORITHM:
	1024	The following is a modular algorithm for computing the echelon
	1025	form. Define the height of a matrix to be the max of the
	1026	absolute values of the entries.
	1027
	1028	Given Matrix A with n columns (self).
	1029
	1030	0. Rescale input matrix A to have integer entries. This does
	1031	not change echelon form and makes reduction modulo lots of
	1032	primes significantly easier if there were denominators.
	1033	Henceforth we assume A has integer entries.
	1034
	1035	1. Let c be a guess for the height of the echelon form. E.g.,
	1036	c=1000, e.g., if matrix is very sparse and application is to
	1037	computing modular symbols.
	1038
	1039	2. Let M = n * c * H(A) + 1,
	1040	where n is the number of columns of A.
	1041
	1042	3. List primes p_1, p_2, ..., such that the product of
	1043	the p_i is at least M.
	1044
	1045	4. Try to compute the rational reconstruction CRT echelon form
	1046	of A mod the product of the p_i. If rational
	1047	reconstruction fails, compute 1 more echelon forms mod the
	1048	next prime, and attempt again. Make sure to keep the
	1049	result of CRT on the primes from before, so we don't have
	1050	to do that computation again. Let E be this matrix.
	1051
	1052	5. Compute the denominator d of E.
	1053	Attempt to prove that result is correct by checking that
	1054
	1055	H(dE)ncols(A)*H(A) < (prod of reduction primes)
	1056
	1057	where H denotes the height. If this fails, do step 4 with
	1058	a few more primes.
	1059	"""
[1593]	1060	cdef Matrix_rational_sparse E
[11]	1061
	1062	if self.nr == 0 or self.nc == 0:
	1063	return self
	1064
[1319]	1065	cdef sage.rings.integer.Integer dd
[11]	1066	dd = self.clear_denom()
	1067	#print "dd = ", dd
	1068
	1069	hA = long(self.height())
[807]	1070	if height_guess is None:
[11]	1071	height_guess = 100000hA*4
	1072	tm = sage.misc.all.verbose("height_guess = %s"%height_guess, level=2)
	1073
	1074	if proof:
	1075	M = self.nc * height_guess * hA + 1
	1076	else:
	1077	M = height_guess + 1
	1078
	1079	p = START_PRIME
	1080	X = []
	1081	best_pivots = []
	1082	prod = 1
	1083	problem = 0
	1084	while True:
	1085	while prod < M:
	1086	problem = problem + 1
	1087	if problem > 50:
	1088	sage.misc.all.verbose("sparse_matrix multi-modular reduce not converging?")
	1089	t = sage.misc.all.verbose("echelon modulo p=%s (%.2f%% done)"%(
	1090	p, 100*float(len(str(prod))) / len(str(M))), level=2)
	1091
	1092	# We use denoms=False, since we made self integral by calling clear_denom above.
	1093	A = self.matrix_modint(p, denoms=False)
	1094	t = sage.misc.all.verbose("time to reduce matrix mod p:",t, level=2)
	1095	A.echelon()
	1096	t = sage.misc.all.verbose("time to put reduced matrix in echelon form:",t, level=2)
[1663]	1097	c = matrix_rational_dense.cmp_pivots(best_pivots, A.pivots())
[11]	1098	if c <= 0:
	1099	best_pivots = A.pivots()
	1100	X.append(A)
	1101	prod = prod * p
	1102	else:
	1103	# do not save A since it is bad.
	1104	if sage.misc.all.LEVEL > 1:
	1105	sage.misc.all.verbose("Excluding this prime (bad pivots).")
	1106	p = sage.rings.arith.next_prime(p)
	1107	t = sage.misc.all.verbose("time for pivot compare", t, level=2)
	1108	# Find set of best matrices.
	1109	Y = []
	1110	# recompute product, since may drop bad matrices
	1111	prod = 1
	1112	for i in range(len(X)):
[1663]	1113	if matrix_rational_dense.cmp_pivots(best_pivots, X[i].pivots()) <= 0:
[11]	1114	Y.append(X[i])
	1115	prod = prod * X[i].prime()
	1116	try:
	1117	t = sage.misc.all.verbose("start crt and rr", level=2)
	1118	E = lift_matrices_modint(Y)
	1119	#print "E = ", E # debug
	1120	sage.misc.all.verbose("crt and rr time is",t, level=2)
	1121	except ValueError, msg:
	1122	#print msg # debug
	1123	sage.misc.all.verbose("Redoing with several more primes", level=2)
	1124	for i in range(3):
	1125	M = M * START_PRIME
	1126	continue
	1127
	1128	if not proof:
	1129	sage.misc.all.verbose("Not checking validity of result (since proof=False).", level=2)
	1130	break
	1131	d = E.denom()
	1132	hdE = long(E.height(d))
	1133	if hdE * self.ncols() * hA < prod:
	1134	break
	1135	for i in range(3):
	1136	M = M * START_PRIME
	1137	#end while
	1138	sage.misc.all.verbose("total time",tm, level=2)
	1139	self.__pivots = best_pivots
	1140	E.__pivots = best_pivots
	1141	self.divide_by(dd)
	1142	return E
	1143
[0]	1144	def echelon(self):
	1145	"""
	1146	Replace self by its reduction to reduced row echelon form.
	1147
	1148	ALGORITHM:
[10]	1149	We use Gauss elimination, which is slightly intelligent, in
	1150	these sense that we clear each column using a row with the
	1151	minimum number of nonzero entries.
[0]	1152
	1153	WARNING: There is no reason to use the code below, except
	1154	for testing this class. It is vastly faster to use
	1155	the multi-modular method, which is implemented in
	1156	sparse_matrix.Sparse_matrix_rational
	1157	"""
	1158	cdef int i, r, c, min, min_row, start_row, sgn
	1159	cdef mpq_vector tmp
	1160	cdef mpq_t a, a_inverse, b, minus_b
	1161	mpq_init(a)
	1162	mpq_init(a_inverse)
	1163	mpq_init(b)
	1164	mpq_init(minus_b)
	1165
	1166	start_row = 0
	1167	self.__pivots = []
	1168	for c from 0 <= c < self.nc:
[11]	1169	_sig_check
	1170	if c % 10 == 0:
	1171	sage.misc.all.verbose('clearing column %s of %s'%(c,self.nc),
[1663]	1172	caller_name = 'matrix_rational_sparse echelon')
[0]	1173	min = self.nc + 1
	1174	min_row = -1
	1175	if PyErr_CheckSignals(): raise KeyboardInterrupt
	1176	for r from start_row <= r < self.nr:
	1177	if self.rows[r].num_nonzero > 0 and self.rows[r].num_nonzero < min:
	1178	# Since there is at least one nonzero entry, the first entry
	1179	# of the positions list is defined. It is the first position
	1180	# of a nonzero entry, and it equals c precisely if row r
	1181	# is a row we could use to clear column c.
	1182	if self.rows[r].positions[0] == c:
	1183	min_row = r
	1184	min = self.rows[r].num_nonzero
	1185	#endif
	1186	#endif
	1187	#endfor
	1188	if min_row != -1:
	1189	r = min_row
	1190	self.__pivots.append(c)
	1191	# Since we can use row r to clear column c, the
	1192	# entry in position c in row r must be the first
	1193	# nonzero entry.
	1194	mpq_inv(a_inverse, self.rows[r].entries[0])
	1195	scale_mpq_vector(&self.rows[r], a_inverse)
	1196	self.swap_rows(r, start_row)
	1197	for i from 0 <= i < self.nr:
	1198	if i != start_row:
	1199	mpq_vector_get_entry(&b, &self.rows[i], c)
	1200	if mpq_sgn(b) != 0: # if b is nonzero
	1201	mpq_neg(minus_b, b)
	1202	add_mpq_vector_init(&tmp, &self.rows[i],
	1203	&self.rows[start_row], minus_b)
	1204	clear_mpq_vector(&self.rows[i])
	1205	self.rows[i] = tmp
	1206	start_row = start_row + 1
	1207	mpq_clear(a)
	1208	mpq_clear(a_inverse)
	1209	mpq_clear(b)
	1210	mpq_clear(minus_b)
	1211
[11]	1212
	1213
	1214	#####################################
	1215
	1216
	1217
[1593]	1218	def Matrix_rational_sparse_from_columns(columns):
[0]	1219	"""
[1593]	1220	Create a sparse Matrix_rational_sparse from a list of sparse Vector_mpq's.
[0]	1221	Each vector must have same degree.
	1222	INPUT:
	1223	columns -- a list of Vector_mpq's.
	1224	OUTPUT:
[1593]	1225	A sparse Matrix_rational_sparse whose columns are as given.
[0]	1226	"""
	1227	if not isinstance(columns, list):
	1228	raise TypeError, "columns must be a list"
	1229
	1230	cdef int j, nc, nr
	1231	nc = len(columns)
	1232	if nc == 0:
[1593]	1233	return Matrix_rational_sparse(0,0)
[0]	1234	if not isinstance(columns[0], Vector_mpq):
	1235	raise TypeError, "each column must be of type Vector_mpq"
	1236	nr = columns[0].degree()
	1237	entries = []
	1238	for j from 0 <= j < nc:
	1239	v = columns[j]
	1240	if not isinstance(v, Vector_mpq):
	1241	raise TypeError, "each column must be of type Vector_mpq"
	1242	if v.degree() != nr:
	1243	raise IndexError, "each column must have degree the number of rows of self."
	1244	for i, x in v.list():
	1245	# now the i,j entry of our matrix should be set equal to x."
	1246	entries.append((i,j,x))
[1593]	1247	return Matrix_rational_sparse(nr, nc, entries)
[0]	1248
[11]	1249
	1250
	1251
	1252	def lift_matrices_modint(X):
	1253	"""
	1254	INPUT:
	1255	X -- list of sparse Matrix_modint matrices.
	1256
	1257	OUTPUT:
[1593]	1258	Matrix_rational_sparse -- computed if possible using rational reconstruction
[11]	1259	raises ValueError if there is no valid lift.
	1260
	1261	ALGORITHM:
	1262	0. validate input -- type of elements of X and dimensions match up.
	1263	1. compute CRT basis as array of mpz_t's
	1264	2. allocate new matrix
	1265	3. for each row:
	1266	* find union of nonzero positions as Python int set/list
	1267	* allocate relevant memory in matrix over Q that we're constructing
	1268	* for each nonzero position:
	1269	- use CRT basis to lift to mod prod of moduli
	1270	- use rational reconstruction to list to Q
	1271	(if fail clean up memory and raise exception)
	1272	-- always multiply element by lcm of denominators so far by attempting
	1273	rr (and only do rr if there is a denominator)
	1274	- enter value in output matrix
	1275	4. memory cleanup:
	1276	* clear crt basis
	1277	"""
	1278	cdef int i, r, c, nr, nc, lenX
[1319]	1279	cdef sage.rings.integer.Integer a, b
[11]	1280
	1281	#print "X = ", X # debug
	1282
	1283	# 0. validate input
	1284	lenX = len(X) # save as C int
	1285	# Create C level access to the entries of the Python array X
	1286	for i from 0 <= i < lenX:
	1287	if i == 0:
	1288	nr = X[i].nrows()
	1289	nc = X[i].ncols()
	1290	else:
	1291	if X[i].nrows() != nr or X[i].ncols() != nc:
	1292	raise ValueError, "number of rows and columns of all input matrices must be the same"
	1293	#print "validated input"
	1294
	1295	# 1. crt basis as mpz_t'
[1319]	1296	# We use Python since this doesn't have to be fast.
	1297	import sage.rings.integer # I don't know why this is needed, but it is. (must be a weird pyrex thing)
[11]	1298	P = []
[1319]	1299	for A in X:
	1300	a = sage.rings.integer.Integer(A.prime())
	1301	P.append(a) # no list comprehensions in pyrex
	1302
[11]	1303	import sage.rings.integer_ring
	1304	_B = sage.rings.integer_ring.crt_basis(P)
	1305	#print "computed crt basis = ", _B
	1306
	1307	cdef mpz_t* B
	1308	cdef mpz_t prod # product of moduli
	1309
	1310	mpz_init(prod)
	1311	mpz_set_si(prod, 1)
	1312
[1793]	1313	B = <mpz_t> sage_malloc(sizeof(mpz_t) len(_B))
[11]	1314	if B == <mpz_t*> 0:
	1315	raise MemoryError, "Error allocating memory"
	1316	for i from 0 <= i < len(_B):
	1317	mpz_init(B[i])
	1318	b = _B[i]
	1319	mpz_set(B[i], b.value)
	1320	mpz_mul_si(prod, prod, X[i].prime())
	1321
	1322	#print "converted crt basis to C"
	1323
	1324
	1325	# 2. create un-allocated matrix over Q
[1593]	1326	cdef Matrix_rational_sparse M
	1327	M = Matrix_rational_sparse(nr, nc, init=False)
[11]	1328	#print "created un-allocated matrix"
	1329
	1330	cdef int num_nonzero
	1331	cdef mpz_t crt_lift, tmp
	1332
	1333	mpz_init(crt_lift)
	1334	mpz_init(tmp)
	1335
	1336	cdef mpq_t denom
	1337	mpq_init(denom)
	1338	mpq_set_si(denom, 1, 1)
	1339
	1340	# 3. for each row...
[1593]	1341	cdef matrix_modn_sparse.Matrix_modn_sparse C
[11]	1342
	1343	error = None
	1344
	1345	cdef int max_row_allocated
	1346
	1347	try:
	1348	for r from 0 <= r < nr:
	1349	#print "considering row r=%s"%r
	1350	#_sig_check
	1351	nonzero_positions = []
	1352	for i from 0 <= i < lenX:
	1353	C = X[i]
	1354	for j from 0 <= j < C.rows[r].num_nonzero:
	1355	nonzero_positions.append(C.rows[r].positions[j])
	1356	nonzero_positions = list(set(nonzero_positions))
	1357	nonzero_positions.sort()
	1358	num_nonzero = len(nonzero_positions)
	1359	#print "num_nonzero = ", num_nonzero
	1360
	1361	# allocate space for that many nonzero entries.
	1362	init_mpq_vector(&M.rows[r], nc, num_nonzero)
	1363	max_row_allocated = r
	1364	for j from 0 <= j < num_nonzero:
	1365	M.rows[r].positions[j] = nonzero_positions[j]
	1366
	1367	for j from 0 <= j < num_nonzero:
	1368	# compute each lifted element and insert into M.rows[r]
	1369	mpz_set_si(crt_lift, 0)
	1370	for i from 0 <= i < lenX:
	1371	C = X[i] # this could slow everything down a lot... ?
	1372	# tmp = (ith CRT basis) * (entry of ith modint matrix)
[1593]	1373	#matrix_rational_sparse.pyx:2058:58: Cannot assign type 'c_vector_modint ()' to 'c_vector_modint ()' ?!
	1374	mpz_mul_si(tmp, B[i], get_entry(<c_vector_modint *>&C.rows[r], M.rows[r].positions[j]))
[11]	1375	mpz_add(crt_lift, crt_lift, tmp)
	1376
	1377	mpz_mul(crt_lift, crt_lift, denom)
	1378	mpq_rational_reconstruction(M.rows[r].entries[j],
	1379	crt_lift, prod)
	1380	mpq_div(M.rows[r].entries[j], M.rows[r].entries[j], denom)
	1381	mpz_lcm(mpq_numref(denom), mpq_numref(denom), mpq_denref(M.rows[r].entries[j]))
	1382
	1383
	1384	except ValueError, msg:
	1385	#print "msg = ", msg
	1386	error = ValueError
	1387	# Delete memory in M, since M.is_init isn't true, so M would never
	1388	# get deleted. We do this here, since otherwise we'd have to
	1389	# finish constructing M just so the __dealloc__ method would work correctly.
	1390	for r from 0 <= r <= max_row_allocated:
	1391	clear_mpq_vector(&M.rows[r])
	1392
	1393	# 4. memory cleanup
	1394	#print "cleaning up memory"
	1395	mpz_clear(crt_lift)
	1396	mpz_clear(denom)
	1397	mpz_clear(prod)
	1398	mpz_clear(tmp)
	1399
	1400	for i from 0 <= i < len(_B):
	1401	mpz_clear(B[i])
[1793]	1402	sage_free(B)
[11]	1403
	1404	if not error is None:
	1405	raise error, msg
	1406
	1407	# done.
	1408	M.is_init = True
	1409	return M
	1410
	1411
[1319]	1412	def make_sparse_rational_matrix(nrows, ncols, entries):
[1593]	1413	return Matrix_rational_sparse(nrows, ncols, entries = entries, init=False)

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

Context Navigation

source: sage/matrix/matrix_rational_sparse.pyx @ 3092:227119d39b73

Download in other formats: