1 /*
2 Copyright (C) 2020 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpz_mod_vec.h"
13 #include "fmpz_mod_mpoly.h"
14 #include "fmpz_mod_mpoly_factor.h"
15
16
fmpz_mod_mpolyn_interp_lift_sm_polyu1n(fmpz_mod_mpolyn_t F,fmpz_mod_polyun_t A,const fmpz_mod_mpoly_ctx_t ctx)17 void fmpz_mod_mpolyn_interp_lift_sm_polyu1n(
18 fmpz_mod_mpolyn_t F,
19 fmpz_mod_polyun_t A,
20 const fmpz_mod_mpoly_ctx_t ctx)
21 {
22 slong N = mpoly_words_per_exp_sp(F->bits, ctx->minfo);
23 slong i, j, Fi;
24 slong off0, shift0, off1, shift1;
25
26 mpoly_gen_offset_shift_sp(&off0, &shift0, 0, F->bits, ctx->minfo);
27 mpoly_gen_offset_shift_sp(&off1, &shift1, 1, F->bits, ctx->minfo);
28
29 Fi = 0;
30 for (i = 0; i < A->length; i++)
31 {
32 fmpz * Aicoeffs = A->coeffs[i].coeffs;
33 ulong e0 = A->exps[i] << shift0;
34
35 for (j = A->coeffs[i].length - 1; j >= 0; j--)
36 {
37 if (fmpz_is_zero(Aicoeffs + j))
38 continue;
39
40 fmpz_mod_mpolyn_fit_length(F, Fi + 1, ctx);
41 mpoly_monomial_zero(F->exps + N*Fi, N);
42 (F->exps + N*Fi)[off0] = e0;
43 (F->exps + N*Fi)[off1] += (j << shift1);
44 fmpz_mod_poly_set_fmpz(F->coeffs + Fi, Aicoeffs + j, ctx->ffinfo);
45 Fi++;
46 }
47 }
48
49 F->length = Fi;
50 }
51
fmpz_mod_mpolyn_interp_crt_sm_polyu1n(slong * lastdeg,fmpz_mod_mpolyn_t F,fmpz_mod_mpolyn_t T,fmpz_mod_polyun_t A,fmpz_mod_poly_t modulus,fmpz_mod_poly_t alphapow,const fmpz_mod_mpoly_ctx_t ctx)52 int fmpz_mod_mpolyn_interp_crt_sm_polyu1n(
53 slong * lastdeg,
54 fmpz_mod_mpolyn_t F,
55 fmpz_mod_mpolyn_t T,
56 fmpz_mod_polyun_t A,
57 fmpz_mod_poly_t modulus,
58 fmpz_mod_poly_t alphapow,
59 const fmpz_mod_mpoly_ctx_t ctx)
60 {
61 int changed = 0;
62 slong N = mpoly_words_per_exp(F->bits, ctx->minfo);
63 slong off0, shift0, off1, shift1;
64 fmpz_mod_poly_struct * Acoeffs = A->coeffs;
65 ulong * Aexps = A->exps;
66 slong Fi, Ti, Ai, ai;
67 slong Alen = A->length;
68 slong Flen = F->length;
69 ulong * Fexps = F->exps;
70 fmpz_mod_poly_struct * Fcoeffs = F->coeffs;
71 ulong * Texps = T->exps;
72 fmpz_mod_poly_struct * Tcoeffs = T->coeffs;
73 fmpz_t v;
74 ulong Fexpi, mask;
75
76 fmpz_init(v);
77
78 mask = (-UWORD(1)) >> (FLINT_BITS - F->bits);
79 mpoly_gen_offset_shift_sp(&off0, &shift0, 0, F->bits, ctx->minfo);
80 mpoly_gen_offset_shift_sp(&off1, &shift1, 1, F->bits, ctx->minfo);
81
82 FLINT_ASSERT(T->bits == F->bits);
83
84 *lastdeg = -1;
85
86 Ti = Fi = 0;
87 Ai = 0;
88 ai = 0;
89 if (Ai < Alen)
90 ai = fmpz_mod_poly_degree(Acoeffs + Ai, ctx->ffinfo);
91
92 while (Fi < Flen || Ai < Alen)
93 {
94 if (Ti >= T->alloc)
95 {
96 slong extra = FLINT_MAX(Flen - Fi, Alen - Ai);
97 fmpz_mod_mpolyn_fit_length(T, Ti + extra + 1, ctx);
98 Tcoeffs = T->coeffs;
99 Texps = T->exps;
100 }
101
102 if (Fi < Flen)
103 Fexpi = pack_exp2(((Fexps + N*Fi)[off0]>>shift0)&mask,
104 ((Fexps + N*Fi)[off1]>>shift1)&mask);
105 else
106 Fexpi = 0;
107
108 if (Fi < Flen && Ai < Alen && Fexpi == pack_exp2(Aexps[Ai], ai))
109 {
110 /* F term ok, A term ok */
111 mpoly_monomial_set(Texps + N*Ti, Fexps + N*Fi, N);
112
113 fmpz_mod_poly_eval_pow(v, Fcoeffs + Fi, alphapow, ctx->ffinfo);
114 fmpz_mod_sub(v, Acoeffs[Ai].coeffs + ai, v, ctx->ffinfo);
115 changed |= !fmpz_is_zero(v);
116 fmpz_mod_poly_scalar_addmul_fmpz_mod(Tcoeffs + Ti,
117 Fcoeffs + Fi, modulus, v, ctx->ffinfo);
118 Fi++;
119 do {
120 ai--;
121 } while (ai >= 0 && fmpz_is_zero(Acoeffs[Ai].coeffs + ai));
122 if (ai < 0)
123 {
124 Ai++;
125 if (Ai < Alen)
126 ai = fmpz_mod_poly_degree(Acoeffs + Ai, ctx->ffinfo);
127 }
128 }
129 else if (Ai < Alen && (Fi >= Flen || Fexpi < pack_exp2(Aexps[Ai], ai)))
130 {
131 /* F term missing, A term ok */
132 mpoly_monomial_zero(Texps + N*Ti, N);
133 (Texps + N*Ti)[off0] += (Ai << shift0);
134 (Texps + N*Ti)[off1] += (ai << shift1);
135
136 changed = 1;
137 fmpz_mod_poly_scalar_mul_fmpz(Tcoeffs + Ti, modulus,
138 Acoeffs[Ai].coeffs + ai, ctx->ffinfo);
139
140 do {
141 ai--;
142 } while (ai >= 0 && fmpz_is_zero(Acoeffs[Ai].coeffs + ai));
143 if (ai < 0)
144 {
145 Ai++;
146 if (Ai < Alen)
147 ai = fmpz_mod_poly_degree(Acoeffs + Ai, ctx->ffinfo);
148 }
149 }
150 else
151 {
152 FLINT_ASSERT(Fi < Flen && (Ai >= Alen || Fexpi > pack_exp2(Aexps[Ai], ai)));
153 /* F term ok, Aterm missing */
154 mpoly_monomial_set(Texps + N*Ti, Fexps + N*Fi, N);
155
156 fmpz_mod_poly_eval_pow(v, Fcoeffs + Fi, alphapow, ctx->ffinfo);
157 fmpz_mod_neg(v, v, ctx->ffinfo);
158 changed |= !fmpz_is_zero(v);
159 fmpz_mod_poly_scalar_addmul_fmpz_mod(Tcoeffs + Ti,
160 Fcoeffs + Fi, modulus, v, ctx->ffinfo);
161 Fi++;
162 }
163
164 FLINT_ASSERT(!fmpz_mod_poly_is_zero(Tcoeffs + Ti, ctx->ffinfo));
165 *lastdeg = FLINT_MAX(*lastdeg, fmpz_mod_poly_degree(Tcoeffs + Ti, ctx->ffinfo));
166 Ti++;
167 }
168
169 T->length = Ti;
170
171 if (changed)
172 fmpz_mod_mpolyn_swap(T, F, ctx);
173
174 fmpz_clear(v);
175
176 return changed;
177 }
178
179
_fmpz_mod_mpoly_bidegree(const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_ctx_t ctx)180 static ulong _fmpz_mod_mpoly_bidegree(
181 const fmpz_mod_mpoly_t A,
182 const fmpz_mod_mpoly_ctx_t ctx)
183 {
184 FLINT_ASSERT(A->length > 0);
185 return _mpoly_bidegree(A->exps, A->bits, ctx->minfo);
186 }
187
fmpz_mod_polyun_zip_solve(fmpz_mod_mpoly_t A,fmpz_mod_polyun_t Z,fmpz_mod_polyun_t H,fmpz_mod_polyun_t M,const fmpz_mod_mpoly_ctx_t ctx)188 int fmpz_mod_polyun_zip_solve(
189 fmpz_mod_mpoly_t A,
190 fmpz_mod_polyun_t Z,
191 fmpz_mod_polyun_t H,
192 fmpz_mod_polyun_t M,
193 const fmpz_mod_mpoly_ctx_t ctx)
194 {
195 int success;
196 slong Ai, i, n;
197 fmpz * Acoeffs = A->coeffs;
198 fmpz_mod_poly_t t;
199
200 fmpz_mod_poly_init(t, ctx->ffinfo);
201
202 FLINT_ASSERT(Z->length == H->length);
203 FLINT_ASSERT(Z->length == M->length);
204
205 Ai = 0;
206 for (i = 0; i < H->length; i++)
207 {
208 n = H->coeffs[i].length;
209 FLINT_ASSERT(M->coeffs[i].length == n + 1);
210 FLINT_ASSERT(Z->coeffs[i].length >= n);
211 FLINT_ASSERT(Ai + n <= A->length);
212
213 fmpz_mod_poly_fit_length(t, n, ctx->ffinfo);
214
215 success = _fmpz_mod_zip_vand_solve(Acoeffs + Ai,
216 H->coeffs[i].coeffs, n,
217 Z->coeffs[i].coeffs, Z->coeffs[i].length,
218 M->coeffs[i].coeffs, t->coeffs, ctx->ffinfo);
219 if (success < 1)
220 {
221 fmpz_mod_poly_clear(t, ctx->ffinfo);
222 return success;
223 }
224
225 Ai += n;
226 FLINT_ASSERT(Ai <= A->length);
227
228 }
229
230 FLINT_ASSERT(Ai == A->length);
231
232 fmpz_mod_poly_clear(t, ctx->ffinfo);
233 return 1;
234 }
235
236
237 #define USE_G 1
238 #define USE_ABAR 2
239 #define USE_BBAR 4
240
interp_cost(double degG,double numABgamma,double maxnumci,double totnumci,double degxAB,double degyAB)241 static double interp_cost(
242 double degG,
243 double numABgamma,
244 double maxnumci,
245 double totnumci,
246 double degxAB,
247 double degyAB)
248 {
249 return degG*(degG*totnumci + numABgamma + 0.01*maxnumci*(
250 numABgamma + totnumci + (degxAB*degyAB)*(degxAB*degyAB)));
251 }
252
fmpz_mod_mpoly_gcd_get_use_new(slong rGdeg,slong Adeg,slong Bdeg,slong gammadeg,slong degxAB,slong degyAB,slong numABgamma,const fmpz_mod_polyun_t G,const fmpz_mod_polyun_t Abar,const fmpz_mod_polyun_t Bbar)253 int fmpz_mod_mpoly_gcd_get_use_new(
254 slong rGdeg,
255 slong Adeg,
256 slong Bdeg,
257 slong gammadeg,
258 slong degxAB,
259 slong degyAB,
260 slong numABgamma,
261 const fmpz_mod_polyun_t G,
262 const fmpz_mod_polyun_t Abar,
263 const fmpz_mod_polyun_t Bbar)
264 {
265 int use = 0;
266 slong i, lower = FLINT_MAX(gammadeg, rGdeg);
267 slong upper = gammadeg + FLINT_MIN(FLINT_MIN(Adeg, Bdeg), rGdeg);
268 if (lower <= upper)
269 {
270 slong Gdeg = ((ulong)upper + (ulong)lower)/2;
271 slong maxnumci, totnumci;
272 double Gcost, Abarcost, Bbarcost;
273
274 maxnumci = totnumci = 0;
275 for (i = 0; i < G->length; i++)
276 {
277 maxnumci = FLINT_MAX(maxnumci, G->coeffs[i].length);
278 totnumci += G->coeffs[i].length;
279 }
280 FLINT_ASSERT(Gdeg >= 0);
281 Gcost = interp_cost(Gdeg,
282 numABgamma, maxnumci, totnumci, degxAB, degyAB);
283
284 maxnumci = totnumci = 0;
285 for (i = 0; i < Abar->length; i++)
286 {
287 maxnumci = FLINT_MAX(maxnumci, Abar->coeffs[i].length);
288 totnumci += Abar->coeffs[i].length;
289 }
290 FLINT_ASSERT(gammadeg + Adeg - Gdeg >= 0);
291 Abarcost = interp_cost(gammadeg + Adeg - Gdeg,
292 numABgamma, maxnumci, totnumci, degxAB, degyAB);
293
294 maxnumci = totnumci = 0;
295 for (i = 0; i < Bbar->length; i++)
296 {
297 maxnumci = FLINT_MAX(maxnumci, Bbar->coeffs[i].length);
298 totnumci += Bbar->coeffs[i].length;
299 }
300 FLINT_ASSERT(gammadeg + Bdeg - Gdeg >= 0);
301 Bbarcost = interp_cost(gammadeg + Bdeg - Gdeg,
302 numABgamma, maxnumci, totnumci, degxAB, degyAB);
303
304 if (Gcost <= FLINT_MIN(Abarcost, Bbarcost)*1.125)
305 use |= USE_G;
306
307 if (Abarcost <= FLINT_MIN(Gcost, Bbarcost)*1.125)
308 use |= USE_ABAR;
309
310 if (Bbarcost <= FLINT_MIN(Gcost, Abarcost)*1.125)
311 use |= USE_BBAR;
312 }
313
314 if (use == 0)
315 use = USE_G | USE_ABAR | USE_BBAR;
316
317 return use;
318 }
319
320
fmpz_mod_poly_zip_eval_cur_inc_coeff(fmpz_t eval,fmpz_mod_poly_t cur,const fmpz_mod_poly_t inc,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_ctx_t ctx)321 static void fmpz_mod_poly_zip_eval_cur_inc_coeff(
322 fmpz_t eval,
323 fmpz_mod_poly_t cur,
324 const fmpz_mod_poly_t inc,
325 const fmpz_mod_mpoly_t A,
326 const fmpz_mod_mpoly_ctx_t ctx)
327 {
328 FLINT_ASSERT(A->length == cur->length);
329 FLINT_ASSERT(A->length == inc->length);
330 _fmpz_mod_zip_eval_step(eval, cur->coeffs, inc->coeffs, A->coeffs, A->length, ctx->ffinfo);
331 }
332
333
fmpz_mod_mpoly_mock_eval_coeff(fmpz_mod_polyun_t mock,const fmpz_mod_mpoly_t A,const fmpz_mod_polyun_t Aeh_inc,const fmpz_mod_mpoly_ctx_t ctx)334 void fmpz_mod_mpoly_mock_eval_coeff(
335 fmpz_mod_polyun_t mock,
336 const fmpz_mod_mpoly_t A,
337 const fmpz_mod_polyun_t Aeh_inc,
338 const fmpz_mod_mpoly_ctx_t ctx)
339 {
340 slong i, k;
341
342 if (mock->alloc < Aeh_inc->length)
343 {
344 mock->alloc = FLINT_MAX(Aeh_inc->length, mock->alloc + mock->alloc/2);
345 mock->coeffs = FLINT_ARRAY_REALLOC(mock->coeffs, mock->alloc,
346 fmpz_mod_poly_struct);
347 }
348
349 mock->length = Aeh_inc->length;
350
351 k = 0;
352 for (i = 0; i < Aeh_inc->length; i++)
353 {
354 slong l = Aeh_inc->coeffs[i].length;
355 mock->coeffs[i].coeffs = A->coeffs + k;
356 mock->coeffs[i].alloc = l;
357 mock->coeffs[i].length = l;
358 k += l;
359 }
360
361 FLINT_ASSERT(k == A->length);
362 }
363
364
fmpz_mod_mpoly_monomial_evals(fmpz_mod_poly_t E,const fmpz_mod_mpoly_t A,fmpz_mod_poly_struct * beta_caches,slong start,slong stop,const fmpz_mod_mpoly_ctx_t ctx)365 static void fmpz_mod_mpoly_monomial_evals(
366 fmpz_mod_poly_t E,
367 const fmpz_mod_mpoly_t A,
368 fmpz_mod_poly_struct * beta_caches,
369 slong start,
370 slong stop,
371 const fmpz_mod_mpoly_ctx_t ctx)
372 {
373 mpoly_monomial_evals_fmpz_mod(E, A->exps, A->length, A->bits,
374 beta_caches, start, stop, ctx->minfo, ctx->ffinfo);
375 }
376
fmpz_mod_mpoly_monomial_evals2(fmpz_mod_polyun_t E,const fmpz_mod_mpoly_t A,fmpz_mod_poly_struct * betas,slong m,const fmpz_mod_mpoly_ctx_t ctx,n_poly_t tmark)377 static void fmpz_mod_mpoly_monomial_evals2(
378 fmpz_mod_polyun_t E,
379 const fmpz_mod_mpoly_t A,
380 fmpz_mod_poly_struct * betas,
381 slong m,
382 const fmpz_mod_mpoly_ctx_t ctx,
383 n_poly_t tmark) /* temp space */
384 {
385 mpoly2_fill_marks(&tmark->coeffs, &tmark->length, &tmark->alloc,
386 A->exps, A->length, A->bits, ctx->minfo);
387
388 mpoly2_monomial_evals_fmpz_mod(E, A->exps, A->bits, tmark->coeffs,
389 tmark->length, betas, m, ctx->minfo, ctx->ffinfo);
390 }
391
392
fmpz_mod_polyun_add_zip_must_match(fmpz_mod_polyun_t Z,const fmpz_mod_polyun_t A,slong cur_length)393 static int fmpz_mod_polyun_add_zip_must_match(
394 fmpz_mod_polyun_t Z,
395 const fmpz_mod_polyun_t A,
396 slong cur_length)
397 {
398 slong i, Ai, ai;
399 slong Alen = A->length;
400 ulong * Zexps = Z->exps;
401 fmpz_mod_poly_struct * Zcoeffs = Z->coeffs;
402 ulong * Aexps = A->exps;
403 fmpz_mod_poly_struct * Acoeffs = A->coeffs;
404
405 Ai = 0;
406 ai = 0;
407 if (Ai < Alen)
408 ai = Acoeffs[Ai].length - 1;
409
410 for (i = 0; i < Z->length; i++)
411 {
412 if (Ai < Alen && Zexps[i] == pack_exp2(Aexps[Ai], ai))
413 {
414 /* Z present, A present */
415 fmpz_set(Zcoeffs[i].coeffs + cur_length,
416 Acoeffs[Ai].coeffs + ai);
417 Zcoeffs[i].length = cur_length + 1;
418 do {
419 ai--;
420 } while (ai >= 0 && fmpz_is_zero(Acoeffs[Ai].coeffs + ai));
421 if (ai < 0)
422 {
423 Ai++;
424 if (Ai < Alen)
425 ai = Acoeffs[Ai].length - 1;
426 }
427 }
428 else if (Ai >= Alen || Zexps[i] > pack_exp2(Aexps[Ai], ai))
429 {
430 /* Z present, A missing */
431 fmpz_zero(Zcoeffs[i].coeffs + cur_length);
432 Zcoeffs[i].length = cur_length + 1;
433 }
434 else
435 {
436 /* Z missing, A present */
437 return 0;
438 }
439 }
440
441 return Ai >= Alen;
442 }
443
444
fmpz_mod_polyu2n_zip_eval_cur_inc_coeff(fmpz_mod_polyun_t E,fmpz_mod_polyun_t Acur,const fmpz_mod_polyun_t Ainc,const fmpz_mod_polyun_t Acoeff,const fmpz_mod_ctx_t ctx)445 void fmpz_mod_polyu2n_zip_eval_cur_inc_coeff(
446 fmpz_mod_polyun_t E,
447 fmpz_mod_polyun_t Acur,
448 const fmpz_mod_polyun_t Ainc,
449 const fmpz_mod_polyun_t Acoeff,
450 const fmpz_mod_ctx_t ctx)
451 {
452 slong i, Ei;
453 slong e0, e1;
454 fmpz_t c;
455
456 FLINT_ASSERT(Acur->length > 0);
457 FLINT_ASSERT(Acur->length == Ainc->length);
458 FLINT_ASSERT(Acur->length == Acoeff->length);
459
460 fmpz_init(c);
461
462 e0 = extract_exp(Acur->exps[0], 1, 2);
463 e1 = extract_exp(Acur->exps[0], 0, 2);
464
465 fmpz_mod_polyun_fit_length(E, 4, ctx);
466 Ei = 0;
467 E->exps[Ei] = e1;
468 fmpz_mod_poly_zero(E->coeffs + Ei, ctx);
469
470 for (i = 0; i < Acur->length; i++)
471 {
472 slong this_len = Acur->coeffs[i].length;
473 FLINT_ASSERT(this_len == Ainc->coeffs[i].length);
474 FLINT_ASSERT(this_len == Acoeff->coeffs[i].length);
475
476 _fmpz_mod_zip_eval_step(c, Acur->coeffs[i].coeffs,
477 Ainc->coeffs[i].coeffs, Acoeff->coeffs[i].coeffs, this_len, ctx);
478
479 e0 = extract_exp(Acur->exps[i], 1, 2);
480 e1 = extract_exp(Acur->exps[i], 0, 2);
481
482 if (E->exps[Ei] != e0)
483 {
484 fmpz_mod_polyun_fit_length(E, Ei + 2, ctx);
485 Ei += !fmpz_mod_poly_is_zero(E->coeffs + Ei, ctx);
486 E->exps[Ei] = e0;
487 fmpz_mod_poly_zero(E->coeffs + Ei, ctx);
488 }
489
490 fmpz_mod_poly_set_coeff_fmpz(E->coeffs + Ei, e1, c, ctx);
491 }
492
493 Ei += !fmpz_mod_poly_is_zero(E->coeffs + Ei, ctx);
494 E->length = Ei;
495
496 FLINT_ASSERT(fmpz_mod_polyun_is_canonical(E, ctx));
497
498 fmpz_clear(c);
499 }
500
501
502 /*
503 gamma = gcd(lc(A), lc(B))
504
505 fake answers G, Abar, Bbar with
506
507 G*Abar = gamma*A
508 G*Bbar = gamma*B
509 lc(G) = gamma
510 lc(Abar) = lc(A)
511 lc(Bbar) = lc(B)
512
513 real answers
514
515 rG = pp(G)
516 rAbar = Abar/lc(rG)
517 rBbar = Bbar/lc(rG)
518
519 The degrees of A, B, and gamma wrt the minor vars must be passed in.
520 A guess of the degrees of rG wrt the minor vars can be passed in.
521
522
523 deg(G) = deg(gamma) - deg(lc(rG)) + deg(rG)
524 deg(G) <= deg(gamma) + deg(rG)
525 deg(G) <= deg(gamma) + deg(A)
526 deg(G) <= deg(gamma) + deg(B)
527 deg(G) >= deg(gamma)
528 deg(G) >= deg(rG)
529
530 deg(A) = deg(gamma) + deg(A) - deg(G)
531 deg(B) = deg(gamma) + deg(B) - deg(G)
532 */
fmpz_mod_mpolyl_gcd_zippel2_smprime(fmpz_mod_mpoly_t rG,const slong * rGdegs,fmpz_mod_mpoly_t rAbar,fmpz_mod_mpoly_t rBbar,const fmpz_mod_mpoly_t A,const slong * Adegs,const fmpz_mod_mpoly_t B,const slong * Bdegs,const fmpz_mod_mpoly_t gamma,const slong * gammadegs,const fmpz_mod_mpoly_ctx_t ctx)533 int fmpz_mod_mpolyl_gcd_zippel2_smprime(
534 fmpz_mod_mpoly_t rG, const slong * rGdegs, /* guess at rG degrees, could be NULL */
535 fmpz_mod_mpoly_t rAbar,
536 fmpz_mod_mpoly_t rBbar,
537 const fmpz_mod_mpoly_t A, const slong * Adegs,
538 const fmpz_mod_mpoly_t B, const slong * Bdegs,
539 const fmpz_mod_mpoly_t gamma, const slong * gammadegs,
540 const fmpz_mod_mpoly_ctx_t ctx)
541 {
542 int success, use, alpha_tries_left;
543 slong i, m;
544 slong nvars = ctx->minfo->nvars;
545 flint_bitcnt_t bits = A->bits;
546 fmpz * alphas, * betas;
547 fmpz_mod_poly_struct * alpha_caches, * beta_caches;
548 flint_rand_t state;
549 n_poly_t tmark;
550 fmpz_mod_mpoly_t cont;
551 fmpz_mod_mpoly_t T, G, Abar, Bbar;
552 fmpz_mod_polyun_t HG, HAbar, HBbar, MG, MAbar, MBbar, ZG, ZAbar, ZBbar;
553 fmpz_mod_polyun_t Aev, Bev, Gev, Abarev, Bbarev;
554 fmpz_t gammaev;
555 fmpz_mod_mpolyn_t Tn, Gn, Abarn, Bbarn;
556 slong lastdeg;
557 slong cur_zip_image, req_zip_images, this_length;
558 fmpz_mod_polyun_t Aeh_coeff_mock, Beh_coeff_mock;
559 fmpz_mod_polyun_t Aeh_cur, Aeh_inc, Beh_cur, Beh_inc;
560 fmpz_mod_poly_t gammaeh_cur, gammaeh_inc;
561 fmpz_mod_poly_t modulus, alphapow;
562 fmpz_mod_mpoly_struct * Aevals, * Bevals;
563 fmpz_mod_mpoly_struct * gammaevals;
564 fmpz_mod_poly_polyun_stack_t St;
565 fmpz_t c, start_alpha;
566 ulong GdegboundXY, newdegXY, Abideg, Bbideg;
567 slong degxAB, degyAB;
568
569 FLINT_ASSERT(bits <= FLINT_BITS);
570 FLINT_ASSERT(bits == A->bits);
571 FLINT_ASSERT(bits == B->bits);
572 FLINT_ASSERT(bits == gamma->bits);
573 FLINT_ASSERT(A->length > 0);
574 FLINT_ASSERT(B->length > 0);
575 FLINT_ASSERT(gamma->length > 0);
576
577 fmpz_mod_mpoly_fit_length_reset_bits(rG, 1, bits, ctx);
578 fmpz_mod_mpoly_fit_length_reset_bits(rAbar, 1, bits, ctx);
579 fmpz_mod_mpoly_fit_length_reset_bits(rBbar, 1, bits, ctx);
580
581 fmpz_init(gammaev);
582 fmpz_init(c);
583 fmpz_init(start_alpha);
584
585 #if FLINT_WANT_ASSERT
586 {
587 slong * tmp_degs = FLINT_ARRAY_ALLOC(nvars, slong);
588
589 fmpz_mod_mpoly_degrees_si(tmp_degs, A, ctx);
590 for (i = 0; i < nvars; i++)
591 FLINT_ASSERT(tmp_degs[i] == Adegs[i]);
592
593 fmpz_mod_mpoly_degrees_si(tmp_degs, B, ctx);
594 for (i = 0; i < nvars; i++)
595 FLINT_ASSERT(tmp_degs[i] == Bdegs[i]);
596
597 fmpz_mod_mpoly_degrees_si(tmp_degs, gamma, ctx);
598 for (i = 0; i < nvars; i++)
599 FLINT_ASSERT(tmp_degs[i] == gammadegs[i]);
600
601 flint_free(tmp_degs);
602 }
603 #endif
604
605 FLINT_ASSERT(gammadegs[0] == 0);
606 FLINT_ASSERT(gammadegs[1] == 0);
607
608 n_poly_init(tmark);
609
610 fmpz_mod_polyun_init(HG, ctx->ffinfo);
611 fmpz_mod_polyun_init(HAbar, ctx->ffinfo);
612 fmpz_mod_polyun_init(HBbar, ctx->ffinfo);
613 fmpz_mod_polyun_init(MG, ctx->ffinfo);
614 fmpz_mod_polyun_init(MAbar, ctx->ffinfo);
615 fmpz_mod_polyun_init(MBbar, ctx->ffinfo);
616 fmpz_mod_polyun_init(ZG, ctx->ffinfo);
617 fmpz_mod_polyun_init(ZAbar, ctx->ffinfo);
618 fmpz_mod_polyun_init(ZBbar, ctx->ffinfo);
619 fmpz_mod_polyun_init(Aev, ctx->ffinfo);
620 fmpz_mod_polyun_init(Bev, ctx->ffinfo);
621 fmpz_mod_polyun_init(Gev, ctx->ffinfo);
622 fmpz_mod_polyun_init(Abarev, ctx->ffinfo);
623 fmpz_mod_polyun_init(Bbarev, ctx->ffinfo);
624 fmpz_mod_poly_init2(alphapow, 4, ctx->ffinfo);
625 fmpz_mod_mpoly_init3(cont, 1, bits, ctx);
626 fmpz_mod_mpoly_init3(T, 1, bits, ctx);
627 fmpz_mod_mpoly_init3(G, 1, bits, ctx);
628 fmpz_mod_mpoly_init3(Abar, 1, bits, ctx);
629 fmpz_mod_mpoly_init3(Bbar, 1, bits, ctx);
630 fmpz_mod_mpolyn_init(Tn, bits, ctx);
631 fmpz_mod_mpolyn_init(Gn, bits, ctx);
632 fmpz_mod_mpolyn_init(Abarn, bits, ctx);
633 fmpz_mod_mpolyn_init(Bbarn, bits, ctx);
634 fmpz_mod_polyun_init(Aeh_cur, ctx->ffinfo);
635 fmpz_mod_polyun_init(Aeh_inc, ctx->ffinfo);
636 fmpz_mod_polyun_init(Beh_cur, ctx->ffinfo);
637 fmpz_mod_polyun_init(Beh_inc, ctx->ffinfo);
638 fmpz_mod_poly_init(gammaeh_cur, ctx->ffinfo);
639 fmpz_mod_poly_init(gammaeh_inc, ctx->ffinfo);
640 fmpz_mod_poly_init(modulus, ctx->ffinfo);
641 fmpz_mod_poly_stack_init(St->poly_stack);
642 fmpz_mod_polyun_stack_init(St->polyun_stack);
643
644 Aeh_coeff_mock->exps = NULL;
645 Aeh_coeff_mock->coeffs = NULL;
646 Aeh_coeff_mock->length = 0;
647 Aeh_coeff_mock->alloc = 0;
648
649 Beh_coeff_mock->exps = NULL;
650 Beh_coeff_mock->coeffs = NULL;
651 Beh_coeff_mock->length = 0;
652 Beh_coeff_mock->alloc = 0;
653
654 betas = _fmpz_vec_init(nvars);
655 alphas = _fmpz_vec_init(nvars);
656 flint_randinit(state);
657
658 beta_caches = FLINT_ARRAY_ALLOC(nvars, fmpz_mod_poly_struct);
659 alpha_caches = FLINT_ARRAY_ALLOC(nvars, fmpz_mod_poly_struct);
660 for (i = 0; i < nvars; i++)
661 {
662 fmpz_mod_poly_init(beta_caches + i, ctx->ffinfo);
663 fmpz_mod_poly_init(alpha_caches + i, ctx->ffinfo);
664 }
665
666 Aevals = FLINT_ARRAY_ALLOC(nvars + 1, fmpz_mod_mpoly_struct);
667 Bevals = FLINT_ARRAY_ALLOC(nvars + 1, fmpz_mod_mpoly_struct);
668 gammaevals = FLINT_ARRAY_ALLOC(nvars + 1, fmpz_mod_mpoly_struct);
669 for (i = 0; i < nvars; i++)
670 {
671 fmpz_mod_mpoly_init3(Aevals + i, 0, bits, ctx);
672 fmpz_mod_mpoly_init3(Bevals + i, 0, bits, ctx);
673 fmpz_mod_mpoly_init3(gammaevals + i, 0, bits, ctx);
674 }
675 Aevals[nvars] = *A;
676 Bevals[nvars] = *B;
677 gammaevals[nvars] = *gamma;
678
679 Abideg = _fmpz_mod_mpoly_bidegree(A, ctx);
680 Bbideg = _fmpz_mod_mpoly_bidegree(B, ctx);
681
682 degxAB = FLINT_MAX(Adegs[0], Bdegs[0]);
683 degyAB = FLINT_MAX(Adegs[1], Bdegs[1]);
684
685 GdegboundXY = pack_exp2(FLINT_MIN(Adegs[0], Bdegs[0]),
686 FLINT_MIN(Adegs[1], Bdegs[1]));
687 if (GdegboundXY == 0)
688 goto gcd_is_trivial;
689
690 alpha_tries_left = 20;
691
692 choose_alphas:
693
694 if (--alpha_tries_left < 0)
695 {
696 success = 0;
697 goto cleanup;
698 }
699
700 for (i = 2; i < nvars; i++)
701 fmpz_mod_rand_not_zero(alphas + i, state, ctx->ffinfo);
702
703 for (i = nvars - 1; i >= 2; i--)
704 {
705 fmpz_mod_mpoly_evaluate_one_fmpz(Aevals + i, Aevals + i + 1, i, alphas + i, ctx);
706 fmpz_mod_mpoly_evaluate_one_fmpz(Bevals + i, Bevals + i + 1, i, alphas + i, ctx);
707 fmpz_mod_mpoly_evaluate_one_fmpz(gammaevals + i, gammaevals + i + 1, i, alphas + i, ctx);
708 if (fmpz_mod_mpoly_is_zero(gammaevals + i, ctx))
709 goto choose_alphas;
710 if (Aevals[i].length < 1 || _fmpz_mod_mpoly_bidegree(Aevals + i, ctx) != Abideg)
711 goto choose_alphas;
712 if (Bevals[i].length < 1 || _fmpz_mod_mpoly_bidegree(Bevals + i, ctx) != Bbideg)
713 goto choose_alphas;
714 }
715
716 m = 2;
717
718 if (rGdegs == NULL)
719 use = USE_G | USE_ABAR | USE_BBAR;
720 else
721 use = mpoly_gcd_get_use_first(rGdegs[m], Adegs[m], Bdegs[m], gammadegs[m]);
722
723 fmpz_mod_mpoly_get_polyu1n(Aev, Aevals + m, 0, 1, ctx);
724 fmpz_mod_mpoly_get_polyu1n(Bev, Bevals + m, 0, 1, ctx);
725
726 success = fmpz_mod_polyu1n_gcd_brown_smprime(Gev, Abarev, Bbarev, Aev, Bev,
727 ctx->ffinfo, St->poly_stack, St->polyun_stack);
728 if (!success)
729 goto cleanup;
730
731 newdegXY = fmpz_mod_polyu1n_bidegree(Gev);
732 if (newdegXY > GdegboundXY)
733 goto choose_alphas;
734 if (newdegXY < GdegboundXY)
735 {
736 GdegboundXY = newdegXY;
737 if (GdegboundXY == 0)
738 goto gcd_is_trivial;
739 }
740
741 fmpz_mod_mpoly_get_fmpz(gammaev, gammaevals + m, ctx);
742 _fmpz_mod_poly_vec_mul_fmpz_mod(Gev->coeffs, Gev->length, gammaev, ctx->ffinfo);
743
744 fmpz_mod_mpolyn_interp_lift_sm_polyu1n(Gn, Gev, ctx);
745 fmpz_mod_mpolyn_interp_lift_sm_polyu1n(Abarn, Abarev, ctx);
746 fmpz_mod_mpolyn_interp_lift_sm_polyu1n(Bbarn, Bbarev, ctx);
747
748 fmpz_mod_poly_one(modulus, ctx->ffinfo);
749 fmpz_mod_neg(c, alphas + m, ctx->ffinfo);
750 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 1, c, ctx->ffinfo);
751
752 fmpz_set(start_alpha, alphas + m);
753 while (1)
754 {
755 choose_alpha_2:
756
757 if (fmpz_cmp_ui(alphas + m, 2) < 0)
758 fmpz_set(alphas + m, fmpz_mod_ctx_modulus(ctx->ffinfo));
759 fmpz_sub_ui(alphas + m, alphas + m, 1);
760 if (fmpz_equal(alphas + m, start_alpha))
761 goto choose_alphas;
762
763 FLINT_ASSERT(alphapow->alloc >= 2);
764 fmpz_one(alphapow->coeffs + 0);
765 fmpz_set(alphapow->coeffs + 1, alphas + m);
766 alphapow->length = 2;
767
768 fmpz_mod_mpoly_evaluate_one_fmpz(Aevals + m, Aevals + m + 1, m, alphas + m, ctx);
769 fmpz_mod_mpoly_evaluate_one_fmpz(Bevals + m, Bevals + m + 1, m, alphas + m, ctx);
770
771 fmpz_mod_mpoly_evaluate_one_fmpz(gammaevals + m, gammaevals + m + 1, m, alphas + m, ctx);
772 if (fmpz_mod_mpoly_is_zero(gammaevals + m, ctx))
773 goto choose_alpha_2;
774 if (Aevals[m].length < 1 || _fmpz_mod_mpoly_bidegree(Aevals + m, ctx) != Abideg)
775 goto choose_alpha_2;
776 if (Bevals[m].length < 1 || _fmpz_mod_mpoly_bidegree(Bevals + m, ctx) != Bbideg)
777 goto choose_alpha_2;
778
779 fmpz_mod_mpoly_get_polyu1n(Aev, Aevals + m, 0, 1, ctx);
780 fmpz_mod_mpoly_get_polyu1n(Bev, Bevals + m, 0, 1, ctx);
781
782 success = fmpz_mod_polyu1n_gcd_brown_smprime(Gev, Abarev, Bbarev,
783 Aev, Bev, ctx->ffinfo, St->poly_stack, St->polyun_stack);
784 if (!success)
785 goto cleanup;
786
787 newdegXY = fmpz_mod_polyu1n_bidegree(Gev);
788 if (newdegXY > GdegboundXY)
789 goto choose_alpha_2;
790 if (newdegXY < GdegboundXY)
791 {
792 GdegboundXY = newdegXY;
793 if (GdegboundXY == 0)
794 goto gcd_is_trivial;
795 goto choose_alphas;
796 }
797
798 fmpz_mod_mpoly_get_fmpz(gammaev, gammaevals + m, ctx);
799 _fmpz_mod_poly_vec_mul_fmpz_mod(Gev->coeffs, Gev->length, gammaev, ctx->ffinfo);
800
801 fmpz_mod_poly_eval_pow(c, modulus, alphapow, ctx->ffinfo);
802 fmpz_mod_inv(c, c, ctx->ffinfo);
803 fmpz_mod_poly_scalar_mul_fmpz(modulus, modulus, c, ctx->ffinfo);
804
805 if ((use & USE_G) && !fmpz_mod_mpolyn_interp_crt_sm_polyu1n(
806 &lastdeg, Gn, Tn, Gev, modulus, alphapow, ctx))
807 {
808 if (m == nvars - 1)
809 {
810 fmpz_mod_mpoly_cvtfrom_mpolyn(rG, Gn, m, ctx);
811 success = fmpz_mod_mpolyl_content(cont, rG, 2, ctx);
812 if (!success)
813 goto cleanup;
814 fmpz_mod_mpoly_divides(rG, rG, cont, ctx);
815 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
816 if (fmpz_mod_mpoly_divides(rAbar, A, rG, ctx) &&
817 fmpz_mod_mpoly_divides(rBbar, B, rG, ctx))
818 {
819 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
820 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
821 break;
822 }
823 }
824 else
825 {
826 fmpz_mod_mpoly_cvtfrom_mpolyn(G, Gn, m, ctx);
827 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
828 if (fmpz_mod_mpoly_divides(Abar, T, G, ctx))
829 {
830 fmpz_mod_mpoly_repack_bits_inplace(Abar, bits, ctx);
831 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
832 if (fmpz_mod_mpoly_divides(Bbar, T, G, ctx))
833 {
834 fmpz_mod_mpoly_repack_bits_inplace(Bbar, bits, ctx);
835 break;
836 }
837 }
838 }
839 }
840
841 if ((use & USE_ABAR) && !fmpz_mod_mpolyn_interp_crt_sm_polyu1n(
842 &lastdeg, Abarn, Tn, Abarev, modulus, alphapow, ctx))
843 {
844 if (m == nvars - 1)
845 {
846 fmpz_mod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, ctx);
847 success = fmpz_mod_mpolyl_content(cont, rAbar, 2, ctx);
848 if (!success)
849 goto cleanup;
850 fmpz_mod_mpoly_divides(rAbar, rAbar, cont, ctx);
851 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
852 if (fmpz_mod_mpoly_divides(rG, A, rAbar, ctx) &&
853 fmpz_mod_mpoly_divides(rBbar, B, rG, ctx))
854 {
855 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
856 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
857 break;
858 }
859 }
860 else
861 {
862 fmpz_mod_mpoly_cvtfrom_mpolyn(Abar, Abarn, m, ctx);
863 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
864 if (fmpz_mod_mpoly_divides(G, T, Abar, ctx))
865 {
866 fmpz_mod_mpoly_repack_bits_inplace(G, bits, ctx);
867 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
868 if (fmpz_mod_mpoly_divides(Bbar, T, G, ctx))
869 {
870 fmpz_mod_mpoly_repack_bits_inplace(Bbar, bits, ctx);
871 break;
872 }
873 }
874 }
875 }
876
877 if ((use & USE_BBAR) && !fmpz_mod_mpolyn_interp_crt_sm_polyu1n(
878 &lastdeg, Bbarn, Tn, Bbarev, modulus, alphapow, ctx))
879 {
880 if (m == nvars - 1)
881 {
882 fmpz_mod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, ctx);
883 success = fmpz_mod_mpolyl_content(cont, rBbar, 2, ctx);
884 if (!success)
885 goto cleanup;
886 fmpz_mod_mpoly_divides(rBbar, rBbar, cont, ctx);
887 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
888 if (fmpz_mod_mpoly_divides(rG, B, rBbar, ctx) &&
889 fmpz_mod_mpoly_divides(rAbar, A, rG, ctx))
890 {
891 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
892 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
893 break;
894 }
895 }
896 else
897 {
898 fmpz_mod_mpoly_cvtfrom_mpolyn(Bbar, Bbarn, m, ctx);
899 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
900 if (fmpz_mod_mpoly_divides(G, T, Bbar, ctx))
901 {
902 fmpz_mod_mpoly_repack_bits_inplace(G, bits, ctx);
903 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
904 if (fmpz_mod_mpoly_divides(Abar, T, G, ctx))
905 {
906 fmpz_mod_mpoly_repack_bits_inplace(Abar, bits, ctx);
907 break;
908 }
909 }
910 }
911 }
912
913 if (fmpz_mod_poly_degree(modulus, ctx->ffinfo) > gammadegs[m] + Adegs[m] &&
914 fmpz_mod_poly_degree(modulus, ctx->ffinfo) > gammadegs[m] + Bdegs[m])
915 {
916 goto choose_alphas;
917 }
918
919 fmpz_mod_neg(c, alphas + m, ctx->ffinfo);
920 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 1, c, ctx->ffinfo);
921 }
922
923 for (m = 3; m < nvars; m++)
924 {
925 /* G, Abar, Bbar are in Fp[gen(0), ..., gen(m - 1)] */
926 fmpz_mod_mpolyn_interp_lift_sm_mpoly(Gn, G, ctx);
927 fmpz_mod_mpolyn_interp_lift_sm_mpoly(Abarn, Abar, ctx);
928 fmpz_mod_mpolyn_interp_lift_sm_mpoly(Bbarn, Bbar, ctx);
929
930 if (rGdegs == NULL)
931 use = USE_G | USE_ABAR | USE_BBAR;
932 else
933 use = 0;
934
935 fmpz_mod_poly_one(modulus, ctx->ffinfo);
936 fmpz_mod_neg(c, alphas + m, ctx->ffinfo);
937 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 1, c, ctx->ffinfo);
938
939 fmpz_set(start_alpha, alphas + m);
940
941 choose_betas:
942
943 /* only beta[2], beta[1], ..., beta[m - 1] will be used */
944 for (i = 2; i < nvars; i++)
945 {
946 fmpz_mod_rand_not_zero(betas + i, state, ctx->ffinfo);
947 fmpz_mod_pow_cache_start(betas + i, beta_caches + i, ctx->ffinfo);
948 }
949
950 fmpz_mod_mpoly_monomial_evals2(HG, G, beta_caches + 2, m, ctx, tmark);
951 fmpz_mod_mpoly_monomial_evals2(HAbar, Abar, beta_caches + 2, m, ctx, tmark);
952 fmpz_mod_mpoly_monomial_evals2(HBbar, Bbar, beta_caches + 2, m, ctx, tmark);
953
954 if (use == 0)
955 {
956 this_length = gammaevals[m + 1].length + Aevals[m + 1].length +
957 Bevals[m + 1].length;
958
959 use = fmpz_mod_mpoly_gcd_get_use_new(rGdegs[m], Adegs[m], Bdegs[m],
960 gammadegs[m], degxAB, degyAB, this_length, HG, HAbar, HBbar);
961 }
962
963 req_zip_images = 1;
964 if (use & USE_G)
965 {
966 this_length = fmpz_mod_polyun_product_roots(MG, HG, ctx->ffinfo);
967 req_zip_images = FLINT_MAX(req_zip_images, this_length);
968 }
969 if (use & USE_ABAR)
970 {
971 this_length = fmpz_mod_polyun_product_roots(MAbar, HAbar, ctx->ffinfo);
972 req_zip_images = FLINT_MAX(req_zip_images, this_length);
973 }
974 if (use & USE_BBAR)
975 {
976 this_length = fmpz_mod_polyun_product_roots(MBbar, HBbar, ctx->ffinfo);
977 req_zip_images = FLINT_MAX(req_zip_images, this_length);
978 }
979
980 while (1)
981 {
982 choose_alpha_m:
983
984 if (fmpz_cmp_ui(alphas + m, 2) < 0)
985 fmpz_set(alphas + m, fmpz_mod_ctx_modulus(ctx->ffinfo));
986 fmpz_sub_ui(alphas + m, alphas + m, 1);
987 if (fmpz_equal(alphas + m, start_alpha))
988 goto choose_alphas;
989
990 FLINT_ASSERT(alphapow->alloc >= 2);
991 fmpz_one(alphapow->coeffs + 0);
992 fmpz_set(alphapow->coeffs + 1, alphas + m);
993 alphapow->length = 2;
994
995 fmpz_mod_mpoly_evaluate_one_fmpz(Aevals + m, Aevals + m + 1, m, alphas + m, ctx);
996 fmpz_mod_mpoly_evaluate_one_fmpz(Bevals + m, Bevals + m + 1, m, alphas + m, ctx);
997 fmpz_mod_mpoly_evaluate_one_fmpz(gammaevals + m, gammaevals + m + 1, m, alphas + m, ctx);
998 if (fmpz_mod_mpoly_is_zero(gammaevals + m, ctx))
999 goto choose_alpha_m;
1000 if (Aevals[m].length < 1 || _fmpz_mod_mpoly_bidegree(Aevals + m, ctx) != Abideg)
1001 goto choose_alpha_m;
1002 if (Bevals[m].length < 1 || _fmpz_mod_mpoly_bidegree(Bevals + m, ctx) != Bbideg)
1003 goto choose_alpha_m;
1004
1005 fmpz_mod_mpoly_monomial_evals2(Aeh_inc, Aevals + m, beta_caches + 2, m, ctx, tmark);
1006 fmpz_mod_mpoly_monomial_evals2(Beh_inc, Bevals + m, beta_caches + 2, m, ctx, tmark);
1007 fmpz_mod_mpoly_monomial_evals(gammaeh_inc, gammaevals + m, beta_caches + 2, 2, m, ctx);
1008 fmpz_mod_polyun_set(Aeh_cur, Aeh_inc, ctx->ffinfo);
1009 fmpz_mod_polyun_set(Beh_cur, Beh_inc, ctx->ffinfo);
1010 fmpz_mod_poly_set(gammaeh_cur, gammaeh_inc, ctx->ffinfo);
1011 fmpz_mod_mpoly_mock_eval_coeff(Aeh_coeff_mock, Aevals + m, Aeh_inc, ctx);
1012 fmpz_mod_mpoly_mock_eval_coeff(Beh_coeff_mock, Bevals + m, Beh_inc, ctx);
1013
1014 fmpz_mod_polyun_zip_start(ZG, HG, req_zip_images, ctx->ffinfo);
1015 fmpz_mod_polyun_zip_start(ZAbar, HAbar, req_zip_images, ctx->ffinfo);
1016 fmpz_mod_polyun_zip_start(ZBbar, HBbar, req_zip_images, ctx->ffinfo);
1017 for (cur_zip_image = 0; cur_zip_image < req_zip_images; cur_zip_image++)
1018 {
1019 fmpz_mod_polyu2n_zip_eval_cur_inc_coeff(Aev, Aeh_cur, Aeh_inc, Aeh_coeff_mock, ctx->ffinfo);
1020 fmpz_mod_polyu2n_zip_eval_cur_inc_coeff(Bev, Beh_cur, Beh_inc, Beh_coeff_mock, ctx->ffinfo);
1021 fmpz_mod_poly_zip_eval_cur_inc_coeff(gammaev, gammaeh_cur, gammaeh_inc, gammaevals + m, ctx);
1022 if (fmpz_is_zero(gammaev))
1023 goto choose_betas;
1024 if (Aev->length < 1 || fmpz_mod_polyu1n_bidegree(Aev) != Abideg)
1025 goto choose_betas;
1026 if (Bev->length < 1 || fmpz_mod_polyu1n_bidegree(Bev) != Bbideg)
1027 goto choose_betas;
1028
1029 success = fmpz_mod_polyu1n_gcd_brown_smprime(Gev, Abarev, Bbarev,
1030 Aev, Bev, ctx->ffinfo, St->poly_stack, St->polyun_stack);
1031 if (!success)
1032 goto cleanup;
1033
1034 newdegXY = fmpz_mod_polyu1n_bidegree(Gev);
1035 if (newdegXY > GdegboundXY)
1036 goto choose_betas;
1037 if (newdegXY < GdegboundXY)
1038 {
1039 GdegboundXY = newdegXY;
1040 if (GdegboundXY == 0)
1041 goto gcd_is_trivial;
1042 goto choose_alphas;
1043 }
1044
1045 _fmpz_mod_poly_vec_mul_fmpz_mod(Gev->coeffs, Gev->length, gammaev, ctx->ffinfo);
1046 if ((use & USE_G) && !fmpz_mod_polyun_add_zip_must_match(ZG, Gev, cur_zip_image))
1047 goto choose_alphas;
1048 if ((use & USE_ABAR) && !fmpz_mod_polyun_add_zip_must_match(ZAbar, Abarev, cur_zip_image))
1049 goto choose_alphas;
1050 if ((use & USE_BBAR) && !fmpz_mod_polyun_add_zip_must_match(ZBbar, Bbarev, cur_zip_image))
1051 goto choose_alphas;
1052 }
1053
1054 if ((use & USE_G) && fmpz_mod_polyun_zip_solve(G, ZG, HG, MG, ctx) < 1)
1055 goto choose_alphas;
1056 if ((use & USE_ABAR) && fmpz_mod_polyun_zip_solve(Abar, ZAbar, HAbar, MAbar, ctx) < 1)
1057 goto choose_alphas;
1058 if ((use & USE_BBAR) && fmpz_mod_polyun_zip_solve(Bbar, ZBbar, HBbar, MBbar, ctx) < 1)
1059 goto choose_alphas;
1060
1061 FLINT_ASSERT(fmpz_mod_poly_degree(modulus, ctx->ffinfo) > 0);
1062 fmpz_mod_poly_eval_pow(c, modulus, alphapow, ctx->ffinfo);
1063 fmpz_mod_inv(c, c, ctx->ffinfo);
1064 fmpz_mod_poly_scalar_mul_fmpz(modulus, modulus, c, ctx->ffinfo);
1065
1066 if ((use & USE_G) && !fmpz_mod_mpolyn_interp_mcrt_sm_mpoly(
1067 &lastdeg, Gn, G, modulus, alphapow, ctx))
1068 {
1069 fmpz_mod_mpoly_cvtfrom_mpolyn(rG, Gn, m, ctx);
1070 if (m == nvars - 1)
1071 {
1072 success = fmpz_mod_mpolyl_content(cont, rG, 2, ctx);
1073 if (!success)
1074 goto cleanup;
1075 fmpz_mod_mpoly_divides(rG, rG, cont, ctx);
1076 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
1077 if (fmpz_mod_mpoly_divides(rAbar, A, rG, ctx) &&
1078 fmpz_mod_mpoly_divides(rBbar, B, rG, ctx))
1079 {
1080 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1081 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1082 break;
1083 }
1084 }
1085 else
1086 {
1087 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1088 if (fmpz_mod_mpoly_divides(rAbar, T, rG, ctx))
1089 {
1090 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1091 if (fmpz_mod_mpoly_divides(rBbar, T, rG, ctx))
1092 {
1093 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1094 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1095 fmpz_mod_mpoly_swap(G, rG, ctx);
1096 fmpz_mod_mpoly_swap(Abar, rAbar, ctx);
1097 fmpz_mod_mpoly_swap(Bbar, rBbar, ctx);
1098 break;
1099 }
1100 }
1101 }
1102 }
1103 if ((use & USE_ABAR) && !fmpz_mod_mpolyn_interp_mcrt_sm_mpoly(
1104 &lastdeg, Abarn, Abar, modulus, alphapow, ctx))
1105 {
1106 fmpz_mod_mpoly_cvtfrom_mpolyn(rAbar, Abarn, m, ctx);
1107 if (m == nvars - 1)
1108 {
1109 success = fmpz_mod_mpolyl_content(cont, rAbar, 2, ctx);
1110 if (!success)
1111 goto cleanup;
1112 fmpz_mod_mpoly_divides(rAbar, rAbar, cont, ctx);
1113 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1114 if (fmpz_mod_mpoly_divides(rG, A, rAbar, ctx) &&
1115 fmpz_mod_mpoly_divides(rBbar, B, rG, ctx))
1116 {
1117 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
1118 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1119 break;
1120 }
1121 }
1122 else
1123 {
1124 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1125 if (fmpz_mod_mpoly_divides(rG, T, rAbar, ctx))
1126 {
1127 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1128 if (fmpz_mod_mpoly_divides(rBbar, T, rG, ctx))
1129 {
1130 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
1131 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1132 fmpz_mod_mpoly_swap(G, rG, ctx);
1133 fmpz_mod_mpoly_swap(Abar, rAbar, ctx);
1134 fmpz_mod_mpoly_swap(Bbar, rBbar, ctx);
1135 break;
1136 }
1137 }
1138 }
1139 }
1140 if ((use & USE_BBAR) && !fmpz_mod_mpolyn_interp_mcrt_sm_mpoly(
1141 &lastdeg, Bbarn, Bbar, modulus, alphapow, ctx))
1142 {
1143 fmpz_mod_mpoly_cvtfrom_mpolyn(rBbar, Bbarn, m, ctx);
1144 if (m == nvars - 1)
1145 {
1146 success = fmpz_mod_mpolyl_content(cont, rBbar, 2, ctx);
1147 if (!success)
1148 goto cleanup;
1149 fmpz_mod_mpoly_divides(rBbar, rBbar, cont, ctx);
1150 fmpz_mod_mpoly_repack_bits_inplace(rBbar, bits, ctx);
1151 if (fmpz_mod_mpoly_divides(rG, B, rBbar, ctx) &&
1152 fmpz_mod_mpoly_divides(rAbar, A, rG, ctx))
1153 {
1154 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
1155 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1156 break;
1157 }
1158 }
1159 else
1160 {
1161 fmpz_mod_mpoly_mul(T, Bevals + m + 1, gammaevals + m + 1, ctx);
1162 if (fmpz_mod_mpoly_divides(rG, T, rBbar, ctx))
1163 {
1164 fmpz_mod_mpoly_mul(T, Aevals + m + 1, gammaevals + m + 1, ctx);
1165 if (fmpz_mod_mpoly_divides(rAbar, T, rG, ctx))
1166 {
1167 fmpz_mod_mpoly_repack_bits_inplace(rG, bits, ctx);
1168 fmpz_mod_mpoly_repack_bits_inplace(rAbar, bits, ctx);
1169 fmpz_mod_mpoly_swap(G, rG, ctx);
1170 fmpz_mod_mpoly_swap(Abar, rAbar, ctx);
1171 fmpz_mod_mpoly_swap(Bbar, rBbar, ctx);
1172 break;
1173 }
1174 }
1175 }
1176 }
1177
1178 if (_fmpz_mod_poly_degree(modulus) > gammadegs[m] + Adegs[m] &&
1179 _fmpz_mod_poly_degree(modulus) > gammadegs[m] + Bdegs[m])
1180 {
1181 goto choose_alphas;
1182 }
1183
1184 fmpz_mod_neg(c, alphas + m, ctx->ffinfo);
1185 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 1, c, ctx->ffinfo);
1186 }
1187 }
1188
1189 success = 1;
1190
1191 cleanup:
1192
1193 n_poly_clear(tmark);
1194
1195 fmpz_mod_polyun_clear(HG, ctx->ffinfo);
1196 fmpz_mod_polyun_clear(HAbar, ctx->ffinfo);
1197 fmpz_mod_polyun_clear(HBbar, ctx->ffinfo);
1198 fmpz_mod_polyun_clear(MG, ctx->ffinfo);
1199 fmpz_mod_polyun_clear(MAbar, ctx->ffinfo);
1200 fmpz_mod_polyun_clear(MBbar, ctx->ffinfo);
1201 fmpz_mod_polyun_clear(ZG, ctx->ffinfo);
1202 fmpz_mod_polyun_clear(ZAbar, ctx->ffinfo);
1203 fmpz_mod_polyun_clear(ZBbar, ctx->ffinfo);
1204 fmpz_mod_polyun_clear(Aev, ctx->ffinfo);
1205 fmpz_mod_polyun_clear(Bev, ctx->ffinfo);
1206 fmpz_mod_polyun_clear(Gev, ctx->ffinfo);
1207 fmpz_mod_polyun_clear(Abarev, ctx->ffinfo);
1208 fmpz_mod_polyun_clear(Bbarev, ctx->ffinfo);
1209 fmpz_mod_poly_clear(alphapow, ctx->ffinfo);
1210 fmpz_mod_mpoly_clear(cont, ctx);
1211 fmpz_mod_mpoly_clear(T, ctx);
1212 fmpz_mod_mpoly_clear(G, ctx);
1213 fmpz_mod_mpoly_clear(Abar, ctx);
1214 fmpz_mod_mpoly_clear(Bbar, ctx);
1215 fmpz_mod_mpolyn_clear(Tn, ctx);
1216 fmpz_mod_mpolyn_clear(Gn, ctx);
1217 fmpz_mod_mpolyn_clear(Abarn, ctx);
1218 fmpz_mod_mpolyn_clear(Bbarn, ctx);
1219 fmpz_mod_polyun_clear(Aeh_cur, ctx->ffinfo);
1220 fmpz_mod_polyun_clear(Aeh_inc, ctx->ffinfo);
1221 fmpz_mod_polyun_clear(Beh_cur, ctx->ffinfo);
1222 fmpz_mod_polyun_clear(Beh_inc, ctx->ffinfo);
1223 fmpz_mod_poly_clear(gammaeh_cur, ctx->ffinfo);
1224 fmpz_mod_poly_clear(gammaeh_inc, ctx->ffinfo);
1225 fmpz_mod_poly_clear(modulus, ctx->ffinfo);
1226 fmpz_mod_poly_stack_clear(St->poly_stack);
1227 fmpz_mod_polyun_stack_clear(St->polyun_stack);
1228
1229 flint_free(Aeh_coeff_mock->coeffs);
1230 flint_free(Beh_coeff_mock->coeffs);
1231
1232 for (i = 0; i < nvars; i++)
1233 {
1234 fmpz_mod_poly_clear(beta_caches + i, ctx->ffinfo);
1235 fmpz_mod_poly_clear(alpha_caches + i, ctx->ffinfo);
1236 }
1237 flint_free(beta_caches);
1238 flint_free(alpha_caches);
1239
1240 _fmpz_vec_clear(betas, nvars);
1241 _fmpz_vec_clear(alphas, nvars);
1242
1243 flint_randclear(state);
1244
1245 for (i = 0; i < nvars; i++)
1246 {
1247 fmpz_mod_mpoly_clear(Aevals + i, ctx);
1248 fmpz_mod_mpoly_clear(Bevals + i, ctx);
1249 fmpz_mod_mpoly_clear(gammaevals + i, ctx);
1250 }
1251 flint_free(Aevals);
1252 flint_free(Bevals);
1253 flint_free(gammaevals);
1254
1255 fmpz_clear(gammaev);
1256 fmpz_clear(c);
1257 fmpz_clear(start_alpha);
1258
1259 FLINT_ASSERT(!success || rG->bits == bits);
1260 FLINT_ASSERT(!success || rAbar->bits == bits);
1261 FLINT_ASSERT(!success || rBbar->bits == bits);
1262
1263 FLINT_ASSERT(success || fmpz_abs_fits_ui(fmpz_mod_mpoly_ctx_modulus(ctx)));
1264
1265 return success;
1266
1267 gcd_is_trivial:
1268
1269 fmpz_mod_mpoly_one(rG, ctx);
1270 fmpz_mod_mpoly_set(rAbar, A, ctx);
1271 fmpz_mod_mpoly_set(rBbar, B, ctx);
1272
1273 success = 1;
1274
1275 goto cleanup;
1276 }
1277
1278