1 /*
2 Copyright (C) 2018 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 "nmod_mpoly.h"
13 #include "fq_nmod_mpoly.h"
14
fq_nmod_next(fq_nmod_t alpha,const fq_nmod_ctx_t fqctx)15 int fq_nmod_next(fq_nmod_t alpha, const fq_nmod_ctx_t fqctx)
16 {
17 slong i;
18 slong deg = nmod_poly_degree(fqctx->modulus);
19
20 for (i = 0; i < deg; i++)
21 {
22 ulong c = nmod_poly_get_coeff_ui(alpha, i);
23 c += UWORD(1);
24 if (c < fqctx->mod.n)
25 {
26 nmod_poly_set_coeff_ui(alpha, i, c);
27 return 1;
28 }
29 nmod_poly_set_coeff_ui(alpha, i, 0);
30 }
31
32 return 0;
33 }
34
fq_nmod_next_not_zero(fq_nmod_t alpha,const fq_nmod_ctx_t fqctx)35 void fq_nmod_next_not_zero(fq_nmod_t alpha, const fq_nmod_ctx_t fqctx)
36 {
37 slong i;
38 slong deg = fqctx->modulus->length - 1;
39
40 for (i = 0; i < deg; i++) {
41 ulong c = nmod_poly_get_coeff_ui(alpha, i);
42 c += UWORD(1);
43 if (c < fqctx->mod.n) {
44 nmod_poly_set_coeff_ui(alpha, i, c);
45 return;
46 }
47 nmod_poly_set_coeff_ui(alpha, i, UWORD(0));
48 }
49
50 /* we hit zero, so skip to 1 */
51 nmod_poly_set_coeff_ui(alpha, 0, UWORD(1));
52 }
53
54
55 /* store in each coefficient the evaluation of the corresponding monomial */
fq_nmod_mpoly_evalsk(fq_nmod_mpoly_t A,fq_nmod_mpoly_t B,slong entries,slong * offs,ulong * masks,fq_nmod_struct * powers,const fq_nmod_mpoly_ctx_t ctx)56 void fq_nmod_mpoly_evalsk(fq_nmod_mpoly_t A, fq_nmod_mpoly_t B,
57 slong entries, slong * offs, ulong * masks, fq_nmod_struct * powers,
58 const fq_nmod_mpoly_ctx_t ctx)
59 {
60 slong i, j;
61 slong N;
62
63 FLINT_ASSERT(A->bits == B->bits);
64 fq_nmod_mpoly_fit_length(A, B->length, ctx);
65 N = mpoly_words_per_exp(B->bits, ctx->minfo);
66 for (i = 0; i < B->length; i++)
67 {
68 fq_nmod_one(A->coeffs + i, ctx->fqctx);
69
70 for (j = 0; j < entries; j++)
71 {
72 if ((B->exps + N*i)[offs[j]] & masks[j])
73 {
74 fq_nmod_mul(A->coeffs + i, A->coeffs + i, powers + j, ctx->fqctx);
75 }
76 }
77
78 mpoly_monomial_set(A->exps + N*i, B->exps + N*i, N);
79 }
80 A->length = B->length;
81 }
82
fq_nmod_mpolyu_evalsk(fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,slong entries,slong * offs,ulong * masks,fq_nmod_struct * powers,const fq_nmod_mpoly_ctx_t ctx)83 void fq_nmod_mpolyu_evalsk(fq_nmod_mpolyu_t A, fq_nmod_mpolyu_t B,
84 slong entries, slong * offs, ulong * masks, fq_nmod_struct * powers,
85 const fq_nmod_mpoly_ctx_t ctx)
86 {
87 slong i;
88
89 fq_nmod_mpolyu_fit_length(A, B->length, ctx);
90 for (i = 0; i < B->length; i++)
91 {
92 A->exps[i] = B->exps[i];
93 fq_nmod_mpoly_evalsk(A->coeffs + i, B->coeffs + i,
94 entries, offs, masks, powers, ctx);
95 }
96 A->length = B->length;
97 }
98
99
100 /* multiply the coefficients of A pointwise by those of B */
fq_nmod_mpolyu_mulsk(fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,const fq_nmod_mpoly_ctx_t ctx)101 void fq_nmod_mpolyu_mulsk(fq_nmod_mpolyu_t A, fq_nmod_mpolyu_t B,
102 const fq_nmod_mpoly_ctx_t ctx)
103 {
104 slong i, j;
105
106 FLINT_ASSERT(A->length == B->length);
107 for (i = 0; i < A->length; i++)
108 {
109 FLINT_ASSERT(A->exps[i] == B->exps[i]);
110
111 FLINT_ASSERT((A->coeffs + i)->length == (B->coeffs + i)->length);
112 for (j = 0; j < (A->coeffs + i)->length; j++)
113 {
114 fq_nmod_mul((A->coeffs + i)->coeffs + j,
115 (A->coeffs + i)->coeffs + j,
116 (B->coeffs + i)->coeffs + j, ctx->fqctx);
117 }
118 }
119 }
120
121 /*
122 return 0 if the leading coeff of A vanishes
123 else return 1
124 */
fq_nmod_mpolyu_evalfromsk(fq_nmod_poly_t e,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t SK,const fq_nmod_mpoly_ctx_t ctx)125 int fq_nmod_mpolyu_evalfromsk(fq_nmod_poly_t e, fq_nmod_mpolyu_t A,
126 fq_nmod_mpolyu_t SK, const fq_nmod_mpoly_ctx_t ctx)
127 {
128 slong i, j;
129 int ret = 0;
130 fq_nmod_t acc, pp;
131
132 FLINT_ASSERT(A->length == SK->length);
133
134 fq_nmod_init(pp, ctx->fqctx);
135 fq_nmod_init(acc, ctx->fqctx);
136
137 fq_nmod_poly_zero(e, ctx->fqctx);
138 for (i = 0; i < A->length; i++)
139 {
140 FLINT_ASSERT((A->coeffs + i)->length == (SK->coeffs + i)->length);
141
142 fq_nmod_zero(acc, ctx->fqctx);
143
144 for (j = 0; j < (A->coeffs + i)->length; j++)
145 {
146 fq_nmod_mul(pp, (A->coeffs + i)->coeffs + j,
147 (SK->coeffs + i)->coeffs + j, ctx->fqctx);
148 fq_nmod_add(acc, acc, pp, ctx->fqctx);
149 }
150
151 fq_nmod_poly_set_coeff(e, A->exps[i], acc, ctx->fqctx);
152 ret |= (i == 0 && !fq_nmod_is_zero(acc, ctx->fqctx));
153 }
154
155 fq_nmod_clear(pp, ctx->fqctx);
156 fq_nmod_clear(acc, ctx->fqctx);
157
158 return ret;
159 }
160
161
fq_nmod_poly_product_roots(fq_nmod_poly_t P,fq_nmod_struct * r,slong n,const fq_nmod_ctx_t fqctx)162 void fq_nmod_poly_product_roots(fq_nmod_poly_t P, fq_nmod_struct * r,
163 slong n, const fq_nmod_ctx_t fqctx)
164 {
165 slong i;
166 fq_nmod_poly_t B;
167 fq_nmod_t a;
168 fq_nmod_init(a, fqctx);
169 fq_nmod_poly_init(B, fqctx);
170 fq_nmod_poly_one(P, fqctx);
171 fq_nmod_poly_gen(B, fqctx);
172 for (i = 0; i < n; i++)
173 {
174 fq_nmod_neg(a, r + i, fqctx);
175 fq_nmod_poly_set_coeff(B, 0, a, fqctx);
176 fq_nmod_poly_mul(P, P, B, fqctx);
177 }
178 fq_nmod_clear(a, fqctx);
179 fq_nmod_poly_clear(B, fqctx);
180 }
181
182 /*
183 solve
184
185 [ a[0] a[1] ... a[n-1] ] [ x[0] ] [ b[0] ]
186 [ a[0]^2 a[1]^2 ... a[n-1]^2 ] [ x[1] ] [ b[1] ]
187 [ ... ... ] . [ ... ] = [ ... ]
188 [ a[0]^n a[1]^n ... a[n-1]^n ] [ x[n-1] ] [ b[n-1] ]
189
190 for x
191 */
fq_nmod_vandsolve(fq_nmod_struct * x,fq_nmod_struct * a,fq_nmod_struct * b,slong n,const fq_nmod_ctx_t fqctx)192 int fq_nmod_vandsolve(fq_nmod_struct * x, fq_nmod_struct * a, fq_nmod_struct * b,
193 slong n, const fq_nmod_ctx_t fqctx)
194 {
195 int success = 0;
196 slong i, j;
197 fq_nmod_t t, s;
198 fq_nmod_t Dinv;
199 fq_nmod_poly_t Q, P, R, u;
200
201 fq_nmod_init(t, fqctx);
202 fq_nmod_init(s, fqctx);
203 fq_nmod_init(Dinv, fqctx);
204
205 for (i = 0; i < n; i++)
206 fq_nmod_zero(x + i, fqctx);
207
208 fq_nmod_poly_init(Q, fqctx);
209 fq_nmod_poly_init(P, fqctx);
210 fq_nmod_poly_init(R, fqctx);
211 fq_nmod_poly_init(u, fqctx);
212 fq_nmod_poly_gen(u, fqctx);
213 fq_nmod_poly_product_roots(P, a, n, fqctx);
214 for (i = 0; i < n; i++)
215 {
216 if (fq_nmod_is_zero(a + i, fqctx))
217 goto cleanup;
218
219 fq_nmod_neg(t, a + i, fqctx);
220 fq_nmod_poly_set_coeff(u, 0, t, fqctx);
221 fq_nmod_poly_divrem(Q, R, P, u, fqctx);
222 fq_nmod_poly_evaluate_fq_nmod(t, Q, a + i, fqctx);
223 fq_nmod_mul(t, a + i, t, fqctx);
224 if (fq_nmod_is_zero(t, fqctx))
225 goto cleanup;
226
227 fq_nmod_inv(Dinv, t, fqctx);
228 for (j = 0; j < n; j++)
229 {
230 fq_nmod_mul(t, b + j, Dinv, fqctx);
231 fq_nmod_poly_get_coeff(s, Q, j, fqctx);
232 fq_nmod_mul(t, t, s, fqctx);
233 fq_nmod_add(x + i, x + i, t, fqctx);
234 }
235 }
236 success = 1;
237
238 cleanup:
239 fq_nmod_poly_clear(Q, fqctx);
240 fq_nmod_poly_clear(P, fqctx);
241 fq_nmod_poly_clear(R, fqctx);
242 fq_nmod_poly_clear(u, fqctx);
243
244 fq_nmod_clear(t, fqctx);
245 fq_nmod_clear(s, fqctx);
246 fq_nmod_clear(Dinv, fqctx);
247
248 return success;
249 }
250
251 /*
252 Try to set G to the gcd of A and B given the form f of G.
253 return codes as enumerated in nmod_mpoly.h:
254
255 nmod_mpoly_gcds_success,
256 nmod_mpoly_gcds_form_wrong,
257 nmod_mpoly_gcds_no_solution,
258 nmod_mpoly_gcds_scales_not_found,
259 nmod_mpoly_gcds_eval_point_not_found,
260 nmod_mpoly_gcds_eval_gcd_deg_too_high
261 */
fq_nmod_mpolyu_gcds_zippel(fq_nmod_mpolyu_t G,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,fq_nmod_mpolyu_t f,slong var,const fq_nmod_mpoly_ctx_t ctx,flint_rand_t randstate,slong * degbound)262 nmod_gcds_ret_t fq_nmod_mpolyu_gcds_zippel(fq_nmod_mpolyu_t G,
263 fq_nmod_mpolyu_t A, fq_nmod_mpolyu_t B, fq_nmod_mpolyu_t f, slong var,
264 const fq_nmod_mpoly_ctx_t ctx, flint_rand_t randstate, slong * degbound)
265 {
266 int eval_points_tried;
267 nmod_gcds_ret_t success;
268 fq_nmod_mpolyu_t Aevalsk1, Bevalsk1, fevalsk1, Aevalski, Bevalski, fevalski;
269 fq_nmod_poly_t Aeval, Beval, Geval;
270 fq_nmod_struct * alpha, * b;
271 fq_nmod_mat_struct * M, * ML;
272 fq_nmod_mat_t MF, Msol;
273 int lc_ok;
274 int underdeterminedcount = 0;
275 int exceededcount = 0;
276 int * ML_is_initialized;
277 slong i, j, k, s, S, nullity;
278 slong * d;
279 slong l;
280 fq_nmod_struct * W;
281 slong entries;
282 slong * offs;
283 ulong * masks;
284 fq_nmod_struct * powers;
285 fq_nmod_t ck;
286 TMP_INIT;
287
288 FLINT_ASSERT(A->length > 0);
289 FLINT_ASSERT(B->length > 0);
290 FLINT_ASSERT(f->length > 0);
291
292 FLINT_ASSERT(A->bits == B->bits);
293 FLINT_ASSERT(A->bits == G->bits);
294 FLINT_ASSERT(A->bits == f->bits);
295 FLINT_ASSERT(var >= 0);
296
297 FLINT_ASSERT(*degbound == f->exps[0]);
298
299 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(A, ctx));
300 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(B, ctx));
301
302 FLINT_ASSERT(var > 0);
303
304 if (f->length == 1)
305 {
306 if ((f->coeffs + 0)->length > 1)
307 {
308 /* impossible to find scale factors in this case */
309 return nmod_gcds_scales_not_found;
310 }
311 else
312 {
313 /* otherwise set the coeff of the monomial to one */
314 nmod_gcds_ret_t ret;
315 FLINT_ASSERT((f->coeffs + 0)->length == 1);
316 fq_nmod_mpolyu_set(G, f, ctx);
317 fq_nmod_one((G->coeffs + 0)->coeffs + 0, ctx->fqctx);
318 fq_nmod_mpolyu_init(Aevalsk1, f->bits, ctx);
319 ret = nmod_gcds_form_wrong;
320 if ( fq_nmod_mpolyuu_divides(Aevalsk1, A, G, 1, ctx)
321 && fq_nmod_mpolyuu_divides(Aevalsk1, B, G, 1, ctx))
322 {
323 ret = nmod_gcds_success;
324 }
325 fq_nmod_mpolyu_clear(Aevalsk1, ctx);
326 return ret;
327 }
328 }
329
330 TMP_START;
331
332 fq_nmod_init(ck, ctx->fqctx);
333
334 fq_nmod_mpolyu_init(Aevalsk1, f->bits, ctx);
335 fq_nmod_mpolyu_init(Bevalsk1, f->bits, ctx);
336 fq_nmod_mpolyu_init(fevalsk1, f->bits, ctx);
337 fq_nmod_mpolyu_init(Aevalski, f->bits, ctx);
338 fq_nmod_mpolyu_init(Bevalski, f->bits, ctx);
339 fq_nmod_mpolyu_init(fevalski, f->bits, ctx);
340 fq_nmod_poly_init(Aeval, ctx->fqctx);
341 fq_nmod_poly_init(Beval, ctx->fqctx);
342 fq_nmod_poly_init(Geval, ctx->fqctx);
343
344 d = (slong *) TMP_ALLOC(f->length*sizeof(slong));
345 for (i = 0; i < f->length; i++)
346 {
347 d[i] = i;
348 }
349
350 /*
351 make d sort the coeffs so that
352 (f->coeffs + d[j-1])->length <= (f->coeffs + d[j-0])->length
353 for all j
354 */
355 for (i = 1; i<f->length; i++)
356 {
357 for (j=i; j > 0 && (f->coeffs + d[j-1])->length
358 > (f->coeffs + d[j-0])->length; j--)
359 {
360 slong temp = d[j-1];
361 d[j-1] = d[j-0];
362 d[j-0] = temp;
363 }
364 }
365
366 /* l is the number of images we will try to construct */
367 l = f->length - 3;
368 for (i = 0; i < f->length; i++)
369 {
370 l += (f->coeffs + i)->length;
371 }
372 l = l / (f->length - 1);
373 l = FLINT_MAX(l, (f->coeffs + d[f->length - 1])->length);
374 /* one extra test image */
375 l += 1;
376
377 alpha = (fq_nmod_struct *) TMP_ALLOC(var*sizeof(fq_nmod_struct));
378 ML = (fq_nmod_mat_struct *) TMP_ALLOC(f->length*sizeof(fq_nmod_mat_struct));
379 b = (fq_nmod_struct *) TMP_ALLOC((f->coeffs + d[f->length - 1])->length
380 *sizeof(fq_nmod_struct));
381
382 for (i = 0; i < var; i++)
383 fq_nmod_init(alpha + i, ctx->fqctx);
384 for (i = 0; i < (f->coeffs + d[f->length - 1])->length; i++)
385 fq_nmod_init(b + i, ctx->fqctx);
386
387 fq_nmod_mat_init(MF, 0, l, ctx->fqctx);
388
389 M = (fq_nmod_mat_struct *) TMP_ALLOC(f->length*sizeof(fq_nmod_mat_struct));
390 ML_is_initialized = (int *) TMP_ALLOC(f->length*sizeof(int));
391 for (i = 0; i < f->length; i++)
392 {
393 fq_nmod_mat_init(M + i, l, (f->coeffs + i)->length, ctx->fqctx);
394 ML_is_initialized[i] = 0;
395 }
396
397 W = (fq_nmod_struct *) flint_malloc(l*f->length*sizeof(fq_nmod_struct));
398 for (i = 0; i < l*f->length; i++)
399 fq_nmod_init(W + i, ctx->fqctx);
400
401
402 fq_nmod_mat_init(Msol, l, 1, ctx->fqctx);
403
404 /* compute how many masks are needed */
405 entries = f->bits * var;
406 offs = (slong *) TMP_ALLOC(entries*sizeof(slong));
407 masks = (ulong *) TMP_ALLOC(entries*sizeof(slong));
408 powers = (fq_nmod_struct *) TMP_ALLOC(entries*sizeof(fq_nmod_struct));
409 for (i = 0; i < entries; i++)
410 fq_nmod_init(powers + i, ctx->fqctx);
411
412
413 /***** evaluation loop head *******/
414 eval_points_tried = 0;
415 pick_evaluation_point:
416
417 if (++eval_points_tried > 10)
418 {
419 success = nmod_gcds_eval_point_not_found;
420 goto finished;
421 }
422
423 /* avoid 0 for the evaluation points */
424 for (i = 0; i < var; i++)
425 fq_nmod_randtest_not_zero(alpha + i, randstate, ctx->fqctx);
426
427 /* store bit masks for each power of two of the non-main variables */
428 for (i = 0; i < var; i++)
429 {
430 slong shift, off;
431 mpoly_gen_offset_shift_sp(&off, &shift, i, f->bits, ctx->minfo);
432 for (j = 0; j < f->bits; j++)
433 {
434 offs[f->bits*i + j] = off;
435 masks[f->bits*i + j] = UWORD(1) << (j + shift);
436 if (j == 0)
437 fq_nmod_set(powers + f->bits*i + j, alpha + i, ctx->fqctx);
438 else
439 fq_nmod_mul(powers + f->bits*i + j, powers + f->bits*i + j - 1,
440 powers + f->bits*i + j - 1, ctx->fqctx);
441 }
442 }
443
444 fq_nmod_mpolyu_evalsk(Aevalsk1, A, entries, offs, masks, powers, ctx);
445 fq_nmod_mpolyu_evalsk(Bevalsk1, B, entries, offs, masks, powers, ctx);
446 fq_nmod_mpolyu_evalsk(fevalsk1, f, entries, offs, masks, powers, ctx);
447
448 for (i = 0; i < l*f->length; i++)
449 {
450 fq_nmod_zero(W + i, ctx->fqctx);
451 }
452
453
454 for (i = 0; i < l; i++)
455 {
456 if (i == 0)
457 {
458 fq_nmod_mpolyu_set(Aevalski, Aevalsk1, ctx);
459 fq_nmod_mpolyu_set(Bevalski, Bevalsk1, ctx);
460 fq_nmod_mpolyu_set(fevalski, fevalsk1, ctx);
461 } else
462 {
463 fq_nmod_mpolyu_mulsk(Aevalski, Aevalsk1, ctx);
464 fq_nmod_mpolyu_mulsk(Bevalski, Bevalsk1, ctx);
465 fq_nmod_mpolyu_mulsk(fevalski, fevalsk1, ctx);
466 }
467
468 for (j = 0; j < f->length; j++)
469 {
470 for (k = 0; k < (f->coeffs + j)->length; k++)
471 {
472 fq_nmod_set((M + j)->rows[i] + k,
473 (fevalski->coeffs + j)->coeffs + k, ctx->fqctx);
474 }
475 }
476
477 lc_ok = 1;
478 lc_ok = lc_ok && fq_nmod_mpolyu_evalfromsk(Aeval, A, Aevalski, ctx);
479 lc_ok = lc_ok && fq_nmod_mpolyu_evalfromsk(Beval, B, Bevalski, ctx);
480 if (!lc_ok)
481 {
482 /* lc of A or B vanished */
483 goto pick_evaluation_point;
484 }
485
486 fq_nmod_poly_gcd(Geval, Aeval, Beval, ctx->fqctx);
487
488 if (f->exps[0] < fq_nmod_poly_degree(Geval, ctx->fqctx))
489 {
490 ++exceededcount;
491 if (exceededcount < 2)
492 goto pick_evaluation_point;
493
494 success = nmod_gcds_eval_gcd_deg_too_high;
495 goto finished;
496 }
497
498 if (f->exps[0] > fq_nmod_poly_degree(Geval, ctx->fqctx))
499 {
500 success = nmod_gcds_form_main_degree_too_high;
501 *degbound = fq_nmod_poly_degree(Geval, ctx->fqctx);
502 goto finished;
503 }
504
505 k = fq_nmod_poly_length(Geval, ctx->fqctx);
506 j = WORD(0);
507 while ((--k) >= 0)
508 {
509 fq_nmod_poly_get_coeff(ck, Geval, k, ctx->fqctx);
510 if (!fq_nmod_is_zero(ck, ctx->fqctx))
511 {
512 while (j < f->length && f->exps[j] > k)
513 {
514 j++;
515 }
516 if (j >= f->length || f->exps[j] != k)
517 {
518 success = nmod_gcds_form_wrong;
519 goto finished;
520 }
521 fq_nmod_set(W + l*j + i, ck, ctx->fqctx);
522 }
523 }
524 }
525
526 nullity = -1;
527 fq_nmod_mat_clear(MF, ctx->fqctx);
528 fq_nmod_mat_init(MF, 0, l, ctx->fqctx);
529
530 for (S = 0; S < f->length; S++)
531 {
532 s = d[S];
533
534 if (!ML_is_initialized[s])
535 {
536 fq_nmod_mat_init(ML + s, l, (f->coeffs + s)->length + l, ctx->fqctx);
537 ML_is_initialized[s] = 1;
538 for (i = 0; i < l; i++)
539 {
540 for (j = 0; j < (f->coeffs + s)->length; j++)
541 {
542 fq_nmod_set((ML + s)->rows[i] + j,
543 (M + s)->rows[i] + j, ctx->fqctx);
544 }
545 fq_nmod_set((ML + s)->rows[i] + (f->coeffs + s)->length + i,
546 W + l*s + i, ctx->fqctx);
547 }
548 } else {
549 for (i = 0; i < l; i++)
550 {
551 for (j = 0; j < (f->coeffs + s)->length; j++)
552 {
553 fq_nmod_set((ML + s)->rows[i] + j,
554 (M + s)->rows[i] + j, ctx->fqctx);
555 }
556 for (j = 0; j < l; j++) {
557 if (j == i)
558 {
559 fq_nmod_set((ML + s)->rows[i]
560 + (f->coeffs + s)->length + j,
561 W + l*s + i, ctx->fqctx);
562 } else {
563 fq_nmod_zero((ML + s)->rows[i]
564 + (f->coeffs + s)->length + j,
565 ctx->fqctx);
566 }
567 }
568 }
569
570 }
571 fq_nmod_mat_rref(ML + s, ctx->fqctx);
572
573 for (i = 0; i < (f->coeffs + s)->length; i++)
574 {
575 if (!fq_nmod_is_one((ML + s)->rows[i] + i, ctx->fqctx))
576 {
577 /* evaluation points produced a singular vandermonde matrix */
578 goto pick_evaluation_point;
579 }
580 }
581
582 {
583 /* appends rows to MF */
584 fq_nmod_mat_t MFtemp;
585 fq_nmod_mat_t Mwindow;
586
587 fq_nmod_mat_window_init(Mwindow, ML + s,
588 (f->coeffs + s)->length, (f->coeffs + s)->length,
589 l, (f->coeffs + s)->length + l, ctx->fqctx);
590 fq_nmod_mat_init(MFtemp,
591 fq_nmod_mat_nrows(MF, ctx->fqctx)
592 + l - (f->coeffs + s)->length,
593 l, ctx->fqctx);
594 fq_nmod_mat_concat_vertical(MFtemp, MF, Mwindow, ctx->fqctx);
595 fq_nmod_mat_swap(MFtemp, MF, ctx->fqctx);
596 fq_nmod_mat_clear(MFtemp, ctx->fqctx);
597 fq_nmod_mat_window_clear(Mwindow, ctx->fqctx);
598 }
599
600 nullity = l - fq_nmod_mat_rref(MF, ctx->fqctx);
601
602 if (nullity == 0)
603 {
604 /* There is no solution for scale factors. Form f must be wrong */
605 success = nmod_gcds_form_wrong;
606 goto finished;
607 }
608 if (nullity == 1)
609 {
610 /*
611 There is one solution for scale factors based on equations
612 considered thus far. Accept this as a solution and perform
613 checks of the remaining equations at the end.
614 */
615 break;
616 }
617 }
618
619 if (nullity != 1)
620 {
621 ++underdeterminedcount;
622 if (underdeterminedcount < 2)
623 goto pick_evaluation_point;
624
625 success = nmod_gcds_scales_not_found;
626 goto finished;
627 }
628
629 nullity = fq_nmod_mat_nullspace(Msol, MF, ctx->fqctx);
630 FLINT_ASSERT(nullity == 1);
631
632 fq_nmod_mpolyu_set(G, f, ctx);
633
634 for (i = 0; i < f->length; i++)
635 {
636 for (j = 0; j < (f->coeffs + i)->length; j++)
637 {
638 FLINT_ASSERT((f->coeffs + i)->length <= l);
639 fq_nmod_mul(b + j, W + l*i + j, Msol->rows[j] + 0, ctx->fqctx);
640 }
641 success = fq_nmod_vandsolve((G->coeffs + i)->coeffs,
642 (fevalsk1->coeffs + i)->coeffs, b,
643 (f->coeffs + i)->length, ctx->fqctx);
644 if (!success)
645 {
646 /* evaluation points produced a singular vandermonde matrix */
647 goto pick_evaluation_point;
648 }
649 }
650
651 /* check solution */
652 for (s = 0; s < f->length; s++)
653 {
654 fq_nmod_t acc, pp, u;
655 fq_nmod_init(acc, ctx->fqctx);
656 fq_nmod_init(pp, ctx->fqctx);
657 fq_nmod_init(u, ctx->fqctx);
658
659 for (i = 0; i < l; i++)
660 {
661 fq_nmod_zero(acc, ctx->fqctx);
662 for (j = 0; j < (f->coeffs + s)->length; j++)
663 {
664 fq_nmod_mul(pp, (M + s)->rows[i] + j,
665 (G->coeffs + s)->coeffs + j, ctx->fqctx);
666 fq_nmod_add(acc, acc, pp, ctx->fqctx);
667 }
668
669 fq_nmod_mul(u, W + l*s + i, Msol->rows[i] + 0, ctx->fqctx);
670 if (!fq_nmod_equal(acc, u, ctx->fqctx))
671 {
672 fq_nmod_clear(acc, ctx->fqctx);
673 fq_nmod_clear(pp, ctx->fqctx);
674 fq_nmod_clear(u, ctx->fqctx);
675 success = nmod_gcds_no_solution;
676 goto finished;
677 }
678 }
679
680 fq_nmod_clear(acc, ctx->fqctx);
681 fq_nmod_clear(pp, ctx->fqctx);
682 fq_nmod_clear(u, ctx->fqctx);
683 }
684
685 success = nmod_gcds_success;
686
687 finished:
688
689 for (i = 0; i < entries; i++)
690 fq_nmod_clear(powers + i, ctx->fqctx);
691
692 for (i = 0; i < var; i++)
693 fq_nmod_clear(alpha + i, ctx->fqctx);
694 for (i = 0; i < (f->coeffs + d[f->length - 1])->length; i++)
695 fq_nmod_clear(b + i, ctx->fqctx);
696 for (i = 0; i < l*f->length; i++)
697 fq_nmod_clear(W + i, ctx->fqctx);
698
699 flint_free(W);
700 fq_nmod_mat_clear(MF, ctx->fqctx);
701 fq_nmod_mat_clear(Msol, ctx->fqctx);
702 for (i = 0; i < f->length; i++)
703 {
704 fq_nmod_mat_clear(M + i, ctx->fqctx);
705 if (ML_is_initialized[i])
706 {
707 fq_nmod_mat_clear(ML + i, ctx->fqctx);
708 }
709 }
710
711 fq_nmod_clear(ck, ctx->fqctx);
712
713 fq_nmod_mpolyu_clear(Aevalsk1, ctx);
714 fq_nmod_mpolyu_clear(Bevalsk1, ctx);
715 fq_nmod_mpolyu_clear(fevalsk1, ctx);
716 fq_nmod_mpolyu_clear(Aevalski, ctx);
717 fq_nmod_mpolyu_clear(Bevalski, ctx);
718 fq_nmod_mpolyu_clear(fevalski, ctx);
719 fq_nmod_poly_clear(Aeval, ctx->fqctx);
720 fq_nmod_poly_clear(Beval, ctx->fqctx);
721 fq_nmod_poly_clear(Geval, ctx->fqctx);
722
723 TMP_END;
724 return success;
725 }
726
727
728 /* setform copies the exponents and zeros the coefficients */
fq_nmod_mpoly_setform_mpolyn(fq_nmod_mpoly_t A,fq_nmod_mpolyn_t B,const fq_nmod_mpoly_ctx_t ctx)729 void fq_nmod_mpoly_setform_mpolyn(fq_nmod_mpoly_t A, fq_nmod_mpolyn_t B,
730 const fq_nmod_mpoly_ctx_t ctx)
731 {
732 slong i;
733 slong N;
734
735 FLINT_ASSERT(A->bits == B->bits);
736
737 fq_nmod_mpoly_fit_length(A, B->length, ctx);
738 N = mpoly_words_per_exp(B->bits, ctx->minfo);
739 for (i = 0; i < B->length; i++)
740 {
741 fq_nmod_zero(A->coeffs + i, ctx->fqctx);
742 mpoly_monomial_set(A->exps + N*i, B->exps + N*i, N);
743 }
744 A->length = B->length;
745 }
746
fq_nmod_mpolyu_setform_mpolyun(fq_nmod_mpolyu_t A,fq_nmod_mpolyun_t B,const fq_nmod_mpoly_ctx_t ctx)747 void fq_nmod_mpolyu_setform_mpolyun(fq_nmod_mpolyu_t A, fq_nmod_mpolyun_t B,
748 const fq_nmod_mpoly_ctx_t ctx)
749 {
750 slong i;
751
752 fq_nmod_mpolyu_fit_length(A, B->length, ctx);
753 for (i = 0; i < B->length; i++)
754 {
755 FLINT_ASSERT((B->coeffs + i)->bits == B->bits);
756 fq_nmod_mpoly_setform_mpolyn(A->coeffs + i, B->coeffs + i, ctx);
757 A->exps[i] = B->exps[i];
758 }
759 A->length = B->length;
760 }
761
762
fq_nmod_mpolyu_gcdp_zippel_univar(fq_nmod_mpolyu_t G,fq_nmod_mpolyu_t Abar,fq_nmod_mpolyu_t Bbar,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,const fq_nmod_mpoly_ctx_t ctx)763 int fq_nmod_mpolyu_gcdp_zippel_univar(
764 fq_nmod_mpolyu_t G,
765 fq_nmod_mpolyu_t Abar,
766 fq_nmod_mpolyu_t Bbar,
767 fq_nmod_mpolyu_t A,
768 fq_nmod_mpolyu_t B,
769 const fq_nmod_mpoly_ctx_t ctx)
770 {
771 fq_nmod_poly_t a, b, g, t, r;
772 FLINT_ASSERT(A->bits == B->bits);
773 fq_nmod_poly_init(a, ctx->fqctx);
774 fq_nmod_poly_init(b, ctx->fqctx);
775 fq_nmod_poly_init(g, ctx->fqctx);
776 fq_nmod_poly_init(t, ctx->fqctx);
777 fq_nmod_poly_init(r, ctx->fqctx);
778 fq_nmod_mpolyu_cvtto_poly(a, A, ctx);
779 fq_nmod_mpolyu_cvtto_poly(b, B, ctx);
780 fq_nmod_poly_gcd(g, a, b, ctx->fqctx);
781 fq_nmod_mpolyu_cvtfrom_poly(G, g, ctx);
782 fq_nmod_poly_divrem(t, r, a, g, ctx->fqctx);
783 fq_nmod_mpolyu_cvtfrom_poly(Abar, t, ctx);
784 fq_nmod_poly_divrem(t, r, b, g, ctx->fqctx);
785 fq_nmod_mpolyu_cvtfrom_poly(Bbar, t, ctx);
786 fq_nmod_poly_clear(a, ctx->fqctx);
787 fq_nmod_poly_clear(b, ctx->fqctx);
788 fq_nmod_poly_clear(g, ctx->fqctx);
789 fq_nmod_poly_clear(t, ctx->fqctx);
790 fq_nmod_poly_clear(r, ctx->fqctx);
791 return 1;
792 }
793
fq_nmod_mpolyu_gcdp_zippel_univar_no_cofactors(fq_nmod_mpolyu_t G,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,const fq_nmod_mpoly_ctx_t ctx)794 int fq_nmod_mpolyu_gcdp_zippel_univar_no_cofactors(
795 fq_nmod_mpolyu_t G,
796 fq_nmod_mpolyu_t A,
797 fq_nmod_mpolyu_t B,
798 const fq_nmod_mpoly_ctx_t ctx)
799 {
800 fq_nmod_poly_t a, b, g;
801 FLINT_ASSERT(A->bits == B->bits);
802
803 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(A,ctx));
804 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(B,ctx));
805
806 fq_nmod_poly_init(a, ctx->fqctx);
807 fq_nmod_poly_init(b, ctx->fqctx);
808 fq_nmod_poly_init(g, ctx->fqctx);
809 fq_nmod_mpolyu_cvtto_poly(a, A, ctx);
810 fq_nmod_mpolyu_cvtto_poly(b, B, ctx);
811
812 fq_nmod_poly_gcd(g, a, b, ctx->fqctx);
813
814 fq_nmod_mpolyu_cvtfrom_poly(G, g, ctx);
815 fq_nmod_poly_clear(a, ctx->fqctx);
816 fq_nmod_poly_clear(b, ctx->fqctx);
817 fq_nmod_poly_clear(g, ctx->fqctx);
818 return 1;
819 }
820
821
fq_nmod_mpolyu_gcdp_zippel_bivar(fq_nmod_mpolyu_t G,fq_nmod_mpolyu_t Abar,fq_nmod_mpolyu_t Bbar,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,const fq_nmod_mpoly_ctx_t ctx,mpoly_zipinfo_t zinfo,flint_rand_t randstate)822 int fq_nmod_mpolyu_gcdp_zippel_bivar(
823 fq_nmod_mpolyu_t G,
824 fq_nmod_mpolyu_t Abar,
825 fq_nmod_mpolyu_t Bbar,
826 fq_nmod_mpolyu_t A,
827 fq_nmod_mpolyu_t B,
828 const fq_nmod_mpoly_ctx_t ctx,
829 mpoly_zipinfo_t zinfo,
830 flint_rand_t randstate)
831 {
832 slong var = 0;
833 slong Alastdeg;
834 slong Blastdeg;
835 ulong Ashift, Bshift, Gshift;
836 slong lastdeg;
837 slong bound;
838 int success = 0, changed, have_enough;
839 fq_nmod_poly_t a, b, c, g, modulus, tempmod;
840 fq_nmod_mpolyu_t Aeval, Beval, Geval;
841 fq_nmod_mpolyun_t An, Bn, H, Ht;
842 fq_nmod_t geval, temp, alpha, alphastart;
843 fmpz_t minusone;
844
845 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
846 FLINT_ASSERT(var >= -WORD(1));
847
848 FLINT_ASSERT(G->bits == A->bits);
849 FLINT_ASSERT(G->bits == B->bits);
850
851 fmpz_init(minusone);
852 fmpz_set_si(minusone, -WORD(1));
853
854 fq_nmod_mpolyun_init(An, A->bits, ctx);
855 fq_nmod_mpolyun_init(Bn, A->bits, ctx);
856 fq_nmod_mpolyu_cvtto_mpolyun(An, A, var, ctx);
857 fq_nmod_mpolyu_cvtto_mpolyun(Bn, B, var, ctx);
858
859 FLINT_ASSERT(An->length > 0);
860 FLINT_ASSERT(Bn->length > 0);
861
862 Ashift = A->exps[A->length - 1];
863 Bshift = B->exps[B->length - 1];
864 Gshift = FLINT_MIN(Ashift, Bshift);
865 fq_nmod_mpolyun_shift_right(An, Ashift);
866 fq_nmod_mpolyun_shift_right(Bn, Bshift);
867
868 fq_nmod_poly_init(a, ctx->fqctx);
869 fq_nmod_poly_init(b, ctx->fqctx);
870 fq_nmod_poly_init(c, ctx->fqctx);
871 fq_nmod_poly_init(g, ctx->fqctx);
872
873 /* if the gcd has content wrt last variable, we are going to fail */
874 fq_nmod_mpolyun_content_poly(a, An, ctx);
875 fq_nmod_mpolyun_content_poly(b, Bn, ctx);
876 fq_nmod_mpolyun_divexact_poly(An, An, a, ctx);
877 fq_nmod_mpolyun_divexact_poly(Bn, Bn, b, ctx);
878 fq_nmod_poly_gcd(c, a, b, ctx->fqctx);
879
880 fq_nmod_poly_gcd(g, fq_nmod_mpolyun_leadcoeff_poly(An, ctx),
881 fq_nmod_mpolyun_leadcoeff_poly(Bn, ctx), ctx->fqctx);
882 Alastdeg = fq_nmod_mpolyun_lastdeg(An, ctx);
883 Blastdeg = fq_nmod_mpolyun_lastdeg(Bn, ctx);
884
885 /* bound of the number of images required */
886 bound = 1 + FLINT_MIN(Alastdeg, Blastdeg)
887 + fq_nmod_poly_degree(g, ctx->fqctx);
888
889 fq_nmod_poly_init(modulus, ctx->fqctx);
890 fq_nmod_poly_init(tempmod, ctx->fqctx);
891 fq_nmod_poly_set_coeff_fmpz(tempmod, 1, minusone, ctx->fqctx);
892 fq_nmod_mpolyu_init(Aeval, A->bits, ctx);
893 fq_nmod_mpolyu_init(Beval, A->bits, ctx);
894 fq_nmod_mpolyu_init(Geval, A->bits, ctx);
895 fq_nmod_mpolyun_init(H, A->bits, ctx);
896 fq_nmod_mpolyun_init(Ht, A->bits, ctx);
897
898 fq_nmod_init(temp, ctx->fqctx);
899 fq_nmod_init(geval, ctx->fqctx);
900 fq_nmod_init(alpha, ctx->fqctx);
901 fq_nmod_init(alphastart, ctx->fqctx);
902
903 /* fail if the gcd has content wrt last variable */
904 if (fq_nmod_poly_degree(c, ctx->fqctx) > 0)
905 {
906 success = 0;
907 goto finished;
908 }
909
910 fq_nmod_poly_one(modulus, ctx->fqctx);
911 fq_nmod_mpolyun_zero(H, ctx);
912
913 fq_nmod_randtest_not_zero(alphastart, randstate, ctx->fqctx);
914 fq_nmod_set(alpha, alphastart, ctx->fqctx);
915
916 while (1)
917 {
918 /* get new evaluation point */
919 fq_nmod_next_not_zero(alpha, ctx->fqctx);
920
921 if (fq_nmod_equal(alpha, alphastart, ctx->fqctx))
922 {
923 success = 0;
924 goto finished;
925 }
926
927 /* make sure evaluation point does not kill both lc(A) and lc(B) */
928 fq_nmod_poly_evaluate_fq_nmod(geval, g, alpha, ctx->fqctx);
929
930 if (fq_nmod_is_zero(geval, ctx->fqctx))
931 goto outer_continue;
932
933 /* make sure evaluation point does not kill either A or B */
934 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Aeval, An, alpha, ctx);
935 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Beval, Bn, alpha, ctx);
936
937 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(Aeval, ctx));
938 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(Beval, ctx));
939
940 if (Aeval->length == 0 || Beval->length == 0)
941 goto outer_continue;
942
943 fq_nmod_mpolyu_gcdp_zippel_univar_no_cofactors(Geval, Aeval, Beval, ctx);
944
945 if (fq_nmod_mpolyu_is_one(Geval, ctx))
946 {
947 fq_nmod_mpolyu_one(G, ctx);
948 fq_nmod_mpolyu_swap(Abar, A);
949 fq_nmod_mpolyu_swap(Bbar, B);
950 fq_nmod_mpolyu_shift_left(G, Gshift);
951 fq_nmod_mpolyu_shift_left(Abar, Ashift - Gshift);
952 fq_nmod_mpolyu_shift_left(Bbar, Bshift - Gshift);
953 success = 1;
954 goto finished;
955 }
956
957 FLINT_ASSERT(Geval->length > 0);
958
959 if (fq_nmod_poly_degree(modulus, ctx->fqctx) > 0)
960 {
961 if (Geval->exps[0] > H->exps[0])
962 {
963 goto outer_continue;
964 }
965 else if (Geval->exps[0] < H->exps[0])
966 {
967 fq_nmod_poly_one(modulus, ctx->fqctx);
968 }
969 }
970
971 /* update interpolant H */
972 fq_nmod_inv(temp, fq_nmod_mpolyu_leadcoeff(Geval, ctx), ctx->fqctx);
973 fq_nmod_mul(temp, geval, temp, ctx->fqctx);
974 fq_nmod_mpolyu_scalar_mul_fq_nmod(Geval, temp, ctx);
975
976 if (fq_nmod_poly_degree(modulus, ctx->fqctx) > 0)
977 {
978 fq_nmod_poly_evaluate_fq_nmod(temp, modulus, alpha, ctx->fqctx);
979 fq_nmod_inv(temp, temp, ctx->fqctx);
980 fq_nmod_poly_scalar_mul_fq_nmod(modulus, modulus, temp, ctx->fqctx);
981 changed = fq_nmod_mpolyun_interp_crt_sm_mpolyu(&lastdeg, H, Ht, Geval,
982 modulus, alpha, ctx);
983 fq_nmod_poly_set_coeff(tempmod, 0, alpha, ctx->fqctx);
984 fq_nmod_poly_mul(modulus, modulus, tempmod, ctx->fqctx);
985
986 have_enough = fq_nmod_poly_degree(modulus, ctx->fqctx) >= bound;
987
988 if (changed && !have_enough)
989 {
990 goto outer_continue;
991 }
992
993 if (!changed || have_enough)
994 {
995 fq_nmod_mpolyun_content_poly(a, H, ctx);
996 fq_nmod_mpolyun_divexact_poly(Ht, H, a, ctx);
997 fq_nmod_mpolyun_shift_left(Ht, Gshift);
998 fq_nmod_mpolyu_cvtfrom_mpolyun(G, Ht, var, ctx);
999
1000 if ( fq_nmod_mpolyuu_divides(Abar, A, G, 1, ctx)
1001 && fq_nmod_mpolyuu_divides(Bbar, B, G, 1, ctx))
1002 {
1003 success = 1;
1004 goto finished;
1005 }
1006 }
1007
1008 if (have_enough)
1009 {
1010 fq_nmod_poly_one(modulus, ctx->fqctx);
1011 goto outer_continue;
1012 }
1013 }
1014 else
1015 {
1016 fq_nmod_mpolyun_interp_lift_sm_mpolyu(H, Geval, ctx);
1017 lastdeg = WORD(0);
1018 fq_nmod_poly_set_coeff(tempmod, 0, alpha, ctx->fqctx);
1019 fq_nmod_poly_mul(modulus, modulus, tempmod, ctx->fqctx);
1020 }
1021
1022 outer_continue:
1023 NULL;
1024 }
1025
1026 success = 0;
1027
1028 finished:
1029
1030 fmpz_clear(minusone);
1031 fq_nmod_clear(temp, ctx->fqctx);
1032 fq_nmod_clear(geval, ctx->fqctx);
1033 fq_nmod_clear(alpha, ctx->fqctx);
1034 fq_nmod_clear(alphastart, ctx->fqctx);
1035 fq_nmod_poly_clear(a, ctx->fqctx);
1036 fq_nmod_poly_clear(b, ctx->fqctx);
1037 fq_nmod_poly_clear(c, ctx->fqctx);
1038 fq_nmod_poly_clear(g, ctx->fqctx);
1039 fq_nmod_poly_clear(modulus, ctx->fqctx);
1040 fq_nmod_poly_clear(tempmod, ctx->fqctx);
1041 fq_nmod_mpolyu_clear(Aeval, ctx);
1042 fq_nmod_mpolyu_clear(Beval, ctx);
1043 fq_nmod_mpolyu_clear(Geval, ctx);
1044 fq_nmod_mpolyun_clear(An, ctx);
1045 fq_nmod_mpolyun_clear(Bn, ctx);
1046 fq_nmod_mpolyun_clear(H, ctx);
1047 fq_nmod_mpolyun_clear(Ht, ctx);
1048
1049 return success;
1050 }
1051
1052
1053
fq_nmod_mpolyu_gcdp_zippel(fq_nmod_mpolyu_t G,fq_nmod_mpolyu_t Abar,fq_nmod_mpolyu_t Bbar,fq_nmod_mpolyu_t A,fq_nmod_mpolyu_t B,slong var,const fq_nmod_mpoly_ctx_t ctx,mpoly_zipinfo_t zinfo,flint_rand_t randstate)1054 int fq_nmod_mpolyu_gcdp_zippel(
1055 fq_nmod_mpolyu_t G,
1056 fq_nmod_mpolyu_t Abar,
1057 fq_nmod_mpolyu_t Bbar,
1058 fq_nmod_mpolyu_t A,
1059 fq_nmod_mpolyu_t B,
1060 slong var,
1061 const fq_nmod_mpoly_ctx_t ctx,
1062 mpoly_zipinfo_t zinfo,
1063 flint_rand_t randstate)
1064 {
1065 slong lastdeg;
1066 slong Alastdeg;
1067 slong Blastdeg;
1068 ulong Ashift, Bshift, Gshift;
1069 slong degbound;
1070 slong bound;
1071 int success = 0, changed, have_enough;
1072 fq_nmod_mpolyun_t An, Bn;
1073 fq_nmod_poly_t a, b, c, g;
1074 fq_nmod_poly_t modulus, tempmod;
1075 fq_nmod_mpolyu_t Aeval, Beval, Geval, Abareval, Bbareval, Gform;
1076 fq_nmod_mpolyun_t H, Ht;
1077 fq_nmod_t geval, temp;
1078 fq_nmod_t alpha, alphastart;
1079 fmpz_t minusone;
1080
1081 FLINT_ASSERT(ctx->minfo->ord == ORD_LEX);
1082 FLINT_ASSERT(var >= -WORD(1));
1083
1084 FLINT_ASSERT(G->bits == A->bits);
1085 FLINT_ASSERT(G->bits == B->bits);
1086 FLINT_ASSERT(A->length > 0);
1087 FLINT_ASSERT(B->length > 0);
1088
1089 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(A,ctx));
1090 FLINT_ASSERT(fq_nmod_mpolyu_is_canonical(B,ctx));
1091
1092 if (var == -WORD(1))
1093 {
1094 /* no more variables left to interpolate */
1095 return fq_nmod_mpolyu_gcdp_zippel_univar(G, Abar, Bbar, A, B, ctx);
1096 }
1097
1098 if (var == WORD(0))
1099 {
1100 /* bivariate is more comfortable separated */
1101 return fq_nmod_mpolyu_gcdp_zippel_bivar(G, Abar, Bbar, A, B, ctx, zinfo, randstate);
1102 }
1103
1104 fq_nmod_mpolyun_init(An, A->bits, ctx);
1105 fq_nmod_mpolyun_init(Bn, A->bits, ctx);
1106 fq_nmod_mpolyu_cvtto_mpolyun(An, A, var, ctx);
1107 fq_nmod_mpolyu_cvtto_mpolyun(Bn, B, var, ctx);
1108
1109 Ashift = A->exps[A->length - 1];
1110 Bshift = B->exps[B->length - 1];
1111 Gshift = FLINT_MIN(Ashift, Bshift);
1112 fq_nmod_mpolyun_shift_right(An, Ashift);
1113 fq_nmod_mpolyun_shift_right(Bn, Bshift);
1114
1115 fq_nmod_poly_init(a, ctx->fqctx);
1116 fq_nmod_poly_init(b, ctx->fqctx);
1117 fq_nmod_poly_init(c, ctx->fqctx);
1118 fq_nmod_poly_init(g, ctx->fqctx);
1119
1120 /* if the gcd has content wrt last variable, we are going to fail */
1121 fq_nmod_mpolyun_content_poly(a, An, ctx);
1122 fq_nmod_mpolyun_content_poly(b, Bn, ctx);
1123 fq_nmod_mpolyun_divexact_poly(An, An, a, ctx);
1124 fq_nmod_mpolyun_divexact_poly(Bn, Bn, b, ctx);
1125 fq_nmod_poly_gcd(c, a, b, ctx->fqctx);
1126 fq_nmod_poly_gcd(g, fq_nmod_mpolyun_leadcoeff_poly(An, ctx),
1127 fq_nmod_mpolyun_leadcoeff_poly(Bn, ctx), ctx->fqctx);
1128 Alastdeg = fq_nmod_mpolyun_lastdeg(An, ctx);
1129 Blastdeg = fq_nmod_mpolyun_lastdeg(Bn, ctx);
1130
1131 /* bound of the number of images required */
1132 bound = 1 + FLINT_MIN(Alastdeg, Blastdeg)
1133 + fq_nmod_poly_degree(g, ctx->fqctx);
1134
1135 /* degree bound on the gcd */
1136 degbound = FLINT_MIN(A->exps[0], B->exps[0]);
1137
1138 fq_nmod_poly_init(modulus, ctx->fqctx);
1139 fq_nmod_poly_init(tempmod, ctx->fqctx);
1140 fmpz_init(minusone);
1141 fmpz_set_si(minusone, WORD(-1));
1142 fq_nmod_poly_set_coeff_fmpz(tempmod, 1, minusone, ctx->fqctx);
1143 fq_nmod_mpolyu_init(Aeval, A->bits, ctx);
1144 fq_nmod_mpolyu_init(Beval, A->bits, ctx);
1145 fq_nmod_mpolyu_init(Geval, A->bits, ctx);
1146 fq_nmod_mpolyu_init(Abareval, A->bits, ctx);
1147 fq_nmod_mpolyu_init(Bbareval, A->bits, ctx);
1148 fq_nmod_mpolyu_init(Gform, A->bits, ctx);
1149 fq_nmod_mpolyun_init(H, A->bits, ctx);
1150 fq_nmod_mpolyun_init(Ht, A->bits, ctx);
1151
1152 fq_nmod_init(geval, ctx->fqctx);
1153 fq_nmod_init(temp, ctx->fqctx);
1154 fq_nmod_init(alpha, ctx->fqctx);
1155 fq_nmod_init(alphastart, ctx->fqctx);
1156
1157 /* fail if the gcd has content wrt last variable */
1158 if (fq_nmod_poly_degree(c, ctx->fqctx) > 0)
1159 {
1160 success = 0;
1161 goto finished;
1162 }
1163
1164 /* we don't expect this function to work over F_p */
1165 if (nmod_poly_degree(ctx->fqctx->modulus) < WORD(2))
1166 {
1167 success = 0;
1168 goto finished;
1169 }
1170
1171 fq_nmod_poly_one(modulus, ctx->fqctx);
1172 fq_nmod_mpolyun_zero(H, ctx);
1173
1174 fq_nmod_randtest_not_zero(alphastart, randstate, ctx->fqctx);
1175 fq_nmod_set(alpha, alphastart, ctx->fqctx);
1176
1177 while (1)
1178 {
1179 /* get new evaluation point */
1180 fq_nmod_next_not_zero(alpha, ctx->fqctx);
1181 if (fq_nmod_equal(alpha, alphastart, ctx->fqctx))
1182 {
1183 success = 0;
1184 goto finished;
1185 }
1186
1187 /* make sure evaluation point does not kill both lc(A) and lc(B) */
1188 fq_nmod_poly_evaluate_fq_nmod(geval, g, alpha, ctx->fqctx);
1189 if (fq_nmod_is_zero(geval, ctx->fqctx))
1190 {
1191 goto outer_continue;
1192 }
1193
1194 /* make sure evaluation point does not kill either A or B */
1195 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Aeval, An, alpha, ctx);
1196 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Beval, Bn, alpha, ctx);
1197 if (Aeval->length == 0 || Beval->length == 0)
1198 {
1199 goto outer_continue;
1200 }
1201
1202 success = fq_nmod_mpolyu_gcdp_zippel(Geval, Abareval, Bbareval,
1203 Aeval, Beval, var - 1, ctx, zinfo, randstate);
1204 if (!success || Geval->exps[0] > degbound)
1205 {
1206 success = 0;
1207 goto finished;
1208 }
1209
1210 degbound = Geval->exps[0];
1211
1212 if (fq_nmod_mpolyu_is_one(Geval, ctx))
1213 {
1214 fq_nmod_mpolyu_one(G, ctx);
1215 fq_nmod_mpolyu_swap(Abar, A);
1216 fq_nmod_mpolyu_swap(Bbar, B);
1217 fq_nmod_mpolyu_shift_left(G, Gshift);
1218 fq_nmod_mpolyu_shift_left(Abar, Ashift - Gshift);
1219 fq_nmod_mpolyu_shift_left(Bbar, Bshift - Gshift);
1220 success = 1;
1221 goto finished;
1222 }
1223
1224 if (fq_nmod_poly_degree(modulus, ctx->fqctx) > 0)
1225 {
1226 if (Geval->exps[0] > H->exps[0])
1227 {
1228 goto outer_continue;
1229 }
1230 else if (Geval->exps[0] < H->exps[0])
1231 {
1232 fq_nmod_poly_one(modulus, ctx->fqctx);
1233 }
1234 }
1235
1236 /* update interpolant H */
1237 fq_nmod_inv(temp, fq_nmod_mpolyu_leadcoeff(Geval, ctx), ctx->fqctx);
1238 fq_nmod_mul(temp, geval, temp, ctx->fqctx);
1239 fq_nmod_mpolyu_scalar_mul_fq_nmod(Geval, temp, ctx);
1240 if (fq_nmod_poly_degree(modulus, ctx->fqctx) > 0)
1241 {
1242 fq_nmod_poly_evaluate_fq_nmod(temp, modulus, alpha, ctx->fqctx);
1243 fq_nmod_inv(temp, temp, ctx->fqctx);
1244 fq_nmod_poly_scalar_mul_fq_nmod(modulus, modulus, temp, ctx->fqctx);
1245 changed = fq_nmod_mpolyun_interp_crt_sm_mpolyu(&lastdeg, H, Ht, Geval,
1246 modulus, alpha, ctx);
1247 if (!changed)
1248 {
1249 fq_nmod_mpolyun_content_poly(a, H, ctx);
1250 fq_nmod_mpolyun_divexact_poly(Ht, H, a, ctx);
1251 fq_nmod_mpolyun_shift_left(Ht, Gshift);
1252 fq_nmod_mpolyu_cvtfrom_mpolyun(G, Ht, var, ctx);
1253 if ( !fq_nmod_mpolyuu_divides(Abar, A, G, 1, ctx)
1254 || !fq_nmod_mpolyuu_divides(Bbar, B, G, 1, ctx))
1255 {
1256 FLINT_ASSERT(0);
1257 goto outer_continue;
1258 }
1259 success = 1;
1260 goto finished;
1261 }
1262 }
1263 else
1264 {
1265 fq_nmod_mpolyun_interp_lift_sm_mpolyu(H, Geval, ctx);
1266 }
1267 fq_nmod_poly_set_coeff(tempmod, 0, alpha, ctx->fqctx);
1268 fq_nmod_poly_mul(modulus, modulus, tempmod, ctx->fqctx);
1269
1270 fq_nmod_mpolyu_setform_mpolyun(Gform, H, ctx);
1271
1272 while (1)
1273 {
1274 /* get new evaluation point */
1275 fq_nmod_next_not_zero(alpha, ctx->fqctx);
1276 if (fq_nmod_equal(alpha, alphastart, ctx->fqctx))
1277 {
1278 success = 0;
1279 goto finished;
1280 }
1281
1282 /* make sure evaluation point does not kill both lc(A) and lc(B) */
1283 fq_nmod_poly_evaluate_fq_nmod(geval, g, alpha, ctx->fqctx);
1284 if (fq_nmod_is_zero(geval, ctx->fqctx))
1285 {
1286 goto inner_continue;
1287 }
1288
1289 /* make sure evaluation point does not kill either A or B */
1290 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Aeval, An, alpha, ctx);
1291 fq_nmod_mpolyun_interp_reduce_sm_mpolyu(Beval, Bn, alpha, ctx);
1292 if (Aeval->length == 0 || Beval->length == 0)
1293 {
1294 goto inner_continue;
1295 }
1296
1297 switch (fq_nmod_mpolyu_gcds_zippel(Geval, Aeval, Beval, Gform, var,
1298 ctx, randstate, °bound))
1299 {
1300 default:
1301 FLINT_ASSERT(0);
1302 case nmod_gcds_form_main_degree_too_high:
1303 /* fq_nmod_mpolyu_gcds_zippel has updated degbound */
1304 fq_nmod_poly_one(modulus, ctx->fqctx);
1305 goto outer_continue;
1306 case nmod_gcds_form_wrong:
1307 case nmod_gcds_no_solution:
1308 success = 0;
1309 goto finished;
1310 case nmod_gcds_scales_not_found:
1311 case nmod_gcds_eval_point_not_found:
1312 case nmod_gcds_eval_gcd_deg_too_high:
1313 goto inner_continue;
1314 case nmod_gcds_success:
1315 (void)(NULL);
1316 }
1317
1318 if (fq_nmod_is_zero(fq_nmod_mpolyu_leadcoeff(Geval, ctx),
1319 ctx->fqctx))
1320 {
1321 goto inner_continue;
1322 }
1323
1324 /* update interpolant H */
1325 FLINT_ASSERT(fq_nmod_poly_degree(modulus, ctx->fqctx) > 0);
1326
1327 fq_nmod_inv(temp, fq_nmod_mpolyu_leadcoeff(Geval, ctx),
1328 ctx->fqctx);
1329 fq_nmod_mul(temp, temp, geval, ctx->fqctx);
1330 fq_nmod_mpolyu_scalar_mul_fq_nmod(Geval, temp, ctx);
1331
1332 fq_nmod_poly_evaluate_fq_nmod(temp, modulus, alpha, ctx->fqctx);
1333 fq_nmod_inv(temp, temp, ctx->fqctx);
1334 fq_nmod_poly_scalar_mul_fq_nmod(modulus, modulus, temp, ctx->fqctx);
1335 changed = fq_nmod_mpolyun_interp_crt_sm_mpolyu(&lastdeg, H, Ht, Geval, modulus,
1336 alpha, ctx);
1337
1338 fq_nmod_poly_set_coeff(tempmod, 0, alpha, ctx->fqctx);
1339 fq_nmod_poly_mul(modulus, modulus, tempmod, ctx->fqctx);
1340
1341 have_enough = fq_nmod_poly_degree(modulus, ctx->fqctx) >= bound;
1342
1343 if (changed && !have_enough)
1344 {
1345 goto inner_continue;
1346 }
1347
1348 if (!changed || have_enough)
1349 {
1350 fq_nmod_mpolyun_content_poly(a, H, ctx);
1351 fq_nmod_mpolyun_divexact_poly(Ht, H, a, ctx);
1352 fq_nmod_mpolyun_shift_left(Ht, Gshift);
1353 fq_nmod_mpolyu_cvtfrom_mpolyun(G, Ht, var, ctx);
1354 if ( fq_nmod_mpolyuu_divides(Abar, A, G, 1, ctx)
1355 && fq_nmod_mpolyuu_divides(Bbar, B, G, 1, ctx))
1356 {
1357 success = 1;
1358 goto finished;
1359 }
1360 }
1361
1362 if (have_enough)
1363 {
1364 fq_nmod_poly_one(modulus, ctx->fqctx);
1365 goto outer_continue;
1366 }
1367
1368 inner_continue:
1369 NULL;
1370 }
1371 FLINT_ASSERT(0 && "not reachable");
1372
1373 outer_continue:
1374 NULL;
1375 }
1376 FLINT_ASSERT(0 && "not reachable");
1377
1378
1379 finished:
1380
1381 fmpz_clear(minusone);
1382 fq_nmod_clear(geval, ctx->fqctx);
1383 fq_nmod_clear(temp, ctx->fqctx);
1384 fq_nmod_clear(alpha, ctx->fqctx);
1385 fq_nmod_clear(alphastart, ctx->fqctx);
1386 fq_nmod_poly_clear(a, ctx->fqctx);
1387 fq_nmod_poly_clear(b, ctx->fqctx);
1388 fq_nmod_poly_clear(c, ctx->fqctx);
1389 fq_nmod_poly_clear(g, ctx->fqctx);
1390 fq_nmod_poly_clear(modulus, ctx->fqctx);
1391 fq_nmod_poly_clear(tempmod, ctx->fqctx);
1392 fq_nmod_mpolyu_clear(Aeval, ctx);
1393 fq_nmod_mpolyu_clear(Beval, ctx);
1394 fq_nmod_mpolyu_clear(Geval, ctx);
1395 fq_nmod_mpolyu_clear(Abareval, ctx);
1396 fq_nmod_mpolyu_clear(Bbareval, ctx);
1397 fq_nmod_mpolyu_clear(Gform, ctx);
1398 fq_nmod_mpolyun_clear(An, ctx);
1399 fq_nmod_mpolyun_clear(Bn, ctx);
1400 fq_nmod_mpolyun_clear(H, ctx);
1401 fq_nmod_mpolyun_clear(Ht, ctx);
1402
1403 return success;
1404 }
1405