1 #include "cado.h" // IWYU pragma: keep
2 // IWYU pragma: no_include <sys/param.h>
3 #include <cstring> // for memmove
4 #include <algorithm> // for min, max
5 #include <utility> // for move, swap
6 #include <gmp.h>
7 #include "lingen_matpoly_select.hpp" // for matpoly, matpoly::const_view_t
8 #include "mpfq_layer.h"
9 #include "omp_proxy.h"
10 #include "macros.h"
11 #include "lingen_matpoly.hpp"
12 #include "lingen_polymat.hpp"
13
14 matpoly::memory_pool_type matpoly::memory;
15
16 /* with the exception of matpoly_realloc, all functions here are exactly
17 * identical to those in lingen-polymat.c */
18 /* {{{ init/zero/clear interface for matpoly */
matpoly(abdst_field ab,unsigned int m,unsigned int n,int len)19 matpoly::matpoly(abdst_field ab, unsigned int m, unsigned int n, int len) : ab(ab), m(m), n(n), alloc(len) {
20 /* As a special case, we allow a pre-init state with m==n==len==0 */
21 /* Note that because we want to handle homogenous and non-homogenous
22 * cases the same way, we support matrices of size 0*n, so that is
23 * not technically a pre-init state */
24 if (!m && !n) {
25 ASSERT_ALWAYS(!len);
26 return;
27 }
28 if (alloc) {
29 if (data_alloc_size_in_bytes()) {
30 x = (abdst_vec) memory.alloc(data_alloc_size_in_bytes());
31 abvec_set_zero(ab, x, m*n*alloc);
32 } else {
33 x = NULL;
34 }
35 }
36 }
37
~matpoly()38 matpoly::~matpoly() {
39 if (x) {
40 memory.free(x, data_alloc_size_in_bytes());
41 // abvec_clear(ab, &(x), m*n*alloc);
42 }
43 }
matpoly(matpoly && a)44 matpoly::matpoly(matpoly && a)
45 : ab(a.ab), m(a.m), n(a.n), alloc(a.alloc)
46 {
47 size=a.size;
48 x=a.x;
49 a.x=NULL;
50 a.m=a.n=a.size=a.alloc=0;
51 a.ab=NULL;
52 }
operator =(matpoly && a)53 matpoly& matpoly::operator=(matpoly&& a)
54 {
55 if (x) {
56 memory.free(x, data_alloc_size_in_bytes());
57 // abvec_clear(ab, &(x), m*n*alloc);
58 }
59 ab = a.ab;
60 m = a.m;
61 n = a.n;
62 alloc = a.alloc;
63 size = a.size;
64 x=a.x;
65 a.x=NULL;
66 a.m=a.n=a.size=a.alloc=0;
67 a.ab=NULL;
68 return *this;
69 }
set(matpoly const & a)70 matpoly& matpoly::set(matpoly const& a)
71 {
72 if (x) {
73 memory.free(x, data_alloc_size_in_bytes());
74 // abvec_clear(ab, &(x), m*n*alloc);
75 }
76 ab = a.ab;
77 m = a.m;
78 n = a.n;
79 alloc = a.alloc;
80 size = a.size;
81 // abvec_init(ab, &(x), m*n*alloc);
82 x = (abdst_vec) memory.alloc(data_alloc_size_in_bytes());
83 abvec_set(ab, x, a.x, m*n*alloc);
84 return *this;
85 }
86
87 /* For structures in pre_init state, this is equivalent to an initial
88 * allocation.
89 *
90 * If data is grown, old data is kept, and the 'size' field is unchanged.
91 *
92 * If data is shrunk below the previous value of 'size', then 'size' is
93 * set to zero.
94 *
95 * The contents of the data area above 'size' on return is unspecified.
96 */
realloc(size_t newalloc)97 void matpoly::realloc(size_t newalloc) {
98 ASSERT_ALWAYS(size <= alloc);
99 size_t oldmem = data_alloc_size_in_bytes();
100 size_t newmem = data_alloc_size_in_bytes(newalloc);
101
102 /* zero out the newly added data */
103 if (newalloc > alloc) {
104 /* allocate new space, then inflate */
105 // abvec_reinit(ab, &(x), m*n*alloc, m*n*newalloc);
106 x = (abdst_vec) memory.realloc(x, oldmem, newmem);
107 abvec rhead = abvec_subvec(ab, x, m*n*alloc);
108 abvec whead = abvec_subvec(ab, x, m*n*newalloc);
109 if (size)
110 for(unsigned int i = m ; i-- ; ) {
111 for(unsigned int j = n ; j-- ; ) {
112 whead = abvec_subvec(ab, whead, -newalloc);
113 rhead = abvec_subvec(ab, rhead, -alloc);
114 abvec_set(ab, whead, rhead, size);
115 abvec_set_zero(ab,
116 abvec_subvec(ab, whead, alloc),
117 newalloc - alloc);
118 }
119 }
120 } else {
121 if (size > newalloc)
122 size = 0;
123 /* deflate, then free space */
124 ASSERT_ALWAYS(size <= newalloc);
125 abvec rhead = x;
126 abvec whead = x;
127 if (size)
128 for(unsigned int i = 0 ; i < m ; i++) {
129 for(unsigned int j = 0 ; j < n ; j++) {
130 abvec_set(ab, whead, rhead, size);
131 whead = abvec_subvec(ab, whead, newalloc);
132 rhead = abvec_subvec(ab, rhead, alloc);
133 }
134 }
135 // abvec_reinit(ab, &(x), m*n*alloc, m*n*newalloc);
136 x =(abdst_vec) memory.realloc(x, oldmem, newmem);
137 }
138 // if (!size) abvec_set_zero(ab, x, m*n*alloc);
139 alloc = newalloc;
140 }
141
get_true_nonzero_size() const142 size_t matpoly::get_true_nonzero_size() const
143 {
144 size_t lb = 0;
145 size_t ub = get_size();
146 for(unsigned int ij = 0 ; ij < m*n && lb < ub ; ij++) {
147 unsigned int i = ij / n;
148 unsigned int j = ij % n;
149 /* Find the last nonzero in the range [lb, ub[ */
150 for(unsigned int k = ub ; k > lb ; k--) {
151 if (!abis_zero(ab, coeff(i, j, k-1))) {
152 lb = k;
153 break;
154 }
155 }
156 }
157 return lb;
158 }
159
zero()160 void matpoly::zero() {
161 size = 0;
162 abvec_set_zero(ab, x, m*n*alloc);
163 }
164
set_constant_ui(unsigned long e)165 void matpoly::set_constant_ui(unsigned long e) {
166 ASSERT_ALWAYS(m == n);
167 size = 0;
168 if (alloc == 0 && e)
169 realloc(1);
170 zero();
171 if (!e) return;
172 size = 1;
173 for(unsigned int i = 0 ; i < m ; ++i)
174 abset_ui(ab, coeff(i, i, 0), e);
175 }
set_constant(absrc_elt e)176 void matpoly::set_constant(absrc_elt e) {
177 ASSERT_ALWAYS(m == n);
178 size = 0;
179 if (alloc == 0 && abcmp_ui(ab, e, 0) != 0)
180 realloc(1);
181 zero();
182 if (abcmp_ui(ab, e, 0) == 0) return;
183 size = 1;
184 for(unsigned int i = 0 ; i < m ; ++i)
185 abset(ab, coeff(i, i, 0), e);
186 }
187 /* }}} */
188
fill_random(unsigned int k0,unsigned int k1,gmp_randstate_t rstate)189 void matpoly::fill_random(unsigned int k0, unsigned int k1, gmp_randstate_t rstate)
190 {
191 ASSERT_ALWAYS(k1 <= alloc);
192 if (k0 == 0 && k1 == alloc) {
193 abvec_random(ab, x, m*n*k1, rstate);
194 } else if (k0 < k1) {
195 for(unsigned int i = 0 ; i < m ; i++) {
196 for(unsigned int j = 0 ; j < n ; j++) {
197 abvec_random(ab, part_head(i, j, k0), k1 - k0, rstate);
198 }
199 }
200 }
201 }
202
cmp(matpoly const & b) const203 int matpoly::cmp(matpoly const& b) const
204 {
205 ASSERT_ALWAYS(n == b.n);
206 ASSERT_ALWAYS(m == b.m);
207 if (size != b.size) {
208 return (size > b.size) - (b.size > size);
209 } else if (size == 0 && b.size == 0) {
210 /* This is for the "intermediary pre-init" state. size is zero,
211 * but a priori the dimension fields are ok. x might be NULL or
212 * not, it depends on the previous life of the object */
213 return 0;
214 } else if (size == alloc && b.size == b.alloc) {
215 return abvec_cmp(ab, x, b.x, m*n*size);
216 } else {
217 for(unsigned int i = 0 ; i < m ; i++) {
218 for(unsigned int j = 0 ; j < n ; j++) {
219 int r = abvec_cmp(ab, part(i, j), b.part(i, j), size);
220 if (r) return r;
221 }
222 }
223 return 0;
224 }
225 }
226
227 /* shift by a multiplication by x all coefficients of degree less than
228 * colsize in column j. This results in a new column of length colsize+1.
229 * Allocation must be sufficient for this length to fit. What happens to
230 * coefficients of degree larger than colsize in the result is unspecified.
231 *
232 * It is often relevant to "colsize++" right after this call, since the
233 * coefficient of degree colsize is well-defined on output
234 */
multiply_column_by_x(unsigned int j,unsigned int colsize)235 void matpoly::multiply_column_by_x(unsigned int j, unsigned int colsize)/*{{{*/
236 {
237 ASSERT_ALWAYS((colsize + 1) <= alloc);
238 for(unsigned int i = 0 ; i < m ; i++) {
239 memmove(part_head(i, j, 1), part(i, j), colsize * abvec_elt_stride(ab, 1));
240 abset_ui(ab, coeff(i, j, 0), 0);
241 }
242 }/*}}}*/
243
244 /* shift right (by a division by x) all coefficients of degree less than
245 * colsize in column j. This results in a new column of length colsize-1,
246 * with a zero coefficient shifted in at degree colsize-1.
247 *
248 * It is often relevant to "colsize--" right after this call.
249 */
divide_column_by_x(unsigned int j,unsigned int colsize)250 void matpoly::divide_column_by_x(unsigned int j, unsigned int colsize)/*{{{*/
251 {
252 ASSERT_ALWAYS(colsize <= alloc);
253 for(unsigned int i = 0 ; i < m ; i++) {
254 memmove(part(i, j), part_head(i, j, 1),
255 (colsize-1) * abvec_elt_stride(ab, 1));
256 abset_ui(ab, coeff(i, j, colsize-1), 0);
257 }
258 }/*}}}*/
259
truncate(matpoly const & src,unsigned int nsize)260 void matpoly::truncate(matpoly const & src, unsigned int nsize)/*{{{*/
261 {
262 ASSERT_ALWAYS(nsize <= src.alloc);
263 if (check_pre_init()) {
264 /* Need to call ctor... */
265 *this = matpoly(src.ab, src.m, src.n, nsize);
266 }
267 ASSERT_ALWAYS(m == src.m);
268 ASSERT_ALWAYS(n == src.n);
269 ASSERT_ALWAYS(nsize <= alloc);
270 ASSERT_ALWAYS(nsize <= src.size);
271 size = nsize;
272 if (this == &src) return;
273 /* XXX Much more cumbersome here than for polymat, of course */
274 for(unsigned int i = 0 ; i < src.m ; i++) {
275 for(unsigned int j = 0 ; j < src.n ; j++) {
276 abvec_set(ab, part(i, j), src.part(i, j), nsize);
277 }
278 }
279 }/*}}}*/
tail_is_zero(unsigned int size0)280 int matpoly::tail_is_zero(unsigned int size0)/*{{{*/
281 {
282 ASSERT_ALWAYS(size0 <= size);
283 for(unsigned int i = 0 ; i < m ; i++) {
284 for(unsigned int j = 0 ; j < n ; j++) {
285 if (!abvec_is_zero(ab, part_head(i, j, size0), size - size0))
286 return 0;
287 }
288 }
289 return 1;
290 }/*}}}*/
zero_pad(unsigned int nsize)291 void matpoly::zero_pad(unsigned int nsize)/*{{{*/
292 {
293 ASSERT_ALWAYS(nsize >= size);
294 if (check_pre_init() || nsize > alloc)
295 realloc(nsize);
296 for(unsigned int i = 0 ; i < m ; i++) {
297 for(unsigned int j = 0 ; j < n ; j++) {
298 abvec_set_zero(ab, part_head(i, j, size), nsize - size);
299 }
300 }
301 size = nsize;
302 }/*}}}*/
303
304
305 /* This takes coefficient ksrc of column jsrc, and copies it to
306 * coefficient kdst of column jdst
307 */
308 /* XXX compared to polymat, our diffferent stride has a consequence,
309 * clearly ! */
extract_column(unsigned int jdst,unsigned int kdst,matpoly const & src,unsigned int jsrc,unsigned int ksrc)310 void matpoly::extract_column( /*{{{*/
311 unsigned int jdst, unsigned int kdst,
312 matpoly const & src, unsigned int jsrc, unsigned int ksrc)
313 {
314 ASSERT_ALWAYS(m == src.m);
315 for(unsigned int i = 0 ; i < src.m ; i++)
316 abset(ab, coeff(i, jdst, kdst), src.coeff(i, jsrc, ksrc));
317 }/*}}}*/
318
319 #if 0
320 void matpoly::transpose_dumb(matpoly const & src) /*{{{*/
321 {
322 if (this == &src) {
323 matpoly tmp;
324 tmp.transpose_dumb(src);
325 *this = std::move(tmp);
326 return;
327 }
328 if (!check_pre_init()) {
329 *this = matpoly();
330 }
331 matpoly tmp(src.ab, src.n, src.m, src.size);
332 tmp.size = src.size;
333 for(unsigned int i = 0 ; i < src.m ; i++) {
334 for(unsigned int j = 0 ; j < src.n ; j++) {
335 for(unsigned int k = 0 ; k < src.size ; k++) {
336 abset(ab, tmp.coeff(j, i, k), src.coeff(i, j, k));
337 }
338 }
339 }
340 *this = std::move(tmp);
341 }/*}}}*/
342 #endif
343
zero_column(unsigned int jdst,unsigned int kdst)344 void matpoly::zero_column(unsigned int jdst, unsigned int kdst) /*{{{*/
345 {
346 for(unsigned int i = 0 ; i < m ; i++)
347 abset_zero(ab, coeff(i, jdst, kdst));
348 }/*}}}*/
349
350 #if 0
351 void matpoly::extract_row_fragment(/*{{{*/
352 unsigned int i1, unsigned int j1,
353 matpoly const & src, unsigned int i0, unsigned int j0,
354 unsigned int n)
355 {
356 ASSERT_ALWAYS(src.size <= alloc);
357 ASSERT_ALWAYS(size == src.size);
358 for(unsigned int k = 0 ; k < n ; k++)
359 abvec_set(ab, part(i1, j1 + k, 0), src.part(i0, j0 + k, 0), size);
360 }/*}}}*/
361 #endif
362
zero()363 void matpoly::view_t::zero() { /*{{{*/
364 unsigned int nrows = this->nrows();
365 unsigned int ncols = this->ncols();
366 #ifdef HAVE_OPENMP
367 #pragma omp parallel for collapse(2)
368 #endif
369 for(unsigned int i = 0 ; i < nrows ; i++) {
370 for(unsigned int j = 0 ; j < ncols ; j++) {
371 abvec_set_zero(M.ab, part(i,j), M.get_size());
372 }
373 }
374 }/*}}}*/
375
rshift(matpoly const & src,unsigned int k)376 void matpoly::rshift(matpoly const & src, unsigned int k)/*{{{*/
377 {
378 ASSERT_ALWAYS(k <= src.size);
379 unsigned int newsize = src.size - k;
380 if (check_pre_init()) {
381 *this = matpoly(src.ab, src.m, src.n, newsize);
382 }
383 ASSERT_ALWAYS(m == src.m);
384 ASSERT_ALWAYS(n == src.n);
385 ASSERT_ALWAYS(newsize <= alloc);
386 size = newsize;
387 for(unsigned int i = 0 ; i < src.m ; i++) {
388 for(unsigned int j = 0 ; j < src.n ; j++) {
389 abvec_set(ab, part(i, j), src.part_head(i, j, k), newsize);
390 }
391 }
392 }/*}}}*/
rshift(unsigned int k)393 void matpoly::rshift(unsigned int k)/*{{{*/
394 {
395 ASSERT_ALWAYS(k <= size);
396 unsigned int newsize = size - k;
397 ASSERT_ALWAYS(newsize <= alloc);
398 size = newsize;
399 for(unsigned int i = 0 ; i < m ; i++) {
400 for(unsigned int j = 0 ; j < n ; j++) {
401 /* can't use abvec_set because memcpy does not accept overlap */
402 for(unsigned s = 0 ; s < newsize ; s++) {
403 abset(ab, coeff(i, j, s), coeff(i, j, s + k));
404 }
405 }
406 }
407 }/*}}}*/
408
add(matpoly const & a,matpoly const & b)409 void matpoly::add(matpoly const & a, matpoly const & b)/*{{{*/
410 {
411 size_t csize = std::max(a.size, b.size);
412 ASSERT_ALWAYS(a.m == b.m);
413 ASSERT_ALWAYS(a.n == b.n);
414 if (check_pre_init()) {
415 *this = matpoly(a.ab, a.m, a.n, csize);
416 }
417 if (alloc < csize)
418 realloc(csize);
419 ASSERT_ALWAYS(alloc >= csize);
420
421 for(unsigned int i = 0 ; i < m ; ++i) {
422 for(unsigned int j = 0 ; j < n ; ++j) {
423 size_t s0 = std::min(a.size, b.size);
424 abvec_add(ab, part(i, j), a.part(i, j), b.part(i, j), s0);
425 if (a.size > s0)
426 abvec_set(ab, part_head(i, j, s0), a.part_head(i, j, s0), a.size - s0);
427 if (b.size > s0)
428 abvec_set(ab, part_head(i, j, s0), b.part_head(i, j, s0), b.size - s0);
429 }
430 }
431 size = csize;
432 }/*}}}*/
sub(matpoly const & a,matpoly const & b)433 void matpoly::sub(matpoly const & a, matpoly const & b)/*{{{*/
434 {
435 size_t csize = std::max(a.size, b.size);
436 ASSERT_ALWAYS(a.m == b.m);
437 ASSERT_ALWAYS(a.n == b.n);
438 if (check_pre_init()) {
439 *this = matpoly(a.ab, a.m, a.n, csize);
440 }
441 if (alloc < csize)
442 realloc(csize);
443 ASSERT_ALWAYS(alloc >= csize);
444
445 for(unsigned int i = 0 ; i < m ; ++i) {
446 for(unsigned int j = 0 ; j < n ; ++j) {
447 size_t s0 = std::min(a.size, b.size);
448 abvec_sub(ab, part(i, j), a.part(i, j), b.part(i, j), s0);
449 if (a.size > s0)
450 abvec_set(ab, part_head(i, j, s0), a.part_head(i, j, s0), a.size - s0);
451 if (b.size > s0)
452 abvec_neg(ab, part_head(i, j, s0), b.part_head(i, j, s0), b.size - s0);
453 }
454 }
455 size = csize;
456 }/*}}}*/
addmul(matpoly const & a,matpoly const & b)457 void matpoly::addmul(matpoly const & a, matpoly const & b)/*{{{*/
458 {
459 size_t csize = a.size + b.size; csize -= (csize > 0);
460 if (this == &a || this == &b) {
461 matpoly tc(a.ab, a.m, b.n, csize);
462 tc.addmul(a, b);
463 *this = std::move(tc);
464 return;
465 }
466 ASSERT_ALWAYS(a.n == b.m);
467 if (check_pre_init()) {
468 *this = matpoly(a.ab, a.m, b.n, csize);
469 }
470 ASSERT_ALWAYS(m == a.m);
471 ASSERT_ALWAYS(n == b.n);
472 ASSERT_ALWAYS(alloc >= csize);
473
474 if (a.size == 0 || b.size == 0)
475 return;
476
477 for(unsigned int i = 0 ; i < m ; ++i) {
478 for(unsigned int j = 0 ; j < n ; ++j) {
479 abvec_set_zero(ab, part_head(i, j, size), csize - size);
480 }
481 }
482
483 if (csize >= size)
484 zero_pad(csize);
485
486 matpoly::addmul(*this, a, b);
487
488 size = csize;
489 }/*}}}*/
490
copy(matpoly::view_t t,matpoly::const_view_t a)491 void matpoly::copy(matpoly::view_t t, matpoly::const_view_t a)/*{{{*/
492 {
493 unsigned int nrows = a.nrows();
494 unsigned int ncols = a.ncols();
495 ASSERT_ALWAYS(t.nrows() == nrows);
496 ASSERT_ALWAYS(t.ncols() == ncols);
497 #ifdef HAVE_OPENMP
498 unsigned int T = std::min((unsigned int) omp_get_max_threads(), nrows*ncols);
499 #pragma omp parallel num_threads(T)
500 #endif
501 {
502 #ifdef HAVE_OPENMP
503 #pragma omp for collapse(2)
504 #endif
505 for(unsigned int i = 0 ; i < nrows ; i++) {
506 for(unsigned int j = 0 ; j < ncols ; j++) {
507 ptr tij = t.part(i, j);
508 srcptr aij = a.part(i, j);
509 abvec_set(a.M.ab, tij, aij, a.M.get_size());
510 }
511 }
512 }
513 }/*}}}*/
514
addmul(matpoly::view_t t,matpoly::const_view_t t0,matpoly::const_view_t t1)515 void matpoly::addmul(matpoly::view_t t, matpoly::const_view_t t0, matpoly::const_view_t t1)/*{{{*/
516 {
517 unsigned int nrows = t.nrows();
518 unsigned int ncols = t.ncols();
519 unsigned int nadd = t0.ncols();
520 ASSERT_ALWAYS(&t.M != &t0.M);
521 ASSERT_ALWAYS(&t.M != &t1.M);
522 ASSERT_ALWAYS(&t0.M != &t1.M);
523 ASSERT_ALWAYS(t0.nrows() == nrows);
524 ASSERT_ALWAYS(t1.ncols() == ncols);
525 ASSERT_ALWAYS(t0.ncols() == nadd);
526 ASSERT_ALWAYS(t1.nrows() == nadd);
527 // ASSERT(t0.check());
528 // ASSERT(t1.check());
529 // ASSERT(t.check());
530 size_t csize = t0.M.size + t1.M.size; csize -= (csize > 0);
531 abdst_field ab = t0.M.ab;
532 if (t0.M.size == 0 || t1.M.size == 0)
533 return;
534 /* Attention: we don't want to count on t.M.size being set yet. */
535 #ifdef HAVE_OPENMP
536 unsigned int T = std::min((unsigned int) omp_get_max_threads(), nrows*ncols);
537 #pragma omp parallel num_threads(T)
538 #endif
539 {
540 abvec_ur tmp[2];
541 abvec_ur_init(ab, &tmp[0], csize);
542 abvec_ur_init(ab, &tmp[1], csize);
543
544 #ifdef HAVE_OPENMP
545 #pragma omp for collapse(2)
546 #endif
547 for(unsigned int i = 0 ; i < nrows ; i++) {
548 for(unsigned int j = 0 ; j < ncols ; j++) {
549 abvec_ur_set_vec(ab, tmp[1], t.part(i, j), csize);
550 for(unsigned int k = 0 ; k < nadd ; k++) {
551 abvec_conv_ur(ab, tmp[0],
552 t0.part(i, k), t0.M.size,
553 t1.part(k, j), t1.M.size);
554 abvec_ur_add(ab, tmp[1], tmp[1], tmp[0], csize);
555 }
556 abvec_reduce(ab, t.part(i, j), tmp[1], csize);
557 }
558 }
559 abvec_ur_clear(ab, &tmp[0], csize);
560 abvec_ur_clear(ab, &tmp[1], csize);
561 }
562 }/*}}}*/
addmp(matpoly::view_t t,matpoly::const_view_t t0,matpoly::const_view_t t1)563 void matpoly::addmp(matpoly::view_t t, matpoly::const_view_t t0, matpoly::const_view_t t1)/*{{{*/
564 {
565 unsigned int nrows = t.nrows();
566 unsigned int ncols = t.ncols();
567 unsigned int nadd = t0.ncols();
568 ASSERT_ALWAYS(&t.M != &t0.M);
569 ASSERT_ALWAYS(&t.M != &t1.M);
570 ASSERT_ALWAYS(&t0.M != &t1.M);
571 ASSERT_ALWAYS(t0.nrows() == nrows);
572 ASSERT_ALWAYS(t1.ncols() == ncols);
573 ASSERT_ALWAYS(t0.ncols() == nadd);
574 ASSERT_ALWAYS(t1.nrows() == nadd);
575 // ASSERT(t0.check());
576 // ASSERT(t1.check());
577 // ASSERT(t.check());
578 size_t fullsize = t0.M.size + t1.M.size; fullsize -= (fullsize > 0);
579 unsigned int nb = MAX(t0.M.size, t1.M.size) - MIN(t0.M.size, t1.M.size) + 1;
580 abdst_field ab = t0.M.ab;
581 if (t0.M.size == 0 || t1.M.size == 0)
582 return;
583 /* Attention: we don't want to count on t.M.size being set yet. */
584
585 /* We are going to make it completely stupid for a beginning. */
586 /* Note that the caching code that is used for larger sizes does a
587 * _real_ middle product.
588 */
589 #ifdef HAVE_OPENMP
590 unsigned int T = std::min((unsigned int) omp_get_max_threads(), nrows*ncols);
591 #pragma omp parallel num_threads(T)
592 #endif
593 {
594 abvec_ur tmp[2];
595 abvec_ur_init(ab, &tmp[0], fullsize);
596 abvec_ur_init(ab, &tmp[1], nb);
597
598 #ifdef HAVE_OPENMP
599 #pragma omp for collapse(2)
600 #endif
601 for(unsigned int i = 0 ; i < nrows ; i++) {
602 for(unsigned int j = 0 ; j < ncols ; j++) {
603 abvec_ur_set_vec(ab, tmp[1], t.part(i, j), nb);
604 for(unsigned int k = 0 ; k < nadd ; k++) {
605 abvec_conv_ur(ab, tmp[0],
606 t0.part(i, k), t0.M.size,
607 t1.part(k, j), t1.M.size);
608 abvec_ur_add(ab, tmp[1], tmp[1],
609 abvec_ur_subvec(ab, tmp[0], MIN(t0.M.size, t1.M.size) - 1), nb);
610 }
611 abvec_reduce(ab, t.part(i, j), tmp[1], nb);
612 }
613 }
614 abvec_ur_clear(ab, &tmp[0], fullsize);
615 abvec_ur_clear(ab, &tmp[1], nb);
616 }
617 }/*}}}*/
618
619
mul(matpoly const & a,matpoly const & b)620 matpoly matpoly::mul(matpoly const & a, matpoly const & b)/*{{{*/
621 {
622 size_t csize = a.size + b.size; csize -= (csize > 0);
623
624 matpoly tc(a.ab, a.m, b.n, csize);
625
626 ASSERT_ALWAYS(a.n == b.m);
627 tc.set_size(csize);
628 tc.zero();
629 tc.addmul(a, b);
630 return tc;
631 }/*}}}*/
632
addmp(matpoly const & a,matpoly const & c)633 void matpoly::addmp(matpoly const & a, matpoly const & c)/*{{{*/
634 {
635 size_t fullsize = a.size + c.size; fullsize -= (fullsize > 0);
636 unsigned int nb = MAX(a.size, c.size) - MIN(a.size, c.size) + 1;
637 ASSERT_ALWAYS(a.n == c.m);
638 if (check_pre_init()) {
639 *this = matpoly(a.ab, a.m, c.n, nb);
640 }
641 ASSERT_ALWAYS(m == a.m);
642 ASSERT_ALWAYS(n == c.n);
643 ASSERT_ALWAYS(alloc >= nb);
644
645 if (a.size == 0 || c.size == 0)
646 return;
647
648 if (nb >= size)
649 zero_pad(nb);
650
651 matpoly::addmp(*this, a, c);
652 }/*}}}*/
653
mp(matpoly const & a,matpoly const & c)654 matpoly matpoly::mp(matpoly const & a, matpoly const & c)/*{{{*/
655 {
656 unsigned int nb = MAX(a.size, c.size) - MIN(a.size, c.size) + 1;
657 ASSERT_ALWAYS(a.n == c.m);
658 matpoly b(a.ab, a.m, c.n, nb);
659 b.zero();
660 b.set_size(nb);
661 b.addmp(a, c);
662 return b;
663 }/*}}}*/
664
665
set_polymat(polymat const & src)666 void matpoly::set_polymat(polymat const & src)
667 {
668 *this = matpoly(src.ab, src.m, src.n, src.get_size());
669 set_size(src.get_size());
670
671 for(unsigned int i = 0 ; i < src.m ; i++) {
672 for(unsigned int j = 0 ; j < src.n ; j++) {
673 for(unsigned int k = 0 ; k < size ; k++) {
674 abset(ab, coeff(i, j, k), src.coeff(i, j, k));
675 }
676 }
677 }
678 }
679
coeff_is_zero(unsigned int k) const680 int matpoly::coeff_is_zero(unsigned int k) const
681 {
682 for(unsigned int j = 0; j < m; j++)
683 for(unsigned int i = 0 ; i < n ; i++)
684 if (!abis_zero(ab, coeff(i, j, k)))
685 return 0;
686 return 1;
687 }
coeff_set_zero(unsigned int k)688 void matpoly::coeff_set_zero(unsigned int k)
689 {
690 for(unsigned int j = 0; j < m; j++)
691 for(unsigned int i = 0 ; i < n ; i++)
692 abset_zero(ab, coeff(i, j, k));
693 }
694
truncate_and_rshift(unsigned int truncated_size,unsigned int shiftcount)695 matpoly matpoly::truncate_and_rshift(unsigned int truncated_size, unsigned int shiftcount)
696 {
697 matpoly other(ab, m, n, size - shiftcount);
698 other.rshift(*this, shiftcount);
699 truncate(*this, truncated_size);
700 shrink_to_fit();
701 std::swap(*this, other);
702 return other;
703 }
704