1 /*
2 Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
10 */
11
12 #include "fmpz_mpoly.h"
13 #include "nmod_mpoly.h"
14
15 /* max = max abs coeffs(A) */
fmpz_mpoly_height(fmpz_t max,const fmpz_mpoly_t A,const fmpz_mpoly_ctx_t ctx)16 void fmpz_mpoly_height(
17 fmpz_t max,
18 const fmpz_mpoly_t A,
19 const fmpz_mpoly_ctx_t ctx)
20 {
21 slong i;
22 fmpz_t t;
23
24 fmpz_init(t);
25 fmpz_zero(max);
26
27 for (i = 0; i < A->length; i++)
28 {
29 fmpz_abs(t, A->coeffs + i);
30 if (fmpz_cmp(max, t) < 0)
31 fmpz_set(max, t);
32 }
33
34 fmpz_clear(t);
35 }
36
37 /* max = max abs coeffs(A), sum = sum abs coeffs(A) */
fmpz_mpoly_heights(fmpz_t max,fmpz_t sum,const fmpz_mpoly_t A,const fmpz_mpoly_ctx_t ctx)38 void fmpz_mpoly_heights(
39 fmpz_t max,
40 fmpz_t sum,
41 const fmpz_mpoly_t A,
42 const fmpz_mpoly_ctx_t ctx)
43 {
44 slong i;
45 fmpz_t t;
46
47 fmpz_init(t);
48 fmpz_zero(max);
49 fmpz_zero(sum);
50
51 for (i = 0; i < A->length; i++)
52 {
53 fmpz_abs(t, A->coeffs + i);
54 fmpz_add(sum, sum, t);
55 if (fmpz_cmp(max, t) < 0)
56 fmpz_set(max, t);
57 }
58
59 fmpz_clear(t);
60 }
61
62
63 /* The inputs A and B are modified */
fmpz_mpolyl_gcd_brown(fmpz_mpoly_t G,fmpz_mpoly_t Abar,fmpz_mpoly_t Bbar,fmpz_mpoly_t A,fmpz_mpoly_t B,const fmpz_mpoly_ctx_t ctx,const mpoly_gcd_info_t I)64 int fmpz_mpolyl_gcd_brown(
65 fmpz_mpoly_t G,
66 fmpz_mpoly_t Abar,
67 fmpz_mpoly_t Bbar,
68 fmpz_mpoly_t A,
69 fmpz_mpoly_t B,
70 const fmpz_mpoly_ctx_t ctx,
71 const mpoly_gcd_info_t I)
72 {
73 int success;
74 fmpz_t bound;
75 slong offset, shift;
76 mp_limb_t p, gammared;
77 fmpz_t gamma, modulus;
78 fmpz_t gnm, gns, anm, ans, bnm, bns;
79 fmpz_t cA, cB, cG, cAbar, cBbar;
80 fmpz_t temp;
81 fmpz_mpoly_t T;
82 nmod_mpolyn_t Gp, Abarp, Bbarp, Ap, Bp;
83 nmod_mpoly_ctx_t pctx;
84 flint_bitcnt_t bits = G->bits;
85 slong N = mpoly_words_per_exp_sp(bits, ctx->minfo);
86 nmod_poly_stack_t Sp;
87
88 FLINT_ASSERT(bits <= FLINT_BITS);
89 FLINT_ASSERT(bits == A->bits);
90 FLINT_ASSERT(bits == B->bits);
91 FLINT_ASSERT(bits == G->bits);
92 FLINT_ASSERT(bits == Abar->bits);
93 FLINT_ASSERT(bits == Bbar->bits);
94 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
95
96 mpoly_gen_offset_shift_sp(&offset, &shift, ctx->minfo->nvars - 1, bits, ctx->minfo);
97
98 fmpz_init(gamma);
99 fmpz_init(gnm);
100 fmpz_init(gns);
101 fmpz_init(anm);
102 fmpz_init(ans);
103 fmpz_init(bnm);
104 fmpz_init(bns);
105 fmpz_init(bound);
106 fmpz_init(temp);
107 fmpz_init_set_si(modulus, 1);
108
109 /* compute contents of G, Abar, Bbar, A, B */
110 fmpz_init(cA);
111 fmpz_init(cB);
112 fmpz_init(cG);
113 fmpz_init(cAbar);
114 fmpz_init(cBbar);
115 _fmpz_vec_content(cA, A->coeffs, A->length);
116 _fmpz_vec_content(cB, B->coeffs, B->length);
117 fmpz_gcd(cG, cA, cB);
118 fmpz_divexact(cAbar, cA, cG);
119 fmpz_divexact(cBbar, cB, cG);
120
121 /* remove content from inputs */
122 fmpz_mpoly_scalar_divexact_fmpz(A, A, cA, ctx);
123 fmpz_mpoly_scalar_divexact_fmpz(B, B, cB, ctx);
124
125 fmpz_gcd(gamma, A->coeffs + 0, B->coeffs + 0);
126
127 fmpz_mpoly_height(bound, A, ctx);
128 fmpz_mpoly_height(temp, B, ctx);
129 if (fmpz_cmp(bound, temp) < 0)
130 fmpz_swap(bound, temp);
131 fmpz_mul(bound, bound, gamma);
132 fmpz_add(bound, bound, bound);
133
134 fmpz_mpoly_init3(T, 0, bits, ctx);
135
136 nmod_mpoly_ctx_init(pctx, ctx->minfo->nvars, ORD_LEX, 2);
137 nmod_poly_stack_init(Sp, bits, pctx);
138 nmod_mpolyn_init(Ap, bits, pctx);
139 nmod_mpolyn_init(Bp, bits, pctx);
140 nmod_mpolyn_init(Gp, bits, pctx);
141 nmod_mpolyn_init(Abarp, bits, pctx);
142 nmod_mpolyn_init(Bbarp, bits, pctx);
143
144 p = UWORD(1) << (FLINT_BITS - 1);
145
146 choose_prime:
147
148 if (p >= UWORD_MAX_PRIME)
149 {
150 /* ran out of primes - absolute failure */
151 success = 0;
152 goto cleanup;
153 }
154 p = n_nextprime(p, 1);
155
156 /* make sure reduction does not kill both lc(A) and lc(B) */
157 gammared = fmpz_fdiv_ui(gamma, p);
158 if (gammared == 0)
159 {
160 goto choose_prime;
161 }
162
163 nmod_mpoly_ctx_set_modulus(pctx, p);
164 /* the unfortunate nmod poly's store their own context :( */
165 nmod_poly_stack_set_ctx(Sp, pctx);
166 nmod_mpolyn_set_mod(Ap, pctx->ffinfo->mod);
167 nmod_mpolyn_set_mod(Bp, pctx->ffinfo->mod);
168 nmod_mpolyn_set_mod(Gp, pctx->ffinfo->mod);
169 nmod_mpolyn_set_mod(Abarp, pctx->ffinfo->mod);
170 nmod_mpolyn_set_mod(Bbarp, pctx->ffinfo->mod);
171
172 /* reduction should kill neither A nor B */
173 fmpz_mpoly_interp_reduce_p_mpolyn(Ap, pctx, A, ctx);
174 fmpz_mpoly_interp_reduce_p_mpolyn(Bp, pctx, B, ctx);
175 FLINT_ASSERT(Ap->length > 0);
176 FLINT_ASSERT(Bp->length > 0);
177
178 success = nmod_mpolyn_gcd_brown_smprime(Gp, Abarp, Bbarp,
179 Ap, Bp, ctx->minfo->nvars - 1, pctx, I, Sp);
180 if (!success)
181 {
182 goto choose_prime;
183 }
184
185 if (nmod_mpolyn_is_nonzero_nmod(Gp, pctx))
186 {
187 fmpz_mpoly_one(G, ctx);
188 fmpz_mpoly_swap(Abar, A, ctx);
189 fmpz_mpoly_swap(Bbar, B, ctx);
190 goto successful_put_content;
191 }
192
193 if (!fmpz_is_one(modulus))
194 {
195 int cmp;
196 slong k;
197 FLINT_ASSERT(G->length > 0);
198
199 k = nmod_poly_degree(Gp->coeffs + 0);
200 cmp = mpoly_monomial_cmp_nomask_extra(G->exps + N*0,
201 Gp->exps + N*0, N, offset, k << shift);
202 if (cmp < 0)
203 {
204 goto choose_prime;
205 }
206 else if (cmp > 0)
207 {
208 fmpz_one(modulus);
209 }
210 }
211
212 FLINT_ASSERT(1 == nmod_mpolyn_leadcoeff(Gp, pctx));
213 nmod_mpolyn_scalar_mul_nmod(Gp, gammared, pctx);
214
215 if (!fmpz_is_one(modulus))
216 {
217 fmpz_mpoly_interp_crt_p_mpolyn(G, T, ctx, modulus, Gp, pctx);
218 fmpz_mpoly_interp_crt_p_mpolyn(Abar, T, ctx, modulus, Abarp, pctx);
219 fmpz_mpoly_interp_crt_p_mpolyn(Bbar, T, ctx, modulus, Bbarp, pctx);
220 }
221 else
222 {
223 fmpz_mpoly_interp_lift_p_mpolyn(G, ctx, Gp, pctx);
224 fmpz_mpoly_interp_lift_p_mpolyn(Abar, ctx, Abarp, pctx);
225 fmpz_mpoly_interp_lift_p_mpolyn(Bbar, ctx, Bbarp, pctx);
226 }
227
228 fmpz_mul_ui(modulus, modulus, p);
229
230 if (fmpz_cmp(modulus, bound) <= 0)
231 {
232 goto choose_prime;
233 }
234
235 fmpz_mpoly_heights(gnm, gns, G, ctx);
236 fmpz_mpoly_heights(anm, ans, Abar, ctx);
237 fmpz_mpoly_heights(bnm, bns, Bbar, ctx);
238 fmpz_mul(ans, ans, gnm);
239 fmpz_mul(anm, anm, gns);
240 fmpz_mul(bns, bns, gnm);
241 fmpz_mul(bnm, bnm, gns);
242
243 if (fmpz_cmp(ans, anm) > 0)
244 fmpz_swap(ans, anm);
245 if (fmpz_cmp(bns, bnm) > 0)
246 fmpz_swap(bns, bnm);
247 fmpz_add(ans, ans, ans);
248 fmpz_add(bns, bns, bns);
249 if (fmpz_cmp(ans, modulus) < 0 && fmpz_cmp(bns, modulus) < 0)
250 {
251 goto successful;
252 }
253
254 /* do not reset modulus to 1 */
255 goto choose_prime;
256
257 successful:
258
259 FLINT_ASSERT(fmpz_equal(gamma, G->coeffs + 0));
260
261 _fmpz_vec_content(temp, G->coeffs, G->length);
262 fmpz_mpoly_scalar_divexact_fmpz(G, G, temp, ctx);
263 fmpz_mpoly_scalar_divexact_fmpz(Abar, Abar, G->coeffs + 0, ctx);
264 fmpz_mpoly_scalar_divexact_fmpz(Bbar, Bbar, G->coeffs + 0, ctx);
265
266 successful_put_content:
267
268 fmpz_mpoly_scalar_mul_fmpz(G, G, cG, ctx);
269 fmpz_mpoly_scalar_mul_fmpz(Abar, Abar, cAbar, ctx);
270 fmpz_mpoly_scalar_mul_fmpz(Bbar, Bbar, cBbar, ctx);
271
272 success = 1;
273
274 cleanup:
275
276 fmpz_clear(cA);
277 fmpz_clear(cB);
278 fmpz_clear(cG);
279 fmpz_clear(cAbar);
280 fmpz_clear(cBbar);
281
282 fmpz_clear(gamma);
283 fmpz_clear(gnm);
284 fmpz_clear(gns);
285 fmpz_clear(anm);
286 fmpz_clear(ans);
287 fmpz_clear(bnm);
288 fmpz_clear(bns);
289 fmpz_clear(bound);
290 fmpz_clear(temp);
291 fmpz_clear(modulus);
292
293 nmod_mpolyn_clear(Gp, pctx);
294 nmod_mpolyn_clear(Abarp, pctx);
295 nmod_mpolyn_clear(Bbarp, pctx);
296 nmod_mpolyn_clear(Ap, pctx);
297 nmod_mpolyn_clear(Bp, pctx);
298 nmod_poly_stack_clear(Sp);
299 nmod_mpoly_ctx_clear(pctx);
300
301 fmpz_mpoly_clear(T, ctx);
302
303 return success;
304 }
305
306
fmpz_mpoly_gcd_brown(fmpz_mpoly_t G,const fmpz_mpoly_t A,const fmpz_mpoly_t B,const fmpz_mpoly_ctx_t ctx)307 int fmpz_mpoly_gcd_brown(
308 fmpz_mpoly_t G,
309 const fmpz_mpoly_t A,
310 const fmpz_mpoly_t B,
311 const fmpz_mpoly_ctx_t ctx)
312 {
313 int success;
314 slong * perm;
315 ulong * shift, * stride;
316 slong i;
317 flint_bitcnt_t wbits;
318 fmpz_mpoly_ctx_t lctx;
319 fmpz_mpoly_t Al, Bl, Gl, Abarl, Bbarl;
320
321 if (fmpz_mpoly_is_zero(A, ctx))
322 {
323 if (fmpz_mpoly_is_zero(B, ctx))
324 {
325 fmpz_mpoly_zero(G, ctx);
326 return 1;
327 }
328 if (fmpz_sgn(B->coeffs + 0) < 0)
329 {
330 fmpz_mpoly_neg(G, B, ctx);
331 return 1;
332 }
333 else
334 {
335 fmpz_mpoly_set(G, B, ctx);
336 return 1;
337 }
338 }
339
340 if (fmpz_mpoly_is_zero(B, ctx))
341 {
342 if (fmpz_sgn(A->coeffs + 0) < 0)
343 {
344 fmpz_mpoly_neg(G, A, ctx);
345 return 1;
346 }
347 else
348 {
349 fmpz_mpoly_set(G, A, ctx);
350 return 1;
351 }
352 }
353
354 if (A->bits > FLINT_BITS || B->bits > FLINT_BITS)
355 {
356 return 0;
357 }
358
359 perm = (slong *) flint_malloc(ctx->minfo->nvars*sizeof(slong));
360 shift = (ulong *) flint_malloc(ctx->minfo->nvars*sizeof(ulong));
361 stride = (ulong *) flint_malloc(ctx->minfo->nvars*sizeof(ulong));
362 for (i = 0; i < ctx->minfo->nvars; i++)
363 {
364 perm[i] = i;
365 shift[i] = 0;
366 stride[i] = 1;
367 }
368
369 if (ctx->minfo->nvars == 1)
370 {
371 fmpz_poly_t a, b, g;
372 fmpz_poly_init(a);
373 fmpz_poly_init(b);
374 fmpz_poly_init(g);
375 _fmpz_mpoly_to_fmpz_poly_deflate(a, A, 0, shift, stride, ctx);
376 _fmpz_mpoly_to_fmpz_poly_deflate(b, B, 0, shift, stride, ctx);
377 fmpz_poly_gcd(g, a, b);
378 _fmpz_mpoly_from_fmpz_poly_inflate(G, A->bits, g, 0, shift, stride, ctx);
379 fmpz_poly_clear(a);
380 fmpz_poly_clear(b);
381 fmpz_poly_clear(g);
382 success = 1;
383 goto cleanup1;
384 }
385
386 wbits = FLINT_MAX(A->bits, B->bits);
387
388 fmpz_mpoly_ctx_init(lctx, ctx->minfo->nvars, ORD_LEX);
389 fmpz_mpoly_init3(Al, 0, wbits, lctx);
390 fmpz_mpoly_init3(Bl, 0, wbits, lctx);
391 fmpz_mpoly_init3(Gl, 0, wbits, lctx);
392 fmpz_mpoly_init3(Abarl, 0, wbits, lctx);
393 fmpz_mpoly_init3(Bbarl, 0, wbits, lctx);
394
395 fmpz_mpoly_to_mpoly_perm_deflate_threaded_pool(Al, lctx, A, ctx,
396 perm, shift, stride, NULL, 0);
397 fmpz_mpoly_to_mpoly_perm_deflate_threaded_pool(Bl, lctx, B, ctx,
398 perm, shift, stride, NULL, 0);
399
400 success = fmpz_mpolyl_gcd_brown(Gl, Abarl, Bbarl, Al, Bl, lctx, NULL);
401 if (!success)
402 goto cleanup;
403
404 fmpz_mpoly_from_mpoly_perm_inflate(G, FLINT_MIN(A->bits, B->bits), ctx,
405 Gl, lctx, perm, shift, stride);
406 if (fmpz_sgn(G->coeffs + 0) < 0)
407 fmpz_mpoly_neg(G, G, ctx);
408
409 cleanup:
410
411 fmpz_mpoly_clear(Al, lctx);
412 fmpz_mpoly_clear(Bl, lctx);
413 fmpz_mpoly_clear(Gl, lctx);
414 fmpz_mpoly_clear(Abarl, lctx);
415 fmpz_mpoly_clear(Bbarl, lctx);
416 fmpz_mpoly_ctx_clear(lctx);
417
418 cleanup1:
419
420 flint_free(perm);
421 flint_free(shift);
422 flint_free(stride);
423
424 return success;
425 }
426
427