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 "fmpz_mod_vec.h"
13 #include "fmpz_mod_mpoly.h"
14 #include "fmpz_mod_mpoly_factor.h"
15
16 /*
17 For each j, set out[j] to the evaluation of A at x_i = alpha[i] (i != j)
18 i.e. if nvars = 3
19 out[0] = A(x, alpha[1], alpha[2])
20 out[1] = A(alpha[0], x, alpha[2])
21 out[2] = A(alpha[0], alpha[1], x)
22
23 If ignore[j] is nonzero, then out[j] need not be calculated, probably
24 because we shouldn't calculate it in dense form.
25 */
fmpz_mod_mpoly_evals(slong * Atdeg,fmpz_mod_poly_struct * out,const int * ignore,const fmpz_mod_mpoly_t A,ulong * Amin_exp,ulong * Amax_exp,ulong * Astride,const fmpz * alphas,const fmpz_mod_mpoly_ctx_t ctx)26 static void fmpz_mod_mpoly_evals(
27 slong * Atdeg, /* total degree of deflated A, or -1 for overflow */
28 fmpz_mod_poly_struct * out,
29 const int * ignore,
30 const fmpz_mod_mpoly_t A,
31 ulong * Amin_exp,
32 ulong * Amax_exp,
33 ulong * Astride,
34 const fmpz * alphas,
35 const fmpz_mod_mpoly_ctx_t ctx)
36 {
37 slong i, j;
38 slong nvars = ctx->minfo->nvars;
39 ulong mask = (-UWORD(1)) >> (FLINT_BITS - A->bits);
40 slong * offsets, * shifts;
41 slong N = mpoly_words_per_exp_sp(A->bits, ctx->minfo);
42 ulong * varexps;
43 ulong varexp;
44 slong total_degree, lo, hi;
45 fmpz_t meval, t, t1;
46
47 fmpz_init(meval);
48 fmpz_init(t);
49 fmpz_init(t1);
50
51 offsets = FLINT_ARRAY_ALLOC(2*nvars, slong);
52 shifts = offsets + nvars;
53 varexps = FLINT_ARRAY_ALLOC(nvars, ulong);
54 for (j = 0; j < nvars; j++)
55 {
56 fmpz_mod_poly_zero(out + j, ctx->ffinfo);
57 mpoly_gen_offset_shift_sp(offsets + j, shifts + j, j, A->bits, ctx->minfo);
58 }
59
60 total_degree = 0;
61 for (i = 0; i < A->length; i++)
62 {
63 fmpz * s = A->coeffs + i;
64
65 hi = lo = 0;
66 for (j = 0; j < nvars; j++)
67 {
68 varexp = ((A->exps + N*i)[offsets[j]]>>shifts[j])&mask;
69
70 FLINT_ASSERT((Astride[j] == 0 && varexp == Amin_exp[j]) ||
71 (varexp - Amin_exp[j]) % Astride[j] == 0);
72
73 varexps[j] = Astride[j] < 2 ? varexp - Amin_exp[j] :
74 (varexp - Amin_exp[j])/Astride[j];
75
76 add_ssaaaa(hi, lo, hi, lo, 0, varexps[j]);
77
78 fmpz_mod_pow_ui(t1, alphas + j, varexps[j], ctx->ffinfo);
79 fmpz_mod_mul(meval, s, t1, ctx->ffinfo);
80 s = meval;
81 }
82
83 if (hi == 0 && FLINT_SIGN_EXT(lo) == 0 && total_degree >= 0)
84 total_degree = FLINT_MAX(total_degree, lo);
85 else
86 total_degree = -1;
87
88 for (j = 0; j < nvars; j++)
89 {
90 varexp = varexps[j];
91
92 if (ignore[j])
93 continue;
94
95 fmpz_mod_poly_fit_length(out + j, varexp + 1, ctx->ffinfo);
96
97 while (out[j].length <= varexp)
98 {
99 fmpz_zero(out[j].coeffs + out[j].length);
100 out[j].length++;
101 }
102
103 fmpz_mod_inv(t1, alphas + j, ctx->ffinfo);
104 fmpz_mod_pow_ui(t1, t1, varexps[j], ctx->ffinfo);
105 fmpz_mod_mul(t, meval, t1, ctx->ffinfo);
106 fmpz_mod_add(out[j].coeffs + varexp, out[j].coeffs + varexp, t,
107 ctx->ffinfo);
108 }
109 }
110
111 *Atdeg = total_degree;
112
113 for (j = 0; j < nvars; j++)
114 _fmpz_mod_poly_normalise(out + j);
115
116 flint_free(offsets);
117 flint_free(varexps);
118
119 fmpz_clear(meval);
120 fmpz_clear(t);
121 fmpz_clear(t1);
122 }
123
124
_set_estimates(mpoly_gcd_info_t I,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx)125 static void _set_estimates(
126 mpoly_gcd_info_t I,
127 const fmpz_mod_mpoly_t A,
128 const fmpz_mod_mpoly_t B,
129 const fmpz_mod_mpoly_ctx_t ctx)
130 {
131 int try_count = 0;
132 slong nvars = ctx->minfo->nvars;
133 slong i, j;
134 fmpz_mod_poly_t Geval;
135 fmpz_mod_poly_struct * Aevals, * Bevals;
136 fmpz * alphas;
137 flint_rand_t state;
138 slong ignore_limit;
139 int * ignore;
140
141 flint_randinit(state);
142
143 ignore = FLINT_ARRAY_ALLOC(nvars, int);
144 alphas = _fmpz_vec_init(nvars);
145 Aevals = FLINT_ARRAY_ALLOC(nvars, fmpz_mod_poly_struct);
146 Bevals = FLINT_ARRAY_ALLOC(nvars, fmpz_mod_poly_struct);
147
148 fmpz_mod_poly_init(Geval, ctx->ffinfo);
149 for (j = 0; j < nvars; j++)
150 {
151 fmpz_mod_poly_init(Aevals + j, ctx->ffinfo);
152 fmpz_mod_poly_init(Bevals + j, ctx->ffinfo);
153 }
154
155 ignore_limit = (A->length + B->length)/4096;
156 ignore_limit = FLINT_MAX(WORD(9999), ignore_limit);
157 I->Gdeflate_deg_bounds_are_nice = 1;
158 for (j = 0; j < nvars; j++)
159 {
160 if (I->Adeflate_deg[j] > ignore_limit ||
161 I->Bdeflate_deg[j] > ignore_limit)
162 {
163 ignore[j] = 1;
164 I->Gdeflate_deg_bounds_are_nice = 0;
165 }
166 else
167 {
168 ignore[j] = 0;
169 }
170 }
171
172 try_again:
173
174 if (++try_count > 10)
175 {
176 I->Gdeflate_deg_bounds_are_nice = 0;
177 for (j = 0; j < nvars; j++)
178 {
179 I->Gdeflate_deg_bound[j] = FLINT_MIN(I->Adeflate_deg[j],
180 I->Bdeflate_deg[j]);
181 I->Gterm_count_est[j] = 1 + I->Gdeflate_deg_bound[j]/2;
182 }
183
184 goto cleanup;
185 }
186
187 for (j = 0; j < nvars; j++)
188 fmpz_mod_rand_not_zero(alphas + j, state, ctx->ffinfo);
189
190 fmpz_mod_mpoly_evals(&I->Adeflate_tdeg, Aevals, ignore, A,
191 I->Amin_exp, I->Amax_exp, I->Gstride, alphas, ctx);
192 fmpz_mod_mpoly_evals(&I->Bdeflate_tdeg, Bevals, ignore, B,
193 I->Bmin_exp, I->Bmax_exp, I->Gstride, alphas, ctx);
194
195 for (j = 0; j < nvars; j++)
196 {
197 if (ignore[j])
198 {
199 I->Gdeflate_deg_bound[j] = FLINT_MIN(I->Adeflate_deg[j],
200 I->Bdeflate_deg[j]);
201 I->Gterm_count_est[j] = 1 + I->Gdeflate_deg_bound[j]/2;
202 }
203 else
204 {
205 if (I->Adeflate_deg[j] != fmpz_mod_poly_degree(Aevals + j, ctx->ffinfo) ||
206 I->Bdeflate_deg[j] != fmpz_mod_poly_degree(Bevals + j, ctx->ffinfo))
207 {
208 goto try_again;
209 }
210
211 fmpz_mod_poly_gcd(Geval, Aevals + j, Bevals + j, ctx->ffinfo);
212
213 I->Gterm_count_est[j] = 0;
214 I->Gdeflate_deg_bound[j] = fmpz_mod_poly_degree(Geval, ctx->ffinfo);
215 for (i = I->Gdeflate_deg_bound[j]; i >= 0; i--)
216 I->Gterm_count_est[j] += (Geval->coeffs[i] != 0);
217 }
218 }
219
220 cleanup:
221
222 fmpz_mod_poly_clear(Geval, ctx->ffinfo);
223 for (j = 0; j < nvars; j++)
224 {
225 fmpz_mod_poly_clear(Aevals + j, ctx->ffinfo);
226 fmpz_mod_poly_clear(Bevals + j, ctx->ffinfo);
227 }
228
229 flint_free(ignore);
230 _fmpz_vec_clear(alphas, nvars);
231 flint_free(Aevals);
232 flint_free(Bevals);
233
234 flint_randclear(state);
235
236 return;
237 }
238
239
240 /* (Abar, Bbar) = (A, B) */
_parallel_set(fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx)241 static void _parallel_set(
242 fmpz_mod_mpoly_t Abar, /* could be NULL */
243 fmpz_mod_mpoly_t Bbar, /* could be NULL */
244 const fmpz_mod_mpoly_t A,
245 const fmpz_mod_mpoly_t B,
246 const fmpz_mod_mpoly_ctx_t ctx)
247 {
248 if (Abar == B && Bbar == A)
249 {
250 FLINT_ASSERT(Abar != NULL && Bbar != NULL);
251 fmpz_mod_mpoly_set(Abar, B, ctx);
252 fmpz_mod_mpoly_set(Bbar, A, ctx);
253 fmpz_mod_mpoly_swap(Abar, Bbar, ctx);
254 }
255 else if (Abar == B && Bbar != A)
256 {
257 FLINT_ASSERT(Abar != NULL);
258 if (Bbar != NULL)
259 fmpz_mod_mpoly_set(Bbar, B, ctx);
260 fmpz_mod_mpoly_set(Abar, A, ctx);
261 }
262 else
263 {
264 if (Abar != NULL)
265 fmpz_mod_mpoly_set(Abar, A, ctx);
266 if (Bbar != NULL)
267 fmpz_mod_mpoly_set(Bbar, B, ctx);
268 }
269 }
270
271
272 /* The variables in ess(A) and ess(B) are disjoint. gcd is trivial to compute */
_do_trivial(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)273 static void _do_trivial(
274 fmpz_mod_mpoly_t G,
275 fmpz_mod_mpoly_t Abar, /* could be NULL */
276 fmpz_mod_mpoly_t Bbar, /* could be NULL */
277 const fmpz_mod_mpoly_t A,
278 const fmpz_mod_mpoly_t B,
279 const mpoly_gcd_info_t I,
280 const fmpz_mod_mpoly_ctx_t ctx)
281 {
282 _parallel_set(Abar, Bbar, A, B, ctx);
283
284 if (Abar != NULL)
285 mpoly_monomials_shift_right_ui(Abar->exps, Abar->bits, Abar->length,
286 I->Gmin_exp, ctx->minfo);
287
288 if (Bbar != NULL)
289 mpoly_monomials_shift_right_ui(Bbar->exps, Bbar->bits, Bbar->length,
290 I->Gmin_exp, ctx->minfo);
291
292 fmpz_mod_mpoly_fit_length_reset_bits(G, 1, I->Gbits, ctx);
293 mpoly_set_monomial_ui(G->exps, I->Gmin_exp, I->Gbits, ctx->minfo);
294 fmpz_one(G->coeffs + 0);
295 _fmpz_mod_mpoly_set_length(G, 1, ctx);
296 }
297
298 /*********************** Easy when B is a monomial ***************************/
_do_monomial_gcd(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx)299 static int _do_monomial_gcd(
300 fmpz_mod_mpoly_t G,
301 fmpz_mod_mpoly_t Abar, /* could be NULL */
302 fmpz_mod_mpoly_t Bbar, /* could be NULL */
303 const fmpz_mod_mpoly_t A,
304 const fmpz_mod_mpoly_t B,
305 const fmpz_mod_mpoly_ctx_t ctx)
306 {
307 slong i;
308 flint_bitcnt_t Gbits = FLINT_MIN(A->bits, B->bits);
309 fmpz * minAfields, * minAdegs, * minBdegs;
310 TMP_INIT;
311
312 FLINT_ASSERT(A->length > 0);
313 FLINT_ASSERT(B->length == 1);
314
315 TMP_START;
316
317 /* get the field-wise minimum of A */
318 minAfields = (fmpz *) TMP_ALLOC(ctx->minfo->nfields*sizeof(fmpz));
319 for (i = 0; i < ctx->minfo->nfields; i++)
320 fmpz_init(minAfields + i);
321 mpoly_min_fields_fmpz(minAfields, A->exps, A->length, A->bits, ctx->minfo);
322
323 /* unpack to get the min degrees of each variable in A */
324 minAdegs = (fmpz *) TMP_ALLOC(ctx->minfo->nvars*sizeof(fmpz));
325 for (i = 0; i < ctx->minfo->nvars; i++)
326 fmpz_init(minAdegs + i);
327 mpoly_get_monomial_ffmpz_unpacked_ffmpz(minAdegs, minAfields, ctx->minfo);
328
329 /* get the degree of each variable in B */
330 minBdegs = (fmpz *) TMP_ALLOC(ctx->minfo->nvars*sizeof(fmpz));
331 for (i = 0; i < ctx->minfo->nvars; i++)
332 fmpz_init(minBdegs + i);
333 mpoly_get_monomial_ffmpz(minBdegs, B->exps, B->bits, ctx->minfo);
334
335 /* compute the degree of each variable in G */
336 _fmpz_vec_min_inplace(minBdegs, minAdegs, ctx->minfo->nvars);
337
338 _parallel_set(Abar, Bbar, A, B, ctx);
339
340 if (Abar != NULL)
341 mpoly_monomials_shift_right_ffmpz(Abar->exps, Abar->bits, Abar->length,
342 minBdegs, ctx->minfo);
343
344 if (Bbar != NULL)
345 mpoly_monomials_shift_right_ffmpz(Bbar->exps, Bbar->bits, Bbar->length,
346 minBdegs, ctx->minfo);
347
348 fmpz_mod_mpoly_fit_length_reset_bits(G, 1, Gbits, ctx);
349 mpoly_set_monomial_ffmpz(G->exps, minBdegs, Gbits, ctx->minfo);
350 fmpz_one(G->coeffs + 0);
351 _fmpz_mod_mpoly_set_length(G, 1, ctx);
352
353 for (i = 0; i < ctx->minfo->nfields; i++)
354 {
355 fmpz_clear(minAfields + i);
356 }
357 for (i = 0; i < ctx->minfo->nvars; i++)
358 {
359 fmpz_clear(minAdegs + i);
360 fmpz_clear(minBdegs + i);
361 }
362
363 TMP_END;
364
365 return 1;
366 }
367
368
369 /********************** See if cofactors are monomials ***********************/
_try_monomial_cofactors(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx)370 static int _try_monomial_cofactors(
371 fmpz_mod_mpoly_t G,
372 fmpz_mod_mpoly_t Abar, /* could be NULL */
373 fmpz_mod_mpoly_t Bbar, /* could be NULL */
374 const fmpz_mod_mpoly_t A,
375 const fmpz_mod_mpoly_t B,
376 const fmpz_mod_mpoly_ctx_t ctx)
377 {
378 int success;
379 slong i, j;
380 slong NA, NG;
381 slong nvars = ctx->minfo->nvars;
382 fmpz * Abarexps, * Bbarexps, * Texps;
383 fmpz_t a0, b0, t1, t2;
384 fmpz_mod_mpoly_t T;
385 flint_bitcnt_t Gbits = FLINT_MIN(A->bits, B->bits);
386 flint_bitcnt_t Abarbits = A->bits;
387 flint_bitcnt_t Bbarbits = B->bits;
388 TMP_INIT;
389
390 FLINT_ASSERT(A->length > 0);
391 FLINT_ASSERT(B->length > 0);
392
393 if (A->length != B->length)
394 return 0;
395
396 fmpz_init(t1);
397 fmpz_init(t2);
398 fmpz_init(a0);
399 fmpz_init(b0);
400
401 for (i = A->length - 1; i > 0; i--)
402 {
403 fmpz_mod_mul(t1, A->coeffs + 0, B->coeffs + i, ctx->ffinfo);
404 fmpz_mod_mul(t2, B->coeffs + 0, A->coeffs + i, ctx->ffinfo);
405 success = fmpz_equal(t1, t2);
406 if (!success)
407 goto cleanup_less;
408 }
409
410 fmpz_set(a0, A->coeffs + 0);
411 fmpz_set(b0, B->coeffs + 0);
412
413 TMP_START;
414
415 Abarexps = (fmpz *) TMP_ALLOC(3*nvars*sizeof(fmpz));
416 Bbarexps = Abarexps + 1*nvars;
417 Texps = Abarexps + 2*nvars;
418 for (j = 0; j < nvars; j++)
419 {
420 fmpz_init(Abarexps + j);
421 fmpz_init(Bbarexps + j);
422 fmpz_init(Texps + j);
423 }
424
425 success = mpoly_monomial_cofactors(Abarexps, Bbarexps, A->exps, A->bits,
426 B->exps, B->bits, A->length, ctx->minfo);
427 if (!success)
428 goto cleanup_more;
429
430 fmpz_mod_mpoly_init3(T, A->length, Gbits, ctx);
431 NG = mpoly_words_per_exp(Gbits, ctx->minfo);
432 NA = mpoly_words_per_exp(A->bits, ctx->minfo);
433 fmpz_mod_inv(t1, A->coeffs + 0, ctx->ffinfo);
434 T->length = A->length;
435 for (i = 0; i < A->length; i++)
436 {
437 mpoly_get_monomial_ffmpz(Texps, A->exps + NA*i, A->bits, ctx->minfo);
438 _fmpz_vec_sub(Texps, Texps, Abarexps, nvars);
439 mpoly_set_monomial_ffmpz(T->exps + NG*i, Texps, Gbits, ctx->minfo);
440 fmpz_mod_mul(T->coeffs + i, A->coeffs + i, t1, ctx->ffinfo);
441 }
442 fmpz_mod_mpoly_swap(G, T, ctx);
443 fmpz_mod_mpoly_clear(T, ctx);
444
445 if (Abar != NULL)
446 {
447 fmpz_mod_mpoly_fit_length_reset_bits(Abar, 1, Abarbits, ctx);
448 mpoly_set_monomial_ffmpz(Abar->exps, Abarexps, Abarbits, ctx->minfo);
449 fmpz_set(Abar->coeffs + 0, a0);
450 _fmpz_mod_mpoly_set_length(Abar, 1, ctx);
451 }
452
453 if (Bbar != NULL)
454 {
455 fmpz_mod_mpoly_fit_length_reset_bits(Bbar, 1, Bbarbits, ctx);
456 mpoly_set_monomial_ffmpz(Bbar->exps, Bbarexps, Bbarbits, ctx->minfo);
457 fmpz_set(Bbar->coeffs + 0, b0);
458 _fmpz_mod_mpoly_set_length(Bbar, 1, ctx);
459 }
460
461 success = 1;
462
463 cleanup_more:
464
465 for (j = 0; j < nvars; j++)
466 {
467 fmpz_clear(Abarexps + j);
468 fmpz_clear(Bbarexps + j);
469 fmpz_clear(Texps + j);
470 }
471
472 TMP_END;
473
474 cleanup_less:
475
476 fmpz_clear(t1);
477 fmpz_clear(t2);
478 fmpz_clear(a0);
479 fmpz_clear(b0);
480
481 return success;
482 }
483
484
485 /*** ess(A) and ess(B) depend on only one variable v_in_both ****************/
_do_univar(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,slong v_in_both,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)486 static int _do_univar(
487 fmpz_mod_mpoly_t G,
488 fmpz_mod_mpoly_t Abar,
489 fmpz_mod_mpoly_t Bbar,
490 const fmpz_mod_mpoly_t A,
491 const fmpz_mod_mpoly_t B,
492 slong v_in_both,
493 const mpoly_gcd_info_t I,
494 const fmpz_mod_mpoly_ctx_t ctx)
495 {
496 fmpz_mod_poly_t a, b, g, t, r;
497
498 fmpz_mod_poly_init(a, ctx->ffinfo);
499 fmpz_mod_poly_init(b, ctx->ffinfo);
500 fmpz_mod_poly_init(g, ctx->ffinfo);
501 fmpz_mod_poly_init(t, ctx->ffinfo);
502 fmpz_mod_poly_init(r, ctx->ffinfo);
503
504 _fmpz_mod_mpoly_to_fmpz_mod_poly_deflate(a, A, v_in_both, I->Amin_exp, I->Gstride, ctx);
505 _fmpz_mod_mpoly_to_fmpz_mod_poly_deflate(b, B, v_in_both, I->Bmin_exp, I->Gstride, ctx);
506 fmpz_mod_poly_gcd(g, a, b, ctx->ffinfo);
507 _fmpz_mod_mpoly_from_fmpz_mod_poly_inflate(G, I->Gbits, g, v_in_both,
508 I->Gmin_exp, I->Gstride, ctx);
509 if (Abar != NULL)
510 {
511 fmpz_mod_poly_divrem(t, r, a, g, ctx->ffinfo);
512 _fmpz_mod_mpoly_from_fmpz_mod_poly_inflate(Abar, I->Abarbits, t, v_in_both,
513 I->Abarmin_exp, I->Gstride, ctx);
514 }
515
516 if (Bbar != NULL)
517 {
518 fmpz_mod_poly_divrem(t, r, b, g, ctx->ffinfo);
519 _fmpz_mod_mpoly_from_fmpz_mod_poly_inflate(Bbar, I->Bbarbits, t, v_in_both,
520 I->Bbarmin_exp, I->Gstride, ctx);
521 }
522
523 fmpz_mod_poly_clear(a, ctx->ffinfo);
524 fmpz_mod_poly_clear(b, ctx->ffinfo);
525 fmpz_mod_poly_clear(g, ctx->ffinfo);
526 fmpz_mod_poly_clear(t, ctx->ffinfo);
527 fmpz_mod_poly_clear(r, ctx->ffinfo);
528
529 return 1;
530 }
531
532
533 /********* Assume B has length one when converted to univar format ***********/
_try_missing_var(fmpz_mod_mpoly_t G,flint_bitcnt_t Gbits,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,slong var,const fmpz_mod_mpoly_t A,ulong Ashift,const fmpz_mod_mpoly_t B,ulong Bshift,const fmpz_mod_mpoly_ctx_t ctx)534 static int _try_missing_var(
535 fmpz_mod_mpoly_t G, flint_bitcnt_t Gbits,
536 fmpz_mod_mpoly_t Abar,
537 fmpz_mod_mpoly_t Bbar,
538 slong var,
539 const fmpz_mod_mpoly_t A, ulong Ashift,
540 const fmpz_mod_mpoly_t B, ulong Bshift,
541 const fmpz_mod_mpoly_ctx_t ctx)
542 {
543 int success;
544 fmpz_mod_mpoly_univar_t Au;
545
546 fmpz_mod_mpoly_univar_init(Au, ctx);
547
548 #if FLINT_WANT_ASSERT
549 fmpz_mod_mpoly_to_univar(Au, B, var, ctx);
550 FLINT_ASSERT(Au->length == 1);
551 #endif
552
553 fmpz_mod_mpoly_to_univar(Au, A, var, ctx);
554
555 fmpz_mod_mpoly_univar_fit_length(Au, Au->length + 1, ctx);
556 fmpz_mod_mpoly_set(Au->coeffs + Au->length, B, ctx);
557 Au->length++;
558
559 if (Abar == NULL && Bbar == NULL)
560 {
561 success = _fmpz_mod_mpoly_vec_content_mpoly(G, Au->coeffs, Au->length, ctx);
562 if (!success)
563 goto cleanup;
564
565 fmpz_mod_mpoly_repack_bits_inplace(G, Gbits, ctx);
566 _mpoly_gen_shift_left(G->exps, G->bits, G->length,
567 var, FLINT_MIN(Ashift, Bshift), ctx->minfo);
568 }
569 else
570 {
571 fmpz_mod_mpoly_t tG, tAbar, tBbar;
572
573 fmpz_mod_mpoly_init(tG, ctx);
574 fmpz_mod_mpoly_init(tAbar, ctx);
575 fmpz_mod_mpoly_init(tBbar, ctx);
576
577 success = _fmpz_mod_mpoly_vec_content_mpoly(tG, Au->coeffs, Au->length, ctx);
578 if (!success)
579 goto cleanup;
580
581 fmpz_mod_mpoly_repack_bits_inplace(tG, Gbits, ctx);
582 _mpoly_gen_shift_left(tG->exps, tG->bits, tG->length,
583 var, FLINT_MIN(Ashift, Bshift), ctx->minfo);
584
585 if (Abar != NULL)
586 {
587 success = fmpz_mod_mpoly_divides(tAbar, A, tG, ctx);
588 FLINT_ASSERT(success);
589 }
590
591 if (Bbar != NULL)
592 {
593 success = fmpz_mod_mpoly_divides(tBbar, B, tG, ctx);
594 FLINT_ASSERT(success);
595 }
596
597 fmpz_mod_mpoly_swap(G, tG, ctx);
598
599 if (Abar != NULL)
600 fmpz_mod_mpoly_swap(Abar, tAbar, ctx);
601
602 if (Bbar != NULL)
603 fmpz_mod_mpoly_swap(Bbar, tBbar, ctx);
604
605 fmpz_mod_mpoly_clear(tG, ctx);
606 fmpz_mod_mpoly_clear(tAbar, ctx);
607 fmpz_mod_mpoly_clear(tBbar, ctx);
608 }
609
610 success = 1;
611
612 cleanup:
613
614 fmpz_mod_mpoly_univar_clear(Au, ctx);
615
616 return success;
617 }
618
619
620 /************************ See if B divides A ********************************/
_try_divides(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t BB,const fmpz_mod_mpoly_ctx_t ctx)621 static int _try_divides(
622 fmpz_mod_mpoly_t G,
623 fmpz_mod_mpoly_t Abar,
624 fmpz_mod_mpoly_t Bbar,
625 const fmpz_mod_mpoly_t A,
626 const fmpz_mod_mpoly_t BB,
627 const fmpz_mod_mpoly_ctx_t ctx)
628 {
629 int success = 0;
630 fmpz_mod_mpoly_t Q, B, M;
631
632 fmpz_mod_mpoly_init(Q, ctx);
633 fmpz_mod_mpoly_init(B, ctx);
634 fmpz_mod_mpoly_init(M, ctx);
635
636 /* BB = M*B */
637 fmpz_mod_mpoly_term_content(M, BB, ctx);
638 fmpz_mod_mpoly_divides(B, BB, M, ctx);
639
640 success = fmpz_mod_mpoly_divides(Q, A, B, ctx);
641 if (success)
642 {
643 _do_monomial_gcd(G, Abar, Bbar, Q, M, ctx);
644 fmpz_mod_mpoly_mul(G, G, B, ctx);
645 }
646
647 fmpz_mod_mpoly_clear(Q, ctx);
648 fmpz_mod_mpoly_clear(B, ctx);
649 fmpz_mod_mpoly_clear(M, ctx);
650
651 return success;
652 }
653
654
655 /********************* Hit A and B with a prs ********************************/
656
_try_prs(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)657 static int _try_prs(
658 fmpz_mod_mpoly_t G,
659 fmpz_mod_mpoly_t Abar,
660 fmpz_mod_mpoly_t Bbar,
661 const fmpz_mod_mpoly_t A,
662 const fmpz_mod_mpoly_t B,
663 const mpoly_gcd_info_t I,
664 const fmpz_mod_mpoly_ctx_t ctx)
665 {
666 int success;
667 slong j, var = -WORD(1);
668 fmpz_mod_mpoly_t Ac, Bc, Gc, s, t;
669 fmpz_mod_mpoly_univar_t Ax, Bx, Gx;
670
671 for (j = 0; j < ctx->minfo->nvars; j++)
672 {
673 if (I->Amax_exp[j] <= I->Amin_exp[j] ||
674 I->Bmax_exp[j] <= I->Bmin_exp[j] /* not needed */)
675 {
676 FLINT_ASSERT(I->Amax_exp[j] == I->Amin_exp[j]);
677 FLINT_ASSERT(I->Bmax_exp[j] == I->Bmin_exp[j]);
678 continue;
679 }
680
681 FLINT_ASSERT(I->Gstride[j] != UWORD(0));
682 FLINT_ASSERT((I->Amax_exp[j] - I->Amin_exp[j]) % I->Gstride[j] == 0);
683 FLINT_ASSERT((I->Bmax_exp[j] - I->Bmin_exp[j]) % I->Gstride[j] == 0);
684
685 var = j;
686 break;
687 }
688
689 if (var < 0)
690 return 0;
691
692 fmpz_mod_mpoly_init(Ac, ctx);
693 fmpz_mod_mpoly_init(Bc, ctx);
694 fmpz_mod_mpoly_init(Gc, ctx);
695 fmpz_mod_mpoly_init(s, ctx);
696 fmpz_mod_mpoly_init(t, ctx);
697 fmpz_mod_mpoly_univar_init(Ax, ctx);
698 fmpz_mod_mpoly_univar_init(Bx, ctx);
699 fmpz_mod_mpoly_univar_init(Gx, ctx);
700
701 fmpz_mod_mpoly_to_univar(Ax, A, var, ctx);
702 fmpz_mod_mpoly_to_univar(Bx, B, var, ctx);
703
704 success = _fmpz_mod_mpoly_vec_content_mpoly(Ac, Ax->coeffs, Ax->length, ctx) &&
705 _fmpz_mod_mpoly_vec_content_mpoly(Bc, Bx->coeffs, Bx->length, ctx) &&
706 fmpz_mod_mpoly_gcd(Gc, Ac, Bc, ctx);
707 if (!success)
708 goto cleanup;
709
710 _fmpz_mod_mpoly_vec_divexact_mpoly(Ax->coeffs, Ax->length, Ac, ctx);
711 _fmpz_mod_mpoly_vec_divexact_mpoly(Bx->coeffs, Bx->length, Bc, ctx);
712
713 success = fmpz_cmp(Ax->exps + 0, Bx->exps + 0) > 0 ?
714 fmpz_mod_mpoly_univar_pseudo_gcd(Gx, Ax, Bx, ctx) :
715 fmpz_mod_mpoly_univar_pseudo_gcd(Gx, Bx, Ax, ctx);
716 if (!success)
717 goto cleanup;
718
719 if (fmpz_mod_mpoly_gcd(t, Ax->coeffs + 0, Bx->coeffs + 0, ctx) &&
720 t->length == 1)
721 {
722 fmpz_mod_mpoly_term_content(s, Gx->coeffs + 0, ctx);
723 fmpz_mod_mpoly_divexact(t, Gx->coeffs + 0, s, ctx);
724 _fmpz_mod_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
725 }
726 else if (fmpz_mod_mpoly_gcd(t, Ax->coeffs + Ax->length - 1,
727 Bx->coeffs + Bx->length - 1, ctx) && t->length == 1)
728 {
729 fmpz_mod_mpoly_term_content(s, Gx->coeffs + Gx->length - 1, ctx);
730 fmpz_mod_mpoly_divexact(t, Gx->coeffs + Gx->length - 1, s, ctx);
731 _fmpz_mod_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
732 }
733
734 success = _fmpz_mod_mpoly_vec_content_mpoly(t, Gx->coeffs, Gx->length, ctx);
735 if (!success)
736 goto cleanup;
737
738 _fmpz_mod_mpoly_vec_divexact_mpoly(Gx->coeffs, Gx->length, t, ctx);
739 _fmpz_mod_mpoly_vec_mul_mpoly(Gx->coeffs, Gx->length, Gc, ctx);
740 _fmpz_mod_mpoly_from_univar(Gc, I->Gbits, Gx, var, ctx);
741
742 if (Abar != NULL)
743 fmpz_mod_mpoly_divexact(s, A, Gc, ctx);
744
745 if (Bbar != NULL)
746 fmpz_mod_mpoly_divexact(t, B, Gc, ctx);
747
748 fmpz_mod_mpoly_swap(G, Gc, ctx);
749
750 if (Abar != NULL)
751 fmpz_mod_mpoly_swap(Abar, s, ctx);
752
753 if (Bbar != NULL)
754 fmpz_mod_mpoly_swap(Bbar, t, ctx);
755
756 success = 1;
757
758 cleanup:
759
760 fmpz_mod_mpoly_clear(Ac, ctx);
761 fmpz_mod_mpoly_clear(Bc, ctx);
762 fmpz_mod_mpoly_clear(Gc, ctx);
763 fmpz_mod_mpoly_clear(s, ctx);
764 fmpz_mod_mpoly_clear(t, ctx);
765 fmpz_mod_mpoly_univar_clear(Ax, ctx);
766 fmpz_mod_mpoly_univar_clear(Bx, ctx);
767 fmpz_mod_mpoly_univar_clear(Gx, ctx);
768
769 return success;
770 }
771
772
773 /********************** Hit A and B with zippel ******************************/
774
_try_zippel(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)775 static int _try_zippel(
776 fmpz_mod_mpoly_t G,
777 fmpz_mod_mpoly_t Abar,
778 fmpz_mod_mpoly_t Bbar,
779 const fmpz_mod_mpoly_t A,
780 const fmpz_mod_mpoly_t B,
781 const mpoly_gcd_info_t I,
782 const fmpz_mod_mpoly_ctx_t ctx)
783 {
784 slong i, k;
785 slong m = I->mvars;
786 int success;
787 flint_bitcnt_t wbits;
788 fmpz_mod_mpoly_ctx_t lctx;
789 fmpz_mod_mpoly_t Al, Bl, Gl, Abarl, Bbarl;
790 fmpz_mod_mpoly_t Ac, Bc, Gc, Abarc, Bbarc;
791 slong max_deg;
792 flint_rand_t state;
793
794 FLINT_ASSERT(A->bits <= FLINT_BITS);
795 FLINT_ASSERT(B->bits <= FLINT_BITS);
796 FLINT_ASSERT(A->length > 0);
797 FLINT_ASSERT(B->length > 0);
798
799 if (!(I->can_use & MPOLY_GCD_USE_ZIPPEL))
800 return 0;
801
802 FLINT_ASSERT(m >= 2);
803
804 flint_randinit(state);
805
806 fmpz_mod_mpoly_ctx_init(lctx, m, ORD_LEX, fmpz_mod_ctx_modulus(ctx->ffinfo));
807
808 max_deg = 0;
809 for (i = 0; i < m; i++)
810 {
811 k = I->zippel_perm[i];
812 max_deg = FLINT_MAX(max_deg, I->Adeflate_deg[k]);
813 max_deg = FLINT_MAX(max_deg, I->Bdeflate_deg[k]);
814 }
815
816 wbits = 1 + FLINT_BIT_COUNT(max_deg);
817 wbits = mpoly_fix_bits(wbits, lctx->minfo);
818 FLINT_ASSERT(wbits <= FLINT_BITS);
819
820 fmpz_mod_mpoly_init3(Al, 0, wbits, lctx);
821 fmpz_mod_mpoly_init3(Bl, 0, wbits, lctx);
822 fmpz_mod_mpoly_init3(Gl, 0, wbits, lctx);
823 fmpz_mod_mpoly_init3(Abarl, 0, wbits, lctx);
824 fmpz_mod_mpoly_init3(Bbarl, 0, wbits, lctx);
825 fmpz_mod_mpoly_init3(Ac, 0, wbits, lctx);
826 fmpz_mod_mpoly_init3(Bc, 0, wbits, lctx);
827 fmpz_mod_mpoly_init3(Gc, 0, wbits, lctx);
828 fmpz_mod_mpoly_init3(Abarc, 0, wbits, lctx);
829 fmpz_mod_mpoly_init3(Bbarc, 0, wbits, lctx);
830
831 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Al, lctx, A, ctx,
832 I->zippel_perm, I->Amin_exp, I->Gstride);
833 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Bl, lctx, B, ctx,
834 I->zippel_perm, I->Bmin_exp, I->Gstride);
835
836 success = fmpz_mod_mpolyl_content(Ac, Al, 1, lctx) &&
837 fmpz_mod_mpolyl_content(Bc, Bl, 1, lctx);
838 if (!success)
839 goto cleanup;
840
841 success = _fmpz_mod_mpoly_gcd_algo(Gc, Abar == NULL ? NULL : Abarc,
842 Bbar == NULL ? NULL : Bbarc, Ac, Bc, lctx, MPOLY_GCD_USE_ALL);
843 if (!success)
844 goto cleanup;
845
846 success = fmpz_mod_mpoly_divides(Al, Al, Ac, lctx);
847 FLINT_ASSERT(success);
848
849 success = fmpz_mod_mpoly_divides(Bl, Bl, Bc, lctx);
850 FLINT_ASSERT(success);
851
852 fmpz_mod_mpoly_repack_bits_inplace(Al, wbits, lctx);
853 fmpz_mod_mpoly_repack_bits_inplace(Bl, wbits, lctx);
854
855 success = fmpz_mod_mpolyl_gcdp_zippel(Gl, Abarl, Bbarl, Al, Bl,
856 m - 1, lctx, state);
857 if (!success)
858 goto cleanup;
859
860 fmpz_mod_mpoly_mul(Gl, Gl, Gc, lctx);
861 fmpz_mod_mpoly_from_mpolyl_perm_inflate(G, I->Gbits, ctx, Gl, lctx,
862 I->zippel_perm, I->Gmin_exp, I->Gstride);
863 if (Abar != NULL)
864 {
865 fmpz_mod_mpoly_mul(Abarl, Abarl, Abarc, lctx);
866 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Abar, I->Abarbits, ctx, Abarl, lctx,
867 I->zippel_perm, I->Abarmin_exp, I->Gstride);
868 }
869
870 if (Bbar != NULL)
871 {
872 fmpz_mod_mpoly_mul(Bbarl, Bbarl, Bbarc, lctx);
873 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Bbar, I->Bbarbits, ctx, Bbarl, lctx,
874 I->zippel_perm, I->Bbarmin_exp, I->Gstride);
875 }
876
877 success = 1;
878
879 cleanup:
880
881 fmpz_mod_mpoly_clear(Al, lctx);
882 fmpz_mod_mpoly_clear(Bl, lctx);
883 fmpz_mod_mpoly_clear(Gl, lctx);
884 fmpz_mod_mpoly_clear(Abarl, lctx);
885 fmpz_mod_mpoly_clear(Bbarl, lctx);
886 fmpz_mod_mpoly_clear(Ac, lctx);
887 fmpz_mod_mpoly_clear(Bc, lctx);
888 fmpz_mod_mpoly_clear(Gc, lctx);
889 fmpz_mod_mpoly_clear(Abarc, lctx);
890 fmpz_mod_mpoly_clear(Bbarc, lctx);
891
892 fmpz_mod_mpoly_ctx_clear(lctx);
893 flint_randclear(state);
894
895 return success;
896 }
897
_try_zippel2(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)898 static int _try_zippel2(
899 fmpz_mod_mpoly_t G,
900 fmpz_mod_mpoly_t Abar,
901 fmpz_mod_mpoly_t Bbar,
902 const fmpz_mod_mpoly_t A,
903 const fmpz_mod_mpoly_t B,
904 const mpoly_gcd_info_t I,
905 const fmpz_mod_mpoly_ctx_t ctx)
906 {
907 slong i, k;
908 slong m = I->mvars;
909 int success;
910 flint_bitcnt_t wbits;
911 fmpz_mod_mpoly_ctx_t lctx;
912 fmpz_mod_mpoly_t Al, Bl, Gl, Abarl, Bbarl;
913 fmpz_mod_mpoly_t Al_lc, Bl_lc, Ac, Bc, Gc, Abarc, Bbarc, Gamma;
914 slong * tmp, * Gl_degs, * Al_degs, * Bl_degs, * Gamma_degs, * Gguess;
915 slong max_degree;
916
917 FLINT_ASSERT(A->bits <= FLINT_BITS);
918 FLINT_ASSERT(B->bits <= FLINT_BITS);
919 FLINT_ASSERT(A->length > 0);
920 FLINT_ASSERT(B->length > 0);
921
922 if (!(I->can_use & MPOLY_GCD_USE_ZIPPEL2))
923 return 0;
924
925 FLINT_ASSERT(m >= 3);
926
927 tmp = FLINT_ARRAY_ALLOC(5*m, slong);
928 Al_degs = tmp + 1*m;
929 Bl_degs = tmp + 2*m;
930 Gl_degs = tmp + 3*m;
931 Gamma_degs = tmp + 4*m;
932
933 fmpz_mod_mpoly_ctx_init(lctx, m, ORD_LEX, fmpz_mod_ctx_modulus(ctx->ffinfo));
934
935 max_degree = 0;
936 for (i = 0; i < m; i++)
937 {
938 k = I->zippel2_perm[i];
939
940 Gl_degs[i] = I->Gdeflate_deg_bound[k];
941
942 Al_degs[i] = I->Adeflate_deg[k];
943 max_degree = FLINT_MAX(max_degree, Al_degs[i]);
944
945 Bl_degs[i] = I->Bdeflate_deg[k];
946 max_degree = FLINT_MAX(max_degree, Bl_degs[i]);
947 }
948
949 wbits = 1 + FLINT_BIT_COUNT(max_degree);
950 wbits = mpoly_fix_bits(wbits, lctx->minfo);
951 FLINT_ASSERT(wbits <= FLINT_BITS);
952
953 fmpz_mod_mpoly_init3(Al, 0, wbits, lctx);
954 fmpz_mod_mpoly_init3(Bl, 0, wbits, lctx);
955 fmpz_mod_mpoly_init3(Gl, 0, wbits, lctx);
956 fmpz_mod_mpoly_init3(Abarl, 0, wbits, lctx);
957 fmpz_mod_mpoly_init3(Bbarl, 0, wbits, lctx);
958 fmpz_mod_mpoly_init3(Ac, 0, wbits, lctx);
959 fmpz_mod_mpoly_init3(Bc, 0, wbits, lctx);
960 fmpz_mod_mpoly_init3(Gc, 0, wbits, lctx);
961 fmpz_mod_mpoly_init3(Abarc, 0, wbits, lctx);
962 fmpz_mod_mpoly_init3(Bbarc, 0, wbits, lctx);
963 fmpz_mod_mpoly_init3(Gamma, 0, wbits, lctx);
964 fmpz_mod_mpoly_init3(Al_lc, 0, wbits, lctx);
965 fmpz_mod_mpoly_init3(Bl_lc, 0, wbits, lctx);
966
967 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Al, lctx, A, ctx,
968 I->zippel2_perm, I->Amin_exp, I->Gstride);
969 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Bl, lctx, B, ctx,
970 I->zippel2_perm, I->Bmin_exp, I->Gstride);
971
972 success = fmpz_mod_mpolyl_content(Ac, Al, 2, lctx) &&
973 fmpz_mod_mpolyl_content(Bc, Bl, 2, lctx);
974 if (!success)
975 goto cleanup;
976
977 success = _fmpz_mod_mpoly_gcd_algo(Gc, Abar == NULL ? NULL : Abarc,
978 Bbar == NULL ? NULL : Bbarc, Ac, Bc, lctx, MPOLY_GCD_USE_ALL);
979 if (!success)
980 goto cleanup;
981
982 fmpz_mod_mpoly_degrees_si(tmp, Ac, lctx);
983 for (i = 0; i < m; i++)
984 Al_degs[i] -= tmp[i];
985
986 success = fmpz_mod_mpoly_divides(Al, Al, Ac, lctx);
987 FLINT_ASSERT(success);
988
989 fmpz_mod_mpoly_degrees_si(tmp, Bc, lctx);
990 for (i = 0; i < m; i++)
991 Bl_degs[i] -= tmp[i];
992
993 success = fmpz_mod_mpoly_divides(Bl, Bl, Bc, lctx);
994 FLINT_ASSERT(success);
995
996 fmpz_mod_mpoly_degrees_si(tmp, Gc, lctx);
997 for (i = 0; i < m; i++)
998 Gl_degs[i] -= tmp[i];
999
1000 fmpz_mod_mpoly_repack_bits_inplace(Al, wbits, lctx);
1001 fmpz_mod_mpoly_repack_bits_inplace(Bl, wbits, lctx);
1002 fmpz_mod_mpolyl_lead_coeff(Al_lc, Al, 2, lctx);
1003 fmpz_mod_mpolyl_lead_coeff(Bl_lc, Bl, 2, lctx);
1004 success = fmpz_mod_mpoly_gcd(Gamma, Al_lc, Bl_lc, lctx);
1005 if (!success)
1006 goto cleanup;
1007 fmpz_mod_mpoly_repack_bits_inplace(Gamma, wbits, lctx);
1008
1009 fmpz_mod_mpoly_degrees_si(Gamma_degs, Gamma, lctx);
1010
1011 Gguess = I->Gdeflate_deg_bounds_are_nice ? Gl_degs : NULL;
1012
1013 success = fmpz_mod_mpolyl_gcd_zippel2_smprime(Gl, Gguess, Abarl, Bbarl,
1014 Al, Al_degs, Bl, Bl_degs, Gamma, Gamma_degs, lctx);
1015 if (!success)
1016 goto cleanup;
1017
1018 fmpz_mod_mpoly_mul(Gl, Gl, Gc, lctx);
1019 fmpz_mod_mpoly_from_mpolyl_perm_inflate(G, I->Gbits, ctx, Gl, lctx,
1020 I->zippel2_perm, I->Gmin_exp, I->Gstride);
1021 if (Abar != NULL)
1022 {
1023 fmpz_mod_mpoly_mul(Abarl, Abarl, Abarc, lctx);
1024 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Abar, I->Abarbits, ctx, Abarl, lctx,
1025 I->zippel2_perm, I->Abarmin_exp, I->Gstride);
1026 }
1027
1028 if (Bbar != NULL)
1029 {
1030 fmpz_mod_mpoly_mul(Bbarl, Bbarl, Bbarc, lctx);
1031 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Bbar, I->Bbarbits, ctx, Bbarl, lctx,
1032 I->zippel2_perm, I->Bbarmin_exp, I->Gstride);
1033 }
1034
1035 success = 1;
1036
1037 cleanup:
1038
1039 fmpz_mod_mpoly_clear(Al, lctx);
1040 fmpz_mod_mpoly_clear(Bl, lctx);
1041 fmpz_mod_mpoly_clear(Gl, lctx);
1042 fmpz_mod_mpoly_clear(Abarl, lctx);
1043 fmpz_mod_mpoly_clear(Bbarl, lctx);
1044 fmpz_mod_mpoly_clear(Ac, lctx);
1045 fmpz_mod_mpoly_clear(Bc, lctx);
1046 fmpz_mod_mpoly_clear(Gc, lctx);
1047 fmpz_mod_mpoly_clear(Abarc, lctx);
1048 fmpz_mod_mpoly_clear(Bbarc, lctx);
1049 fmpz_mod_mpoly_clear(Gamma, lctx);
1050 fmpz_mod_mpoly_clear(Al_lc, lctx);
1051 fmpz_mod_mpoly_clear(Bl_lc, lctx);
1052
1053 fmpz_mod_mpoly_ctx_clear(lctx);
1054
1055 flint_free(tmp);
1056
1057 return success;
1058 }
1059
1060
1061 /******************** Hit A and B with hensel lifting ************************/
_try_hensel(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)1062 static int _try_hensel(
1063 fmpz_mod_mpoly_t G,
1064 fmpz_mod_mpoly_t Abar,
1065 fmpz_mod_mpoly_t Bbar,
1066 const fmpz_mod_mpoly_t A,
1067 const fmpz_mod_mpoly_t B,
1068 const mpoly_gcd_info_t I,
1069 const fmpz_mod_mpoly_ctx_t ctx)
1070 {
1071 slong i, k;
1072 slong m = I->mvars;
1073 int success;
1074 flint_bitcnt_t wbits;
1075 fmpz_mod_mpoly_ctx_t lctx;
1076 fmpz_mod_mpoly_t Al, Bl, Gl, Abarl, Bbarl;
1077 fmpz_mod_mpoly_t Ac, Bc, Gc, Abarc, Bbarc;
1078 slong max_deg;
1079
1080 FLINT_ASSERT(A->bits <= FLINT_BITS);
1081 FLINT_ASSERT(B->bits <= FLINT_BITS);
1082 FLINT_ASSERT(A->length > 0);
1083 FLINT_ASSERT(B->length > 0);
1084
1085 if (!(I->can_use & MPOLY_GCD_USE_HENSEL))
1086 return 0;
1087
1088 FLINT_ASSERT(m >= 2);
1089
1090 fmpz_mod_mpoly_ctx_init(lctx, m, ORD_LEX, fmpz_mod_ctx_modulus(ctx->ffinfo));
1091
1092 max_deg = 0;
1093 for (i = 0; i < m; i++)
1094 {
1095 k = I->hensel_perm[i];
1096 max_deg = FLINT_MAX(max_deg, I->Adeflate_deg[k]);
1097 max_deg = FLINT_MAX(max_deg, I->Bdeflate_deg[k]);
1098 }
1099
1100 wbits = 1 + FLINT_BIT_COUNT(max_deg);
1101 wbits = mpoly_fix_bits(wbits, lctx->minfo);
1102 FLINT_ASSERT(wbits <= FLINT_BITS);
1103
1104 fmpz_mod_mpoly_init3(Al, 0, wbits, lctx);
1105 fmpz_mod_mpoly_init3(Bl, 0, wbits, lctx);
1106 fmpz_mod_mpoly_init3(Gl, 0, wbits, lctx);
1107 fmpz_mod_mpoly_init3(Abarl, 0, wbits, lctx);
1108 fmpz_mod_mpoly_init3(Bbarl, 0, wbits, lctx);
1109 fmpz_mod_mpoly_init3(Ac, 0, wbits, lctx);
1110 fmpz_mod_mpoly_init3(Bc, 0, wbits, lctx);
1111 fmpz_mod_mpoly_init3(Gc, 0, wbits, lctx);
1112 fmpz_mod_mpoly_init3(Abarc, 0, wbits, lctx);
1113 fmpz_mod_mpoly_init3(Bbarc, 0, wbits, lctx);
1114
1115 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Al, lctx, A, ctx,
1116 I->hensel_perm, I->Amin_exp, I->Gstride);
1117 fmpz_mod_mpoly_to_mpolyl_perm_deflate(Bl, lctx, B, ctx,
1118 I->hensel_perm, I->Bmin_exp, I->Gstride);
1119
1120 success = fmpz_mod_mpolyl_content(Ac, Al, 1, lctx) &&
1121 fmpz_mod_mpolyl_content(Bc, Bl, 1, lctx);
1122 if (!success)
1123 goto cleanup;
1124
1125 success = _fmpz_mod_mpoly_gcd_algo(Gc, Abar == NULL ? NULL : Abarc,
1126 Bbar == NULL ? NULL : Bbarc, Ac, Bc, lctx, MPOLY_GCD_USE_ALL);
1127 if (!success)
1128 goto cleanup;
1129
1130 success = fmpz_mod_mpoly_divides(Al, Al, Ac, lctx);
1131 FLINT_ASSERT(success);
1132
1133 success = fmpz_mod_mpoly_divides(Bl, Bl, Bc, lctx);
1134 FLINT_ASSERT(success);
1135
1136 fmpz_mod_mpoly_repack_bits_inplace(Al, wbits, lctx);
1137 fmpz_mod_mpoly_repack_bits_inplace(Bl, wbits, lctx);
1138
1139 max_deg = I->Gdeflate_deg_bound[I->hensel_perm[0]];
1140 success = fmpz_mod_mpolyl_gcd_hensel_smprime(Gl, max_deg, Abarl, Bbarl, Al, Bl, lctx);
1141 if (!success)
1142 goto cleanup;
1143
1144 fmpz_mod_mpoly_mul(Gl, Gl, Gc, lctx);
1145 fmpz_mod_mpoly_from_mpolyl_perm_inflate(G, I->Gbits, ctx, Gl, lctx,
1146 I->hensel_perm, I->Gmin_exp, I->Gstride);
1147 if (Abar != NULL)
1148 {
1149 fmpz_mod_mpoly_mul(Abarl, Abarl, Abarc, lctx);
1150 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Abar, I->Abarbits, ctx, Abarl, lctx,
1151 I->hensel_perm, I->Abarmin_exp, I->Gstride);
1152 }
1153
1154 if (Bbar != NULL)
1155 {
1156 fmpz_mod_mpoly_mul(Bbarl, Bbarl, Bbarc, lctx);
1157 fmpz_mod_mpoly_from_mpolyl_perm_inflate(Bbar, I->Bbarbits, ctx, Bbarl, lctx,
1158 I->hensel_perm, I->Bbarmin_exp, I->Gstride);
1159 }
1160
1161 success = 1;
1162
1163 cleanup:
1164
1165 fmpz_mod_mpoly_clear(Al, lctx);
1166 fmpz_mod_mpoly_clear(Bl, lctx);
1167 fmpz_mod_mpoly_clear(Gl, lctx);
1168 fmpz_mod_mpoly_clear(Abarl, lctx);
1169 fmpz_mod_mpoly_clear(Bbarl, lctx);
1170 fmpz_mod_mpoly_clear(Ac, lctx);
1171 fmpz_mod_mpoly_clear(Bc, lctx);
1172 fmpz_mod_mpoly_clear(Gc, lctx);
1173 fmpz_mod_mpoly_clear(Abarc, lctx);
1174 fmpz_mod_mpoly_clear(Bbarc, lctx);
1175
1176 fmpz_mod_mpoly_ctx_clear(lctx);
1177
1178 return success;
1179 }
1180
1181
1182 /*********************** Hit A and B with brown ******************************/
1183
_try_brown(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,mpoly_gcd_info_t I,const fmpz_mod_mpoly_ctx_t ctx)1184 static int _try_brown(
1185 fmpz_mod_mpoly_t G,
1186 fmpz_mod_mpoly_t Abar,
1187 fmpz_mod_mpoly_t Bbar,
1188 const fmpz_mod_mpoly_t A,
1189 const fmpz_mod_mpoly_t B,
1190 mpoly_gcd_info_t I,
1191 const fmpz_mod_mpoly_ctx_t ctx)
1192 {
1193 int success;
1194 slong m = I->mvars;
1195 flint_bitcnt_t wbits;
1196 fmpz_mod_mpoly_ctx_t nctx;
1197 fmpz_mod_mpolyn_t An, Bn, Gn, Abarn, Bbarn;
1198 fmpz_mod_poly_polyun_mpolyn_stack_t St;
1199
1200 if (!(I->can_use & MPOLY_GCD_USE_BROWN))
1201 return 0;
1202
1203 FLINT_ASSERT(m >= 2);
1204 FLINT_ASSERT(A->bits <= FLINT_BITS);
1205 FLINT_ASSERT(B->bits <= FLINT_BITS);
1206 FLINT_ASSERT(A->length > 0);
1207 FLINT_ASSERT(B->length > 0);
1208
1209 wbits = FLINT_MAX(A->bits, B->bits);
1210
1211 fmpz_mod_mpoly_ctx_init(nctx, m, ORD_LEX, fmpz_mod_ctx_modulus(ctx->ffinfo));
1212 fmpz_mod_mpolyn_init(An, wbits, nctx);
1213 fmpz_mod_mpolyn_init(Bn, wbits, nctx);
1214 fmpz_mod_mpolyn_init(Gn, wbits, nctx);
1215 fmpz_mod_mpolyn_init(Abarn, wbits, nctx);
1216 fmpz_mod_mpolyn_init(Bbarn, wbits, nctx);
1217 fmpz_mod_poly_stack_init(St->poly_stack);
1218 fmpz_mod_polyun_stack_init(St->polyun_stack);
1219 fmpz_mod_mpolyn_stack_init(St->mpolyn_stack, wbits, nctx);
1220
1221 fmpz_mod_mpoly_to_mpolyn_perm_deflate(An, nctx, A, ctx,
1222 I->brown_perm, I->Amin_exp, I->Gstride);
1223 fmpz_mod_mpoly_to_mpolyn_perm_deflate(Bn, nctx, B, ctx,
1224 I->brown_perm, I->Bmin_exp, I->Gstride);
1225
1226 FLINT_ASSERT(An->bits == wbits);
1227 FLINT_ASSERT(Bn->bits == wbits);
1228 FLINT_ASSERT(An->length > 0);
1229 FLINT_ASSERT(Bn->length > 0);
1230
1231 success = fmpz_mod_mpolyn_gcd_brown_smprime(Gn, Abarn, Bbarn, An, Bn,
1232 m - 1, nctx, I, St);
1233 if (!success)
1234 goto cleanup;
1235
1236 fmpz_mod_mpoly_from_mpolyn_perm_inflate(G, I->Gbits, ctx, Gn, nctx,
1237 I->brown_perm, I->Gmin_exp, I->Gstride);
1238
1239 if (Abar != NULL)
1240 fmpz_mod_mpoly_from_mpolyn_perm_inflate(Abar, I->Abarbits, ctx, Abarn, nctx,
1241 I->brown_perm, I->Abarmin_exp, I->Gstride);
1242
1243 if (Bbar != NULL)
1244 fmpz_mod_mpoly_from_mpolyn_perm_inflate(Bbar, I->Bbarbits, ctx, Bbarn, nctx,
1245 I->brown_perm, I->Bbarmin_exp, I->Gstride);
1246
1247 success = 1;
1248
1249 cleanup:
1250
1251 fmpz_mod_poly_stack_clear(St->poly_stack);
1252 fmpz_mod_polyun_stack_clear(St->polyun_stack);
1253 fmpz_mod_mpolyn_stack_clear(St->mpolyn_stack, nctx);
1254 fmpz_mod_mpolyn_clear(An, nctx);
1255 fmpz_mod_mpolyn_clear(Bn, nctx);
1256 fmpz_mod_mpolyn_clear(Gn, nctx);
1257 fmpz_mod_mpolyn_clear(Abarn, nctx);
1258 fmpz_mod_mpolyn_clear(Bbarn, nctx);
1259 fmpz_mod_mpoly_ctx_clear(nctx);
1260
1261 return success;
1262 }
1263
_try_nmod(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx,unsigned int algo)1264 static int _try_nmod(
1265 fmpz_mod_mpoly_t G,
1266 fmpz_mod_mpoly_t Abar,
1267 fmpz_mod_mpoly_t Bbar,
1268 const fmpz_mod_mpoly_t A,
1269 const fmpz_mod_mpoly_t B,
1270 const fmpz_mod_mpoly_ctx_t ctx,
1271 unsigned int algo)
1272 {
1273 int success;
1274 nmod_mpoly_ctx_t nctx;
1275 nmod_mpoly_t nG, nAbar, nBbar, nA, nB;
1276
1277 FLINT_ASSERT(fmpz_abs_fits_ui(fmpz_mod_ctx_modulus(ctx->ffinfo)));
1278
1279 *nctx->minfo = *ctx->minfo;
1280 nmod_init(&nctx->mod, fmpz_get_ui(fmpz_mod_ctx_modulus(ctx->ffinfo)));
1281 nmod_mpoly_init(nG, nctx);
1282 nmod_mpoly_init(nAbar, nctx);
1283 nmod_mpoly_init(nBbar, nctx);
1284 nmod_mpoly_init(nA, nctx);
1285 nmod_mpoly_init(nB, nctx);
1286
1287 _fmpz_mod_mpoly_get_nmod_mpoly(nA, nctx, A, ctx);
1288 _fmpz_mod_mpoly_get_nmod_mpoly(nB, nctx, B, ctx);
1289
1290 success = _nmod_mpoly_gcd_algo_small(nG, Abar == NULL ? NULL : nAbar,
1291 Bbar == NULL ? NULL : nBbar, nA, nB, nctx, algo);
1292
1293 _fmpz_mod_mpoly_set_nmod_mpoly(G, ctx, nG, nctx);
1294 if (Abar != NULL)
1295 _fmpz_mod_mpoly_set_nmod_mpoly(Abar, ctx, nAbar, nctx);
1296 if (Bbar != NULL)
1297 _fmpz_mod_mpoly_set_nmod_mpoly(Bbar, ctx, nBbar, nctx);
1298
1299 nmod_mpoly_clear(nG, nctx);
1300 nmod_mpoly_clear(nAbar, nctx);
1301 nmod_mpoly_clear(nBbar, nctx);
1302 nmod_mpoly_clear(nA, nctx);
1303 nmod_mpoly_clear(nB, nctx);
1304
1305 return success;
1306 }
1307
1308
1309 /*
1310 Both A and B have to be packed into bits <= FLINT_BITS
1311 return is 1 for success, 0 for failure.
1312 */
_fmpz_mod_mpoly_gcd_algo_small(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx,unsigned int algo)1313 static int _fmpz_mod_mpoly_gcd_algo_small(
1314 fmpz_mod_mpoly_t G,
1315 fmpz_mod_mpoly_t Abar, /* could be NULL */
1316 fmpz_mod_mpoly_t Bbar, /* could be NULL */
1317 const fmpz_mod_mpoly_t A,
1318 const fmpz_mod_mpoly_t B,
1319 const fmpz_mod_mpoly_ctx_t ctx,
1320 unsigned int algo)
1321 {
1322 int success;
1323 slong v_in_both;
1324 slong v_in_either;
1325 slong v_in_A_only;
1326 slong v_in_B_only;
1327 slong j;
1328 slong nvars = ctx->minfo->nvars;
1329 mpoly_gcd_info_t I;
1330 #if FLINT_WANT_ASSERT
1331 fmpz_mod_mpoly_t T, Asave, Bsave;
1332 #endif
1333
1334 FLINT_ASSERT(A->length > 0);
1335 FLINT_ASSERT(B->length > 0);
1336 FLINT_ASSERT(A->bits <= FLINT_BITS);
1337 FLINT_ASSERT(B->bits <= FLINT_BITS);
1338
1339 if (A->length == 1)
1340 return _do_monomial_gcd(G, Bbar, Abar, B, A, ctx);
1341
1342 if (B->length == 1)
1343 return _do_monomial_gcd(G, Abar, Bbar, A, B, ctx);
1344
1345 #if FLINT_WANT_ASSERT
1346 fmpz_mod_mpoly_init(T, ctx);
1347 fmpz_mod_mpoly_init(Asave, ctx);
1348 fmpz_mod_mpoly_init(Bsave, ctx);
1349 fmpz_mod_mpoly_set(Asave, A, ctx);
1350 fmpz_mod_mpoly_set(Bsave, B, ctx);
1351 #endif
1352
1353 mpoly_gcd_info_init(I, nvars);
1354
1355 /* entries of I are all now invalid */
1356
1357 I->Gbits = FLINT_MIN(A->bits, B->bits);
1358 I->Abarbits = A->bits;
1359 I->Bbarbits = B->bits;
1360
1361 mpoly_gcd_info_limits(I->Amax_exp, I->Amin_exp, I->Alead_count,
1362 I->Atail_count, A->exps, A->bits, A->length, ctx->minfo);
1363 mpoly_gcd_info_limits(I->Bmax_exp, I->Bmin_exp, I->Blead_count,
1364 I->Btail_count, B->exps, B->bits, B->length, ctx->minfo);
1365
1366 /* set ess(p) := p/term_content(p) */
1367
1368 /* check if the cofactors could be monomials, i.e. ess(A) == ess(B) */
1369 for (j = 0; j < nvars; j++)
1370 {
1371 if (I->Amax_exp[j] - I->Amin_exp[j] != I->Bmax_exp[j] - I->Bmin_exp[j])
1372 goto skip_monomial_cofactors;
1373 }
1374
1375 if (_try_monomial_cofactors(G, Abar, Bbar, A, B, ctx))
1376 goto successful;
1377
1378 skip_monomial_cofactors:
1379
1380 if (fmpz_abs_fits_ui(fmpz_mod_ctx_modulus(ctx->ffinfo)))
1381 {
1382 success = _try_nmod(G, Abar, Bbar, A, B, ctx, algo);
1383 goto cleanup;
1384 }
1385
1386 mpoly_gcd_info_stride(I->Gstride,
1387 A->exps, A->bits, A->length, I->Amax_exp, I->Amin_exp,
1388 B->exps, B->bits, B->length, I->Bmax_exp, I->Bmin_exp, ctx->minfo);
1389
1390 for (j = 0; j < nvars; j++)
1391 {
1392 ulong t = I->Gstride[j];
1393
1394 if (t == 0)
1395 {
1396 FLINT_ASSERT(I->Amax_exp[j] == I->Amin_exp[j] ||
1397 I->Bmax_exp[j] == I->Bmin_exp[j]);
1398 }
1399 else
1400 {
1401 FLINT_ASSERT((I->Amax_exp[j] - I->Amin_exp[j]) % t == 0);
1402 FLINT_ASSERT((I->Bmax_exp[j] - I->Bmin_exp[j]) % t == 0);
1403 }
1404
1405 I->Adeflate_deg[j] = t == 0 ? 0 : (I->Amax_exp[j] - I->Amin_exp[j])/t;
1406 I->Bdeflate_deg[j] = t == 0 ? 0 : (I->Bmax_exp[j] - I->Bmin_exp[j])/t;
1407
1408 t = FLINT_MIN(I->Amin_exp[j], I->Bmin_exp[j]);
1409 I->Gmin_exp[j] = t;
1410 I->Abarmin_exp[j] = I->Amin_exp[j] - t;
1411 I->Bbarmin_exp[j] = I->Bmin_exp[j] - t;
1412 }
1413
1414 /*
1415 The following are now valid:
1416 I->Amax_exp, I->Amin_exp, I->Alead_count, I->Atail_count,
1417 I->Bmax_exp, I->Bmin_exp, I->Blead_count, I->Btail_count,
1418 I->Gstride
1419 I->Adeflate_deg
1420 I->Bdeflate_deg
1421 I->Gmin_exp
1422 */
1423
1424 /* check if ess(A) and ess(B) have a variable v_in_both in common */
1425 v_in_both = -WORD(1);
1426 for (j = 0; j < nvars; j++)
1427 {
1428 if (I->Amax_exp[j] > I->Amin_exp[j] && I->Bmax_exp[j] > I->Bmin_exp[j])
1429 {
1430 v_in_both = j;
1431 break;
1432 }
1433 }
1434 if (v_in_both == -WORD(1))
1435 {
1436 _do_trivial(G, Abar, Bbar, A, B, I, ctx);
1437 goto successful;
1438 }
1439
1440 /* check if ess(A) and ess(B) depend on another variable v_in_either */
1441 FLINT_ASSERT(0 <= v_in_both && v_in_both < nvars);
1442
1443 v_in_either = -WORD(1);
1444 for (j = 0; j < nvars; j++)
1445 {
1446 if (j == v_in_both)
1447 continue;
1448
1449 if (I->Amax_exp[j] > I->Amin_exp[j] || I->Bmax_exp[j] > I->Bmin_exp[j])
1450 {
1451 v_in_either = j;
1452 break;
1453 }
1454 }
1455
1456 if (v_in_either == -WORD(1))
1457 {
1458 _do_univar(G, Abar, Bbar, A, B, v_in_both, I, ctx);
1459 goto successful;
1460 }
1461
1462 /* check if there is a variable in ess(A) that is not in ess(B) */
1463 v_in_A_only = -WORD(1);
1464 v_in_B_only = -WORD(1);
1465 for (j = 0; j < nvars; j++)
1466 {
1467 if (I->Amax_exp[j] > I->Amin_exp[j] && I->Bmax_exp[j] == I->Bmin_exp[j])
1468 {
1469 v_in_A_only = j;
1470 break;
1471 }
1472 if (I->Bmax_exp[j] > I->Bmin_exp[j] && I->Amax_exp[j] == I->Amin_exp[j])
1473 {
1474 v_in_B_only = j;
1475 break;
1476 }
1477 }
1478 if (v_in_A_only != -WORD(1))
1479 {
1480 success = _try_missing_var(G, I->Gbits, Abar, Bbar, v_in_A_only,
1481 A, I->Amin_exp[v_in_A_only],
1482 B, I->Bmin_exp[v_in_A_only],
1483 ctx);
1484 goto cleanup;
1485 }
1486 if (v_in_B_only != -WORD(1))
1487 {
1488 success = _try_missing_var(G, I->Gbits, Bbar, Abar, v_in_B_only,
1489 B, I->Bmin_exp[v_in_B_only],
1490 A, I->Amin_exp[v_in_B_only],
1491 ctx);
1492 goto cleanup;
1493 }
1494
1495 /*
1496 all variable are now either
1497 missing from both ess(A) and ess(B), or
1498 present in both ess(A) and ess(B)
1499 and there are at least two in the latter case
1500 */
1501
1502 mpoly_gcd_info_set_perm(I, A->length, B->length, ctx->minfo);
1503
1504 /* _set_estimates will probably calculate the correct total degrees */
1505 I->Adeflate_tdeg = I->Bdeflate_tdeg = -1;
1506
1507 _set_estimates(I, A, B, ctx);
1508
1509 /* everything in I is valid now */
1510
1511 /* check divisibility A/B and B/A */
1512 {
1513 int gcd_is_trivial = 1;
1514 int try_a = I->Gdeflate_deg_bounds_are_nice;
1515 int try_b = I->Gdeflate_deg_bounds_are_nice;
1516
1517 for (j = 0; j < nvars; j++)
1518 {
1519 if (I->Gdeflate_deg_bound[j] != 0)
1520 gcd_is_trivial = 0;
1521
1522 if (I->Adeflate_deg[j] != I->Gdeflate_deg_bound[j])
1523 try_a = 0;
1524
1525 if (I->Bdeflate_deg[j] != I->Gdeflate_deg_bound[j])
1526 try_b = 0;
1527 }
1528
1529 if (gcd_is_trivial)
1530 {
1531 _do_trivial(G, Abar, Bbar, A, B, I, ctx);
1532 goto successful;
1533 }
1534
1535 if (try_a && _try_divides(G, Bbar, Abar, B, A, ctx))
1536 goto successful;
1537
1538 if (try_b && _try_divides(G, Abar, Bbar, A, B, ctx))
1539 goto successful;
1540 }
1541
1542 if (algo == MPOLY_GCD_USE_PRS)
1543 {
1544 success = _try_prs(G, Abar, Bbar, A, B, I, ctx);
1545 goto cleanup;
1546 }
1547
1548 if (I->mvars < 3)
1549 {
1550 mpoly_gcd_info_measure_brown(I, A->length, B->length, ctx->minfo);
1551 mpoly_gcd_info_measure_hensel(I, A->length, B->length, ctx->minfo);
1552
1553 algo &= (MPOLY_GCD_USE_BROWN | MPOLY_GCD_USE_HENSEL);
1554
1555 if (algo == MPOLY_GCD_USE_BROWN)
1556 {
1557 success = _try_brown(G, Abar, Bbar, A, B, I, ctx);
1558 }
1559 else if (algo == MPOLY_GCD_USE_HENSEL)
1560 {
1561 success = _try_hensel(G, Abar, Bbar, A, B, I, ctx);
1562 }
1563 else
1564 {
1565 if (I->Adensity + I->Bdensity > 0.05)
1566 {
1567 success = _try_brown(G, Abar, Bbar, A, B, I, ctx) ||
1568 _try_hensel(G, Abar, Bbar, A, B, I, ctx);
1569 }
1570 else
1571 {
1572 success = _try_hensel(G, Abar, Bbar, A, B, I, ctx) ||
1573 _try_brown(G, Abar, Bbar, A, B, I, ctx);
1574 }
1575 }
1576
1577 goto cleanup;
1578 }
1579 else if (algo == MPOLY_GCD_USE_HENSEL)
1580 {
1581 mpoly_gcd_info_measure_hensel(I, A->length, B->length, ctx->minfo);
1582 success = _try_hensel(G, Abar, Bbar, A, B, I, ctx);
1583 goto cleanup;
1584 }
1585 else if (algo == MPOLY_GCD_USE_BROWN)
1586 {
1587 mpoly_gcd_info_measure_brown(I, A->length, B->length, ctx->minfo);
1588 success = _try_brown(G, Abar, Bbar, A, B, I, ctx);
1589 goto cleanup;
1590 }
1591 else if (algo == MPOLY_GCD_USE_ZIPPEL)
1592 {
1593 mpoly_gcd_info_measure_zippel(I, A->length, B->length, ctx->minfo);
1594 success = _try_zippel(G, Abar, Bbar, A, B, I, ctx);
1595 goto cleanup;
1596 }
1597 else if (algo == MPOLY_GCD_USE_ZIPPEL2)
1598 {
1599 mpoly_gcd_info_measure_zippel2(I, A->length, B->length, ctx->minfo);
1600 success = _try_zippel2(G, Abar, Bbar, A, B, I, ctx);
1601 goto cleanup;
1602 }
1603 else
1604 {
1605 slong k;
1606 int deg_is_small = 1;
1607 double density = I->Adensity + I->Bdensity;
1608
1609 /*
1610 mpoly gcd case.
1611 Only rule is that measure_X must be called before
1612 try_X is called or I->X_perm is accessed.
1613 */
1614
1615 mpoly_gcd_info_measure_hensel(I, A->length, B->length, ctx->minfo);
1616 mpoly_gcd_info_measure_brown(I, A->length, B->length, ctx->minfo);
1617 mpoly_gcd_info_measure_zippel(I, A->length, B->length, ctx->minfo);
1618 mpoly_gcd_info_measure_zippel2(I, A->length, B->length, ctx->minfo);
1619
1620 if (density > 0.08)
1621 {
1622 if (!deg_is_small && _try_hensel(G, Abar, Bbar, A, B, I, ctx))
1623 goto successful;
1624
1625 if (_try_brown(G, Abar, Bbar, A, B, I, ctx))
1626 goto successful;
1627 }
1628
1629 if (I->Adeflate_tdeg > 0 && I->Bdeflate_tdeg > 0)
1630 {
1631 fmpz_t x;
1632 double tdensity;
1633 fmpz_init(x);
1634 fmpz_bin_uiui(x, (ulong)I->Adeflate_tdeg + I->mvars, I->mvars);
1635 tdensity = A->length/fmpz_get_d(x);
1636 fmpz_bin_uiui(x, (ulong)I->Bdeflate_tdeg + I->mvars, I->mvars);
1637 tdensity += B->length/fmpz_get_d(x);
1638 density = FLINT_MAX(density, tdensity);
1639 fmpz_clear(x);
1640 }
1641
1642 if (density > (deg_is_small ? 0.05 : 0.001))
1643 {
1644 if (_try_hensel(G, Abar, Bbar, A, B, I, ctx))
1645 goto successful;
1646 }
1647
1648 k = I->zippel2_perm[1];
1649 k = FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
1650 if ((A->length + B->length)/64 < k)
1651 {
1652 if (_try_zippel(G, Abar, Bbar, A, B, I, ctx))
1653 goto successful;
1654 if (_try_zippel2(G, Abar, Bbar, A, B, I, ctx))
1655 goto successful;
1656 }
1657 else
1658 {
1659 if (_try_zippel2(G, Abar, Bbar, A, B, I, ctx))
1660 goto successful;
1661 if (_try_zippel(G, Abar, Bbar, A, B, I, ctx))
1662 goto successful;
1663 }
1664
1665 success = _try_hensel(G, Abar, Bbar, A, B, I, ctx) ||
1666 _try_brown(G, Abar, Bbar, A, B, I, ctx);
1667 goto cleanup;
1668 }
1669
1670 success = 0;
1671 goto cleanup;
1672
1673 successful:
1674
1675 success = 1;
1676
1677 cleanup:
1678
1679 mpoly_gcd_info_clear(I);
1680
1681 if (success)
1682 {
1683 FLINT_ASSERT(G->length > 0);
1684
1685 if (!fmpz_is_one(G->coeffs + 0))
1686 {
1687 if (Abar != NULL)
1688 _fmpz_mod_vec_scalar_mul_fmpz_mod(Abar->coeffs, Abar->coeffs,
1689 Abar->length, G->coeffs + 0, ctx->ffinfo);
1690 if (Bbar != NULL)
1691 _fmpz_mod_vec_scalar_mul_fmpz_mod(Bbar->coeffs, Bbar->coeffs,
1692 Bbar->length, G->coeffs + 0, ctx->ffinfo);
1693
1694 _fmpz_mod_vec_scalar_div_fmpz_mod(G->coeffs, G->coeffs, G->length,
1695 G->coeffs + 0, ctx->ffinfo);
1696 }
1697
1698 FLINT_ASSERT(fmpz_mod_mpoly_divides(T, Asave, G, ctx));
1699 FLINT_ASSERT(Abar == NULL || fmpz_mod_mpoly_equal(T, Abar, ctx));
1700
1701 FLINT_ASSERT(fmpz_mod_mpoly_divides(T, Bsave, G, ctx));
1702 FLINT_ASSERT(Bbar == NULL || fmpz_mod_mpoly_equal(T, Bbar, ctx));
1703 }
1704
1705 #if FLINT_WANT_ASSERT
1706 fmpz_mod_mpoly_clear(T, ctx);
1707 fmpz_mod_mpoly_clear(Asave, ctx);
1708 fmpz_mod_mpoly_clear(Bsave, ctx);
1709 #endif
1710
1711 return success;
1712 }
1713
1714
1715 /*
1716 The gcd calculation is unusual.
1717 First see if both inputs fit into FLINT_BITS.
1718 Then, try deflation as a last resort.
1719 */
_fmpz_mod_mpoly_gcd_algo_large(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx,unsigned int algo)1720 static int _fmpz_mod_mpoly_gcd_algo_large(
1721 fmpz_mod_mpoly_t G,
1722 fmpz_mod_mpoly_t Abar,
1723 fmpz_mod_mpoly_t Bbar,
1724 const fmpz_mod_mpoly_t A,
1725 const fmpz_mod_mpoly_t B,
1726 const fmpz_mod_mpoly_ctx_t ctx,
1727 unsigned int algo)
1728 {
1729 int success;
1730 slong k;
1731 fmpz * Ashift, * Astride;
1732 fmpz * Bshift, * Bstride;
1733 fmpz * Gshift, * Gstride;
1734 fmpz_mod_mpoly_t Anew, Bnew;
1735 const fmpz_mod_mpoly_struct * Ause, * Buse;
1736
1737 if (A->length == 1)
1738 return _do_monomial_gcd(G, Bbar, Abar, B, A, ctx);
1739
1740 if (B->length == 1)
1741 return _do_monomial_gcd(G, Abar, Bbar, A, B, ctx);
1742
1743 if (_try_monomial_cofactors(G, Abar, Bbar, A, B, ctx))
1744 return 1;
1745
1746 fmpz_mod_mpoly_init(Anew, ctx);
1747 fmpz_mod_mpoly_init(Bnew, ctx);
1748
1749 Ause = A;
1750 if (A->bits > FLINT_BITS)
1751 {
1752 if (!fmpz_mod_mpoly_repack_bits(Anew, A, FLINT_BITS, ctx))
1753 goto could_not_repack;
1754 Ause = Anew;
1755 }
1756
1757 Buse = B;
1758 if (B->bits > FLINT_BITS)
1759 {
1760 if (!fmpz_mod_mpoly_repack_bits(Bnew, B, FLINT_BITS, ctx))
1761 goto could_not_repack;
1762 Buse = Bnew;
1763 }
1764
1765 success = _fmpz_mod_mpoly_gcd_algo_small(G, Abar, Bbar, Ause, Buse, ctx, algo);
1766 goto cleanup;
1767
1768 could_not_repack:
1769
1770 /*
1771 One of A or B could not be repacked into FLINT_BITS. See if
1772 they both fit into FLINT_BITS after deflation.
1773 */
1774
1775 Ashift = _fmpz_vec_init(ctx->minfo->nvars);
1776 Astride = _fmpz_vec_init(ctx->minfo->nvars);
1777 Bshift = _fmpz_vec_init(ctx->minfo->nvars);
1778 Bstride = _fmpz_vec_init(ctx->minfo->nvars);
1779 Gshift = _fmpz_vec_init(ctx->minfo->nvars);
1780 Gstride = _fmpz_vec_init(ctx->minfo->nvars);
1781
1782 fmpz_mod_mpoly_deflation(Ashift, Astride, A, ctx);
1783 fmpz_mod_mpoly_deflation(Bshift, Bstride, B, ctx);
1784 _fmpz_vec_min(Gshift, Ashift, Bshift, ctx->minfo->nvars);
1785 for (k = 0; k < ctx->minfo->nvars; k++)
1786 fmpz_gcd(Gstride + k, Astride + k, Bstride + k);
1787
1788 fmpz_mod_mpoly_deflate(Anew, A, Ashift, Gstride, ctx);
1789 if (Anew->bits > FLINT_BITS)
1790 {
1791 success = fmpz_mod_mpoly_repack_bits(Anew, Anew, FLINT_BITS, ctx);
1792 if (!success)
1793 goto deflate_cleanup;
1794 }
1795
1796 fmpz_mod_mpoly_deflate(Bnew, B, Bshift, Gstride, ctx);
1797 if (Bnew->bits > FLINT_BITS)
1798 {
1799 success = fmpz_mod_mpoly_repack_bits(Bnew, Bnew, FLINT_BITS, ctx);
1800 if (!success)
1801 goto deflate_cleanup;
1802 }
1803
1804 success = _fmpz_mod_mpoly_gcd_algo_small(G, Abar, Bbar, Anew, Bnew, ctx, algo);
1805 if (!success)
1806 goto deflate_cleanup;
1807
1808 for (k = 0; k < ctx->minfo->nvars; k++)
1809 {
1810 fmpz_sub(Ashift + k, Ashift + k, Gshift + k);
1811 fmpz_sub(Bshift + k, Bshift + k, Gshift + k);
1812 FLINT_ASSERT(fmpz_sgn(Ashift + k) >= 0);
1813 FLINT_ASSERT(fmpz_sgn(Bshift + k) >= 0);
1814 }
1815
1816 fmpz_mod_mpoly_inflate(G, G, Gshift, Gstride, ctx);
1817 if (Abar != NULL)
1818 fmpz_mod_mpoly_inflate(Abar, Abar, Ashift, Gstride, ctx);
1819 if (Bbar != NULL)
1820 fmpz_mod_mpoly_inflate(Bbar, Bbar, Bshift, Gstride, ctx);
1821
1822 FLINT_ASSERT(G->length > 0);
1823 if (!fmpz_is_one(G->coeffs + 0))
1824 {
1825 if (Abar != NULL)
1826 _fmpz_mod_vec_scalar_mul_fmpz_mod(Abar->coeffs, Abar->coeffs,
1827 Abar->length, G->coeffs + 0, ctx->ffinfo);
1828
1829 if (Bbar != NULL)
1830 _fmpz_mod_vec_scalar_mul_fmpz_mod(Bbar->coeffs, Bbar->coeffs,
1831 Bbar->length, G->coeffs + 0, ctx->ffinfo);
1832
1833 _fmpz_mod_vec_scalar_div_fmpz_mod(G->coeffs, G->coeffs, G->length,
1834 G->coeffs + 0, ctx->ffinfo);
1835 }
1836
1837 deflate_cleanup:
1838
1839 _fmpz_vec_clear(Ashift, ctx->minfo->nvars);
1840 _fmpz_vec_clear(Astride, ctx->minfo->nvars);
1841 _fmpz_vec_clear(Bshift, ctx->minfo->nvars);
1842 _fmpz_vec_clear(Bstride, ctx->minfo->nvars);
1843 _fmpz_vec_clear(Gshift, ctx->minfo->nvars);
1844 _fmpz_vec_clear(Gstride, ctx->minfo->nvars);
1845
1846 cleanup:
1847
1848 fmpz_mod_mpoly_clear(Anew, ctx);
1849 fmpz_mod_mpoly_clear(Bnew, ctx);
1850
1851 return success;
1852 }
1853
1854
_fmpz_mod_mpoly_gcd_algo(fmpz_mod_mpoly_t G,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx,unsigned int algo)1855 int _fmpz_mod_mpoly_gcd_algo(
1856 fmpz_mod_mpoly_t G,
1857 fmpz_mod_mpoly_t Abar,
1858 fmpz_mod_mpoly_t Bbar,
1859 const fmpz_mod_mpoly_t A,
1860 const fmpz_mod_mpoly_t B,
1861 const fmpz_mod_mpoly_ctx_t ctx,
1862 unsigned int algo)
1863 {
1864 FLINT_ASSERT(!fmpz_mod_mpoly_is_zero(A, ctx));
1865 FLINT_ASSERT(!fmpz_mod_mpoly_is_zero(B, ctx));
1866
1867 if (A->bits <= FLINT_BITS && B->bits <= FLINT_BITS)
1868 return _fmpz_mod_mpoly_gcd_algo_small(G, Abar, Bbar, A, B, ctx, algo);
1869 else
1870 return _fmpz_mod_mpoly_gcd_algo_large(G, Abar, Bbar, A, B, ctx, algo);
1871 }
1872
1873