source: sage/matrix/matrix_rational_dense.pyx @ 2654:237abdf81fb6

Revision 2654:237abdf81fb6, 16.8 KB checked in by Robert Bradshaw <robertwb@…>, 6 years ago (diff)

enable multi-modular for rationals
Line 
1"""
2Dense matrices over the rational field.
3"""
4
5##############################################################################
6#       Copyright (C) 2004,2005,2006 William Stein <wstein@gmail.com>
7#  Distributed under the terms of the GNU General Public License (GPL)
8#  The full text of the GPL is available at:
9#                  http://www.gnu.org/licenses/
10##############################################################################
11
12include "../ext/interrupt.pxi"
13include "../ext/stdsage.pxi"
14include "../ext/cdefs.pxi"
15
16from sage.rings.rational cimport Rational
17from matrix cimport Matrix
18from matrix_integer_dense cimport Matrix_integer_dense
19import sage.structure.coerce
20from sage.structure.element cimport ModuleElement
21from sage.rings.integer cimport Integer
22from sage.rings.integer_ring import ZZ
23
24cdef class Matrix_rational_dense(matrix_dense.Matrix_dense): 
25
26    ########################################################################
27    # LEVEL 1 functionality
28    # x * __new__ 
29    # x * __dealloc__   
30    # x * __init__     
31    # x * set_unsafe
32    # x * get_unsafe
33    # x * cdef _pickle
34    # x * cdef _unpickle
35    ########################################################################
36    def __new__(self, parent, entries, copy, coerce):
37        """
38        Create and allocate memory for the matrix.
39
40        Unlike over matrix_integer_dense, mpq_init() is called (as there is no mpq_init_set function).
41       
42        INPUT:
43            parent, entries, coerce, copy -- as for __init__.
44
45        EXAMPLES:
46            sage: from sage.matrix.matrix_rational_dense import Matrix_rational_dense
47            sage: a = Matrix_rational_dense.__new__(Matrix_rational_dense, Mat(ZZ,3), 0,0,0)
48            sage: type(a)
49            <type 'sage.matrix.matrix_rational_dense.Matrix_rational_dense'>
50
51        WARNING: This is for internal use only, or if you really know what you're doing.
52        """
53        matrix_dense.Matrix_dense.__init__(self, parent)
54
55        cdef Py_ssize_t i, k
56
57        self._entries = <mpq_t *> sage_malloc(sizeof(mpq_t)*(self._nrows * self._ncols))
58        if self._entries == NULL:
59            raise MemoryError, "out of memory allocating a matrix"
60
61        self._matrix =  <mpq_t **> sage_malloc(sizeof(mpq_t*) * self._nrows)
62        if self._matrix == NULL:
63            raise MemoryError, "out of memory allocating a matrix"
64
65        # store pointers to the starts of the rows
66        k = 0
67        for i from 0 <= i < self._nrows:
68            self._matrix[i] = self._entries + k
69            k = k + self._ncols
70
71        for i from 0 <= i < self._nrows * self._ncols:
72            mpq_init(self._entries[i])
73
74    def  __dealloc__(self):
75        cdef Py_ssize_t i
76        for i from 0 <= i < self._nrows * self._ncols:
77            mpq_clear(self._entries[i])
78        sage_free(self._entries)
79        sage_free(self._matrix)
80   
81    def __init__(self, parent, entries=0, coerce=True, copy=True):
82   
83        cdef Py_ssize_t i
84        cdef Rational z
85
86        if isinstance(entries, list):
87            if len(entries) != self._nrows * self._ncols:
88                raise TypeError, "entries has the wrong length"
89
90            _sig_on
91            if coerce:
92                for i from 0 <= i < self._nrows * self._ncols:
93                    # TODO: Should use an unsafe un-bounds-checked array access here.
94                    z = Rational(entries[i])
95                    mpq_set(self._entries[i], z.value)
96            else:
97                for i from 0 <= i < self._nrows * self._ncols:
98                    # TODO: Should use an unsafe un-bounds-checked array access here.
99                    mpq_set(self._entries[i], (<Rational> entries[i]).value)
100            _sig_off
101                   
102        else:
103            # is it a scalar?
104            try:
105                # Try to coerce entries to a scalar (an integer)
106                z = Rational(entries)
107                is_list = False
108            except TypeError:
109                raise TypeError, "entries must be coercible to a list or integer"
110
111            if not z.is_zero():
112                if self._nrows != self._ncols:
113                    raise TypeError, "nonzero scalar matrix must be square"
114                for i from 0 <= i < self._nrows:
115                    mpq_set(self._entries[i*self._ncols+i], z.value)
116
117
118    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
119        cdef Rational y
120        y = value
121        mpq_set(self._matrix[i][j], y.value)
122
123
124    cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
125        cdef Rational x
126        x = Rational.__new__(Rational)
127        mpq_set(x.value, self._matrix[i][j])
128        return x
129
130    def _pickle(self):
131        return self._pickle_version0(), 0
132   
133    def _unpickle(self, data, int version):
134        if version == 0:
135            self._unpickle_version0(data)
136        else:
137            raise RuntimeError, "unknown matrix version (=%s)"%version
138       
139    cdef _pickle_version0(self):
140        cdef Py_ssize_t i, j, len_so_far, m, n
141        cdef char *a
142        cdef char *s, *t, *tmp
143
144        if self._nrows == 0 or self._ncols == 0:
145            data = ''
146        else:
147            n = self._nrows*self._ncols*10
148            s = <char*> sage_malloc(n * sizeof(char))
149            t = s
150            len_so_far = 0
151
152            _sig_on
153            for i from 0 <= i < self._nrows:
154                for j from 0 <= j < self._ncols:
155                    m = mpz_sizeinbase (mpq_numref(self._matrix[i][j]), 32) + \
156                        mpz_sizeinbase (mpq_denref(self._matrix[i][j]), 32) + 3
157                    if len_so_far + m + 1 >= n:
158                        # copy to new string with double the size
159                        n = 2*n + m + 1
160                        tmp = <char*> sage_malloc(n * sizeof(char))
161                        strcpy(tmp, s)
162                        sage_free(s)
163                        s = tmp
164                        t = s + len_so_far
165                    #endif
166                    mpq_get_str(t, 32, self._matrix[i][j])
167                    m = strlen(t)
168                    len_so_far = len_so_far + m + 1
169                    t = t + m
170                    t[0] = <char>32
171                    t[1] = <char>0
172                    t = t + 1
173            _sig_off
174            data = str(s)[:-1]
175            free(s)
176        return data
177
178    cdef _unpickle_version0(self, data):
179        cdef Py_ssize_t i, n
180        data = data.split()
181        n = self._nrows * self._ncols
182        if len(data) != n:
183            raise RuntimeError, "invalid pickle data."
184        for i from 0 <= i < n:
185            s = data[i]
186            if mpq_set_str(self._entries[i], s, 32):
187                raise RuntimeError, "invalid pickle data"
188       
189    def __richcmp__(Matrix self, right, int op):
190        return self._richcmp(right, op)
191    def __hash__(self):
192        return self._hash()
193
194    ########################################################################
195    # LEVEL 2 functionality
196    # x * cdef _add_c_impl         
197    #   * cdef _mul_c_impl
198    #   * cdef _cmp_c_impl
199    # x * __neg__
200    #   * __invert__
201    # x * __copy__
202    #   * _multiply_classical
203    #   * _list -- list of underlying elements (need not be a copy)
204    #   * _dict -- sparse dictionary of underlying elements (need not be a copy)
205    ########################################################################
206   
207    cdef ModuleElement _add_c_impl(self, ModuleElement right):
208        """
209        Add two dense matrices over QQ.
210
211        EXAMPLES:
212        sage: a = MatrixSpace(QQ,3)(range(9))
213        sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)])
214        sage: a+b
215        [   1  3/2  7/3]
216        [13/4 21/5 31/6]
217        [43/7 57/8 73/9]
218        sage: b.swap_rows(1,2)
219        sage: #a+b
220       
221        """
222        cdef Py_ssize_t i, j
223        cdef Matrix_rational_dense M
224        M = Matrix_rational_dense.__new__(Matrix_rational_dense, self._parent, None, None, None)
225
226        cdef mpq_t *M_row
227        cdef mpq_t *self_row
228        cdef mpq_t *right_row
229        _sig_on
230        for i from 0 <= i < self._nrows:
231            M_row = M._matrix[i]
232            self_row = self._matrix[i]
233            right_row = (<Matrix_rational_dense>right)._matrix[i]
234            for j from 0 <= j < self._ncols:
235                mpq_add(M_row[0], self_row[0], right_row[0])
236                M_row = M_row + 1
237                self_row = self_row + 1
238                right_row = right_row + 1
239        _sig_off
240        return M
241       
242    cdef ModuleElement _sub_c_impl(self, ModuleElement right):
243        """
244        Add two dense matrices over QQ.
245
246        EXAMPLES:
247        sage: a = MatrixSpace(QQ,3)(range(9))
248        sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)])
249        sage: a-b
250        [  -1  1/2  5/3]
251        [11/4 19/5 29/6]
252        [41/7 55/8 71/9]
253        """
254        cdef Py_ssize_t i, j
255        cdef Matrix_rational_dense M
256        M = Matrix_rational_dense.__new__(Matrix_rational_dense, self._parent, None, None, None)
257
258        cdef mpq_t *M_row
259        cdef mpq_t *self_row
260        cdef mpq_t *right_row
261        _sig_on
262        for i from 0 <= i < self._nrows:
263            M_row = M._matrix[i]
264            self_row = self._matrix[i]
265            right_row = (<Matrix_rational_dense>right)._matrix[i]
266            for j from 0 <= j < self._ncols:
267                mpq_sub(M_row[0], self_row[0], right_row[0])
268                M_row = M_row + 1
269                self_row = self_row + 1
270                right_row = right_row + 1
271        _sig_off
272        return M
273
274    def __neg__(self):
275        """
276        Negate a matrix over QQ.
277
278        EXAMPLES:
279        sage: a = MatrixSpace(QQ,3)([1/n for n in range(1,10)])
280        sage: -a
281        [  -1 -1/2 -1/3]
282        [-1/4 -1/5 -1/6]
283        [-1/7 -1/8 -1/9]
284        """
285        cdef Py_ssize_t i, j
286        cdef Matrix_rational_dense M
287        M = Matrix_rational_dense.__new__(Matrix_rational_dense, self._parent, None, None, None)
288
289        cdef mpq_t *M_row
290        cdef mpq_t *self_row
291        _sig_on
292        for i from 0 <= i < self._nrows:
293            M_row = M._matrix[i]
294            self_row = self._matrix[i]
295            for j from 0 <= j < self._ncols:
296                mpq_neg(M_row[0], self_row[0])
297                M_row = M_row + 1
298                self_row = self_row + 1
299        _sig_off
300        return M
301       
302    def __copy__(self):
303        """
304        Negate a matrix over QQ.
305
306        EXAMPLES:
307        sage: a = MatrixSpace(QQ,3)([1/n for n in range(1,10)])
308        sage: -a
309        [  -1 -1/2 -1/3]
310        [-1/4 -1/5 -1/6]
311        [-1/7 -1/8 -1/9]
312        """
313        cdef Py_ssize_t i, j
314        cdef Matrix_rational_dense M
315        M = Matrix_rational_dense.__new__(Matrix_rational_dense, self._parent, None, None, None)
316
317        cdef mpq_t *M_row
318        cdef mpq_t *self_row
319        _sig_on
320        for i from 0 <= i < self._nrows:
321            M_row = M._matrix[i]
322            self_row = self._matrix[i]
323            for j from 0 <= j < self._ncols:
324                mpq_set(M_row[0], self_row[0])
325                M_row = M_row + 1
326                self_row = self_row + 1
327        _sig_off
328        return M
329
330
331   
332    # cdef _mul_c_impl(self, Matrix right):
333    # cdef int _cmp_c_impl(self, Matrix right) except -2:
334    # def __invert__(self):
335    # def _multiply_classical(left, matrix.Matrix _right):
336    # def _list(self):
337    # def _dict(self):   
338   
339
340    ########################################################################
341    # LEVEL 3 functionality (Optional)
342    # x * cdef _sub_c_impl
343    #   * __deepcopy__
344    #   * __invert__
345    #   * Matrix windows -- only if you need strassen for that base
346    #   * Other functions (list them here):
347    # x * denom(self):
348    # x * mpz_denom(self, mpz_t d):
349    # x * _clear_denom(self):
350    # x * _multiply_multi_modular(self, Matrix_rational_dense right):
351    ########################################################################
352    def denom(self):
353        """
354        Return the denominator of this matrix.
355
356        OUTPUT:
357            -- SAGE Integer
358
359        EXAMPLES:
360            sage: b = matrix(QQ,2,range(6)); b[0,0]=-5007/293; b
361            [-5007/293         1         2]
362            [        3         4         5]
363            sage: b.denom()
364            293
365        """
366        cdef Integer z
367        z = Integer.__new__(Integer)
368        self.mpz_denom(z.value)
369        return z
370
371    cdef int mpz_denom(self, mpz_t d) except -1:
372        mpz_set_si(d,1)
373        cdef int i, j
374        cdef mpq_t *self_row
375        _sig_on
376        for i from 0 <= i < self._nrows:
377            self_row = self._matrix[i]
378            for j from 0 <= j < self._ncols:
379                mpz_lcm(d, d, mpq_denref(self_row[0]))
380                self_row = self_row + 1
381        _sig_off
382        return 0
383       
384    def _clear_denom(self):
385        """
386        INPUT:
387            self -- a matrix
388        OUTPUT:
389            D*self, D
390
391        The product is a matrix over ZZ
392        """
393        cdef Integer D
394        cdef Py_ssize_t i, j
395        cdef Matrix_integer_dense A
396        cdef mpq_t *self_row
397        cdef mpz_t *A_row
398        D = <Integer>Integer.__new__(Integer)
399        self.mpz_denom(D.value)
400        MZ = sage.matrix.matrix_space.MatrixSpace(ZZ, self._nrows, self._ncols, sparse=self.is_sparse())
401        A = Matrix_integer_dense.__new__(Matrix_integer_dense, MZ, 0, 0, 0)
402        _sig_on
403        for i from 0 <= i < self._nrows:
404            A_row = A._matrix[i]
405            self_row = self._matrix[i]
406            for j from 0 <= j < self._ncols:
407                mpz_init(A_row[0])
408                mpz_divexact(A_row[0], D.value, mpq_denref(self_row[0]))
409                mpz_mul(A_row[0], A_row[0], mpq_numref(self_row[0]))
410                A_row = A_row + 1
411                self_row = self_row + 1
412        _sig_off
413        return A, D
414
415    def _multiply_multi_modular(left, Matrix_rational_dense right):
416        """
417        Multiply this matrix by right using a multimodular algorithm
418        and return the result.
419       
420        EXAMPLES:
421            sage: a = MatrixSpace(QQ,3)(range(9))
422            sage: b = MatrixSpace(QQ,3)([1/n for n in range(1,10)])
423            sage: a._multiply_multi_modular(b)
424            [ 15/28   9/20   7/18]
425            [  33/7 117/40   20/9]
426            [249/28   27/5  73/18]
427            sage: a = MatrixSpace(QQ,10,5)(range(50))
428            sage: b = MatrixSpace(QQ,5,12)([1/n for n in range(1,61)])
429            sage: a._multiply_multi_modular(b) == a._multiply_classical(b)
430            True
431
432        """
433        cdef Matrix_integer_dense A, B, AB
434        cdef Matrix_rational_dense res
435        cdef Integer D
436        cdef mpz_t* AB_row,
437        cdef mpq_t* res_row
438        A_denom, B_denom
439        A, A_denom = left._clear_denom()
440        B, B_denom = right._clear_denom()
441        AB = A._multiply_multi_modular(B)
442        D = A_denom * B_denom
443        res = Matrix_rational_dense.__new__(Matrix_rational_dense, left.matrix_space(AB._nrows, AB._ncols), 0, 0, 0)
444        for i from 0 <= i < res._nrows:
445            AB_row = AB._matrix[i]
446            res_row = res._matrix[i]
447            for j from 0 <= j < res._ncols:
448                mpz_set(mpq_numref(res_row[0]), AB_row[0])
449                mpz_set(mpq_denref(res_row[0]), D.value)
450                mpq_canonicalize(res_row[0])
451                AB_row = AB_row + 1
452                res_row = res_row + 1
453        _sig_off
454        return res
455   
456    cdef int mpz_height(self, mpz_t height) except -1:
457        cdef mpz_t x, h
458        mpz_init(x)
459        mpz_init_set_si(h, 0)
460        cdef int i, j
461        _sig_on
462        for i from 0 <= i < self._nrows:
463            for j from 0 <= j < self._ncols:
464                mpq_get_num(x,self._matrix[i][j])
465                mpz_abs(x, x)
466                if mpz_cmp(h,x) < 0:
467                    mpz_set(h,x)
468                mpq_get_den(x,self._matrix[i][j])
469                mpz_abs(x, x)               
470                if mpz_cmp(h,x) < 0:
471                    mpz_set(h,x)
472        _sig_off
473        mpz_set(height, h)
474        mpz_clear(h)
475        mpz_clear(x)
476        return 0
477
478    def height(self):
479        """
480        Return the height of this matrix, which is the least common
481        multiple of all numerators and denominators of elements of
482        this matrix.
483
484        OUTPUT:
485            -- SAGE Integer
486
487        EXAMPLES:
488            sage: b = matrix(QQ,2,range(6)); b[0,0]=-5007/293; b
489            [-5007/293         1         2]
490            [        3         4         5]
491            sage: b.height()
492            5007
493        """
494        cdef Integer z
495        z = Integer.__new__(Integer)
496        self.mpz_height(z.value)
497        return z
498
499    cdef int _rescale(self, mpq_t a) except -1:
500        cdef int i, j
501        _sig_on
502        for i from 0 <= i < self._nrows:
503            for j from 0 <= j < self._ncols:
504                mpq_mul(self._matrix[i][j], self._matrix[i][j], a)
505        _sig_off       
506
507    def _adjoint(self):
508        """
509        Return the adjoint of this matrix.
510       
511        Assumes self is a square matrix (checked in adjoint).
512        """
513        return self.parent()(self._pari_().matadjoint().python())
514
515###########################
516
Note: See TracBrowser for help on using the repository browser.