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