1 /*
2 Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
10 */
11
12 #include "fq_nmod_mpoly.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 */
fq_nmod_mpoly_evals(fq_nmod_poly_struct * out,const int * ignore,const fq_nmod_mpoly_t A,ulong * Amin_exp,ulong * Amax_exp,ulong * Astride,fq_nmod_struct * alpha,const fq_nmod_mpoly_ctx_t ctx)24 void fq_nmod_mpoly_evals(
25 fq_nmod_poly_struct * out,
26 const int * ignore,
27 const fq_nmod_mpoly_t A,
28 ulong * Amin_exp,
29 ulong * Amax_exp,
30 ulong * Astride,
31 fq_nmod_struct * alpha,
32 const fq_nmod_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 fq_nmod_struct * Acoeff = A->coeffs;
44 fq_nmod_t meval, t, t2;
45
46 FLINT_ASSERT(A->bits <= FLINT_BITS);
47
48 fq_nmod_init(meval, ctx->fqctx);
49 fq_nmod_init(t, ctx->fqctx);
50 fq_nmod_init(t2, ctx->fqctx);
51
52 mask = (-UWORD(1)) >> (FLINT_BITS - A->bits);
53 offsets = (slong *) flint_malloc(ctx->minfo->nvars*sizeof(slong));
54 shifts = (slong *) flint_malloc(ctx->minfo->nvars*sizeof(slong));
55
56 for (j = 0; j < ctx->minfo->nvars; j++)
57 {
58 fq_nmod_poly_zero(out + j, ctx->fqctx);
59 mpoly_gen_offset_shift_sp(offsets + j, shifts + j, j, A->bits, ctx->minfo);
60 }
61
62 /*
63 two cases:
64 (1) the Amax_exp[j] are small enough to calculate a direct LUT
65 (2) use a LUT for exponents that are powers of two
66 */
67
68 total_limit = A->length/256;
69 total_limit = FLINT_MAX(WORD(9999), total_limit);
70 total_length = 0;
71 use_direct_LUT = 1;
72 for (j = 0; j < ctx->minfo->nvars; j++)
73 {
74 total_length += Amax_exp[j] + 1;
75 if ((ulong) total_length > (ulong) total_limit)
76 use_direct_LUT = 0;
77 }
78
79 if (use_direct_LUT)
80 {
81 slong off;
82 fq_nmod_struct * LUT, ** LUTvalue, ** LUTvalueinv;
83
84 /* value of powers of alpha[j] */
85 LUT = (fq_nmod_struct *) flint_malloc(2*total_length*sizeof(fq_nmod_struct));
86 for (j = 0; j < 2*total_length; j++)
87 fq_nmod_init(LUT + j, ctx->fqctx);
88
89 /* pointers into LUT */
90 LUTvalue = (fq_nmod_struct **) flint_malloc(nvars*sizeof(fq_nmod_struct *));
91 LUTvalueinv = (fq_nmod_struct **) flint_malloc(nvars*sizeof(fq_nmod_struct *));
92
93 off = 0;
94 for (j = 0; j < nvars; j++)
95 {
96 ulong k;
97
98 fq_nmod_inv(t2, alpha + j, ctx->fqctx);
99
100 LUTvalue[j] = LUT + off;
101 LUTvalueinv[j] = LUT + total_length + off;
102 fq_nmod_one(LUTvalue[j] + 0, ctx->fqctx);
103 fq_nmod_one(LUTvalueinv[j] + 0, ctx->fqctx);
104 for (k = 0; k < Amax_exp[j]; k++)
105 {
106 fq_nmod_mul(LUTvalue[j] + k + 1, LUTvalue[j] + k,
107 alpha + j, ctx->fqctx);
108 fq_nmod_mul(LUTvalueinv[j] + k + 1, LUTvalueinv[j] + k,
109 t2, ctx->fqctx);
110 }
111
112 off += Amax_exp[j] + 1;
113 }
114 FLINT_ASSERT(off == total_length);
115
116 for (i = 0; i < A->length; i++)
117 {
118 fq_nmod_set(meval, Acoeff + i, ctx->fqctx);
119
120 for (j = 0; j < nvars; j++)
121 {
122 varexp = ((Aexp + N*i)[offsets[j]]>>shifts[j])&mask;
123 FLINT_ASSERT(varexp <= Amax_exp[j]);
124 fq_nmod_mul(t, meval, LUTvalue[j] + varexp, ctx->fqctx);
125 fq_nmod_swap(meval, t, ctx->fqctx);
126 }
127
128 for (j = 0; j < nvars; j++)
129 {
130 varexp = ((Aexp + N*i)[offsets[j]]>>shifts[j])&mask;
131
132 if (ignore[j])
133 continue;
134
135 fq_nmod_mul(t, meval, LUTvalueinv[j] + varexp, ctx->fqctx);
136
137 FLINT_ASSERT((Astride[j] == 0 && varexp == Amin_exp[j])
138 || (varexp - Amin_exp[j]) % Astride[j] == 0);
139
140 varexp = Astride[j] < 2 ? varexp - Amin_exp[j] :
141 (varexp - Amin_exp[j])/Astride[j];
142
143 fq_nmod_poly_get_coeff(t2, out + j, varexp, ctx->fqctx);
144 fq_nmod_add(t, t, t2, ctx->fqctx);
145 fq_nmod_poly_set_coeff(out + j, varexp, t, ctx->fqctx);
146 }
147 }
148
149 for (j = 0; j < 2*total_length; j++)
150 fq_nmod_clear(LUT + j, ctx->fqctx);
151
152 flint_free(LUT);
153 flint_free(LUTvalue);
154 flint_free(LUTvalueinv);
155 }
156 else
157 {
158 slong LUTlen;
159 ulong * LUTmask;
160 slong * LUToffset, * LUTvar;
161 fq_nmod_struct * LUTvalue, * LUTvalueinv;
162 fq_nmod_struct * vieval;
163 fq_nmod_t xpoweval, xinvpoweval;
164
165 fq_nmod_init(xpoweval, ctx->fqctx);
166 fq_nmod_init(xinvpoweval, ctx->fqctx);
167
168 LUToffset = (slong *) flint_malloc(N*FLINT_BITS*sizeof(slong));
169 LUTmask = (ulong *) flint_malloc(N*FLINT_BITS*sizeof(ulong));
170 LUTvalue = (fq_nmod_struct *) flint_malloc(N*FLINT_BITS*sizeof(fq_nmod_struct));
171 LUTvar = (slong *) flint_malloc(N*FLINT_BITS*sizeof(slong));
172 LUTvalueinv = (fq_nmod_struct *) flint_malloc(N*FLINT_BITS*sizeof(fq_nmod_struct));
173 for (j = 0; j < N*FLINT_BITS; j++)
174 {
175 fq_nmod_init(LUTvalue + j, ctx->fqctx);
176 fq_nmod_init(LUTvalueinv + j, ctx->fqctx);
177 }
178
179 vieval = (fq_nmod_struct *) flint_malloc(nvars*sizeof(fq_nmod_struct));
180 for (j = 0; j < nvars; j++)
181 {
182 fq_nmod_init(vieval + j, ctx->fqctx);
183 }
184
185 LUTlen = 0;
186 for (j = nvars - 1; j >= 0; j--)
187 {
188 flint_bitcnt_t bits = FLINT_BIT_COUNT(Amax_exp[j]);
189 fq_nmod_set(xpoweval, alpha + j, ctx->fqctx); /* xpoweval = alpha[j]^(2^i) */
190 fq_nmod_inv(xinvpoweval, xpoweval, ctx->fqctx); /* alpha[j]^(-2^i) */
191 for (i = 0; i < bits; i++)
192 {
193 LUToffset[LUTlen] = offsets[j];
194 LUTmask[LUTlen] = (UWORD(1) << (shifts[j] + i));
195 fq_nmod_set(LUTvalue + LUTlen, xpoweval, ctx->fqctx);
196 fq_nmod_set(LUTvalueinv + LUTlen, xinvpoweval, ctx->fqctx);
197 LUTvar[LUTlen] = j;
198 LUTlen++;
199 fq_nmod_mul(xpoweval, xpoweval, xpoweval, ctx->fqctx);
200 fq_nmod_mul(xinvpoweval, xinvpoweval, xinvpoweval, ctx->fqctx);
201 }
202
203 fq_nmod_one(vieval + j, ctx->fqctx);
204 }
205 FLINT_ASSERT(LUTlen < N*FLINT_BITS);
206
207 for (i = 0; i < A->length; i++)
208 {
209 fq_nmod_set(meval, Acoeff + i, ctx->fqctx);
210
211 for (j = 0; j < LUTlen; j++)
212 {
213 if (((Aexp + N*i)[LUToffset[j]] & LUTmask[j]) != 0)
214 {
215 fq_nmod_mul(meval, meval, LUTvalue + j, ctx->fqctx);
216 fq_nmod_mul(vieval + LUTvar[j], vieval + LUTvar[j],
217 LUTvalueinv + j, ctx->fqctx);
218 }
219 }
220
221 for (j = 0; j < nvars; j++)
222 {
223 varexp = ((Aexp + N*i)[offsets[j]]>>shifts[j])&mask;
224
225 FLINT_ASSERT((Astride[j] == 0 && varexp == Amin_exp[j])
226 || (varexp - Amin_exp[j]) % Astride[j] == 0);
227
228 varexp = Astride[j] < 2 ? varexp - Amin_exp[j] :
229 (varexp - Amin_exp[j])/Astride[j];
230
231 fq_nmod_mul(t, meval, vieval + j, ctx->fqctx);
232 fq_nmod_poly_get_coeff(t2, out + j, varexp, ctx->fqctx);
233 fq_nmod_add(t, t, t2, ctx->fqctx);
234 fq_nmod_poly_set_coeff(out + j, varexp, t, ctx->fqctx);
235 fq_nmod_one(vieval + j, ctx->fqctx);
236 }
237 }
238
239 for (j = 0; j < N*FLINT_BITS; j++)
240 {
241 fq_nmod_clear(LUTvalue + j, ctx->fqctx);
242 fq_nmod_clear(LUTvalueinv + j, ctx->fqctx);
243 }
244 flint_free(LUToffset);
245 flint_free(LUTmask);
246 flint_free(LUTvalue);
247 flint_free(LUTvar);
248 flint_free(LUTvalueinv);
249
250 for (j = 0; j < nvars; j++)
251 {
252 fq_nmod_clear(vieval + j, ctx->fqctx);
253 }
254 flint_free(vieval);
255
256 fq_nmod_clear(xpoweval, ctx->fqctx);
257 fq_nmod_clear(xinvpoweval, ctx->fqctx);
258 }
259
260 flint_free(offsets);
261 flint_free(shifts);
262
263 fq_nmod_clear(meval, ctx->fqctx);
264 fq_nmod_clear(t, ctx->fqctx);
265 fq_nmod_clear(t2, ctx->fqctx);
266 }
267
268
mpoly_gcd_info_set_estimates_fq_nmod_mpoly(mpoly_gcd_info_t I,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)269 void mpoly_gcd_info_set_estimates_fq_nmod_mpoly(
270 mpoly_gcd_info_t I,
271 const fq_nmod_mpoly_t A,
272 const fq_nmod_mpoly_t B,
273 const fq_nmod_mpoly_ctx_t ctx)
274 {
275 int try_count = 0;
276 slong i, j;
277 fq_nmod_poly_t Geval;
278 fq_nmod_poly_struct * Aevals, * Bevals;
279 fq_nmod_struct * alpha;
280 flint_rand_t randstate;
281 slong ignore_limit;
282 int * ignore;
283
284 flint_randinit(randstate);
285
286 ignore = (int *) flint_malloc(ctx->minfo->nvars*sizeof(int));
287 alpha = (fq_nmod_struct *) flint_malloc(
288 ctx->minfo->nvars*sizeof(fq_nmod_struct));
289 Aevals = (fq_nmod_poly_struct *) flint_malloc(
290 ctx->minfo->nvars*sizeof(fq_nmod_poly_struct));
291 Bevals = (fq_nmod_poly_struct *) flint_malloc(
292 ctx->minfo->nvars*sizeof(fq_nmod_poly_struct));
293
294 fq_nmod_poly_init(Geval, ctx->fqctx);
295 for (j = 0; j < ctx->minfo->nvars; j++)
296 {
297 fq_nmod_init(alpha + j, ctx->fqctx);
298 fq_nmod_poly_init(Aevals + j, ctx->fqctx);
299 fq_nmod_poly_init(Bevals + j, ctx->fqctx);
300 }
301
302 ignore_limit = A->length/4096 + B->length/4096;
303 ignore_limit = FLINT_MAX(WORD(9999), ignore_limit);
304 I->Gdeflate_deg_bounds_are_nice = 1;
305 for (j = 0; j < ctx->minfo->nvars; j++)
306 {
307 if ( I->Adeflate_deg[j] > ignore_limit
308 || I->Bdeflate_deg[j] > ignore_limit)
309 {
310 ignore[j] = 1;
311 I->Gdeflate_deg_bounds_are_nice = 0;
312 }
313 else
314 {
315 ignore[j] = 0;
316 }
317 }
318
319 try_again:
320
321 if (++try_count > 10)
322 {
323 I->Gdeflate_deg_bounds_are_nice = 0;
324 for (j = 0; j < ctx->minfo->nvars; j++)
325 {
326 I->Gdeflate_deg_bound[j] = FLINT_MIN(I->Adeflate_deg[j],
327 I->Bdeflate_deg[j]);
328 I->Gterm_count_est[j] = 1 + I->Gdeflate_deg_bound[j]/2;
329 }
330
331 goto cleanup;
332 }
333
334 for (j = 0; j < ctx->minfo->nvars; j++)
335 {
336 fq_nmod_randtest_not_zero(alpha + j, randstate, ctx->fqctx);
337 }
338
339
340 fq_nmod_mpoly_evals(Aevals, ignore, A, I->Amin_exp, I->Amax_exp,
341 I->Gstride, alpha, ctx);
342 fq_nmod_mpoly_evals(Bevals, ignore, B, I->Bmin_exp, I->Bmax_exp,
343 I->Gstride, alpha, ctx);
344
345 for (j = 0; j < ctx->minfo->nvars; j++)
346 {
347 if (ignore[j])
348 {
349 I->Gdeflate_deg_bound[j] = FLINT_MIN(I->Adeflate_deg[j],
350 I->Bdeflate_deg[j]);
351 I->Gterm_count_est[j] = 1 + I->Gdeflate_deg_bound[j]/2;
352 }
353 else
354 {
355 if ( I->Adeflate_deg[j] != fq_nmod_poly_degree(Aevals + j, ctx->fqctx)
356 || I->Bdeflate_deg[j] != fq_nmod_poly_degree(Bevals + j, ctx->fqctx))
357 {
358 goto try_again;
359 }
360
361 fq_nmod_poly_gcd(Geval, Aevals + j, Bevals + j, ctx->fqctx);
362
363 I->Gterm_count_est[j] = 0;
364 I->Gdeflate_deg_bound[j] = fq_nmod_poly_degree(Geval, ctx->fqctx);
365 for (i = I->Gdeflate_deg_bound[j]; i >= 0; i--)
366 {
367 I->Gterm_count_est[j] += !fq_nmod_is_zero(Geval->coeffs + i, ctx->fqctx);
368 }
369 }
370 }
371
372 cleanup:
373
374 fq_nmod_poly_clear(Geval, ctx->fqctx);
375 for (j = 0; j < ctx->minfo->nvars; j++)
376 {
377 fq_nmod_clear(alpha + j, ctx->fqctx);
378 fq_nmod_poly_clear(Aevals + j, ctx->fqctx);
379 fq_nmod_poly_clear(Bevals + j, ctx->fqctx);
380 }
381
382 flint_free(ignore);
383 flint_free(alpha);
384 flint_free(Aevals);
385 flint_free(Bevals);
386
387 flint_randclear(randstate);
388
389 return;
390 }
391
392 /*********************** Easy when B is a monomial ***************************/
_try_monomial_gcd(fq_nmod_mpoly_t G,flint_bitcnt_t Gbits,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)393 static int _try_monomial_gcd(
394 fq_nmod_mpoly_t G, flint_bitcnt_t Gbits,
395 const fq_nmod_mpoly_t A,
396 const fq_nmod_mpoly_t B,
397 const fq_nmod_mpoly_ctx_t ctx)
398 {
399 slong i;
400 fmpz * minAfields, * minAdegs, * minBdegs;
401 TMP_INIT;
402
403 FLINT_ASSERT(A->length > 0);
404 FLINT_ASSERT(B->length == 1);
405
406 TMP_START;
407
408 /* get the field-wise minimum of A */
409 minAfields = (fmpz *) TMP_ALLOC(ctx->minfo->nfields*sizeof(fmpz));
410 for (i = 0; i < ctx->minfo->nfields; i++)
411 fmpz_init(minAfields + i);
412 mpoly_min_fields_fmpz(minAfields, A->exps, A->length, A->bits, ctx->minfo);
413
414 /* unpack to get the min degrees of each variable in A */
415 minAdegs = (fmpz *) TMP_ALLOC(ctx->minfo->nvars*sizeof(fmpz));
416 for (i = 0; i < ctx->minfo->nvars; i++)
417 fmpz_init(minAdegs + i);
418 mpoly_get_monomial_ffmpz_unpacked_ffmpz(minAdegs, minAfields, ctx->minfo);
419
420 /* get the degree of each variable in B */
421 minBdegs = (fmpz *) TMP_ALLOC(ctx->minfo->nvars*sizeof(fmpz));
422 for (i = 0; i < ctx->minfo->nvars; i++)
423 fmpz_init(minBdegs + i);
424 mpoly_get_monomial_ffmpz(minBdegs, B->exps, B->bits, ctx->minfo);
425
426 /* compute the degree of each variable in G */
427 _fmpz_vec_min_inplace(minBdegs, minAdegs, ctx->minfo->nvars);
428
429 fq_nmod_mpoly_fit_length(G, 1, ctx);
430 fq_nmod_mpoly_fit_bits(G, Gbits, ctx);
431 G->bits = Gbits;
432 mpoly_set_monomial_ffmpz(G->exps, minBdegs, Gbits, ctx->minfo);
433 fq_nmod_one(G->coeffs + 0, ctx->fqctx);
434 G->length = 1;
435
436 for (i = 0; i < ctx->minfo->nfields; i++)
437 {
438 fmpz_clear(minAfields + i);
439 }
440 for (i = 0; i < ctx->minfo->nvars; i++)
441 {
442 fmpz_clear(minAdegs + i);
443 fmpz_clear(minBdegs + i);
444 }
445
446 TMP_END;
447
448 return 1;
449 }
450
451
452 /********************** See if cofactors are monomials ***********************/
_try_monomial_cofactors(fq_nmod_mpoly_t G,flint_bitcnt_t Gbits,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)453 static int _try_monomial_cofactors(
454 fq_nmod_mpoly_t G, flint_bitcnt_t Gbits,
455 const fq_nmod_mpoly_t A,
456 const fq_nmod_mpoly_t B,
457 const fq_nmod_mpoly_ctx_t ctx)
458 {
459 int success;
460 slong i, j;
461 slong NA, NG;
462 slong nvars = ctx->minfo->nvars;
463 fmpz * Abarexps, * Bbarexps, * Texps;
464 fq_nmod_t t1, t2;
465 fq_nmod_mpoly_t T;
466 TMP_INIT;
467
468 FLINT_ASSERT(A->length > 0);
469 FLINT_ASSERT(B->length > 0);
470
471 if (A->length != B->length)
472 return 0;
473
474 fq_nmod_init(t1, ctx->fqctx);
475 fq_nmod_init(t2, ctx->fqctx);
476
477 for (i = A->length - 1; i > 0; i--)
478 {
479 fq_nmod_mul(t1, A->coeffs + 0, B->coeffs + i, ctx->fqctx);
480 fq_nmod_mul(t2, B->coeffs + 0, A->coeffs + i, ctx->fqctx);
481 success = fq_nmod_equal(t1, t2, ctx->fqctx);
482 if (!success)
483 goto cleanup;
484 }
485
486 TMP_START;
487
488 Abarexps = (fmpz *) TMP_ALLOC(3*nvars*sizeof(fmpz));
489 Bbarexps = Abarexps + 1*nvars;
490 Texps = Abarexps + 2*nvars;
491 for (j = 0; j < nvars; j++)
492 {
493 fmpz_init(Abarexps + j);
494 fmpz_init(Bbarexps + j);
495 fmpz_init(Texps + j);
496 }
497
498 success = mpoly_monomial_cofactors(Abarexps, Bbarexps, A->exps, A->bits,
499 B->exps, B->bits, A->length, ctx->minfo);
500 if (!success)
501 goto cleanup_tmp;
502
503 fq_nmod_mpoly_init3(T, A->length, Gbits, ctx);
504 NG = mpoly_words_per_exp(Gbits, ctx->minfo);
505 NA = mpoly_words_per_exp(A->bits, ctx->minfo);
506 fq_nmod_inv(t1, A->coeffs + 0, ctx->fqctx);
507 T->length = A->length;
508 for (i = 0; i < A->length; i++)
509 {
510 mpoly_get_monomial_ffmpz(Texps, A->exps + NA*i, A->bits, ctx->minfo);
511 _fmpz_vec_sub(Texps, Texps, Abarexps, nvars);
512 mpoly_set_monomial_ffmpz(T->exps + NG*i, Texps, Gbits, ctx->minfo);
513 fq_nmod_mul(T->coeffs + i, A->coeffs + i, t1, ctx->fqctx);
514 }
515 fq_nmod_mpoly_swap(G, T, ctx);
516 fq_nmod_mpoly_clear(T, ctx);
517
518 success = 1;
519
520 cleanup_tmp:
521
522 for (j = 0; j < nvars; j++)
523 {
524 fmpz_clear(Abarexps + j);
525 fmpz_clear(Bbarexps + j);
526 fmpz_clear(Texps + j);
527 }
528
529 TMP_END;
530
531 cleanup:
532
533 fq_nmod_clear(t1, ctx->fqctx);
534 fq_nmod_clear(t2, ctx->fqctx);
535
536 return success;
537 }
538
539
540 /********* Assume B has length one when converted to univar format ***********/
_try_missing_var(fq_nmod_mpoly_t G,flint_bitcnt_t Gbits,slong var,const fq_nmod_mpoly_t A,ulong Ashift,const fq_nmod_mpoly_t B,ulong Bshift,const fq_nmod_mpoly_ctx_t ctx)541 static int _try_missing_var(
542 fq_nmod_mpoly_t G, flint_bitcnt_t Gbits,
543 slong var,
544 const fq_nmod_mpoly_t A, ulong Ashift,
545 const fq_nmod_mpoly_t B, ulong Bshift,
546 const fq_nmod_mpoly_ctx_t ctx)
547 {
548 int success;
549 slong i;
550 fq_nmod_mpoly_t tG;
551 fq_nmod_mpoly_univar_t Ax;
552
553 fq_nmod_mpoly_init(tG, ctx);
554 fq_nmod_mpoly_univar_init(Ax, ctx);
555
556 fq_nmod_mpoly_to_univar(Ax, A, var, ctx);
557
558 FLINT_ASSERT(Ax->length > 0);
559 success = _fq_nmod_mpoly_gcd(tG, Gbits, B, Ax->coeffs + 0, ctx);
560 if (!success)
561 goto cleanup;
562
563 for (i = 1; i < Ax->length; i++)
564 {
565 success = _fq_nmod_mpoly_gcd(tG, Gbits, tG, Ax->coeffs + i, ctx);
566 if (!success)
567 goto cleanup;
568 }
569
570 fq_nmod_mpoly_swap(G, tG, ctx);
571 _mpoly_gen_shift_left(G->exps, G->bits, G->length,
572 var, FLINT_MIN(Ashift, Bshift), ctx->minfo);
573
574 cleanup:
575
576 fq_nmod_mpoly_clear(tG, ctx);
577 fq_nmod_mpoly_univar_clear(Ax, ctx);
578
579 return success;
580 }
581
582
583 /******************* Test if B divides A or A divides B **********************/
584 /*
585 Test if B divides A or A divides B
586 TODO: incorporate deflation
587 */
_try_divides(fq_nmod_mpoly_t G,const fq_nmod_mpoly_t A,int try_a,const fq_nmod_mpoly_t B,int try_b,const fq_nmod_mpoly_ctx_t ctx)588 static int _try_divides(
589 fq_nmod_mpoly_t G,
590 const fq_nmod_mpoly_t A, int try_a,
591 const fq_nmod_mpoly_t B, int try_b,
592 const fq_nmod_mpoly_ctx_t ctx)
593 {
594 int success;
595 fq_nmod_mpoly_t Q;
596
597 fq_nmod_mpoly_init(Q, ctx);
598
599 if (try_b && fq_nmod_mpoly_divides(Q, A, B, ctx))
600 {
601 fq_nmod_mpoly_set(G, B, ctx);
602 success = 1;
603 goto cleanup;
604 }
605
606 if (try_a && fq_nmod_mpoly_divides(Q, B, A, ctx))
607 {
608 fq_nmod_mpoly_set(G, A, ctx);
609 success = 1;
610 goto cleanup;
611 }
612
613 success = 0;
614
615 cleanup:
616
617 fq_nmod_mpoly_clear(Q, ctx);
618
619 return success;
620 }
621
622
623 /********************** Hit A and B with zippel ******************************/
_try_zippel(fq_nmod_mpoly_t G,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const mpoly_gcd_info_t I,const fq_nmod_mpoly_ctx_t ctx)624 static int _try_zippel(
625 fq_nmod_mpoly_t G,
626 const fq_nmod_mpoly_t A,
627 const fq_nmod_mpoly_t B,
628 const mpoly_gcd_info_t I,
629 const fq_nmod_mpoly_ctx_t ctx)
630 {
631 slong i, k;
632 slong m = I->mvars;
633 int success;
634 mpoly_zipinfo_t zinfo;
635 flint_bitcnt_t wbits;
636 flint_rand_t randstate;
637 fq_nmod_mpoly_ctx_t uctx;
638 fq_nmod_mpolyu_t Au, Bu, Gu, Abaru, Bbaru;
639 fq_nmod_mpoly_t Ac, Bc, Gc;
640
641 FLINT_ASSERT(A->bits <= FLINT_BITS);
642 FLINT_ASSERT(B->bits <= FLINT_BITS);
643
644 if (!I->can_use_zippel)
645 return 0;
646
647 FLINT_ASSERT(m >= WORD(2));
648 FLINT_ASSERT(A->length > 0);
649 FLINT_ASSERT(B->length > 0);
650
651 flint_randinit(randstate);
652
653 /* interpolation will continue in m variables */
654 mpoly_zipinfo_init(zinfo, m);
655
656 /* uctx is context for Fq[y_1,...,y_{m-1}]*/
657 fq_nmod_mpoly_ctx_init(uctx, m - 1, ORD_LEX, ctx->fqctx);
658
659 /* fill in a valid zinfo->perm and degrees */
660 for (i = 0; i < m; i++)
661 {
662 k = I->zippel_perm[i];
663 zinfo->perm[i] = k;
664 zinfo->Adegs[i] = I->Adeflate_deg[k];
665 zinfo->Bdegs[i] = I->Bdeflate_deg[k];
666 FLINT_ASSERT(I->Adeflate_deg[k] != 0);
667 FLINT_ASSERT(I->Bdeflate_deg[k] != 0);
668 }
669
670 wbits = FLINT_MAX(A->bits, B->bits);
671
672 fq_nmod_mpolyu_init(Au, wbits, uctx);
673 fq_nmod_mpolyu_init(Bu, wbits, uctx);
674 fq_nmod_mpolyu_init(Gu, wbits, uctx);
675 fq_nmod_mpolyu_init(Abaru, wbits, uctx);
676 fq_nmod_mpolyu_init(Bbaru, wbits, uctx);
677 fq_nmod_mpoly_init3(Ac, 0, wbits, uctx);
678 fq_nmod_mpoly_init3(Bc, 0, wbits, uctx);
679 fq_nmod_mpoly_init3(Gc, 0, wbits, uctx);
680
681 fq_nmod_mpoly_to_mpolyu_perm_deflate(Au, uctx, A, ctx, zinfo->perm,
682 I->Amin_exp, I->Gstride);
683 fq_nmod_mpoly_to_mpolyu_perm_deflate(Bu, uctx, B, ctx, zinfo->perm,
684 I->Bmin_exp, I->Gstride);
685
686 success = fq_nmod_mpolyu_content_mpoly(Ac, Au, uctx);
687 success = success && fq_nmod_mpolyu_content_mpoly(Bc, Bu, uctx);
688 if (!success)
689 goto cleanup;
690
691 fq_nmod_mpolyu_divexact_mpoly_inplace(Au, Ac, uctx);
692 fq_nmod_mpolyu_divexact_mpoly_inplace(Bu, Bc, uctx);
693
694 /* after removing content, degree bounds in zinfo are still valid bounds */
695 success = fq_nmod_mpolyu_gcdm_zippel(Gu, Abaru, Bbaru, Au, Bu,
696 uctx, zinfo, randstate);
697 if (!success)
698 goto cleanup;
699
700 success = _fq_nmod_mpoly_gcd(Gc, wbits, Ac, Bc, uctx);
701 if (!success)
702 goto cleanup;
703
704 fq_nmod_mpolyu_mul_mpoly_inplace(Gu, Gc, uctx);
705
706 fq_nmod_mpoly_from_mpolyu_perm_inflate(G, I->Gbits, ctx, Gu, uctx,
707 zinfo->perm, I->Gmin_exp, I->Gstride);
708 success = 1;
709
710 cleanup:
711
712 fq_nmod_mpolyu_clear(Au, uctx);
713 fq_nmod_mpolyu_clear(Bu, uctx);
714 fq_nmod_mpolyu_clear(Gu, uctx);
715 fq_nmod_mpolyu_clear(Abaru, uctx);
716 fq_nmod_mpolyu_clear(Bbaru, uctx);
717 fq_nmod_mpoly_clear(Ac, uctx);
718 fq_nmod_mpoly_clear(Bc, uctx);
719 fq_nmod_mpoly_clear(Gc, uctx);
720
721 fq_nmod_mpoly_ctx_clear(uctx);
722
723 mpoly_zipinfo_clear(zinfo);
724
725 flint_randclear(randstate);
726
727 return success;
728 }
729
730
731 /*********************** Hit A and B with brown ******************************/
_try_brown(fq_nmod_mpoly_t G,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,mpoly_gcd_info_t I,const fq_nmod_mpoly_ctx_t ctx)732 static int _try_brown(
733 fq_nmod_mpoly_t G,
734 const fq_nmod_mpoly_t A,
735 const fq_nmod_mpoly_t B,
736 mpoly_gcd_info_t I,
737 const fq_nmod_mpoly_ctx_t ctx)
738 {
739 int success;
740 slong m = I->mvars;
741 flint_bitcnt_t wbits;
742 fq_nmod_mpoly_ctx_t nctx;
743 fq_nmod_mpolyn_t An, Bn, Gn, Abarn, Bbarn;
744
745 if (!I->can_use_brown)
746 return 0;
747
748 FLINT_ASSERT(m >= 2);
749 FLINT_ASSERT(A->bits <= FLINT_BITS);
750 FLINT_ASSERT(B->bits <= FLINT_BITS);
751 FLINT_ASSERT(A->length > 0);
752 FLINT_ASSERT(B->length > 0);
753
754 wbits = FLINT_MAX(A->bits, B->bits);
755
756 fq_nmod_mpoly_ctx_init(nctx, m, ORD_LEX, ctx->fqctx);
757 fq_nmod_mpolyn_init(An, wbits, nctx);
758 fq_nmod_mpolyn_init(Bn, wbits, nctx);
759 fq_nmod_mpolyn_init(Gn, wbits, nctx);
760 fq_nmod_mpolyn_init(Abarn, wbits, nctx);
761 fq_nmod_mpolyn_init(Bbarn, wbits, nctx);
762
763 fq_nmod_mpoly_to_mpolyn_perm_deflate(An, nctx, A, ctx,
764 I->brown_perm, I->Amin_exp, I->Gstride);
765 fq_nmod_mpoly_to_mpolyn_perm_deflate(Bn, nctx, B, ctx,
766 I->brown_perm, I->Bmin_exp, I->Gstride);
767
768 FLINT_ASSERT(An->bits == wbits);
769 FLINT_ASSERT(Bn->bits == wbits);
770 FLINT_ASSERT(An->length > 1);
771 FLINT_ASSERT(Bn->length > 1);
772
773 success = fq_nmod_mpolyn_gcd_brown_smprime(Gn, Abarn, Bbarn, An, Bn,
774 m - 1, nctx);
775 if (!success)
776 {
777 fq_nmod_mpoly_to_mpolyn_perm_deflate(An, nctx, A, ctx,
778 I->brown_perm, I->Amin_exp, I->Gstride);
779 fq_nmod_mpoly_to_mpolyn_perm_deflate(Bn, nctx, B, ctx,
780 I->brown_perm, I->Bmin_exp, I->Gstride);
781 success = fq_nmod_mpolyn_gcd_brown_lgprime(Gn, Abarn, Bbarn, An, Bn,
782 m - 1, nctx);
783 }
784
785 if (!success)
786 goto cleanup;
787
788 fq_nmod_mpoly_from_mpolyn_perm_inflate(G, I->Gbits, ctx, Gn, nctx,
789 I->brown_perm, I->Gmin_exp, I->Gstride);
790 success = 1;
791
792 cleanup:
793
794 fq_nmod_mpolyn_clear(An, nctx);
795 fq_nmod_mpolyn_clear(Bn, nctx);
796 fq_nmod_mpolyn_clear(Gn, nctx);
797 fq_nmod_mpolyn_clear(Abarn, nctx);
798 fq_nmod_mpolyn_clear(Bbarn, nctx);
799 fq_nmod_mpoly_ctx_clear(nctx);
800
801 return success;
802 }
803
804
805 /*
806 The function must pack its answer into bits = Gbits <= FLINT_BITS
807 Both A and B have to be packed into bits <= FLINT_BITS
808
809 return is 1 for success, 0 for failure.
810 */
_fq_nmod_mpoly_gcd(fq_nmod_mpoly_t G,flint_bitcnt_t Gbits,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)811 int _fq_nmod_mpoly_gcd(
812 fq_nmod_mpoly_t G, flint_bitcnt_t Gbits,
813 const fq_nmod_mpoly_t A,
814 const fq_nmod_mpoly_t B,
815 const fq_nmod_mpoly_ctx_t ctx)
816 {
817 int success;
818 slong v_in_both;
819 slong v_in_either;
820 slong v_in_A_only;
821 slong v_in_B_only;
822 slong j;
823 slong nvars = ctx->minfo->nvars;
824 mpoly_gcd_info_t I;
825
826 if (A->length == 1)
827 {
828 return _try_monomial_gcd(G, Gbits, B, A, ctx);
829 }
830 else if (B->length == 1)
831 {
832 return _try_monomial_gcd(G, Gbits, A, B, ctx);
833 }
834
835 mpoly_gcd_info_init(I, nvars);
836
837 /* entries of I are all now invalid */
838
839 I->Gbits = Gbits;
840
841 mpoly_gcd_info_limits(I->Amax_exp, I->Amin_exp, I->Alead_count,
842 I->Atail_count, A->exps, A->bits, A->length, ctx->minfo);
843 mpoly_gcd_info_limits(I->Bmax_exp, I->Bmin_exp, I->Blead_count,
844 I->Btail_count, B->exps, B->bits, B->length, ctx->minfo);
845
846 /* set ess(p) := p/term_content(p) */
847
848 /* check if the cofactors could be monomials, i.e. ess(A) == ess(B) */
849 for (j = 0; j < nvars; j++)
850 {
851 if (I->Amax_exp[j] - I->Amin_exp[j] != I->Bmax_exp[j] - I->Bmin_exp[j])
852 goto skip_monomial_cofactors;
853 }
854 if (_try_monomial_cofactors(G, I->Gbits, A, B, ctx))
855 {
856 goto successful;
857 }
858
859 skip_monomial_cofactors:
860
861 mpoly_gcd_info_stride(I->Gstride,
862 A->exps, A->bits, A->length, I->Amax_exp, I->Amin_exp,
863 B->exps, B->bits, B->length, I->Bmax_exp, I->Bmin_exp, ctx->minfo);
864
865 for (j = 0; j < nvars; j++)
866 {
867 ulong t = I->Gstride[j];
868
869 if (t == 0)
870 {
871 FLINT_ASSERT( I->Amax_exp[j] == I->Amin_exp[j]
872 || I->Bmax_exp[j] == I->Bmin_exp[j]);
873 }
874 else
875 {
876 FLINT_ASSERT((I->Amax_exp[j] - I->Amin_exp[j]) % t == 0);
877 FLINT_ASSERT((I->Bmax_exp[j] - I->Bmin_exp[j]) % t == 0);
878 }
879
880 I->Adeflate_deg[j] = t == 0 ? 0 : (I->Amax_exp[j] - I->Amin_exp[j])/t;
881 I->Bdeflate_deg[j] = t == 0 ? 0 : (I->Bmax_exp[j] - I->Bmin_exp[j])/t;
882 I->Gmin_exp[j] = FLINT_MIN(I->Amin_exp[j], I->Bmin_exp[j]);
883 }
884
885 /*
886 The following are now valid:
887 I->Amax_exp, I->Amin_exp, I->Alead_count, I->Atail_count,
888 I->Bmax_exp, I->Bmin_exp, I->Blead_count, I->Btail_count,
889 I->Gstride
890 I->Adeflate_deg
891 I->Bdeflate_deg
892 I->Gmin_exp
893 */
894
895 /* check if ess(A) and ess(B) have a variable v_in_both in common */
896 v_in_both = -WORD(1);
897 for (j = 0; j < nvars; j++)
898 {
899 if (I->Amax_exp[j] > I->Amin_exp[j] && I->Bmax_exp[j] > I->Bmin_exp[j])
900 {
901 v_in_both = j;
902 break;
903 }
904 }
905 if (v_in_both == -WORD(1))
906 {
907 /*
908 The variables in ess(A) and ess(B) are disjoint.
909 gcd is trivial to compute.
910 */
911
912 calculate_trivial_gcd:
913
914 fq_nmod_mpoly_fit_length(G, 1, ctx);
915 fq_nmod_mpoly_fit_bits(G, Gbits, ctx);
916 G->bits = Gbits;
917 mpoly_set_monomial_ui(G->exps, I->Gmin_exp, Gbits, ctx->minfo);
918 fq_nmod_one(G->coeffs + 0, ctx->fqctx);
919 _fq_nmod_mpoly_set_length(G, 1, ctx);
920
921 goto successful;
922 }
923
924 /* check if ess(A) and ess(B) depend on another variable v_in_either */
925 FLINT_ASSERT(0 <= v_in_both);
926 FLINT_ASSERT(v_in_both < nvars);
927
928 v_in_either = -WORD(1);
929 for (j = 0; j < nvars; j++)
930 {
931 if (j == v_in_both)
932 continue;
933
934 if (I->Amax_exp[j] > I->Amin_exp[j] || I->Bmax_exp[j] > I->Bmin_exp[j])
935 {
936 v_in_either = j;
937 break;
938 }
939 }
940
941 if (v_in_either == -WORD(1))
942 {
943 /*
944 The ess(A) and ess(B) depend on only one variable v_in_both
945 Calculate gcd using univariates
946 */
947 fq_nmod_poly_t a, b, g;
948
949 fq_nmod_poly_init(a, ctx->fqctx);
950 fq_nmod_poly_init(b, ctx->fqctx);
951 fq_nmod_poly_init(g, ctx->fqctx);
952 _fq_nmod_mpoly_to_fq_nmod_poly_deflate(a, A, v_in_both,
953 I->Amin_exp, I->Gstride, ctx);
954 _fq_nmod_mpoly_to_fq_nmod_poly_deflate(b, B, v_in_both,
955 I->Bmin_exp, I->Gstride, ctx);
956 fq_nmod_poly_gcd(g, a, b, ctx->fqctx);
957 _fq_nmod_mpoly_from_fq_nmod_poly_inflate(G, Gbits, g, v_in_both,
958 I->Gmin_exp, I->Gstride, ctx);
959 fq_nmod_poly_clear(a, ctx->fqctx);
960 fq_nmod_poly_clear(b, ctx->fqctx);
961 fq_nmod_poly_clear(g, ctx->fqctx);
962
963 goto successful;
964 }
965
966 /* check if there is a variable in ess(A) that is not in ess(B) */
967 v_in_A_only = -WORD(1);
968 v_in_B_only = -WORD(1);
969 for (j = 0; j < nvars; j++)
970 {
971 if (I->Amax_exp[j] > I->Amin_exp[j] && I->Bmax_exp[j] == I->Bmin_exp[j])
972 {
973 v_in_A_only = j;
974 break;
975 }
976 if (I->Bmax_exp[j] > I->Bmin_exp[j] && I->Amax_exp[j] == I->Amin_exp[j])
977 {
978 v_in_B_only = j;
979 break;
980 }
981 }
982 if (v_in_A_only != -WORD(1))
983 {
984 success = _try_missing_var(G, I->Gbits,
985 v_in_A_only,
986 A, I->Amin_exp[v_in_A_only],
987 B, I->Bmin_exp[v_in_A_only],
988 ctx);
989 goto cleanup;
990 }
991 if (v_in_B_only != -WORD(1))
992 {
993 success = _try_missing_var(G, I->Gbits,
994 v_in_B_only,
995 B, I->Bmin_exp[v_in_B_only],
996 A, I->Amin_exp[v_in_B_only],
997 ctx);
998 goto cleanup;
999 }
1000
1001 /* all variable are now either
1002 missing from both ess(A) and ess(B), or
1003 present in both ess(A) and ess(B)
1004 and there are at least two in the latter case
1005 */
1006
1007 mpoly_gcd_info_set_estimates_fq_nmod_mpoly(I, A, B, ctx);
1008 mpoly_gcd_info_set_perm(I, A->length, B->length, ctx->minfo);
1009
1010 /* everything in I is valid now */
1011
1012 /* check divisibility A/B and B/A */
1013 {
1014 int gcd_is_trivial = 1;
1015 int try_a = I->Gdeflate_deg_bounds_are_nice;
1016 int try_b = I->Gdeflate_deg_bounds_are_nice;
1017 for (j = 0; j < nvars; j++)
1018 {
1019 if (I->Gdeflate_deg_bound[j] != 0)
1020 {
1021 gcd_is_trivial = 0;
1022 }
1023
1024 if (I->Adeflate_deg[j] != I->Gdeflate_deg_bound[j]
1025 || I->Amin_exp[j] > I->Bmin_exp[j])
1026 {
1027 try_a = 0;
1028 }
1029
1030 if (I->Bdeflate_deg[j] != I->Gdeflate_deg_bound[j]
1031 || I->Bmin_exp[j] > I->Amin_exp[j])
1032 {
1033 try_b = 0;
1034 }
1035 }
1036
1037 if (gcd_is_trivial)
1038 goto calculate_trivial_gcd;
1039
1040 if ((try_a || try_b) && _try_divides(G, A, try_a, B, try_b, ctx))
1041 goto successful;
1042 }
1043
1044 mpoly_gcd_info_measure_brown(I, A->length, B->length, ctx->minfo);
1045 mpoly_gcd_info_measure_zippel(I, A->length, B->length, ctx->minfo);
1046
1047 if (I->zippel_time_est < I->brown_time_est)
1048 {
1049 if (_try_zippel(G, A, B, I, ctx))
1050 goto successful;
1051
1052 if (_try_brown(G, A, B, I, ctx))
1053 goto successful;
1054 }
1055 else
1056 {
1057 if (_try_brown(G, A, B, I, ctx))
1058 goto successful;
1059
1060 if (_try_zippel(G, A, B, I, ctx))
1061 goto successful;
1062 }
1063
1064 success = 0;
1065 goto cleanup;
1066
1067 successful:
1068
1069 success = 1;
1070
1071 cleanup:
1072
1073 mpoly_gcd_info_clear(I);
1074
1075 if (success)
1076 {
1077 fq_nmod_mpoly_repack_bits_inplace(G, Gbits, ctx);
1078 fq_nmod_mpoly_make_monic(G, G, ctx);
1079 }
1080
1081 return success;
1082 }
1083
fq_nmod_mpoly_gcd(fq_nmod_mpoly_t G,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)1084 int fq_nmod_mpoly_gcd(fq_nmod_mpoly_t G, const fq_nmod_mpoly_t A,
1085 const fq_nmod_mpoly_t B, const fq_nmod_mpoly_ctx_t ctx)
1086 {
1087 flint_bitcnt_t Gbits;
1088
1089 if (fq_nmod_mpoly_is_zero(A, ctx))
1090 {
1091 if (fq_nmod_mpoly_is_zero(B, ctx))
1092 fq_nmod_mpoly_zero(G, ctx);
1093 else
1094 fq_nmod_mpoly_make_monic(G, B, ctx);
1095 return 1;
1096 }
1097
1098 if (fq_nmod_mpoly_is_zero(B, ctx))
1099 {
1100 fq_nmod_mpoly_make_monic(G, A, ctx);
1101 return 1;
1102 }
1103
1104 Gbits = FLINT_MIN(A->bits, B->bits);
1105
1106 if (A->bits <= FLINT_BITS && B->bits <= FLINT_BITS)
1107 {
1108 /* usual gcd's go right down here */
1109 return _fq_nmod_mpoly_gcd(G, Gbits, A, B, ctx);
1110 }
1111
1112 if (A->length == 1)
1113 {
1114 return _try_monomial_gcd(G, Gbits, B, A, ctx);
1115 }
1116 else if (B->length == 1)
1117 {
1118 return _try_monomial_gcd(G, Gbits, A, B, ctx);
1119 }
1120 else if (_try_monomial_cofactors(G, Gbits, A, B, ctx))
1121 {
1122 return 1;
1123 }
1124 else
1125 {
1126 /*
1127 The gcd calculation is unusual.
1128 First see if both inputs fit into FLINT_BITS.
1129 Then, try deflation as a last resort.
1130 */
1131
1132 int success;
1133 int useAnew = 0;
1134 int useBnew = 0;
1135 slong k;
1136 fmpz * Ashift, * Astride;
1137 fmpz * Bshift, * Bstride;
1138 fmpz * Gshift, * Gstride;
1139 fq_nmod_mpoly_t Anew;
1140 fq_nmod_mpoly_t Bnew;
1141
1142 fq_nmod_mpoly_init(Anew, ctx);
1143 fq_nmod_mpoly_init(Bnew, ctx);
1144
1145 if (A->bits > FLINT_BITS)
1146 {
1147 useAnew = fq_nmod_mpoly_repack_bits(Anew, A, FLINT_BITS, ctx);
1148 if (!useAnew)
1149 goto could_not_repack;
1150 }
1151
1152 if (B->bits > FLINT_BITS)
1153 {
1154 useBnew = fq_nmod_mpoly_repack_bits(Bnew, B, FLINT_BITS, ctx);
1155 if (!useBnew)
1156 goto could_not_repack;
1157 }
1158
1159 success = _fq_nmod_mpoly_gcd(G, FLINT_BITS, useAnew ? Anew : A,
1160 useBnew ? Bnew : B, ctx);
1161 goto cleanup;
1162
1163 could_not_repack:
1164
1165 /*
1166 One of A or B could not be repacked into FLINT_BITS. See if
1167 they both fit into FLINT_BITS after deflation.
1168 */
1169
1170 Ashift = _fmpz_vec_init(ctx->minfo->nvars);
1171 Astride = _fmpz_vec_init(ctx->minfo->nvars);
1172 Bshift = _fmpz_vec_init(ctx->minfo->nvars);
1173 Bstride = _fmpz_vec_init(ctx->minfo->nvars);
1174 Gshift = _fmpz_vec_init(ctx->minfo->nvars);
1175 Gstride = _fmpz_vec_init(ctx->minfo->nvars);
1176
1177 fq_nmod_mpoly_deflation(Ashift, Astride, A, ctx);
1178 fq_nmod_mpoly_deflation(Bshift, Bstride, B, ctx);
1179 _fmpz_vec_min(Gshift, Ashift, Bshift, ctx->minfo->nvars);
1180 for (k = 0; k < ctx->minfo->nvars; k++)
1181 {
1182 fmpz_gcd(Gstride + k, Astride + k, Bstride + k);
1183 }
1184
1185 success = 0;
1186
1187 fq_nmod_mpoly_deflate(Anew, A, Ashift, Gstride, ctx);
1188 if (Anew->bits > FLINT_BITS)
1189 {
1190 if (!fq_nmod_mpoly_repack_bits(Anew, Anew, FLINT_BITS, ctx))
1191 goto deflate_cleanup;
1192 }
1193
1194 fq_nmod_mpoly_deflate(Bnew, B, Bshift, Gstride, ctx);
1195 if (Bnew->bits > FLINT_BITS)
1196 {
1197 if (!fq_nmod_mpoly_repack_bits(Bnew, Bnew, FLINT_BITS, ctx))
1198 goto deflate_cleanup;
1199 }
1200
1201 success = _fq_nmod_mpoly_gcd(G, FLINT_BITS, Anew, Bnew, ctx);
1202
1203 if (success)
1204 {
1205 fq_nmod_mpoly_inflate(G, G, Gshift, Gstride, ctx);
1206 fq_nmod_mpoly_make_monic(G, G, ctx);
1207 }
1208
1209 deflate_cleanup:
1210
1211 _fmpz_vec_clear(Ashift, ctx->minfo->nvars);
1212 _fmpz_vec_clear(Astride, ctx->minfo->nvars);
1213 _fmpz_vec_clear(Bshift, ctx->minfo->nvars);
1214 _fmpz_vec_clear(Bstride, ctx->minfo->nvars);
1215 _fmpz_vec_clear(Gshift, ctx->minfo->nvars);
1216 _fmpz_vec_clear(Gstride, ctx->minfo->nvars);
1217
1218 cleanup:
1219
1220 fq_nmod_mpoly_clear(Anew, ctx);
1221 fq_nmod_mpoly_clear(Bnew, ctx);
1222
1223 return success;
1224 }
1225 }
1226