source: sage/libs/linbox/linbox.pyx @ 4367:ce011a5e756c

Revision 4367:ce011a5e756c, 9.7 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

new ver
Line 
1## NOTE: The _sig_on/_sig_off stuff can't go in here -- it has to be in the
2## code that calls these functions.  Otherwise strangely objects get left
3## in an incorrect state.
4
5from sage.rings.integer cimport Integer
6from sage.misc.misc import verbose, get_verbose, cputime, UNAME
7
8##########################################################################
9## Dense matrices modulo p
10##########################################################################
11
12# LinBox bugs to address:
13#  * charpoly and minpoly don't work randomly
14
15cdef extern from "linbox_wrap.h":
16    void linbox_modn_dense_delete_array(mod_int *f)
17
18    int linbox_modn_dense_echelonize(unsigned long modulus, 
19                              mod_int **matrix, size_t nrows, size_t ncols)
20    void linbox_modn_dense_minpoly(unsigned long modulus, mod_int **mp, size_t* degree, size_t n,
21                                   mod_int **matrix, int do_minpoly)
22
23    int  linbox_modn_dense_matrix_matrix_multiply(unsigned long modulus, mod_int **ans,
24                                                  mod_int **A, mod_int **B,
25                                                  size_t A_nr, size_t A_nc,
26                                                  size_t B_nr, size_t B_nc)
27
28    int linbox_modn_dense_rank(unsigned long modulus, 
29                               mod_int** matrix, size_t nrows, size_t ncols)
30
31
32cdef class Linbox_modn_dense:
33    cdef set(self, mod_int n, mod_int** matrix,
34             size_t nrows, size_t ncols):
35        self.n = n
36        self.nrows = nrows
37        self.ncols = ncols
38        self.matrix = matrix
39
40    cdef int echelonize(self):
41        cdef int r
42        r = linbox_modn_dense_echelonize(self.n, self.matrix,
43                                         self.nrows, self.ncols)
44        return r
45   
46    def minpoly(self):
47        return self._poly(True)
48
49    def charpoly(self):
50        return self._poly(False)
51
52    def _poly(self, minpoly):
53        """
54        INPUT:
55            as given
56           
57        OUTPUT:
58            coefficients of charpoly or minpoly as a Python list
59        """
60        cdef mod_int *f
61        cdef size_t degree
62        linbox_modn_dense_minpoly(self.n, &f, &degree,
63                                  self.nrows, self.matrix,
64                                  minpoly)
65        v = []
66        cdef Py_ssize_t i
67        for i from 0 <= i <= degree:
68            v.append(f[i])
69        linbox_modn_dense_delete_array(f)
70        return v
71       
72    cdef matrix_matrix_multiply(self,
73                                mod_int **ans,
74                                mod_int **B,
75                                size_t B_nr, size_t B_nc):
76        cdef int e
77        e = linbox_modn_dense_matrix_matrix_multiply(self.n, ans,
78                                                     self.matrix,  B,
79                                                     self.nrows, self.ncols,
80                                                     B_nr, B_nc)
81        if e:
82            raise RuntimeError, "error doing matrix matrix multiply modn using linbox"
83   
84
85    cdef unsigned long rank(self) except -1:
86        cdef unsigned long r
87        r = linbox_modn_dense_rank(self.n,   self.matrix, self.nrows, self.ncols)
88        return r
89
90
91
92
93
94
95##########################################################################
96## Dense matrices GF(2)
97##########################################################################
98
99## cdef extern from "linbox_wrap.h":
100##     ctypedef struct packedmatrix
101
102##     cdef int linbox_mod2_dense_echelonize(packedmatrix *m)
103   
104##     cdef void linbox_mod2_dense_minpoly(mod_int **mp, size_t* degree, packedmatrix *matrix, int do_minpoly)
105   
106##     cdef int linbox_mod2_dense_matrix_matrix_multiply(packedmatrix *ans, packedmatrix *A, packedmatrix *B)
107
108##     cdef int linbox_mod2_dense_rank(packedmatrix *m)
109
110## cdef class Linbox_mod2_dense:
111##     cdef set(self, packedmatrix *matrix):
112##         self.matrix = matrix
113
114##     cdef int echelonize(self):
115##         cdef int r
116##         r = linbox_mod2_dense_echelonize(self.matrix)
117##         return r
118   
119##     def minpoly(self):
120##         return self._poly(True)
121
122##     def charpoly(self):
123##         return self._poly(False)
124
125##     def _poly(self, minpoly):
126##         """
127##         INPUT:
128##             as given
129           
130##         OUTPUT:
131##             coefficients of charpoly or minpoly as a Python list
132##         """
133##         cdef mod_int *f
134##         cdef size_t degree
135##         linbox_mod2_dense_minpoly(&f, &degree, self.matrix, minpoly)
136
137##         v = []
138##         cdef Py_ssize_t i
139##         for i from 0 <= i <= degree:
140##             v.append(f[i])
141##         linbox_modn_dense_delete_array(f)
142##         return v
143       
144##     cdef matrix_matrix_multiply(self,
145##                                 packedmatrix *ans,
146##                                 packedmatrix *B):
147##         cdef int e
148
149##         e = linbox_mod2_dense_matrix_matrix_multiply(ans, self.matrix,  B)
150##         if e:
151##             raise RuntimeError, "error doing matrix matrix multiply mod2 using linbox"
152
153##     cdef unsigned long rank(self) except -1:
154##         cdef unsigned long r
155##         r = linbox_mod2_dense_rank(self.matrix)
156##         return r
157       
158
159##########################################################################
160## Sparse matices modulo p.
161##########################################################################
162cdef class Linbox_modn_sparse:
163    pass
164
165
166##########################################################################
167## Sparse matrices over ZZ
168##########################################################################
169
170
171##########################################################################
172## Dense matrices over ZZ
173##########################################################################
174
175cdef extern from "linbox_wrap.h":
176    void linbox_integer_dense_minpoly_hacked(mpz_t* *minpoly, size_t* degree,
177                                      size_t n, mpz_t** matrix, int do_minpoly)
178
179    void linbox_integer_dense_minpoly(mpz_t* *minpoly, size_t* degree,
180                                      size_t n, mpz_t** matrix)
181
182    void linbox_integer_dense_charpoly(mpz_t* *charpoly, size_t* degree,
183                                       size_t n, mpz_t** matrix)
184
185    void linbox_integer_dense_delete_array(mpz_t* f)
186
187    int linbox_integer_dense_matrix_matrix_multiply(mpz_t** ans, mpz_t **A, mpz_t **B,
188                                      size_t A_nr, size_t A_nc, size_t B_nr, size_t B_nc)
189
190    unsigned long linbox_integer_dense_rank(mpz_t** matrix, size_t nrows,
191                                            size_t ncols)
192
193    void linbox_integer_dense_det(mpz_t ans, mpz_t** matrix,
194                             size_t nrows, size_t ncols)
195   
196    void linbox_integer_dense_smithform(mpz_t **v,
197                                        mpz_t **matrix,
198                                        size_t nrows, size_t ncols)
199   
200cdef class Linbox_integer_dense:
201    cdef set(self, mpz_t** matrix, size_t nrows, size_t ncols):
202        self.nrows = nrows
203        self.ncols = ncols
204        self.matrix = matrix
205
206    def minpoly(self):
207        return self._poly(True)
208
209    def charpoly(self):
210        return self._poly(False)
211   
212    def _poly(self, do_minpoly):
213        """
214        INPUT:
215            as given
216           
217        OUTPUT:
218            coefficients of charpoly or minpoly as a Python list
219        """
220        cdef mpz_t* poly
221        cdef size_t degree
222        if self.nrows % 4 == 0 and UNAME == "Darwin":
223            verbose("using hack to get around bug in linbox on OS X since n is divisible by 4")
224            if do_minpoly:
225                linbox_integer_dense_minpoly_hacked(&poly, &degree, self.nrows, self.matrix, 1)
226            else:
227                linbox_integer_dense_minpoly_hacked(&poly, &degree, self.nrows, self.matrix, 0)
228        else:
229            verbose("using linbox poly comp")
230            if do_minpoly:
231                linbox_integer_dense_minpoly(&poly, &degree, self.nrows, self.matrix)
232            else:
233                linbox_integer_dense_charpoly(&poly, &degree, self.nrows, self.matrix)
234            verbose("computed poly -- now converting back to SAGE")
235           
236        v = []
237        cdef Integer k
238        cdef size_t n
239        for n from 0 <= n <= degree:
240            k = Integer()
241            mpz_set(k.value, poly[n])
242            mpz_clear(poly[n])
243            v.append(k)
244        linbox_integer_dense_delete_array(poly)
245        return v
246       
247    cdef matrix_matrix_multiply(self,
248                                mpz_t **ans,
249                                mpz_t **B,
250                                size_t B_nr, size_t B_nc):
251        cdef int e
252        e = linbox_integer_dense_matrix_matrix_multiply(ans,
253                                                        self.matrix,  B,
254                                                        self.nrows, self.ncols,
255                                                        B_nr, B_nc)
256        if e:
257            raise RuntimeError, "error doing matrix matrix multiply over ZZ using linbox"
258   
259
260    cdef unsigned long rank(self) except -1:
261        return linbox_integer_dense_rank(self.matrix, self.nrows, self.ncols)
262
263    def det(self):
264        cdef Integer z
265        z = Integer()
266        linbox_integer_dense_det(z.value, self.matrix, self.nrows, self.ncols)
267        return z
268
269    def smithform(self):
270        raise NotImplementedError
271        #cdef mpz_t* v
272        #linbox_integer_dense_smithform(&v, self.matrix, self.nrows, self.ncols)
273        #z = []
274        #cdef Integer k
275        #cdef size_t n
276        #for n from 0 <= n < self.ncols:
277        #    k = Integer()
278        #    mpz_set(k.value, v[n])
279        #    mpz_clear(v[n])
280        #    z.append(k)
281        #linbox_integer_dense_delete_array(v)
282        #return z
283       
Note: See TracBrowser for help on using the repository browser.