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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "nmod_mpoly.h"
13 #include "nmod_mpoly_factor.h"
14
nmod_mpolyu_init(nmod_mpolyu_t A,flint_bitcnt_t bits,const nmod_mpoly_ctx_t ctx)15 void nmod_mpolyu_init(nmod_mpolyu_t A, flint_bitcnt_t bits,
16 const nmod_mpoly_ctx_t ctx)
17 {
18 A->coeffs = NULL;
19 A->exps = NULL;
20 A->alloc = 0;
21 A->length = 0;
22 A->bits = bits;
23 }
24
25
nmod_mpolyu_clear(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx)26 void nmod_mpolyu_clear(nmod_mpolyu_t A, const nmod_mpoly_ctx_t uctx)
27 {
28 slong i;
29 for (i = 0; i < A->alloc; i++)
30 nmod_mpoly_clear(A->coeffs + i, uctx);
31 flint_free(A->coeffs);
32 flint_free(A->exps);
33 }
34
35
nmod_mpolyu_swap(nmod_mpolyu_t A,nmod_mpolyu_t B,const nmod_mpoly_ctx_t uctx)36 void nmod_mpolyu_swap(nmod_mpolyu_t A, nmod_mpolyu_t B,
37 const nmod_mpoly_ctx_t uctx)
38 {
39 nmod_mpolyu_struct t = *A;
40 *A = *B;
41 *B = t;
42 }
43
nmod_mpolyu_zero(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx)44 void nmod_mpolyu_zero(nmod_mpolyu_t A, const nmod_mpoly_ctx_t uctx)
45 {
46 A->length = 0;
47 }
48
nmod_mpolyu_is_one(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx)49 int nmod_mpolyu_is_one(nmod_mpolyu_t A, const nmod_mpoly_ctx_t uctx)
50 {
51 if (A->length != 1 || A->exps[0] != UWORD(0))
52 return 0;
53
54 return nmod_mpoly_is_one(A->coeffs + 0, uctx);
55 }
56
nmod_mpolyu_print_pretty(const nmod_mpolyu_t poly,const char ** x,const nmod_mpoly_ctx_t ctx)57 void nmod_mpolyu_print_pretty(const nmod_mpolyu_t poly,
58 const char ** x, const nmod_mpoly_ctx_t ctx)
59 {
60 slong i;
61 if (poly->length == 0)
62 flint_printf("0");
63 for (i = 0; i < poly->length; i++)
64 {
65 if (i != 0)
66 flint_printf(" + ");
67 flint_printf("(");
68 nmod_mpoly_print_pretty(poly->coeffs + i,x,ctx);
69 flint_printf(")*X^%wd",poly->exps[i]);
70 }
71 }
72
nmod_mpolyu_fit_length(nmod_mpolyu_t A,slong length,const nmod_mpoly_ctx_t uctx)73 void nmod_mpolyu_fit_length(nmod_mpolyu_t A, slong length,
74 const nmod_mpoly_ctx_t uctx)
75 {
76 slong i;
77 slong old_alloc = A->alloc;
78 slong new_alloc = FLINT_MAX(length, 2*A->alloc);
79
80 if (length > old_alloc)
81 {
82 A->exps = (ulong *) flint_realloc(A->exps, new_alloc*sizeof(ulong));
83 A->coeffs = (nmod_mpoly_struct *) flint_realloc(A->coeffs,
84 new_alloc*sizeof(nmod_mpoly_struct));
85
86 for (i = old_alloc; i < new_alloc; i++)
87 nmod_mpoly_init3(A->coeffs + i, 0, A->bits, uctx);
88
89 A->alloc = new_alloc;
90 }
91 }
92
nmod_mpolyu_one(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx)93 void nmod_mpolyu_one(nmod_mpolyu_t A, const nmod_mpoly_ctx_t uctx)
94 {
95 nmod_mpolyu_fit_length(A, WORD(1), uctx);
96 A->exps[0] = UWORD(0);
97 nmod_mpoly_one(A->coeffs + 0, uctx);
98 A->length = WORD(1);
99 }
100
101
nmod_mpolyu_degrees_si(slong * degs,const nmod_mpolyu_t A,const nmod_mpoly_ctx_t ctx)102 void nmod_mpolyu_degrees_si(
103 slong * degs,
104 const nmod_mpolyu_t A,
105 const nmod_mpoly_ctx_t ctx)
106 {
107 slong i, j;
108 flint_bitcnt_t bits = A->bits;
109 slong N = mpoly_words_per_exp_sp(bits, ctx->minfo);
110 ulong * pmax, mask;
111 TMP_INIT;
112
113 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
114
115 if (A->length < 1)
116 {
117 for (j = 0; j < ctx->minfo->nvars; j++)
118 degs[j] = -1;
119 }
120
121 TMP_START;
122
123 mask = mpoly_overflow_mask_sp(bits);
124
125 pmax = (ulong *) TMP_ALLOC(N*sizeof(ulong));
126 mpoly_monomial_zero(pmax, N);
127
128 for (i = 0; i < A->length; i++)
129 {
130 ulong * Aiexps = A->coeffs[i].exps;
131 FLINT_ASSERT(A->coeffs[i].bits == bits);
132
133 for (j = 0; j < A->coeffs[i].length; j++)
134 mpoly_monomial_max(pmax, pmax, Aiexps + N*j, bits, N, mask);
135 }
136
137 mpoly_unpack_vec_ui((ulong *) degs, pmax, bits, ctx->minfo->nvars, 1);
138
139 for (i = 0; i < ctx->minfo->nvars/2; i++)
140 SLONG_SWAP(degs[i], degs[ctx->minfo->nvars - i - 1]);
141
142 TMP_END;
143 }
144
nmod_mpolyu_repack_bits_inplace(nmod_mpolyu_t A,flint_bitcnt_t bits,const nmod_mpoly_ctx_t ctx)145 void nmod_mpolyu_repack_bits_inplace(
146 nmod_mpolyu_t A,
147 flint_bitcnt_t bits,
148 const nmod_mpoly_ctx_t ctx)
149 {
150 slong i;
151
152 if (bits == A->bits)
153 return;
154
155 A->bits = bits;
156
157 for (i = 0; i < A->alloc; i++)
158 nmod_mpoly_repack_bits_inplace(A->coeffs + i, bits, ctx);
159 }
160
161 /* if the coefficient doesn't exist, a new one is created (and set to zero) */
_nmod_mpolyu_get_coeff(nmod_mpolyu_t A,ulong pow,const nmod_mpoly_ctx_t uctx)162 nmod_mpoly_struct * _nmod_mpolyu_get_coeff(nmod_mpolyu_t A,
163 ulong pow, const nmod_mpoly_ctx_t uctx)
164 {
165 slong i, j;
166 nmod_mpoly_struct * xk;
167
168 for (i = 0; i < A->length && A->exps[i] >= pow; i++)
169 {
170 if (A->exps[i] == pow)
171 {
172 return A->coeffs + i;
173 }
174 }
175
176 nmod_mpolyu_fit_length(A, A->length + 1, uctx);
177
178 for (j = A->length; j > i; j--)
179 {
180 A->exps[j] = A->exps[j - 1];
181 nmod_mpoly_swap(A->coeffs + j, A->coeffs + j - 1, uctx);
182 }
183
184 A->length++;
185 A->exps[i] = pow;
186 xk = A->coeffs + i;
187 xk->length = 0;
188 FLINT_ASSERT(xk->bits == A->bits);
189
190 return xk;
191 }
192
193
194 typedef struct
195 {
196 volatile slong index;
197 #if FLINT_USES_PTHREAD
198 pthread_mutex_t mutex;
199 #endif
200 slong length;
201 nmod_mpoly_struct * coeffs;
202 const nmod_mpoly_ctx_struct * ctx;
203 }
204 _sort_arg_struct;
205
206 typedef _sort_arg_struct _sort_arg_t[1];
207
208
_worker_sort(void * varg)209 static void _worker_sort(void * varg)
210 {
211 _sort_arg_struct * arg = (_sort_arg_struct *) varg;
212 slong i;
213
214 get_next_index:
215
216 #if FLINT_USES_PTHREAD
217 pthread_mutex_lock(&arg->mutex);
218 #endif
219 i = arg->index;
220 arg->index++;
221 #if FLINT_USES_PTHREAD
222 pthread_mutex_unlock(&arg->mutex);
223 #endif
224
225 if (i >= arg->length)
226 goto cleanup;
227
228 nmod_mpoly_sort_terms(arg->coeffs + i, arg->ctx);
229
230 goto get_next_index;
231
232 cleanup:
233
234 return;
235 }
236
237
238 /** 1 variables *************************************/
239
240
241 /*
242 Convert B to A using the variable permutation perm.
243 The uctx (m vars) should be the context of the coefficients of A.
244 The ctx (n vars) should be the context of B.
245
246 operation on each term:
247
248 for 0 <= k < m + 1
249 l = perm[k]
250 Aexp[k] = (Bexp[l] - shift[l])/stride[l]
251
252 the most significant main variable uses k = 0
253 the coefficients of A use variables k = 1 ... m
254 */
nmod_mpoly_to_mpolyu_perm_deflate_threaded_pool(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx,const slong * perm,const ulong * shift,const ulong * stride,const thread_pool_handle * handles,slong num_handles)255 void nmod_mpoly_to_mpolyu_perm_deflate_threaded_pool(
256 nmod_mpolyu_t A,
257 const nmod_mpoly_ctx_t uctx,
258 const nmod_mpoly_t B,
259 const nmod_mpoly_ctx_t ctx,
260 const slong * perm,
261 const ulong * shift,
262 const ulong * stride,
263 const thread_pool_handle * handles,
264 slong num_handles)
265 {
266 slong i, j, k, l;
267 slong n = ctx->minfo->nvars;
268 slong m = uctx->minfo->nvars;
269 slong NA, NB;
270 ulong * uexps;
271 ulong * Bexps;
272 nmod_mpoly_struct * Ac;
273 TMP_INIT;
274
275 FLINT_ASSERT(A->bits <= FLINT_BITS);
276 FLINT_ASSERT(B->bits <= FLINT_BITS);
277 FLINT_ASSERT(m + 1 <= n);
278
279 TMP_START;
280
281 uexps = (ulong *) TMP_ALLOC((m + 1)*sizeof(ulong));
282 Bexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
283
284 nmod_mpolyu_zero(A, uctx);
285
286 NA = mpoly_words_per_exp(A->bits, uctx->minfo);
287 NB = mpoly_words_per_exp(B->bits, ctx->minfo);
288
289 for (j = 0; j < B->length; j++)
290 {
291 mpoly_get_monomial_ui(Bexps, B->exps + NB*j, B->bits, ctx->minfo);
292 for (k = 0; k < m + 1; k++)
293 {
294 l = perm[k];
295 FLINT_ASSERT(stride[l] != UWORD(0));
296 FLINT_ASSERT(((Bexps[l] - shift[l]) % stride[l]) == UWORD(0));
297 uexps[k] = (Bexps[l] - shift[l]) / stride[l];
298 }
299 Ac = _nmod_mpolyu_get_coeff(A, uexps[0], uctx);
300 FLINT_ASSERT(Ac->bits == A->bits);
301
302 nmod_mpoly_fit_length(Ac, Ac->length + 1, uctx);
303 Ac->coeffs[Ac->length] = B->coeffs[j];
304 mpoly_set_monomial_ui(Ac->exps + NA*Ac->length, uexps + 1, A->bits, uctx->minfo);
305 Ac->length++;
306 }
307
308 if (num_handles > 0)
309 {
310 _sort_arg_t arg;
311
312 #if FLINT_USES_PTHREAD
313 pthread_mutex_init(&arg->mutex, NULL);
314 #endif
315 arg->index = 0;
316 arg->coeffs = A->coeffs;
317 arg->length = A->length;
318 arg->ctx = uctx;
319
320 for (i = 0; i < num_handles; i++)
321 {
322 thread_pool_wake(global_thread_pool, handles[i], 0, _worker_sort, arg);
323 }
324 _worker_sort(arg);
325 for (i = 0; i < num_handles; i++)
326 {
327 thread_pool_wait(global_thread_pool, handles[i]);
328 }
329
330 #if FLINT_USES_PTHREAD
331 pthread_mutex_destroy(&arg->mutex);
332 #endif
333 }
334 else
335 {
336 for (i = 0; i < A->length; i++)
337 {
338 nmod_mpoly_sort_terms(A->coeffs + i, uctx);
339 }
340 }
341
342 TMP_END;
343 }
344
345 /*
346 Convert B to A using the variable permutation vector perm.
347 This function inverts nmod_mpoly_to_mpolyu_perm_deflate.
348 A will be constructed with bits = Abits.
349
350 operation on each term:
351
352 for 0 <= l < n
353 Aexp[l] = shift[l]
354
355 for 0 <= k < m + 1
356 l = perm[k]
357 Aexp[l] += scale[l]*Bexp[k]
358 */
nmod_mpoly_from_mpolyu_perm_inflate(nmod_mpoly_t A,flint_bitcnt_t Abits,const nmod_mpoly_ctx_t ctx,const nmod_mpolyu_t B,const nmod_mpoly_ctx_t uctx,const slong * perm,const ulong * shift,const ulong * stride)359 void nmod_mpoly_from_mpolyu_perm_inflate(
360 nmod_mpoly_t A, flint_bitcnt_t Abits,
361 const nmod_mpoly_ctx_t ctx,
362 const nmod_mpolyu_t B,
363 const nmod_mpoly_ctx_t uctx,
364 const slong * perm,
365 const ulong * shift,
366 const ulong * stride)
367 {
368 slong n = ctx->minfo->nvars;
369 slong m = uctx->minfo->nvars;
370 slong i, j, k, l;
371 slong NA, NB;
372 slong Alen;
373 mp_limb_t * Acoeff;
374 ulong * Aexp;
375 ulong * uexps;
376 ulong * Aexps;
377 TMP_INIT;
378
379 FLINT_ASSERT(B->length > 0);
380 FLINT_ASSERT(Abits <= FLINT_BITS);
381 FLINT_ASSERT(B->bits <= FLINT_BITS);
382 FLINT_ASSERT(m + 1 <= n);
383 TMP_START;
384
385 uexps = (ulong *) TMP_ALLOC((m + 1)*sizeof(ulong));
386 Aexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
387
388 NA = mpoly_words_per_exp(Abits, ctx->minfo);
389 NB = mpoly_words_per_exp(B->bits, uctx->minfo);
390
391 nmod_mpoly_fit_length_reset_bits(A, B->length, Abits, ctx);
392
393 Acoeff = A->coeffs;
394 Aexp = A->exps;
395 Alen = 0;
396 for (i = 0; i < B->length; i++)
397 {
398 nmod_mpoly_struct * Bc = B->coeffs + i;
399 FLINT_ASSERT(Bc->bits == B->bits);
400
401 _nmod_mpoly_fit_length(&Acoeff, &A->coeffs_alloc,
402 &Aexp, &A->exps_alloc, NA, Alen + Bc->length);
403
404 for (j = 0; j < Bc->length; j++)
405 {
406 Acoeff[Alen + j] = Bc->coeffs[j];
407 mpoly_get_monomial_ui(uexps + 1, Bc->exps + NB*j, Bc->bits, uctx->minfo);
408 uexps[0] = B->exps[i];
409 for (l = 0; l < n; l++)
410 {
411 Aexps[l] = shift[l];
412 }
413 for (k = 0; k < m + 1; k++)
414 {
415 l = perm[k];
416 Aexps[l] += stride[l]*uexps[k];
417 }
418 mpoly_set_monomial_ui(Aexp + NA*(Alen + j), Aexps, Abits, ctx->minfo);
419 }
420 Alen += Bc->length;
421 }
422 A->coeffs = Acoeff;
423 A->exps = Aexp;
424 _nmod_mpoly_set_length(A, Alen, ctx);
425
426 TMP_END;
427
428 nmod_mpoly_sort_terms(A, ctx);
429 }
430
431
432 /** 0 variables *************************************/
433
434
nmod_mpoly_to_mpolyl_perm_deflate(nmod_mpoly_t A,const nmod_mpoly_ctx_t lctx,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx,const slong * perm,const ulong * shift,const ulong * stride)435 void nmod_mpoly_to_mpolyl_perm_deflate(
436 nmod_mpoly_t A,
437 const nmod_mpoly_ctx_t lctx,
438 const nmod_mpoly_t B,
439 const nmod_mpoly_ctx_t ctx,
440 const slong * perm,
441 const ulong * shift,
442 const ulong * stride)
443 {
444 slong j, k, l;
445 slong m = lctx->minfo->nvars;
446 slong n = ctx->minfo->nvars;
447 slong NA, NB;
448 ulong * lexps;
449 ulong * Bexps;
450 TMP_INIT;
451
452 FLINT_ASSERT(A->bits <= FLINT_BITS);
453 FLINT_ASSERT(B->bits <= FLINT_BITS);
454 FLINT_ASSERT(m <= n);
455 TMP_START;
456
457 nmod_mpoly_fit_length(A, B->length, ctx);
458 A->length = B->length;
459
460 lexps = (ulong *) TMP_ALLOC(m*sizeof(ulong));
461 Bexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
462 NA = mpoly_words_per_exp(A->bits, lctx->minfo);
463 NB = mpoly_words_per_exp(B->bits, ctx->minfo);
464
465 for (j = 0; j < B->length; j++)
466 {
467 A->coeffs[j] = B->coeffs[j];
468
469 mpoly_get_monomial_ui(Bexps, B->exps + NB*j, B->bits, ctx->minfo);
470 for (k = 0; k < m; k++)
471 {
472 l = perm[k];
473 if (stride[l] == 1)
474 {
475 lexps[k] = (Bexps[l] - shift[l]);
476 }
477 else
478 {
479 FLINT_ASSERT(stride[l] != 0);
480 FLINT_ASSERT(((Bexps[l] - shift[l]) % stride[l]) == 0);
481 lexps[k] = (Bexps[l] - shift[l]) / stride[l];
482 }
483 }
484
485 mpoly_set_monomial_ui(A->exps + NA*j, lexps, A->bits, lctx->minfo);
486 }
487
488 TMP_END;
489
490 nmod_mpoly_sort_terms(A, lctx);
491
492 FLINT_ASSERT(nmod_mpoly_is_canonical(A, lctx));
493 }
494
495
496 /*
497 Convert B to A using the variable permutation vector perm.
498 A must be constructed with bits = Abits.
499
500 operation on each term:
501
502 for 0 <= l < n
503 Aexp[l] = shift[l]
504
505 for 0 <= k < m + 2
506 l = perm[k]
507 Aexp[l] += scale[l]*Bexp[k]
508 */
nmod_mpoly_from_mpolyl_perm_inflate(nmod_mpoly_t A,flint_bitcnt_t Abits,const nmod_mpoly_ctx_t ctx,const nmod_mpoly_t B,const nmod_mpoly_ctx_t lctx,const slong * perm,const ulong * shift,const ulong * stride)509 void nmod_mpoly_from_mpolyl_perm_inflate(
510 nmod_mpoly_t A,
511 flint_bitcnt_t Abits,
512 const nmod_mpoly_ctx_t ctx,
513 const nmod_mpoly_t B,
514 const nmod_mpoly_ctx_t lctx,
515 const slong * perm,
516 const ulong * shift,
517 const ulong * stride)
518 {
519 slong n = ctx->minfo->nvars;
520 slong m = lctx->minfo->nvars;
521 slong i, k, l;
522 slong NA, NB;
523 ulong * Bexps;
524 ulong * Aexps;
525 TMP_INIT;
526
527 FLINT_ASSERT(B->length > 0);
528 FLINT_ASSERT(Abits <= FLINT_BITS);
529 FLINT_ASSERT(B->bits <= FLINT_BITS);
530 FLINT_ASSERT(m <= n);
531 TMP_START;
532
533 Bexps = (ulong *) TMP_ALLOC(m*sizeof(ulong));
534 Aexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
535
536 NA = mpoly_words_per_exp(Abits, ctx->minfo);
537 NB = mpoly_words_per_exp(B->bits, lctx->minfo);
538
539 nmod_mpoly_fit_length_reset_bits(A, B->length, Abits, ctx);
540 A->length = B->length;
541
542 for (i = 0; i < B->length; i++)
543 {
544 A->coeffs[i] = B->coeffs[i];
545 mpoly_get_monomial_ui(Bexps, B->exps + NB*i, B->bits, lctx->minfo);
546
547 for (l = 0; l < n; l++)
548 {
549 Aexps[l] = shift[l];
550 }
551 for (k = 0; k < m; k++)
552 {
553 l = perm[k];
554 Aexps[l] += stride[l]*Bexps[k];
555 }
556
557 mpoly_set_monomial_ui(A->exps + NA*i, Aexps, Abits, ctx->minfo);
558 }
559
560 TMP_END;
561
562 nmod_mpoly_sort_terms(A, ctx);
563
564 FLINT_ASSERT(nmod_mpoly_is_canonical(A, ctx));
565 }
566
567
568
569 /** 2 variables *************************************/
570
571 /*
572 Convert B to A using the variable permutation perm.
573 The uctx should be the context of the coefficients of A.
574 The ctx should be the context of B.
575
576 operation on each term:
577
578 for 0 <= k < m + 2
579 l = perm[k]
580 Aexp[k] = (Bexp[l] - shift[l])/stride[l]
581
582 the most significant main variable uses Aexp[0]
583 the least significant main variable uses Aexp[1]
584 the coefficients of A use variables Aexp[2], ..., Aexp[m + 1]
585 maxexps if it exists is supposed to be a degree bound on B
586 */
nmod_mpoly_to_mpolyuu_perm_deflate_threaded_pool(nmod_mpolyu_t A,const nmod_mpoly_ctx_t uctx,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx,const slong * perm,const ulong * shift,const ulong * stride,const ulong * maxexps,const thread_pool_handle * handles,slong num_handles)587 void nmod_mpoly_to_mpolyuu_perm_deflate_threaded_pool(
588 nmod_mpolyu_t A,
589 const nmod_mpoly_ctx_t uctx,
590 const nmod_mpoly_t B,
591 const nmod_mpoly_ctx_t ctx,
592 const slong * perm,
593 const ulong * shift,
594 const ulong * stride,
595 const ulong * maxexps, /* nullptr is ok */
596 const thread_pool_handle * handles,
597 slong num_handles)
598 {
599 slong i, j, k, l;
600 slong n = ctx->minfo->nvars;
601 slong m = uctx->minfo->nvars;
602 slong NA, NB;
603 ulong * uexps;
604 ulong * Bexps;
605 nmod_mpoly_struct * Ac;
606 TMP_INIT;
607
608 FLINT_ASSERT(A->bits <= FLINT_BITS);
609 FLINT_ASSERT(B->bits <= FLINT_BITS);
610 FLINT_ASSERT(m + 2 <= n);
611
612 nmod_mpolyu_zero(A, uctx);
613
614 {
615 TMP_START;
616
617 uexps = (ulong *) TMP_ALLOC((m + 2)*sizeof(ulong));
618 Bexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
619 NA = mpoly_words_per_exp(A->bits, uctx->minfo);
620 NB = mpoly_words_per_exp(B->bits, ctx->minfo);
621
622 for (j = 0; j < B->length; j++)
623 {
624 mpoly_get_monomial_ui(Bexps, B->exps + NB*j, B->bits, ctx->minfo);
625 for (k = 0; k < m + 2; k++)
626 {
627 l = perm[k];
628 if (stride[l] == 1)
629 {
630 uexps[k] = (Bexps[l] - shift[l]);
631 }
632 else
633 {
634 FLINT_ASSERT(stride[l] != 0);
635 FLINT_ASSERT(((Bexps[l] - shift[l]) % stride[l]) == 0);
636 uexps[k] = (Bexps[l] - shift[l]) / stride[l];
637 }
638 }
639 FLINT_ASSERT(FLINT_BIT_COUNT(uexps[0]) < FLINT_BITS/2);
640 FLINT_ASSERT(FLINT_BIT_COUNT(uexps[1]) < FLINT_BITS/2);
641 Ac = _nmod_mpolyu_get_coeff(A, (uexps[0] << (FLINT_BITS/2)) + uexps[1], uctx);
642 FLINT_ASSERT(Ac->bits == A->bits);
643
644 nmod_mpoly_fit_length(Ac, Ac->length + 1, uctx);
645 Ac->coeffs[Ac->length] = B->coeffs[j];
646 mpoly_set_monomial_ui(Ac->exps + NA*Ac->length, uexps + 2, A->bits, uctx->minfo);
647 Ac->length++;
648 }
649
650
651 for (i = 0; i < A->length; i++)
652 {
653 nmod_mpoly_sort_terms(A->coeffs + i, uctx);
654 }
655
656 TMP_END;
657 }
658 }
659
660
661 /*
662 Convert B to A using the variable permutation vector perm.
663 A must be constructed with bits = Abits.
664
665 operation on each term:
666
667 for 0 <= l < n
668 Aexp[l] = shift[l]
669
670 for 0 <= k < m + 2
671 l = perm[k]
672 Aexp[l] += scale[l]*Bexp[k]
673 */
nmod_mpoly_from_mpolyuu_perm_inflate(nmod_mpoly_t A,flint_bitcnt_t Abits,const nmod_mpoly_ctx_t ctx,const nmod_mpolyu_t B,const nmod_mpoly_ctx_t uctx,const slong * perm,const ulong * shift,const ulong * stride)674 void nmod_mpoly_from_mpolyuu_perm_inflate( /* only for 2 main vars */
675 nmod_mpoly_t A,
676 flint_bitcnt_t Abits,
677 const nmod_mpoly_ctx_t ctx,
678 const nmod_mpolyu_t B,
679 const nmod_mpoly_ctx_t uctx,
680 const slong * perm,
681 const ulong * shift,
682 const ulong * stride)
683 {
684 slong n = ctx->minfo->nvars;
685 slong m = uctx->minfo->nvars;
686 slong i, j, k, l;
687 slong NA, NB;
688 slong Alen;
689 mp_limb_t * Acoeff;
690 ulong * Aexp;
691 ulong * uexps;
692 ulong * Aexps;
693 TMP_INIT;
694
695 FLINT_ASSERT(B->length > 0);
696 FLINT_ASSERT(Abits <= FLINT_BITS);
697 FLINT_ASSERT(B->bits <= FLINT_BITS);
698 FLINT_ASSERT(m + 2 <= n);
699 TMP_START;
700
701 uexps = (ulong *) TMP_ALLOC((m + 2)*sizeof(ulong));
702 Aexps = (ulong *) TMP_ALLOC(n*sizeof(ulong));
703
704 NA = mpoly_words_per_exp(Abits, ctx->minfo);
705 NB = mpoly_words_per_exp(B->bits, uctx->minfo);
706
707 nmod_mpoly_fit_length_reset_bits(A, B->length, Abits, ctx);
708
709 Acoeff = A->coeffs;
710 Aexp = A->exps;
711 Alen = 0;
712 for (i = 0; i < B->length; i++)
713 {
714 nmod_mpoly_struct * Bc = B->coeffs + i;
715 FLINT_ASSERT(Bc->bits == B->bits);
716
717 _nmod_mpoly_fit_length(&Acoeff, &A->coeffs_alloc,
718 &Aexp, &A->exps_alloc, NA, Alen + Bc->length);
719
720 for (j = 0; j < Bc->length; j++)
721 {
722 Acoeff[Alen + j] = Bc->coeffs[j];
723 mpoly_get_monomial_ui(uexps + 2, Bc->exps + NB*j, Bc->bits, uctx->minfo);
724 uexps[0] = B->exps[i] >> (FLINT_BITS/2);
725 uexps[1] = B->exps[i] & ((-UWORD(1)) >> (FLINT_BITS - FLINT_BITS/2));
726 for (l = 0; l < n; l++)
727 {
728 Aexps[l] = shift[l];
729 }
730 for (k = 0; k < m + 2; k++)
731 {
732 l = perm[k];
733 Aexps[l] += stride[l]*uexps[k];
734 }
735 mpoly_set_monomial_ui(Aexp + NA*(Alen + j), Aexps, Abits, ctx->minfo);
736 }
737 Alen += Bc->length;
738 }
739 A->coeffs = Acoeff;
740 A->exps = Aexp;
741 _nmod_mpoly_set_length(A, Alen, ctx);
742
743 nmod_mpoly_sort_terms(A, ctx);
744 TMP_END;
745 }
746
747
748
nmod_mpolyu_shift_right(nmod_mpolyu_t A,ulong s)749 void nmod_mpolyu_shift_right(nmod_mpolyu_t A, ulong s)
750 {
751 slong i;
752 for (i = 0; i < A->length; i++)
753 {
754 FLINT_ASSERT(A->exps[i] >= s);
755 A->exps[i] -= s;
756 }
757 }
758
nmod_mpolyu_shift_left(nmod_mpolyu_t A,ulong s)759 void nmod_mpolyu_shift_left(nmod_mpolyu_t A, ulong s)
760 {
761 slong i;
762 for (i = 0; i < A->length; i++)
763 {
764 A->exps[i] += s;
765 }
766 }
767
nmod_mpolyu_scalar_mul_nmod(nmod_mpolyu_t A,mp_limb_t c,const nmod_mpoly_ctx_t ctx)768 void nmod_mpolyu_scalar_mul_nmod(nmod_mpolyu_t A, mp_limb_t c,
769 const nmod_mpoly_ctx_t ctx)
770 {
771 slong i, j;
772
773 for (i = 0; i < A->length; i++)
774 {
775 for (j = 0; j < (A->coeffs + i)->length; j++)
776 {
777 A->coeffs[i].coeffs[j] = nmod_mul(A->coeffs[i].coeffs[j], c, ctx->mod);
778 }
779 }
780 }
781
782
nmod_mpolyu_set(nmod_mpolyu_t A,const nmod_mpolyu_t B,const nmod_mpoly_ctx_t uctx)783 void nmod_mpolyu_set(nmod_mpolyu_t A, const nmod_mpolyu_t B,
784 const nmod_mpoly_ctx_t uctx)
785 {
786 slong i;
787 nmod_mpoly_struct * Acoeff, * Bcoeff;
788 ulong * Aexp, * Bexp;
789 slong Alen, Blen;
790
791 Alen = 0;
792 Blen = B->length;
793 nmod_mpolyu_fit_length(A, Blen, uctx);
794 Acoeff = A->coeffs;
795 Bcoeff = B->coeffs;
796 Aexp = A->exps;
797 Bexp = B->exps;
798
799 for (i = 0; i < Blen; i++)
800 {
801 nmod_mpoly_set(Acoeff + Alen, Bcoeff + i, uctx);
802 Aexp[Alen++] = Bexp[i];
803 }
804 Alen = Blen;
805
806 /* demote remaining coefficients */
807 for (i = Alen; i < A->length; i++)
808 {
809 nmod_mpoly_clear(Acoeff + i, uctx);
810 nmod_mpoly_init(Acoeff + i, uctx);
811 }
812 A->length = Alen;
813 }
814
815
816
nmod_mpoly_cvtfrom_poly_notmain(nmod_mpoly_t A,nmod_poly_t a,slong var,const nmod_mpoly_ctx_t ctx)817 void nmod_mpoly_cvtfrom_poly_notmain(nmod_mpoly_t A, nmod_poly_t a,
818 slong var, const nmod_mpoly_ctx_t ctx)
819 {
820 slong i;
821 slong k;
822 ulong * oneexp;
823 slong N;
824 TMP_INIT;
825 TMP_START;
826
827 FLINT_ASSERT(A->bits <= FLINT_BITS);
828
829 N = mpoly_words_per_exp_sp(A->bits, ctx->minfo);
830
831 oneexp = (ulong *)TMP_ALLOC(N*sizeof(ulong));
832 mpoly_gen_monomial_sp(oneexp, var, A->bits, ctx->minfo);
833
834 nmod_mpoly_fit_length(A, nmod_poly_length(a), ctx);
835
836 k = 0;
837 for (i = nmod_poly_length(a) - 1; i >= 0; i--)
838 {
839 mp_limb_t c = nmod_poly_get_coeff_ui(a, i);
840 if (c != UWORD(0))
841 {
842 A->coeffs[k] = c;
843 mpoly_monomial_mul_ui(A->exps + N*k, oneexp, N, i);
844 k++;
845 }
846 }
847 A->length = k;
848 TMP_END;
849 }
850 /*
851 Set "A" to "a" where "a" is a polynomial in a non-main variable "var"
852 */
nmod_mpolyu_cvtfrom_poly_notmain(nmod_mpolyu_t A,nmod_poly_t a,slong var,const nmod_mpoly_ctx_t ctx)853 void nmod_mpolyu_cvtfrom_poly_notmain(nmod_mpolyu_t A, nmod_poly_t a,
854 slong var, const nmod_mpoly_ctx_t ctx)
855 {
856 nmod_mpolyu_fit_length(A, 1, ctx);
857 A->exps[0] = 0;
858 nmod_mpoly_cvtfrom_poly_notmain(A->coeffs + 0, a, var, ctx);
859 A->length = !nmod_mpoly_is_zero(A->coeffs + 0, ctx);
860 }
861
862
863
864 /*
865 Assuming that "A" depends only on the main variable,
866 convert it to a poly "a".
867 */
nmod_mpolyu_cvtto_poly(nmod_poly_t a,nmod_mpolyu_t A,const nmod_mpoly_ctx_t ctx)868 void nmod_mpolyu_cvtto_poly(nmod_poly_t a, nmod_mpolyu_t A,
869 const nmod_mpoly_ctx_t ctx)
870 {
871 slong i;
872 nmod_poly_zero(a);
873 for (i = 0; i < A->length; i++)
874 {
875 FLINT_ASSERT((A->coeffs + i)->length == 1);
876 FLINT_ASSERT(mpoly_monomial_is_zero((A->coeffs + i)->exps,
877 mpoly_words_per_exp((A->coeffs + i)->bits, ctx->minfo)));
878 nmod_poly_set_coeff_ui(a, A->exps[i], (A->coeffs + i)->coeffs[0]);
879 }
880 }
881
882 /*
883 Convert a poly "a" to "A" in the main variable,
884 */
nmod_mpolyu_cvtfrom_poly(nmod_mpolyu_t A,nmod_poly_t a,const nmod_mpoly_ctx_t ctx)885 void nmod_mpolyu_cvtfrom_poly(nmod_mpolyu_t A, nmod_poly_t a,
886 const nmod_mpoly_ctx_t ctx)
887 {
888 slong i;
889 slong k;
890 slong N = mpoly_words_per_exp(A->bits, ctx->minfo);
891
892 nmod_mpolyu_zero(A, ctx);
893 k = 0;
894 for (i = nmod_poly_length(a) - 1; i >= 0; i--)
895 {
896 mp_limb_t c = nmod_poly_get_coeff_ui(a, i);
897 if (c != UWORD(0))
898 {
899 nmod_mpolyu_fit_length(A, k + 1, ctx);
900 A->exps[k] = i;
901 nmod_mpoly_fit_length_reset_bits(A->coeffs + k, 1, A->bits, ctx);
902 A->coeffs[k].coeffs[0] = c;
903 A->coeffs[k].length = 1;
904 mpoly_monomial_zero(A->coeffs[k].exps + N*0, N);
905 k++;
906 }
907 }
908 A->length = k;
909 }
910
911
nmod_mpolyu_msub(nmod_mpolyu_t R,nmod_mpolyu_t A,nmod_mpolyu_t B,nmod_mpoly_t c,slong e,const nmod_mpoly_ctx_t ctx)912 void nmod_mpolyu_msub(nmod_mpolyu_t R, nmod_mpolyu_t A, nmod_mpolyu_t B,
913 nmod_mpoly_t c, slong e, const nmod_mpoly_ctx_t ctx)
914 {
915 slong i, j, k;
916 nmod_mpoly_t T;
917
918 nmod_mpolyu_fit_length(R, A->length + B->length, ctx);
919
920 nmod_mpoly_init(T, ctx);
921
922 i = j = k = 0;
923 while (i < A->length || j < B->length)
924 {
925 if (i < A->length && (j >= B->length || A->exps[i] > B->exps[j] + e))
926 {
927 /* only A ok */
928 nmod_mpoly_set(R->coeffs + k, A->coeffs + i, ctx);
929 R->exps[k] = A->exps[i];
930 k++;
931 i++;
932 }
933 else if (j < B->length && (i >= A->length || B->exps[j] + e > A->exps[i]))
934 {
935 /* only B ok */
936 nmod_mpoly_mul(R->coeffs + k, B->coeffs + j, c, ctx);
937 nmod_mpoly_neg(R->coeffs + k, R->coeffs + k, ctx);
938 R->exps[k] = B->exps[j] + e;
939 k++;
940 j++;
941 }
942 else if (i < A->length && j < B->length && (A->exps[i] == B->exps[j] + e))
943 {
944 nmod_mpoly_mul(T, B->coeffs + j, c, ctx);
945 nmod_mpoly_sub(R->coeffs + k, A->coeffs + i, T, ctx);
946 R->exps[k] = A->exps[i];
947 k += !nmod_mpoly_is_zero(R->coeffs + k, ctx);
948 i++;
949 j++;
950 } else
951 {
952 FLINT_ASSERT(0);
953 }
954 }
955
956 nmod_mpoly_clear(T, ctx);
957 R->length = k;
958 }
959
960
nmod_mpolyu_divexact_mpoly_inplace(nmod_mpolyu_t A,nmod_mpoly_t c,const nmod_mpoly_ctx_t ctx)961 void nmod_mpolyu_divexact_mpoly_inplace(nmod_mpolyu_t A, nmod_mpoly_t c,
962 const nmod_mpoly_ctx_t ctx)
963 {
964 slong i;
965 flint_bitcnt_t bits = A->bits;
966 slong N = mpoly_words_per_exp(bits, ctx->minfo);
967 ulong * cmpmask;
968 nmod_mpoly_t t;
969 TMP_INIT;
970
971 FLINT_ASSERT(bits == c->bits);
972 FLINT_ASSERT(c->length > 0);
973
974 if (nmod_mpoly_is_ui(c, ctx))
975 {
976 if (c->coeffs[0] == 1)
977 return;
978
979 for (i = 0; i < A->length; i++)
980 {
981 nmod_mpoly_struct * Ai = A->coeffs + i;
982 _nmod_vec_scalar_mul_nmod(Ai->coeffs, Ai->coeffs, Ai->length,
983 nmod_inv(c->coeffs[0], ctx->mod), ctx->mod);
984 }
985
986 return;
987 }
988
989 nmod_mpoly_init3(t, 0, bits, ctx);
990
991 TMP_START;
992
993 cmpmask = (ulong *) TMP_ALLOC(N*sizeof(ulong));
994 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
995
996 for (i = A->length - 1; i >= 0; i--)
997 {
998 FLINT_ASSERT(A->coeffs[i].bits == bits);
999 _nmod_mpoly_divides_monagan_pearce(t, A->coeffs[i].coeffs,
1000 A->coeffs[i].exps, A->coeffs[i].length, c->coeffs, c->exps,
1001 c->length, bits, N, cmpmask, ctx->mod);
1002 nmod_mpoly_swap(A->coeffs + i, t, ctx);
1003 FLINT_ASSERT(A->coeffs[i].length > 0);
1004 }
1005
1006 TMP_END;
1007
1008 nmod_mpoly_clear(t, ctx);
1009 }
1010
1011
1012 /*
1013 A = B * c and preserve the bit packing
1014 */
nmod_mpolyu_mul_mpoly(nmod_mpolyu_t A,nmod_mpolyu_t B,nmod_mpoly_t c,const nmod_mpoly_ctx_t ctx)1015 void nmod_mpolyu_mul_mpoly(
1016 nmod_mpolyu_t A,
1017 nmod_mpolyu_t B,
1018 nmod_mpoly_t c,
1019 const nmod_mpoly_ctx_t ctx)
1020 {
1021 slong i;
1022 flint_bitcnt_t bits = B->bits;
1023 slong N = mpoly_words_per_exp(bits, ctx->minfo);
1024 ulong * cmpmask;
1025 TMP_INIT;
1026
1027 TMP_START;
1028
1029 bits = B->bits;
1030 FLINT_ASSERT(A->bits == B->bits);
1031 FLINT_ASSERT(A->bits == c->bits);
1032 FLINT_ASSERT(A != B);
1033
1034 nmod_mpolyu_fit_length(A, B->length, ctx);
1035
1036 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
1037 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
1038
1039 for (i = 0; i < B->length; i++)
1040 {
1041 nmod_mpoly_fit_length(A->coeffs + i,
1042 B->coeffs[i].length + c->length + 1, ctx);
1043 _nmod_mpoly_mul_johnson(A->coeffs + i, B->coeffs[i].coeffs,
1044 B->coeffs[i].exps, B->coeffs[i].length, c->coeffs,
1045 c->exps, c->length, bits, N, cmpmask, ctx->mod);
1046 FLINT_ASSERT(A->coeffs[i].length > 0);
1047 A->exps[i] = B->exps[i];
1048 }
1049 A->length = B->length;
1050 TMP_END;
1051 }
1052
1053
nmod_mpolyu_mul_mpoly_inplace(nmod_mpolyu_t A,nmod_mpoly_t c,const nmod_mpoly_ctx_t ctx)1054 void nmod_mpolyu_mul_mpoly_inplace(nmod_mpolyu_t A, nmod_mpoly_t c,
1055 const nmod_mpoly_ctx_t ctx)
1056 {
1057 slong i;
1058 flint_bitcnt_t bits = A->bits;
1059 slong N = mpoly_words_per_exp(bits, ctx->minfo);
1060 ulong * cmpmask;
1061 nmod_mpoly_t t;
1062 TMP_INIT;
1063
1064 FLINT_ASSERT(bits == c->bits);
1065 FLINT_ASSERT(c->length > 0);
1066
1067 if (nmod_mpoly_is_ui(c, ctx))
1068 {
1069 if (c->coeffs[0] == 1)
1070 return;
1071
1072 for (i = 0; i < A->length; i++)
1073 {
1074 nmod_mpoly_struct * Ai = A->coeffs + i;
1075 _nmod_vec_scalar_mul_nmod(Ai->coeffs, Ai->coeffs, Ai->length,
1076 c->coeffs[0], ctx->mod);
1077 }
1078
1079 return;
1080 }
1081
1082 nmod_mpoly_init3(t, 0, bits, ctx);
1083
1084 N = mpoly_words_per_exp(bits, ctx->minfo);
1085
1086 TMP_START;
1087
1088 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
1089 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
1090
1091 for (i = A->length - 1; i >= 0; i--)
1092 {
1093 FLINT_ASSERT(A->coeffs[i].bits == bits);
1094 _nmod_mpoly_mul_johnson(t, A->coeffs[i].coeffs,
1095 A->coeffs[i].exps, A->coeffs[i].length, c->coeffs, c->exps,
1096 c->length, bits, N, cmpmask, ctx->mod);
1097 nmod_mpoly_swap(A->coeffs + i, t, ctx);
1098 FLINT_ASSERT(A->coeffs[i].length > 0);
1099 }
1100
1101 TMP_END;
1102
1103 nmod_mpoly_clear(t, ctx);
1104 }
1105
1106
1107
nmod_mpolyu_content_mpoly(nmod_mpoly_t g,const nmod_mpolyu_t A,const nmod_mpoly_ctx_t ctx)1108 int nmod_mpolyu_content_mpoly(
1109 nmod_mpoly_t g,
1110 const nmod_mpolyu_t A,
1111 const nmod_mpoly_ctx_t ctx)
1112 {
1113 int success;
1114
1115 success = _nmod_mpoly_vec_content_mpoly(g, A->coeffs, A->length, ctx);
1116 if (!success)
1117 nmod_mpoly_zero(g, ctx);
1118
1119 nmod_mpoly_repack_bits_inplace(g, A->bits, ctx);
1120
1121 return success;
1122 }
1123