1 /*
2 Copyright (C) 2021 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpz_mod_mpoly_factor.h"
13
14
_fmpz_mod_poly_gcd_cofactors(fmpz_mod_poly_t G,fmpz_mod_poly_t Abar,fmpz_mod_poly_t Bbar,const fmpz_mod_poly_t A,const fmpz_mod_poly_t B,const fmpz_mod_ctx_t ctx,fmpz_mod_poly_t t)15 static void _fmpz_mod_poly_gcd_cofactors(
16 fmpz_mod_poly_t G,
17 fmpz_mod_poly_t Abar,
18 fmpz_mod_poly_t Bbar,
19 const fmpz_mod_poly_t A,
20 const fmpz_mod_poly_t B,
21 const fmpz_mod_ctx_t ctx,
22 fmpz_mod_poly_t t)
23 {
24 fmpz_mod_poly_gcd(G, A, B, ctx);
25 fmpz_mod_poly_divrem(Abar, t, A, G, ctx);
26 fmpz_mod_poly_divrem(Bbar, t, B, G, ctx);
27 }
28
29
fmpz_mod_polyu1n_gcd_brown_smprime(fmpz_mod_polyun_t G,fmpz_mod_polyun_t Abar,fmpz_mod_polyun_t Bbar,fmpz_mod_polyun_t A,fmpz_mod_polyun_t B,const fmpz_mod_ctx_t ctx,fmpz_mod_poly_stack_t St_poly,fmpz_mod_polyun_stack_t St_polyun)30 int fmpz_mod_polyu1n_gcd_brown_smprime(
31 fmpz_mod_polyun_t G,
32 fmpz_mod_polyun_t Abar,
33 fmpz_mod_polyun_t Bbar,
34 fmpz_mod_polyun_t A,
35 fmpz_mod_polyun_t B,
36 const fmpz_mod_ctx_t ctx,
37 fmpz_mod_poly_stack_t St_poly,
38 fmpz_mod_polyun_stack_t St_polyun)
39 {
40 int success;
41 slong bound;
42 fmpz_t alpha, temp, gammaevalp, gammaevalm;
43 fmpz_mod_poly_struct * Aevalp, * Bevalp, * Gevalp, * Abarevalp, * Bbarevalp;
44 fmpz_mod_poly_struct * Aevalm, * Bevalm, * Gevalm, * Abarevalm, * Bbarevalm;
45 fmpz_mod_polyun_struct * T;
46 slong deggamma, ldegG, ldegAbar, ldegBbar, ldegA, ldegB;
47 fmpz_mod_poly_struct * cA, * cB, * cG, * cAbar, * cBbar, * gamma;
48 fmpz_mod_poly_struct * modulus, * alphapow, * r;
49 int gstab, astab, bstab, use_stab;
50 #if FLINT_WANT_ASSERT
51 fmpz_mod_poly_t leadA, leadB;
52 #endif
53
54 fmpz_init(alpha);
55 fmpz_init(temp);
56 fmpz_init(gammaevalp);
57 fmpz_init(gammaevalm);
58
59 #if FLINT_WANT_ASSERT
60 fmpz_mod_poly_init(leadA, ctx);
61 fmpz_mod_poly_init(leadB, ctx);
62 fmpz_mod_poly_set(leadA, A->coeffs + 0, ctx);
63 fmpz_mod_poly_set(leadB, B->coeffs + 0, ctx);
64 #endif
65
66 fmpz_mod_poly_stack_fit_request(St_poly, 19);
67 cA = fmpz_mod_poly_stack_take_top(St_poly);
68 cB = fmpz_mod_poly_stack_take_top(St_poly);
69 cG = fmpz_mod_poly_stack_take_top(St_poly);
70 cAbar = fmpz_mod_poly_stack_take_top(St_poly);
71 cBbar = fmpz_mod_poly_stack_take_top(St_poly);
72 gamma = fmpz_mod_poly_stack_take_top(St_poly);
73 Aevalp = fmpz_mod_poly_stack_take_top(St_poly);
74 Bevalp = fmpz_mod_poly_stack_take_top(St_poly);
75 Gevalp = fmpz_mod_poly_stack_take_top(St_poly);
76 Abarevalp = fmpz_mod_poly_stack_take_top(St_poly);
77 Bbarevalp = fmpz_mod_poly_stack_take_top(St_poly);
78 Aevalm = fmpz_mod_poly_stack_take_top(St_poly);
79 Bevalm = fmpz_mod_poly_stack_take_top(St_poly);
80 Gevalm = fmpz_mod_poly_stack_take_top(St_poly);
81 Abarevalm = fmpz_mod_poly_stack_take_top(St_poly);
82 Bbarevalm = fmpz_mod_poly_stack_take_top(St_poly);
83 r = fmpz_mod_poly_stack_take_top(St_poly);
84 alphapow = fmpz_mod_poly_stack_take_top(St_poly);
85 modulus = fmpz_mod_poly_stack_take_top(St_poly);
86
87 fmpz_mod_polyun_stack_fit_request(St_polyun, 1);
88 T = fmpz_mod_polyun_stack_take_top(St_polyun);
89
90 _fmpz_mod_poly_vec_remove_content(cA, A->coeffs, A->length, ctx);
91 _fmpz_mod_poly_vec_remove_content(cB, B->coeffs, B->length, ctx);
92
93 _fmpz_mod_poly_gcd_cofactors(cG, cAbar, cBbar, cA, cB, ctx, r);
94
95 fmpz_mod_poly_gcd(gamma, A->coeffs + 0, B->coeffs + 0, ctx);
96
97 ldegA = _fmpz_mod_poly_vec_max_degree(A->coeffs, A->length, ctx);
98 ldegB = _fmpz_mod_poly_vec_max_degree(B->coeffs, B->length, ctx);
99 deggamma = fmpz_mod_poly_degree(gamma, ctx);
100 bound = 1 + deggamma + FLINT_MAX(ldegA, ldegB);
101
102 fmpz_mod_poly_fit_length(alphapow, FLINT_MAX(WORD(3), bound + 1), ctx);
103 fmpz_mod_poly_one(modulus, ctx);
104
105 use_stab = 1;
106 gstab = bstab = astab = 0;
107
108 fmpz_fdiv_q_2exp(alpha, fmpz_mod_ctx_modulus(ctx), 1);
109
110 choose_prime:
111
112 fmpz_sub_ui(alpha, alpha, 1);
113 if (fmpz_sgn(alpha) <= 0)
114 {
115 success = 0;
116 goto cleanup;
117 }
118
119 FLINT_ASSERT(alphapow->alloc >= 2);
120 alphapow->length = 2;
121 fmpz_one(alphapow->coeffs + 0);
122 fmpz_set(alphapow->coeffs + 1, alpha);
123
124 /* make sure evaluation point does not kill both lc(A) and lc(B) */
125 fmpz_mod_poly_eval2_pow(gammaevalp, gammaevalm, gamma, alphapow, ctx);
126 if (fmpz_is_zero(gammaevalp) || fmpz_is_zero(gammaevalm))
127 {
128 goto choose_prime;
129 }
130
131 /* evaluation point should kill neither A nor B */
132 fmpz_mod_polyu1n_interp_reduce_2sm_poly(Aevalp, Aevalm, A, alphapow, ctx);
133 fmpz_mod_polyu1n_interp_reduce_2sm_poly(Bevalp, Bevalm, B, alphapow, ctx);
134 FLINT_ASSERT(Aevalp->length > 0);
135 FLINT_ASSERT(Aevalm->length > 0);
136 FLINT_ASSERT(Bevalp->length > 0);
137 FLINT_ASSERT(Bevalm->length > 0);
138
139 if (use_stab && gstab)
140 {
141 slong Gdeg;
142 fmpz_mod_polyu1n_interp_reduce_2sm_poly(Gevalp, Gevalm, G, alphapow, ctx);
143 Gdeg = G->exps[0];
144 success = 1;
145 success = success && _fmpz_mod_poly_degree(Gevalp) == Gdeg;
146 success = success && _fmpz_mod_poly_degree(Gevalm) == Gdeg;
147 success = success && fmpz_equal(Gevalp->coeffs + Gdeg, gammaevalp);
148 success = success && fmpz_equal(Gevalm->coeffs + Gdeg, gammaevalm);
149 fmpz_mod_poly_divrem(Abarevalp, r, Aevalp, Gevalp, ctx);
150 success = success && (r->length == 0);
151 fmpz_mod_poly_divrem(Abarevalm, r, Aevalm, Gevalm, ctx);
152 success = success && (r->length == 0);
153 fmpz_mod_poly_divrem(Bbarevalp, r, Bevalp, Gevalp, ctx);
154 success = success && (r->length == 0);
155 fmpz_mod_poly_divrem(Bbarevalm, r, Bevalm, Gevalm, ctx);
156 success = success && (r->length == 0);
157
158 if (!success)
159 {
160 use_stab = 0;
161 fmpz_mod_poly_one(modulus, ctx);
162 goto choose_prime;
163 }
164
165 fmpz_mod_poly_scalar_mul_fmpz(Abarevalp, Abarevalp, gammaevalp, ctx);
166 fmpz_mod_poly_scalar_mul_fmpz(Abarevalm, Abarevalm, gammaevalm, ctx);
167 fmpz_mod_poly_scalar_mul_fmpz(Bbarevalp, Bbarevalp, gammaevalp, ctx);
168 fmpz_mod_poly_scalar_mul_fmpz(Bbarevalm, Bbarevalm, gammaevalm, ctx);
169 }
170 else
171 {
172 _fmpz_mod_poly_gcd_cofactors(Gevalp, Abarevalp, Bbarevalp,
173 Aevalp, Bevalp, ctx, r);
174 _fmpz_mod_poly_gcd_cofactors(Gevalm, Abarevalm, Bbarevalm,
175 Aevalm, Bevalm, ctx, r);
176 }
177
178 FLINT_ASSERT(Gevalp->length > 0);
179 FLINT_ASSERT(Abarevalp->length > 0);
180 FLINT_ASSERT(Bbarevalp->length > 0);
181 FLINT_ASSERT(Gevalm->length > 0);
182 FLINT_ASSERT(Abarevalm->length > 0);
183 FLINT_ASSERT(Bbarevalm->length > 0);
184
185 if (_fmpz_mod_poly_degree(Gevalp) == 0 || _fmpz_mod_poly_degree(Gevalm) == 0)
186 {
187 fmpz_mod_polyun_one(G, ctx);
188 fmpz_mod_polyun_swap(Abar, A);
189 fmpz_mod_polyun_swap(Bbar, B);
190 goto successful_put_content;
191 }
192
193 if (_fmpz_mod_poly_degree(Gevalp) != _fmpz_mod_poly_degree(Gevalm))
194 {
195 goto choose_prime;
196 }
197
198 if (_fmpz_mod_poly_degree(modulus) > 0)
199 {
200 FLINT_ASSERT(G->length > 0);
201 if (_fmpz_mod_poly_degree(Gevalp) > G->exps[0])
202 {
203 goto choose_prime;
204 }
205 else if (_fmpz_mod_poly_degree(Gevalp) < G->exps[0])
206 {
207 fmpz_mod_poly_one(modulus, ctx);
208 }
209 }
210
211 fmpz_mod_poly_scalar_mul_fmpz(Gevalp, Gevalp, gammaevalp, ctx);
212 fmpz_mod_poly_scalar_mul_fmpz(Gevalm, Gevalm, gammaevalm, ctx);
213 if (fmpz_mod_poly_degree(modulus, ctx) > 0)
214 {
215 fmpz_mod_poly_evaluate_fmpz(temp, modulus, alpha, ctx);
216 fmpz_mod_mul(temp, temp, alpha, ctx);
217 fmpz_mod_add(temp, temp, temp, ctx);
218 fmpz_mod_poly_scalar_div_fmpz(modulus, modulus, temp, ctx);
219 gstab = gstab || !fmpz_mod_polyu1n_interp_crt_2sm_poly(&ldegG, G, T, Gevalp, Gevalm, modulus, alphapow, ctx);
220 fmpz_mod_polyu1n_interp_crt_2sm_poly(&ldegAbar, Abar, T, Abarevalp, Abarevalm, modulus, alphapow, ctx);
221 fmpz_mod_polyu1n_interp_crt_2sm_poly(&ldegBbar, Bbar, T, Bbarevalp, Bbarevalm, modulus, alphapow, ctx);
222 }
223 else
224 {
225 fmpz_mod_polyu1n_interp_lift_2sm_poly(&ldegG, G, Gevalp, Gevalm, alpha, ctx);
226 fmpz_mod_polyu1n_interp_lift_2sm_poly(&ldegAbar, Abar, Abarevalp, Abarevalm, alpha, ctx);
227 fmpz_mod_polyu1n_interp_lift_2sm_poly(&ldegBbar, Bbar, Bbarevalp, Bbarevalm, alpha, ctx);
228 gstab = astab = bstab = 0;
229 }
230
231 fmpz_mod_mul(temp, alpha, alpha, ctx);
232 fmpz_mod_neg(temp, temp, ctx);
233 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 2, temp, ctx);
234
235 if (fmpz_mod_poly_degree(modulus, ctx) < bound)
236 {
237 goto choose_prime;
238 }
239
240 FLINT_ASSERT(ldegG >= 0);
241 FLINT_ASSERT(ldegAbar >= 0);
242 FLINT_ASSERT(ldegBbar >= 0);
243
244 if (deggamma + ldegA == ldegG + ldegAbar &&
245 deggamma + ldegB == ldegG + ldegBbar)
246 {
247 goto successful;
248 }
249
250 fmpz_mod_poly_one(modulus, ctx);
251 goto choose_prime;
252
253 successful:
254
255 _fmpz_mod_poly_vec_remove_content(r, G->coeffs, G->length, ctx);
256 _fmpz_mod_poly_vec_divexact_poly(Abar->coeffs, Abar->length, G->coeffs + 0, ctx);
257 _fmpz_mod_poly_vec_divexact_poly(Bbar->coeffs, Bbar->length, G->coeffs + 0, ctx);
258
259 successful_put_content:
260
261 _fmpz_mod_poly_vec_mul_poly(G->coeffs, G->length, cG, ctx);
262 _fmpz_mod_poly_vec_mul_poly(Abar->coeffs, Abar->length, cAbar, ctx);
263 _fmpz_mod_poly_vec_mul_poly(Bbar->coeffs, Bbar->length, cBbar, ctx);
264
265 success = 1;
266
267 cleanup:
268
269 #if FLINT_WANT_ASSERT
270 if (success)
271 {
272 FLINT_ASSERT(fmpz_is_one(fmpz_mod_poly_lead(G->coeffs + 0, ctx)));
273 fmpz_mod_poly_mul(modulus, G->coeffs + 0, Abar->coeffs + 0, ctx);
274 FLINT_ASSERT(fmpz_mod_poly_equal(modulus, leadA, ctx));
275 fmpz_mod_poly_mul(modulus, G->coeffs + 0, Bbar->coeffs + 0, ctx);
276 FLINT_ASSERT(fmpz_mod_poly_equal(modulus, leadB, ctx));
277 }
278 fmpz_mod_poly_clear(leadA, ctx);
279 fmpz_mod_poly_clear(leadB, ctx);
280 #endif
281
282 fmpz_mod_poly_stack_give_back(St_poly, 19);
283 fmpz_mod_polyun_stack_give_back(St_polyun, 1);
284
285 fmpz_clear(alpha);
286 fmpz_clear(temp);
287 fmpz_clear(gammaevalp);
288 fmpz_clear(gammaevalm);
289
290 return success;
291 }
292
293
fmpz_mod_mpolyn_set_polyun_swap(fmpz_mod_mpolyn_t A,fmpz_mod_polyun_t B,const fmpz_mod_mpoly_ctx_t ctx)294 static void fmpz_mod_mpolyn_set_polyun_swap(
295 fmpz_mod_mpolyn_t A,
296 fmpz_mod_polyun_t B,
297 const fmpz_mod_mpoly_ctx_t ctx)
298 {
299 slong i, N, off, shift;
300 flint_bitcnt_t bits = A->bits;
301 N = mpoly_words_per_exp_sp(bits, ctx->minfo);
302 mpoly_gen_offset_shift_sp(&off, &shift, 0, bits, ctx->minfo);
303
304 fmpz_mod_mpolyn_fit_length(A, B->length, ctx);
305 for (i = 0; i < B->length; i++)
306 {
307 mpoly_monomial_zero(A->exps + N*i, N);
308 (A->exps + N*i)[off] = B->exps[i] << shift;
309 fmpz_mod_poly_swap(A->coeffs + i, B->coeffs + i, ctx->ffinfo);
310 }
311
312 A->length = B->length;
313 }
314
fmpz_mod_mpolyn_get_polyun_swap(fmpz_mod_polyun_t B,fmpz_mod_mpolyn_t A,const fmpz_mod_mpoly_ctx_t ctx)315 static void fmpz_mod_mpolyn_get_polyun_swap(
316 fmpz_mod_polyun_t B,
317 fmpz_mod_mpolyn_t A,
318 const fmpz_mod_mpoly_ctx_t ctx)
319 {
320 slong i, N, off, shift;
321 flint_bitcnt_t bits = A->bits;
322 N = mpoly_words_per_exp_sp(bits, ctx->minfo);
323 mpoly_gen_offset_shift_sp(&off, &shift, 0, bits, ctx->minfo);
324
325 fmpz_mod_polyun_fit_length(B, A->length, ctx->ffinfo);
326 for (i = 0; i < A->length; i++)
327 {
328 B->exps[i] = (A->exps + N*i)[off] >> shift;
329 fmpz_mod_poly_swap(B->coeffs + i, A->coeffs + i, ctx->ffinfo);
330 }
331
332 B->length = A->length;
333 }
334
335
fmpz_mod_mpolyn_gcd_brown_smprime(fmpz_mod_mpolyn_t G,fmpz_mod_mpolyn_t Abar,fmpz_mod_mpolyn_t Bbar,fmpz_mod_mpolyn_t A,fmpz_mod_mpolyn_t B,slong var,const fmpz_mod_mpoly_ctx_t ctx,const mpoly_gcd_info_t I,fmpz_mod_poly_polyun_mpolyn_stack_t St)336 int fmpz_mod_mpolyn_gcd_brown_smprime(
337 fmpz_mod_mpolyn_t G,
338 fmpz_mod_mpolyn_t Abar,
339 fmpz_mod_mpolyn_t Bbar,
340 fmpz_mod_mpolyn_t A,
341 fmpz_mod_mpolyn_t B,
342 slong var,
343 const fmpz_mod_mpoly_ctx_t ctx,
344 const mpoly_gcd_info_t I,
345 fmpz_mod_poly_polyun_mpolyn_stack_t St)
346 {
347 int success, changed;
348 slong bound, upper_limit;
349 slong offset, shift;
350 fmpz_t alpha, temp, gammaevalp, gammaevalm;
351 fmpz_mod_mpolyn_struct * Aevalp, * Bevalp, * Gevalp, * Abarevalp, * Bbarevalp;
352 fmpz_mod_mpolyn_struct * Aevalm, * Bevalm, * Gevalm, * Abarevalm, * Bbarevalm;
353 fmpz_mod_mpolyn_struct * T1, * T2;
354 slong deggamma, ldegG, ldegAbar, ldegBbar, ldegA, ldegB;
355 fmpz_mod_poly_struct * cA, * cB, * cG, * cAbar, * cBbar, * gamma;
356 fmpz_mod_poly_struct * modulus, * alphapow, * t1;
357 flint_bitcnt_t bits = A->bits;
358 slong N = mpoly_words_per_exp_sp(bits, ctx->minfo);
359 #if FLINT_WANT_ASSERT
360 fmpz_mod_mpolyn_t Aorg, Borg;
361 fmpz_mod_poly_t leadA, leadB;
362 #endif
363
364 FLINT_ASSERT(bits == St->mpolyn_stack->bits);
365 FLINT_ASSERT(bits == A->bits);
366 FLINT_ASSERT(bits == B->bits);
367 FLINT_ASSERT(bits == G->bits);
368 FLINT_ASSERT(bits == Abar->bits);
369 FLINT_ASSERT(bits == Bbar->bits);
370
371 FLINT_ASSERT(var > 0);
372
373 if (var == 1)
374 {
375 fmpz_mod_polyun_struct * mA, * mB, * mG, * mAbar, * mBbar;
376
377 fmpz_mod_polyun_stack_fit_request(St->polyun_stack, 5);
378 mA = fmpz_mod_polyun_stack_take_top(St->polyun_stack);
379 mB = fmpz_mod_polyun_stack_take_top(St->polyun_stack);
380 mG = fmpz_mod_polyun_stack_take_top(St->polyun_stack);
381 mAbar = fmpz_mod_polyun_stack_take_top(St->polyun_stack);
382 mBbar = fmpz_mod_polyun_stack_take_top(St->polyun_stack);
383
384 fmpz_mod_mpolyn_get_polyun_swap(mA, A, ctx);
385 fmpz_mod_mpolyn_get_polyun_swap(mB, B, ctx);
386
387 success = fmpz_mod_polyu1n_gcd_brown_smprime(mG, mAbar, mBbar, mA, mB,
388 ctx->ffinfo, St->poly_stack, St->polyun_stack);
389 fmpz_mod_mpolyn_set_polyun_swap(G, mG, ctx);
390 fmpz_mod_mpolyn_set_polyun_swap(Abar, mAbar, ctx);
391 fmpz_mod_mpolyn_set_polyun_swap(Bbar, mBbar, ctx);
392
393 fmpz_mod_polyun_stack_give_back(St->polyun_stack, 5);
394
395 return success;
396 }
397
398 mpoly_gen_offset_shift_sp(&offset, &shift, var - 1, G->bits, ctx->minfo);
399
400 #if FLINT_WANT_ASSERT
401 fmpz_mod_poly_init(leadA, ctx->ffinfo);
402 fmpz_mod_poly_init(leadB, ctx->ffinfo);
403 fmpz_mod_poly_set(leadA, A->coeffs + 0, ctx->ffinfo);
404 fmpz_mod_poly_set(leadB, B->coeffs + 0, ctx->ffinfo);
405 fmpz_mod_mpolyn_init(Aorg, bits, ctx);
406 fmpz_mod_mpolyn_init(Borg, bits, ctx);
407 fmpz_mod_mpolyn_set(Aorg, A, ctx);
408 fmpz_mod_mpolyn_set(Borg, B, ctx);
409 #endif
410
411 fmpz_init(alpha);
412 fmpz_init(temp);
413 fmpz_init(gammaevalp);
414 fmpz_init(gammaevalm);
415
416 fmpz_mod_poly_stack_fit_request(St->poly_stack, 9);
417 cA = fmpz_mod_poly_stack_take_top(St->poly_stack);
418 cB = fmpz_mod_poly_stack_take_top(St->poly_stack);
419 cG = fmpz_mod_poly_stack_take_top(St->poly_stack);
420 cAbar = fmpz_mod_poly_stack_take_top(St->poly_stack);
421 cBbar = fmpz_mod_poly_stack_take_top(St->poly_stack);
422 gamma = fmpz_mod_poly_stack_take_top(St->poly_stack);
423 alphapow = fmpz_mod_poly_stack_take_top(St->poly_stack);
424 modulus = fmpz_mod_poly_stack_take_top(St->poly_stack);
425 t1 = fmpz_mod_poly_stack_take_top(St->poly_stack);
426
427 fmpz_mod_mpolyn_stack_fit_request(St->mpolyn_stack, 12, ctx);
428 Aevalp = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
429 Bevalp = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
430 Gevalp = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
431 Abarevalp = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
432 Bbarevalp = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
433 Aevalm = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
434 Bevalm = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
435 Gevalm = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
436 Abarevalm = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
437 Bbarevalm = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
438 T1 = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
439 T2 = fmpz_mod_mpolyn_stack_take_top(St->mpolyn_stack);
440
441 _fmpz_mod_poly_vec_remove_content(cA, A->coeffs, A->length, ctx->ffinfo);
442 _fmpz_mod_poly_vec_remove_content(cB, B->coeffs, B->length, ctx->ffinfo);
443 _fmpz_mod_poly_gcd_cofactors(cG, cAbar, cBbar, cA, cB, ctx->ffinfo, t1);
444 fmpz_mod_poly_gcd(gamma, A->coeffs + 0, B->coeffs + 0, ctx->ffinfo);
445
446 ldegA = fmpz_mod_mpolyn_lastdeg(A, ctx);
447 ldegB = fmpz_mod_mpolyn_lastdeg(B, ctx);
448 deggamma = fmpz_mod_poly_degree(gamma, ctx->ffinfo);
449 bound = 1 + deggamma + FLINT_MAX(ldegA, ldegB);
450
451 upper_limit = mpoly_gcd_info_get_brown_upper_limit(I, var, bound);
452
453 fmpz_mod_poly_fit_length(alphapow, FLINT_MAX(WORD(3), bound + 1), ctx->ffinfo);
454
455 fmpz_mod_poly_one(modulus, ctx->ffinfo);
456
457 fmpz_fdiv_q_2exp(alpha, fmpz_mod_ctx_modulus(ctx->ffinfo), 1);
458
459 choose_prime:
460
461 fmpz_sub_ui(alpha, alpha, 1);
462 if (fmpz_sgn(alpha) <= 0)
463 {
464 success = 0;
465 goto cleanup;
466 }
467
468 FLINT_ASSERT(alphapow->alloc >= 2);
469 alphapow->length = 2;
470 fmpz_one(alphapow->coeffs + 0);
471 fmpz_set(alphapow->coeffs + 1, alpha);
472
473 /* make sure evaluation point does not kill both lc(A) and lc(B) */
474 fmpz_mod_poly_eval2_pow(gammaevalp, gammaevalm, gamma, alphapow, ctx->ffinfo);
475 if (fmpz_is_zero(gammaevalp) || fmpz_is_zero(gammaevalm))
476 {
477 goto choose_prime;
478 }
479
480 /* evaluation point should kill neither A nor B */
481 fmpz_mod_mpolyn_interp_reduce_2sm_mpolyn(Aevalp, Aevalm, A, var, alphapow, ctx);
482 fmpz_mod_mpolyn_interp_reduce_2sm_mpolyn(Bevalp, Bevalm, B, var, alphapow, ctx);
483 FLINT_ASSERT(Aevalp->length > 0);
484 FLINT_ASSERT(Aevalm->length > 0);
485 FLINT_ASSERT(Bevalp->length > 0);
486 FLINT_ASSERT(Bevalm->length > 0);
487
488 success = fmpz_mod_mpolyn_gcd_brown_smprime(Gevalp, Abarevalp, Bbarevalp,
489 Aevalp, Bevalp, var - 1, ctx, I, St);
490 success = success &&
491 fmpz_mod_mpolyn_gcd_brown_smprime(Gevalm, Abarevalm, Bbarevalm,
492 Aevalm, Bevalm, var - 1, ctx, I, St);
493 if (success == 0)
494 {
495 goto choose_prime;
496 }
497
498 FLINT_ASSERT(Gevalp->length > 0);
499 FLINT_ASSERT(Abarevalp->length > 0);
500 FLINT_ASSERT(Bbarevalp->length > 0);
501 FLINT_ASSERT(Gevalm->length > 0);
502 FLINT_ASSERT(Abarevalm->length > 0);
503 FLINT_ASSERT(Bbarevalm->length > 0);
504
505 if (fmpz_mod_mpolyn_is_nonzero_fmpz(Gevalp, ctx) ||
506 fmpz_mod_mpolyn_is_nonzero_fmpz(Gevalm, ctx))
507 {
508 fmpz_mod_mpolyn_one(G, ctx);
509 fmpz_mod_mpolyn_swap(Abar, A, ctx);
510 fmpz_mod_mpolyn_swap(Bbar, B, ctx);
511 goto successful_put_content;
512 }
513
514 if (Gevalp->coeffs[0].length != Gevalm->coeffs[0].length ||
515 !mpoly_monomial_equal(Gevalp->exps + N*0, Gevalm->exps + N*0, N))
516 {
517 goto choose_prime;
518 }
519
520 /* the Geval have matching degrees */
521 if (fmpz_mod_poly_degree(modulus, ctx->ffinfo) > 0)
522 {
523 slong k = fmpz_mod_poly_degree(Gevalp->coeffs + 0, ctx->ffinfo);
524 int cmp = mpoly_monomial_cmp_nomask_extra(G->exps + N*0,
525 Gevalp->exps + N*0, N, offset, k << shift);
526 if (cmp < 0)
527 {
528 goto choose_prime;
529 }
530 else if (cmp > 0)
531 {
532 fmpz_mod_poly_one(modulus, ctx->ffinfo);
533 }
534 }
535
536 /* update interpolants */
537 fmpz_mod_mpolyn_scalar_mul_fmpz_mod(Gevalp, gammaevalp, ctx);
538 fmpz_mod_mpolyn_scalar_mul_fmpz_mod(Gevalm, gammaevalm, ctx);
539 if (fmpz_mod_poly_degree(modulus, ctx->ffinfo) > 0)
540 {
541 fmpz_mod_poly_evaluate_fmpz(temp, modulus, alpha, ctx->ffinfo);
542 fmpz_mod_mul(temp, temp, alpha, ctx->ffinfo);
543 fmpz_mod_add(temp, temp, temp, ctx->ffinfo);
544 fmpz_mod_poly_scalar_div_fmpz(modulus, modulus, temp, ctx->ffinfo);
545
546 changed = fmpz_mod_mpolyn_interp_crt_2sm_mpolyn(&ldegG, G, T1,
547 Gevalp, Gevalm, var, modulus, alphapow, ctx);
548 if (!changed && _fmpz_mod_poly_degree(modulus) < upper_limit)
549 {
550 _fmpz_mod_poly_vec_remove_content(t1, G->coeffs, G->length, ctx->ffinfo);
551 if (fmpz_mod_mpolyn_divides(T1, A, G, ctx) &&
552 fmpz_mod_mpolyn_divides(T2, B, G, ctx))
553 {
554 fmpz_mod_mpolyn_swap(Abar, T1, ctx);
555 fmpz_mod_mpolyn_swap(Bbar, T2, ctx);
556 goto successful_fix_lc;
557 }
558 _fmpz_mod_poly_vec_mul_poly(G->coeffs, G->length, t1, ctx->ffinfo);
559 }
560
561 changed = fmpz_mod_mpolyn_interp_crt_2sm_mpolyn(&ldegAbar, Abar, T1,
562 Abarevalp, Abarevalm, var, modulus, alphapow, ctx);
563 if (!changed && _fmpz_mod_poly_degree(modulus) < upper_limit)
564 {
565 _fmpz_mod_poly_vec_remove_content(t1, Abar->coeffs, Abar->length, ctx->ffinfo);
566 if (fmpz_mod_mpolyn_divides(T1, A, Abar, ctx) &&
567 fmpz_mod_mpolyn_divides(T2, B, T1, ctx))
568 {
569 fmpz_mod_mpolyn_swap(G, T1, ctx);
570 fmpz_mod_mpolyn_swap(Bbar, T2, ctx);
571 goto successful_fix_lc;
572 }
573 _fmpz_mod_poly_vec_mul_poly(Abar->coeffs, Abar->length, t1, ctx->ffinfo);
574 }
575
576 changed = fmpz_mod_mpolyn_interp_crt_2sm_mpolyn(&ldegBbar, Bbar, T1,
577 Bbarevalp, Bbarevalm, var, modulus, alphapow, ctx);
578 if (!changed && _fmpz_mod_poly_degree(modulus) < upper_limit)
579 {
580 _fmpz_mod_poly_vec_remove_content(t1, Bbar->coeffs, Bbar->length, ctx->ffinfo);
581 if (fmpz_mod_mpolyn_divides(T1, B, Bbar, ctx) &&
582 fmpz_mod_mpolyn_divides(T2, A, T1, ctx))
583 {
584 fmpz_mod_mpolyn_swap(G, T1, ctx);
585 fmpz_mod_mpolyn_swap(Abar, T2, ctx);
586 goto successful_fix_lc;
587 }
588 _fmpz_mod_poly_vec_mul_poly(Bbar->coeffs, Bbar->length, t1, ctx->ffinfo);
589 }
590 }
591 else
592 {
593 fmpz_mod_mpolyn_interp_lift_2sm_mpolyn(&ldegG, G,
594 Gevalp, Gevalm, var, alpha, ctx);
595 fmpz_mod_mpolyn_interp_lift_2sm_mpolyn(&ldegAbar, Abar,
596 Abarevalp, Abarevalm, var, alpha, ctx);
597 fmpz_mod_mpolyn_interp_lift_2sm_mpolyn(&ldegBbar, Bbar,
598 Bbarevalp, Bbarevalm, var, alpha, ctx);
599 }
600
601 fmpz_mod_mul(temp, alpha, alpha, ctx->ffinfo);
602 fmpz_mod_neg(temp, temp, ctx->ffinfo);
603 fmpz_mod_poly_shift_left_scalar_addmul_fmpz_mod(modulus, 2, temp, ctx->ffinfo);
604
605 if (fmpz_mod_poly_degree(modulus, ctx->ffinfo) < bound)
606 {
607 goto choose_prime;
608 }
609
610 FLINT_ASSERT(ldegG >= 0);
611 FLINT_ASSERT(ldegAbar >= 0);
612 FLINT_ASSERT(ldegBbar >= 0);
613
614 if (deggamma + ldegA == ldegG + ldegAbar &&
615 deggamma + ldegB == ldegG + ldegBbar)
616 {
617 goto successful;
618 }
619
620 fmpz_mod_poly_one(modulus, ctx->ffinfo);
621 goto choose_prime;
622
623 successful_fix_lc:
624
625 fmpz_set(temp, fmpz_mod_poly_lead(G->coeffs + 0, ctx->ffinfo));
626 if (fmpz_is_one(temp))
627 goto successful_put_content;
628
629 _fmpz_mod_poly_vec_mul_fmpz_mod(Abar->coeffs, Abar->length, temp, ctx->ffinfo);
630 _fmpz_mod_poly_vec_mul_fmpz_mod(Bbar->coeffs, Bbar->length, temp, ctx->ffinfo);
631 fmpz_mod_inv(temp, temp, ctx->ffinfo);
632 _fmpz_mod_poly_vec_mul_fmpz_mod(G->coeffs, G->length, temp, ctx->ffinfo);
633
634 goto successful_put_content;
635
636 successful:
637
638 _fmpz_mod_poly_vec_remove_content(t1, G->coeffs, G->length, ctx->ffinfo);
639 _fmpz_mod_poly_vec_divexact_poly(Abar->coeffs, Abar->length, G->coeffs + 0, ctx->ffinfo);
640 _fmpz_mod_poly_vec_divexact_poly(Bbar->coeffs, Bbar->length, G->coeffs + 0, ctx->ffinfo);
641
642 successful_put_content:
643
644 _fmpz_mod_poly_vec_mul_poly(G->coeffs, G->length, cG, ctx->ffinfo);
645 _fmpz_mod_poly_vec_mul_poly(Abar->coeffs, Abar->length, cAbar, ctx->ffinfo);
646 _fmpz_mod_poly_vec_mul_poly(Bbar->coeffs, Bbar->length, cBbar, ctx->ffinfo);
647
648 success = 1;
649
650 cleanup:
651
652 #if FLINT_WANT_ASSERT
653 if (success)
654 {
655 FLINT_ASSERT(fmpz_is_one(fmpz_mod_mpolyn_leadcoeff(G)));
656 fmpz_mod_poly_mul(modulus, G->coeffs + 0, Abar->coeffs + 0, ctx->ffinfo);
657 FLINT_ASSERT(fmpz_mod_poly_equal(modulus, leadA, ctx->ffinfo));
658 fmpz_mod_poly_mul(modulus, G->coeffs + 0, Bbar->coeffs + 0, ctx->ffinfo);
659 FLINT_ASSERT(fmpz_mod_poly_equal(modulus, leadB, ctx->ffinfo));
660 }
661 fmpz_mod_poly_clear(leadA, ctx->ffinfo);
662 fmpz_mod_poly_clear(leadB, ctx->ffinfo);
663 fmpz_mod_mpolyn_clear(Aorg, ctx);
664 fmpz_mod_mpolyn_clear(Borg, ctx);
665 #endif
666
667 fmpz_clear(alpha);
668 fmpz_clear(temp);
669 fmpz_clear(gammaevalp);
670 fmpz_clear(gammaevalm);
671
672 fmpz_mod_poly_stack_give_back(St->poly_stack, 9);
673 fmpz_mod_mpolyn_stack_give_back(St->mpolyn_stack, 12);
674
675 FLINT_ASSERT(success || fmpz_bits(fmpz_mod_mpoly_ctx_modulus(ctx)) < 60);
676
677 return success;
678 }
679
680