source: sage/matrix/matrix_integer_dense.pyx @ 3110:91db895cb507

Revision 3110:91db895cb507, 52.5 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

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