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