1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
|
diff --git a/src/sage/libs/linbox/fflas.pxd b/src/sage/libs/linbox/fflas.pxd
index ca23d71..03bbee9 100644
--- a/src/sage/libs/linbox/fflas.pxd
+++ b/src/sage/libs/linbox/fflas.pxd
@@ -20,19 +20,29 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "std":
iterator end()
size_t size()
- cdef cppclass list[T]:
- cppclass iterator:
- T operator*()
- iterator operator++()
- bint operator==(iterator)
- bint operator!=(iterator)
- void push_back(T&)
- void pop_front()
- T& front()
- iterator begin()
- iterator end()
+cdef extern from "givaro/givpoly1.h":
+ ## template < typename T, typename A=std::allocator<T> >
+ ## class givvector : public __GIV_STANDARD_VECTOR<T,A>
+ cdef cppclass givvector "Givaro::givvector" [T,ALLOCATOR=*]:
+ T& operator[](size_t i)
size_t size()
- void clear()
+
+ctypedef givvector[ModDoubleFieldElement] ModDoubleDensePolynomial
+ctypedef givvector[ModFloatFieldElement] ModFloatDensePolynomial
+
+cdef extern from "givaro/givpoly1.h":
+ ## template <class Domain, class StorageTag= Givaro::Dense>
+ ## class GivPolynomialRing : public Givaro::Poly1FactorDom< Domain,StorageTag>
+ cdef cppclass ModDoublePolynomialRing "Givaro::Poly1Dom<Givaro::Modular<double>, Givaro::Dense>":
+ ctypedef givvector[ModDoubleField] Element
+ ctypedef givvector[ModDoubleField] Polynomial
+ ModDoublePolynomialRing(ModDoubleField& F)
+ ## template <class Domain, class StorageTag= Givaro::Dense>
+ ## class GivPolynomialRing : public Givaro::Poly1FactorDom< Domain,StorageTag>
+ cdef cppclass ModFloatPolynomialRing "Givaro::Poly1Dom<Givaro::Modular<float>, Givaro::Dense>":
+ ctypedef givvector[ModFloatField] Element
+ ctypedef givvector[ModFloatField] Polynomial
+ ModFloatPolynomialRing(ModFloatField& F)
cdef extern from "fflas-ffpack/fflas-ffpack.h":
ctypedef enum fflas_trans_enum "FFLAS::FFLAS_TRANSPOSE":
@@ -105,13 +115,12 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h":
size_t nr, size_t foo, size_t r,
ModDoubleFieldElement* matrix, size_t nc, size_t* Q)
- void ModDouble_MinPoly "FFPACK::MinPoly" ( ModDoubleField F,
+ void ModDouble_MinPoly "FFPACK::MinPoly" ( ModDoubleField& F,
vector[ModDoubleFieldElement] minP, size_t N,
- ModDoubleFieldElement* A, size_t lda,
- ModDoubleFieldElement* X, size_t ldx, size_t* P)
+ ModDoubleFieldElement*A, size_t lda)
- void ModDouble_CharPoly "FFPACK::CharPoly" ( ModDoubleField F,
- list[vector[ModDoubleFieldElement]] charp, size_t N,
+ void ModDouble_CharPoly "FFPACK::CharPoly" ( ModDoublePolynomialRing& R,
+ ModDoubleDensePolynomial& charp, size_t N,
ModDoubleFieldElement* A, size_t lda)
# float
@@ -142,9 +151,8 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h":
void ModFloat_MinPoly "FFPACK::MinPoly" ( ModFloatField F,
vector[ModFloatFieldElement] minP, size_t N,
- ModFloatFieldElement* A, size_t lda,
- ModFloatFieldElement* X, size_t ldx, size_t* P)
+ ModFloatFieldElement* A, size_t lda)
- void ModFloat_CharPoly "FFPACK::CharPoly" ( ModFloatField F,
- list[vector[ModFloatFieldElement]] charp, size_t N,
+ void ModFloat_CharPoly "FFPACK::CharPoly" ( ModFloatPolynomialRing& F,
+ ModFloatDensePolynomial& charp, size_t N,
ModFloatFieldElement* A, size_t lda )
diff --git a/src/sage/libs/linbox/linbox_flint_interface.pyx b/src/sage/libs/linbox/linbox_flint_interface.pyx
index 02a6151..a00d303 100644
--- a/src/sage/libs/linbox/linbox_flint_interface.pyx
+++ b/src/sage/libs/linbox/linbox_flint_interface.pyx
@@ -76,19 +76,23 @@ cdef extern from "givaro/zring.h":
ctypedef GivaroInteger Element
-cdef extern from "givaro/givpoly1.h":
- ## template < typename T, typename A=std::allocator<T> >
- ## class givvector : public __GIV_STANDARD_VECTOR<T,A>
- cdef cppclass givvector "Givaro::givvector" [T,ALLOCATOR=*]:
- T& operator[](size_t i)
+cdef extern from "linbox/polynomial/dense-polynomial.h":
+ ## template<class Field>
+ ## class DensePolynomial : public Givaro::Poly1FactorDom<Field, Givaro::Dense>::Element
+ cdef cppclass LinBoxIntegerDensePolynomial "LinBox::DensePolynomial<Givaro::ZRing<Givaro::Integer> >":
+ ctypedef GivaroIntegerRing BaseRing
+ ctypedef GivaroInteger BaseRingElement
+ LinBoxIntegerDensePolynomial(BaseRing &F)
+ LinBoxIntegerDensePolynomial(BaseRing &F, size_t s)
+ BaseRingElement& operator[](size_t i)
size_t size()
-cdef extern from "linbox/ring/givaro-polynomial.h":
- ## template <class Domain, class StorageTag= Givaro::Dense>
- ## class GivPolynomialRing : public Givaro::Poly1FactorDom< Domain,StorageTag>
- cdef cppclass LinBoxIntegerPolynomialRing "LinBox::GivPolynomialRing<Givaro::ZRing<Givaro::Integer>, Givaro::Dense>":
- ctypedef givvector[GivaroInteger] Element
- ctypedef givvector[GivaroInteger] Polynomial
+cdef extern from "linbox/ring/polynomial-ring.h":
+ ## template <class BaseRing, class StorageTag= Givaro::Dense>
+ ## class PolynomialRing : public Givaro::Poly1FactorDom<BaseRing,StorageTag>
+ cdef cppclass LinBoxIntegerPolynomialRing "LinBox::PolynomialRing<Givaro::ZRing<Givaro::Integer>, Givaro::Dense>":
+ ctypedef LinBoxIntegerDensePolynomial Element
+ ctypedef LinBoxIntegerDensePolynomial Polynomial
cdef extern from "linbox/matrix/matrix-domain.h":
## template <class Field_ >
@@ -190,40 +194,32 @@ cdef void linbox_fmpz_mat_mul(fmpz_mat_t C, fmpz_mat_t A, fmpz_mat_t B):
cdef void linbox_fmpz_mat_charpoly(fmpz_poly_t cp, fmpz_mat_t A):
cdef GivaroIntegerRing ZZ
cdef LinBoxIntegerDenseMatrix * LBA
- cdef LinBoxIntegerPolynomialRing.Element m_A
-
- # FIXME: bug in LinBox
- # see https://github.com/linbox-team/linbox/issues/51
- if fmpz_mat_nrows(A) == 0:
- fmpz_poly_one(cp)
- return
+ cdef LinBoxIntegerDensePolynomial * m_A
LBA = new LinBoxIntegerDenseMatrix(ZZ, fmpz_mat_nrows(A), fmpz_mat_ncols(A))
fmpz_mat_get_linbox(LBA[0], A)
- LinBoxIntegerDense_charpoly(m_A, LBA[0])
- fmpz_poly_set_linbox(cp, m_A)
+ m_A = new LinBoxIntegerDensePolynomial(ZZ, fmpz_mat_nrows(A))
+ LinBoxIntegerDense_charpoly(m_A[0], LBA[0])
+ fmpz_poly_set_linbox(cp, m_A[0])
del LBA
+ del m_A
# set mp to the minimal polynomial of A
cdef void linbox_fmpz_mat_minpoly(fmpz_poly_t mp, fmpz_mat_t A):
cdef GivaroIntegerRing ZZ
cdef LinBoxIntegerDenseMatrix * LBA
- cdef LinBoxIntegerPolynomialRing.Element m_A
-
- # FIXME: bug in LinBox
- # see https://github.com/linbox-team/linbox/issues/51
- if fmpz_mat_nrows(A) == 0:
- fmpz_poly_one(mp)
- return
+ cdef LinBoxIntegerDensePolynomial * m_A
LBA = new LinBoxIntegerDenseMatrix(ZZ, fmpz_mat_nrows(A), fmpz_mat_ncols(A))
+ m_A = new LinBoxIntegerDensePolynomial(ZZ)
fmpz_mat_get_linbox(LBA[0], A)
- LinBoxIntegerDense_minpoly(m_A, LBA[0])
- fmpz_poly_set_linbox(mp, m_A)
+ LinBoxIntegerDense_minpoly(m_A[0], LBA[0])
+ fmpz_poly_set_linbox(mp, m_A[0])
del LBA
+ del m_A
# return the rank of A
diff --git a/src/sage/matrix/matrix_modn_dense_double.pyx b/src/sage/matrix/matrix_modn_dense_double.pyx
index 15c35cb..b1f6154 100644
--- a/src/sage/matrix/matrix_modn_dense_double.pyx
+++ b/src/sage/matrix/matrix_modn_dense_double.pyx
@@ -27,7 +27,9 @@ from sage.libs.linbox.fflas cimport ModDouble_fgemm as Mod_fgemm, ModDouble_fgem
ModDoubleRank as ModRank, ModDouble_echelon as Mod_echelon, \
ModDouble_applyp as Mod_applyp, \
ModDouble_MinPoly as Mod_MinPoly, \
- ModDouble_CharPoly as Mod_CharPoly
+ ModDouble_CharPoly as Mod_CharPoly, \
+ ModDoublePolynomialRing as ModDensePolyRing,\
+ ModDoubleDensePolynomial as ModDensePoly
# Limit for LinBox Modular<double>
MAX_MODULUS = 2**23
diff --git a/src/sage/matrix/matrix_modn_dense_float.pyx b/src/sage/matrix/matrix_modn_dense_float.pyx
index 61457ce..eda8f0b 100644
--- a/src/sage/matrix/matrix_modn_dense_float.pyx
+++ b/src/sage/matrix/matrix_modn_dense_float.pyx
@@ -25,7 +25,10 @@ from sage.libs.linbox.fflas cimport ModFloat_fgemm as Mod_fgemm, ModFloat_fgemv
ModFloatRank as ModRank, ModFloat_echelon as Mod_echelon, \
ModFloat_applyp as Mod_applyp, \
ModFloat_MinPoly as Mod_MinPoly, \
- ModFloat_CharPoly as Mod_CharPoly
+ ModFloat_CharPoly as Mod_CharPoly, \
+ ModFloatPolynomialRing as ModDensePolyRing,\
+ ModFloatDensePolynomial as ModDensePoly
+
# LinBox supports up to 2^11 using float but that's double dog slow,
# so we pick a smaller value for crossover
diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi
index 3da3e31..819ef2d 100644
--- a/src/sage/matrix/matrix_modn_dense_template.pxi
+++ b/src/sage/matrix/matrix_modn_dense_template.pxi
@@ -306,23 +306,15 @@ cdef inline linbox_minpoly(celement modulus, Py_ssize_t nrows, celement* entries
cdef Py_ssize_t i
cdef ModField *F = new ModField(<long>modulus)
cdef vector[ModFieldElement] *minP = new vector[ModFieldElement]()
- cdef ModFieldElement *X = <ModFieldElement*>check_allocarray(nrows * (nrows+1), sizeof(ModFieldElement))
- cdef size_t *P = <size_t*>check_allocarray(nrows, sizeof(size_t))
-
- cdef celement *cpy = linbox_copy(modulus, entries, nrows, nrows)
if nrows*nrows > 1000: sig_on()
- Mod_MinPoly(F[0], minP[0], nrows, <ModFieldElement*>cpy, nrows, X, nrows, P)
+ Mod_MinPoly(F[0], minP[0], nrows, <ModFieldElement*>entries, nrows)
if nrows*nrows > 1000: sig_off()
- sig_free(cpy)
-
l = []
for i in range(minP.size()):
l.append( <celement>minP.at(i) )
- sig_free(P)
- sig_free(X)
del F
return l
@@ -332,27 +324,23 @@ cdef inline linbox_charpoly(celement modulus, Py_ssize_t nrows, celement* entrie
"""
cdef Py_ssize_t i
cdef ModField *F = new ModField(<long>modulus)
- cdef std_list[vector[ModFieldElement]] P_list
- P_list.clear()
+ cdef ModDensePolyRing * R = new ModDensePolyRing(F[0])
+ cdef ModDensePoly P
cdef celement *cpy = linbox_copy(modulus, entries, nrows, nrows)
if nrows*nrows > 1000: sig_on()
- Mod_CharPoly(F[0], P_list, nrows, <ModFieldElement*>cpy, nrows)
+ Mod_CharPoly(R[0], P, nrows, <ModFieldElement*>cpy, nrows)
if nrows*nrows > 1000: sig_off()
sig_free(cpy)
- cdef vector[ModFieldElement] tmp
l = []
- while P_list.size():
- l.append([])
- tmp = P_list.front()
- for i in range(tmp.size()):
- l[-1].append(<celement>tmp.at(i))
- P_list.pop_front()
+ for i in range(P.size()):
+ l.append(<celement>P[i])
del F
+ del R
return l
@@ -1731,9 +1719,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
R = self._base_ring[var]
# call linbox for charpoly
v = linbox_charpoly(self.p, self._nrows, self._entries)
- r = R(1)
- for e in v:
- r *= R(e)
+ r = R(v)
return r
def echelonize(self, algorithm="linbox", **kwds):
|