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_mpoly_factor.h"
13
fmpz_mod_mpolyl_gcd_hensel_smprime(fmpz_mod_mpoly_t G,slong Gdeg,fmpz_mod_mpoly_t Abar,fmpz_mod_mpoly_t Bbar,const fmpz_mod_mpoly_t A,const fmpz_mod_mpoly_t B,const fmpz_mod_mpoly_ctx_t ctx)14 int fmpz_mod_mpolyl_gcd_hensel_smprime(
15 fmpz_mod_mpoly_t G, slong Gdeg, /* upperbound on deg_X(G) */
16 fmpz_mod_mpoly_t Abar,
17 fmpz_mod_mpoly_t Bbar,
18 const fmpz_mod_mpoly_t A,
19 const fmpz_mod_mpoly_t B,
20 const fmpz_mod_mpoly_ctx_t ctx)
21 {
22 int success, alphas_tries_remaining, gamma_is_one;
23 const slong n = ctx->minfo->nvars - 1;
24 slong i, k;
25 flint_bitcnt_t bits = A->bits;
26 fmpz * alphas, * prev_alphas;
27 fmpz_t q, mu1, mu2;
28 fmpz_mod_mpoly_struct * Aevals, * Bevals, * Hevals;
29 fmpz_mod_mpoly_struct * H; /* points to A, B, or Hevals + n */
30 fmpz_mod_mpoly_struct * Glcs, * Hlcs;
31 fmpz_mod_mpoly_struct Hfac[2], Htfac[2];
32 slong * Hdegs;
33 slong Adegx, Bdegx, gdegx;
34 fmpz_mod_mpoly_t t1, t2, g, abar, bbar, hbar;
35 flint_rand_t state;
36
37 FLINT_ASSERT(n > 0);
38 FLINT_ASSERT(A->length > 0);
39 FLINT_ASSERT(B->length > 0);
40 FLINT_ASSERT(bits <= FLINT_BITS);
41 FLINT_ASSERT(A->bits == bits);
42 FLINT_ASSERT(B->bits == bits);
43 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
44
45 flint_randinit(state);
46
47 Hdegs = FLINT_ARRAY_ALLOC(n + 1, slong);
48
49 Glcs = FLINT_ARRAY_ALLOC(3*(n + 1), fmpz_mod_mpoly_struct);
50 Hlcs = Glcs + (n + 1);
51 Hevals = Hlcs + (n + 1);
52 for (i = 0; i < n + 1; i++)
53 {
54 fmpz_mod_mpoly_init(Glcs + i, ctx);
55 fmpz_mod_mpoly_init(Hlcs + i, ctx);
56 fmpz_mod_mpoly_init(Hevals + i, ctx);
57 }
58
59 alphas = _fmpz_vec_init(2*n);
60 prev_alphas = alphas + n;
61 Aevals = FLINT_ARRAY_ALLOC(2*(n + 1), fmpz_mod_mpoly_struct);
62 Bevals = Aevals + (n + 1);
63 for (i = 0; i < n; i++)
64 {
65 fmpz_mod_mpoly_init(Aevals + i, ctx);
66 fmpz_mod_mpoly_init(Bevals + i, ctx);
67 }
68
69 fmpz_init(q);
70 fmpz_init(mu1);
71 fmpz_init(mu2);
72
73 fmpz_mod_mpoly_init(t1, ctx);
74 fmpz_mod_mpoly_init(t2, ctx);
75 fmpz_mod_mpoly_init(g, ctx);
76 fmpz_mod_mpoly_init(abar, ctx);
77 fmpz_mod_mpoly_init(bbar, ctx);
78 fmpz_mod_mpoly_init(hbar, ctx);
79
80 fmpz_mod_mpoly_init(Hfac + 0, ctx);
81 fmpz_mod_mpoly_init(Hfac + 1, ctx);
82 fmpz_mod_mpoly_init(Htfac + 0, ctx);
83 fmpz_mod_mpoly_init(Htfac + 1, ctx);
84
85 /* init done */
86
87 alphas_tries_remaining = 10;
88
89 /* try all zeros first */
90 for (i = 0; i < n; i++)
91 {
92 /* no previous at this point */
93 fmpz_set(prev_alphas + i, fmpz_mod_ctx_modulus(ctx->ffinfo));
94 fmpz_zero(alphas + i);
95 }
96
97 goto got_alpha;
98
99 next_alpha:
100
101 if (--alphas_tries_remaining < 0)
102 {
103 success = 0;
104 goto cleanup;
105 }
106
107 for (i = 0; i < n; i++)
108 {
109 do {
110 fmpz_mod_rand(alphas + i, state, ctx->ffinfo);
111 } while (fmpz_equal(alphas + i, prev_alphas + i));
112 }
113
114 got_alpha:
115
116 /* ensure deg_X do not drop under evaluation */
117 Adegx = fmpz_mod_mpoly_degree_si(A, 0, ctx);
118 Bdegx = fmpz_mod_mpoly_degree_si(B, 0, ctx);
119 for (i = n - 1; i >= 0; i--)
120 {
121 fmpz_mod_mpoly_evaluate_one_fmpz(Aevals + i, i == n - 1 ? A :
122 Aevals + i + 1, i + 1, alphas + i, ctx);
123 fmpz_mod_mpoly_evaluate_one_fmpz(Bevals + i, i == n - 1 ? B :
124 Bevals + i + 1, i + 1, alphas + i, ctx);
125 if (Adegx != fmpz_mod_mpoly_degree_si(Aevals + i, 0, ctx) ||
126 Bdegx != fmpz_mod_mpoly_degree_si(Bevals + i, 0, ctx))
127 {
128 goto next_alpha;
129 }
130 }
131
132 /* univariate gcd */
133 success = fmpz_mod_mpoly_gcd_cofactors(g, abar, bbar,
134 Aevals + 0, Bevals + 0, ctx) &&
135 fmpz_mod_mpoly_gcd(t1, g, abar, ctx) &&
136 fmpz_mod_mpoly_gcd(t2, g, bbar, ctx);
137 if (!success)
138 goto cleanup;
139
140 gdegx = fmpz_mod_mpoly_degree_si(g, 0, ctx);
141
142 if (gdegx == 0)
143 {
144 /* G is trivial */
145 fmpz_mod_mpoly_set(Abar, A, ctx);
146 fmpz_mod_mpoly_set(Bbar, B, ctx);
147 fmpz_mod_mpoly_one(G, ctx);
148 success = 1;
149 goto cleanup;
150 }
151 else if (gdegx > Gdeg)
152 {
153 goto next_alpha;
154 }
155 else if (gdegx < Gdeg)
156 {
157 Gdeg = gdegx;
158 for (i = 0; i < n; i++)
159 fmpz_set(prev_alphas + i, alphas + i);
160
161 goto next_alpha;
162 }
163
164 /* the degbound gdegx (== Gdeg) has at least two witnesses now */
165
166 if (gdegx == Adegx)
167 {
168 if (fmpz_mod_mpoly_divides(Bbar, B, A, ctx))
169 {
170 fmpz_mod_mpoly_set(G, A, ctx);
171 fmpz_mod_mpoly_one(Abar, ctx);
172 success = 1;
173 goto cleanup;
174 }
175
176 goto next_alpha;
177 }
178 else if (gdegx == Bdegx)
179 {
180 if (fmpz_mod_mpoly_divides(Abar, A, B, ctx))
181 {
182 fmpz_mod_mpoly_set(G, B, ctx);
183 fmpz_mod_mpoly_one(Bbar, ctx);
184 success = 1;
185 goto cleanup;
186 }
187
188 goto next_alpha;
189 }
190
191 FLINT_ASSERT(0 < gdegx && gdegx < FLINT_MIN(Adegx, Bdegx));
192
193 /* set Hlcs[n], Glcs[n] (gamma), H, and Hevals */
194 if (fmpz_mod_mpoly_is_one(t1, ctx))
195 {
196 fmpz_one(mu1);
197 fmpz_zero(mu2);
198
199 fmpz_mod_mpoly_swap(hbar, abar, ctx);
200
201 fmpz_mod_mpolyl_lead_coeff(Hlcs + n, A, 1, ctx);
202 fmpz_mod_mpolyl_lead_coeff(t2, B, 1, ctx);
203 success = fmpz_mod_mpoly_gcd(Glcs + n, Hlcs + n, t2, ctx);
204 if (!success)
205 goto cleanup;
206
207 H = (fmpz_mod_mpoly_struct *) A;
208
209 gamma_is_one = fmpz_mod_mpoly_is_one(Glcs + n, ctx);
210 if (gamma_is_one)
211 for (i = 0; i < n; i++)
212 fmpz_mod_mpoly_swap(Hevals + i, Aevals + i, ctx);
213 }
214 else if (fmpz_mod_mpoly_is_one(t2, ctx))
215 {
216 fmpz_zero(mu1);
217 fmpz_one(mu2);
218
219 fmpz_mod_mpoly_swap(hbar, bbar, ctx);
220
221 fmpz_mod_mpolyl_lead_coeff(Hlcs + n, B, 1, ctx);
222 fmpz_mod_mpolyl_lead_coeff(t2, A, 1, ctx);
223 success = fmpz_mod_mpoly_gcd(Glcs + n, Hlcs + n, t2, ctx);
224 if (!success)
225 goto cleanup;
226
227 H = (fmpz_mod_mpoly_struct *) B;
228
229 gamma_is_one = fmpz_mod_mpoly_is_one(Glcs + n, ctx);
230 if (gamma_is_one)
231 for (i = 0; i < n; i++)
232 fmpz_mod_mpoly_swap(Hevals + i, Bevals + i, ctx);
233 }
234 else
235 {
236 int mu_tries_remaining = 10;
237
238 next_mu:
239
240 if (--mu_tries_remaining < 0)
241 {
242 success = 0;
243 goto cleanup;
244 }
245
246 fmpz_one(mu1);
247 fmpz_mod_rand_not_zero(mu2, state, ctx->ffinfo);
248 fmpz_mod_mpoly_scalar_addmul_fmpz(hbar, abar, bbar, mu2, ctx);
249
250 /* make sure the linear combo did not drop degree */
251 if (fmpz_mod_mpoly_degree_si(hbar, 0, ctx) != FLINT_MAX(Adegx, Bdegx) - gdegx)
252 goto next_mu;
253
254 /* make sure the linear combo is prime to g */
255 success = fmpz_mod_mpoly_gcd(t1, hbar, g, ctx);
256 if (!success)
257 goto cleanup;
258
259 if (!fmpz_mod_mpoly_is_one(t1, ctx))
260 goto next_mu;
261
262 fmpz_mod_mpolyl_lead_coeff(t1, A, 1, ctx);
263 fmpz_mod_mpolyl_lead_coeff(t2, B, 1, ctx);
264 success = fmpz_mod_mpoly_gcd(Glcs + n, t1, t2, ctx);
265 if (!success)
266 goto cleanup;
267
268 H = Hevals + n;
269 fmpz_mod_mpoly_scalar_addmul_fmpz(H, A, B, mu2, ctx);
270 fmpz_mod_mpolyl_lead_coeff(Hlcs + n, H, 1, ctx);
271
272 gamma_is_one = fmpz_mod_mpoly_is_one(Glcs + n, ctx);
273 if (gamma_is_one)
274 for (i = 0; i < n; i++)
275 fmpz_mod_mpoly_scalar_addmul_fmpz(Hevals + i, Aevals + i,
276 Bevals + i, mu2, ctx);
277 }
278
279 if (!gamma_is_one)
280 {
281 fmpz_mod_mpoly_mul(Hevals + n, H, Glcs + n, ctx);
282 H = Hevals + n;
283 for (i = n - 1; i >= 0; i--)
284 fmpz_mod_mpoly_evaluate_one_fmpz(Hevals + i, Hevals + i + 1,
285 i + 1, alphas + i, ctx);
286 }
287
288 success = H->bits <= FLINT_BITS ||
289 fmpz_mod_mpoly_repack_bits_inplace(H, FLINT_BITS, ctx);
290 if (!success)
291 goto cleanup;
292
293 /* the evals should all fit in H->bits */
294 for (i = 0; i < n; i++)
295 fmpz_mod_mpoly_repack_bits_inplace(Hevals + i, H->bits, ctx);
296
297 fmpz_mod_mpoly_degrees_si(Hdegs, H, ctx);
298
299 /* computed evaluated leading coeffs */
300 for (i = n - 1; i >= 0; i--)
301 {
302 fmpz_mod_mpoly_evaluate_one_fmpz(Glcs + i, Glcs + i + 1, i + 1, alphas + i, ctx);
303 fmpz_mod_mpoly_evaluate_one_fmpz(Hlcs + i, Hlcs + i + 1, i + 1, alphas + i, ctx);
304 /* evaluation could have killed gamma */
305 if (fmpz_mod_mpoly_is_zero(Glcs + i, ctx) ||
306 fmpz_mod_mpoly_is_zero(Hlcs + i, ctx))
307 {
308 goto next_alpha;
309 }
310 }
311
312 /* make the leading coefficients match Glcs[0], Hlcs[0] */
313
314 FLINT_ASSERT(fmpz_mod_mpoly_is_fmpz(Glcs + 0, ctx) && Glcs[0].length == 1);
315 FLINT_ASSERT(fmpz_mod_mpoly_is_fmpz(Hlcs + 0, ctx) && Hlcs[0].length == 1);
316
317 fmpz_mod_inv(q, g->coeffs + 0, ctx->ffinfo);
318 fmpz_mod_mul(q, q, Glcs[0].coeffs + 0, ctx->ffinfo);
319 fmpz_mod_mpoly_scalar_mul_fmpz_mod_invertible(Hfac + 0, g, q, ctx);
320
321 fmpz_mod_inv(q, hbar->coeffs + 0, ctx->ffinfo);
322 fmpz_mod_mul(q, q, Hlcs[0].coeffs + 0, ctx->ffinfo);
323 fmpz_mod_mpoly_scalar_mul_fmpz_mod_invertible(Hfac + 1, hbar, q, ctx);
324
325 for (k = 1; k <= n; k++)
326 {
327 _fmpz_mod_mpoly_set_lead0(Htfac + 0, Hfac + 0, Glcs + k, ctx);
328 _fmpz_mod_mpoly_set_lead0(Htfac + 1, Hfac + 1, Hlcs + k, ctx);
329 success = fmpz_mod_mpoly_hlift(k, Htfac, 2, alphas,
330 k < n ? Hevals + k : H, Hdegs, ctx);
331 if (!success)
332 goto next_alpha;
333
334 fmpz_mod_mpoly_swap(Hfac + 0, Htfac + 0, ctx);
335 fmpz_mod_mpoly_swap(Hfac + 1, Htfac + 1, ctx);
336 }
337
338 success = fmpz_mod_mpolyl_content(t1, Hfac + 0, 1, ctx);
339 if (!success)
340 goto cleanup;
341
342 success = fmpz_mod_mpoly_divides(G, Hfac + 0, t1, ctx);
343 FLINT_ASSERT(success);
344
345 if (fmpz_is_zero(mu2))
346 {
347 FLINT_ASSERT(fmpz_is_one(mu1));
348 /* the division by t1 should succeed, but let's be careful */
349 fmpz_mod_mpolyl_lead_coeff(t1, G, 1, ctx);
350 success = fmpz_mod_mpoly_divides(Abar, Hfac + 1, t1, ctx) &&
351 fmpz_mod_mpoly_divides(Bbar, B, G, ctx);
352 }
353 else if (fmpz_is_zero(mu1))
354 {
355 FLINT_ASSERT(fmpz_is_one(mu2));
356 /* ditto */
357 fmpz_mod_mpolyl_lead_coeff(t1, G, 1, ctx);
358 success = fmpz_mod_mpoly_divides(Bbar, Hfac + 1, t1, ctx) &&
359 fmpz_mod_mpoly_divides(Abar, A, G, ctx);
360 }
361 else
362 {
363 FLINT_ASSERT(fmpz_is_one(mu1));
364 success = fmpz_mod_mpoly_divides(Abar, A, G, ctx) &&
365 fmpz_mod_mpoly_divides(Bbar, B, G, ctx);
366 }
367
368 if (!success)
369 goto next_alpha;
370
371 success = 1;
372
373 cleanup:
374
375 flint_randclear(state);
376
377 flint_free(Hdegs);
378
379 for (i = 0; i < n + 1; i++)
380 {
381 fmpz_mod_mpoly_clear(Glcs + i, ctx);
382 fmpz_mod_mpoly_clear(Hlcs + i, ctx);
383 fmpz_mod_mpoly_clear(Hevals + i, ctx);
384 }
385 flint_free(Glcs);
386
387 _fmpz_vec_clear(alphas, 2*n);
388
389 for (i = 0; i < n; i++)
390 {
391 fmpz_mod_mpoly_clear(Aevals + i, ctx);
392 fmpz_mod_mpoly_clear(Bevals + i, ctx);
393 }
394 flint_free(Aevals);
395
396 fmpz_clear(q);
397 fmpz_clear(mu1);
398 fmpz_clear(mu2);
399
400 fmpz_mod_mpoly_clear(t1, ctx);
401 fmpz_mod_mpoly_clear(t2, ctx);
402 fmpz_mod_mpoly_clear(g, ctx);
403 fmpz_mod_mpoly_clear(abar, ctx);
404 fmpz_mod_mpoly_clear(bbar, ctx);
405 fmpz_mod_mpoly_clear(hbar, ctx);
406
407 fmpz_mod_mpoly_clear(Hfac + 0, ctx);
408 fmpz_mod_mpoly_clear(Hfac + 1, ctx);
409 fmpz_mod_mpoly_clear(Htfac + 0, ctx);
410 fmpz_mod_mpoly_clear(Htfac + 1, ctx);
411
412 if (success)
413 {
414 fmpz_mod_mpoly_repack_bits_inplace(G, bits, ctx);
415 fmpz_mod_mpoly_repack_bits_inplace(Abar, bits, ctx);
416 fmpz_mod_mpoly_repack_bits_inplace(Bbar, bits, ctx);
417
418 FLINT_ASSERT(G->length > 0);
419 FLINT_ASSERT(Abar->length > 0);
420 FLINT_ASSERT(Bbar->length > 0);
421 }
422
423 return success;
424 }
425
426