1 /*
2 Copyright (C) 2017 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <gmp.h>
13 #include <stdlib.h>
14 #include "flint.h"
15 #include "fmpz.h"
16 #include "fmpz_mpoly.h"
17 #include "longlong.h"
18
19 /*
20 As for divrem_monagan_pearce1 except that an array of divisor polynomials is
21 passed and an array of quotient polynomials is returned. These are not in
22 low level format.
23 */
_fmpz_mpoly_divrem_ideal_monagan_pearce1(fmpz_mpoly_struct ** polyq,fmpz ** polyr,ulong ** expr,slong * allocr,const fmpz * poly2,const ulong * exp2,slong len2,fmpz_mpoly_struct * const * poly3,ulong * const * exp3,slong len,slong bits,const fmpz_mpoly_ctx_t ctx,ulong maskhi)24 slong _fmpz_mpoly_divrem_ideal_monagan_pearce1(fmpz_mpoly_struct ** polyq,
25 fmpz ** polyr, ulong ** expr, slong * allocr, const fmpz * poly2,
26 const ulong * exp2, slong len2, fmpz_mpoly_struct * const * poly3,
27 ulong * const * exp3, slong len, slong bits, const fmpz_mpoly_ctx_t ctx,
28 ulong maskhi)
29 {
30 slong i, j, p, l, w;
31 slong next_loc;
32 slong * store, * store_base;
33 slong len3;
34 slong heap_len = 2; /* heap zero index unused */
35 mpoly_heap1_s * heap;
36 mpoly_nheap_t ** chains;
37 slong ** hinds;
38 mpoly_nheap_t * x;
39 fmpz * p2 = *polyr;
40 ulong * e2 = *expr;
41 ulong exp, texp;
42 ulong c[3]; /* for accumulating coefficients */
43 ulong mask;
44 ulong * ub;
45 slong * k, * s;
46 fmpz_t qc, q;
47 fmpz * mb;
48 int small;
49 slong bits2, bits3;
50 int d1, d2, div_flag;
51 TMP_INIT;
52
53 TMP_START;
54
55
56 fmpz_init(q);
57 fmpz_init(qc);
58
59 bits2 = _fmpz_vec_max_bits(poly2, len2);
60
61 chains = (mpoly_nheap_t **) TMP_ALLOC(len*sizeof(mpoly_nheap_t *));
62 hinds = (slong **) TMP_ALLOC(len*sizeof(mpoly_heap_t *));
63
64 bits3 = 0;
65 len3 = 0;
66 for (w = 0; w < len; w++)
67 {
68 chains[w] = (mpoly_nheap_t *) TMP_ALLOC((poly3[w]->length)*sizeof(mpoly_nheap_t));
69 hinds[w] = (slong *) TMP_ALLOC((poly3[w]->length)*sizeof(slong));
70 bits3 = FLINT_MAX(bits3, FLINT_ABS(fmpz_mpoly_max_bits(poly3[w])));
71 len3 += poly3[w]->length;
72 for (i = 0; i < poly3[w]->length; i++)
73 hinds[w][i] = 1;
74 }
75
76 /* allow one bit for sign, one bit for subtraction */
77 small = FLINT_ABS(bits2) <= (FLINT_ABS(bits3) + FLINT_BIT_COUNT(len3) + FLINT_BITS - 2)
78 && FLINT_ABS(bits3) <= FLINT_BITS - 2;
79
80 next_loc = len3 + 4; /* something bigger than heap can ever be */
81 heap = (mpoly_heap1_s *) TMP_ALLOC((len3 + 1)*sizeof(mpoly_heap1_s));
82 store = store_base = (slong *) TMP_ALLOC(3*len3*sizeof(mpoly_nheap_t *));
83
84 k = (slong *) TMP_ALLOC(len*sizeof(slong));
85 s = (slong *) TMP_ALLOC(len*sizeof(slong));
86 ub = (ulong *) TMP_ALLOC(len*sizeof(ulong));
87 mb = (fmpz * ) TMP_ALLOC(len*sizeof(fmpz));
88
89 mask = 0;
90 for (i = 0; i < FLINT_BITS/bits; i++)
91 mask = (mask << bits) + (UWORD(1) << (bits - 1));
92
93 for (w = 0; w < len; w++)
94 {
95 k[w] = -WORD(1);
96 s[w] = poly3[w]->length;
97 }
98 l = -WORD(1);
99
100 x = chains[0] + 0;
101 x->i = -WORD(1);
102 x->j = 0;
103 x->p = -WORD(1);
104 x->next = NULL;
105 HEAP_ASSIGN(heap[1], exp2[0], x);
106
107 for (i = 0; i < len; i++)
108 {
109 fmpz_init(mb + i);
110 fmpz_neg(mb + i, poly3[i]->coeffs);
111 ub[i] = ((ulong) FLINT_ABS(mb[i])) >> 1; /* abs(poly3[0])/2 */
112 }
113
114 while (heap_len > 1)
115 {
116 exp = heap[1].exp;
117 if (mpoly_monomial_overflows1(exp, mask))
118 goto exp_overflow;
119
120 c[0] = c[1] = c[2] = 0;
121 fmpz_zero(qc);
122 while (heap_len > 1 && heap[1].exp == exp)
123 {
124 x = _mpoly_heap_pop1(heap, &heap_len, maskhi);
125 do
126 {
127 *store++ = x->i;
128 *store++ = x->j;
129 *store++ = x->p;
130 if (x->i != -WORD(1))
131 hinds[x->p][x->i] |= WORD(1);
132
133 if (small)
134 {
135 if (x->i == -WORD(1))
136 _fmpz_mpoly_sub_uiuiui_fmpz(c, poly2 + x->j);
137 else
138 _fmpz_mpoly_addmul_uiuiui_fmpz(c,
139 poly3[x->p]->coeffs[x->i], polyq[x->p]->coeffs[x->j]);
140 } else
141 {
142 if (x->i == -WORD(1))
143 fmpz_sub(qc, qc, poly2 + x->j);
144 else
145 fmpz_addmul(qc, poly3[x->p]->coeffs + x->i,
146 polyq[x->p]->coeffs + x->j);
147 }
148 } while ((x = x->next) != NULL);
149 }
150
151 while (store > store_base)
152 {
153 p = *--store;
154 j = *--store;
155 i = *--store;
156 if (i == -WORD(1))
157 {
158 if (j + 1 < len2)
159 {
160 x = chains[0] + 0;
161 x->i = -WORD(1);
162 x->j = j + 1;
163 x->p = p;
164 x->next = NULL;
165 _mpoly_heap_insert1(heap, exp2[x->j], x, &next_loc, &heap_len, maskhi);
166 }
167 } else
168 {
169 if ( (i + 1 < poly3[p]->length)
170 && (hinds[p][i + 1] == 2*j + 1)
171 )
172 {
173 x = chains[p] + i + 1;
174 x->i = i + 1;
175 x->j = j;
176 x->p = p;
177 x->next = NULL;
178 hinds[p][x->i] = 2*(x->j + 1) + 0;
179 _mpoly_heap_insert1(heap, exp3[x->p][x->i] +
180 polyq[x->p]->exps[x->j], x, &next_loc, &heap_len, maskhi);
181 }
182 if (j == k[p])
183 {
184 s[p]++;
185 } else if ( ((hinds[p][i] & 1) == 1)
186 && ((i == 1) || (hinds[p][i - 1] >= 2*(j + 2) + 1))
187 )
188 {
189 x = chains[p] + i;
190 x->i = i;
191 x->j = j + 1;
192 x->p = p;
193 x->next = NULL;
194 hinds[p][x->i] = 2*(x->j + 1) + 0;
195 _mpoly_heap_insert1(heap, exp3[x->p][x->i] +
196 polyq[x->p]->exps[x->j], x, &next_loc, &heap_len, maskhi);
197 }
198 }
199 }
200
201 if ((small && (c[2] != 0 || c[1] != 0 || c[0] != 0)) ||
202 (!small && !fmpz_is_zero(qc)))
203 {
204 div_flag = 0;
205 for (w = 0; w < len; w++)
206 {
207 d1 = mpoly_monomial_divides1(&texp, exp, exp3[w][0], mask);
208 if (d1)
209 {
210 d2 = 0;
211 if (small)
212 {
213 ulong d[3];
214 if (0 > (slong) c[2])
215 mpn_neg(d, c, 3);
216 else
217 flint_mpn_copyi(d, c, 3);
218
219 if (d[2] != 0 || ub[w] <= d[1] ||
220 (ub[w] == 0 && 0 > (slong) d[0])) /* quotient not a small */
221 {
222 fmpz_set_signed_uiuiui(qc, c[2], c[1], c[0]);
223 small = 0;
224 } else /* quotient fits a small */
225 {
226 slong r1;
227 slong tq;
228 sdiv_qrnnd(tq, r1, c[1], c[0], mb[w]);
229 if (tq > COEFF_MAX || tq < COEFF_MIN)
230 {
231 small = 0;
232 fmpz_set_signed_uiuiui(qc, c[2], c[1], c[0]);
233 } else
234 {
235 div_flag = (r1 == 0);
236 d2 = tq != 0;
237 if (d2)
238 {
239 k[w]++;
240 fmpz_mpoly_fit_length(polyq[w], k[w] + 1, ctx);
241 fmpz_set_si(polyq[w]->coeffs + k[w], tq);
242 polyq[w]->exps[k[w]] = texp;
243 }
244 c[0] = r1;
245 c[2] = c[1] = -(slong)(r1 < 0);
246 }
247 }
248 }
249 /* quotient non-small case */
250 if (!small)
251 {
252 fmpz_fdiv_qr(q, qc, qc, mb + w);
253 div_flag = fmpz_is_zero(qc);
254 d2 = !fmpz_is_zero(q);
255 if (d2)
256 {
257 k[w]++;
258 fmpz_mpoly_fit_length(polyq[w], k[w] + 1, ctx);
259 fmpz_set(polyq[w]->coeffs + k[w], q);
260 polyq[w]->exps[k[w]] = texp;
261 }
262 }
263 if (d2)
264 {
265 if (s[w] > 1)
266 {
267 i = 1;
268 x = chains[w] + i;
269 x->i = i;
270 x->j = k[w];
271 x->p = w;
272 x->next = NULL;
273 hinds[w][x->i] = 2*(x->j + 1) + 0;
274 _mpoly_heap_insert1(heap, exp3[w][i] +
275 polyq[w]->exps[k[w]], x,
276 &next_loc, &heap_len, maskhi);
277 }
278 s[w] = 1;
279 }
280 }
281 }
282 if (!div_flag)
283 {
284 l++;
285 _fmpz_mpoly_fit_length(&p2, &e2, allocr, l + 1, 1);
286
287 if (small)
288 {
289 fmpz_set_signed_uiuiui(p2 + l, c[2], c[1], c[0]);
290 fmpz_neg(p2 + l, p2 + l);
291 } else
292 {
293 fmpz_neg(p2 + l, qc);
294 }
295 e2[l] = exp;
296 }
297 }
298 }
299
300 cleanup:
301
302 for (i = 0; i < len; i++)
303 _fmpz_mpoly_set_length(polyq[i], k[i] + 1, ctx);
304 for (i = 0; i < len; i++)
305 fmpz_clear(mb + i);
306 fmpz_clear(qc);
307 fmpz_clear(q);
308
309 (*polyr) = p2;
310 (*expr) = e2;
311
312 TMP_END;
313 return l + 1;
314
315 exp_overflow:
316 for (i = 0; i < l; i++)
317 _fmpz_demote(p2 + i);
318 for (w = 0; w < len; w++)
319 {
320 for (i = 0; i < k[w]; i++)
321 _fmpz_demote(polyq[w]->coeffs + i);
322 k[w] = -WORD(1); /* we add 1 to this before exit */
323 }
324 l = -WORD(2); /* we add 1 to this upon return */
325 goto cleanup;
326 }
327
328 /*
329 As for divrem_monagan_pearce except that an array of divisor polynomials is
330 passed and an array of quotient polynomials is returned. These are not in
331 low level format.
332 */
_fmpz_mpoly_divrem_ideal_monagan_pearce(fmpz_mpoly_struct ** polyq,fmpz ** polyr,ulong ** expr,slong * allocr,const fmpz * poly2,const ulong * exp2,slong len2,fmpz_mpoly_struct * const * poly3,ulong * const * exp3,slong len,slong N,slong bits,const fmpz_mpoly_ctx_t ctx,const ulong * cmpmask)333 slong _fmpz_mpoly_divrem_ideal_monagan_pearce(fmpz_mpoly_struct ** polyq,
334 fmpz ** polyr, ulong ** expr, slong * allocr, const fmpz * poly2,
335 const ulong * exp2, slong len2, fmpz_mpoly_struct * const * poly3,
336 ulong * const * exp3, slong len, slong N, slong bits,
337 const fmpz_mpoly_ctx_t ctx, const ulong * cmpmask)
338 {
339 slong i, j, p, l, w;
340 slong next_loc;
341 slong * store, * store_base;
342 slong len3;
343 slong heap_len = 2; /* heap zero index unused */
344 mpoly_heap_s * heap;
345 mpoly_nheap_t ** chains;
346 slong ** hinds;
347 mpoly_nheap_t * x;
348 fmpz * p2 = *polyr;
349 ulong * e2 = *expr;
350 ulong * exp, * exps, * texp;
351 ulong ** exp_list;
352 ulong c[3]; /* for accumulating coefficients */
353 slong exp_next;
354 ulong mask = 0;
355 ulong * ub;
356 slong * k, * s;
357 fmpz_t qc, q;
358 fmpz * mb;
359 int small;
360 slong bits2, bits3;
361 int d1, d2, div_flag;
362 TMP_INIT;
363
364 if (N == 1)
365 return _fmpz_mpoly_divrem_ideal_monagan_pearce1(polyq, polyr, expr,
366 allocr, poly2, exp2, len2, poly3, exp3, len, bits, ctx, cmpmask[0]);
367
368 TMP_START;
369
370 fmpz_init(q);
371 fmpz_init(qc);
372
373 bits2 = _fmpz_vec_max_bits(poly2, len2);
374
375 chains = (mpoly_nheap_t **) TMP_ALLOC(len*sizeof(mpoly_nheap_t *));
376 hinds = (slong **) TMP_ALLOC(len*sizeof(mpoly_heap_t *));
377 bits3 = 0;
378 len3 = 0;
379 for (w = 0; w < len; w++)
380 {
381 chains[w] = (mpoly_nheap_t *) TMP_ALLOC((poly3[w]->length)*sizeof(mpoly_nheap_t));
382 hinds[w] = (slong *) TMP_ALLOC((poly3[w]->length)*sizeof(slong));
383 bits3 = FLINT_MAX(bits3, FLINT_ABS(fmpz_mpoly_max_bits(poly3[w])));
384 len3 += poly3[w]->length;
385 for (i = 0; i < poly3[w]->length; i++)
386 hinds[w][i] = 1;
387 }
388
389 /* allow one bit for sign, one bit for subtraction */
390 small = FLINT_ABS(bits2) <= (FLINT_ABS(bits3) +
391 FLINT_BIT_COUNT(len3) + FLINT_BITS - 2) &&
392 FLINT_ABS(bits3) <= FLINT_BITS - 2;
393
394 next_loc = len3 + 4; /* something bigger than heap can ever be */
395 heap = (mpoly_heap_s *) TMP_ALLOC((len3 + 1)*sizeof(mpoly_heap_s));
396 store = store_base = (slong *) TMP_ALLOC(3*len3*sizeof(mpoly_nheap_t *));
397
398 exps = (ulong *) TMP_ALLOC(len3*N*sizeof(ulong));
399 exp_list = (ulong **) TMP_ALLOC(len3*sizeof(ulong *));
400 texp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
401 exp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
402 k = (slong *) TMP_ALLOC(len*sizeof(slong));
403 s = (slong *) TMP_ALLOC(len*sizeof(slong));
404 ub = (ulong *) TMP_ALLOC(len*sizeof(ulong));
405 mb = (fmpz * ) TMP_ALLOC(len*sizeof(fmpz));
406
407 exp_next = 0;
408 for (i = 0; i < len3; i++)
409 exp_list[i] = exps + i*N;
410
411 mask = 0;
412 for (i = 0; i < FLINT_BITS/bits; i++)
413 mask = (mask << bits) + (UWORD(1) << (bits - 1));
414
415 for (w = 0; w < len; w++)
416 {
417 k[w] = -WORD(1);
418 s[w] = poly3[w]->length;
419 }
420 l = -WORD(1);
421
422 x = chains[0] + 0;
423 x->i = -WORD(1);
424 x->j = 0;
425 x->p = -WORD(1);
426 x->next = NULL;
427 heap[1].next = x;
428 heap[1].exp = exp_list[exp_next++];
429 mpoly_monomial_set(heap[1].exp, exp2, N);
430
431 for (i = 0; i < len; i++)
432 {
433 fmpz_init(mb + i);
434 fmpz_neg(mb + i, poly3[i]->coeffs);
435 ub[i] = ((ulong) FLINT_ABS(mb[i])) >> 1; /* abs(poly3[0])/2 */
436 }
437
438 while (heap_len > 1)
439 {
440 mpoly_monomial_set(exp, heap[1].exp, N);
441
442 c[0] = c[1] = c[2] = 0;
443 fmpz_zero(qc);
444
445 if (bits <= FLINT_BITS)
446 {
447 if (mpoly_monomial_overflows(exp, N, mask))
448 goto exp_overflow;
449 }
450 else
451 {
452 if (mpoly_monomial_overflows_mp(exp, N, bits))
453 goto exp_overflow;
454 }
455
456 while (heap_len > 1 && mpoly_monomial_equal(heap[1].exp, exp, N))
457 {
458 exp_list[--exp_next] = heap[1].exp;
459 x = _mpoly_heap_pop(heap, &heap_len, N, cmpmask);
460 do
461 {
462 *store++ = x->i;
463 *store++ = x->j;
464 *store++ = x->p;
465 if (x->i != -WORD(1))
466 hinds[x->p][x->i] |= WORD(1);
467
468 if (small)
469 {
470 if (x->i == -WORD(1))
471 _fmpz_mpoly_sub_uiuiui_fmpz(c, poly2 + x->j);
472 else
473 _fmpz_mpoly_addmul_uiuiui_fmpz(c,
474 poly3[x->p]->coeffs[x->i], polyq[x->p]->coeffs[x->j]);
475 } else
476 {
477 if (x->i == -WORD(1))
478 fmpz_sub(qc, qc, poly2 + x->j);
479 else
480 fmpz_addmul(qc, poly3[x->p]->coeffs + x->i,
481 polyq[x->p]->coeffs + x->j);
482 }
483 } while ((x = x->next) != NULL);
484 }
485
486 while (store > store_base)
487 {
488 p = *--store;
489 j = *--store;
490 i = *--store;
491 if (i == -WORD(1))
492 {
493 if (j + 1 < len2)
494 {
495 x = chains[0] + 0;
496 x->i = -WORD(1);
497 x->j = j + 1;
498 x->p = p;
499 x->next = NULL;
500 mpoly_monomial_set(exp_list[exp_next], exp2 + x->j*N, N);
501 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
502 &next_loc, &heap_len, N, cmpmask);
503 }
504 } else
505 {
506 /* should we go right */
507 if ( (i + 1 < poly3[p]->length)
508 && (hinds[p][i + 1] == 2*j + 1)
509 )
510 {
511 x = chains[p] + i + 1;
512 x->i = i + 1;
513 x->j = j;
514 x->p = p;
515 x->next = NULL;
516 hinds[p][x->i] = 2*(x->j + 1) + 0;
517 mpoly_monomial_add_mp(exp_list[exp_next], exp3[x->p] + x->i*N,
518 polyq[x->p]->exps + x->j*N, N);
519 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
520 &next_loc, &heap_len, N, cmpmask);
521 }
522
523 /* should we go up? */
524 if (j == k[p])
525 {
526 s[p]++;
527 } else if ( ((hinds[p][i] & 1) == 1)
528 && ((i == 1) || (hinds[p][i - 1] >= 2*(j + 2) + 1))
529 )
530 {
531 x = chains[p] + i;
532 x->i = i;
533 x->j = j + 1;
534 x->p = p;
535 x->next = NULL;
536 hinds[p][x->i] = 2*(x->j + 1) + 0;
537 mpoly_monomial_add_mp(exp_list[exp_next], exp3[x->p] + x->i*N,
538 polyq[x->p]->exps + x->j*N, N);
539 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
540 &next_loc, &heap_len, N, cmpmask);
541 }
542 }
543 }
544
545 if ((small && (c[2] != 0 || c[1] != 0 || c[0] != 0)) ||
546 (!small && !fmpz_is_zero(qc)))
547 {
548 div_flag = 0;
549 for (w = 0; w < len; w++)
550 {
551 if (bits <= FLINT_BITS)
552 d1 = mpoly_monomial_divides(texp, exp, exp3[w], N, mask);
553 else
554 d1 = mpoly_monomial_divides_mp(texp, exp, exp3[w], N, bits);
555
556 if (d1)
557 {
558 d2 = 0;
559 if (small)
560 {
561 ulong d[3];
562 if (0 > (slong) c[2])
563 mpn_neg(d, c, 3);
564 else
565 flint_mpn_copyi(d, c, 3);
566
567 if (d[2] != 0 || ub[w] <= d[1] ||
568 (ub[w] == 0 && 0 > (slong) d[0])) /* quotient not a small */
569 {
570 small = 0;
571 fmpz_set_signed_uiuiui(qc, c[2], c[1], c[0]);
572 } else /* quotient fits a small */
573 {
574 slong r1;
575 slong tq;
576 sdiv_qrnnd(tq, r1, c[1], c[0], mb[w]);
577 if (tq > COEFF_MAX || tq < COEFF_MIN)
578 {
579 small = 0;
580 fmpz_set_signed_uiuiui(qc, c[2], c[1], c[0]);
581 } else
582 {
583 d2 = tq != 0;
584 div_flag = r1 == 0;
585 if (d2)
586 {
587 k[w]++;
588 fmpz_mpoly_fit_length(polyq[w], k[w] + 1, ctx);
589 fmpz_set_si(polyq[w]->coeffs + k[w], tq);
590 mpoly_monomial_set(polyq[w]->exps + k[w]*N, texp, N);
591 }
592 c[0] = r1;
593 c[2] = c[1] = -(slong)(r1 < 0);
594 }
595 }
596 }
597 /* quotient non-small case */
598 if (!small)
599 {
600 fmpz_fdiv_qr(q, qc, qc, mb + w);
601 d2 = !fmpz_is_zero(q);
602 div_flag = fmpz_is_zero(qc);
603 if (d2)
604 {
605 k[w]++;
606 fmpz_mpoly_fit_length(polyq[w], k[w] + 1, ctx);
607 fmpz_set(polyq[w]->coeffs + k[w], q);
608 mpoly_monomial_set(polyq[w]->exps + k[w]*N, texp, N);
609 }
610 }
611 if (d2)
612 {
613 if (s[w] > 1)
614 {
615 i = 1;
616 x = chains[w] + i;
617 x->i = i;
618 x->j = k[w];
619 x->p = w;
620 x->next = NULL;
621 hinds[w][x->i] = 2*(x->j + 1) + 0;
622 mpoly_monomial_add_mp(exp_list[exp_next], exp3[w] + i*N,
623 polyq[w]->exps + k[w]*N, N);
624 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
625 &next_loc, &heap_len, N, cmpmask);
626 }
627 s[w] = 1;
628 }
629 }
630 }
631 if (!div_flag)
632 {
633 l++;
634 _fmpz_mpoly_fit_length(&p2, &e2, allocr, l + 1, N);
635 if (small)
636 {
637 fmpz_set_signed_uiuiui(p2 + l, c[2], c[1], c[0]);
638 fmpz_neg(p2 + l, p2 + l);
639 } else
640 {
641 fmpz_neg(p2 + l, qc);
642 }
643 mpoly_monomial_set(e2 + l*N, exp, N);
644 }
645 }
646 }
647
648 cleanup2:
649
650 for (i = 0; i < len; i++)
651 _fmpz_mpoly_set_length(polyq[i], k[i] + 1, ctx);
652
653 for (i = 0; i < len; i++)
654 fmpz_clear(mb + i);
655 fmpz_clear(qc);
656 fmpz_clear(q);
657
658 (*polyr) = p2;
659 (*expr) = e2;
660
661 TMP_END;
662
663 return l + 1;
664
665 exp_overflow:
666 for (i = 0; i < l; i++)
667 _fmpz_demote(p2 + i);
668 for (w = 0; w < len; w++)
669 {
670 for (i = 0; i < k[w]; i++)
671 _fmpz_demote(polyq[w]->coeffs + i);
672 k[w] = -WORD(1);
673 }
674 l = -WORD(2);
675 goto cleanup2;
676
677 }
678
679 /* Assumes divisor polys don't alias any output polys */
fmpz_mpoly_divrem_ideal_monagan_pearce(fmpz_mpoly_struct ** q,fmpz_mpoly_t r,const fmpz_mpoly_t poly2,fmpz_mpoly_struct * const * poly3,slong len,const fmpz_mpoly_ctx_t ctx)680 void fmpz_mpoly_divrem_ideal_monagan_pearce(fmpz_mpoly_struct ** q, fmpz_mpoly_t r,
681 const fmpz_mpoly_t poly2, fmpz_mpoly_struct * const * poly3, slong len,
682 const fmpz_mpoly_ctx_t ctx)
683 {
684 slong i, exp_bits, N, lenr = 0;
685 slong len3 = 0;
686 ulong * cmpmask;
687 ulong * exp2;
688 ulong ** exp3;
689 int free2 = 0;
690 int * free3;
691 fmpz_mpoly_t temp2;
692 fmpz_mpoly_struct * tr;
693 TMP_INIT;
694
695 /* check none of the divisor polynomials is zero */
696 for (i = 0; i < len; i++)
697 {
698 if (poly3[i]->length == 0)
699 flint_throw(FLINT_DIVZERO,
700 "Divide by zero in fmpz_mpoly_divrem_ideal_monagan_pearce");
701
702 len3 = FLINT_MAX(len3, poly3[i]->length);
703 }
704
705 /* dividend is zero, write out quotients and remainder */
706 if (poly2->length == 0)
707 {
708 for (i = 0; i < len; i++)
709 {
710 fmpz_mpoly_zero(q[i], ctx);
711 }
712
713 fmpz_mpoly_zero(r, ctx);
714
715 return;
716 }
717
718
719 TMP_START;
720
721 free3 = (int *) TMP_ALLOC(len*sizeof(int));
722
723 exp3 = (ulong **) TMP_ALLOC(len*sizeof(ulong *));
724
725 /* compute maximum degrees that can occur in any input or output polys */
726 exp_bits = poly2->bits;
727 for (i = 0; i < len; i++)
728 exp_bits = FLINT_MAX(exp_bits, poly3[i]->bits);
729
730 exp_bits = mpoly_fix_bits(exp_bits, ctx->minfo);
731
732 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
733 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
734 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
735
736 /* ensure input exponents packed to same size as output exponents */
737 exp2 = poly2->exps;
738 free2 = 0;
739 if (exp_bits > poly2->bits)
740 {
741 free2 = 1;
742 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
743 mpoly_repack_monomials(exp2, exp_bits, poly2->exps, poly2->bits,
744 poly2->length, ctx->minfo);
745 }
746
747 for (i = 0; i < len; i++)
748 {
749 exp3[i] = poly3[i]->exps;
750 free3[i] = 0;
751 if (exp_bits > poly3[i]->bits)
752 {
753 free3[i] = 1;
754 exp3[i] = (ulong *) flint_malloc(N*poly3[i]->length*sizeof(ulong));
755 mpoly_repack_monomials(exp3[i], exp_bits, poly3[i]->exps,
756 poly3[i]->bits, poly3[i]->length, ctx->minfo);
757 }
758 fmpz_mpoly_fit_length(q[i], 1, ctx);
759 fmpz_mpoly_fit_bits(q[i], exp_bits, ctx);
760 q[i]->bits = exp_bits;
761 }
762
763 /* check leading mon. of at least one divisor is at most that of dividend */
764 for (i = 0; i < len; i++)
765 {
766 if (!mpoly_monomial_lt(exp2, exp3[i], N, cmpmask))
767 break;
768 }
769
770 if (i == len)
771 {
772 fmpz_mpoly_set(r, poly2, ctx);
773 for (i = 0; i < len; i++)
774 fmpz_mpoly_zero(q[i], ctx);
775
776 goto cleanup3;
777 }
778
779 /* take care of aliasing */
780 if (r == poly2)
781 {
782 fmpz_mpoly_init2(temp2, len3, ctx);
783 fmpz_mpoly_fit_bits(temp2, exp_bits, ctx);
784 temp2->bits = exp_bits;
785
786 tr = temp2;
787 } else
788 {
789 fmpz_mpoly_fit_length(r, len3, ctx);
790 fmpz_mpoly_fit_bits(r, exp_bits, ctx);
791 r->bits = exp_bits;
792
793 tr = r;
794 }
795
796 /* do division with remainder */
797 while (1)
798 {
799 slong old_exp_bits = exp_bits;
800 ulong * old_exp2 = exp2, * old_exp3;
801
802 lenr = _fmpz_mpoly_divrem_ideal_monagan_pearce(q, &tr->coeffs, &tr->exps,
803 &tr->alloc, poly2->coeffs, exp2, poly2->length,
804 poly3, exp3, len, N, exp_bits, ctx, cmpmask);
805
806 if (lenr >= 0) /* check if division was successful */
807 break;
808
809 exp_bits = mpoly_fix_bits(exp_bits + 1, ctx->minfo);
810 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
811 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
812 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
813
814 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
815 mpoly_repack_monomials(exp2, exp_bits, old_exp2, old_exp_bits,
816 poly2->length, ctx->minfo);
817
818 if (free2)
819 flint_free(old_exp2);
820
821 free2 = 1;
822
823 fmpz_mpoly_fit_bits(tr, exp_bits, ctx);
824 tr->bits = exp_bits;
825
826 for (i = 0; i < len; i++)
827 {
828 old_exp3 = exp3[i];
829
830 exp3[i] = (ulong *) flint_malloc(N*poly3[i]->length*sizeof(ulong));
831 mpoly_repack_monomials(exp3[i], exp_bits, old_exp3, old_exp_bits,
832 poly3[i]->length, ctx->minfo);
833
834 if (free3[i])
835 flint_free(old_exp3);
836
837 free3[i] = 1;
838
839 fmpz_mpoly_fit_bits(q[i], exp_bits, ctx);
840 q[i]->bits = exp_bits;
841 }
842
843 }
844
845 /* take care of aliasing */
846 if (r == poly2)
847 {
848 fmpz_mpoly_swap(temp2, r, ctx);
849 fmpz_mpoly_clear(temp2, ctx);
850 }
851
852 _fmpz_mpoly_set_length(r, lenr, ctx);
853
854 cleanup3:
855
856 if (free2)
857 flint_free(exp2);
858
859 for (i = 0; i < len; i++)
860 {
861 if (free3[i])
862 flint_free(exp3[i]);
863 }
864
865 TMP_END;
866 }
867