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