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

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

matrix coding sprint -- snapshot.

RevLine
[1775]1r"""nodoctest
[1687]2Sparse matrices over \$\Q\$
[6]3
[3092]4This is an optimzed implementation of sparse matrices over
5the rational numbers.
[771]6
[6]7AUTHOR:
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
13EXAMPLES:
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>
45#  The full text of the GPL is available at:
[10]47#############################################################################
[0]48
49import random
50
[11]51import sage.misc.all
[1593]52import sage.matrix.matrix_rational_sparse
[11]53import sage.rings.integer_ring
[823]54import sage.rings.finite_field
[11]55import sage.rings.arith
[0]56
[1319]57cimport sage.rings.rational
58import  sage.rings.rational
[11]59
[1319]60cimport sage.rings.integer
61import  sage.rings.integer
[11]62
[1319]63cimport sage.ext.arith
64import sage.ext.arith
65cdef sage.ext.arith.arith_int ai
66ai = sage.ext.arith.arith_int()
[0]67
[1664]68cimport matrix_field_sparse
69import matrix_field_sparse
[0]70
[1593]71cimport matrix_modn_sparse
72import matrix_modn_sparse
73
74cimport matrix_rational_dense
75import matrix_rational_dense
76
[1319]77include "../ext/gmp.pxi"
[3092]78include '../ext/interrupt.pxi'
[0]79
[10]80START_PRIME = 20011  # used for multi-modular algorithms
[0]81
82cdef class Vector_mpq
83
84cdef void Vector_mpq_rescale(Vector_mpq w, mpq_t x):
85    scale_mpq_vector(&w.v, x)
86
87cdef 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
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]228cdef 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(nrows*sizeof(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])
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 = 2*n + 20*ln
[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]
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])
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(d*E)*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 = 100000*hA**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)
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]1218def 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
1252def 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]))
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]1412def 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.