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 "fq_nmod_mpoly.h"
13 #include "fq_nmod_mpoly_factor.h"
14 
15 
_fq_nmod_mpoly_monomial_evals_cache(n_fq_poly_t E,const ulong * Aexps,flint_bitcnt_t Abits,slong Alen,const fq_nmod_struct * betas,slong start,slong stop,const fq_nmod_mpoly_ctx_t ctx)16 void _fq_nmod_mpoly_monomial_evals_cache(
17     n_fq_poly_t E,
18     const ulong * Aexps,
19     flint_bitcnt_t Abits,
20     slong Alen,
21     const fq_nmod_struct * betas,
22     slong start,
23     slong stop,
24     const fq_nmod_mpoly_ctx_t ctx)
25 {
26     slong d = fq_nmod_ctx_degree(ctx->fqctx);
27     slong i, Ai;
28     ulong mask = (-UWORD(1)) >> (FLINT_BITS - Abits);
29     slong N = mpoly_words_per_exp_sp(Abits, ctx->minfo);
30     slong * off, * shift;
31     n_poly_struct * caches;
32     mp_limb_t * c;
33     slong num = stop - start;
34 
35     FLINT_ASSERT(Abits <= FLINT_BITS);
36     FLINT_ASSERT(Alen > 0);
37     FLINT_ASSERT(num > 0);
38 
39     caches = FLINT_ARRAY_ALLOC(3*num, n_fq_poly_struct);
40     off = FLINT_ARRAY_ALLOC(2*num, slong);
41     shift = off + num;
42     for (i = 0; i < num; i++)
43     {
44         mpoly_gen_offset_shift_sp(&off[i], &shift[i], i + start, Abits, ctx->minfo);
45         n_fq_poly_init(caches + 3*i + 0);
46         n_fq_poly_init(caches + 3*i + 1);
47         n_fq_poly_init(caches + 3*i + 2);
48         n_fq_pow_cache_start_fq_nmod(betas + i, caches + 3*i + 0,
49                                                 caches + 3*i + 1,
50                                                 caches + 3*i + 2, ctx->fqctx);
51     }
52 
53     n_poly_fit_length(E, d*Alen);
54     E->length = Alen;
55 
56     for (Ai = 0; Ai < Alen; Ai++)
57     {
58         c = E->coeffs + d*Ai;
59         _n_fq_one(c, d);
60         for (i = 0; i < num; i++)
61         {
62             ulong ei = (Aexps[N*Ai + off[i]] >> shift[i]) & mask;
63             n_fq_pow_cache_mulpow_ui(c, c, ei, caches + 3*i + 0,
64                                                caches + 3*i + 1,
65                                                caches + 3*i + 2, ctx->fqctx);
66         }
67     }
68 
69     for (i = 0; i < num; i++)
70     {
71         n_fq_poly_clear(caches + 3*i + 0);
72         n_fq_poly_clear(caches + 3*i + 1);
73         n_fq_poly_clear(caches + 3*i + 2);
74     }
75     flint_free(caches);
76     flint_free(off);
77 }
78 
79 
80 
_fq_nmod_mpoly_monomial_evals2_cache(n_fq_polyun_t E,const ulong * Aexps,flint_bitcnt_t Abits,slong Alen,const fq_nmod_struct * betas,slong m,const fq_nmod_mpoly_ctx_t ctx)81 void _fq_nmod_mpoly_monomial_evals2_cache(
82     n_fq_polyun_t E,
83     const ulong * Aexps,
84     flint_bitcnt_t Abits,
85     slong Alen,
86     const fq_nmod_struct * betas,
87     slong m,
88     const fq_nmod_mpoly_ctx_t ctx)
89 {
90     slong d = fq_nmod_ctx_degree(ctx->fqctx);
91     slong i, Ai, Ei;
92     ulong e0, e1, e01;
93     ulong mask = (-UWORD(1)) >> (FLINT_BITS - Abits);
94     slong N = mpoly_words_per_exp_sp(Abits, ctx->minfo);
95     slong * off, * shift;
96     n_fq_poly_struct * caches;
97     mp_limb_t * c;
98 
99     FLINT_ASSERT(Abits <= FLINT_BITS);
100     FLINT_ASSERT(Alen > 0);
101     FLINT_ASSERT(m > 2);
102 
103     caches = FLINT_ARRAY_ALLOC(3*(m - 2), n_fq_poly_struct);
104     off = FLINT_ARRAY_ALLOC(2*m, slong);
105     shift = off + m;
106     for (i = 0; i < m; i++)
107     {
108         mpoly_gen_offset_shift_sp(&off[i], &shift[i], i, Abits, ctx->minfo);
109         if (i >= 2)
110         {
111             n_fq_poly_init(caches + 3*(i - 2) + 0);
112             n_fq_poly_init(caches + 3*(i - 2) + 1);
113             n_fq_poly_init(caches + 3*(i - 2) + 2);
114             n_fq_pow_cache_start_fq_nmod(betas + i - 2, caches + 3*(i - 2) + 0,
115                                            caches + 3*(i - 2) + 1,
116                                            caches + 3*(i - 2) + 2, ctx->fqctx);
117         }
118     }
119 
120     Ai = 0;
121     Ei = 0;
122 
123     e0 = (Aexps[N*Ai + off[0]] >> shift[0]) & mask;
124     e1 = (Aexps[N*Ai + off[1]] >> shift[1]) & mask;
125     e01 = pack_exp2(e0, e1);
126     n_polyun_fit_length(E, Ei + 1);
127     E->exps[Ei] = e01;
128     n_poly_fit_length(E->coeffs + Ei, d*1);
129     c = E->coeffs[Ei].coeffs + d*0;
130     E->coeffs[Ei].length = 1;
131     Ei++;
132     _n_fq_one(c, d);
133     for (i = 2; i < m; i++)
134     {
135         ulong ei = (Aexps[N*Ai + off[i]] >> shift[i]) & mask;
136         n_fq_pow_cache_mulpow_ui(c, c, ei, caches + 3*(i - 2) + 0,
137                                            caches + 3*(i - 2) + 1,
138                                            caches + 3*(i - 2) + 2, ctx->fqctx);
139     }
140 
141     for (Ai++; Ai < Alen; Ai++)
142     {
143         e0 = (Aexps[N*Ai + off[0]] >> shift[0]) & mask;
144         e1 = (Aexps[N*Ai + off[1]] >> shift[1]) & mask;
145         e01 = pack_exp2(e0, e1);
146         if (e01 == E->exps[Ei - 1])
147         {
148             slong len = E->coeffs[Ei - 1].length;
149             n_poly_fit_length(E->coeffs + Ei - 1, d*(len + 1));
150             c = E->coeffs[Ei - 1].coeffs + d*len;
151             E->coeffs[Ei - 1].length = len + 1;
152         }
153         else
154         {
155             n_polyun_fit_length(E, Ei + 1);
156             E->exps[Ei] = e01;
157             n_poly_fit_length(E->coeffs + Ei, d*1);
158             c = E->coeffs[Ei].coeffs + d*0;
159             E->coeffs[Ei].length = 1;
160             Ei++;
161         }
162 
163         _n_fq_one(c, d);
164         for (i = 2; i < m; i++)
165         {
166             ulong ei = (Aexps[N*Ai + off[i]] >> shift[i]) & mask;
167             n_fq_pow_cache_mulpow_ui(c, c, ei, caches + 3*(i - 2) + 0,
168                                                caches + 3*(i - 2) + 1,
169                                                caches + 3*(i - 2) + 2, ctx->fqctx);
170         }
171     }
172 
173     E->length = Ei;
174 
175     for (i = 0; i < m - 2; i++)
176     {
177         n_fq_poly_clear(caches + 3*i + 0);
178         n_fq_poly_clear(caches + 3*i + 1);
179         n_fq_poly_clear(caches + 3*i + 2);
180     }
181     flint_free(caches);
182     flint_free(off);
183 
184 #if FLINT_WANT_ASSERT
185     Ai = 0;
186     for (i = 0; i < E->length; i++)
187         Ai += E->coeffs[i].length;
188     FLINT_ASSERT(Ai == Alen);
189 #endif
190 }
191 
192 
193 
n_fq_bpoly_eval_step_sep(n_fq_bpoly_t E,n_fq_polyun_t cur,const n_fq_polyun_t inc,const fq_nmod_mpoly_t A,const fq_nmod_ctx_t ctx)194 void n_fq_bpoly_eval_step_sep(
195     n_fq_bpoly_t E,
196     n_fq_polyun_t cur,
197     const n_fq_polyun_t inc,
198     const fq_nmod_mpoly_t A,
199     const fq_nmod_ctx_t ctx)
200 {
201     slong d = fq_nmod_ctx_degree(ctx);
202     slong i, Ai;
203     slong e0, e1;
204     mp_limb_t * c = FLINT_ARRAY_ALLOC(d, mp_limb_t);
205 
206     n_bpoly_zero(E);
207 
208     Ai = 0;
209     for (i = 0; i < cur->length; i++)
210     {
211         slong this_len = cur->coeffs[i].length;
212 
213         _n_fq_zip_eval_step(c, cur->coeffs[i].coeffs, inc->coeffs[i].coeffs,
214                                   A->coeffs + d*Ai, this_len, ctx);
215 
216         Ai += this_len;
217 
218         e0 = extract_exp(cur->exps[i], 1, 2);
219         e1 = extract_exp(cur->exps[i], 0, 2);
220         if (_n_fq_is_zero(c, d))
221             continue;
222 
223         n_fq_bpoly_set_coeff_n_fq(E, e0, e1, c, ctx);
224     }
225 
226     FLINT_ASSERT(Ai == A->length);
227 
228     flint_free(c);
229 }
230 
n_fq_poly_eval_step_sep(mp_limb_t * res,n_fq_poly_t cur,const n_fq_poly_t inc,const fq_nmod_mpoly_t A,const fq_nmod_ctx_t ctx)231 static void n_fq_poly_eval_step_sep(
232     mp_limb_t * res,
233     n_fq_poly_t cur,
234     const n_fq_poly_t inc,
235     const fq_nmod_mpoly_t A,
236     const fq_nmod_ctx_t ctx)
237 {
238     FLINT_ASSERT(A->length == cur->length);
239     FLINT_ASSERT(A->length == inc->length);
240     _n_fq_zip_eval_step(res, cur->coeffs, inc->coeffs, A->coeffs,
241                                                                A->length, ctx);
242 }
243 
n_fq_bpoly_evalp_step_sep(n_fq_bpoly_t E,n_polyun_t cur,const n_polyun_t inc,const fq_nmod_mpoly_t A,const fq_nmod_ctx_t ctx)244 static void n_fq_bpoly_evalp_step_sep(
245     n_fq_bpoly_t E,
246     n_polyun_t cur,
247     const n_polyun_t inc,
248     const fq_nmod_mpoly_t A,
249     const fq_nmod_ctx_t ctx)
250 {
251     slong d = fq_nmod_ctx_degree(ctx);
252     slong i, Ai;
253     slong e0, e1;
254     mp_limb_t * c = FLINT_ARRAY_ALLOC(d, mp_limb_t);
255 
256     n_bpoly_zero(E);
257 
258     Ai = 0;
259     for (i = 0; i < cur->length; i++)
260     {
261         slong this_len = cur->coeffs[i].length;
262 
263         _n_fqp_zip_eval_step(c, cur->coeffs[i].coeffs, inc->coeffs[i].coeffs,
264                                   A->coeffs + d*Ai, this_len, d, ctx->mod);
265 
266         Ai += this_len;
267 
268         e0 = extract_exp(cur->exps[i], 1, 2);
269         e1 = extract_exp(cur->exps[i], 0, 2);
270         if (_n_fq_is_zero(c, d))
271             continue;
272 
273         n_fq_bpoly_set_coeff_n_fq(E, e0, e1, c, ctx);
274     }
275 
276     FLINT_ASSERT(Ai == A->length);
277 
278     flint_free(c);
279 }
280 
281 
n_fq_poly_evalp_step_sep(mp_limb_t * res,n_poly_t cur,const n_poly_t inc,const fq_nmod_mpoly_t A,const fq_nmod_ctx_t ctx)282 static void n_fq_poly_evalp_step_sep(
283     mp_limb_t * res,
284     n_poly_t cur,
285     const n_poly_t inc,
286     const fq_nmod_mpoly_t A,
287     const fq_nmod_ctx_t ctx)
288 {
289     FLINT_ASSERT(A->length == cur->length);
290     FLINT_ASSERT(A->length == inc->length);
291     _n_fqp_zip_eval_step(res, cur->coeffs, inc->coeffs, A->coeffs,
292                                  A->length, fq_nmod_ctx_degree(ctx), ctx->mod);
293 }
294 
_fq_nmod_mpoly_bidegree(const fq_nmod_mpoly_t A,const fq_nmod_mpoly_ctx_t ctx)295 static ulong _fq_nmod_mpoly_bidegree(
296     const fq_nmod_mpoly_t A,
297     const fq_nmod_mpoly_ctx_t ctx)
298 {
299     FLINT_ASSERT(A->length > 0);
300     return _mpoly_bidegree(A->exps, A->bits, ctx->minfo);
301 }
302 
fq_nmod_mpoly_monomial_evals2(n_fq_polyun_t E,const fq_nmod_mpoly_t A,const fq_nmod_struct * betas,slong m,const fq_nmod_mpoly_ctx_t ctx)303 static void fq_nmod_mpoly_monomial_evals2(
304     n_fq_polyun_t E,
305     const fq_nmod_mpoly_t A,
306     const fq_nmod_struct * betas,
307     slong m,
308     const fq_nmod_mpoly_ctx_t ctx)
309 {
310     _fq_nmod_mpoly_monomial_evals2_cache(E, A->exps, A->bits, A->length, betas, m, ctx);
311 }
312 
fq_nmod_mpoly_monomial_evalsp2(n_polyun_t E,const fq_nmod_mpoly_t A,const mp_limb_t * betas,slong m,const fq_nmod_mpoly_ctx_t ctx)313 static void fq_nmod_mpoly_monomial_evalsp2(
314     n_polyun_t E,
315     const fq_nmod_mpoly_t A,
316     const mp_limb_t * betas,
317     slong m,
318     const fq_nmod_mpoly_ctx_t ctx)
319 {
320     _nmod_mpoly_monomial_evals2_cache(E, A->exps, A->bits, A->length, betas, m,
321                                                   ctx->minfo, ctx->fqctx->mod);
322 }
323 
fq_nmod_mpoly_monomial_evalsp(n_poly_t E,const fq_nmod_mpoly_t A,const mp_limb_t * betas,slong start,slong stop,const fq_nmod_mpoly_ctx_t ctx)324 static void fq_nmod_mpoly_monomial_evalsp(
325     n_poly_t E,
326     const fq_nmod_mpoly_t A,
327     const mp_limb_t * betas,
328     slong start,
329     slong stop,
330     const fq_nmod_mpoly_ctx_t ctx)
331 {
332     _nmod_mpoly_monomial_evals_cache(E, A->exps, A->bits, A->length, betas,
333                                      start, stop, ctx->minfo, ctx->fqctx->mod);
334 }
335 
fq_nmod_mpoly_monomial_evals(n_fq_poly_t E,const fq_nmod_mpoly_t A,const fq_nmod_struct * betas,slong start,slong stop,const fq_nmod_mpoly_ctx_t ctx)336 static void fq_nmod_mpoly_monomial_evals(
337     n_fq_poly_t E,
338     const fq_nmod_mpoly_t A,
339     const fq_nmod_struct * betas,
340     slong start,
341     slong stop,
342     const fq_nmod_mpoly_ctx_t ctx)
343 {
344     _fq_nmod_mpoly_monomial_evals_cache(E, A->exps, A->bits, A->length,
345                                                       betas, start, stop, ctx);
346 }
347 
348 
n_fq_polyun_zip_solve(fq_nmod_mpoly_t A,n_fq_polyun_t Z,n_fq_polyun_t H,n_fq_polyun_t M,const fq_nmod_mpoly_ctx_t ctx)349 int n_fq_polyun_zip_solve(
350     fq_nmod_mpoly_t A,
351     n_fq_polyun_t Z,
352     n_fq_polyun_t H,
353     n_fq_polyun_t M,
354     const fq_nmod_mpoly_ctx_t ctx)
355 {
356     slong d = fq_nmod_ctx_degree(ctx->fqctx);
357     int success;
358     slong Ai, i, n;
359     n_poly_t t;
360 
361     n_poly_init(t);
362 
363     FLINT_ASSERT(Z->length == H->length);
364     FLINT_ASSERT(Z->length == M->length);
365 
366     /* A has the right length, but possibly from a smaller ctx */
367     if (A->length*d > A->coeffs_alloc)
368     {
369         slong new_alloc = FLINT_MAX(A->coeffs_alloc + A->coeffs_alloc/2, A->length*d);
370         A->coeffs = (mp_limb_t *) flint_realloc(A->coeffs, new_alloc*sizeof(mp_limb_t));
371         A->coeffs_alloc = new_alloc;
372     }
373 
374     Ai = 0;
375     for (i = 0; i < H->length; i++)
376     {
377         n = H->coeffs[i].length;
378         FLINT_ASSERT(M->coeffs[i].length == n + 1);
379         FLINT_ASSERT(Z->coeffs[i].length >= n);
380         FLINT_ASSERT(Ai + n <= A->length);
381 
382         n_poly_fit_length(t, d*n);
383 
384         success = _n_fq_zip_vand_solve(A->coeffs + d*Ai,
385                          H->coeffs[i].coeffs, n,
386                          Z->coeffs[i].coeffs, Z->coeffs[i].length,
387                          M->coeffs[i].coeffs, t->coeffs, ctx->fqctx);
388         if (success < 1)
389         {
390             n_poly_clear(t);
391             return success;
392         }
393 
394         Ai += n;
395         FLINT_ASSERT(Ai <= A->length);
396     }
397 
398     FLINT_ASSERT(Ai == A->length);
399 
400     n_poly_clear(t);
401     return 1;
402 }
403 
404 
405 
n_fq_polyun_zip_solvep(fq_nmod_mpoly_t A,n_fq_polyun_t Z,n_fq_polyun_t H,n_fq_polyun_t M,const fq_nmod_mpoly_ctx_t ctx)406 static int n_fq_polyun_zip_solvep(
407     fq_nmod_mpoly_t A,
408     n_fq_polyun_t Z,
409     n_fq_polyun_t H,
410     n_fq_polyun_t M,
411     const fq_nmod_mpoly_ctx_t ctx)
412 {
413     slong d = fq_nmod_ctx_degree(ctx->fqctx);
414     int success;
415     slong Ai, i, n;
416     n_poly_t t;
417 
418     n_poly_init(t);
419 
420     FLINT_ASSERT(Z->length == H->length);
421     FLINT_ASSERT(Z->length == M->length);
422 
423     /* A has the right length, but possibly from a smaller ctx */
424     if (A->length*d > A->coeffs_alloc)
425     {
426         slong new_alloc = FLINT_MAX(A->coeffs_alloc + A->coeffs_alloc/2, A->length*d);
427         A->coeffs = FLINT_ARRAY_REALLOC(A->coeffs, new_alloc, mp_limb_t);
428         A->coeffs_alloc = new_alloc;
429     }
430 
431     Ai = 0;
432     for (i = 0; i < H->length; i++)
433     {
434         n = H->coeffs[i].length;
435         FLINT_ASSERT(M->coeffs[i].length == n + 1);
436         FLINT_ASSERT(Z->coeffs[i].length >= n);
437         FLINT_ASSERT(Ai + n <= A->length);
438 
439         n_poly_fit_length(t, n);
440 
441         success = _n_fqp_zip_vand_solve(A->coeffs + d*Ai,
442                          H->coeffs[i].coeffs, n,
443                          Z->coeffs[i].coeffs, Z->coeffs[i].length,
444                          M->coeffs[i].coeffs, t->coeffs, ctx->fqctx);
445         if (success < 1)
446         {
447             n_poly_clear(t);
448             return success;
449         }
450 
451         Ai += n;
452         FLINT_ASSERT(Ai <= A->length);
453     }
454 
455     FLINT_ASSERT(Ai == A->length);
456 
457     n_poly_clear(t);
458     return 1;
459 }
460 
461 
n_fq_polyun_zip_start(n_fq_polyun_t Z,n_fq_polyun_t H,slong req_images,const fq_nmod_ctx_t ctx)462 void n_fq_polyun_zip_start(
463     n_fq_polyun_t Z,
464     n_fq_polyun_t H,
465     slong req_images,
466     const fq_nmod_ctx_t ctx)
467 {
468     slong d = fq_nmod_ctx_degree(ctx);
469     slong j;
470     n_polyun_fit_length(Z, H->length);
471     Z->length = H->length;
472     for (j = 0; j < H->length; j++)
473     {
474         Z->exps[j] = H->exps[j];
475         n_poly_fit_length(Z->coeffs + j, d*req_images);
476         Z->coeffs[j].length = 0;
477     }
478 }
479 
480 
n_fq_polyu2n_add_zip_must_match(n_fq_polyun_t Z,const n_fq_bpoly_t A,slong cur_length,const fq_nmod_ctx_t ctx)481 int n_fq_polyu2n_add_zip_must_match(
482     n_fq_polyun_t Z,
483     const n_fq_bpoly_t A,
484     slong cur_length,
485     const fq_nmod_ctx_t ctx)
486 {
487     slong d = fq_nmod_ctx_degree(ctx);
488     slong i, Ai, ai;
489     n_poly_struct * Zcoeffs = Z->coeffs;
490     ulong * Zexps = Z->exps;
491     const n_poly_struct * Acoeffs = A->coeffs;
492 
493     Ai = A->length - 1;
494     ai = (Ai < 0) ? 0 : n_fq_poly_degree(A->coeffs + Ai);
495 
496     for (i = 0; i < Z->length; i++)
497     {
498         if (Ai >= 0 && Zexps[i] == pack_exp2(Ai, ai))
499         {
500             /* Z present, A present */
501             _n_fq_set(Zcoeffs[i].coeffs + d*cur_length, Acoeffs[Ai].coeffs + d*ai, d);
502             Zcoeffs[i].length = cur_length + 1;
503             do {
504                 ai--;
505             } while (ai >= 0 && _n_fq_is_zero(Acoeffs[Ai].coeffs + d*ai, d));
506             if (ai < 0)
507             {
508                 do {
509                     Ai--;
510                 } while (Ai >= 0 && Acoeffs[Ai].length == 0);
511                 if (Ai >= 0)
512                     ai = n_fq_poly_degree(Acoeffs + Ai);
513             }
514         }
515         else if (Ai < 0 || Zexps[i] > pack_exp2(Ai, ai))
516         {
517             /* Z present, A missing */
518             _n_fq_zero(Zcoeffs[i].coeffs + d*cur_length, d);
519             Zcoeffs[i].length = cur_length + 1;
520         }
521         else
522         {
523             /* Z missing, A present */
524             return 0;
525         }
526     }
527 
528     return 1;
529 }
530 
531 
532 
533 #define USE_G    1
534 #define USE_ABAR 2
535 #define USE_BBAR 4
536 
537 /*
538     gamma = gcd(lc(A), lc(B))
539 
540     fake answers G, Abar, Bbar with
541 
542         G*Abar = gamma*A
543         G*Bbar = gamma*B
544         lc(G) = gamma
545         lc(Abar) = lc(A)
546         lc(Bbar) = lc(B)
547 
548     real answers
549 
550         rG = pp(G)
551         rAbar = Abar/lc(rG)
552         rBbar = Bbar/lc(rG)
553 */
fq_nmod_mpolyl_gcd_zippel_smprime(fq_nmod_mpoly_t rG,const slong * rGdegs,fq_nmod_mpoly_t rAbar,fq_nmod_mpoly_t rBbar,const fq_nmod_mpoly_t A,const slong * Adegs,const fq_nmod_mpoly_t B,const slong * Bdegs,const fq_nmod_mpoly_t gamma,const slong * gammadegs,const fq_nmod_mpoly_ctx_t ctx)554 int fq_nmod_mpolyl_gcd_zippel_smprime(
555     fq_nmod_mpoly_t rG, const slong * rGdegs,
556     fq_nmod_mpoly_t rAbar,
557     fq_nmod_mpoly_t rBbar,
558     const fq_nmod_mpoly_t A, const slong * Adegs,
559     const fq_nmod_mpoly_t B, const slong * Bdegs,
560     const fq_nmod_mpoly_t gamma, const slong * gammadegs,
561     const fq_nmod_mpoly_ctx_t ctx)
562 {
563     slong d = fq_nmod_ctx_degree(ctx->fqctx);
564     int success, use, betas_in_fp, main_tries_left;
565     slong i, m;
566     slong nvars = ctx->minfo->nvars;
567     flint_bitcnt_t bits = A->bits;
568     fq_nmod_struct * alphas, * betas;
569     mp_limb_t * betasp;
570     flint_rand_t state;
571     fq_nmod_mpoly_t cont;
572     fq_nmod_mpoly_t T, G, Abar, Bbar;
573     n_fq_polyun_t HG, HAbar, HBbar, MG, MAbar, MBbar, ZG, ZAbar, ZBbar;
574     n_fq_bpoly_t Aev, Bev, Gev, Abarev, Bbarev;
575     const mp_limb_t * gammaev;
576     fq_nmod_mpolyn_t Tn, Gn, Abarn, Bbarn;
577     slong lastdeg;
578     slong cur_zip_image, req_zip_images, this_length;
579     n_fq_polyun_t Aeh_cur, Aeh_inc, Beh_cur, Beh_inc;
580     n_fq_poly_t gammaeh_cur, gammaeh_inc;
581     n_fq_poly_t modulus, alphapow;
582     fq_nmod_mpoly_struct * Aevals, * Bevals;
583     fq_nmod_mpoly_struct * gammaevals;
584     n_poly_bpoly_stack_t St;
585     mp_limb_t * c;
586     fq_nmod_t start_alpha;
587     ulong GdegboundXY, newdegXY, Abideg, Bbideg;
588     slong degxAB, degyAB;
589 
590     FLINT_ASSERT(bits <= FLINT_BITS);
591     FLINT_ASSERT(bits == A->bits);
592     FLINT_ASSERT(bits == B->bits);
593     FLINT_ASSERT(bits == gamma->bits);
594     FLINT_ASSERT(A->length > 0);
595     FLINT_ASSERT(B->length > 0);
596     FLINT_ASSERT(gamma->length > 0);
597 
598     fq_nmod_mpoly_fit_length_reset_bits(rG, 1, bits, ctx);
599     fq_nmod_mpoly_fit_length_reset_bits(rAbar, 1, bits, ctx);
600     fq_nmod_mpoly_fit_length_reset_bits(rBbar, 1, bits, ctx);
601 
602 #if FLINT_WANT_ASSERT
603     {
604         slong * tmp_degs = FLINT_ARRAY_ALLOC(nvars, slong);
605 
606         fq_nmod_mpoly_degrees_si(tmp_degs, A, ctx);
607         for (i = 0; i < nvars; i++)
608             FLINT_ASSERT(tmp_degs[i] == Adegs[i]);
609 
610         fq_nmod_mpoly_degrees_si(tmp_degs, B, ctx);
611         for (i = 0; i < nvars; i++)
612             FLINT_ASSERT(tmp_degs[i] == Bdegs[i]);
613 
614         fq_nmod_mpoly_degrees_si(tmp_degs, gamma, ctx);
615         for (i = 0; i < nvars; i++)
616             FLINT_ASSERT(tmp_degs[i] == gammadegs[i]);
617 
618         flint_free(tmp_degs);
619     }
620 #endif
621 
622     FLINT_ASSERT(gammadegs[0] == 0);
623     FLINT_ASSERT(gammadegs[1] == 0);
624 
625     fq_nmod_init(start_alpha, ctx->fqctx);
626     c = FLINT_ARRAY_ALLOC(d, mp_limb_t);
627 
628     n_fq_polyun_init(HG);
629     n_fq_polyun_init(HAbar);
630     n_fq_polyun_init(HBbar);
631     n_fq_polyun_init(MG);
632     n_fq_polyun_init(MAbar);
633     n_fq_polyun_init(MBbar);
634     n_fq_polyun_init(ZG);
635     n_fq_polyun_init(ZAbar);
636     n_fq_polyun_init(ZBbar);
637     n_fq_bpoly_init(Aev);
638     n_fq_bpoly_init(Bev);
639     n_fq_bpoly_init(Gev);
640     n_fq_bpoly_init(Abarev);
641     n_fq_bpoly_init(Bbarev);
642     n_fq_poly_init2(alphapow, 4, ctx->fqctx);
643     fq_nmod_mpoly_init3(cont, 1, bits, ctx);
644     fq_nmod_mpoly_init3(T, 0, bits, ctx);
645     fq_nmod_mpoly_init3(G, 0, bits, ctx);
646     fq_nmod_mpoly_init3(Abar, 0, bits, ctx);
647     fq_nmod_mpoly_init3(Bbar, 0, bits, ctx);
648     fq_nmod_mpolyn_init(Tn, bits, ctx);
649     fq_nmod_mpolyn_init(Gn, bits, ctx);
650     fq_nmod_mpolyn_init(Abarn, bits, ctx);
651     fq_nmod_mpolyn_init(Bbarn, bits, ctx);
652     n_fq_polyun_init(Aeh_cur);
653     n_fq_polyun_init(Aeh_inc);
654     n_fq_polyun_init(Beh_cur);
655     n_fq_polyun_init(Beh_inc);
656     n_fq_poly_init(gammaeh_cur);
657     n_fq_poly_init(gammaeh_inc);
658     n_fq_poly_init(modulus);
659     n_poly_stack_init(St->poly_stack);
660     n_bpoly_stack_init(St->bpoly_stack);
661 
662     betasp = FLINT_ARRAY_ALLOC(nvars, mp_limb_t);
663     betas = FLINT_ARRAY_ALLOC(nvars, fq_nmod_struct);
664     alphas = FLINT_ARRAY_ALLOC(nvars, fq_nmod_struct);
665     for (i = 0; i < nvars; i++)
666     {
667         fq_nmod_init(betas + i, ctx->fqctx);
668         fq_nmod_init(alphas + i, ctx->fqctx);
669     }
670 
671     flint_randinit(state);
672 
673     Aevals = FLINT_ARRAY_ALLOC(nvars + 1, fq_nmod_mpoly_struct);
674     Bevals = FLINT_ARRAY_ALLOC(nvars + 1, fq_nmod_mpoly_struct);
675     gammaevals = FLINT_ARRAY_ALLOC(nvars + 1, fq_nmod_mpoly_struct);
676     for (i = 0; i < nvars; i++)
677     {
678         fq_nmod_mpoly_init3(Aevals + i, 0, bits, ctx);
679         fq_nmod_mpoly_init3(Bevals + i, 0, bits, ctx);
680         fq_nmod_mpoly_init3(gammaevals + i, 0, bits, ctx);
681     }
682     Aevals[nvars] = *A;
683     Bevals[nvars] = *B;
684     gammaevals[nvars] = *gamma;
685 
686     Abideg = _fq_nmod_mpoly_bidegree(A, ctx);
687     Bbideg = _fq_nmod_mpoly_bidegree(B, ctx);
688 
689     degxAB = FLINT_MAX(Adegs[0], Bdegs[0]);
690     degyAB = FLINT_MAX(Adegs[1], Bdegs[1]);
691 
692     GdegboundXY = pack_exp2(FLINT_MIN(Adegs[0], Bdegs[0]),
693                             FLINT_MIN(Adegs[1], Bdegs[1]));
694     if (GdegboundXY == 0)
695         goto gcd_is_trivial;
696 
697     main_tries_left = 3;
698 
699 choose_main:
700 
701     if (--main_tries_left < 0)
702     {
703         success = 0;
704         goto cleanup;
705     }
706 
707     for (i = 2; i < nvars; i++)
708         fq_nmod_rand_not_zero(alphas + i, state, ctx->fqctx);
709 
710     for (i = nvars - 1; i >= 2; i--)
711     {
712         fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + i, Aevals + i + 1, i, alphas + i, ctx);
713         fq_nmod_mpoly_repack_bits_inplace(Aevals + i, bits, ctx);
714         fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + i, Bevals + i + 1, i, alphas + i, ctx);
715         fq_nmod_mpoly_repack_bits_inplace(Bevals + i, bits, ctx);
716         fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + i, gammaevals + i + 1, i, alphas + i, ctx);
717         fq_nmod_mpoly_repack_bits_inplace(gammaevals + i, bits, ctx);
718         if (fq_nmod_mpoly_is_zero(gammaevals + i, ctx))
719             goto choose_main;
720         if (Aevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + i, ctx) != Abideg)
721             goto choose_main;
722         if (Bevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + i, ctx) != Bbideg)
723             goto choose_main;
724     }
725 
726     m = 2;
727     if (rGdegs == NULL)
728         use = USE_G | USE_ABAR | USE_BBAR;
729     else
730         use = mpoly_gcd_get_use_first(rGdegs[m], Adegs[m], Bdegs[m], gammadegs[m]);
731 
732     fq_nmod_mpoly_get_n_fq_bpoly(Aev, Aevals + m, 0, 1, ctx);
733     fq_nmod_mpoly_get_n_fq_bpoly(Bev, Bevals + m, 0, 1, ctx);
734 
735     success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
736                                                      Aev, Bev, ctx->fqctx, St);
737     if (!success)
738         goto cleanup;
739 
740     newdegXY = n_bpoly_bidegree(Gev);
741     if (newdegXY > GdegboundXY)
742         goto choose_main;
743     if (newdegXY < GdegboundXY)
744     {
745         GdegboundXY = newdegXY;
746         if (GdegboundXY == 0)
747             goto gcd_is_trivial;
748     }
749 
750     gammaev = fq_nmod_mpoly_get_nonzero_n_fq(gammaevals + m, ctx);
751     n_fq_bpoly_scalar_mul_n_fq(Gev, gammaev, ctx->fqctx);
752 
753     fq_nmod_mpolyn_interp_lift_sm_bpoly(Gn, Gev, ctx);
754     fq_nmod_mpolyn_interp_lift_sm_bpoly(Abarn, Abarev, ctx);
755     fq_nmod_mpolyn_interp_lift_sm_bpoly(Bbarn, Bbarev, ctx);
756 
757     n_fq_poly_one(modulus, ctx->fqctx);
758     n_fq_set_fq_nmod(c, alphas + m, ctx->fqctx);
759     n_fq_poly_shift_left_scalar_submul(modulus, 1, c, ctx->fqctx);
760 
761     fq_nmod_set(start_alpha, alphas + m, ctx->fqctx);
762     while (1)
763     {
764     choose_alpha_2:
765 
766         fq_nmod_next_not_zero(alphas + m, ctx->fqctx);
767 
768         if (fq_nmod_equal(alphas + m, start_alpha, ctx->fqctx))
769             goto choose_main;
770 
771         FLINT_ASSERT(alphapow->alloc >= d*2);
772         alphapow->length = 2;
773         _n_fq_one(alphapow->coeffs + d*0, d);
774         n_fq_set_fq_nmod(alphapow->coeffs + d*1, alphas + m, ctx->fqctx);
775 
776         fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + m, Aevals + m + 1, m, alphas + m, ctx);
777         fq_nmod_mpoly_repack_bits_inplace(Aevals + m, bits, ctx);
778         fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + m, Bevals + m + 1, m, alphas + m, ctx);
779         fq_nmod_mpoly_repack_bits_inplace(Bevals + m, bits, ctx);
780         fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + m, gammaevals + m + 1, m, alphas + m, ctx);
781         fq_nmod_mpoly_repack_bits_inplace(gammaevals + m, bits, ctx);
782         if (fq_nmod_mpoly_is_zero(gammaevals + m, ctx))
783             goto choose_alpha_2;
784         if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, ctx) != Abideg)
785             goto choose_alpha_2;
786         if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, ctx) != Bbideg)
787             goto choose_alpha_2;
788 
789         fq_nmod_mpoly_get_n_fq_bpoly(Aev, Aevals + m, 0, 1, ctx);
790         fq_nmod_mpoly_get_n_fq_bpoly(Bev, Bevals + m, 0, 1, ctx);
791 
792         success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
793                                                      Aev, Bev, ctx->fqctx, St);
794         if (!success)
795             goto cleanup;
796 
797         newdegXY = n_bpoly_bidegree(Gev);
798         if (newdegXY > GdegboundXY)
799             goto choose_alpha_2;
800         if (newdegXY < GdegboundXY)
801         {
802             GdegboundXY = newdegXY;
803             if (GdegboundXY == 0)
804                 goto gcd_is_trivial;
805             goto choose_main;
806         }
807 
808         gammaev = fq_nmod_mpoly_get_nonzero_n_fq(gammaevals + m, ctx);
809         n_fq_bpoly_scalar_mul_n_fq(Gev, gammaev, ctx->fqctx);
810 
811         n_fq_poly_eval_pow(c, modulus, alphapow, ctx->fqctx);
812         n_fq_inv(c, c, ctx->fqctx);
813         n_fq_poly_scalar_mul_n_fq(modulus, modulus, c, ctx->fqctx);
814 
815         if ((use & USE_G) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
816                                 &lastdeg, Gn, Tn, Gev, modulus, alphapow, ctx))
817         {
818             if (m == nvars - 1)
819             {
820                 fq_nmod_mpoly_cvtfrom_mpolyn(rG, Gn, m, ctx);
821                 success = fq_nmod_mpolyl_content(cont, rG, 2, ctx);
822                 if (!success)
823                     goto cleanup;
824                 fq_nmod_mpoly_divides(rG, rG, cont, ctx);
825                 fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
826                 if (fq_nmod_mpoly_divides(rAbar, A, rG, ctx) &&
827                     fq_nmod_mpoly_divides(rBbar, B, rG, ctx))
828                 {
829                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
830                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
831                     break;
832                 }
833             }
834             else
835             {
836                 fq_nmod_mpoly_cvtfrom_mpolyn(G, Gn, m, ctx);
837                 fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
838                 if (fq_nmod_mpoly_divides(Abar, T, G, ctx))
839                 {
840                     fq_nmod_mpoly_repack_bits_inplace(Abar, bits, ctx);
841                     fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
842                     if (fq_nmod_mpoly_divides(Bbar, T, G, ctx))
843                     {
844                         fq_nmod_mpoly_repack_bits_inplace(Bbar, bits, ctx);
845                         break;
846                     }
847                 }
848             }
849         }
850 
851         if ((use & USE_ABAR) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
852                           &lastdeg, Abarn, Tn, Abarev, modulus, alphapow, ctx))
853         {
854             if (m == nvars - 1)
855             {
856                 fq_nmod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, ctx);
857                 success = fq_nmod_mpolyl_content(cont, rAbar, 2, ctx);
858                 if (!success)
859                     goto cleanup;
860                 fq_nmod_mpoly_divides(rAbar, rAbar, cont, ctx);
861                 fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
862                 if (fq_nmod_mpoly_divides(rG, A, rAbar, ctx) &&
863                     fq_nmod_mpoly_divides(rBbar, B, rG, ctx))
864                 {
865                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
866                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
867                     break;
868                 }
869             }
870             else
871             {
872                 fq_nmod_mpoly_cvtfrom_mpolyn(Abar, Abarn, m, ctx);
873                 fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
874                 if (fq_nmod_mpoly_divides(G, T, Abar, ctx))
875                 {
876                     fq_nmod_mpoly_repack_bits_inplace(G, bits, ctx);
877                     fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
878                     if (fq_nmod_mpoly_divides(Bbar, T, G, ctx))
879                     {
880                         fq_nmod_mpoly_repack_bits_inplace(Bbar, bits, ctx);
881                         break;
882                     }
883                 }
884             }
885         }
886 
887         if ((use & USE_BBAR) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
888                           &lastdeg, Bbarn, Tn, Bbarev, modulus, alphapow, ctx))
889         {
890             if (m == nvars - 1)
891             {
892                 fq_nmod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, ctx);
893                 success = fq_nmod_mpolyl_content(cont, rBbar, 2, ctx);
894                 if (!success)
895                     goto cleanup;
896                 fq_nmod_mpoly_divides(rBbar, rBbar, cont, ctx);
897                 fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
898                 if (fq_nmod_mpoly_divides(rG, B, rBbar, ctx) &&
899                     fq_nmod_mpoly_divides(rAbar, A, rG, ctx))
900                 {
901                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
902                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
903                     break;
904                 }
905             }
906             else
907             {
908                 fq_nmod_mpoly_cvtfrom_mpolyn(Bbar, Bbarn, m, ctx);
909                 fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
910                 if (fq_nmod_mpoly_divides(G, T, Bbar, ctx))
911                 {
912                     fq_nmod_mpoly_repack_bits_inplace(G, bits, ctx);
913                     fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
914                     if (fq_nmod_mpoly_divides(Abar, T, G, ctx))
915                     {
916                         fq_nmod_mpoly_repack_bits_inplace(Abar, bits, ctx);
917                         break;
918                     }
919                 }
920             }
921         }
922 
923         if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
924             n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
925         {
926             goto choose_main;
927         }
928 
929         FLINT_ASSERT(n_fq_equal_fq_nmod(alphapow->coeffs + d*1, alphas + m, ctx->fqctx));
930         n_fq_poly_shift_left_scalar_submul(modulus, 1, alphapow->coeffs + d*1, ctx->fqctx);
931     }
932 
933     for (m = 3; m < nvars; m++)
934     {
935         /* G, Abar, Bbar are in Fq[gen(0), ..., gen(m - 1)] */
936         fq_nmod_mpolyn_interp_lift_sm_mpoly(Gn, G, ctx);
937         fq_nmod_mpolyn_interp_lift_sm_mpoly(Abarn, Abar, ctx);
938         fq_nmod_mpolyn_interp_lift_sm_mpoly(Bbarn, Bbar, ctx);
939 
940         if (rGdegs == NULL)
941             use = USE_G | USE_ABAR | USE_BBAR;
942         else
943             use = 0;
944 
945         n_fq_poly_one(modulus, ctx->fqctx);
946         n_fq_set_fq_nmod(c, alphas + m, ctx->fqctx);
947         n_fq_poly_shift_left_scalar_submul(modulus, 1, c, ctx->fqctx);
948 
949         fq_nmod_set(start_alpha, alphas + m, ctx->fqctx);
950 
951     choose_betas:
952 
953         /* only beta[0], beta[1], ..., beta[m - 1] will be used */
954         betas_in_fp = (ctx->fqctx->mod.norm < FLINT_BITS/4);
955         if (betas_in_fp)
956         {
957             for (i = 2; i < ctx->minfo->nvars; i++)
958                 betasp[i] = n_urandint(state, ctx->fqctx->mod.n - 3) + 2;
959 
960             fq_nmod_mpoly_monomial_evalsp2(HG, G, betasp + 2, m, ctx);
961             fq_nmod_mpoly_monomial_evalsp2(HAbar, Abar, betasp + 2, m, ctx);
962             fq_nmod_mpoly_monomial_evalsp2(HBbar, Bbar, betasp + 2, m, ctx);
963         }
964         else
965         {
966             for (i = 2; i < ctx->minfo->nvars; i++)
967                 fq_nmod_rand_not_zero(betas + i, state, ctx->fqctx);
968 
969             FLINT_ASSERT(G->bits == bits);
970             FLINT_ASSERT(Abar->bits == bits);
971             FLINT_ASSERT(Bbar->bits == bits);
972 
973             fq_nmod_mpoly_monomial_evals2(HG, G, betas + 2, m, ctx);
974             fq_nmod_mpoly_monomial_evals2(HAbar, Abar, betas + 2, m, ctx);
975             fq_nmod_mpoly_monomial_evals2(HBbar, Bbar, betas + 2, m, ctx);
976         }
977 
978         if (use == 0)
979         {
980             this_length = gammaevals[m + 1].length + Aevals[m + 1].length +
981                                                      Bevals[m + 1].length;
982 
983             use = nmod_mpoly_gcd_get_use_new(rGdegs[m], Adegs[m], Bdegs[m],
984                   gammadegs[m], degxAB, degyAB, this_length, HG, HAbar, HBbar);
985         }
986 
987         req_zip_images = 1;
988         if (use & USE_G)
989         {
990             this_length = betas_in_fp ?
991                          n_polyun_product_roots(MG, HG, ctx->fqctx->mod) :
992                          n_fq_polyun_product_roots(MG, HG, ctx->fqctx, St->poly_stack);
993             req_zip_images = FLINT_MAX(req_zip_images, this_length);
994         }
995         if (use & USE_ABAR)
996         {
997             this_length = betas_in_fp ?
998                          n_polyun_product_roots(MAbar, HAbar, ctx->fqctx->mod) :
999                          n_fq_polyun_product_roots(MAbar, HAbar, ctx->fqctx, St->poly_stack);
1000             req_zip_images = FLINT_MAX(req_zip_images, this_length);
1001         }
1002         if (use & USE_BBAR)
1003         {
1004             this_length = betas_in_fp ?
1005                          n_polyun_product_roots(MBbar, HBbar, ctx->fqctx->mod) :
1006                          n_fq_polyun_product_roots(MBbar, HBbar, ctx->fqctx, St->poly_stack);
1007             req_zip_images = FLINT_MAX(req_zip_images, this_length);
1008         }
1009 
1010         while (1)
1011         {
1012         choose_alpha_m:
1013 
1014             fq_nmod_next_not_zero(alphas + m, ctx->fqctx);
1015 
1016             if (fq_nmod_equal(alphas + m, start_alpha, ctx->fqctx))
1017                 goto choose_main;
1018 
1019             FLINT_ASSERT(alphapow->alloc >= d*2);
1020             alphapow->length = 2;
1021             _n_fq_one(alphapow->coeffs + d*0, d);
1022             n_fq_set_fq_nmod(alphapow->coeffs + d*1, alphas + m, ctx->fqctx);
1023 
1024             fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + m, Aevals + m + 1, m, alphas + m, ctx);
1025             fq_nmod_mpoly_repack_bits_inplace(Aevals + m, bits, ctx);
1026             fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + m, Bevals + m + 1, m, alphas + m, ctx);
1027             fq_nmod_mpoly_repack_bits_inplace(Bevals + m, bits, ctx);
1028             fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + m, gammaevals + m + 1, m, alphas + m, ctx);
1029             fq_nmod_mpoly_repack_bits_inplace(gammaevals + m, bits, ctx);
1030             if (fq_nmod_mpoly_is_zero(gammaevals + m, ctx))
1031                 goto choose_alpha_m;
1032             if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, ctx) != Abideg)
1033                 goto choose_alpha_m;
1034             if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, ctx) != Bbideg)
1035                 goto choose_alpha_m;
1036 
1037             if (betas_in_fp)
1038             {
1039                 fq_nmod_mpoly_monomial_evalsp2(Aeh_inc, Aevals + m, betasp + 2, m, ctx);
1040                 fq_nmod_mpoly_monomial_evalsp2(Beh_inc, Bevals + m, betasp + 2, m, ctx);
1041                 fq_nmod_mpoly_monomial_evalsp(gammaeh_inc, gammaevals + m, betasp + 2, 2, m, ctx);
1042                 n_polyun_set(Aeh_cur, Aeh_inc);
1043                 n_polyun_set(Beh_cur, Beh_inc);
1044                 n_poly_set(gammaeh_cur, gammaeh_inc);
1045             }
1046             else
1047             {
1048                 fq_nmod_mpoly_monomial_evals2(Aeh_inc, Aevals + m, betas + 2, m, ctx);
1049                 fq_nmod_mpoly_monomial_evals2(Beh_inc, Bevals + m, betas + 2, m, ctx);
1050                 fq_nmod_mpoly_monomial_evals(gammaeh_inc, gammaevals + m, betas + 2, 2, m, ctx);
1051                 n_fq_polyun_set(Aeh_cur, Aeh_inc, ctx->fqctx);
1052                 n_fq_polyun_set(Beh_cur, Beh_inc, ctx->fqctx);
1053                 n_fq_poly_set(gammaeh_cur, gammaeh_inc, ctx->fqctx);
1054             }
1055 
1056             n_fq_polyun_zip_start(ZG, HG, req_zip_images, ctx->fqctx);
1057             n_fq_polyun_zip_start(ZAbar, HAbar, req_zip_images, ctx->fqctx);
1058             n_fq_polyun_zip_start(ZBbar, HBbar, req_zip_images, ctx->fqctx);
1059             for (cur_zip_image = 0; cur_zip_image < req_zip_images; cur_zip_image++)
1060             {
1061                 if (betas_in_fp)
1062                 {
1063                     n_fq_bpoly_evalp_step_sep(Aev, Aeh_cur, Aeh_inc, Aevals + m, ctx->fqctx);
1064                     n_fq_bpoly_evalp_step_sep(Bev, Beh_cur, Beh_inc, Bevals + m, ctx->fqctx);
1065                     n_fq_poly_evalp_step_sep(c, gammaeh_cur, gammaeh_inc, gammaevals + m, ctx->fqctx);
1066                 }
1067                 else
1068                 {
1069                     n_fq_bpoly_eval_step_sep(Aev, Aeh_cur, Aeh_inc, Aevals + m, ctx->fqctx);
1070                     n_fq_bpoly_eval_step_sep(Bev, Beh_cur, Beh_inc, Bevals + m, ctx->fqctx);
1071                     n_fq_poly_eval_step_sep(c, gammaeh_cur, gammaeh_inc, gammaevals + m, ctx->fqctx);
1072                 }
1073 
1074                 if (_n_fq_is_zero(c, d))
1075                     goto choose_betas;
1076                 if (Aev->length < 1 || n_bpoly_bidegree(Aev) != Abideg)
1077                     goto choose_betas;
1078                 if (Bev->length < 1 || n_bpoly_bidegree(Bev) != Bbideg)
1079                     goto choose_betas;
1080 
1081                 success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
1082                                                      Aev, Bev, ctx->fqctx, St);
1083                 if (!success)
1084                     goto cleanup;
1085 
1086                 newdegXY = n_bpoly_bidegree(Gev);
1087                 if (newdegXY > GdegboundXY)
1088                     goto choose_betas;
1089                 if (newdegXY < GdegboundXY)
1090                 {
1091                     GdegboundXY = newdegXY;
1092                     if (GdegboundXY == 0)
1093                         goto gcd_is_trivial;
1094                     goto choose_main;
1095                 }
1096 
1097                 n_fq_bpoly_scalar_mul_n_fq(Gev, c, ctx->fqctx);
1098                 if ((use & USE_G) && !n_fq_polyu2n_add_zip_must_match(ZG, Gev, cur_zip_image, ctx->fqctx))
1099                     goto choose_main;
1100                 if ((use & USE_ABAR) && !n_fq_polyu2n_add_zip_must_match(ZAbar, Abarev, cur_zip_image, ctx->fqctx))
1101                     goto choose_main;
1102                 if ((use & USE_BBAR) && !n_fq_polyu2n_add_zip_must_match(ZBbar, Bbarev, cur_zip_image, ctx->fqctx))
1103                     goto choose_main;
1104             }
1105 
1106             if (betas_in_fp)
1107             {
1108                 if ((use & USE_G) && n_fq_polyun_zip_solvep(G, ZG, HG, MG, ctx) < 1)
1109                     goto choose_main;
1110                 if ((use & USE_ABAR) && n_fq_polyun_zip_solvep(Abar, ZAbar, HAbar, MAbar, ctx) < 1)
1111                     goto choose_main;
1112                 if ((use & USE_BBAR) && n_fq_polyun_zip_solvep(Bbar, ZBbar, HBbar, MBbar, ctx) < 1)
1113                     goto choose_main;
1114             }
1115             else
1116             {
1117                 if ((use & USE_G) && n_fq_polyun_zip_solve(G, ZG, HG, MG, ctx) < 1)
1118                     goto choose_main;
1119                 if ((use & USE_ABAR) && n_fq_polyun_zip_solve(Abar, ZAbar, HAbar, MAbar, ctx) < 1)
1120                     goto choose_main;
1121                 if ((use & USE_BBAR) && n_fq_polyun_zip_solve(Bbar, ZBbar, HBbar, MBbar, ctx) < 1)
1122                     goto choose_main;
1123             }
1124 
1125             n_fq_poly_eval_pow(c, modulus, alphapow, ctx->fqctx);
1126             n_fq_inv(c, c, ctx->fqctx);
1127             n_fq_poly_scalar_mul_n_fq(modulus, modulus, c, ctx->fqctx);
1128 
1129             if ((use & USE_G) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
1130                                       &lastdeg, Gn, G, modulus, alphapow, ctx))
1131             {
1132                 fq_nmod_mpoly_cvtfrom_mpolyn(rG, Gn, m, ctx);
1133                 if (m == nvars - 1)
1134                 {
1135                     success = fq_nmod_mpolyl_content(cont, rG, 2, ctx);
1136                     if (!success)
1137                         goto cleanup;
1138                     fq_nmod_mpoly_divides(rG, rG, cont, ctx);
1139                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
1140                     if (fq_nmod_mpoly_divides(rAbar, A, rG, ctx) &&
1141                         fq_nmod_mpoly_divides(rBbar, B, rG,  ctx))
1142                     {
1143                         fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1144                         fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1145                         break;
1146                     }
1147                 }
1148                 else
1149                 {
1150                     fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1151                     if (fq_nmod_mpoly_divides(rAbar, T, rG, ctx))
1152                     {
1153                         fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1154                         if (fq_nmod_mpoly_divides(rBbar, T, rG, ctx))
1155                         {
1156                             fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1157                             fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1158                             fq_nmod_mpoly_swap(G, rG, ctx);
1159                             fq_nmod_mpoly_swap(Abar, rAbar, ctx);
1160                             fq_nmod_mpoly_swap(Bbar, rBbar, ctx);
1161                             break;
1162                         }
1163                     }
1164                 }
1165             }
1166             if ((use & USE_ABAR) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
1167                                 &lastdeg, Abarn, Abar, modulus, alphapow, ctx))
1168             {
1169                 fq_nmod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, ctx);
1170                 if (m == nvars - 1)
1171                 {
1172                     success = fq_nmod_mpolyl_content(cont, rAbar, 2, ctx);
1173                     if (!success)
1174                         goto cleanup;
1175                     success = fq_nmod_mpoly_divides(rAbar, rAbar, cont, ctx);
1176                     FLINT_ASSERT(success);
1177                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1178                     if (fq_nmod_mpoly_divides(rG, A, rAbar, ctx) &&
1179                         fq_nmod_mpoly_divides(rBbar, B, rG, ctx))
1180                     {
1181                         fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
1182                         fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1183                         break;
1184                     }
1185                 }
1186                 else
1187                 {
1188                     fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1189                     if (fq_nmod_mpoly_divides(rG, T, rAbar, ctx))
1190                     {
1191                         fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1192                         if (fq_nmod_mpoly_divides(rBbar, T, rG, ctx))
1193                         {
1194                             fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
1195                             fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1196                             fq_nmod_mpoly_swap(G, rG, ctx);
1197                             fq_nmod_mpoly_swap(Abar, rAbar, ctx);
1198                             fq_nmod_mpoly_swap(Bbar, rBbar, ctx);
1199                             break;
1200                         }
1201                     }
1202                 }
1203             }
1204             if ((use & USE_BBAR) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
1205                                 &lastdeg, Bbarn, Bbar, modulus, alphapow, ctx))
1206             {
1207                 fq_nmod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, ctx);
1208                 if (m == nvars - 1)
1209                 {
1210                     success = fq_nmod_mpolyl_content(cont, rBbar, 2, ctx);
1211                     if (!success)
1212                         goto cleanup;
1213                     fq_nmod_mpoly_divides(rBbar, rBbar, cont, ctx);
1214                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1215                     if (fq_nmod_mpoly_divides(rG, B, rBbar, ctx) &&
1216                         fq_nmod_mpoly_divides(rAbar, A, rG, ctx))
1217                     {
1218                         fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
1219                         fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1220                         break;
1221                     }
1222                 }
1223                 else
1224                 {
1225                     fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1226                     if (fq_nmod_mpoly_divides(rG, T, rBbar, ctx))
1227                     {
1228                         fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1229                         if (fq_nmod_mpoly_divides(rAbar, T, rG, ctx))
1230                         {
1231                             fq_nmod_mpoly_repack_bits_inplace(rG, bits, ctx);
1232                             fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1233                             fq_nmod_mpoly_swap(G, rG, ctx);
1234                             fq_nmod_mpoly_swap(Abar, rAbar, ctx);
1235                             fq_nmod_mpoly_swap(Bbar, rBbar, ctx);
1236                             break;
1237                         }
1238                     }
1239                 }
1240             }
1241 
1242             if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
1243                 n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
1244             {
1245                 goto choose_main;
1246             }
1247 
1248             FLINT_ASSERT(n_fq_equal_fq_nmod(alphapow->coeffs + d*1, alphas + m, ctx->fqctx));
1249             n_fq_poly_shift_left_scalar_submul(modulus, 1, alphapow->coeffs + d*1, ctx->fqctx);
1250         }
1251     }
1252 
1253     success = 1;
1254 
1255 cleanup:
1256 
1257     n_fq_polyun_clear(HG);
1258     n_fq_polyun_clear(HAbar);
1259     n_fq_polyun_clear(HBbar);
1260     n_fq_polyun_clear(MG);
1261     n_fq_polyun_clear(MAbar);
1262     n_fq_polyun_clear(MBbar);
1263     n_fq_polyun_clear(ZG);
1264     n_fq_polyun_clear(ZAbar);
1265     n_fq_polyun_clear(ZBbar);
1266     n_fq_bpoly_clear(Aev);
1267     n_fq_bpoly_clear(Bev);
1268     n_fq_bpoly_clear(Gev);
1269     n_fq_bpoly_clear(Abarev);
1270     n_fq_bpoly_clear(Bbarev);
1271     n_fq_poly_clear(alphapow);
1272     fq_nmod_mpoly_clear(cont, ctx);
1273     fq_nmod_mpoly_clear(T, ctx);
1274     fq_nmod_mpoly_clear(G, ctx);
1275     fq_nmod_mpoly_clear(Abar, ctx);
1276     fq_nmod_mpoly_clear(Bbar, ctx);
1277     fq_nmod_mpolyn_clear(Tn, ctx);
1278     fq_nmod_mpolyn_clear(Gn, ctx);
1279     fq_nmod_mpolyn_clear(Abarn, ctx);
1280     fq_nmod_mpolyn_clear(Bbarn, ctx);
1281     n_fq_polyun_clear(Aeh_cur);
1282     n_fq_polyun_clear(Aeh_inc);
1283     n_fq_polyun_clear(Beh_cur);
1284     n_fq_polyun_clear(Beh_inc);
1285     n_fq_poly_clear(gammaeh_cur);
1286     n_fq_poly_clear(gammaeh_inc);
1287     n_fq_poly_clear(modulus);
1288     n_poly_stack_clear(St->poly_stack);
1289     n_bpoly_stack_clear(St->bpoly_stack);
1290 
1291     fq_nmod_clear(start_alpha, ctx->fqctx);
1292     flint_free(c);
1293 
1294     flint_free(betasp);
1295 
1296     for (i = 0; i < nvars; i++)
1297     {
1298         fq_nmod_clear(betas + i, ctx->fqctx);
1299         fq_nmod_clear(alphas + i, ctx->fqctx);
1300     }
1301     flint_free(betas);
1302     flint_free(alphas);
1303     flint_randclear(state);
1304 
1305     for (i = 0; i < nvars; i++)
1306     {
1307         fq_nmod_mpoly_clear(Aevals + i, ctx);
1308         fq_nmod_mpoly_clear(Bevals + i, ctx);
1309         fq_nmod_mpoly_clear(gammaevals + i, ctx);
1310     }
1311     flint_free(Aevals);
1312     flint_free(Bevals);
1313     flint_free(gammaevals);
1314 
1315     FLINT_ASSERT(!success || rG->bits == bits);
1316     FLINT_ASSERT(!success || rAbar->bits == bits);
1317     FLINT_ASSERT(!success || rBbar->bits == bits);
1318 
1319     return success;
1320 
1321 gcd_is_trivial:
1322 
1323     fq_nmod_mpoly_one(rG, ctx);
1324     fq_nmod_mpoly_set(rAbar, A, ctx);
1325     fq_nmod_mpoly_set(rBbar, B, ctx);
1326 
1327     success = 1;
1328 
1329     goto cleanup;
1330 }
1331 
1332 
1333 /* TODO: make this the real fq_nmod_mpolyn_interp_mcrt_lg_mpoly */
newfq_nmod_mpolyn_interp_mcrt_lg_mpoly(slong * lastdeg,fq_nmod_mpolyn_t H,const fq_nmod_mpoly_ctx_t ctx,const n_fq_poly_t m,const mp_limb_t * inv_m_eval,fq_nmod_mpoly_t A,const fq_nmod_mpoly_ctx_t ectx,const bad_fq_nmod_embed_t emb)1334 static int newfq_nmod_mpolyn_interp_mcrt_lg_mpoly(
1335     slong * lastdeg,
1336     fq_nmod_mpolyn_t H,
1337     const fq_nmod_mpoly_ctx_t ctx,
1338     const n_fq_poly_t m,
1339     const mp_limb_t * inv_m_eval,
1340     fq_nmod_mpoly_t A,
1341     const fq_nmod_mpoly_ctx_t ectx,
1342     const bad_fq_nmod_embed_t emb)
1343 {
1344     slong lgd = fq_nmod_ctx_degree(ectx->fqctx);
1345     slong i;
1346 #if FLINT_WANT_ASSERT
1347     slong N = mpoly_words_per_exp(A->bits, ctx->minfo);
1348 #endif
1349     int changed = 0;
1350     mp_limb_t * u, * v, * tmp;
1351     n_fq_poly_struct * w, * u_sm;
1352     n_poly_stack_t St;
1353 
1354     FLINT_ASSERT(H->length == A->length);
1355     FLINT_ASSERT(H->bits == A->bits);
1356 
1357     n_poly_stack_init(St);
1358     n_poly_stack_fit_request(St, 3);
1359 
1360     w = n_poly_stack_take_top(St);
1361     u_sm = n_poly_stack_take_top(St);
1362 
1363     tmp = n_poly_stack_vec_init(St, lgd*(2 + N_FQ_MUL_ITCH));
1364     u = tmp + lgd*N_FQ_MUL_ITCH;
1365     v = u + lgd;
1366 
1367     for (i = 0; i < A->length; i++)
1368     {
1369         FLINT_ASSERT(mpoly_monomial_equal(H->exps + N*i, A->exps + N*i, N));
1370         bad_n_fq_embed_sm_to_lg(u, H->coeffs + i, emb);
1371         _n_fq_sub(v, A->coeffs + lgd*i, u, lgd, ectx->fqctx->mod);
1372         if (!_n_fq_is_zero(v, lgd))
1373         {
1374             changed = 1;
1375             _n_fq_mul(u, v, inv_m_eval, ectx->fqctx, tmp);
1376             bad_n_fq_embed_lg_to_sm(u_sm, u, emb);
1377             n_fq_poly_mul_(w, u_sm, m, ctx->fqctx, St);
1378             n_fq_poly_add(H->coeffs + i, H->coeffs + i, w, ctx->fqctx);
1379         }
1380 
1381         *lastdeg = FLINT_MAX(*lastdeg, n_fq_poly_degree(H->coeffs + i));
1382 
1383         FLINT_ASSERT(n_fq_poly_degree(H->coeffs + i) < n_fq_poly_degree(m) +
1384                                       fq_nmod_poly_degree(emb->h, ctx->fqctx));
1385     }
1386 
1387     n_poly_stack_vec_clear(St);
1388 
1389     n_poly_stack_give_back(St, 2);
1390 
1391     n_poly_stack_clear(St);
1392 
1393     return changed;
1394 }
1395 
1396 
1397 
fq_nmod_mpolyl_gcd_zippel_lgprime(fq_nmod_mpoly_t rG,const slong * rGdegs,fq_nmod_mpoly_t rAbar,fq_nmod_mpoly_t rBbar,const fq_nmod_mpoly_t A,const slong * Adegs,const fq_nmod_mpoly_t B,const slong * Bdegs,const fq_nmod_mpoly_t gamma,const slong * gammadegs,const fq_nmod_mpoly_ctx_t smctx)1398 int fq_nmod_mpolyl_gcd_zippel_lgprime(
1399     fq_nmod_mpoly_t rG, const slong * rGdegs,
1400     fq_nmod_mpoly_t rAbar,
1401     fq_nmod_mpoly_t rBbar,
1402     const fq_nmod_mpoly_t A, const slong * Adegs,
1403     const fq_nmod_mpoly_t B, const slong * Bdegs,
1404     const fq_nmod_mpoly_t gamma, const slong * gammadegs,
1405     const fq_nmod_mpoly_ctx_t smctx)
1406 {
1407     slong lgd;
1408     int success, use;
1409     slong i, m;
1410     slong nvars = smctx->minfo->nvars;
1411     flint_bitcnt_t bits = A->bits;
1412     fq_nmod_struct * alphas, * betas;
1413     flint_rand_t state;
1414     fq_nmod_mpoly_t cont;
1415     fq_nmod_mpoly_t T, G, Abar, Bbar;
1416     n_fq_polyun_t HG, HAbar, HBbar, MG, MAbar, MBbar, ZG, ZAbar, ZBbar;
1417     n_fq_bpoly_t Aev, Bev, Gev, Abarev, Bbarev;
1418     const mp_limb_t * gammaev;
1419     fq_nmod_mpolyn_t Tn, Gn, Abarn, Bbarn;
1420     slong lastdeg;
1421     slong cur_zip_image, req_zip_images, this_length;
1422     n_fq_polyun_t Aeh_cur, Aeh_inc, Beh_cur, Beh_inc;
1423     n_fq_poly_t gammaeh_cur, gammaeh_inc;
1424     fq_nmod_mpoly_struct * Aevals, * Bevals;
1425     fq_nmod_mpoly_struct * gammaevals;
1426     n_fq_poly_t modulus, alphapow;
1427     n_poly_bpoly_stack_t St;
1428     fq_nmod_t start_alpha;
1429     n_poly_t tmp;  /* tmp arithmetic space */
1430     ulong GdegboundXY, newdegXY, Abideg, Bbideg;
1431     slong degxAB, degyAB;
1432     bad_fq_nmod_mpoly_embed_chooser_t embc;
1433     bad_fq_nmod_embed_struct * cur_emb;
1434     fq_nmod_mpoly_ctx_t lgctx;
1435     fq_nmod_mpolyn_t gamman;
1436     fq_nmod_mpolyn_t An, Bn;
1437 
1438     FLINT_ASSERT(bits <= FLINT_BITS);
1439     FLINT_ASSERT(bits == A->bits);
1440     FLINT_ASSERT(bits == B->bits);
1441     FLINT_ASSERT(bits == gamma->bits);
1442     FLINT_ASSERT(A->length > 0);
1443     FLINT_ASSERT(B->length > 0);
1444     FLINT_ASSERT(gamma->length > 0);
1445 
1446     fq_nmod_mpoly_fit_length_reset_bits(rG, 1, bits, smctx);
1447     fq_nmod_mpoly_fit_length_reset_bits(rAbar, 1, bits, smctx);
1448     fq_nmod_mpoly_fit_length_reset_bits(rBbar, 1, bits, smctx);
1449 
1450 #if FLINT_WANT_ASSERT
1451     {
1452         slong * tmp_degs = FLINT_ARRAY_ALLOC(nvars, slong);
1453 
1454         fq_nmod_mpoly_degrees_si(tmp_degs, A, smctx);
1455         for (i = 0; i < nvars; i++)
1456             FLINT_ASSERT(tmp_degs[i] == Adegs[i]);
1457 
1458         fq_nmod_mpoly_degrees_si(tmp_degs, B, smctx);
1459         for (i = 0; i < nvars; i++)
1460             FLINT_ASSERT(tmp_degs[i] == Bdegs[i]);
1461 
1462         fq_nmod_mpoly_degrees_si(tmp_degs, gamma, smctx);
1463         for (i = 0; i < nvars; i++)
1464             FLINT_ASSERT(tmp_degs[i] == gammadegs[i]);
1465 
1466         flint_free(tmp_degs);
1467     }
1468 #endif
1469 
1470     FLINT_ASSERT(gammadegs[0] == 0);
1471     FLINT_ASSERT(gammadegs[1] == 0);
1472 
1473     n_fq_polyun_init(HG);
1474     n_fq_polyun_init(HAbar);
1475     n_fq_polyun_init(HBbar);
1476     n_fq_polyun_init(MG);
1477     n_fq_polyun_init(MAbar);
1478     n_fq_polyun_init(MBbar);
1479     n_fq_polyun_init(ZG);
1480     n_fq_polyun_init(ZAbar);
1481     n_fq_polyun_init(ZBbar);
1482     n_fq_bpoly_init(Aev);
1483     n_fq_bpoly_init(Bev);
1484     n_fq_bpoly_init(Gev);
1485     n_fq_bpoly_init(Abarev);
1486     n_fq_bpoly_init(Bbarev);
1487     fq_nmod_mpoly_init3(cont, 1, bits, smctx);
1488     fq_nmod_mpoly_init3(T, 0, bits, smctx);
1489     fq_nmod_mpoly_init3(G, 0, bits, smctx);
1490     fq_nmod_mpoly_init3(Abar, 0, bits, smctx);
1491     fq_nmod_mpoly_init3(Bbar, 0, bits, smctx);
1492     fq_nmod_mpolyn_init(Tn, bits, smctx);
1493     fq_nmod_mpolyn_init(Gn, bits, smctx);
1494     fq_nmod_mpolyn_init(Abarn, bits, smctx);
1495     fq_nmod_mpolyn_init(Bbarn, bits, smctx);
1496     n_fq_polyun_init(Aeh_cur);
1497     n_fq_polyun_init(Aeh_inc);
1498     n_fq_polyun_init(Beh_cur);
1499     n_fq_polyun_init(Beh_inc);
1500     n_fq_poly_init(gammaeh_cur);
1501     n_fq_poly_init(gammaeh_inc);
1502     n_fq_poly_init(modulus);
1503     n_poly_stack_init(St->poly_stack);
1504     n_bpoly_stack_init(St->bpoly_stack);
1505     fq_nmod_mpolyn_init(An, bits, smctx);
1506     fq_nmod_mpolyn_init(Bn, bits, smctx);
1507     fq_nmod_mpolyn_init(gamman, bits, smctx);
1508 
1509     /* alphas[nvars - 1] not used - it is replaced by cur_emb */
1510     betas = FLINT_ARRAY_ALLOC(nvars, fq_nmod_struct);
1511     alphas = FLINT_ARRAY_ALLOC(nvars, fq_nmod_struct);
1512     for (i = 0; i < nvars; i++)
1513     {
1514         fq_nmod_init(betas + i, smctx->fqctx);
1515         fq_nmod_init(alphas + i, smctx->fqctx);
1516     }
1517 
1518     flint_randinit(state);
1519     cur_emb = bad_fq_nmod_mpoly_embed_chooser_init(embc, lgctx, smctx, state);
1520     lgd = fq_nmod_ctx_degree(lgctx->fqctx);
1521     n_poly_init2(tmp, lgd);
1522     n_poly_init2(alphapow, 2*lgd);
1523     fq_nmod_init(start_alpha, lgctx->fqctx);
1524 
1525     /* Aevals[nvars] does not exist - it is replaced by An */
1526     Aevals = FLINT_ARRAY_ALLOC(nvars, fq_nmod_mpoly_struct);
1527     Bevals = FLINT_ARRAY_ALLOC(nvars, fq_nmod_mpoly_struct);
1528     gammaevals = FLINT_ARRAY_ALLOC(nvars, fq_nmod_mpoly_struct);
1529     for (i = 0; i < nvars; i++)
1530     {
1531         fq_nmod_mpoly_init3(Aevals + i, 0, bits, smctx);
1532         fq_nmod_mpoly_init3(Bevals + i, 0, bits, smctx);
1533         fq_nmod_mpoly_init3(gammaevals + i, 0, bits, smctx);
1534     }
1535     fq_nmod_mpoly_cvtto_mpolyn(An, A, nvars - 1, smctx);
1536     fq_nmod_mpoly_cvtto_mpolyn(Bn, B, nvars - 1, smctx);
1537     fq_nmod_mpoly_cvtto_mpolyn(gamman, gamma, nvars - 1, smctx);
1538 
1539     Abideg = _fq_nmod_mpoly_bidegree(A, smctx);
1540     Bbideg = _fq_nmod_mpoly_bidegree(B, smctx);
1541 
1542     degxAB = FLINT_MAX(Adegs[0], Bdegs[0]);
1543     degyAB = FLINT_MAX(Adegs[1], Bdegs[1]);
1544 
1545     GdegboundXY = pack_exp2(FLINT_MIN(Adegs[0], Bdegs[0]),
1546                             FLINT_MIN(Adegs[1], Bdegs[1]));
1547     if (GdegboundXY == 0)
1548         goto gcd_is_trivial;
1549 
1550     goto got_alpha_m;
1551 
1552 increase_degree:
1553 
1554     do {
1555         cur_emb = bad_fq_nmod_mpoly_embed_chooser_next(embc, lgctx, smctx, state);
1556         if (cur_emb == NULL)
1557         {
1558             /* ran out of primes */
1559             success = 0;
1560             goto cleanup;
1561         }
1562     } while (fq_nmod_ctx_degree(lgctx->fqctx) <= lgd);
1563 
1564     lgd = fq_nmod_ctx_degree(lgctx->fqctx);
1565 
1566     goto got_alpha_m;
1567 
1568 choose_main:
1569 
1570     cur_emb = bad_fq_nmod_mpoly_embed_chooser_next(embc, lgctx, smctx, state);
1571     if (cur_emb == NULL)
1572     {
1573         /* ran out of primes */
1574         success = 0;
1575         goto cleanup;
1576     }
1577 
1578     lgd = fq_nmod_ctx_degree(lgctx->fqctx);
1579 
1580 got_alpha_m:
1581 
1582     n_poly_fit_length(tmp, lgd);
1583     n_poly_fit_length(alphapow, 2*lgd);
1584 
1585     for (i = 2; i < nvars - 1; i++)
1586         fq_nmod_rand_not_zero(alphas + i, state, lgctx->fqctx);
1587 
1588     i = nvars - 1;
1589     fq_nmod_mpolyn_interp_reduce_lg_mpoly(Aevals + i, An, lgctx, smctx, cur_emb);
1590     fq_nmod_mpolyn_interp_reduce_lg_mpoly(Bevals + i, Bn, lgctx, smctx, cur_emb);
1591     fq_nmod_mpolyn_interp_reduce_lg_mpoly(gammaevals + i, gamman, lgctx, smctx, cur_emb);
1592     if (fq_nmod_mpoly_is_zero(gammaevals + i, lgctx))
1593         goto choose_main;
1594     if (Aevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + i, lgctx) != Abideg)
1595         goto choose_main;
1596     if (Bevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + i, lgctx) != Bbideg)
1597         goto choose_main;
1598     for (i--; i >= 2; i--)
1599     {
1600         fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + i, Aevals + i + 1, i, alphas + i, lgctx);
1601         fq_nmod_mpoly_repack_bits_inplace(Aevals + i, bits, lgctx);
1602         fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + i, Bevals + i + 1, i, alphas + i, lgctx);
1603         fq_nmod_mpoly_repack_bits_inplace(Bevals + i, bits, lgctx);
1604         fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + i, gammaevals + i + 1, i, alphas + i, lgctx);
1605         fq_nmod_mpoly_repack_bits_inplace(gammaevals + i, bits, lgctx);
1606         if (fq_nmod_mpoly_is_zero(gammaevals + i, lgctx))
1607             goto choose_main;
1608         if (Aevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + i, lgctx) != Abideg)
1609             goto choose_main;
1610         if (Bevals[i].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + i, lgctx) != Bbideg)
1611             goto choose_main;
1612     }
1613 
1614     m = 2;
1615 
1616     if (rGdegs == NULL)
1617         use = USE_G | USE_ABAR | USE_BBAR;
1618     else
1619         use = mpoly_gcd_get_use_first(rGdegs[m], Adegs[m], Bdegs[m], gammadegs[m]);
1620 
1621     fq_nmod_mpoly_get_n_fq_bpoly(Aev, Aevals + m, 0, 1, lgctx);
1622     fq_nmod_mpoly_get_n_fq_bpoly(Bev, Bevals + m, 0, 1, lgctx);
1623 
1624     success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
1625                                                    Aev, Bev, lgctx->fqctx, St);
1626     if (!success)
1627         goto increase_degree;
1628 
1629     newdegXY = n_bpoly_bidegree(Gev);
1630     if (newdegXY > GdegboundXY)
1631         goto choose_main;
1632     if (newdegXY < GdegboundXY)
1633     {
1634         GdegboundXY = newdegXY;
1635         if (GdegboundXY == 0)
1636             goto gcd_is_trivial;
1637     }
1638 
1639     gammaev = fq_nmod_mpoly_get_nonzero_n_fq(gammaevals + m, lgctx);
1640     n_fq_bpoly_scalar_mul_n_fq(Gev, gammaev, lgctx->fqctx);
1641 
1642     if (nvars == 3)
1643     {
1644         fq_nmod_mpolyn_interp_lift_lg_bpoly(&lastdeg, Gn, smctx, Gev, lgctx, cur_emb);
1645         fq_nmod_mpolyn_interp_lift_lg_bpoly(&lastdeg, Abarn, smctx, Abarev, lgctx, cur_emb);
1646         fq_nmod_mpolyn_interp_lift_lg_bpoly(&lastdeg, Bbarn, smctx, Bbarev, lgctx, cur_emb);
1647 
1648         n_fq_poly_set(modulus, cur_emb->h_as_n_fq_poly, smctx->fqctx);
1649 
1650         while (1)
1651         {
1652         choose_alpha_2_last:
1653 
1654             cur_emb = bad_fq_nmod_mpoly_embed_chooser_next(embc, lgctx, smctx, state);
1655             if (cur_emb == NULL)
1656             {
1657                 /* ran out of primes */
1658                 success = 0;
1659                 goto cleanup;
1660             }
1661 
1662             lgd = fq_nmod_ctx_degree(lgctx->fqctx);
1663             n_poly_fit_length(tmp, lgd);
1664             n_poly_fit_length(alphapow, 2*lgd);
1665 
1666             fq_nmod_mpolyn_interp_reduce_lg_mpoly(Aevals + m, An, lgctx, smctx, cur_emb);
1667             fq_nmod_mpolyn_interp_reduce_lg_mpoly(Bevals + m, Bn, lgctx, smctx, cur_emb);
1668             fq_nmod_mpolyn_interp_reduce_lg_mpoly(gammaevals + m, gamman, lgctx, smctx, cur_emb);
1669             if (fq_nmod_mpoly_is_zero(gammaevals + m, lgctx))
1670                 goto choose_alpha_2_last;
1671             if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, lgctx) != Abideg)
1672                 goto choose_alpha_2_last;
1673             if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, lgctx) != Bbideg)
1674                 goto choose_alpha_2_last;
1675 
1676             fq_nmod_mpoly_get_n_fq_bpoly(Aev, Aevals + m, 0, 1, lgctx);
1677             fq_nmod_mpoly_get_n_fq_bpoly(Bev, Bevals + m, 0, 1, lgctx);
1678 
1679             success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
1680                                                    Aev, Bev, lgctx->fqctx, St);
1681             if (!success)
1682                 goto increase_degree;
1683 
1684             newdegXY = n_bpoly_bidegree(Gev);
1685             if (newdegXY > GdegboundXY)
1686                 goto choose_alpha_2_last;
1687             if (newdegXY < GdegboundXY)
1688             {
1689                 GdegboundXY = newdegXY;
1690                 if (GdegboundXY == 0)
1691                     goto gcd_is_trivial;
1692                 goto choose_main;
1693             }
1694 
1695             gammaev = fq_nmod_mpoly_get_nonzero_n_fq(gammaevals + m, lgctx);
1696             n_fq_bpoly_scalar_mul_n_fq(Gev, gammaev, lgctx->fqctx);
1697 
1698             if ((use & USE_G) && !fq_nmod_mpolyn_interp_crt_lg_bpoly(
1699                         &lastdeg, Gn, Tn, modulus, smctx, Gev, lgctx, cur_emb))
1700             {
1701                 fq_nmod_mpoly_cvtfrom_mpolyn(rG, Gn, m, smctx);
1702                 success = fq_nmod_mpolyl_content(cont, rG, 2, smctx);
1703                 if (!success)
1704                     goto cleanup;
1705                 fq_nmod_mpoly_divides(rG, rG, cont, smctx);
1706                 fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
1707                 if (fq_nmod_mpoly_divides(rAbar, A, rG, smctx) &&
1708                     fq_nmod_mpoly_divides(rBbar, B, rG, smctx))
1709                 {
1710                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
1711                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
1712                     break;
1713                 }
1714             }
1715 
1716             if ((use & USE_ABAR) && !fq_nmod_mpolyn_interp_crt_lg_bpoly(
1717                   &lastdeg, Abarn, Tn, modulus, smctx, Abarev, lgctx, cur_emb))
1718             {
1719                 fq_nmod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, smctx);
1720                 success = fq_nmod_mpolyl_content(cont, rAbar, 2, smctx);
1721                 if (!success)
1722                     goto cleanup;
1723                 fq_nmod_mpoly_divides(rAbar, rAbar, cont, smctx);
1724                 fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
1725                 if (fq_nmod_mpoly_divides(rG, A, rAbar, smctx) &&
1726                     fq_nmod_mpoly_divides(rBbar, B, rG, smctx))
1727                 {
1728                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
1729                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
1730                     break;
1731                 }
1732             }
1733 
1734             if ((use & USE_BBAR) && !fq_nmod_mpolyn_interp_crt_lg_bpoly(
1735                   &lastdeg, Bbarn, Tn, modulus, smctx, Bbarev, lgctx, cur_emb))
1736             {
1737                 fq_nmod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, smctx);
1738                 success = fq_nmod_mpolyl_content(cont, rBbar, 2, smctx);
1739                 if (!success)
1740                     goto cleanup;
1741                 fq_nmod_mpoly_divides(rBbar, rBbar, cont, smctx);
1742                 fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
1743                 if (fq_nmod_mpoly_divides(rG, B, rBbar, smctx) &&
1744                     fq_nmod_mpoly_divides(rAbar, A, rG, smctx))
1745                 {
1746                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
1747                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
1748                     break;
1749                 }
1750             }
1751 
1752             if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
1753                 n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
1754             {
1755                 goto choose_main;
1756             }
1757 
1758             n_fq_poly_mul_(modulus, modulus, cur_emb->h_as_n_fq_poly,
1759                                                  smctx->fqctx, St->poly_stack);
1760         }
1761 
1762         success = 1;
1763         goto cleanup;
1764     }
1765 
1766     fq_nmod_mpolyn_interp_lift_sm_bpoly(Gn, Gev, lgctx);
1767     fq_nmod_mpolyn_interp_lift_sm_bpoly(Abarn, Abarev, lgctx);
1768     fq_nmod_mpolyn_interp_lift_sm_bpoly(Bbarn, Bbarev, lgctx);
1769 
1770     FLINT_ASSERT(tmp->alloc >= lgd);
1771     n_fq_poly_one(modulus, lgctx->fqctx);
1772     n_fq_set_fq_nmod(tmp->coeffs, alphas + m, lgctx->fqctx);
1773     n_fq_poly_shift_left_scalar_submul(modulus, 1, tmp->coeffs, lgctx->fqctx);
1774 
1775     fq_nmod_set(start_alpha, alphas + m, lgctx->fqctx);
1776 
1777     while (1)
1778     {
1779     choose_alpha_2:
1780 
1781         fq_nmod_next_not_zero(alphas + m, lgctx->fqctx);
1782         if (fq_nmod_equal(alphas + m, start_alpha, lgctx->fqctx))
1783             goto increase_degree;
1784 
1785         FLINT_ASSERT(alphapow->alloc >= lgd*2);
1786         alphapow->length = 2;
1787         _n_fq_one(alphapow->coeffs + lgd*0, lgd);
1788         n_fq_set_fq_nmod(alphapow->coeffs + lgd*1, alphas + m, lgctx->fqctx);
1789 
1790         fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + m, Aevals + m + 1, m, alphas + m, lgctx);
1791         fq_nmod_mpoly_repack_bits_inplace(Aevals + m, bits, lgctx);
1792         fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + m, Bevals + m + 1, m, alphas + m, lgctx);
1793         fq_nmod_mpoly_repack_bits_inplace(Bevals + m, bits, lgctx);
1794         fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + m, gammaevals + m + 1, m, alphas + m, lgctx);
1795         fq_nmod_mpoly_repack_bits_inplace(gammaevals + m, bits, lgctx);
1796         if (fq_nmod_mpoly_is_zero(gammaevals + m, lgctx))
1797             goto choose_alpha_2;
1798         if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, lgctx) != Abideg)
1799             goto choose_alpha_2;
1800         if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, lgctx) != Bbideg)
1801             goto choose_alpha_2;
1802 
1803         fq_nmod_mpoly_get_n_fq_bpoly(Aev, Aevals + m, 0, 1, lgctx);
1804         fq_nmod_mpoly_get_n_fq_bpoly(Bev, Bevals + m, 0, 1, lgctx);
1805 
1806         success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
1807                                                    Aev, Bev, lgctx->fqctx, St);
1808         if (!success)
1809             goto increase_degree;
1810 
1811         newdegXY = n_bpoly_bidegree(Gev);
1812         if (newdegXY > GdegboundXY)
1813             goto choose_alpha_2;
1814         if (newdegXY < GdegboundXY)
1815         {
1816             GdegboundXY = newdegXY;
1817             if (GdegboundXY == 0)
1818                 goto gcd_is_trivial;
1819             goto choose_main;
1820         }
1821 
1822         gammaev = fq_nmod_mpoly_get_nonzero_n_fq(gammaevals + m, lgctx);
1823         n_fq_bpoly_scalar_mul_n_fq(Gev, gammaev, lgctx->fqctx);
1824 
1825         FLINT_ASSERT(tmp->alloc >= lgd);
1826         n_fq_poly_eval_pow(tmp->coeffs, modulus, alphapow, lgctx->fqctx);
1827         n_fq_inv(tmp->coeffs, tmp->coeffs, lgctx->fqctx);
1828         n_fq_poly_scalar_mul_n_fq(modulus, modulus, tmp->coeffs, lgctx->fqctx);
1829 
1830         if ((use & USE_G) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
1831                               &lastdeg, Gn, Tn, Gev, modulus, alphapow, lgctx))
1832         {
1833             fq_nmod_mpoly_cvtfrom_mpolyn(G, Gn, m, lgctx);
1834             fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
1835             if (fq_nmod_mpoly_divides(Abar, T, G, lgctx))
1836             {
1837                 fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
1838                 if (fq_nmod_mpoly_divides(Bbar, T, G, lgctx))
1839                 {
1840                     fq_nmod_mpoly_repack_bits_inplace(Abar, bits, lgctx);
1841                     fq_nmod_mpoly_repack_bits_inplace(Bbar, bits, lgctx);
1842                     break;
1843                 }
1844             }
1845         }
1846 
1847         if ((use & USE_ABAR) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
1848                         &lastdeg, Abarn, Tn, Abarev, modulus, alphapow, lgctx))
1849         {
1850             fq_nmod_mpoly_cvtfrom_mpolyn(Abar, Abarn, m, lgctx);
1851             fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
1852             if (fq_nmod_mpoly_divides(G, T, Abar, lgctx))
1853             {
1854                 fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
1855                 if (fq_nmod_mpoly_divides(Bbar, T, G, lgctx))
1856                 {
1857                     fq_nmod_mpoly_repack_bits_inplace(G, bits, lgctx);
1858                     fq_nmod_mpoly_repack_bits_inplace(Bbar, bits, lgctx);
1859                     break;
1860                 }
1861             }
1862         }
1863 
1864         if ((use & USE_BBAR) && !fq_nmod_mpolyn_interp_crt_sm_bpoly(
1865                         &lastdeg, Bbarn, Tn, Bbarev, modulus, alphapow, lgctx))
1866         {
1867             fq_nmod_mpoly_cvtfrom_mpolyn(Bbar, Bbarn, m, lgctx);
1868             fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
1869             if (fq_nmod_mpoly_divides(G, T, Bbar, lgctx))
1870             {
1871                 fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
1872                 if (fq_nmod_mpoly_divides(Abar, T, G, lgctx))
1873                 {
1874                     fq_nmod_mpoly_repack_bits_inplace(G, bits, lgctx);
1875                     fq_nmod_mpoly_repack_bits_inplace(Abar, bits, lgctx);
1876                     break;
1877                 }
1878             }
1879         }
1880 
1881         if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
1882             n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
1883         {
1884             goto choose_main;
1885         }
1886 
1887         FLINT_ASSERT(n_fq_equal_fq_nmod(alphapow->coeffs + lgd*1, alphas + m, lgctx->fqctx));
1888         n_fq_poly_shift_left_scalar_submul(modulus, 1, alphapow->coeffs + lgd*1, lgctx->fqctx);
1889     }
1890 
1891     for (m = 3; m < nvars - 1; m++)
1892     {
1893         /* G, Abar, Bbar are in Fq[gen(0), ..., gen(m - 1)] */
1894         fq_nmod_mpolyn_interp_lift_sm_mpoly(Gn, G, lgctx);
1895         fq_nmod_mpolyn_interp_lift_sm_mpoly(Abarn, Abar, lgctx);
1896         fq_nmod_mpolyn_interp_lift_sm_mpoly(Bbarn, Bbar, lgctx);
1897 
1898         if (rGdegs == NULL)
1899             use = USE_G | USE_ABAR | USE_BBAR;
1900         else
1901             use = 0;
1902 
1903         FLINT_ASSERT(tmp->alloc >= lgd);
1904         n_fq_poly_one(modulus, lgctx->fqctx);
1905         n_fq_set_fq_nmod(tmp->coeffs, alphas + m, lgctx->fqctx);
1906         n_fq_poly_shift_left_scalar_submul(modulus, 1, tmp->coeffs, lgctx->fqctx);
1907 
1908         fq_nmod_set(start_alpha, alphas + m, lgctx->fqctx);
1909 
1910     choose_betas_m:
1911 
1912         /* only beta[2], beta[1], ..., beta[m - 1] will be used */
1913         for (i = 2; i < nvars; i++)
1914             fq_nmod_rand_not_zero(betas + i, state, lgctx->fqctx);
1915 
1916         FLINT_ASSERT(G->bits == bits);
1917         FLINT_ASSERT(Abar->bits == bits);
1918         FLINT_ASSERT(Bbar->bits == bits);
1919 
1920         fq_nmod_mpoly_monomial_evals2(HG, G, betas + 2, m, lgctx);
1921         fq_nmod_mpoly_monomial_evals2(HAbar, Abar, betas + 2, m, lgctx);
1922         fq_nmod_mpoly_monomial_evals2(HBbar, Bbar, betas + 2, m, lgctx);
1923         if (use == 0)
1924         {
1925             this_length = gammaevals[m + 1].length + Aevals[m + 1].length +
1926                                                      Bevals[m + 1].length;
1927 
1928             use = nmod_mpoly_gcd_get_use_new(rGdegs[m], Adegs[m], Bdegs[m],
1929                    gammadegs[m], degxAB, degyAB, this_length, HG, HAbar, HBbar);
1930         }
1931 
1932         req_zip_images = 1;
1933         if (use & USE_G)
1934         {
1935             this_length = n_fq_polyun_product_roots(MG, HG, lgctx->fqctx, St->poly_stack);
1936             req_zip_images = FLINT_MAX(req_zip_images, this_length);
1937         }
1938         if (use & USE_ABAR)
1939         {
1940             this_length = n_fq_polyun_product_roots(MAbar, HAbar, lgctx->fqctx, St->poly_stack);
1941             req_zip_images = FLINT_MAX(req_zip_images, this_length);
1942         }
1943         if (use & USE_BBAR)
1944         {
1945             this_length = n_fq_polyun_product_roots(MBbar, HBbar, lgctx->fqctx, St->poly_stack);
1946             req_zip_images = FLINT_MAX(req_zip_images, this_length);
1947         }
1948 
1949         while (1)
1950         {
1951         choose_alpha_m:
1952 
1953             fq_nmod_next_not_zero(alphas + m, lgctx->fqctx);
1954             if (fq_nmod_equal(alphas + m, start_alpha, lgctx->fqctx))
1955                 goto choose_main;
1956 
1957             FLINT_ASSERT(alphapow->alloc >= lgd*2);
1958             alphapow->length = 2;
1959             _n_fq_one(alphapow->coeffs + lgd*0, lgd);
1960             n_fq_set_fq_nmod(alphapow->coeffs + lgd*1, alphas + m, lgctx->fqctx);
1961 
1962             fq_nmod_mpoly_evaluate_one_fq_nmod(Aevals + m, Aevals + m + 1, m, alphas + m, lgctx);
1963             fq_nmod_mpoly_repack_bits_inplace(Aevals + m, bits, lgctx);
1964             fq_nmod_mpoly_evaluate_one_fq_nmod(Bevals + m, Bevals + m + 1, m, alphas + m, lgctx);
1965             fq_nmod_mpoly_repack_bits_inplace(Bevals + m, bits, lgctx);
1966             fq_nmod_mpoly_evaluate_one_fq_nmod(gammaevals + m, gammaevals + m + 1, m, alphas + m, lgctx);
1967             fq_nmod_mpoly_repack_bits_inplace(gammaevals + m, bits, lgctx);
1968             if (fq_nmod_mpoly_is_zero(gammaevals + m, lgctx))
1969                 goto choose_alpha_m;
1970             if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, lgctx) != Abideg)
1971                 goto choose_alpha_m;
1972             if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, lgctx) != Bbideg)
1973                 goto choose_alpha_m;
1974 
1975             fq_nmod_mpoly_monomial_evals2(Aeh_inc, Aevals + m, betas + 2, m, lgctx);
1976             fq_nmod_mpoly_monomial_evals2(Beh_inc, Bevals + m, betas + 2, m, lgctx);
1977             fq_nmod_mpoly_monomial_evals(gammaeh_inc, gammaevals + m, betas + 2, 2, m, lgctx);
1978             n_fq_polyun_set(Aeh_cur, Aeh_inc, lgctx->fqctx);
1979             n_fq_polyun_set(Beh_cur, Beh_inc, lgctx->fqctx);
1980             n_fq_poly_set(gammaeh_cur, gammaeh_inc, lgctx->fqctx);
1981 
1982             n_fq_polyun_zip_start(ZG, HG, req_zip_images, lgctx->fqctx);
1983             n_fq_polyun_zip_start(ZAbar, HAbar, req_zip_images, lgctx->fqctx);
1984             n_fq_polyun_zip_start(ZBbar, HBbar, req_zip_images, lgctx->fqctx);
1985             for (cur_zip_image = 0; cur_zip_image < req_zip_images; cur_zip_image++)
1986             {
1987                 n_fq_bpoly_eval_step_sep(Aev, Aeh_cur, Aeh_inc, Aevals + m, lgctx->fqctx);
1988                 n_fq_bpoly_eval_step_sep(Bev, Beh_cur, Beh_inc, Bevals + m, lgctx->fqctx);
1989                 FLINT_ASSERT(tmp->alloc >= lgd);
1990                 n_fq_poly_eval_step_sep(tmp->coeffs, gammaeh_cur, gammaeh_inc, gammaevals + m, lgctx->fqctx);
1991                 if (_n_fq_is_zero(tmp->coeffs, lgd))
1992                     goto choose_betas_m;
1993                 if (Aev->length < 1 || n_bpoly_bidegree(Aev) != Abideg)
1994                     goto choose_betas_m;
1995                 if (Bev->length < 1 || n_bpoly_bidegree(Bev) != Bbideg)
1996                     goto choose_betas_m;
1997 
1998                 success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
1999                                                    Aev, Bev, lgctx->fqctx, St);
2000                 if (!success)
2001                     goto increase_degree;
2002 
2003                 newdegXY = n_bpoly_bidegree(Gev);
2004                 if (newdegXY > GdegboundXY)
2005                     goto choose_betas_m;
2006                 if (newdegXY < GdegboundXY)
2007                 {
2008                     GdegboundXY = newdegXY;
2009                     if (GdegboundXY == 0)
2010                         goto gcd_is_trivial;
2011                     goto choose_main;
2012                 }
2013 
2014                 n_fq_bpoly_scalar_mul_n_fq(Gev, tmp->coeffs, lgctx->fqctx);
2015 
2016                 if ((use & USE_G) && !n_fq_polyu2n_add_zip_must_match(ZG, Gev, cur_zip_image, lgctx->fqctx))
2017                     goto choose_main;
2018                 if ((use & USE_ABAR) && !n_fq_polyu2n_add_zip_must_match(ZAbar, Abarev, cur_zip_image, lgctx->fqctx))
2019                     goto choose_main;
2020                 if ((use & USE_BBAR) && !n_fq_polyu2n_add_zip_must_match(ZBbar, Bbarev, cur_zip_image, lgctx->fqctx))
2021                     goto choose_main;
2022             }
2023 
2024             if ((use & USE_G) && n_fq_polyun_zip_solve(G, ZG, HG, MG, lgctx) < 1)
2025                 goto choose_main;
2026             if ((use & USE_ABAR) && n_fq_polyun_zip_solve(Abar, ZAbar, HAbar, MAbar, lgctx) < 1)
2027                 goto choose_main;
2028             if ((use & USE_BBAR) && n_fq_polyun_zip_solve(Bbar, ZBbar, HBbar, MBbar, lgctx) < 1)
2029                 goto choose_main;
2030 
2031             FLINT_ASSERT(n_fq_poly_degree(modulus) > 0);
2032 
2033             FLINT_ASSERT(tmp->alloc >= lgd);
2034             n_fq_poly_eval_pow(tmp->coeffs, modulus, alphapow, lgctx->fqctx);
2035             n_fq_inv(tmp->coeffs, tmp->coeffs, lgctx->fqctx);
2036             n_fq_poly_scalar_mul_n_fq(modulus, modulus, tmp->coeffs, lgctx->fqctx);
2037 
2038             if ((use & USE_G) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
2039                                     &lastdeg, Gn, G, modulus, alphapow, lgctx))
2040             {
2041                 fq_nmod_mpoly_cvtfrom_mpolyn(rG, Gn, m, lgctx);
2042                 fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
2043                 if (fq_nmod_mpoly_divides(rAbar, T, rG, lgctx))
2044                 {
2045                     fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
2046                     if (fq_nmod_mpoly_divides(rBbar, T, rG, lgctx))
2047                     {
2048                         fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, lgctx);
2049                         fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, lgctx);
2050                         fq_nmod_mpoly_swap(G, rG, lgctx);
2051                         fq_nmod_mpoly_swap(Abar, rAbar, lgctx);
2052                         fq_nmod_mpoly_swap(Bbar, rBbar, lgctx);
2053                         break;
2054                     }
2055                 }
2056             }
2057 
2058             if ((use & USE_ABAR) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
2059                               &lastdeg, Abarn, Abar, modulus, alphapow, lgctx))
2060             {
2061                 fq_nmod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, lgctx);
2062                 fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
2063                 if (fq_nmod_mpoly_divides(rG, T, rAbar, lgctx))
2064                 {
2065                     fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
2066                     if (fq_nmod_mpoly_divides(rBbar, T, rG, lgctx))
2067                     {
2068                         fq_nmod_mpoly_repack_bits_inplace(rG, bits, lgctx);
2069                         fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, lgctx);
2070                         fq_nmod_mpoly_swap(G, rG, lgctx);
2071                         fq_nmod_mpoly_swap(Abar, rAbar, lgctx);
2072                         fq_nmod_mpoly_swap(Bbar, rBbar, lgctx);
2073                         break;
2074                     }
2075                 }
2076             }
2077 
2078             if ((use & USE_BBAR) && !fq_nmod_mpolyn_interp_mcrt_sm_mpoly(
2079                               &lastdeg, Bbarn, Bbar, modulus, alphapow, lgctx))
2080             {
2081                 fq_nmod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, lgctx);
2082                 fq_nmod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, lgctx);
2083                 if (fq_nmod_mpoly_divides(rG, T, rBbar, lgctx))
2084                 {
2085                     fq_nmod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, lgctx);
2086                     if (fq_nmod_mpoly_divides(rAbar, T, rG, lgctx))
2087                     {
2088                         fq_nmod_mpoly_repack_bits_inplace(rG, bits, lgctx);
2089                         fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, lgctx);
2090                         fq_nmod_mpoly_swap(G, rG, lgctx);
2091                         fq_nmod_mpoly_swap(Abar, rAbar, lgctx);
2092                         fq_nmod_mpoly_swap(Bbar, rBbar, lgctx);
2093                         break;
2094                     }
2095                 }
2096             }
2097 
2098             if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
2099                 n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
2100             {
2101                 goto choose_main;
2102             }
2103 
2104             FLINT_ASSERT(n_fq_equal_fq_nmod(alphapow->coeffs + lgd*1, alphas + m, lgctx->fqctx));
2105             n_fq_poly_shift_left_scalar_submul(modulus, 1, alphapow->coeffs + lgd*1, lgctx->fqctx);
2106         }
2107     }
2108 
2109     m = nvars - 1;
2110     {
2111         /* G, Abar, Bbar are in Fq/alpha(gen(m-1))[gen(0), ..., gen(m - 1)] */
2112 
2113         fq_nmod_mpolyn_interp_lift_lg_mpoly(Gn, smctx, G, lgctx, cur_emb);
2114         fq_nmod_mpolyn_interp_lift_lg_mpoly(Abarn, smctx, Abar, lgctx, cur_emb);
2115         fq_nmod_mpolyn_interp_lift_lg_mpoly(Bbarn, smctx, Bbar, lgctx, cur_emb);
2116 
2117         if (rGdegs == NULL)
2118             use = USE_G | USE_ABAR | USE_BBAR;
2119         else
2120             use = 0;
2121 
2122         n_fq_poly_set(modulus, cur_emb->h_as_n_fq_poly, smctx->fqctx);
2123 
2124         while (1)
2125         {
2126         choose_alpha_last:
2127 
2128             cur_emb = bad_fq_nmod_mpoly_embed_chooser_next(embc, lgctx, smctx, state);
2129             if (cur_emb == NULL)
2130             {
2131                 /* ran out of primes */
2132                 success = 0;
2133                 goto cleanup;
2134             }
2135 
2136             lgd = fq_nmod_ctx_degree(lgctx->fqctx);
2137             n_poly_fit_length(tmp, lgd);
2138             n_poly_fit_length(alphapow, 2*lgd);
2139 
2140             fq_nmod_mpolyn_interp_reduce_lg_mpoly(Aevals + m, An, lgctx, smctx, cur_emb);
2141             fq_nmod_mpolyn_interp_reduce_lg_mpoly(Bevals + m, Bn, lgctx, smctx, cur_emb);
2142             fq_nmod_mpolyn_interp_reduce_lg_mpoly(gammaevals + m, gamman, lgctx, smctx, cur_emb);
2143             if (fq_nmod_mpoly_is_zero(gammaevals + m, lgctx))
2144                 goto choose_alpha_last;
2145             if (Aevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Aevals + m, lgctx) != Abideg)
2146                 goto choose_alpha_last;
2147             if (Bevals[m].length < 1 || _fq_nmod_mpoly_bidegree(Bevals + m, lgctx) != Bbideg)
2148                 goto choose_alpha_last;
2149 
2150         choose_betas_last:
2151 
2152             /* only beta[2], ..., beta[m - 1] will be used */
2153             for (i = 2; i < nvars; i++)
2154                 fq_nmod_rand_not_zero(betas + i, state, lgctx->fqctx);
2155 
2156             fq_nmod_mpoly_monomial_evals2(HG, G, betas + 2, m, lgctx);
2157             fq_nmod_mpoly_monomial_evals2(HAbar, Abar, betas + 2, m, lgctx);
2158             fq_nmod_mpoly_monomial_evals2(HBbar, Bbar, betas + 2, m, lgctx);
2159 
2160             if (use == 0)
2161             {
2162                 this_length = gamma->length + A->length + B->length;
2163                 use = nmod_mpoly_gcd_get_use_new(rGdegs[m], Adegs[m], Bdegs[m],
2164                        gammadegs[m], degxAB, degyAB, this_length, HG, HAbar, HBbar);
2165             }
2166 
2167             req_zip_images = 1;
2168             if (use & USE_G)
2169             {
2170                 this_length = n_fq_polyun_product_roots(MG, HG, lgctx->fqctx, St->poly_stack);
2171                 req_zip_images = FLINT_MAX(req_zip_images, this_length);
2172             }
2173             if (use & USE_ABAR)
2174             {
2175                 this_length = n_fq_polyun_product_roots(MAbar, HAbar, lgctx->fqctx, St->poly_stack);
2176                 req_zip_images = FLINT_MAX(req_zip_images, this_length);
2177             }
2178             if (use & USE_BBAR)
2179             {
2180                 this_length = n_fq_polyun_product_roots(MBbar, HBbar, lgctx->fqctx, St->poly_stack);
2181                 req_zip_images = FLINT_MAX(req_zip_images, this_length);
2182             }
2183 
2184             fq_nmod_mpoly_monomial_evals2(Aeh_inc, Aevals + m, betas + 2, m, lgctx);
2185             fq_nmod_mpoly_monomial_evals2(Beh_inc, Bevals + m, betas + 2, m, lgctx);
2186             fq_nmod_mpoly_monomial_evals(gammaeh_inc, gammaevals + m, betas + 2, 2, m, lgctx);
2187             n_fq_polyun_set(Aeh_cur, Aeh_inc, lgctx->fqctx);
2188             n_fq_polyun_set(Beh_cur, Beh_inc, lgctx->fqctx);
2189             n_fq_poly_set(gammaeh_cur, gammaeh_inc, lgctx->fqctx);
2190 
2191             n_fq_polyun_zip_start(ZG, HG, req_zip_images, lgctx->fqctx);
2192             n_fq_polyun_zip_start(ZAbar, HAbar, req_zip_images, lgctx->fqctx);
2193             n_fq_polyun_zip_start(ZBbar, HBbar, req_zip_images, lgctx->fqctx);
2194             for (cur_zip_image = 0; cur_zip_image < req_zip_images; cur_zip_image++)
2195             {
2196                 n_fq_bpoly_eval_step_sep(Aev, Aeh_cur, Aeh_inc, Aevals + m, lgctx->fqctx);
2197                 n_fq_bpoly_eval_step_sep(Bev, Beh_cur, Beh_inc, Bevals + m, lgctx->fqctx);
2198                 FLINT_ASSERT(tmp->alloc >= lgd);
2199                 n_fq_poly_eval_step_sep(tmp->coeffs, gammaeh_cur, gammaeh_inc, gammaevals + m, lgctx->fqctx);
2200                 if (_n_fq_is_zero(tmp->coeffs, lgd))
2201                     goto choose_betas_last;
2202                 if (Aev->length < 1 || n_bpoly_bidegree(Aev) != Abideg)
2203                     goto choose_betas_last;
2204                 if (Bev->length < 1 || n_bpoly_bidegree(Bev) != Bbideg)
2205                     goto choose_betas_last;
2206 
2207                 success = n_fq_bpoly_gcd_brown_smprime(Gev, Abarev, Bbarev,
2208                                                    Aev, Bev, lgctx->fqctx, St);
2209                 if (!success)
2210                     goto cleanup;
2211 
2212                 newdegXY = n_bpoly_bidegree(Gev);
2213                 if (newdegXY > GdegboundXY)
2214                     goto choose_betas_last;
2215                 if (newdegXY < GdegboundXY)
2216                 {
2217                     GdegboundXY = newdegXY;
2218                     if (GdegboundXY == 0)
2219                         goto gcd_is_trivial;
2220                     goto choose_main;
2221                 }
2222 
2223                 n_fq_bpoly_scalar_mul_n_fq(Gev, tmp->coeffs, lgctx->fqctx);
2224 
2225                 if ((use & USE_G) && !n_fq_polyu2n_add_zip_must_match(ZG, Gev, cur_zip_image, lgctx->fqctx))
2226                     goto choose_main;
2227                 if ((use & USE_ABAR) && !n_fq_polyu2n_add_zip_must_match(ZAbar, Abarev, cur_zip_image, lgctx->fqctx))
2228                     goto choose_main;
2229                 if ((use & USE_BBAR) && !n_fq_polyu2n_add_zip_must_match(ZBbar, Bbarev, cur_zip_image, lgctx->fqctx))
2230                     goto choose_main;
2231             }
2232             if ((use & USE_G) && n_fq_polyun_zip_solve(G, ZG, HG, MG, lgctx) < 1)
2233                 goto choose_main;
2234             if ((use & USE_ABAR) && n_fq_polyun_zip_solve(Abar, ZAbar, HAbar, MAbar, lgctx) < 1)
2235                 goto choose_main;
2236             if ((use & USE_BBAR) && n_fq_polyun_zip_solve(Bbar, ZBbar, HBbar, MBbar, lgctx) < 1)
2237                 goto choose_main;
2238 
2239             FLINT_ASSERT(n_fq_poly_degree(modulus) > 0);
2240 
2241             FLINT_ASSERT(tmp->alloc >= lgd);
2242             bad_n_fq_embed_sm_to_lg(tmp->coeffs, modulus, cur_emb);
2243             n_fq_inv(tmp->coeffs, tmp->coeffs, lgctx->fqctx);
2244 
2245             if ((use & USE_G) && !newfq_nmod_mpolyn_interp_mcrt_lg_mpoly(
2246                                     &lastdeg, Gn, smctx, modulus, tmp->coeffs,
2247                                                             G, lgctx, cur_emb))
2248             {
2249                 fq_nmod_mpoly_cvtfrom_mpolyn(rG, Gn, m, smctx);
2250                 success = fq_nmod_mpolyl_content(cont, rG, 2, smctx);
2251                 if (!success)
2252                     goto cleanup;
2253                 fq_nmod_mpoly_divides(rG, rG, cont, smctx);
2254                 fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
2255                 if (fq_nmod_mpoly_divides(rAbar, A, rG, smctx) &&
2256                     fq_nmod_mpoly_divides(rBbar, B, rG, smctx))
2257                 {
2258                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
2259                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
2260                     break;
2261                 }
2262             }
2263 
2264             if ((use & USE_ABAR) && !newfq_nmod_mpolyn_interp_mcrt_lg_mpoly(
2265                                  &lastdeg, Abarn, smctx, modulus, tmp->coeffs,
2266                                                          Abar, lgctx, cur_emb))
2267             {
2268                 fq_nmod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, smctx);
2269                 success = fq_nmod_mpolyl_content(cont, rAbar, 2, smctx);
2270                 if (!success)
2271                     goto cleanup;
2272                 fq_nmod_mpoly_divides(rAbar, rAbar, cont, smctx);
2273                 fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
2274                 if (fq_nmod_mpoly_divides(rG, A, rAbar, smctx) &&
2275                     fq_nmod_mpoly_divides(rBbar, B, rG, smctx))
2276                 {
2277                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
2278                     fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
2279                     break;
2280                 }
2281             }
2282 
2283             if ((use & USE_BBAR) && !newfq_nmod_mpolyn_interp_mcrt_lg_mpoly(
2284                                  &lastdeg, Bbarn, smctx, modulus, tmp->coeffs,
2285                                                          Bbar, lgctx, cur_emb))
2286             {
2287                 fq_nmod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, smctx);
2288                 success = fq_nmod_mpolyl_content(cont, rBbar, 2, smctx);
2289                 if (!success)
2290                     goto cleanup;
2291                 fq_nmod_mpoly_divides(rBbar, rBbar, cont, smctx);
2292                 fq_nmod_mpoly_repack_bits_inplace(rBbar, bits, smctx);
2293                 if (fq_nmod_mpoly_divides(rG, B, rBbar, smctx) &&
2294                     fq_nmod_mpoly_divides(rAbar, A, rG, smctx))
2295                 {
2296                     fq_nmod_mpoly_repack_bits_inplace(rG, bits, smctx);
2297                     fq_nmod_mpoly_repack_bits_inplace(rAbar, bits, smctx);
2298                     break;
2299                 }
2300             }
2301 
2302             if (n_fq_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
2303                 n_fq_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
2304             {
2305                 goto choose_main;
2306             }
2307 
2308             n_fq_poly_mul_(modulus, modulus, cur_emb->h_as_n_fq_poly,
2309                                                  smctx->fqctx, St->poly_stack);
2310         }
2311     }
2312 
2313     success = 1;
2314 
2315 cleanup:
2316 
2317     for (i = 0; i < nvars; i++)
2318     {
2319         fq_nmod_clear(betas + i, smctx->fqctx);
2320         fq_nmod_clear(alphas + i, smctx->fqctx);
2321     }
2322     flint_free(betas);
2323     flint_free(alphas);
2324 
2325     fq_nmod_mpolyn_clear(An, smctx);
2326     fq_nmod_mpolyn_clear(Bn, smctx);
2327     fq_nmod_mpolyn_clear(gamman, smctx);
2328     bad_fq_nmod_mpoly_embed_chooser_clear(embc, lgctx, smctx, state);
2329 
2330     n_fq_polyun_clear(HG);
2331     n_fq_polyun_clear(HAbar);
2332     n_fq_polyun_clear(HBbar);
2333     n_fq_polyun_clear(MG);
2334     n_fq_polyun_clear(MAbar);
2335     n_fq_polyun_clear(MBbar);
2336     n_fq_polyun_clear(ZG);
2337     n_fq_polyun_clear(ZAbar);
2338     n_fq_polyun_clear(ZBbar);
2339     n_fq_bpoly_clear(Aev);
2340     n_fq_bpoly_clear(Bev);
2341     n_fq_bpoly_clear(Gev);
2342     n_fq_bpoly_clear(Abarev);
2343     n_fq_bpoly_clear(Bbarev);
2344     fq_nmod_mpoly_clear(cont, smctx);
2345     fq_nmod_mpoly_clear(T, smctx);
2346     fq_nmod_mpoly_clear(G, smctx);
2347     fq_nmod_mpoly_clear(Abar, smctx);
2348     fq_nmod_mpoly_clear(Bbar, smctx);
2349     fq_nmod_mpolyn_clear(Tn, smctx);
2350     fq_nmod_mpolyn_clear(Gn, smctx);
2351     fq_nmod_mpolyn_clear(Abarn, smctx);
2352     fq_nmod_mpolyn_clear(Bbarn, smctx);
2353     n_fq_polyun_clear(Aeh_cur);
2354     n_fq_polyun_clear(Aeh_inc);
2355     n_fq_polyun_clear(Beh_cur);
2356     n_fq_polyun_clear(Beh_inc);
2357     n_fq_poly_clear(gammaeh_cur);
2358     n_fq_poly_clear(gammaeh_inc);
2359     n_fq_poly_clear(modulus);
2360     n_poly_stack_clear(St->poly_stack);
2361     n_bpoly_stack_clear(St->bpoly_stack);
2362 
2363     n_poly_clear(tmp);
2364     n_fq_poly_clear(alphapow);
2365     fq_nmod_clear(start_alpha, smctx->fqctx);
2366 
2367     flint_randclear(state);
2368 
2369     for (i = 0; i < nvars; i++)
2370     {
2371         fq_nmod_mpoly_clear(Aevals + i, smctx);
2372         fq_nmod_mpoly_clear(Bevals + i, smctx);
2373         fq_nmod_mpoly_clear(gammaevals + i, smctx);
2374     }
2375     flint_free(Aevals);
2376     flint_free(Bevals);
2377     flint_free(gammaevals);
2378 
2379     FLINT_ASSERT(!success || rG->bits == bits);
2380     FLINT_ASSERT(!success || rAbar->bits == bits);
2381     FLINT_ASSERT(!success || rBbar->bits == bits);
2382 
2383     return success;
2384 
2385 gcd_is_trivial:
2386 
2387     fq_nmod_mpoly_one(rG, smctx);
2388     fq_nmod_mpoly_set(rAbar, A, smctx);
2389     fq_nmod_mpoly_set(rBbar, B, smctx);
2390 
2391     success = 1;
2392 
2393     goto cleanup;
2394 }
2395 
2396 
2397 /* should find its way back here in interesting cases */
fq_nmod_mpoly_gcd_zippel2(fq_nmod_mpoly_t G,const fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,const fq_nmod_mpoly_ctx_t ctx)2398 int fq_nmod_mpoly_gcd_zippel2(
2399     fq_nmod_mpoly_t G,
2400     const fq_nmod_mpoly_t A,
2401     const fq_nmod_mpoly_t B,
2402     const fq_nmod_mpoly_ctx_t ctx)
2403 {
2404     if (fq_nmod_mpoly_is_zero(A, ctx) || fq_nmod_mpoly_is_zero(B, ctx))
2405         return fq_nmod_mpoly_gcd(G, A, B, ctx);
2406 
2407     return _fq_nmod_mpoly_gcd_algo(G, NULL, NULL, A, B, ctx, MPOLY_GCD_USE_ZIPPEL2);
2408 }
2409 
2410