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

Revision 2830:bf1870252c66, 45.0 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

linear algebra fixes for linux for now

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