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 "nmod_mpoly_factor.h"
13 #include "fq_nmod_mpoly_factor.h"
14 #include "long_extras.h"
15
16
_deflate(nmod_mpoly_t A,slong tot_deg,const ulong * strides,const slong * perm,const nmod_mpoly_ctx_t ctx)17 static slong _deflate(
18 nmod_mpoly_t A,
19 slong tot_deg,
20 const ulong * strides,
21 const slong * perm,
22 const nmod_mpoly_ctx_t ctx)
23 {
24 slong i, j;
25 slong nvars = ctx->minfo->nvars;
26 flint_bitcnt_t bits = A->bits;
27 slong N = mpoly_words_per_exp_sp(bits, ctx->minfo);
28 ulong * texps, * sexps;
29 TMP_INIT;
30
31 for (j = 0; j < nvars; j++)
32 {
33 if (strides[j] != 1 || perm[j] != j)
34 goto do_it;
35 }
36
37 return tot_deg;
38
39 do_it:
40
41 TMP_START;
42
43 texps = (ulong *) TMP_ALLOC(2*nvars*sizeof(ulong));
44 sexps = texps + nvars;
45
46 tot_deg = 1;
47 for (i = 0; i < A->length; i++)
48 {
49 slong this_deg = 0;
50
51 mpoly_get_monomial_ui(texps, A->exps + N*i, bits, ctx->minfo);
52
53 for (j = 0; j < nvars; j++)
54 {
55 FLINT_ASSERT(0 == texps[j] % strides[j]);
56 texps[j] = texps[j]/strides[j];
57 }
58
59 for (j = 0; j < nvars; j++)
60 {
61 sexps[j] = texps[perm[j]];
62 this_deg += sexps[j];
63 }
64
65 tot_deg = FLINT_MAX(tot_deg, this_deg);
66
67 mpoly_set_monomial_ui(A->exps + N*i, sexps, bits, ctx->minfo);
68 }
69
70 TMP_END;
71
72 nmod_mpoly_sort_terms(A, ctx);
73 nmod_mpoly_make_monic(A, A, ctx);
74
75 return tot_deg;
76 }
77
78
_inflate(nmod_mpoly_t A,flint_bitcnt_t bits,const ulong * strides,const slong * perm,const nmod_mpoly_ctx_t ctx)79 static void _inflate(
80 nmod_mpoly_t A,
81 flint_bitcnt_t bits,
82 const ulong * strides,
83 const slong * perm,
84 const nmod_mpoly_ctx_t ctx)
85 {
86 slong i, j;
87 slong nvars = ctx->minfo->nvars;
88 slong N = mpoly_words_per_exp_sp(bits, ctx->minfo);
89 ulong * texps, * sexps;
90 TMP_INIT;
91
92 for (j = 0; j < nvars; j++)
93 {
94 if (strides[j] != 1 || perm[j] != j)
95 goto do_it;
96 }
97
98 return;
99
100 do_it:
101
102 nmod_mpoly_repack_bits_inplace(A, bits, ctx);
103
104 TMP_START;
105
106 texps = (ulong *) TMP_ALLOC(2*nvars*sizeof(ulong));
107 sexps = texps + nvars;
108
109 for (i = 0; i < A->length; i++)
110 {
111 mpoly_get_monomial_ui(sexps, A->exps + N*i, bits, ctx->minfo);
112
113 for (j = 0; j < nvars; j++)
114 texps[perm[j]] = sexps[j];
115
116 for (j = 0; j < nvars; j++)
117 texps[j] = texps[j]*strides[j];
118
119 mpoly_set_monomial_ui(A->exps + N*i, texps, bits, ctx->minfo);
120 }
121
122 TMP_END;
123
124 nmod_mpoly_sort_terms(A, ctx);
125 nmod_mpoly_make_monic(A, A, ctx);
126
127 return;
128 }
129
130
131 /* A has degree 2 wrt gen(0) */
_apply_quadratic(nmod_mpolyv_t Af,nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx)132 static int _apply_quadratic(
133 nmod_mpolyv_t Af,
134 nmod_mpoly_t A,
135 const nmod_mpoly_ctx_t ctx)
136 {
137 int success;
138 slong i, shift, off, N;
139 flint_bitcnt_t bits = A->bits;
140 ulong mask = (-UWORD(1)) >> (FLINT_BITS - bits);
141 nmod_mpoly_t a_mock, b_mock, c_mock;
142 nmod_mpoly_t t0, t1, t2, t3;
143
144 FLINT_ASSERT(A->length > 1 || A->coeffs[0] != 0);
145 FLINT_ASSERT(A->bits <= FLINT_BITS);
146 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
147
148 nmod_mpoly_init(t0, ctx);
149 nmod_mpoly_init(t1, ctx);
150 nmod_mpoly_init(t2, ctx);
151 nmod_mpoly_init(t3, ctx);
152
153 mpoly_gen_offset_shift_sp(&off, &shift, 0, bits, ctx->minfo);
154 N = mpoly_words_per_exp_sp(bits, ctx->minfo);
155
156 i = 0;
157 a_mock->exps = A->exps + N*i;
158 a_mock->coeffs = A->coeffs + i;
159 a_mock->bits = bits;
160 while (i < A->length && (mask & (A->exps[N*i + off] >> shift)) == 2)
161 i++;
162 a_mock->length = i;
163 a_mock->coeffs_alloc = a_mock->length;
164 a_mock->exps_alloc = N*a_mock->length;
165
166 b_mock->exps = A->exps + N*i;
167 b_mock->coeffs = A->coeffs + i;
168 b_mock->bits = bits;
169 while (i < A->length && (mask & (A->exps[N*i + off] >> shift)) == 1)
170 i++;
171 b_mock->length = i - a_mock->length;
172 b_mock->coeffs_alloc = b_mock->length;
173 b_mock->exps_alloc = N*b_mock->length;
174
175 c_mock->exps = A->exps + N*i;
176 c_mock->coeffs = A->coeffs + i;
177 c_mock->bits = bits;
178 c_mock->length = A->length - i;
179 c_mock->coeffs_alloc = c_mock->length;
180 c_mock->exps_alloc = N*c_mock->length;
181
182 FLINT_ASSERT(a_mock->length > 0);
183 FLINT_ASSERT(c_mock->length > 0);
184
185 nmod_mpoly_mul(t1, a_mock, c_mock, ctx);
186 nmod_mpoly_neg(t1, t1, ctx);
187 if (!nmod_mpoly_quadratic_root(t2, b_mock, t1, ctx))
188 {
189 nmod_mpolyv_fit_length(Af, 1, ctx);
190 Af->length = 1;
191 nmod_mpoly_swap(Af->coeffs + 0, A, ctx);
192 success = 1;
193 goto cleanup;
194 }
195 nmod_mpoly_neg(t2, t2, ctx);
196
197 success = nmod_mpoly_gcd_cofactors(t0, t1, t2, a_mock, t2, ctx);
198 if (!success)
199 goto cleanup;
200
201 nmod_mpoly_divides(t3, c_mock, t2, ctx);
202
203 nmod_mpolyv_fit_length(Af, 2, ctx);
204 Af->length = 2;
205 nmod_mpoly_add(Af->coeffs + 0, t1, t2, ctx);
206 nmod_mpoly_add(Af->coeffs + 1, t0, t3, ctx);
207
208 success = 1;
209
210 cleanup:
211
212 nmod_mpoly_clear(t0, ctx);
213 nmod_mpoly_clear(t1, ctx);
214 nmod_mpoly_clear(t2, ctx);
215 nmod_mpoly_clear(t3, ctx);
216
217 return success;
218 }
219
220 /*
221 The property "sep" used here is that of the returned factors of
222 _nmod_mpoly_factor_separable with sep = 1, namely:
223 (1) monic
224 (2) primitive wrt each variable
225 (3) for all i, derivative(A, gen(i)) = 0, or
226 gcd(A, derivative(A, gen(i))) = 1
227 (4) there is at least one i for which derivative(A, gen(i)) != 0
228
229 Input A is sep and compressed.
230
231 return 1 for success, 0 for failure
232 */
_factor_irred_compressed(nmod_mpolyv_t Af,nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx,unsigned int algo)233 static int _factor_irred_compressed(
234 nmod_mpolyv_t Af,
235 nmod_mpoly_t A,
236 const nmod_mpoly_ctx_t ctx,
237 unsigned int algo)
238 {
239 int success;
240 slong i, j, tot_deg;
241 slong nvars = ctx->minfo->nvars;
242 slong * perm;
243 ulong * strides, * texps;
244 flint_bitcnt_t Abits;
245 flint_rand_t state;
246 #if FLINT_WANT_ASSERT
247 nmod_mpoly_t Aorg;
248
249 nmod_mpoly_init(Aorg, ctx);
250 nmod_mpoly_set(Aorg, A, ctx);
251 #endif
252
253 FLINT_ASSERT(A->length > 0);
254 FLINT_ASSERT(A->coeffs[0] == 1);
255 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
256
257 if (A->length < 2)
258 {
259 FLINT_ASSERT(A->length == 1);
260 FLINT_ASSERT(!nmod_mpoly_is_ui(A, ctx));
261 nmod_mpolyv_fit_length(Af, 1, ctx);
262 nmod_mpoly_swap(Af->coeffs + 0, A, ctx);
263 Af->length = 1;
264 return 1;
265 }
266
267 if (A->bits > FLINT_BITS &&
268 !nmod_mpoly_repack_bits_inplace(A, FLINT_BITS, ctx))
269 {
270 return 0;
271 }
272
273 Abits = A->bits;
274
275 flint_randinit(state);
276
277 strides = FLINT_ARRAY_ALLOC(2*nvars, ulong);
278 texps = strides + nvars;
279
280 perm = FLINT_ARRAY_ALLOC(nvars, slong);
281
282 /* fill perm with id, and fill in strides */
283 {
284 ulong ppowt, ppow = ctx->mod.n;
285 slong N = mpoly_words_per_exp_sp(Abits, ctx->minfo);
286
287 while (!n_mul_checked(&ppowt, ppow, ctx->mod.n))
288 ppow = ppowt;
289
290 for (j = 0; j < nvars; j++)
291 {
292 strides[j] = ppow;
293 perm[j] = j;
294 }
295
296 tot_deg = 1;
297 for (i = 0; i < A->length; i++)
298 {
299 slong this_deg = 0;
300 mpoly_get_monomial_ui(texps, A->exps + N*i, Abits, ctx->minfo);
301 for (j = 0; j < nvars; j++)
302 {
303 if (z_add_checked(&this_deg, this_deg, texps[j]))
304 {
305 success = 0;
306 goto cleanup;
307 }
308 strides[j] = n_gcd(strides[j], texps[j]);
309 }
310
311 tot_deg = FLINT_MAX(tot_deg, this_deg);
312 }
313 }
314
315 /* find permutation with gcd(A, derivative(A, gen(perm[0]))) = 1 */
316 for (i = 0; i < nvars; i++)
317 {
318 if (strides[i] == 1)
319 {
320 SLONG_SWAP(perm[0], perm[i]);
321 break;
322 }
323 }
324
325 if (nvars < 2)
326 {
327 nmod_poly_t Au;
328 nmod_poly_factor_t Auf;
329
330 FLINT_ASSERT(nvars == 1);
331
332 nmod_poly_init_mod(Au, ctx->mod);
333 nmod_poly_factor_init(Auf);
334
335 FLINT_ASSERT(nmod_mpoly_is_nmod_poly(A, perm[0], ctx));
336 success = nmod_mpoly_get_nmod_poly(Au, A, perm[0], ctx);
337 FLINT_ASSERT(success);
338 nmod_poly_factor(Auf, Au);
339
340 nmod_mpolyv_fit_length(Af, Auf->num, ctx);
341 Af->length = Auf->num;
342 for (i = 0; i < Auf->num; i++)
343 {
344 FLINT_ASSERT(Auf->exp[i] == 1);
345 _nmod_mpoly_set_nmod_poly(Af->coeffs + i, Abits,
346 Auf->p[i].coeffs, Auf->p[i].length, perm[0], ctx);
347 }
348
349 nmod_poly_clear(Au);
350 nmod_poly_factor_clear(Auf);
351
352 success = 1;
353 }
354 else if (nvars == 2)
355 {
356 n_poly_t c;
357 n_bpoly_t Ab;
358 n_tpoly_t Abf;
359
360 n_poly_init(c);
361 n_bpoly_init(Ab);
362 n_tpoly_init(Abf);
363
364 nmod_mpoly_get_bpoly(Ab, A, perm[0], perm[1], ctx);
365 success = n_bpoly_mod_factor_smprime(c, Abf, Ab, 1, ctx->mod);
366 if (!success)
367 {
368 nmod_mpoly_get_bpoly(Ab, A, perm[0], perm[1], ctx);
369 n_bpoly_mod_factor_lgprime(c, Abf, Ab, ctx->mod);
370 }
371
372 FLINT_ASSERT(n_poly_degree(c) == 0);
373
374 nmod_mpolyv_fit_length(Af, Abf->length, ctx);
375 Af->length = Abf->length;
376 for (i = 0; i < Abf->length; i++)
377 {
378 nmod_mpoly_set_bpoly(Af->coeffs + i, Abits, Abf->coeffs + i,
379 perm[0], perm[1], ctx);
380 nmod_mpoly_make_monic(Af->coeffs + i, Af->coeffs + i, ctx);
381 }
382
383 n_poly_clear(c);
384 n_bpoly_clear(Ab);
385 n_tpoly_clear(Abf);
386
387 success = 1;
388 }
389 else
390 {
391 slong Adeg0;
392 nmod_mpoly_t lcA;
393 nmod_mpoly_factor_t lcAf;
394
395 nmod_mpoly_init(lcA, ctx);
396 nmod_mpoly_factor_init(lcAf, ctx);
397
398 tot_deg = _deflate(A, tot_deg, strides, perm, ctx);
399
400 #if FLINT_WANT_ASSERT
401 {
402 nmod_mpoly_t g;
403 nmod_mpoly_init(g, ctx);
404 nmod_mpoly_derivative(g, A, 0, ctx);
405 FLINT_ASSERT(nmod_mpoly_gcd(g, g, A, ctx));
406 FLINT_ASSERT(nmod_mpoly_is_one(g, ctx));
407 nmod_mpoly_clear(g, ctx);
408 }
409 #endif
410
411 Adeg0 = nmod_mpoly_degree_si(A, 0, ctx);
412 if (Adeg0 == 1)
413 {
414 nmod_mpolyv_fit_length(Af, 1, ctx);
415 Af->length = 1;
416 nmod_mpoly_swap(Af->coeffs + 0, A, ctx);
417 success = 1;
418 goto cleanup_inflate;
419 }
420 else if (Adeg0 == 2)
421 {
422 success = _apply_quadratic(Af, A, ctx);
423 goto cleanup_inflate;
424 }
425
426 success = 0;
427
428 if (!(algo & (MPOLY_FACTOR_USE_WANG | MPOLY_FACTOR_USE_ZIP)))
429 goto try_zassenhaus;
430
431 /* TODO lcc_kaltofen */
432 _nmod_mpoly_get_lead0(lcA, A, ctx);
433 if (!nmod_mpoly_factor(lcAf, lcA, ctx))
434 goto try_zassenhaus;
435
436 if (!(algo & MPOLY_FACTOR_USE_ZIP))
437 {
438 if (success == 0)
439 success = nmod_mpoly_factor_irred_smprime_wang(
440 Af, A, lcAf, lcA, ctx, state);
441 if (success == 0)
442 success = nmod_mpoly_factor_irred_medprime_wang(
443 Af, A, lcAf, lcA, ctx, state);
444 if (success == 0)
445 success = nmod_mpoly_factor_irred_lgprime_wang(
446 Af, A, lcAf, lcA, ctx, state);
447 }
448 else if (!(algo & MPOLY_FACTOR_USE_WANG))
449 {
450 if (success == 0)
451 success = nmod_mpoly_factor_irred_smprime_zippel(
452 Af, A, lcAf, lcA, ctx, state);
453 if (success == 0)
454 success = nmod_mpoly_factor_irred_medprime_zippel(
455 Af, A, lcAf, lcA, ctx, state);
456 if (success == 0)
457 success = nmod_mpoly_factor_irred_lgprime_zippel(
458 Af, A, lcAf, lcA, ctx, state);
459 }
460 else
461 {
462 double tdensity;
463 fmpz_t x;
464 fmpz_init(x);
465 fmpz_bin_uiui(x, (ulong)tot_deg + nvars, nvars);
466 tdensity = A->length/fmpz_get_d(x);
467 fmpz_clear(x);
468
469 if (tdensity > 0.005)
470 {
471 if (success == 0)
472 success = nmod_mpoly_factor_irred_smprime_wang(
473 Af, A, lcAf, lcA, ctx, state);
474 if (success == 0)
475 success = nmod_mpoly_factor_irred_medprime_wang(
476 Af, A, lcAf, lcA, ctx, state);
477 if (success == 0)
478 success = nmod_mpoly_factor_irred_smprime_zippel(
479 Af, A, lcAf, lcA, ctx, state);
480 if (success == 0)
481 success = nmod_mpoly_factor_irred_medprime_zippel(
482 Af, A, lcAf, lcA, ctx, state);
483 }
484 else
485 {
486 if (success == 0)
487 success = nmod_mpoly_factor_irred_smprime_zippel(
488 Af, A, lcAf, lcA, ctx, state);
489 if (success == 0)
490 success = nmod_mpoly_factor_irred_medprime_zippel(
491 Af, A, lcAf, lcA, ctx, state);
492 if (success == 0)
493 success = nmod_mpoly_factor_irred_smprime_wang(
494 Af, A, lcAf, lcA, ctx, state);
495 if (success == 0)
496 success = nmod_mpoly_factor_irred_medprime_wang(
497 Af, A, lcAf, lcA, ctx, state);
498 }
499
500 if (tdensity > 0.001)
501 {
502 if (success == 0)
503 success = nmod_mpoly_factor_irred_lgprime_wang(
504 Af, A, lcAf, lcA, ctx, state);
505 if (success == 0)
506 success = nmod_mpoly_factor_irred_lgprime_zippel(
507 Af, A, lcAf, lcA, ctx, state);
508 }
509 else
510 {
511 if (success == 0)
512 success = nmod_mpoly_factor_irred_lgprime_zippel(
513 Af, A, lcAf, lcA, ctx, state);
514 if (success == 0)
515 success = nmod_mpoly_factor_irred_lgprime_wang(
516 Af, A, lcAf, lcA, ctx, state);
517 }
518 }
519
520 try_zassenhaus:
521
522 if (algo & MPOLY_FACTOR_USE_ZAS)
523 {
524 if (success == 0)
525 success = nmod_mpoly_factor_irred_smprime_zassenhaus(
526 Af, A, ctx, state);
527 if (success == 0)
528 success = nmod_mpoly_factor_irred_medprime_zassenhaus(
529 Af, A, ctx, state);
530 if (success == 0)
531 success = nmod_mpoly_factor_irred_lgprime_zassenhaus(
532 Af, A, ctx, state);
533 }
534
535 cleanup_inflate:
536
537 success = (success > 0);
538 if (success)
539 {
540 for (i = 0; i < Af->length; i++)
541 _inflate(Af->coeffs + i, Abits, strides, perm, ctx);
542 }
543
544 nmod_mpoly_clear(lcA, ctx);
545 nmod_mpoly_factor_clear(lcAf, ctx);
546 }
547
548 cleanup:
549
550 flint_randclear(state);
551 flint_free(strides);
552 flint_free(perm);
553
554 #if FLINT_WANT_ASSERT
555 if (success)
556 {
557 nmod_mpoly_t prod;
558 nmod_mpoly_init(prod, ctx);
559 nmod_mpoly_one(prod, ctx);
560 for (i = 0; i < Af->length; i++)
561 nmod_mpoly_mul(prod, prod, Af->coeffs + i, ctx);
562
563 FLINT_ASSERT(nmod_mpoly_equal(prod, Aorg, ctx));
564 nmod_mpoly_clear(prod, ctx);
565 nmod_mpoly_clear(Aorg, ctx);
566 }
567 #endif
568
569 FLINT_ASSERT(success == 0 || success == 1);
570 return success;
571 }
572
573
574 /*
575 f is already squarefree
576 make the factors in f have the sep property
577 */
_refine_sep(nmod_mpolyv_t f,const nmod_mpoly_ctx_t ctx,nmod_mpolyv_t g)578 static int _refine_sep(
579 nmod_mpolyv_t f,
580 const nmod_mpoly_ctx_t ctx,
581 nmod_mpolyv_t g) /* temp */
582 {
583 int success;
584 slong v, i;
585 nmod_mpoly_struct * t;
586 nmod_mpoly_univar_t u;
587
588 nmod_mpoly_univar_init(u, ctx);
589
590 /* first make primitive */
591 for (v = 0; v < ctx->minfo->nvars; v++)
592 {
593 g->length = 0;
594 for (i = 0; i < f->length; i++)
595 {
596 nmod_mpoly_to_univar(u, f->coeffs + i, v, ctx);
597 FLINT_ASSERT(u->length > 0);
598 FLINT_ASSERT(fmpz_is_zero(u->exps + u->length - 1));
599
600 nmod_mpolyv_fit_length(g, g->length + 2, ctx);
601 success = _nmod_mpoly_vec_content_mpoly(g->coeffs + g->length,
602 u->coeffs, u->length, ctx);
603 if (!success)
604 goto cleanup;
605
606 if (nmod_mpoly_is_ui(g->coeffs + g->length, ctx))
607 {
608 nmod_mpoly_swap(g->coeffs + g->length, f->coeffs + i, ctx);
609 g->length++;
610 }
611 else
612 {
613 success = nmod_mpoly_divides(g->coeffs + g->length + 1,
614 f->coeffs + i, g->coeffs + g->length, ctx);
615 FLINT_ASSERT(success);
616
617 if (nmod_mpoly_is_ui(g->coeffs + g->length + 1, ctx))
618 g->length += 1;
619 else
620 g->length += 2;
621 }
622 }
623
624 nmod_mpolyv_swap(f, g, ctx);
625 }
626
627 /* now make separable/derivative zero wrt each variable */
628 nmod_mpolyv_fit_length(g, 1, ctx);
629 t = g->coeffs + 0;
630 for (v = 0; v < ctx->minfo->nvars; v++)
631 {
632 i = 0;
633 while (i < f->length)
634 {
635 nmod_mpoly_derivative(t, f->coeffs + i, v, ctx);
636 if (nmod_mpoly_is_zero(t, ctx))
637 {
638 /* f[i] has zero derivative */
639 i++;
640 continue;
641 }
642
643 nmod_mpolyv_fit_length(f, f->length + 1, ctx);
644
645 success = nmod_mpoly_gcd_cofactors(f->coeffs + f->length,
646 f->coeffs + i, t, f->coeffs + i, t, ctx);
647 if (!success)
648 goto cleanup;
649
650 if (nmod_mpoly_is_ui(f->coeffs + f->length, ctx))
651 {
652 /* f[i] is comprime with its derivative */
653 i++;
654 }
655 else
656 {
657 /* f[i] and f[end] at least got smaller */
658 f->length++;
659 }
660 }
661 }
662
663 success = 1;
664
665 cleanup:
666
667 nmod_mpoly_univar_clear(u, ctx);
668
669 return 1;
670 }
671
672
673 /*
674 A is sep.
675
676 return 1 for success, 0 for failure
677 */
_factor_irred(nmod_mpolyv_t Af,nmod_mpoly_t A,const nmod_mpoly_ctx_t Actx,unsigned int algo)678 static int _factor_irred(
679 nmod_mpolyv_t Af,
680 nmod_mpoly_t A,
681 const nmod_mpoly_ctx_t Actx,
682 unsigned int algo)
683 {
684 int success;
685 slong i, j;
686 flint_bitcnt_t Abits;
687 mpoly_compression_t M;
688 #if FLINT_WANT_ASSERT
689 nmod_mpoly_t Aorg;
690
691 nmod_mpoly_init(Aorg, Actx);
692 nmod_mpoly_set(Aorg, A, Actx);
693 #endif
694
695 FLINT_ASSERT(A->length > 0);
696 FLINT_ASSERT(A->coeffs[0] == 1);
697
698 if (A->length < 2)
699 {
700 FLINT_ASSERT(A->length == 1);
701 FLINT_ASSERT(!nmod_mpoly_is_ui(A, Actx));
702 nmod_mpolyv_fit_length(Af, 1, Actx);
703 Af->length = 1;
704 nmod_mpoly_swap(Af->coeffs + 0, A, Actx);
705 success = 1;
706 goto cleanup_less;
707 }
708
709 if (A->bits > FLINT_BITS &&
710 !nmod_mpoly_repack_bits_inplace(A, FLINT_BITS, Actx))
711 {
712 success = 0;
713 goto cleanup_less;
714 }
715
716 Abits = A->bits;
717
718 mpoly_compression_init(M);
719 mpoly_compression_set(M, A->exps, A->bits, A->length, Actx->minfo);
720
721 if (M->is_irred)
722 {
723 nmod_mpolyv_fit_length(Af, 1, Actx);
724 Af->length = 1;
725 nmod_mpoly_swap(Af->coeffs + 0, A, Actx);
726 success = 1;
727 }
728 else if (M->is_trivial)
729 {
730 success = _factor_irred_compressed(Af, A, Actx, algo);
731 }
732 else
733 {
734 nmod_mpoly_ctx_t Lctx;
735 nmod_mpolyv_t Lf, Lft, Lfs;
736
737 nmod_mpoly_ctx_init(Lctx, M->mvars, ORD_LEX, Actx->mod.n);
738 nmod_mpolyv_init(Lf, Lctx);
739 nmod_mpolyv_init(Lft, Lctx);
740 nmod_mpolyv_init(Lfs, Lctx);
741
742 nmod_mpolyv_fit_length(Lft, 1, Lctx);
743 Lft->length = 1;
744 nmod_mpoly_compression_do(Lft->coeffs + 0, Lctx, A->coeffs, A->length, M);
745
746 _refine_sep(Lft, Lctx, Lf);
747
748 if (Lft->length == 1)
749 {
750 success = _factor_irred_compressed(Lf, Lft->coeffs + 0, Lctx, algo);
751 }
752 else
753 {
754 success = 1;
755 Lf->length = 0;
756 for (i = 0; i < Lft->length; i++)
757 {
758 success = _factor_irred(Lfs, Lft->coeffs + i, Lctx, algo);
759 if (!success)
760 break;
761
762 nmod_mpolyv_fit_length(Lf, Lf->length + Lfs->length, Lctx);
763 for (j = 0; j < Lfs->length; j++)
764 {
765 nmod_mpoly_swap(Lf->coeffs + Lf->length, Lfs->coeffs + j, Lctx);
766 Lf->length++;
767 }
768 }
769 }
770
771 if (success)
772 {
773 nmod_mpolyv_fit_length(Af, Lf->length, Actx);
774 Af->length = Lf->length;
775 for (i = 0; i < Lf->length; i++)
776 {
777 nmod_mpoly_compression_undo(Af->coeffs + i, Abits, Actx,
778 Lf->coeffs + i, Lctx, M);
779 }
780 }
781
782 nmod_mpolyv_clear(Lf, Lctx);
783 nmod_mpolyv_clear(Lft, Lctx);
784 nmod_mpolyv_clear(Lfs, Lctx);
785 nmod_mpoly_ctx_clear(Lctx);
786 }
787
788 mpoly_compression_clear(M);
789
790 cleanup_less:
791
792 #if FLINT_WANT_ASSERT
793 if (success)
794 {
795 nmod_mpoly_t prod;
796 nmod_mpoly_init(prod, Actx);
797 nmod_mpoly_one(prod, Actx);
798 for (i = 0; i < Af->length; i++)
799 nmod_mpoly_mul(prod, prod, Af->coeffs + i, Actx);
800 FLINT_ASSERT(nmod_mpoly_equal(prod, Aorg, Actx));
801 nmod_mpoly_clear(prod, Actx);
802 nmod_mpoly_clear(Aorg, Actx);
803 }
804 #endif
805
806 FLINT_ASSERT(success == 0 || success == 1);
807 return success;
808 }
809
810
811 /*
812 Assume each factor in f is sep.
813 Replace f by an irreducible factorization.
814 */
nmod_mpoly_factor_irred(nmod_mpoly_factor_t f,const nmod_mpoly_ctx_t ctx,unsigned int algo)815 int nmod_mpoly_factor_irred(
816 nmod_mpoly_factor_t f,
817 const nmod_mpoly_ctx_t ctx,
818 unsigned int algo)
819 {
820 int success;
821 slong i, j;
822 nmod_mpolyv_t t;
823 nmod_mpoly_factor_t g;
824
825 nmod_mpolyv_init(t, ctx);
826 nmod_mpoly_factor_init(g, ctx);
827
828 g->constant = f->constant;
829 g->num = 0;
830 for (j = 0; j < f->num; j++)
831 {
832 success = _factor_irred(t, f->poly + j, ctx, algo);
833 if (!success)
834 goto cleanup;
835
836 nmod_mpoly_factor_fit_length(g, g->num + t->length, ctx);
837 for (i = 0; i < t->length; i++)
838 {
839 fmpz_set(g->exp + g->num, f->exp + j);
840 nmod_mpoly_swap(g->poly + g->num, t->coeffs + i, ctx);
841 g->num++;
842 }
843 }
844 nmod_mpoly_factor_swap(f, g, ctx);
845
846 success = 1;
847
848 cleanup:
849
850 nmod_mpolyv_clear(t, ctx);
851 nmod_mpoly_factor_clear(g, ctx);
852
853 return success;
854 }
855
856
857 /*
858 append factor(f)^e to g
859 assuming f is compressed and content free
860 */
_compressed_content_to_irred(nmod_mpoly_factor_t g,nmod_mpoly_t f,const fmpz_t e,const nmod_mpoly_ctx_t ctx,unsigned int algo)861 static int _compressed_content_to_irred(
862 nmod_mpoly_factor_t g,
863 nmod_mpoly_t f,
864 const fmpz_t e,
865 const nmod_mpoly_ctx_t ctx,
866 unsigned int algo)
867 {
868 int success;
869 slong j, k;
870 nmod_mpoly_factor_t h;
871 nmod_mpolyv_t v;
872
873 nmod_mpoly_factor_init(h, ctx);
874 nmod_mpolyv_init(v, ctx);
875
876 success = _nmod_mpoly_factor_separable(h, f, ctx, 1);
877 if (!success)
878 goto cleanup;
879
880 for (j = 0; j < h->num; j++)
881 {
882 success = h->num > 1 ? _factor_irred(v, h->poly + j, ctx, algo) :
883 _factor_irred_compressed(v, h->poly + j, ctx, algo);
884 if (!success)
885 goto cleanup;
886
887 nmod_mpoly_factor_fit_length(g, g->num + v->length, ctx);
888 for (k = 0; k < v->length; k++)
889 {
890 fmpz_mul(g->exp + g->num, h->exp + j, e);
891 nmod_mpoly_swap(g->poly + g->num, v->coeffs + k, ctx);
892 g->num++;
893 }
894 }
895
896 cleanup:
897
898 nmod_mpoly_factor_clear(h, ctx);
899 nmod_mpolyv_clear(v, ctx);
900
901 return success;
902 }
903
904
nmod_mpoly_factor_algo(nmod_mpoly_factor_t f,const nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx,unsigned int algo)905 int nmod_mpoly_factor_algo(
906 nmod_mpoly_factor_t f,
907 const nmod_mpoly_t A,
908 const nmod_mpoly_ctx_t ctx,
909 unsigned int algo)
910 {
911 int success;
912 slong i, j;
913 flint_bitcnt_t bits;
914 nmod_mpoly_factor_t g;
915 mpoly_compression_t M;
916
917 if (!nmod_mpoly_factor_content(f, A, ctx))
918 return 0;
919
920 nmod_mpoly_factor_init(g, ctx);
921 mpoly_compression_init(M);
922
923 /* write into g */
924 g->constant = f->constant;
925 g->num = 0;
926 for (i = 0; i < f->num; i++)
927 {
928 if (f->poly[i].length < 2)
929 {
930 nmod_mpoly_factor_fit_length(g, g->num + 1, ctx);
931 nmod_mpoly_swap(g->poly + g->num, f->poly + i, ctx);
932 fmpz_swap(g->exp + g->num, f->exp + i);
933 g->num++;
934 continue;
935 }
936
937 if (f->poly[i].bits > FLINT_BITS &&
938 !nmod_mpoly_repack_bits_inplace(f->poly + i, FLINT_BITS, ctx))
939 {
940 success = 0;
941 goto cleanup;
942 }
943
944 bits = f->poly[i].bits;
945
946 mpoly_compression_set(M, f->poly[i].exps, bits, f->poly[i].length, ctx->minfo);
947 if (M->is_irred)
948 {
949 nmod_mpoly_factor_fit_length(g, g->num + 1, ctx);
950 nmod_mpoly_swap(g->poly + g->num, f->poly + i, ctx);
951 fmpz_swap(g->exp + g->num, f->exp + i);
952 g->num++;
953 }
954 else if (M->is_trivial)
955 {
956 success = _compressed_content_to_irred(g, f->poly + i, f->exp + i, ctx, algo);
957 if (!success)
958 goto cleanup;
959 }
960 else
961 {
962 nmod_mpoly_ctx_t Lctx;
963 nmod_mpoly_t L;
964 nmod_mpoly_factor_t h;
965
966 /* compression may have messed up the content factorization */
967 nmod_mpoly_ctx_init(Lctx, M->mvars, ORD_LEX, ctx->mod.n);
968 nmod_mpoly_init(L, Lctx);
969 nmod_mpoly_factor_init(h, Lctx);
970
971 nmod_mpoly_compression_do(L, Lctx, f->poly[i].coeffs,
972 f->poly[i].length, M);
973 if (M->is_perm)
974 {
975 success = _compressed_content_to_irred(h, L, f->exp + i, Lctx, algo);
976 fmpz_one(f->exp + i);
977 }
978 else
979 {
980 success = nmod_mpoly_factor_separable(h, L, Lctx, 1) &&
981 nmod_mpoly_factor_irred(h, Lctx, algo);
982 }
983
984 if (success)
985 {
986 FLINT_ASSERT(h->constant == 1);
987 nmod_mpoly_factor_fit_length(g, g->num + h->num, ctx);
988 for (j = 0; j < h->num; j++)
989 {
990 fmpz_mul(g->exp + g->num, f->exp + i, h->exp + j);
991 nmod_mpoly_compression_undo(g->poly + g->num, bits, ctx,
992 h->poly + j, Lctx, M);
993 g->num++;
994 }
995 }
996
997 nmod_mpoly_factor_clear(h, Lctx);
998 nmod_mpoly_clear(L, Lctx);
999 nmod_mpoly_ctx_clear(Lctx);
1000
1001 if (!success)
1002 goto cleanup;
1003 }
1004 }
1005
1006 nmod_mpoly_factor_swap(f, g, ctx);
1007
1008 success = 1;
1009
1010 cleanup:
1011
1012 nmod_mpoly_factor_clear(g, ctx);
1013 mpoly_compression_clear(M);
1014
1015 FLINT_ASSERT(!success || nmod_mpoly_factor_matches(A, f, ctx));
1016 return success;
1017 }
1018
1019
nmod_mpoly_factor(nmod_mpoly_factor_t f,const nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx)1020 int nmod_mpoly_factor(
1021 nmod_mpoly_factor_t f,
1022 const nmod_mpoly_t A,
1023 const nmod_mpoly_ctx_t ctx)
1024 {
1025 return nmod_mpoly_factor_algo(f, A, ctx, MPOLY_FACTOR_USE_ALL);
1026 }
1027
1028
nmod_mpoly_factor_zassenhaus(nmod_mpoly_factor_t f,const nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx)1029 int nmod_mpoly_factor_zassenhaus(
1030 nmod_mpoly_factor_t f,
1031 const nmod_mpoly_t A,
1032 const nmod_mpoly_ctx_t ctx)
1033 {
1034 return nmod_mpoly_factor_algo(f, A, ctx, MPOLY_FACTOR_USE_ZAS);
1035 }
1036
1037
nmod_mpoly_factor_wang(nmod_mpoly_factor_t f,const nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx)1038 int nmod_mpoly_factor_wang(
1039 nmod_mpoly_factor_t f,
1040 const nmod_mpoly_t A,
1041 const nmod_mpoly_ctx_t ctx)
1042 {
1043 return nmod_mpoly_factor_algo(f, A, ctx, MPOLY_FACTOR_USE_WANG);
1044 }
1045
1046
nmod_mpoly_factor_zippel(nmod_mpoly_factor_t f,const nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx)1047 int nmod_mpoly_factor_zippel(
1048 nmod_mpoly_factor_t f,
1049 const nmod_mpoly_t A,
1050 const nmod_mpoly_ctx_t ctx)
1051 {
1052 return nmod_mpoly_factor_algo(f, A, ctx, MPOLY_FACTOR_USE_ZIP);
1053 }
1054
1055