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 "n_poly.h"
13 #include "nmod_mpoly_factor.h"
14
n_polyu1n_mod_interp_reduce_2sm_poly(n_poly_t E,n_poly_t F,const n_polyun_t A,n_poly_t alphapow,nmod_t ctx)15 static void n_polyu1n_mod_interp_reduce_2sm_poly(
16 n_poly_t E,
17 n_poly_t F,
18 const n_polyun_t A,
19 n_poly_t alphapow,
20 nmod_t ctx)
21 {
22 slong i;
23 mp_limb_t u, v;
24
25 n_poly_zero(E);
26 n_poly_zero(F);
27 for (i = 0; i < A->length; i++)
28 {
29 n_poly_mod_eval2_pow(&u, &v, A->coeffs + i, alphapow, ctx);
30 n_poly_set_coeff(E, A->exps[i], u);
31 n_poly_set_coeff(F, A->exps[i], v);
32 }
33 }
34
n_polyu1n_mod_interp_lift_2sm_poly(slong * lastdeg,n_polyun_t F,const n_poly_t A,const n_poly_t B,mp_limb_t alpha,nmod_t ctx)35 static void n_polyu1n_mod_interp_lift_2sm_poly(
36 slong * lastdeg,
37 n_polyun_t F,
38 const n_poly_t A,
39 const n_poly_t B,
40 mp_limb_t alpha,
41 nmod_t ctx)
42 {
43 slong lastlen = 0;
44 slong Fi, Aexp, Bexp;
45 const mp_limb_t * Acoeffs = A->coeffs;
46 const mp_limb_t * Bcoeffs = B->coeffs;
47 slong e;
48 mp_limb_t d0 = (1 + ctx.n)/2;
49 mp_limb_t d1 = nmod_inv(nmod_add(alpha, alpha, ctx), ctx);
50 mp_limb_t Avalue, Bvalue, u, v;
51
52 Aexp = n_poly_degree(A);
53 Bexp = n_poly_degree(B);
54
55 n_polyun_fit_length(F, FLINT_MAX(Aexp, Bexp) + 1);
56
57 Fi = 0;
58 while (Aexp >= 0 || Bexp >= 0)
59 {
60 e = Aexp;
61 Avalue = 0;
62 Bvalue = 0;
63 if (Aexp == Bexp)
64 {
65 Avalue = Acoeffs[Aexp];
66 Bvalue = Bcoeffs[Bexp];
67 }
68 else if (Aexp > Bexp)
69 {
70 Avalue = Acoeffs[Aexp];
71 }
72 else
73 {
74 FLINT_ASSERT(Bexp > Aexp);
75 e = Bexp;
76 Bvalue = Bcoeffs[Bexp];
77 }
78 FLINT_ASSERT(Avalue != 0 || Bvalue != 0);
79 u = nmod_add(Avalue, Bvalue, ctx);
80 v = nmod_sub(Avalue, Bvalue, ctx);
81 u = nmod_mul(u, d0, ctx);
82 v = nmod_mul(v, d1, ctx);
83
84 FLINT_ASSERT(Fi < F->alloc);
85
86 F->exps[Fi] = e;
87
88 FLINT_ASSERT(u != 0 || v != 0);
89 n_poly_fit_length(F->coeffs + Fi, 2);
90 F->coeffs[Fi].coeffs[0] = u;
91 F->coeffs[Fi].coeffs[1] = v;
92 F->coeffs[Fi].length = 1 + (v != 0);
93 lastlen = FLINT_MAX(lastlen, F->coeffs[Fi].length);
94 Fi++;
95
96 if (e == Aexp)
97 {
98 do {
99 Aexp--;
100 } while (Aexp >= 0 && Acoeffs[Aexp] == 0);
101 }
102 if (e == Bexp)
103 {
104 do {
105 Bexp--;
106 } while (Bexp >= 0 && Bcoeffs[Bexp] == 0);
107 }
108 }
109 F->length = Fi;
110
111 *lastdeg = lastlen - 1;
112 return;
113 }
114
n_polyu1n_mod_interp_crt_2sm_poly(slong * lastdeg,n_polyun_t F,n_polyun_t T,const n_poly_t A,const n_poly_t B,const n_poly_t modulus,n_poly_t alphapow,nmod_t ctx)115 static int n_polyu1n_mod_interp_crt_2sm_poly(
116 slong * lastdeg,
117 n_polyun_t F,
118 n_polyun_t T,
119 const n_poly_t A,
120 const n_poly_t B,
121 const n_poly_t modulus,
122 n_poly_t alphapow,
123 nmod_t ctx)
124 {
125 int changed = 0, Finc;
126 slong lastlen = 0;
127 n_poly_struct * Fvalue;
128 mp_limb_t u, v, FvalueA, FvalueB;
129 slong Fi, Ti, Aexp, Bexp, e, fexp;
130 const mp_limb_t * Acoeff = A->coeffs;
131 const mp_limb_t * Bcoeff = B->coeffs;
132 slong Flen = F->length;
133 n_poly_t zero;
134
135 zero->alloc = 0;
136 zero->length = 0;
137 zero->coeffs = NULL;
138
139 Fi = 0;
140 Aexp = n_poly_degree(A);
141 Bexp = n_poly_degree(B);
142
143 n_polyun_fit_length(T, Flen + FLINT_MAX(Aexp, Bexp) + 1);
144 Ti = 0;
145
146 #if FLINT_WANT_ASSERT
147 u = n_poly_mod_evaluate_nmod(modulus, alphapow->coeffs[1], ctx);
148 u = nmod_mul(u, alphapow->coeffs[1], ctx);
149 u = nmod_add(u, u, ctx);
150 FLINT_ASSERT(u == 1);
151
152 v = nmod_neg(alphapow->coeffs[1], ctx);
153 u = n_poly_mod_evaluate_nmod(modulus, v, ctx);
154 u = nmod_mul(u, alphapow->coeffs[1], ctx);
155 u = nmod_add(u, u, ctx);
156 FLINT_ASSERT(u == 1);
157 #endif
158
159 while (Fi < Flen || Aexp >= 0 || Bexp >= 0)
160 {
161 FLINT_ASSERT(Ti < T->alloc);
162
163 fexp = e = -WORD(1);
164 if (Fi < Flen)
165 {
166 fexp = e = F->exps[Fi];
167 FLINT_ASSERT(!n_poly_is_zero(F->coeffs + Fi));
168 FLINT_ASSERT(n_poly_degree(F->coeffs + Fi) < n_poly_degree(modulus));
169 }
170 if (Aexp >= 0)
171 {
172 e = FLINT_MAX(e, Aexp);
173 FLINT_ASSERT(Acoeff[Aexp] != 0);
174 }
175 if (Bexp >= 0)
176 {
177 e = FLINT_MAX(e, Bexp);
178 FLINT_ASSERT(Bcoeff[Bexp] != 0);
179 }
180
181 FLINT_ASSERT(e >= 0);
182
183 T->exps[Ti] = e;
184
185 Fvalue = zero;
186 FvalueA = 0;
187 FvalueB = 0;
188 Finc = 0;
189 if (Fi < Flen && e == fexp)
190 {
191 Finc = 1;
192 Fvalue = F->coeffs + Fi;
193 n_poly_mod_eval2_pow(&FvalueA, &FvalueB, Fvalue, alphapow, ctx);
194 }
195
196 if (e == Aexp)
197 {
198 FvalueA = nmod_sub(FvalueA, Acoeff[Aexp], ctx);
199 }
200 if (e == Bexp)
201 {
202 FvalueB = nmod_sub(FvalueB, Bcoeff[Bexp], ctx);
203 }
204
205 u = nmod_sub(FvalueB, FvalueA, ctx);
206 v = nmod_add(FvalueB, FvalueA, ctx);
207 v = nmod_mul(v, alphapow->coeffs[1], ctx);
208 v = nmod_neg(v, ctx);
209 changed |= u != 0 || v != 0;
210 n_poly_mod_addmul_linear(T->coeffs + Ti, Fvalue, modulus, u, v, ctx);
211
212 FLINT_ASSERT(T->coeffs[Ti].length > 0);
213 lastlen = FLINT_MAX(lastlen, T->coeffs[Ti].length);
214 Ti++;
215
216 Fi += Finc;
217 if (e == Aexp)
218 {
219 do {
220 Aexp--;
221 } while (Aexp >= 0 && Acoeff[Aexp] == 0);
222 }
223 if (e == Bexp)
224 {
225 do {
226 Bexp--;
227 } while (Bexp >= 0 && Bcoeff[Bexp] == 0);
228 }
229 }
230 T->length = Ti;
231
232 *lastdeg = lastlen - 1;
233
234 if (changed)
235 n_polyun_swap(T, F);
236
237 return changed;
238 }
239
240
n_polyu1n_mod_gcd_brown_smprime(n_polyun_t G,n_polyun_t Abar,n_polyun_t Bbar,n_polyun_t A,n_polyun_t B,nmod_t ctx,n_poly_polyun_stack_t St)241 int n_polyu1n_mod_gcd_brown_smprime(
242 n_polyun_t G,
243 n_polyun_t Abar,
244 n_polyun_t Bbar,
245 n_polyun_t A,
246 n_polyun_t B,
247 nmod_t ctx,
248 n_poly_polyun_stack_t St)
249 {
250 int success;
251 slong bound;
252 mp_limb_t alpha, temp, gammaevalp, gammaevalm;
253 n_poly_struct * Aevalp, * Bevalp, * Gevalp, * Abarevalp, * Bbarevalp;
254 n_poly_struct * Aevalm, * Bevalm, * Gevalm, * Abarevalm, * Bbarevalm;
255 n_polyun_struct * T;
256 slong deggamma, ldegG, ldegAbar, ldegBbar, ldegA, ldegB;
257 n_poly_struct * cA, * cB, * cG, * cAbar, * cBbar, * gamma;
258 n_poly_struct * modulus, * alphapow, * r;
259 int gstab, astab, bstab, use_stab;
260 #if FLINT_WANT_ASSERT
261 n_poly_t leadA, leadB;
262
263 FLINT_ASSERT(A->length > 0);
264 FLINT_ASSERT(B->length > 0);
265
266 n_poly_init(leadA);
267 n_poly_init(leadB);
268 n_poly_set(leadA, A->coeffs + 0);
269 n_poly_set(leadB, B->coeffs + 0);
270 #endif
271
272 n_poly_stack_fit_request(St->poly_stack, 19);
273 cA = n_poly_stack_take_top(St->poly_stack);
274 cB = n_poly_stack_take_top(St->poly_stack);
275 cG = n_poly_stack_take_top(St->poly_stack);
276 cAbar = n_poly_stack_take_top(St->poly_stack);
277 cBbar = n_poly_stack_take_top(St->poly_stack);
278 gamma = n_poly_stack_take_top(St->poly_stack);
279 Aevalp = n_poly_stack_take_top(St->poly_stack);
280 Bevalp = n_poly_stack_take_top(St->poly_stack);
281 Gevalp = n_poly_stack_take_top(St->poly_stack);
282 Abarevalp = n_poly_stack_take_top(St->poly_stack);
283 Bbarevalp = n_poly_stack_take_top(St->poly_stack);
284 Aevalm = n_poly_stack_take_top(St->poly_stack);
285 Bevalm = n_poly_stack_take_top(St->poly_stack);
286 Gevalm = n_poly_stack_take_top(St->poly_stack);
287 Abarevalm = n_poly_stack_take_top(St->poly_stack);
288 Bbarevalm = n_poly_stack_take_top(St->poly_stack);
289 r = n_poly_stack_take_top(St->poly_stack);
290 alphapow = n_poly_stack_take_top(St->poly_stack);
291 modulus = n_poly_stack_take_top(St->poly_stack);
292
293 n_polyun_stack_fit_request(St->polyun_stack, 1);
294 T = n_polyun_stack_take_top(St->polyun_stack);
295
296 _n_poly_vec_mod_remove_content(cA, A->coeffs, A->length, ctx);
297 _n_poly_vec_mod_remove_content(cB, B->coeffs, B->length, ctx);
298
299 n_poly_mod_gcd(cG, cA, cB, ctx);
300 n_poly_mod_div(cAbar, cA, cG, ctx);
301 n_poly_mod_div(cBbar, cB, cG, ctx);
302
303 n_poly_mod_gcd(gamma, A->coeffs + 0, B->coeffs + 0, ctx);
304
305 ldegA = _n_poly_vec_max_degree(A->coeffs, A->length);
306 ldegB = _n_poly_vec_max_degree(B->coeffs, B->length);
307 deggamma = n_poly_degree(gamma);
308 bound = 1 + deggamma + FLINT_MAX(ldegA, ldegB);
309
310 n_poly_fit_length(alphapow, FLINT_MAX(WORD(3), bound + 1));
311 n_poly_one(modulus);
312
313 use_stab = 1;
314 gstab = bstab = astab = 0;
315
316 alpha = (ctx.n - 1)/2;
317
318 choose_prime: /* primes are v - alpha, v + alpha */
319
320 if (alpha < 2)
321 {
322 success = 0;
323 goto cleanup;
324 }
325
326 alpha--;
327
328 FLINT_ASSERT(0 < alpha && alpha < ctx.n/2);
329 FLINT_ASSERT(alphapow->alloc >= 2);
330 alphapow->length = 2;
331 alphapow->coeffs[0] = 1;
332 alphapow->coeffs[1] = alpha;
333
334 n_poly_mod_eval2_pow(&gammaevalp, &gammaevalm, gamma, alphapow, ctx);
335 if (gammaevalp == 0 || gammaevalm == 0)
336 goto choose_prime;
337
338 /* evaluation point should kill neither A nor B */
339 n_polyu1n_mod_interp_reduce_2sm_poly(Aevalp, Aevalm, A, alphapow, ctx);
340 n_polyu1n_mod_interp_reduce_2sm_poly(Bevalp, Bevalm, B, alphapow, ctx);
341 FLINT_ASSERT(Aevalp->length > 0);
342 FLINT_ASSERT(Aevalm->length > 0);
343 FLINT_ASSERT(Bevalp->length > 0);
344 FLINT_ASSERT(Bevalm->length > 0);
345
346 if (use_stab && gstab)
347 {
348 slong Gdeg;
349 n_polyu1n_mod_interp_reduce_2sm_poly(Gevalp, Gevalm, G, alphapow, ctx);
350 Gdeg = G->exps[0];
351 success = 1;
352 success = success && n_poly_degree(Gevalp) == Gdeg;
353 success = success && n_poly_degree(Gevalm) == Gdeg;
354 success = success && Gevalp->coeffs[Gdeg] == gammaevalp;
355 success = success && Gevalm->coeffs[Gdeg] == gammaevalm;
356 success = success && (n_poly_mod_divrem(Abarevalp, r,
357 Aevalp, Gevalp, ctx), r->length == 0);
358 success = success && (n_poly_mod_divrem(Abarevalm, r,
359 Aevalm, Gevalm, ctx), r->length == 0);
360 success = success && (n_poly_mod_divrem(Bbarevalp, r,
361 Bevalp, Gevalp, ctx), r->length == 0);
362 success = success && (n_poly_mod_divrem(Bbarevalm, r,
363 Bevalm, Gevalm, ctx), r->length == 0);
364 if (!success)
365 {
366 use_stab = 0;
367 n_poly_one(modulus);
368 goto choose_prime;
369 }
370
371 _n_poly_mod_scalar_mul_nmod_inplace(Abarevalp, gammaevalp, ctx);
372 _n_poly_mod_scalar_mul_nmod_inplace(Abarevalm, gammaevalm, ctx);
373 _n_poly_mod_scalar_mul_nmod_inplace(Bbarevalp, gammaevalp, ctx);
374 _n_poly_mod_scalar_mul_nmod_inplace(Bbarevalm, gammaevalm, ctx);
375 }
376 else
377 {
378 n_poly_mod_gcd(Gevalp, Aevalp, Bevalp, ctx);
379 n_poly_mod_div(Abarevalp, Aevalp, Gevalp, ctx);
380 n_poly_mod_div(Bbarevalp, Bevalp, Gevalp, ctx);
381 n_poly_mod_gcd(Gevalm, Aevalm, Bevalm, ctx);
382 n_poly_mod_div(Abarevalm, Aevalm, Gevalm, ctx);
383 n_poly_mod_div(Bbarevalm, Bevalm, Gevalm, ctx);
384 gstab = astab = bstab = 0;
385 }
386
387 FLINT_ASSERT(Gevalp->length > 0);
388 FLINT_ASSERT(Abarevalp->length > 0);
389 FLINT_ASSERT(Bbarevalp->length > 0);
390 FLINT_ASSERT(Gevalm->length > 0);
391 FLINT_ASSERT(Abarevalm->length > 0);
392 FLINT_ASSERT(Bbarevalm->length > 0);
393
394 if (n_poly_degree(Gevalp) == 0 || n_poly_degree(Gevalm) == 0)
395 {
396 n_polyun_one(G);
397 n_polyun_swap(Abar, A);
398 n_polyun_swap(Bbar, B);
399 goto successful_put_content;
400 }
401
402 if (n_poly_degree(Gevalp) != n_poly_degree(Gevalm))
403 {
404 goto choose_prime;
405 }
406
407 /* the Geval have matching degrees */
408 if (n_poly_degree(modulus) > 0)
409 {
410 FLINT_ASSERT(G->length > 0);
411 if (n_poly_degree(Gevalp) > G->exps[0])
412 {
413 goto choose_prime;
414 }
415 else if (n_poly_degree(Gevalp) < G->exps[0])
416 {
417 n_poly_one(modulus);
418 }
419 }
420
421 /* update interpolants */
422 _n_poly_mod_scalar_mul_nmod_inplace(Gevalp, gammaevalp, ctx);
423 _n_poly_mod_scalar_mul_nmod_inplace(Gevalm, gammaevalm, ctx);
424 if (n_poly_degree(modulus) > 0)
425 {
426 temp = n_poly_mod_evaluate_nmod(modulus, alpha, ctx);
427 FLINT_ASSERT(temp == n_poly_mod_evaluate_nmod(modulus, ctx.n - alpha, ctx));
428 temp = nmod_mul(temp, alpha, ctx);
429 temp = nmod_add(temp, temp, ctx);
430 temp = nmod_inv(temp, ctx);
431 _n_poly_mod_scalar_mul_nmod_inplace(modulus, temp, ctx);
432 gstab = gstab || !n_polyu1n_mod_interp_crt_2sm_poly(&ldegG, G, T, Gevalp, Gevalm, modulus, alphapow, ctx);
433 n_polyu1n_mod_interp_crt_2sm_poly(&ldegAbar, Abar, T, Abarevalp, Abarevalm, modulus, alphapow, ctx);
434 n_polyu1n_mod_interp_crt_2sm_poly(&ldegBbar, Bbar, T, Bbarevalp, Bbarevalm, modulus, alphapow, ctx);
435 }
436 else
437 {
438 n_polyu1n_mod_interp_lift_2sm_poly(&ldegG, G, Gevalp, Gevalm, alpha, ctx);
439 n_polyu1n_mod_interp_lift_2sm_poly(&ldegAbar, Abar, Abarevalp, Abarevalm, alpha, ctx);
440 n_polyu1n_mod_interp_lift_2sm_poly(&ldegBbar, Bbar, Bbarevalp, Bbarevalm, alpha, ctx);
441 gstab = astab = bstab = 0;
442 }
443
444 temp = ctx.n - nmod_mul(alpha, alpha, ctx);
445 n_poly_mod_shift_left_scalar_addmul(modulus, 2, temp, ctx);
446
447 if (n_poly_degree(modulus) < bound)
448 goto choose_prime;
449
450 FLINT_ASSERT(ldegG >= 0);
451 FLINT_ASSERT(ldegAbar >= 0);
452 FLINT_ASSERT(ldegBbar >= 0);
453
454 if (deggamma + ldegA == ldegG + ldegAbar &&
455 deggamma + ldegB == ldegG + ldegBbar)
456 {
457 goto successful;
458 }
459
460 n_poly_one(modulus);
461 goto choose_prime;
462
463 successful:
464
465 _n_poly_vec_mod_remove_content(modulus, G->coeffs, G->length, ctx);
466 _n_poly_vec_mod_divexact_poly(Abar->coeffs, Abar->length, G->coeffs + 0, ctx);
467 _n_poly_vec_mod_divexact_poly(Bbar->coeffs, Bbar->length, G->coeffs + 0, ctx);
468
469 successful_put_content:
470
471 _n_poly_vec_mod_mul_poly(G->coeffs, G->length, cG, ctx);
472 _n_poly_vec_mod_mul_poly(Abar->coeffs, Abar->length, cAbar, ctx);
473 _n_poly_vec_mod_mul_poly(Bbar->coeffs, Bbar->length, cBbar, ctx);
474
475 success = 1;
476
477 cleanup:
478
479 #if FLINT_WANT_ASSERT
480 if (success)
481 {
482 FLINT_ASSERT(1 == n_poly_lead(G->coeffs + 0));
483 n_poly_mod_mul(modulus, G->coeffs + 0, Abar->coeffs + 0, ctx);
484 FLINT_ASSERT(n_poly_equal(modulus, leadA));
485 n_poly_mod_mul(modulus, G->coeffs + 0, Bbar->coeffs + 0, ctx);
486 FLINT_ASSERT(n_poly_equal(modulus, leadB));
487 }
488 n_poly_clear(leadA);
489 n_poly_clear(leadB);
490 #endif
491
492 n_poly_stack_give_back(St->poly_stack, 19);
493 n_polyun_stack_give_back(St->polyun_stack, 1);
494
495 return success;
496 }
497
498