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