1 /*
2 Copyright (C) 2020 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fq_zech_mpoly.h"
13
_fq_zech_mpoly_divrem_monagan_pearce(slong * lenr,fq_zech_struct ** polyq,ulong ** expq,slong * allocq,fq_zech_struct ** polyr,ulong ** expr,slong * allocr,const fq_zech_struct * coeff2,const ulong * exp2,slong len2,const fq_zech_struct * coeff3,const ulong * exp3,slong len3,slong bits,slong N,const ulong * cmpmask,const fq_zech_ctx_t fqctx)14 static slong _fq_zech_mpoly_divrem_monagan_pearce(slong * lenr,
15 fq_zech_struct ** polyq, ulong ** expq, slong * allocq,
16 fq_zech_struct ** polyr, ulong ** expr, slong * allocr,
17 const fq_zech_struct * coeff2, const ulong * exp2, slong len2,
18 const fq_zech_struct * coeff3, const ulong * exp3, slong len3,
19 slong bits, slong N, const ulong * cmpmask, const fq_zech_ctx_t fqctx)
20 {
21 slong i, j, q_len, r_len, s;
22 slong next_loc;
23 slong heap_len = 2; /* heap zero index unused */
24 mpoly_heap_s * heap;
25 mpoly_heap_t * chain;
26 slong * store, * store_base;
27 mpoly_heap_t * x;
28 fq_zech_struct * q_coeff = *polyq;
29 fq_zech_struct * r_coeff = *polyr;
30 ulong * q_exp = *expq;
31 ulong * r_exp = *expr;
32 ulong * exp, * exps;
33 ulong ** exp_list;
34 slong exp_next;
35 ulong mask;
36 slong * hind;
37 int lt_divides;
38 fq_zech_t lc_minus_inv, pp;
39 TMP_INIT;
40
41 TMP_START;
42 fq_zech_init(pp, fqctx);
43 fq_zech_init(lc_minus_inv, fqctx);
44
45 /* alloc array of heap nodes which can be chained together */
46 next_loc = len3 + 4; /* something bigger than heap can ever be */
47 heap = (mpoly_heap_s *) TMP_ALLOC((len3 + 1)*sizeof(mpoly_heap_s));
48 chain = (mpoly_heap_t *) TMP_ALLOC(len3*sizeof(mpoly_heap_t));
49 store = store_base = (slong *) TMP_ALLOC(2*len3*sizeof(mpoly_heap_t *));
50
51 /* array of exponent vectors, each of "N" words */
52 exps = (ulong *) TMP_ALLOC(len3*N*sizeof(ulong));
53 /* list of pointers to available exponent vectors */
54 exp_list = (ulong **) TMP_ALLOC(len3*sizeof(ulong *));
55 /* space to save copy of current exponent vector */
56 exp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
57 /* set up list of available exponent vectors */
58 exp_next = 0;
59 for (i = 0; i < len3; i++)
60 exp_list[i] = exps + i*N;
61
62 /* space for flagged heap indicies */
63 hind = (slong *) TMP_ALLOC(len3*sizeof(slong));
64 for (i = 0; i < len3; i++)
65 hind[i] = 1;
66
67 mask = bits <= FLINT_BITS ? mpoly_overflow_mask_sp(bits) : 0;
68
69 q_len = WORD(0);
70 r_len = WORD(0);
71
72 /* s is the number of terms * (latest quotient) we should put into heap */
73 s = len3;
74
75 /* insert (-1, 0, exp2[0]) into heap */
76 x = chain + 0;
77 x->i = -WORD(1);
78 x->j = 0;
79 x->next = NULL;
80 heap[1].next = x;
81 heap[1].exp = exp_list[exp_next++];
82 mpoly_monomial_set(heap[1].exp, exp2, N);
83
84 /* precompute leading cofficient info */
85 fq_zech_inv(lc_minus_inv, coeff3 + 0, fqctx);
86 fq_zech_neg(lc_minus_inv, lc_minus_inv, fqctx);
87
88 while (heap_len > 1)
89 {
90 _fq_zech_mpoly_fit_length(&q_coeff, &q_exp, allocq, q_len + 1, N, fqctx);
91
92 mpoly_monomial_set(exp, heap[1].exp, N);
93
94 if (bits <= FLINT_BITS)
95 {
96 if (mpoly_monomial_overflows(exp, N, mask))
97 goto exp_overflow2;
98 lt_divides = mpoly_monomial_divides(q_exp + q_len*N, exp, exp3, N, mask);
99 }
100 else
101 {
102 if (mpoly_monomial_overflows_mp(exp, N, bits))
103 goto exp_overflow2;
104 lt_divides = mpoly_monomial_divides_mp(q_exp + q_len*N, exp, exp3, N, bits);
105 }
106
107 fq_zech_zero(q_coeff + q_len, fqctx);
108 do {
109 exp_list[--exp_next] = heap[1].exp;
110 x = _mpoly_heap_pop(heap, &heap_len, N, cmpmask);
111 do {
112 *store++ = x->i;
113 *store++ = x->j;
114 if (x->i != -WORD(1))
115 hind[x->i] |= WORD(1);
116
117 if (x->i == -WORD(1))
118 {
119 fq_zech_sub(q_coeff + q_len, q_coeff + q_len, coeff2 + x->j, fqctx);
120 }
121 else
122 {
123 fq_zech_mul(pp, coeff3 + x->i, q_coeff + x->j, fqctx);
124 fq_zech_add(q_coeff + q_len, q_coeff + q_len, pp, fqctx);
125 }
126 } while ((x = x->next) != NULL);
127 } while (heap_len > 1 && mpoly_monomial_equal(heap[1].exp, exp, N));
128
129 /* process nodes taken from the heap */
130 while (store > store_base)
131 {
132 j = *--store;
133 i = *--store;
134
135 if (i == -WORD(1))
136 {
137 /* take next dividend term */
138 if (j + 1 < len2)
139 {
140 x = chain + 0;
141 x->i = i;
142 x->j = j + 1;
143 x->next = NULL;
144 mpoly_monomial_set(exp_list[exp_next], exp2 + x->j*N, N);
145 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
146 &next_loc, &heap_len, N, cmpmask);
147 }
148 }
149 else
150 {
151 /* should we go right? */
152 if ( (i + 1 < len3)
153 && (hind[i + 1] == 2*j + 1)
154 )
155 {
156 x = chain + i + 1;
157 x->i = i + 1;
158 x->j = j;
159 x->next = NULL;
160 hind[x->i] = 2*(x->j + 1) + 0;
161 mpoly_monomial_add_mp(exp_list[exp_next], exp3 + x->i*N,
162 q_exp + x->j*N, N);
163 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
164 &next_loc, &heap_len, N, cmpmask);
165 }
166 /* should we go up? */
167 if (j + 1 == q_len)
168 {
169 s++;
170 }
171 else if ( ((hind[i] & 1) == 1)
172 && ((i == 1) || (hind[i - 1] >= 2*(j + 2) + 1)) )
173 {
174 x = chain + i;
175 x->i = i;
176 x->j = j + 1;
177 x->next = NULL;
178 hind[x->i] = 2*(x->j + 1) + 0;
179 mpoly_monomial_add_mp(exp_list[exp_next], exp3 + x->i*N,
180 q_exp + x->j*N, N);
181 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
182 &next_loc, &heap_len, N, cmpmask);
183 }
184 }
185 }
186
187 /* try to divide accumulated term by leading term */
188 if (fq_zech_is_zero(q_coeff + q_len, fqctx))
189 {
190 continue;
191 }
192 if (!lt_divides)
193 {
194 _fq_zech_mpoly_fit_length(&r_coeff, &r_exp, allocr, r_len + 1, N, fqctx);
195 fq_zech_neg(r_coeff + r_len, q_coeff + q_len, fqctx);
196 mpoly_monomial_set(r_exp + r_len*N, exp, N);
197 r_len++;
198 continue;
199 }
200
201 fq_zech_mul(q_coeff + q_len, q_coeff + q_len, lc_minus_inv, fqctx);
202
203 /* put newly generated quotient term back into the heap if neccesary */
204 if (s > 1)
205 {
206 i = 1;
207 x = chain + i;
208 x->i = i;
209 x->j = q_len;
210 x->next = NULL;
211 hind[x->i] = 2*(x->j + 1) + 0;
212 mpoly_monomial_add_mp(exp_list[exp_next], exp3 + x->i*N,
213 q_exp + x->j*N, N);
214 exp_next += _mpoly_heap_insert(heap, exp_list[exp_next], x,
215 &next_loc, &heap_len, N, cmpmask);
216 }
217 s = 1;
218 q_len++;
219 }
220
221
222 cleanup2:
223
224 (*polyq) = q_coeff;
225 (*expq) = q_exp;
226 (*polyr) = r_coeff;
227 (*expr) = r_exp;
228
229 /* set remainder poly length */
230 (*lenr) = r_len;
231
232 TMP_END;
233 fq_zech_clear(pp, fqctx);
234 fq_zech_clear(lc_minus_inv, fqctx);
235
236 /* return quotient poly length */
237 return q_len;
238
239 exp_overflow2:
240 q_len = 0;
241 r_len = 0;
242 goto cleanup2;
243 }
244
fq_zech_mpoly_divrem_monagan_pearce(fq_zech_mpoly_t q,fq_zech_mpoly_t r,const fq_zech_mpoly_t poly2,const fq_zech_mpoly_t poly3,const fq_zech_mpoly_ctx_t ctx)245 void fq_zech_mpoly_divrem_monagan_pearce(fq_zech_mpoly_t q, fq_zech_mpoly_t r,
246 const fq_zech_mpoly_t poly2, const fq_zech_mpoly_t poly3,
247 const fq_zech_mpoly_ctx_t ctx)
248 {
249 slong exp_bits, N, lenq = 0, lenr = 0;
250 ulong * exp2 = poly2->exps, * exp3 = poly3->exps;
251 ulong * cmpmask;
252 int free2 = 0, free3 = 0;
253 fq_zech_mpoly_t temp1, temp2;
254 fq_zech_mpoly_struct * tq, * tr;
255
256 if (poly3->length == 0)
257 {
258 flint_throw(FLINT_DIVZERO, "Divide by zero in nmod_mpoly_divrem_monagan_pearce");
259 }
260
261 if (poly2->length == 0)
262 {
263 fq_zech_mpoly_zero(q, ctx);
264 fq_zech_mpoly_zero(r, ctx);
265 return;
266 }
267
268 exp_bits = FLINT_MAX(poly2->bits, poly3->bits);
269 exp_bits = mpoly_fix_bits(exp_bits, ctx->minfo);
270
271 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
272 cmpmask = (ulong *) flint_malloc(N*sizeof(ulong));
273 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
274
275 /* ensure input exponents packed to same size as output exponents */
276 if (exp_bits > poly2->bits)
277 {
278 free2 = 1;
279 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
280 mpoly_repack_monomials(exp2, exp_bits, poly2->exps, poly2->bits,
281 poly2->length, ctx->minfo);
282 }
283
284 if (exp_bits > poly3->bits)
285 {
286 free3 = 1;
287 exp3 = (ulong *) flint_malloc(N*poly3->length*sizeof(ulong));
288 mpoly_repack_monomials(exp3, exp_bits, poly3->exps, poly3->bits,
289 poly3->length, ctx->minfo);
290 }
291
292 /* check divisor leading monomial is at most that of the dividend */
293 if (mpoly_monomial_lt(exp2, exp3, N, cmpmask))
294 {
295 fq_zech_mpoly_set(r, poly2, ctx);
296 fq_zech_mpoly_zero(q, ctx);
297 goto cleanup3;
298 }
299
300 /* take care of aliasing */
301 if (q == poly2 || q == poly3)
302 {
303 fq_zech_mpoly_init2(temp1, poly2->length/poly3->length + 1, ctx);
304 fq_zech_mpoly_fit_bits(temp1, exp_bits, ctx);
305 temp1->bits = exp_bits;
306 tq = temp1;
307 }
308 else
309 {
310 fq_zech_mpoly_fit_length(q, poly2->length/poly3->length + 1, ctx);
311 fq_zech_mpoly_fit_bits(q, exp_bits, ctx);
312 q->bits = exp_bits;
313 tq = q;
314 }
315
316 if (r == poly2 || r == poly3)
317 {
318 fq_zech_mpoly_init2(temp2, poly3->length, ctx);
319 fq_zech_mpoly_fit_bits(temp2, exp_bits, ctx);
320 temp2->bits = exp_bits;
321 tr = temp2;
322 }
323 else
324 {
325 fq_zech_mpoly_fit_length(r, poly3->length, ctx);
326 fq_zech_mpoly_fit_bits(r, exp_bits, ctx);
327 r->bits = exp_bits;
328 tr = r;
329 }
330
331 /* do division with remainder */
332 while ((lenq = _fq_zech_mpoly_divrem_monagan_pearce(&lenr, &tq->coeffs, &tq->exps,
333 &tq->alloc, &tr->coeffs, &tr->exps, &tr->alloc, poly2->coeffs, exp2,
334 poly2->length, poly3->coeffs, exp3, poly3->length, exp_bits,
335 N, cmpmask, ctx->fqctx)) == 0
336 && lenr == 0)
337 {
338 ulong * old_exp2 = exp2, * old_exp3 = exp3;
339 slong old_exp_bits = exp_bits;
340
341 exp_bits = mpoly_fix_bits(exp_bits + 1, ctx->minfo);
342
343 N = mpoly_words_per_exp(exp_bits, ctx->minfo);
344 cmpmask = (ulong *) flint_realloc(cmpmask, N*sizeof(ulong));
345 mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
346
347 exp2 = (ulong *) flint_malloc(N*poly2->length*sizeof(ulong));
348 mpoly_repack_monomials(exp2, exp_bits, old_exp2, old_exp_bits,
349 poly2->length, ctx->minfo);
350
351 exp3 = (ulong *) flint_malloc(N*poly3->length*sizeof(ulong));
352 mpoly_repack_monomials(exp3, exp_bits, old_exp3, old_exp_bits,
353 poly3->length, ctx->minfo);
354
355 if (free2)
356 flint_free(old_exp2);
357
358 if (free3)
359 flint_free(old_exp3);
360
361 free2 = free3 = 1;
362
363 fq_zech_mpoly_fit_bits(tq, exp_bits, ctx);
364 tq->bits = exp_bits;
365
366 fq_zech_mpoly_fit_bits(tr, exp_bits, ctx);
367 tr->bits = exp_bits;
368 }
369
370 /* deal with aliasing */
371 if (q == poly2 || q == poly3)
372 {
373 fq_zech_mpoly_swap(temp1, q, ctx);
374 fq_zech_mpoly_clear(temp1, ctx);
375 }
376
377 if (r == poly2 || r == poly3)
378 {
379 fq_zech_mpoly_swap(temp2, r, ctx);
380 fq_zech_mpoly_clear(temp2, ctx);
381 }
382
383 _fq_zech_mpoly_set_length(q, lenq, ctx);
384 _fq_zech_mpoly_set_length(r, lenr, ctx);
385
386 cleanup3:
387
388 if (free2)
389 flint_free(exp2);
390
391 if (free3)
392 flint_free(exp3);
393
394 flint_free(cmpmask);
395 }
396