1 /*
2 Copyright (C) 2019 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 <stdlib.h>
13 #include "fq_nmod_mpoly.h"
14
15 /*
16 As for divrem_monagan_pearce except that an array of divisor polynomials is
17 passed and an array of quotient polynomials is returned. These are not in
18 low level format.
19 */
_fq_nmod_mpoly_divrem_ideal_monagan_pearce(fq_nmod_mpoly_struct ** polyq,fq_nmod_struct ** polyr,ulong ** expr,slong * allocr,const fq_nmod_struct * poly2,const ulong * exp2,slong len2,fq_nmod_mpoly_struct * const * poly3,ulong * const * exp3,slong len,slong N,slong bits,const fq_nmod_mpoly_ctx_t ctx,const ulong * cmpmask)20 slong _fq_nmod_mpoly_divrem_ideal_monagan_pearce(fq_nmod_mpoly_struct ** polyq,
21 fq_nmod_struct ** polyr, ulong ** expr, slong * allocr,
22 const fq_nmod_struct * poly2, const ulong * exp2, slong len2,
23 fq_nmod_mpoly_struct * const * poly3, ulong * const * exp3,
24 slong len, slong N, slong bits,
25 const fq_nmod_mpoly_ctx_t ctx, const ulong * cmpmask)
26 {
27 slong i, j, p, r_len, w;
28 slong next_loc;
29 slong * store, * store_base;
30 slong len3;
31 slong heap_len = 2; /* heap zero index unused */
32 mpoly_heap_s * heap;
33 mpoly_nheap_t ** chains;
34 slong ** hinds;
35 mpoly_nheap_t * x;
36 fq_nmod_struct * r_coeff = *polyr;
37 ulong * r_exp = *expr;
38 ulong * exp, * exps, * texp;
39 ulong ** exp_list;
40 slong exp_next;
41 ulong mask;
42 slong * q_len, * s;
43 fq_nmod_struct * lc_minus_inv;
44 fq_nmod_t acc, pp;
45 TMP_INIT;
46
47 TMP_START;
48
49 chains = (mpoly_nheap_t **) TMP_ALLOC(len*sizeof(mpoly_nheap_t *));
50 hinds = (slong **) TMP_ALLOC(len*sizeof(mpoly_heap_t *));
51 len3 = 0;
52 for (w = 0; w < len; w++)
53 {
54 chains[w] = (mpoly_nheap_t *) TMP_ALLOC((poly3[w]->length)*sizeof(mpoly_nheap_t));
55 hinds[w] = (slong *) TMP_ALLOC((poly3[w]->length)*sizeof(slong));
56 len3 += poly3[w]->length;
57 for (i = 0; i < poly3[w]->length; i++)
58 hinds[w][i] = 1;
59 }
60
61 next_loc = len3 + 4; /* something bigger than heap can ever be */
62 heap = (mpoly_heap_s *) TMP_ALLOC((len3 + 1)*sizeof(mpoly_heap_s));
63 store = store_base = (slong *) TMP_ALLOC(3*len3*sizeof(mpoly_nheap_t *));
64
65 exps = (ulong *) TMP_ALLOC(len3*N*sizeof(ulong));
66 exp_list = (ulong **) TMP_ALLOC(len3*sizeof(ulong *));
67 texp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
68 exp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
69 q_len = (slong *) TMP_ALLOC(len*sizeof(slong));
70 s = (slong *) TMP_ALLOC(len*sizeof(slong));
71
72 exp_next = 0;
73 for (i = 0; i < len3; i++)
74 exp_list[i] = exps + i*N;
75
76 mask = 0;
77 for (i = 0; i < FLINT_BITS/bits; i++)
78 mask = (mask << bits) + (UWORD(1) << (bits - 1));
79
80 for (w = 0; w < len; w++)
81 {
82 q_len[w] = WORD(0);
83 s[w] = poly3[w]->length;
84 }
85 r_len = WORD(0);
86
87 x = chains[0] + 0;
88 x->i = -WORD(1);
89 x->j = 0;
90 x->p = -WORD(1);
91 x->next = NULL;
92 heap[1].next = x;
93 heap[1].exp = exp_list[exp_next++];
94 mpoly_monomial_set(heap[1].exp, exp2, N);
95
96 /* precompute leading coeff info */
97 fq_nmod_init(pp, ctx->fqctx);
98 fq_nmod_init(acc, ctx->fqctx);
99 lc_minus_inv = (fq_nmod_struct *) TMP_ALLOC(len*sizeof(fq_nmod_struct));
100 for (w = 0; w < len; w++)
101 {
102 fq_nmod_init(lc_minus_inv + w, ctx->fqctx);
103 fq_nmod_inv(lc_minus_inv + w, poly3[w]->coeffs + 0, ctx->fqctx);
104 fq_nmod_neg(lc_minus_inv + w, lc_minus_inv + w, ctx->fqctx);
105 }
106
107 while (heap_len > 1)
108 {
109 mpoly_monomial_set(exp, heap[1].exp, N);
110
111 if (bits <= FLINT_BITS)
112 {
113 if (mpoly_monomial_overflows(exp, N, mask))
114 goto exp_overflow;
115 }
116 else
117 {
118 if (mpoly_monomial_overflows_mp(exp, N, bits))
119 goto exp_overflow;
120 }
121
122 fq_nmod_zero(acc, ctx->fqctx);
123 do
124 {
125 exp_list[--exp_next] = heap[1].exp;
126 x = _mpoly_heap_pop(heap, &heap_len, N, cmpmask);
127 do
128 {
129 *store++ = x->i;
130 *store++ = x->j;
131 *store++ = x->p;
132 if (x->i != -WORD(1))
133 {
134 hinds[x->p][x->i] |= WORD(1);
135 }
136
137 if (x->i == -WORD(1))
138 {
139 fq_nmod_sub(acc, acc, poly2 + x->j, ctx->fqctx);
140 }
141 else
142 {
143 fq_nmod_mul(pp, poly3[x->p]->coeffs + x->i,
144 polyq[x->p]->coeffs + x->j, ctx->fqctx);
145 fq_nmod_add(acc, acc, pp, ctx->fqctx);
146 }
147 } while ((x = x->next) != NULL);
148 } while (heap_len > 1 && mpoly_monomial_equal(heap[1].exp, exp, N));
149
150 while (store > store_base)
151 {
152 p = *--store;
153 j = *--store;
154 i = *--store;
155 if (i == -WORD(1))
156 {
157 if (j + 1 < len2)
158 {
159 x = chains[0] + 0;
160 x->i = -WORD(1);
161 x->j = j + 1;
162 x->p = p;
163 x->next = NULL;
164 mpoly_monomial_set(exp_list[exp_next], exp2 + x->j*N, N);
165 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
166 &next_loc, &heap_len, N, cmpmask);
167 }
168 }
169 else
170 {
171 /* should we go right? */
172 if ( (i + 1 < poly3[p]->length)
173 && (hinds[p][i + 1] == 2*j + 1)
174 )
175 {
176 x = chains[p] + i + 1;
177 x->i = i + 1;
178 x->j = j;
179 x->p = p;
180 x->next = NULL;
181 hinds[p][x->i] = 2*(x->j + 1) + 0;
182 mpoly_monomial_add_mp(exp_list[exp_next], exp3[x->p] + x->i*N,
183 polyq[x->p]->exps + x->j*N, N);
184 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
185 &next_loc, &heap_len, N, cmpmask);
186 }
187
188 /* should we go up? */
189 if (j + 1 == q_len[p])
190 {
191 s[p]++;
192 }
193 else if ( ((hinds[p][i] & 1) == 1)
194 && ((i == 1) || (hinds[p][i - 1] >= 2*(j + 2) + 1)) )
195 {
196 x = chains[p] + i;
197 x->i = i;
198 x->j = j + 1;
199 x->p = p;
200 x->next = NULL;
201 hinds[p][x->i] = 2*(x->j + 1) + 0;
202 mpoly_monomial_add_mp(exp_list[exp_next], exp3[x->p] + x->i*N,
203 polyq[x->p]->exps + x->j*N, N);
204 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
205 &next_loc, &heap_len, N, cmpmask);
206 }
207 }
208 }
209
210 if (fq_nmod_is_zero(acc, ctx->fqctx))
211 continue;
212
213 for (w = 0; w < len; w++)
214 {
215 int divides;
216
217 if (bits <= FLINT_BITS)
218 divides = mpoly_monomial_divides(texp, exp, exp3[w] + N*0, N, mask);
219 else
220 divides = mpoly_monomial_divides_mp(texp, exp, exp3[w] + N*0, N, bits);
221
222 if (divides)
223 {
224 fq_nmod_mpoly_fit_length(polyq[w], q_len[w] + 1, ctx);
225 fq_nmod_mul(polyq[w]->coeffs + q_len[w],
226 acc, lc_minus_inv + w, ctx->fqctx);
227 mpoly_monomial_set(polyq[w]->exps + N*q_len[w], texp, N);
228 if (s[w] > 1)
229 {
230 i = 1;
231 x = chains[w] + i;
232 x->i = i;
233 x->j = q_len[w];
234 x->p = w;
235 x->next = NULL;
236 hinds[w][x->i] = 2*(x->j + 1) + 0;
237 mpoly_monomial_add_mp(exp_list[exp_next], exp3[w] + N*x->i,
238 polyq[w]->exps + N*x->j, N);
239 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
240 &next_loc, &heap_len, N, cmpmask);
241 }
242 s[w] = 1;
243 q_len[w]++;
244 /* break out of w for loop and continue in heap loop */
245 goto break_continue;
246 }
247 }
248
249 /* if get here, no leading terms divided */
250 _fq_nmod_mpoly_fit_length(&r_coeff, &r_exp, allocr, r_len + 1, N, ctx->fqctx);
251 fq_nmod_neg(r_coeff + r_len, acc, ctx->fqctx);
252 mpoly_monomial_set(r_exp + r_len*N, exp, N);
253 r_len++;
254
255 break_continue:
256
257 (void)(NULL);
258 }
259
260 cleanup2:
261
262 fq_nmod_clear(pp, ctx->fqctx);
263 fq_nmod_clear(acc, ctx->fqctx);
264 for (w = 0; w < len; w++)
265 {
266 fq_nmod_clear(lc_minus_inv + w, ctx->fqctx);
267 }
268
269 for (i = 0; i < len; i++)
270 {
271 _fq_nmod_mpoly_set_length(polyq[i], q_len[i], ctx);
272 }
273
274 (*polyr) = r_coeff;
275 (*expr) = r_exp;
276
277 TMP_END;
278
279 return r_len;
280
281 exp_overflow:
282 for (w = 0; w < len; w++)
283 q_len[w] = WORD(0);
284
285 r_len = -WORD(1);
286 goto cleanup2;
287
288 }
289
290 /* Assumes divisor polys don't alias any output polys */
fq_nmod_mpoly_divrem_ideal_monagan_pearce(fq_nmod_mpoly_struct ** q,fq_nmod_mpoly_t r,const fq_nmod_mpoly_t poly2,fq_nmod_mpoly_struct * const * poly3,slong len,const fq_nmod_mpoly_ctx_t ctx)291 void fq_nmod_mpoly_divrem_ideal_monagan_pearce(fq_nmod_mpoly_struct ** q,
292 fq_nmod_mpoly_t r, const fq_nmod_mpoly_t poly2, fq_nmod_mpoly_struct * const * poly3,
293 slong len, const fq_nmod_mpoly_ctx_t ctx)
294 {
295 slong i, exp_bits, N, lenr = 0;
296 slong len3 = 0;
297 ulong * cmpmask;
298 ulong * exp2;
299 ulong ** exp3;
300 int free2 = 0;
301 int * free3;
302 fq_nmod_mpoly_t temp2;
303 fq_nmod_mpoly_struct * tr;
304 TMP_INIT;
305
306 for (i = 0; i < len; i++)
307 {
308 len3 = FLINT_MAX(len3, poly3[i]->length);
309 if (poly3[i]->length == 0)
310 {
311 flint_throw(FLINT_DIVZERO,
312 "Divide by zero in fq_nmod_mpoly_divrem_ideal_monagan_pearce");
313 }
314 }
315
316 /* dividend is zero, write out quotients and remainder */
317 if (poly2->length == 0)
318 {
319 for (i = 0; i < len; i++)
320 fq_nmod_mpoly_zero(q[i], ctx);
321 fq_nmod_mpoly_zero(r, ctx);
322 return;
323 }
324
325 TMP_START;
326
327 free3 = (int *) TMP_ALLOC(len*sizeof(int));
328 exp3 = (ulong **) TMP_ALLOC(len*sizeof(ulong *));
329
330 /* compute maximum degrees that can occur in any input or output polys */
331 exp_bits = poly2->bits;
332 for (i = 0; i < len; i++)
333 {
334 exp_bits = FLINT_MAX(exp_bits, poly3[i]->bits);
335 }
336 exp_bits = mpoly_fix_bits(exp_bits, ctx->minfo);
337
338 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
339 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
340 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
341
342 /* ensure input exponents packed to same size as output exponents */
343 exp2 = poly2->exps;
344 free2 = 0;
345 if (exp_bits > poly2->bits)
346 {
347 free2 = 1;
348 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
349 mpoly_repack_monomials(exp2, exp_bits, poly2->exps, poly2->bits,
350 poly2->length, ctx->minfo);
351 }
352
353 for (i = 0; i < len; i++)
354 {
355 exp3[i] = poly3[i]->exps;
356 free3[i] = 0;
357 if (exp_bits > poly3[i]->bits)
358 {
359 free3[i] = 1;
360 exp3[i] = (ulong *) flint_malloc(N*poly3[i]->length*sizeof(ulong));
361 mpoly_repack_monomials(exp3[i], exp_bits, poly3[i]->exps,
362 poly3[i]->bits, poly3[i]->length, ctx->minfo);
363 }
364 fq_nmod_mpoly_fit_length(q[i], 1, ctx);
365 fq_nmod_mpoly_fit_bits(q[i], exp_bits, ctx);
366 q[i]->bits = exp_bits;
367 }
368
369 /* check leading mon. of at least one divisor is at most that of dividend */
370 for (i = 0; i < len; i++)
371 {
372 if (!mpoly_monomial_lt(exp2, exp3[i], N, cmpmask))
373 break;
374 }
375
376 if (i == len)
377 {
378 fq_nmod_mpoly_set(r, poly2, ctx);
379 for (i = 0; i < len; i++)
380 fq_nmod_mpoly_zero(q[i], ctx);
381
382 goto cleanup3;
383 }
384
385 /* take care of aliasing */
386 if (r == poly2)
387 {
388 fq_nmod_mpoly_init2(temp2, len3, ctx);
389 fq_nmod_mpoly_fit_bits(temp2, exp_bits, ctx);
390 temp2->bits = exp_bits;
391 tr = temp2;
392 }
393 else
394 {
395 fq_nmod_mpoly_fit_length(r, len3, ctx);
396 fq_nmod_mpoly_fit_bits(r, exp_bits, ctx);
397 r->bits = exp_bits;
398 tr = r;
399 }
400
401 /* do division with remainder */
402 while (1)
403 {
404 slong old_exp_bits = exp_bits;
405 ulong * old_exp2 = exp2, * old_exp3;
406
407 lenr = _fq_nmod_mpoly_divrem_ideal_monagan_pearce(q, &tr->coeffs, &tr->exps,
408 &tr->alloc, poly2->coeffs, exp2, poly2->length,
409 poly3, exp3, len, N, exp_bits, ctx, cmpmask);
410
411 if (lenr >= 0) /* check if division was successful */
412 break;
413
414 exp_bits = mpoly_fix_bits(exp_bits + 1, ctx->minfo);
415 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
416 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
417 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
418
419 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
420 mpoly_repack_monomials(exp2, exp_bits, old_exp2, old_exp_bits,
421 poly2->length, ctx->minfo);
422
423 if (free2)
424 flint_free(old_exp2);
425
426 free2 = 1;
427
428 fq_nmod_mpoly_fit_bits(tr, exp_bits, ctx);
429 tr->bits = exp_bits;
430
431 for (i = 0; i < len; i++)
432 {
433 old_exp3 = exp3[i];
434
435 exp3[i] = (ulong *) flint_malloc(N*poly3[i]->length*sizeof(ulong));
436 mpoly_repack_monomials(exp3[i], exp_bits, old_exp3, old_exp_bits,
437 poly3[i]->length, ctx->minfo);
438
439 if (free3[i])
440 flint_free(old_exp3);
441
442 free3[i] = 1;
443
444 fq_nmod_mpoly_fit_bits(q[i], exp_bits, ctx);
445 q[i]->bits = exp_bits;
446 }
447 }
448
449 /* take care of aliasing */
450 if (r == poly2)
451 {
452 fq_nmod_mpoly_swap(temp2, r, ctx);
453 fq_nmod_mpoly_clear(temp2, ctx);
454 }
455
456 _fq_nmod_mpoly_set_length(r, lenr, ctx);
457
458 cleanup3:
459
460 if (free2)
461 flint_free(exp2);
462
463 for (i = 0; i < len; i++)
464 {
465 if (free3[i])
466 flint_free(exp3[i]);
467 }
468
469 TMP_END;
470 }
471