Line data Source code
1 : /* Copyright (C) 2000-2005 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /***********************************************************************/
16 : /** **/
17 : /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
18 : /** (third part) **/
19 : /** **/
20 : /***********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_pol
25 :
26 : /************************************************************************
27 : ** **
28 : ** Ring membership **
29 : ** **
30 : ************************************************************************/
31 : struct charact {
32 : GEN q;
33 : int isprime;
34 : };
35 : static void
36 1239 : char_update_prime(struct charact *S, GEN p)
37 : {
38 1239 : if (!S->isprime) { S->isprime = 1; S->q = p; }
39 1239 : if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
40 1232 : }
41 : static void
42 6573 : char_update_int(struct charact *S, GEN n)
43 : {
44 6573 : if (S->isprime)
45 : {
46 7 : if (dvdii(n, S->q)) return;
47 7 : pari_err_MODULUS("characteristic", S->q, n);
48 : }
49 6566 : S->q = gcdii(S->q, n);
50 : }
51 : static void
52 686889 : charact(struct charact *S, GEN x)
53 : {
54 686889 : const long tx = typ(x);
55 : long i, l;
56 686889 : switch(tx)
57 : {
58 5124 : case t_INTMOD:char_update_int(S, gel(x,1)); break;
59 1148 : case t_FFELT: char_update_prime(S, gel(x,4)); break;
60 21665 : case t_COMPLEX: case t_QUAD:
61 : case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
62 : case t_VEC: case t_COL: case t_MAT:
63 21665 : l = lg(x);
64 692048 : for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
65 21651 : break;
66 7 : case t_LIST:
67 7 : x = list_data(x);
68 7 : if (x) charact(S, x);
69 7 : break;
70 : }
71 686861 : }
72 : static void
73 4634 : charact_res(struct charact *S, GEN x)
74 : {
75 4634 : const long tx = typ(x);
76 : long i, l;
77 4634 : switch(tx)
78 : {
79 1449 : case t_INTMOD:char_update_int(S, gel(x,1)); break;
80 0 : case t_FFELT: char_update_prime(S, gel(x,4)); break;
81 91 : case t_PADIC: char_update_prime(S, gel(x,2)); break;
82 1617 : case t_COMPLEX: case t_QUAD:
83 : case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
84 : case t_VEC: case t_COL: case t_MAT:
85 1617 : l = lg(x);
86 5922 : for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
87 1617 : break;
88 0 : case t_LIST:
89 0 : x = list_data(x);
90 0 : if (x) charact_res(S, x);
91 0 : break;
92 : }
93 4634 : }
94 : GEN
95 16492 : characteristic(GEN x)
96 : {
97 : struct charact S;
98 16492 : S.q = gen_0; S.isprime = 0;
99 16492 : charact(&S, x); return S.q;
100 : }
101 : GEN
102 329 : residual_characteristic(GEN x)
103 : {
104 : struct charact S;
105 329 : S.q = gen_0; S.isprime = 0;
106 329 : charact_res(&S, x); return S.q;
107 : }
108 :
109 : int
110 66474320 : Rg_is_Fp(GEN x, GEN *pp)
111 : {
112 : GEN mod;
113 66474320 : switch(typ(x))
114 : {
115 4259766 : case t_INTMOD:
116 4259766 : mod = gel(x,1);
117 4259766 : if (!*pp) *pp = mod;
118 4024244 : else if (mod != *pp && !equalii(mod, *pp))
119 : {
120 0 : if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_Fp");
121 0 : return 0;
122 : }
123 4259766 : return 1;
124 53924442 : case t_INT:
125 53924442 : return 1;
126 8290112 : default: return 0;
127 : }
128 : }
129 :
130 : int
131 24797489 : RgX_is_FpX(GEN x, GEN *pp)
132 : {
133 24797489 : long i, lx = lg(x);
134 82955888 : for (i=2; i<lx; i++)
135 66448503 : if (!Rg_is_Fp(gel(x, i), pp))
136 8290124 : return 0;
137 16507385 : return 1;
138 : }
139 :
140 : int
141 0 : RgV_is_FpV(GEN x, GEN *pp)
142 : {
143 0 : long i, lx = lg(x);
144 0 : for (i=1; i<lx; i++)
145 0 : if (!Rg_is_Fp(gel(x,i), pp)) return 0;
146 0 : return 1;
147 : }
148 :
149 : int
150 0 : RgM_is_FpM(GEN x, GEN *pp)
151 : {
152 0 : long i, lx = lg(x);
153 0 : for (i=1; i<lx; i++)
154 0 : if (!RgV_is_FpV(gel(x, i), pp)) return 0;
155 0 : return 1;
156 : }
157 :
158 : int
159 58485 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
160 : {
161 : GEN pol, mod, p;
162 58485 : switch(typ(x))
163 : {
164 25795 : case t_INTMOD:
165 25795 : return Rg_is_Fp(x, pp);
166 6587 : case t_INT:
167 6587 : return 1;
168 21 : case t_POL:
169 21 : return RgX_is_FpX(x, pp);
170 21350 : case t_FFELT:
171 21350 : mod = x; p = FF_p_i(x);
172 21350 : if (!*pp) *pp = p;
173 21350 : if (!*pT) *pT = mod;
174 19824 : else if (typ(*pT)!=t_FFELT || !FF_samefield(*pT,mod))
175 : {
176 42 : if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
177 42 : return 0;
178 : }
179 21308 : return 1;
180 4578 : case t_POLMOD:
181 4578 : mod = gel(x,1); pol = gel(x, 2);
182 4578 : if (!RgX_is_FpX(mod, pp)) return 0;
183 4578 : if (typ(pol)==t_POL)
184 : {
185 4571 : if (!RgX_is_FpX(pol, pp)) return 0;
186 : }
187 7 : else if (!Rg_is_Fp(pol, pp)) return 0;
188 4578 : if (!*pT) *pT = mod;
189 4431 : else if (mod != *pT && !gequal(mod, *pT))
190 : {
191 0 : if (DEBUGLEVEL) pari_warn(warner,"different moduli in Rg_is_FpXQ");
192 0 : return 0;
193 : }
194 4578 : return 1;
195 :
196 154 : default: return 0;
197 : }
198 : }
199 :
200 : int
201 3115 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
202 : {
203 3115 : long i, lx = lg(x);
204 60956 : for (i = 2; i < lx; i++)
205 57939 : if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
206 3017 : return 1;
207 : }
208 :
209 : /************************************************************************
210 : ** **
211 : ** Ring conversion **
212 : ** **
213 : ************************************************************************/
214 :
215 : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
216 : * If x is an INTMOD, assume modulus is a multiple of p. */
217 : GEN
218 34155664 : Rg_to_Fp(GEN x, GEN p)
219 : {
220 34155664 : if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
221 3102270 : switch(typ(x))
222 : {
223 213392 : case t_INT: return modii(x, p);
224 121 : case t_FRAC: {
225 121 : pari_sp av = avma;
226 121 : GEN z = modii(gel(x,1), p);
227 121 : if (z == gen_0) return gen_0;
228 121 : return gerepileuptoint(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
229 : }
230 0 : case t_PADIC: return padic_to_Fp(x, p);
231 2888778 : case t_INTMOD: {
232 2888778 : GEN q = gel(x,1), a = gel(x,2);
233 2888778 : if (equalii(q, p)) return icopy(a);
234 14 : if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
235 0 : return remii(a, p);
236 : }
237 0 : default: pari_err_TYPE("Rg_to_Fp",x);
238 : return NULL; /* LCOV_EXCL_LINE */
239 : }
240 : }
241 : /* If x is a POLMOD, assume modulus is a multiple of T. */
242 : GEN
243 1292383 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
244 : {
245 1292383 : long ta, tx = typ(x), v = get_FpX_var(T);
246 : GEN a, b;
247 1292383 : if (is_const_t(tx))
248 : {
249 57913 : if (tx == t_FFELT)
250 : {
251 17085 : GEN z = FF_to_FpXQ(x);
252 17085 : setvarn(z, v);
253 17085 : return z;
254 : }
255 40828 : return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
256 : }
257 1234470 : switch(tx)
258 : {
259 1229619 : case t_POLMOD:
260 1229619 : b = gel(x,1);
261 1229619 : a = gel(x,2); ta = typ(a);
262 1229619 : if (is_const_t(ta))
263 4025 : return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
264 1225594 : b = RgX_to_FpX(b, p); if (varn(b) != v) break;
265 1225594 : a = RgX_to_FpX(a, p);
266 1225594 : if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
267 1225594 : return FpX_rem(a, T, p);
268 0 : break;
269 4851 : case t_POL:
270 4851 : if (varn(x) != v) break;
271 4851 : return FpX_rem(RgX_to_FpX(x,p), T, p);
272 0 : case t_RFRAC:
273 0 : a = Rg_to_FpXQ(gel(x,1), T,p);
274 0 : b = Rg_to_FpXQ(gel(x,2), T,p);
275 0 : return FpXQ_div(a,b, T,p);
276 : }
277 0 : pari_err_TYPE("Rg_to_FpXQ",x);
278 : return NULL; /* LCOV_EXCL_LINE */
279 : }
280 : GEN
281 3440007 : RgX_to_FpX(GEN x, GEN p)
282 : {
283 : long i, l;
284 3440007 : GEN z = cgetg_copy(x, &l); z[1] = x[1];
285 15143955 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
286 3440007 : return FpX_renormalize(z, l);
287 : }
288 :
289 : GEN
290 126 : RgV_to_FpV(GEN x, GEN p)
291 385 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
292 :
293 : GEN
294 940222 : RgC_to_FpC(GEN x, GEN p)
295 11883318 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
296 :
297 : GEN
298 132141 : RgM_to_FpM(GEN x, GEN p)
299 1072321 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
300 :
301 : GEN
302 285907 : RgV_to_Flv(GEN x, ulong p)
303 1172280 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
304 :
305 : GEN
306 114124 : RgM_to_Flm(GEN x, ulong p)
307 392639 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
308 :
309 : GEN
310 5000 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
311 : {
312 5000 : long i, l = lg(x);
313 5000 : GEN z = cgetg(l, t_POL); z[1] = x[1];
314 42855 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
315 5000 : return FpXQX_renormalize(z, l);
316 : }
317 : GEN
318 64512 : RgX_to_FqX(GEN x, GEN T, GEN p)
319 : {
320 64512 : long i, l = lg(x);
321 64512 : GEN z = cgetg(l, t_POL); z[1] = x[1];
322 64512 : if (T)
323 10990 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
324 : else
325 1163750 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
326 64512 : return FpXQX_renormalize(z, l);
327 : }
328 :
329 : GEN
330 219184 : RgC_to_FqC(GEN x, GEN T, GEN p)
331 : {
332 219184 : long i, l = lg(x);
333 219184 : GEN z = cgetg(l, t_COL);
334 219184 : if (T)
335 1431710 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
336 : else
337 0 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
338 219184 : return z;
339 : }
340 :
341 : GEN
342 52332 : RgM_to_FqM(GEN x, GEN T, GEN p)
343 271502 : { pari_APPLY_same(RgC_to_FqC(gel(x, i), T, p)) }
344 :
345 : /* lg(V) > 1 */
346 : GEN
347 850311 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
348 : {
349 850311 : pari_sp av = avma;
350 850311 : long i, l = lg(V);
351 850311 : GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
352 4188597 : for(i=2; i<l; i++)
353 : {
354 3338286 : z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
355 3338286 : if ((i & 7) == 0) z = gerepileupto(av, z);
356 : }
357 850311 : return gerepileupto(av, FpX_red(z,p));
358 : }
359 :
360 : GEN
361 52976 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
362 : {
363 52976 : long i, lz = lg(y);
364 : GEN z;
365 52976 : if (!T) return FpX_Fp_add(y, x, p);
366 5012 : if (lz == 2) return scalarpol(x, varn(y));
367 4599 : z = cgetg(lz,t_POL); z[1] = y[1];
368 4599 : gel(z,2) = Fq_add(gel(y,2),x, T, p);
369 4599 : if (lz == 3) z = FpXX_renormalize(z,lz);
370 : else
371 18725 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
372 4599 : return z;
373 : }
374 :
375 : GEN
376 1048 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
377 : {
378 1048 : long i, lz = lg(y);
379 : GEN z;
380 1048 : if (!T) return FpX_Fp_sub(y, x, p);
381 1048 : if (lz == 2) return scalarpol(x, varn(y));
382 1048 : z = cgetg(lz,t_POL); z[1] = y[1];
383 1048 : gel(z,2) = Fq_sub(gel(y,2), x, T, p);
384 1048 : if (lz == 3) z = FpXX_renormalize(z,lz);
385 : else
386 10288 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
387 1048 : return z;
388 : }
389 :
390 : GEN
391 148839 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
392 : {
393 : long i, lP;
394 148839 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
395 918219 : for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
396 148839 : gel(res,lP-1) = gen_1; return res;
397 : }
398 :
399 : GEN
400 37847 : FpXQX_normalize(GEN z, GEN T, GEN p)
401 : {
402 : GEN lc;
403 37847 : if (lg(z) == 2) return z;
404 37833 : lc = leading_coeff(z);
405 37833 : if (typ(lc) == t_POL)
406 : {
407 1982 : if (lg(lc) > 3) /* nonconstant */
408 1712 : return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
409 : /* constant */
410 270 : lc = gel(lc,2);
411 270 : z = shallowcopy(z);
412 270 : gel(z, lg(z)-1) = lc;
413 : }
414 : /* lc a t_INT */
415 36121 : if (equali1(lc)) return z;
416 64 : return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
417 : }
418 :
419 : GEN
420 159698 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
421 : {
422 : pari_sp av;
423 : GEN p1, r;
424 159698 : long j, i=lg(x)-1;
425 159698 : if (i<=2)
426 26327 : return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
427 133371 : av=avma; p1=gel(x,i);
428 : /* specific attention to sparse polynomials (see poleval)*/
429 : /*You've guessed it! It's a copy-paste(tm)*/
430 445700 : for (i--; i>=2; i=j-1)
431 : {
432 312833 : for (j=i; !signe(gel(x,j)); j--)
433 504 : if (j==2)
434 : {
435 315 : if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
436 315 : return gerepileupto(av, Fq_mul(p1,y, T, p));
437 : }
438 312329 : r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
439 312329 : p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
440 : }
441 133056 : return gerepileupto(av, p1);
442 : }
443 :
444 : GEN
445 31598 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
446 : {
447 31598 : long i, lb = lg(Q);
448 : GEN z;
449 31598 : if (!T) return FpXY_evalx(Q, x, p);
450 20993 : z = cgetg(lb, t_POL); z[1] = Q[1];
451 117068 : for (i=2; i<lb; i++)
452 : {
453 96075 : GEN q = gel(Q,i);
454 96075 : gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
455 : }
456 20993 : return FpXQX_renormalize(z, lb);
457 : }
458 :
459 : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
460 : GEN
461 14392 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
462 : {
463 14392 : pari_sp av = avma;
464 14392 : if (!T) return FpXY_eval(Q, y, x, p);
465 588 : return gerepileupto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
466 : }
467 :
468 : /* a X^d */
469 : GEN
470 8810262 : monomial(GEN a, long d, long v)
471 : {
472 : long i, n;
473 : GEN P;
474 8810262 : if (d < 0) {
475 14 : if (isrationalzero(a)) return pol_0(v);
476 14 : retmkrfrac(a, pol_xn(-d, v));
477 : }
478 8810248 : if (gequal0(a))
479 : {
480 19894 : if (isexactzero(a)) return scalarpol_shallow(a,v);
481 0 : n = d+2; P = cgetg(n+1, t_POL);
482 0 : P[1] = evalsigne(0) | evalvarn(v);
483 : }
484 : else
485 : {
486 8790356 : n = d+2; P = cgetg(n+1, t_POL);
487 8790355 : P[1] = evalsigne(1) | evalvarn(v);
488 : }
489 25432324 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
490 8790355 : gel(P,i) = a; return P;
491 : }
492 : GEN
493 1861042 : monomialcopy(GEN a, long d, long v)
494 : {
495 : long i, n;
496 : GEN P;
497 1861042 : if (d < 0) {
498 14 : if (isrationalzero(a)) return pol_0(v);
499 14 : retmkrfrac(gcopy(a), pol_xn(-d, v));
500 : }
501 1861028 : if (gequal0(a))
502 : {
503 7 : if (isexactzero(a)) return scalarpol(a,v);
504 0 : n = d+2; P = cgetg(n+1, t_POL);
505 0 : P[1] = evalsigne(0) | evalvarn(v);
506 : }
507 : else
508 : {
509 1861021 : n = d+2; P = cgetg(n+1, t_POL);
510 1861021 : P[1] = evalsigne(1) | evalvarn(v);
511 : }
512 3505432 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
513 1861021 : gel(P,i) = gcopy(a); return P;
514 : }
515 : GEN
516 322981 : pol_x_powers(long N, long v)
517 : {
518 322981 : GEN L = cgetg(N+1,t_VEC);
519 : long i;
520 780135 : for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
521 322985 : return L;
522 : }
523 :
524 : GEN
525 0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
526 : {
527 0 : return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
528 : }
529 :
530 : GEN
531 0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
532 : {
533 0 : return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
534 : }
535 :
536 : /*******************************************************************/
537 : /* */
538 : /* Fq */
539 : /* */
540 : /*******************************************************************/
541 :
542 : GEN
543 4555245 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
544 : {
545 : (void)T;
546 4555245 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
547 : {
548 541302 : case 0: return Fp_add(x,y,p);
549 205282 : case 1: return FpX_Fp_add(x,y,p);
550 57919 : case 2: return FpX_Fp_add(y,x,p);
551 3750742 : case 3: return FpX_add(x,y,p);
552 : }
553 : return NULL;/*LCOV_EXCL_LINE*/
554 : }
555 :
556 : GEN
557 4781769 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
558 : {
559 : (void)T;
560 4781769 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
561 : {
562 228390 : case 0: return Fp_sub(x,y,p);
563 8601 : case 1: return FpX_Fp_sub(x,y,p);
564 10149 : case 2: return Fp_FpX_sub(x,y,p);
565 4534629 : case 3: return FpX_sub(x,y,p);
566 : }
567 : return NULL;/*LCOV_EXCL_LINE*/
568 : }
569 :
570 : GEN
571 770453 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
572 : {
573 : (void)T;
574 770453 : return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
575 : }
576 :
577 : GEN
578 12829 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
579 : {
580 : (void)T;
581 12829 : return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
582 : }
583 :
584 : /* If T==NULL do not reduce*/
585 : GEN
586 5961700 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
587 : {
588 5961700 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
589 : {
590 665981 : case 0: return Fp_mul(x,y,p);
591 66909 : case 1: return FpX_Fp_mul(x,y,p);
592 155475 : case 2: return FpX_Fp_mul(y,x,p);
593 5073337 : case 3: if (T) return FpXQ_mul(x,y,T,p);
594 2832902 : else return FpX_mul(x,y,p);
595 : }
596 : return NULL;/*LCOV_EXCL_LINE*/
597 : }
598 :
599 : /* If T==NULL do not reduce*/
600 : GEN
601 395643 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
602 : {
603 : (void) T;
604 395643 : return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
605 : }
606 :
607 : /* y t_INT */
608 : GEN
609 33291 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
610 : {
611 : (void)T;
612 5830 : return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
613 39121 : : Fp_mul(x,y,p);
614 : }
615 : /* If T==NULL do not reduce*/
616 : GEN
617 293214 : Fq_sqr(GEN x, GEN T, GEN p)
618 : {
619 293214 : if (typ(x) == t_POL)
620 : {
621 12047 : if (T) return FpXQ_sqr(x,T,p);
622 0 : else return FpX_sqr(x,p);
623 : }
624 : else
625 281167 : return Fp_sqr(x,p);
626 : }
627 :
628 : GEN
629 0 : Fq_neg_inv(GEN x, GEN T, GEN p)
630 : {
631 0 : if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
632 0 : return FpXQ_inv(FpX_neg(x,p),T,p);
633 : }
634 :
635 : GEN
636 0 : Fq_invsafe(GEN x, GEN pol, GEN p)
637 : {
638 0 : if (typ(x) == t_INT) return Fp_invsafe(x,p);
639 0 : return FpXQ_invsafe(x,pol,p);
640 : }
641 :
642 : GEN
643 67971 : Fq_inv(GEN x, GEN pol, GEN p)
644 : {
645 67971 : if (typ(x) == t_INT) return Fp_inv(x,p);
646 61716 : return FpXQ_inv(x,pol,p);
647 : }
648 :
649 : GEN
650 303009 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
651 : {
652 303009 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
653 : {
654 287245 : case 0: return Fp_div(x,y,p);
655 10241 : case 1: return FpX_Fp_div(x,y,p);
656 280 : case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
657 5243 : case 3: return FpXQ_div(x,y,pol,p);
658 : }
659 : return NULL;/*LCOV_EXCL_LINE*/
660 : }
661 :
662 : GEN
663 694995 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
664 : {
665 694995 : if (typ(x) == t_INT) return Fp_pow(x,n,p);
666 130514 : return FpXQ_pow(x,n,pol,p);
667 : }
668 :
669 : GEN
670 14693 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
671 : {
672 14693 : if (typ(x) == t_INT) return Fp_powu(x,n,p);
673 749 : return FpXQ_powu(x,n,pol,p);
674 : }
675 :
676 : GEN
677 709980 : Fq_sqrt(GEN x, GEN T, GEN p)
678 : {
679 709980 : if (typ(x) == t_INT)
680 : {
681 699160 : if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
682 336 : x = scalarpol_shallow(x, get_FpX_var(T));
683 : }
684 11156 : return FpXQ_sqrt(x,T,p);
685 : }
686 : GEN
687 60123 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
688 : {
689 60123 : if (typ(x) == t_INT)
690 : {
691 : long d;
692 59857 : if (!T) return Fp_sqrtn(x,n,p,zeta);
693 77 : d = get_FpX_degree(T);
694 77 : if (ugcdiu(n,d) == 1)
695 : {
696 14 : if (!zeta) return Fp_sqrtn(x,n,p,NULL);
697 : /* gcd(n,p-1)=gcd(n,q-1): same number of solutions in Fp and F_q */
698 7 : if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
699 7 : return Fp_sqrtn(x,n,p,zeta);
700 : }
701 63 : x = scalarpol(x, get_FpX_var(T)); /* left on stack */
702 : }
703 329 : return FpXQ_sqrtn(x,n,T,p,zeta);
704 : }
705 :
706 : struct _Fq_field
707 : {
708 : GEN T, p;
709 : };
710 :
711 : static GEN
712 302463 : _Fq_red(void *E, GEN x)
713 302463 : { struct _Fq_field *s = (struct _Fq_field *)E;
714 302463 : return Fq_red(x, s->T, s->p);
715 : }
716 :
717 : static GEN
718 357525 : _Fq_add(void *E, GEN x, GEN y)
719 : {
720 : (void) E;
721 357525 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
722 : {
723 14 : case 0: return addii(x,y);
724 0 : case 1: return ZX_Z_add(x,y);
725 15918 : case 2: return ZX_Z_add(y,x);
726 341593 : default: return ZX_add(x,y);
727 : }
728 : }
729 :
730 : static GEN
731 315028 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
732 :
733 : static GEN
734 237559 : _Fq_mul(void *E, GEN x, GEN y)
735 : {
736 : (void) E;
737 237559 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
738 : {
739 133 : case 0: return mulii(x,y);
740 1085 : case 1: return ZX_Z_mul(x,y);
741 1043 : case 2: return ZX_Z_mul(y,x);
742 235298 : default: return ZX_mul(x,y);
743 : }
744 : }
745 :
746 : static GEN
747 9331 : _Fq_inv(void *E, GEN x)
748 9331 : { struct _Fq_field *s = (struct _Fq_field *)E;
749 9331 : return Fq_inv(x,s->T,s->p);
750 : }
751 :
752 : static int
753 38388 : _Fq_equal0(GEN x) { return signe(x)==0; }
754 :
755 : static GEN
756 13965 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
757 :
758 : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
759 : _Fq_inv,_Fq_equal0,_Fq_s};
760 :
761 4165 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
762 : {
763 4165 : if (!T)
764 0 : return get_Fp_field(E, p);
765 : else
766 : {
767 4165 : GEN z = new_chunk(sizeof(struct _Fq_field));
768 4165 : struct _Fq_field *e = (struct _Fq_field *) z;
769 4165 : e->T = T; e->p = p; *E = (void*)e;
770 4165 : return &Fq_field;
771 : }
772 : }
773 :
774 : /*******************************************************************/
775 : /* */
776 : /* Fq[X] */
777 : /* */
778 : /*******************************************************************/
779 : /* P(X + c) */
780 : GEN
781 266 : FpX_translate(GEN P, GEN c, GEN p)
782 : {
783 266 : pari_sp av = avma;
784 : GEN Q, *R;
785 : long i, k, n;
786 :
787 266 : if (!signe(P) || !signe(c)) return ZX_copy(P);
788 266 : Q = leafcopy(P);
789 266 : R = (GEN*)(Q+2); n = degpol(P);
790 3738 : for (i=1; i<=n; i++)
791 : {
792 118153 : for (k=n-i; k<n; k++)
793 114681 : R[k] = Fp_add(R[k], Fp_mul(c, R[k+1], p), p);
794 :
795 3472 : if (gc_needed(av,2))
796 : {
797 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FpX_translate, i = %ld/%ld", i,n);
798 0 : Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
799 : }
800 : }
801 266 : return gerepilecopy(av, FpX_renormalize(Q, lg(Q)));
802 : }
803 : /* P(X + c), c an Fq */
804 : GEN
805 33880 : FqX_translate(GEN P, GEN c, GEN T, GEN p)
806 : {
807 33880 : pari_sp av = avma;
808 : GEN Q, *R;
809 : long i, k, n;
810 :
811 : /* signe works for t_(INT|POL) */
812 33880 : if (!signe(P) || !signe(c)) return RgX_copy(P);
813 33880 : Q = leafcopy(P);
814 33880 : R = (GEN*)(Q+2); n = degpol(P);
815 150059 : for (i=1; i<=n; i++)
816 : {
817 433559 : for (k=n-i; k<n; k++)
818 317380 : R[k] = Fq_add(R[k], Fq_mul(c, R[k+1], T, p), T, p);
819 :
820 116179 : if (gc_needed(av,2))
821 : {
822 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FqX_translate, i = %ld/%ld", i,n);
823 0 : Q = gerepilecopy(av, Q); R = (GEN*)Q+2;
824 : }
825 : }
826 33880 : return gerepilecopy(av, FpXQX_renormalize(Q, lg(Q)));
827 : }
828 :
829 : GEN
830 40548 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
831 : {
832 40548 : pari_sp ltop = avma;
833 : long k;
834 : GEN W;
835 40548 : if (lgefint(p) == 3)
836 : {
837 31938 : ulong pp = p[2];
838 31938 : GEN Tl = ZX_to_Flx(T, pp);
839 31937 : GEN Vl = ZXC_to_FlxC(V, pp, get_Flx_var(Tl));
840 31939 : Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
841 31941 : return gerepileupto(ltop, FlxX_to_ZXX(Tl));
842 : }
843 8610 : W = cgetg(lg(V),t_VEC);
844 77390 : for(k=1; k < lg(V); k++)
845 68780 : gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
846 8610 : return gerepileupto(ltop, FpXQXV_prod(W, T, p));
847 : }
848 :
849 : GEN
850 144961 : FqV_red(GEN x, GEN T, GEN p)
851 1185172 : { pari_APPLY_same(Fq_red(gel(x,i), T, p)) }
852 :
853 : GEN
854 0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
855 : {
856 0 : if (!T) return FpC_add(x, y, p);
857 0 : pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
858 : }
859 :
860 : GEN
861 0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
862 : {
863 0 : if (!T) return FpC_sub(x, y, p);
864 0 : pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
865 : }
866 :
867 : GEN
868 0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
869 : {
870 0 : if (!T) return FpC_Fp_mul(x, y, p);
871 0 : pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
872 : }
873 :
874 : GEN
875 105 : FqC_FqV_mul(GEN x, GEN y, GEN T, GEN p)
876 : {
877 105 : long i,j, lx=lg(x), ly=lg(y);
878 : GEN z;
879 105 : if (ly==1) return cgetg(1,t_MAT);
880 105 : z = cgetg(ly,t_MAT);
881 819 : for (j=1; j < ly; j++)
882 : {
883 714 : GEN zj = cgetg(lx,t_COL);
884 4200 : for (i=1; i<lx; i++) gel(zj,i) = Fq_mul(gel(x,i),gel(y,j), T, p);
885 714 : gel(z, j) = zj;
886 : }
887 105 : return z;
888 : }
889 :
890 : GEN
891 5271 : FpXC_center(GEN x, GEN p, GEN pov2)
892 40524 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
893 :
894 : GEN
895 1730 : FpXM_center(GEN x, GEN p, GEN pov2)
896 7001 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
897 :
898 : /*******************************************************************/
899 : /* */
900 : /* GENERIC CRT */
901 : /* */
902 : /*******************************************************************/
903 : static GEN
904 5885986 : primelist(forprime_t *S, long n, GEN dB)
905 : {
906 5885986 : GEN P = cgetg(n+1, t_VECSMALL);
907 5885959 : long i = 1;
908 : ulong p;
909 14458804 : while (i <= n && (p = u_forprime_next(S)))
910 8572844 : if (!dB || umodiu(dB, p)) P[i++] = p;
911 5885958 : return P;
912 : }
913 :
914 : void
915 5483229 : gen_inccrt_i(const char *str, GEN worker, GEN dB, long n, long mmin,
916 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
917 : GEN center(GEN, GEN, GEN))
918 : {
919 5483229 : long m = mmin? minss(mmin, n): usqrt(n);
920 : GEN H, P, mod;
921 : pari_timer ti;
922 5483224 : if (DEBUGLEVEL > 4)
923 : {
924 0 : timer_start(&ti);
925 0 : err_printf("%s: nb primes: %ld\n",str, n);
926 : }
927 5483227 : if (m == 1)
928 : {
929 5279547 : GEN P = primelist(S, n, dB);
930 5279502 : GEN done = closure_callgen1(worker, P);
931 5279482 : H = gel(done,1);
932 5279482 : mod = gel(done,2);
933 5279482 : if (!*pH && center) H = center(H, mod, shifti(mod,-1));
934 5279479 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
935 : }
936 : else
937 : {
938 203680 : long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
939 : struct pari_mt pt;
940 203680 : long pending = 0;
941 203680 : H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
942 203680 : mt_queue_start_lim(&pt, worker, m);
943 854655 : for (i=1; i<=m || pending; i++)
944 : {
945 : GEN done;
946 650975 : GEN pr = i <= m ? mkvec(primelist(S, i<=r ? s: s-1, dB)): NULL;
947 650975 : mt_queue_submit(&pt, i, pr);
948 650976 : done = mt_queue_get(&pt, NULL, &pending);
949 650976 : if (done)
950 : {
951 606451 : di++;
952 606451 : gel(H, di) = gel(done,1);
953 606451 : gel(P, di) = gel(done,2);
954 606451 : if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
955 : }
956 : }
957 203680 : mt_queue_end(&pt);
958 203680 : if (DEBUGLEVEL>5) err_printf("\n");
959 203680 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
960 203680 : H = crt(H, P, &mod);
961 203680 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
962 : }
963 5483159 : if (*pH) H = crt(mkvec2(*pH, H), mkvec2(*pmod, mod), &mod);
964 5483160 : *pH = H; *pmod = mod;
965 5483160 : }
966 : void
967 1920995 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
968 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
969 : GEN center(GEN, GEN, GEN))
970 : {
971 1920995 : pari_sp av = avma;
972 1920995 : gen_inccrt_i(str, worker, dB, n, mmin, S, pH, pmod, crt, center);
973 1920945 : gerepileall(av, 2, pH, pmod);
974 1921070 : }
975 :
976 : GEN
977 1330273 : gen_crt(const char *str, GEN worker, forprime_t *S, GEN dB, ulong bound, long mmin, GEN *pmod,
978 : GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
979 : {
980 1330273 : GEN mod = gen_1, H = NULL;
981 : ulong e;
982 :
983 1330273 : bound++;
984 2660626 : while (bound > (e = expi(mod)))
985 : {
986 1330251 : long n = (bound - e) / expu(S->p) + 1;
987 1330285 : gen_inccrt(str, worker, dB, n, mmin, S, &H, &mod, crt, center);
988 : }
989 1330327 : if (pmod) *pmod = mod;
990 1330327 : return H;
991 : }
992 :
993 : /*******************************************************************/
994 : /* */
995 : /* MODULAR GCD */
996 : /* */
997 : /*******************************************************************/
998 : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
999 : static GEN
1000 5092505 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
1001 : {
1002 5092505 : ulong d, amod = umodiu(a, p);
1003 5092587 : pari_sp av = avma;
1004 : GEN ax;
1005 :
1006 5092587 : if (b == amod) return NULL;
1007 2096855 : d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
1008 2097646 : if (d >= 1 + (p>>1))
1009 1023550 : ax = subii(a, mului(p-d, q));
1010 : else
1011 : {
1012 1074096 : ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
1013 1073557 : if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
1014 : }
1015 2096718 : return gerepileuptoint(av, ax);
1016 : }
1017 : GEN
1018 378 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
1019 : GEN
1020 31477 : ZX_init_CRT(GEN Hp, ulong p, long v)
1021 : {
1022 31477 : long i, l = lg(Hp), lim = (long)(p>>1);
1023 31477 : GEN H = cgetg(l, t_POL);
1024 31477 : H[1] = evalsigne(1) | evalvarn(v);
1025 793616 : for (i=2; i<l; i++)
1026 762137 : gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
1027 31479 : return ZX_renormalize(H,l);
1028 : }
1029 :
1030 : GEN
1031 2828 : ZM_init_CRT(GEN Hp, ulong p)
1032 : {
1033 2828 : long i,j, m, l = lg(Hp), lim = (long)(p>>1);
1034 2828 : GEN c, cp, H = cgetg(l, t_MAT);
1035 2828 : if (l==1) return H;
1036 2828 : m = lgcols(Hp);
1037 10129 : for (j=1; j<l; j++)
1038 : {
1039 7301 : cp = gel(Hp,j);
1040 7301 : c = cgetg(m, t_COL);
1041 7301 : gel(H,j) = c;
1042 79695 : for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
1043 : }
1044 2828 : return H;
1045 : }
1046 :
1047 : int
1048 7581 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
1049 : {
1050 7581 : GEN h, q = *ptq, qp = muliu(q,p);
1051 7581 : ulong qinv = Fl_inv(umodiu(q,p), p);
1052 7581 : int stable = 1;
1053 7581 : h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
1054 7581 : if (h) { *H = h; stable = 0; }
1055 7581 : *ptq = qp; return stable;
1056 : }
1057 :
1058 : static int
1059 147025 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
1060 : {
1061 147025 : GEN H = *ptH, h, qp2 = shifti(qp,-1);
1062 147015 : ulong qinv = Fl_inv(umodiu(q,p), p);
1063 147026 : long i, l = lg(H), lp = lg(Hp);
1064 147026 : int stable = 1;
1065 :
1066 147026 : if (l < lp)
1067 : { /* degree increases */
1068 0 : GEN x = cgetg(lp, t_POL);
1069 0 : for (i=1; i<l; i++) x[i] = H[i];
1070 0 : for ( ; i<lp; i++) gel(x,i) = gen_0;
1071 0 : *ptH = H = x;
1072 0 : stable = 0;
1073 147026 : } else if (l > lp)
1074 : { /* degree decreases */
1075 0 : GEN x = cgetg(l, t_VECSMALL);
1076 0 : for (i=1; i<lp; i++) x[i] = Hp[i];
1077 0 : for ( ; i<l; i++) x[i] = 0;
1078 0 : Hp = x; lp = l;
1079 : }
1080 4926896 : for (i=2; i<lp; i++)
1081 : {
1082 4779955 : h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
1083 4779870 : if (h) { gel(H,i) = h; stable = 0; }
1084 : }
1085 146941 : (void)ZX_renormalize(H,lp);
1086 147028 : return stable;
1087 : }
1088 :
1089 : int
1090 0 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
1091 : {
1092 0 : GEN q = *ptq, qp = muliu(q,p);
1093 0 : int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
1094 0 : *ptq = qp; return stable;
1095 : }
1096 :
1097 : int
1098 4771 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1099 : {
1100 4771 : GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1101 4771 : ulong qinv = Fl_inv(umodiu(q,p), p);
1102 4771 : long i,j, l = lg(H), m = lgcols(H);
1103 4771 : int stable = 1;
1104 17877 : for (j=1; j<l; j++)
1105 140327 : for (i=1; i<m; i++)
1106 : {
1107 127221 : h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
1108 127221 : if (h) { gcoeff(H,i,j) = h; stable = 0; }
1109 : }
1110 4771 : *ptq = qp; return stable;
1111 : }
1112 :
1113 : GEN
1114 623 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
1115 : {
1116 : long i, j, k;
1117 : GEN H;
1118 623 : long m, l = lg(Hp), lim = (long)(p>>1), n;
1119 623 : H = cgetg(l, t_MAT);
1120 623 : if (l==1) return H;
1121 623 : m = lgcols(Hp);
1122 623 : n = deg + 3;
1123 2114 : for (j=1; j<l; j++)
1124 : {
1125 1491 : GEN cp = gel(Hp,j);
1126 1491 : GEN c = cgetg(m, t_COL);
1127 1491 : gel(H,j) = c;
1128 23905 : for (i=1; i<m; i++)
1129 : {
1130 22414 : GEN dp = gel(cp, i);
1131 22414 : long l = lg(dp);
1132 22414 : GEN d = cgetg(n, t_POL);
1133 22414 : gel(c, i) = d;
1134 22414 : d[1] = dp[1] | evalsigne(1);
1135 45647 : for (k=2; k<l; k++)
1136 23233 : gel(d,k) = stoi(Fl_center(dp[k], p, lim));
1137 44457 : for ( ; k<n; k++)
1138 22043 : gel(d,k) = gen_0;
1139 : }
1140 : }
1141 623 : return H;
1142 : }
1143 :
1144 : int
1145 653 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1146 : {
1147 653 : GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1148 653 : ulong qinv = Fl_inv(umodiu(q,p), p);
1149 653 : long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
1150 653 : int stable = 1;
1151 2225 : for (j=1; j<l; j++)
1152 90418 : for (i=1; i<m; i++)
1153 : {
1154 88846 : GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
1155 88846 : long lh = lg(hp);
1156 246641 : for (k=2; k<lh; k++)
1157 : {
1158 157795 : v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
1159 157795 : if (v) { gel(h,k) = v; stable = 0; }
1160 : }
1161 108763 : for (; k<n; k++)
1162 : {
1163 19917 : v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
1164 19917 : if (v) { gel(h,k) = v; stable = 0; }
1165 : }
1166 : }
1167 653 : *ptq = qp; return stable;
1168 : }
1169 :
1170 : /* record the degrees of Euclidean remainders (make them as large as
1171 : * possible : smaller values correspond to a degenerate sequence) */
1172 : static void
1173 23065 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
1174 : {
1175 : long da,db,dc, ind;
1176 23065 : pari_sp av = avma;
1177 :
1178 23065 : if (lgpol(a)==0 || lgpol(b)==0) return;
1179 21798 : da = degpol(a);
1180 21798 : db = degpol(b);
1181 21798 : if (db > da)
1182 0 : { swapspec(a,b, da,db); }
1183 21798 : else if (!da) return;
1184 21798 : ind = 0;
1185 143740 : while (db)
1186 : {
1187 121946 : GEN c = Flx_rem(a,b, p);
1188 121942 : a = b; b = c; dc = degpol(c);
1189 121942 : if (dc < 0) break;
1190 :
1191 121942 : ind++;
1192 121942 : if (dc > dglist[ind]) dglist[ind] = dc;
1193 121942 : if (gc_needed(av,2))
1194 : {
1195 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1196 0 : gerepileall(av, 2, &a,&b);
1197 : }
1198 121942 : db = dc; /* = degpol(b) */
1199 : }
1200 21794 : if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
1201 21794 : set_avma(av);
1202 : }
1203 : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
1204 : * "generic" degree sequence as given by dglist, set *Ci and return
1205 : * resultant(a,b). Modular version of Collins's subresultant */
1206 : static ulong
1207 2080981 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
1208 : {
1209 : long da,db,dc, ind;
1210 2080981 : ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
1211 2080981 : int s = 1;
1212 2080981 : pari_sp av = avma;
1213 :
1214 2080981 : *C0 = 1; *C1 = 0;
1215 2080981 : if (lgpol(a)==0 || lgpol(b)==0) return 0;
1216 2071685 : da = degpol(a);
1217 2071676 : db = degpol(b);
1218 2071744 : if (db > da)
1219 : {
1220 0 : swapspec(a,b, da,db);
1221 0 : if (both_odd(da,db)) s = -s;
1222 : }
1223 2071744 : else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
1224 2071744 : ind = 0;
1225 19765566 : while (db)
1226 : { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
1227 : * da = deg a, db = deg b */
1228 17698654 : GEN c = Flx_rem(a,b, p);
1229 17557539 : long delta = da - db;
1230 :
1231 17557539 : if (both_odd(da,db)) s = -s;
1232 17565934 : lb = Fl_mul(b[db+2], cb, p);
1233 17603059 : a = b; b = c; dc = degpol(c);
1234 17611854 : ind++;
1235 17611854 : if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
1236 17606802 : if (g == h)
1237 : { /* frequent */
1238 17546234 : ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
1239 17634130 : ca = cb;
1240 17634130 : cb = cc;
1241 : }
1242 : else
1243 : {
1244 60568 : ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
1245 60568 : ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
1246 60568 : ca = cb;
1247 60568 : cb = Fl_div(cc, ghdelta, p);
1248 : }
1249 17695594 : da = db; /* = degpol(a) */
1250 17695594 : db = dc; /* = degpol(b) */
1251 :
1252 17695594 : g = lb;
1253 17695594 : if (delta == 1)
1254 17595874 : h = g; /* frequent */
1255 : else
1256 99720 : h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
1257 :
1258 17692613 : if (gc_needed(av,2))
1259 : {
1260 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1261 0 : gerepileall(av, 2, &a,&b);
1262 : }
1263 : }
1264 2066912 : if (da > 1) return 0; /* Failure */
1265 : /* last nonconstant polynomial has degree 1 */
1266 2066912 : *C0 = Fl_mul(ca, a[2], p);
1267 2066872 : *C1 = Fl_mul(ca, a[3], p);
1268 2066883 : res = Fl_mul(cb, b[2], p);
1269 2066865 : if (s == -1) res = p - res;
1270 2066865 : return gc_ulong(av,res);
1271 : }
1272 :
1273 : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
1274 : * Return 0 in case of degree drop. */
1275 : static GEN
1276 2104279 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
1277 : {
1278 : GEN z;
1279 2104279 : long i, lb = lg(Q);
1280 2104279 : ulong leadz = Flx_eval(leading_coeff(Q), x, p);
1281 2103663 : long vs=mael(Q,2,1);
1282 2103663 : if (!leadz) return zero_Flx(vs);
1283 :
1284 2093003 : z = cgetg(lb, t_VECSMALL); z[1] = vs;
1285 20009181 : for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
1286 2086887 : z[i] = leadz; return z;
1287 : }
1288 :
1289 : GEN
1290 1232 : FpXY_FpXQ_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
1291 : {
1292 1232 : pari_sp av = avma;
1293 1232 : long i, lb = lg(Q);
1294 : GEN z;
1295 1232 : if (lb == 2) return pol_0(vx);
1296 1232 : z = gel(Q, lb-1);
1297 1232 : if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1298 :
1299 1232 : if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1300 26572 : for (i=lb-2; i>=2; i--)
1301 : {
1302 25340 : GEN c = gel(Q,i);
1303 25340 : z = FqX_Fq_mul(z, y, T, p);
1304 25340 : z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
1305 : }
1306 1232 : return gerepileupto(av, z);
1307 : }
1308 :
1309 : static GEN
1310 152607 : ZX_norml1(GEN x)
1311 : {
1312 152607 : long i, l = lg(x);
1313 : GEN s;
1314 :
1315 152607 : if (l == 2) return gen_0;
1316 133434 : s = gel(x, l-1); /* != 0 */
1317 480801 : for (i = l-2; i > 1; i--) {
1318 347375 : GEN xi = gel(x,i);
1319 347375 : if (!signe(xi)) continue;
1320 238306 : s = addii_sign(s,1, xi,1);
1321 : }
1322 133426 : return s;
1323 : }
1324 : /* x >= 0, y != 0, return x + |y| */
1325 : static GEN
1326 25577 : addii_abs(GEN x, GEN y)
1327 : {
1328 25577 : if (!signe(x)) return absi_shallow(y);
1329 16043 : return addii_sign(x,1, y,1);
1330 : }
1331 :
1332 : /* x a ZX, return sum_{i >= k} |x[i]| binomial(i, k) */
1333 : static GEN
1334 31712 : ZX_norml1_1(GEN x, long k)
1335 : {
1336 31712 : long i, d = degpol(x);
1337 : GEN s, C; /* = binomial(i, k) */
1338 :
1339 31712 : if (!d || k > d) return gen_0;
1340 31712 : s = absi_shallow(gel(x, k+2)); /* may be 0 */
1341 31711 : C = gen_1;
1342 68171 : for (i = k+1; i <= d; i++) {
1343 36464 : GEN xi = gel(x,i+2);
1344 36464 : if (k) C = diviuexact(muliu(C, i), i-k);
1345 36460 : if (signe(xi)) s = addii_abs(s, mulii(C, xi));
1346 : }
1347 31707 : return s;
1348 : }
1349 : /* x has non-negative real coefficients */
1350 : static GEN
1351 2667 : RgX_norml1_1(GEN x, long k)
1352 : {
1353 2667 : long i, d = degpol(x);
1354 : GEN s, C; /* = binomial(i, k) */
1355 :
1356 2667 : if (!d || k > d) return gen_0;
1357 2667 : s = gel(x, k+2); /* may be 0 */
1358 2667 : C = gen_1;
1359 7784 : for (i = k+1; i <= d; i++) {
1360 5117 : GEN xi = gel(x,i+2);
1361 5117 : if (k) C = diviuexact(muliu(C, i), i-k);
1362 5117 : if (!gequal0(xi)) s = gadd(s, gmul(C, xi));
1363 : }
1364 2667 : return s;
1365 : }
1366 :
1367 : /* N_2(A)^2 */
1368 : static GEN
1369 6485 : sqrN2(GEN A, long prec)
1370 : {
1371 6485 : pari_sp av = avma;
1372 6485 : long i, l = lg(A);
1373 6485 : GEN a = gen_0;
1374 32002 : for (i = 2; i < l; i++)
1375 : {
1376 25517 : a = gadd(a, gabs(gnorm(gel(A,i)), prec));
1377 25517 : if (gc_needed(av,1))
1378 : {
1379 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1380 0 : a = gerepileupto(av, a);
1381 : }
1382 : }
1383 6485 : return a;
1384 : }
1385 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1386 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1387 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1388 : * Return e such that Res(A, B) < 2^e */
1389 : static GEN
1390 5806 : RgX_RgXY_ResBound(GEN A, GEN B, long prec)
1391 : {
1392 5806 : pari_sp av = avma;
1393 5806 : GEN b = gen_0, bnd;
1394 5806 : long i, lB = lg(B);
1395 23102 : for (i=2; i<lB; i++)
1396 : {
1397 17296 : GEN t = gel(B,i);
1398 17296 : if (typ(t) == t_POL) t = gnorml1(t, prec);
1399 17296 : b = gadd(b, gabs(gsqr(t), prec));
1400 17296 : if (gc_needed(av,1))
1401 : {
1402 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1403 0 : b = gerepileupto(av, b);
1404 : }
1405 : }
1406 5806 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1407 : gpowgs(b, degpol(A))), prec);
1408 5806 : return gerepileupto(av, bnd);
1409 : }
1410 : /* A,B in C[X] return RgX_RgXY_ResBound(A, B(x+y)) */
1411 : static GEN
1412 679 : RgX_RgXY_ResBound_1(GEN A, GEN B, long prec)
1413 : {
1414 679 : pari_sp av = avma, av2;
1415 679 : GEN b = gen_0, bnd;
1416 679 : long i, lB = lg(B);
1417 679 : B = shallowcopy(B);
1418 3346 : for (i=2; i<lB; i++) gel(B,i) = gabs(gel(B,i), prec);
1419 679 : av2 = avma;
1420 3346 : for (i=2; i<lB; i++)
1421 : {
1422 2667 : b = gadd(b, gsqr(RgX_norml1_1(B, i-2)));
1423 2667 : if (gc_needed(av2,1))
1424 : {
1425 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1426 0 : b = gerepileupto(av2, b);
1427 : }
1428 : }
1429 679 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1430 : gpowgs(b, degpol(A))), prec);
1431 679 : return gerepileupto(av, bnd);
1432 : }
1433 :
1434 : /* log2 N_2(A)^2 */
1435 : static double
1436 175046 : log2N2(GEN A)
1437 : {
1438 175046 : pari_sp av = avma;
1439 175046 : long i, l = lg(A);
1440 175046 : GEN a = gen_0;
1441 1136256 : for (i=2; i < l; i++)
1442 : {
1443 961221 : a = addii(a, sqri(gel(A,i)));
1444 961211 : if (gc_needed(av,1))
1445 : {
1446 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1447 0 : a = gerepileupto(av, a);
1448 : }
1449 : }
1450 175035 : return gc_double(av, dbllog2(a));
1451 : }
1452 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1453 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1454 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1455 : * Return e such that Res(A, B) < 2^e */
1456 : ulong
1457 164941 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
1458 : {
1459 164941 : pari_sp av = avma;
1460 164941 : GEN b = gen_0;
1461 164941 : long i, lB = lg(B);
1462 : double logb;
1463 929264 : for (i=2; i<lB; i++)
1464 : {
1465 764335 : GEN t = gel(B,i);
1466 764335 : if (typ(t) == t_POL) t = ZX_norml1(t);
1467 764337 : b = addii(b, sqri(t));
1468 764324 : if (gc_needed(av,1))
1469 : {
1470 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1471 0 : b = gerepileupto(av, b);
1472 : }
1473 : }
1474 164929 : logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
1475 164941 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * logb) / 2);
1476 164940 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1477 : }
1478 : /* A,B ZX. Return ZX_ZXY_ResBound(A(x), B(x+y)) */
1479 : static ulong
1480 10103 : ZX_ZXY_ResBound_1(GEN A, GEN B)
1481 : {
1482 10103 : pari_sp av = avma;
1483 10103 : GEN b = gen_0;
1484 10103 : long i, lB = lg(B);
1485 41817 : for (i=2; i<lB; i++)
1486 : {
1487 31712 : b = addii(b, sqri(ZX_norml1_1(B, i-2)));
1488 31714 : if (gc_needed(av,1))
1489 : {
1490 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1491 0 : b = gerepileupto(av, b);
1492 : }
1493 : }
1494 10105 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * dbllog2(b)) / 2);
1495 10108 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1496 : }
1497 : /* special case B = A' */
1498 : static ulong
1499 1080157 : ZX_discbound(GEN A)
1500 : {
1501 1080157 : pari_sp av = avma;
1502 1080157 : GEN a = gen_0, b = gen_0;
1503 1080157 : long i , lA = lg(A), dA = degpol(A);
1504 : double loga, logb;
1505 6500451 : for (i = 2; i < lA; i++)
1506 : {
1507 5420682 : GEN c = sqri(gel(A,i));
1508 5420065 : a = addii(a, c);
1509 5420199 : if (i > 2) b = addii(b, mulii(c, sqru(i-2)));
1510 5420273 : if (gc_needed(av,1))
1511 : {
1512 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_discbound i = %ld",i);
1513 0 : gerepileall(av, 2, &a, &b);
1514 : }
1515 : }
1516 1079769 : loga = dbllog2(a);
1517 1080073 : logb = dbllog2(b); set_avma(av);
1518 1080116 : i = (long)(((dA-1) * loga + dA * logb) / 2);
1519 1080116 : return (i <= 0)? 1: 1 + (ulong)i;
1520 : }
1521 :
1522 : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
1523 : static ulong
1524 1228704 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong la)
1525 : {
1526 1228704 : GEN ev = FlxY_evalx(b, n, p);
1527 1228738 : long drop = lg(b) - lg(ev);
1528 1228738 : ulong r = Flx_resultant(a, ev, p);
1529 1228695 : if (drop && la != 1) r = Fl_mul(r, Fl_powu(la, drop,p),p);
1530 1228695 : return r;
1531 : }
1532 : static GEN
1533 4 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
1534 : {
1535 4 : GEN ev = FpXY_evaly(b, n, p, vX);
1536 4 : long drop = db-degpol(ev);
1537 4 : GEN r = FpX_resultant(a, ev, p);
1538 4 : if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
1539 4 : return r;
1540 : }
1541 :
1542 : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
1543 : /* Return a Fly */
1544 : static GEN
1545 110856 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, long dres, long sx)
1546 : {
1547 : long i;
1548 110856 : ulong n, la = Flx_lead(a);
1549 110856 : GEN x = cgetg(dres+2, t_VECSMALL);
1550 110856 : GEN y = cgetg(dres+2, t_VECSMALL);
1551 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1552 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1553 672483 : for (i=0,n = 1; i < dres; n++)
1554 : {
1555 561624 : x[++i] = n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
1556 561613 : x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
1557 : }
1558 110859 : if (i == dres)
1559 : {
1560 105555 : x[++i] = 0; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,la);
1561 : }
1562 110858 : return Flv_polint(x,y, p, sx);
1563 : }
1564 :
1565 : static GEN
1566 7906 : FlxX_pseudorem(GEN x, GEN y, ulong p)
1567 : {
1568 7906 : long vx = varn(x), dx, dy, dz, i, lx, dp;
1569 7906 : pari_sp av = avma, av2;
1570 :
1571 7906 : if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
1572 7906 : (void)new_chunk(2);
1573 7907 : dx=degpol(x); x = RgX_recip_i(x)+2;
1574 7909 : dy=degpol(y); y = RgX_recip_i(y)+2; dz=dx-dy; dp = dz+1;
1575 7909 : av2 = avma;
1576 : for (;;)
1577 : {
1578 65079 : gel(x,0) = Flx_neg(gel(x,0), p); dp--;
1579 243212 : for (i=1; i<=dy; i++)
1580 176469 : gel(x,i) = Flx_add( Flx_mul(gel(y,0), gel(x,i), p),
1581 178064 : Flx_mul(gel(x,0), gel(y,i), p), p );
1582 1135732 : for ( ; i<=dx; i++)
1583 1071704 : gel(x,i) = Flx_mul(gel(y,0), gel(x,i), p);
1584 69119 : do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
1585 64028 : if (dx < dy) break;
1586 56124 : if (gc_needed(av2,1))
1587 : {
1588 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
1589 0 : gerepilecoeffs(av2,x,dx+1);
1590 : }
1591 : }
1592 7904 : if (dx < 0) return zero_Flx(0);
1593 7904 : lx = dx+3; x -= 2;
1594 7904 : x[0]=evaltyp(t_POL) | evallg(lx);
1595 7902 : x[1]=evalsigne(1) | evalvarn(vx);
1596 7902 : x = RgX_recip_i(x);
1597 7903 : if (dp)
1598 : { /* multiply by y[0]^dp [beware dummy vars from FpX_FpXY_resultant] */
1599 2089 : GEN t = Flx_powu(gel(y,0), dp, p);
1600 8368 : for (i=2; i<lx; i++)
1601 6276 : gel(x,i) = Flx_mul(gel(x,i), t, p);
1602 : }
1603 7906 : return gerepilecopy(av, x);
1604 : }
1605 :
1606 : /* return a Flx */
1607 : GEN
1608 2647 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
1609 : {
1610 2647 : pari_sp av = avma, av2;
1611 : long degq,dx,dy,du,dv,dr,signh;
1612 : GEN z,g,h,r,p1;
1613 :
1614 2647 : dx=degpol(u); dy=degpol(v); signh=1;
1615 2647 : if (dx < dy)
1616 : {
1617 7 : swap(u,v); lswap(dx,dy);
1618 7 : if (both_odd(dx, dy)) signh = -signh;
1619 : }
1620 2647 : if (dy < 0) return zero_Flx(sx);
1621 2647 : if (dy==0) return gerepileupto(av, Flx_powu(gel(v,2),dx,p));
1622 :
1623 2647 : g = h = pol1_Flx(sx); av2 = avma;
1624 : for(;;)
1625 : {
1626 7910 : r = FlxX_pseudorem(u,v,p); dr = lg(r);
1627 7913 : if (dr == 2) { set_avma(av); return zero_Flx(sx); }
1628 7913 : du = degpol(u); dv = degpol(v); degq = du-dv;
1629 7911 : u = v; p1 = g; g = leading_coeff(u);
1630 7911 : switch(degq)
1631 : {
1632 0 : case 0: break;
1633 5806 : case 1:
1634 5806 : p1 = Flx_mul(h,p1, p); h = g; break;
1635 2105 : default:
1636 2105 : p1 = Flx_mul(Flx_powu(h,degq,p), p1, p);
1637 2105 : h = Flx_div(Flx_powu(g,degq,p), Flx_powu(h,degq-1,p), p);
1638 : }
1639 7908 : if (both_odd(du,dv)) signh = -signh;
1640 7903 : v = FlxY_Flx_div(r, p1, p);
1641 7903 : if (dr==3) break;
1642 5256 : if (gc_needed(av2,1))
1643 : {
1644 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
1645 0 : gerepileall(av2,4, &u, &v, &g, &h);
1646 : }
1647 : }
1648 2647 : z = gel(v,2);
1649 2647 : if (dv > 1) z = Flx_div(Flx_powu(z,dv,p), Flx_powu(h,dv-1,p), p);
1650 2647 : if (signh < 0) z = Flx_neg(z,p);
1651 2647 : return gerepileupto(av, z);
1652 : }
1653 :
1654 : /* Warning:
1655 : * This function switches between valid and invalid variable ordering*/
1656 :
1657 : static GEN
1658 6259 : FlxY_to_FlyX(GEN b, long sv)
1659 : {
1660 6259 : long i, n=-1;
1661 6259 : long sw = b[1]&VARNBITS;
1662 21431 : for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
1663 6259 : return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
1664 : }
1665 :
1666 : /* Return a Fly*/
1667 : GEN
1668 6259 : Flx_FlxY_resultant(GEN a, GEN b, ulong pp)
1669 : {
1670 6259 : pari_sp ltop=avma;
1671 6259 : long dres = degpol(a)*degpol(b);
1672 6259 : long sx=a[1], sy=b[1]&VARNBITS;
1673 : GEN z;
1674 6259 : b = FlxY_to_FlyX(b,sx);
1675 6256 : if ((ulong)dres >= pp)
1676 2644 : z = FlxX_resultant(Fly_to_FlxY(a, sy), b, pp, sx);
1677 : else
1678 3612 : z = Flx_FlxY_resultant_polint(a, b, pp, (ulong)dres, sy);
1679 6259 : return gerepileupto(ltop,z);
1680 : }
1681 :
1682 : /* return a t_POL (in variable v >= 0) whose coeffs are the coeffs of b,
1683 : * in variable v. This is an incorrect PARI object if initially varn(b) << v.
1684 : * We could return a vector of coeffs, but it is convenient to have degpol()
1685 : * and friends available. Even in that case, it will behave nicely with all
1686 : * functions treating a polynomial as a vector of coeffs (eg poleval).
1687 : * FOR INTERNAL USE! */
1688 : GEN
1689 110117 : swap_vars(GEN b0, long v)
1690 : {
1691 110117 : long i, n = RgX_degree(b0, v);
1692 : GEN b, x;
1693 110117 : if (n < 0) return pol_0(v);
1694 110117 : b = cgetg(n+3, t_POL); x = b + 2;
1695 110116 : b[1] = evalsigne(1) | evalvarn(v);
1696 502719 : for (i=0; i<=n; i++) gel(x,i) = polcoef_i(b0, i, v);
1697 110116 : return b;
1698 : }
1699 :
1700 : /* assume varn(b) << varn(a) */
1701 : /* return a FpY*/
1702 : GEN
1703 1 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
1704 : {
1705 1 : long i,n,dres, db, vY = varn(b), vX = varn(a);
1706 : GEN la,x,y;
1707 :
1708 1 : if (lgefint(p) == 3)
1709 : {
1710 0 : ulong pp = uel(p,2);
1711 0 : b = ZXX_to_FlxX(b, pp, vX);
1712 0 : a = ZX_to_Flx(a, pp);
1713 0 : x = Flx_FlxY_resultant(a, b, pp);
1714 0 : return Flx_to_ZX(x);
1715 : }
1716 1 : db = RgXY_degreex(b);
1717 1 : dres = degpol(a)*degpol(b);
1718 1 : la = leading_coeff(a);
1719 1 : x = cgetg(dres+2, t_VEC);
1720 1 : y = cgetg(dres+2, t_VEC);
1721 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1722 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1723 3 : for (i=0,n = 1; i < dres; n++)
1724 : {
1725 2 : gel(x,++i) = utoipos(n);
1726 2 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1727 2 : gel(x,++i) = subiu(p,n);
1728 2 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1729 : }
1730 1 : if (i == dres)
1731 : {
1732 0 : gel(x,++i) = gen_0;
1733 0 : gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
1734 : }
1735 1 : return FpV_polint(x,y, p, vY);
1736 : }
1737 :
1738 : static GEN
1739 30 : FpX_composedsum(GEN P, GEN Q, GEN p)
1740 : {
1741 30 : long n = 1+ degpol(P)*degpol(Q);
1742 30 : GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
1743 30 : GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
1744 30 : GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
1745 30 : GEN lead = Fp_mul(Fp_powu(leading_coeff(P),degpol(Q), p),
1746 30 : Fp_powu(leading_coeff(Q),degpol(P), p), p);
1747 30 : GEN R = FpX_fromNewton(L, p);
1748 30 : return FpX_Fp_mul(R, lead, p);
1749 : }
1750 :
1751 : #if 0
1752 : GEN
1753 : FpX_composedprod(GEN P, GEN Q, GEN p)
1754 : {
1755 : long n = 1+ degpol(P)*degpol(Q);
1756 : GEN L=FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
1757 : return FpX_fromNewton(L, p);
1758 : }
1759 : #endif
1760 :
1761 : GEN
1762 30 : FpX_direct_compositum(GEN a, GEN b, GEN p)
1763 : {
1764 30 : if (lgefint(p)==3)
1765 : {
1766 0 : pari_sp av = avma;
1767 0 : ulong pp = p[2];
1768 0 : GEN z = Flx_direct_compositum(ZX_to_Flx(a, pp), ZX_to_Flx(b, pp), pp);
1769 0 : return gerepileupto(av, Flx_to_ZX(z));
1770 : }
1771 30 : return FpX_composedsum(a, b, p);
1772 : }
1773 :
1774 : static GEN
1775 30 : _FpX_direct_compositum(void *E, GEN a, GEN b)
1776 30 : { return FpX_direct_compositum(a,b, (GEN)E); }
1777 :
1778 : GEN
1779 497 : FpXV_direct_compositum(GEN V, GEN p)
1780 : {
1781 497 : if (lgefint(p)==3)
1782 : {
1783 0 : ulong pp = p[2];
1784 0 : return Flx_to_ZX(FlxV_direct_compositum(ZXV_to_FlxV(V, pp), pp));
1785 : }
1786 497 : return gen_product(V, (void *)p, &_FpX_direct_compositum);
1787 : }
1788 :
1789 : /* 0, 1, -1, 2, -2, ... */
1790 : #define next_lambda(a) (a>0 ? -a : 1-a)
1791 :
1792 : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
1793 : * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
1794 : * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
1795 : * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
1796 : * the Last nonconstant polynomial in the Euclidean Remainder Sequence */
1797 : static GEN
1798 20733 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
1799 : {
1800 : ulong bound, dp;
1801 20733 : pari_sp av = avma, av2 = 0;
1802 20733 : long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
1803 : long stable, checksqfree, i,n, cnt, degB;
1804 20733 : long v, vX = varn(B0), vY = varn(A); /* vY < vX */
1805 : GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
1806 : forprime_t S;
1807 :
1808 20733 : if (degA == 1)
1809 : {
1810 966 : GEN a1 = gel(A,3), a0 = gel(A,2);
1811 966 : B = lambda? RgX_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
1812 966 : H = gsubst(B, vY, gdiv(gneg(a0),a1));
1813 966 : if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
1814 966 : *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
1815 966 : return gc_all(av, 2, &H, LERS);
1816 : }
1817 :
1818 19767 : dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
1819 19767 : C0 = cgetg(dres+2, t_VECSMALL);
1820 19767 : C1 = cgetg(dres+2, t_VECSMALL);
1821 19768 : dglist = cgetg(dres+1, t_VECSMALL);
1822 19768 : x = cgetg(dres+2, t_VECSMALL);
1823 19768 : y = cgetg(dres+2, t_VECSMALL);
1824 19768 : B0 = leafcopy(B0);
1825 19768 : A = leafcopy(A);
1826 19767 : B = B0;
1827 19767 : v = fetch_var_higher(); setvarn(A,v);
1828 : /* make sure p large enough */
1829 20369 : INIT:
1830 : /* always except the first time */
1831 20369 : if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
1832 20369 : if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
1833 20369 : B = swap_vars(B, vY); setvarn(B,v);
1834 : /* B0(lambda v + x, v) */
1835 20369 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
1836 20369 : av2 = avma;
1837 :
1838 20369 : if (degA <= 3)
1839 : { /* sub-resultant faster for small degrees */
1840 9779 : H = RgX_resultant_all(A,B,&q);
1841 9779 : if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
1842 9317 : H0 = gel(q,2);
1843 9317 : if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
1844 9317 : H1 = gel(q,3);
1845 9317 : if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
1846 9317 : if (!ZX_is_squarefree(H)) goto INIT;
1847 9275 : goto END;
1848 : }
1849 :
1850 10590 : H = H0 = H1 = NULL;
1851 10590 : degB = degpol(B);
1852 10589 : bound = ZX_ZXY_ResBound(A, B, NULL);
1853 10591 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
1854 10591 : dp = 1;
1855 10591 : init_modular_big(&S);
1856 10591 : for(cnt = 0, checksqfree = 1;;)
1857 49010 : {
1858 59601 : ulong p = u_forprime_next(&S);
1859 : GEN Hi;
1860 59601 : a = ZX_to_Flx(A, p);
1861 59598 : b = ZXX_to_FlxX(B, p, varn(A));
1862 59601 : if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
1863 59601 : if (checksqfree)
1864 : { /* find degree list for generic Euclidean Remainder Sequence */
1865 10591 : long goal = minss(degpol(a), degpol(b)); /* longest possible */
1866 72779 : for (n=1; n <= goal; n++) dglist[n] = 0;
1867 10591 : setlg(dglist, 1);
1868 23449 : for (n=0; n <= dres; n++)
1869 : {
1870 23064 : ev = FlxY_evalx_drop(b, n, p);
1871 23065 : Flx_resultant_set_dglist(a, ev, dglist, p);
1872 23064 : if (lg(dglist)-1 == goal) break;
1873 : }
1874 : /* last pol in ERS has degree > 1 ? */
1875 10591 : goal = lg(dglist)-1;
1876 10591 : if (degpol(B) == 1) { if (!goal) goto INIT; }
1877 : else
1878 : {
1879 10535 : if (goal <= 1) goto INIT;
1880 10493 : if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
1881 : }
1882 10549 : if (DEBUGLEVEL>4)
1883 0 : err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
1884 : }
1885 :
1886 2140772 : for (i=0,n = 0; i <= dres; n++)
1887 : {
1888 2081257 : ev = FlxY_evalx_drop(b, n, p);
1889 2080894 : x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
1890 2081213 : if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
1891 : }
1892 59515 : Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
1893 59559 : Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
1894 59559 : if (!H && degpol(Hp) != dres) continue;
1895 59559 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
1896 59559 : if (checksqfree) {
1897 10549 : if (!Flx_is_squarefree(Hp, p)) goto INIT;
1898 10493 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
1899 10493 : checksqfree = 0;
1900 : }
1901 :
1902 59503 : if (!H)
1903 : { /* initialize */
1904 10493 : q = utoipos(p); stable = 0;
1905 10492 : H = ZX_init_CRT(Hp, p,vX);
1906 10492 : H0= ZX_init_CRT(H0p, p,vX);
1907 10493 : H1= ZX_init_CRT(H1p, p,vX);
1908 : }
1909 : else
1910 : {
1911 49010 : GEN qp = muliu(q,p);
1912 49008 : stable = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
1913 49010 : & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
1914 49010 : & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
1915 49010 : q = qp;
1916 : }
1917 : /* could make it probabilistic for H ? [e.g if stable twice, etc]
1918 : * Probabilistic anyway for H0, H1 */
1919 59503 : if (DEBUGLEVEL>5 && (stable || ++cnt==100))
1920 0 : { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
1921 59503 : if (stable && (ulong)expi(q) >= bound) break; /* DONE */
1922 49010 : if (gc_needed(av,2))
1923 : {
1924 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
1925 0 : gerepileall(av2, 4, &H, &q, &H0, &H1);
1926 : }
1927 : }
1928 19768 : END:
1929 19768 : if (DEBUGLEVEL>5) err_printf(" done\n");
1930 19768 : setvarn(H, vX); (void)delete_var();
1931 19768 : *LERS = mkvec2(H0,H1);
1932 19768 : *plambda = lambda; return gc_all(av, 2, &H, LERS);
1933 : }
1934 :
1935 : GEN
1936 58365 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
1937 : {
1938 58365 : if (LERS)
1939 : {
1940 20733 : if (!plambda)
1941 0 : pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
1942 20733 : return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
1943 : }
1944 37632 : return ZX_ZXY_rnfequation(A, B, plambda);
1945 : }
1946 :
1947 : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
1948 : * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
1949 : * squarefree */
1950 : GEN
1951 3400 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
1952 : {
1953 3400 : pari_sp av = avma;
1954 : GEN R, a;
1955 : long dA;
1956 : int delvar;
1957 :
1958 3400 : if (v < 0) v = 0;
1959 3400 : switch (typ(A))
1960 : {
1961 3400 : case t_POL: dA = degpol(A); if (dA > 0) break;
1962 0 : A = constant_coeff(A);
1963 0 : default:
1964 0 : if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
1965 0 : return gerepileupto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
1966 : }
1967 3400 : delvar = 0;
1968 3400 : if (varn(T) == 0)
1969 : {
1970 3192 : long v0 = fetch_var(); delvar = 1;
1971 3192 : T = leafcopy(T); setvarn(T,v0);
1972 3192 : A = leafcopy(A); setvarn(A,v0);
1973 : }
1974 3400 : R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
1975 3400 : if (delvar) (void)delete_var();
1976 3400 : setvarn(R, v); a = leading_coeff(T);
1977 3400 : if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
1978 3400 : return gerepileupto(av, R);
1979 : }
1980 :
1981 : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
1982 : GEN
1983 107430 : ZXQ_charpoly(GEN A, GEN T, long v)
1984 : {
1985 107430 : return (degpol(T) < 16) ? RgXQ_charpoly(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
1986 : }
1987 :
1988 : GEN
1989 7616 : QXQ_charpoly(GEN A, GEN T, long v)
1990 : {
1991 7616 : pari_sp av = avma;
1992 7616 : GEN den, B = Q_remove_denom(A, &den);
1993 7615 : GEN P = ZXQ_charpoly(B, T, v);
1994 7616 : return gerepilecopy(av, den ? RgX_rescale(P, ginv(den)): P);
1995 : }
1996 :
1997 : static ulong
1998 3848081 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
1999 : {
2000 3848081 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2001 : ulong H, dp;
2002 3847930 : if (dropa && dropb) return 0; /* p | lc(A), p | lc(B) */
2003 3847930 : H = Flx_resultant(a, b, p);
2004 3847800 : if (dropa)
2005 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2006 0 : ulong c = b[degB+2]; /* lc(B) */
2007 0 : if (odd(degB)) c = p - c;
2008 0 : c = Fl_powu(c, dropa, p);
2009 0 : if (c != 1) H = Fl_mul(H, c, p);
2010 : }
2011 3847800 : else if (dropb)
2012 : { /* multiply by lc(A)^(deg B - deg b) */
2013 0 : ulong c = a[degA+2]; /* lc(A) */
2014 0 : c = Fl_powu(c, dropb, p);
2015 0 : if (c != 1) H = Fl_mul(H, c, p);
2016 : }
2017 3847812 : dp = dB ? umodiu(dB, p): 1;
2018 3847812 : if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
2019 3847812 : return H;
2020 : }
2021 :
2022 : /* If B=NULL, assume B=A' */
2023 : static GEN
2024 1579315 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
2025 : {
2026 1579315 : pari_sp av = avma, av2;
2027 1579315 : long degA, degB, i, n = lg(P)-1;
2028 : GEN H, T;
2029 :
2030 1579315 : degA = degpol(A);
2031 1579319 : degB = B? degpol(B): degA - 1;
2032 1579324 : if (n == 1)
2033 : {
2034 898091 : ulong Hp, p = uel(P,1);
2035 898091 : GEN a = ZX_to_Flx(A, p), b = B? ZX_to_Flx(B, p): Flx_deriv(a, p);
2036 898048 : Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
2037 898110 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
2038 : }
2039 681233 : T = ZV_producttree(P);
2040 681224 : A = ZX_nv_mod_tree(A, P, T);
2041 681220 : if (B) B = ZX_nv_mod_tree(B, P, T);
2042 681222 : H = cgetg(n+1, t_VECSMALL); av2 = avma;
2043 3631035 : for(i=1; i <= n; i++, set_avma(av2))
2044 : {
2045 2949828 : ulong p = P[i];
2046 2949828 : GEN a = gel(A,i), b = B? gel(B,i): Flx_deriv(a, p);
2047 2950054 : H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
2048 : }
2049 681207 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
2050 681219 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2051 : }
2052 :
2053 : GEN
2054 1579328 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
2055 : {
2056 1579328 : GEN V = cgetg(3, t_VEC);
2057 1579319 : if (typ(B) == t_INT) B = NULL;
2058 1579319 : if (!signe(dB)) dB = NULL;
2059 1579319 : gel(V,1) = ZX_resultant_slice(A, B, dB, P, &gel(V,2));
2060 1579309 : return V;
2061 : }
2062 :
2063 : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
2064 : * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
2065 : GEN
2066 1253408 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
2067 : {
2068 1253408 : pari_sp av = avma;
2069 : forprime_t S;
2070 : GEN H, worker;
2071 1253408 : if (B)
2072 : {
2073 111438 : long a = degpol(A), b = degpol(B);
2074 111438 : if (a < 0 || b < 0) return gen_0;
2075 111407 : if (!a) return powiu(gel(A,2), b);
2076 111407 : if (!b) return powiu(gel(B,2), a);
2077 109675 : if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
2078 : }
2079 1251638 : worker = snm_closure(is_entry("_ZX_resultant_worker"),
2080 : mkvec3(A, B? B: gen_0, dB? dB: gen_0));
2081 1251688 : init_modular_big(&S);
2082 1251656 : H = gen_crt("ZX_resultant_all", worker, &S, dB, bound, 0, NULL,
2083 : ZV_chinese_center, Fp_center);
2084 1251694 : return gerepileuptoint(av, H);
2085 : }
2086 :
2087 : /* A0 and B0 in Q[X] */
2088 : GEN
2089 56 : QX_resultant(GEN A0, GEN B0)
2090 : {
2091 : GEN s, a, b, A, B;
2092 56 : pari_sp av = avma;
2093 :
2094 56 : A = Q_primitive_part(A0, &a);
2095 56 : B = Q_primitive_part(B0, &b);
2096 56 : s = ZX_resultant(A, B);
2097 56 : if (!signe(s)) { set_avma(av); return gen_0; }
2098 56 : if (a) s = gmul(s, gpowgs(a,degpol(B)));
2099 56 : if (b) s = gmul(s, gpowgs(b,degpol(A)));
2100 56 : return gerepileupto(av, s);
2101 : }
2102 :
2103 : GEN
2104 21637 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
2105 :
2106 : GEN
2107 0 : QXQ_intnorm(GEN A, GEN B)
2108 : {
2109 : GEN c, n, R, lB;
2110 0 : long dA = degpol(A), dB = degpol(B);
2111 0 : pari_sp av = avma;
2112 0 : if (dA < 0) return gen_0;
2113 0 : A = Q_primitive_part(A, &c);
2114 0 : if (!c || typ(c) == t_INT) {
2115 0 : n = c;
2116 0 : R = ZX_resultant(B, A);
2117 : } else {
2118 0 : n = gel(c,1);
2119 0 : R = ZX_resultant_all(B, A, gel(c,2), 0);
2120 : }
2121 0 : if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
2122 0 : lB = leading_coeff(B);
2123 0 : if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
2124 0 : return gerepileuptoint(av, R);
2125 : }
2126 :
2127 : GEN
2128 18249 : QXQ_norm(GEN A, GEN B)
2129 : {
2130 : GEN c, R, lB;
2131 18249 : long dA = degpol(A), dB = degpol(B);
2132 18249 : pari_sp av = avma;
2133 18249 : if (dA < 0) return gen_0;
2134 18249 : A = Q_primitive_part(A, &c);
2135 18249 : R = ZX_resultant(B, A);
2136 18249 : if (c) R = gmul(R, gpowgs(c, dB));
2137 18249 : lB = leading_coeff(B);
2138 18249 : if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
2139 18249 : return gerepileupto(av, R);
2140 : }
2141 :
2142 : /* assume x has integral coefficients */
2143 : GEN
2144 1144888 : ZX_disc_all(GEN x, ulong bound)
2145 : {
2146 1144888 : pari_sp av = avma;
2147 1144888 : long s, d = degpol(x);
2148 : GEN l, R;
2149 :
2150 1144885 : if (d <= 1) return d == 1? gen_1: gen_0;
2151 1142005 : s = (d & 2) ? -1: 1;
2152 1142005 : l = leading_coeff(x);
2153 1142010 : if (!bound) bound = ZX_discbound(x);
2154 1141967 : R = ZX_resultant_all(x, NULL, NULL, bound);
2155 1142007 : if (is_pm1(l))
2156 967061 : { if (signe(l) < 0) s = -s; }
2157 : else
2158 174944 : R = diviiexact(R,l);
2159 1142005 : if (s == -1) togglesign_safe(&R);
2160 1142006 : return gerepileuptoint(av,R);
2161 : }
2162 :
2163 : GEN
2164 1082991 : ZX_disc(GEN x) { return ZX_disc_all(x,0); }
2165 :
2166 : static GEN
2167 6887 : ZXQX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, GEN T, ulong p)
2168 : {
2169 6887 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2170 : GEN H, dp;
2171 6887 : if (dropa && dropb) return pol0_Flx(T[1]); /* p | lc(A), p | lc(B) */
2172 6887 : H = FlxqX_saferesultant(a, b, T, p);
2173 6887 : if (!H) return NULL;
2174 6887 : if (dropa)
2175 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2176 0 : GEN c = gel(b,degB+2); /* lc(B) */
2177 0 : if (odd(degB)) c = Flx_neg(c, p);
2178 0 : c = Flxq_powu(c, dropa, T, p);
2179 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2180 : }
2181 6887 : else if (dropb)
2182 : { /* multiply by lc(A)^(deg B - deg b) */
2183 0 : GEN c = gel(a,degA+2); /* lc(A) */
2184 0 : c = Flxq_powu(c, dropb, T, p);
2185 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2186 : }
2187 6887 : dp = dB ? ZX_to_Flx(dB, p): pol1_Flx(T[1]);
2188 6888 : if (!Flx_equal1(dp))
2189 : {
2190 0 : GEN idp = Flxq_invsafe(dp, T, p);
2191 0 : if (!idp) return NULL;
2192 0 : H = Flxq_mul(H, Flxq_powu(idp, degA, T, p), T, p);
2193 : }
2194 6888 : return H;
2195 : }
2196 :
2197 : /* If B=NULL, assume B=A' */
2198 : static GEN
2199 3306 : ZXQX_resultant_slice(GEN A, GEN B, GEN U, GEN dB, GEN P, GEN *mod)
2200 : {
2201 3306 : pari_sp av = avma;
2202 3306 : long degA, degB, i, n = lg(P)-1;
2203 : GEN H, T;
2204 3306 : long v = varn(U), redo = 0;
2205 :
2206 3306 : degA = degpol(A);
2207 3306 : degB = B? degpol(B): degA - 1;
2208 3306 : if (n == 1)
2209 : {
2210 2144 : ulong p = uel(P,1);
2211 2144 : GEN a = ZXX_to_FlxX(A, p, v), b = B? ZXX_to_FlxX(B, p, v): FlxX_deriv(a, p);
2212 2144 : GEN u = ZX_to_Flx(U, p);
2213 2144 : GEN Hp = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2214 2144 : if (!Hp) { set_avma(av); *mod = gen_1; return pol_0(v); }
2215 2144 : Hp = gerepileupto(av, Flx_to_ZX(Hp)); *mod = utoipos(p); return Hp;
2216 : }
2217 1162 : T = ZV_producttree(P);
2218 1162 : A = ZXX_nv_mod_tree(A, P, T, v);
2219 1162 : if (B) B = ZXX_nv_mod_tree(B, P, T, v);
2220 1162 : U = ZX_nv_mod_tree(U, P, T);
2221 1162 : H = cgetg(n+1, t_VEC);
2222 5905 : for(i=1; i <= n; i++)
2223 : {
2224 4743 : ulong p = P[i];
2225 4743 : GEN a = gel(A,i), b = B? gel(B,i): FlxX_deriv(a, p), u = gel(U, i);
2226 4743 : GEN h = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2227 4744 : if (!h)
2228 : {
2229 0 : gel(H,i) = pol_0(v);
2230 0 : P[i] = 1; redo = 1;
2231 : }
2232 : else
2233 4744 : gel(H,i) = h;
2234 : }
2235 1162 : if (redo) T = ZV_producttree(P);
2236 1162 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2237 1162 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2238 : }
2239 :
2240 : GEN
2241 3306 : ZXQX_resultant_worker(GEN P, GEN A, GEN B, GEN T, GEN dB)
2242 : {
2243 3306 : GEN V = cgetg(3, t_VEC);
2244 3306 : if (isintzero(B)) B = NULL;
2245 3306 : if (!signe(dB)) dB = NULL;
2246 3306 : gel(V,1) = ZXQX_resultant_slice(A, B, T, dB, P, &gel(V,2));
2247 3306 : return V;
2248 : }
2249 :
2250 : static ulong
2251 2999 : ZXQX_resultant_bound_i(GEN nf, GEN A, GEN B, GEN (*f)(GEN,GEN,long))
2252 : {
2253 2999 : pari_sp av = avma;
2254 2999 : GEN r, M = nf_L2_bound(nf, NULL, &r);
2255 2999 : long v = nf_get_varn(nf), i, l = lg(r);
2256 2999 : GEN a = cgetg(l, t_COL);
2257 9484 : for (i = 1; i < l; i++)
2258 6485 : gel(a, i) = f(gsubst(A, v, gel(r,i)), gsubst(B, v, gel(r,i)), DEFAULTPREC);
2259 2999 : return gc_ulong(av, (ulong) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
2260 : }
2261 : static ulong
2262 2747 : ZXQX_resultant_bound(GEN nf, GEN A, GEN B)
2263 2747 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound); }
2264 :
2265 : /* Compute Res(A, B/dB) in Z[X]/T, assuming A,B in Z[X,Y], dB in Z or NULL (= 1)
2266 : * If B=NULL, take B = A' and assume deg A > 1 */
2267 : static GEN
2268 2744 : ZXQX_resultant_all(GEN A, GEN B, GEN T, GEN dB, ulong bound)
2269 : {
2270 2744 : pari_sp av = avma;
2271 : forprime_t S;
2272 : GEN H, worker;
2273 2744 : if (B)
2274 : {
2275 70 : long a = degpol(A), b = degpol(B);
2276 70 : if (a < 0 || b < 0) return gen_0;
2277 70 : if (!a) return gpowgs(gel(A,2), b);
2278 70 : if (!b) return gpowgs(gel(B,2), a);
2279 : } else
2280 2674 : if (!bound) B = RgX_deriv(A);
2281 2723 : if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, B);
2282 2723 : worker = snm_closure(is_entry("_ZXQX_resultant_worker"),
2283 : mkvec4(A, B? B: gen_0, T, dB? dB: gen_0));
2284 2723 : init_modular_big(&S);
2285 2723 : H = gen_crt("ZXQX_resultant_all", worker, &S, dB, bound, 0, NULL,
2286 : nxV_chinese_center, FpX_center);
2287 2723 : if (DEBUGLEVEL)
2288 0 : err_printf("ZXQX_resultant_all: a priori bound: %lu, a posteriori: %lu\n",
2289 : bound, expi(gsupnorm(H, DEFAULTPREC)));
2290 2723 : return gerepileupto(av, H);
2291 : }
2292 :
2293 : GEN
2294 70 : nfX_resultant(GEN nf, GEN x, GEN y)
2295 : {
2296 70 : pari_sp av = avma;
2297 70 : GEN cx, cy, D, T = nf_get_pol(nf);
2298 : ulong bound;
2299 70 : long d = degpol(x), v = varn(T);
2300 70 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2301 70 : x = Q_primitive_part(x, &cx);
2302 70 : y = Q_primitive_part(y, &cy);
2303 70 : bound = ZXQX_resultant_bound(nf, x, y);
2304 70 : D = ZXQX_resultant_all(x, y, T, NULL, bound);
2305 70 : if (cx) D = gmul(D, gpowgs(cx, degpol(y)));
2306 70 : if (cy) D = gmul(D, gpowgs(cy, degpol(x)));
2307 70 : return gerepileupto(av, D);
2308 : }
2309 :
2310 : static GEN
2311 182 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2312 :
2313 : static GEN
2314 2674 : ZXQX_disc_all(GEN x, GEN T, ulong bound)
2315 : {
2316 2674 : pari_sp av = avma;
2317 2674 : long s, d = degpol(x), v = varn(T);
2318 : GEN l, R;
2319 :
2320 2674 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2321 2674 : s = (d & 2) ? -1: 1;
2322 2674 : l = leading_coeff(x);
2323 2674 : R = ZXQX_resultant_all(x, NULL, T, NULL, bound);
2324 2674 : if (!gequal1(l)) R = QXQ_div(R, to_ZX(l,v), T);
2325 2674 : if (s == -1) R = RgX_neg(R);
2326 2674 : return gerepileupto(av, R);
2327 : }
2328 :
2329 : GEN
2330 7 : QX_disc(GEN x)
2331 : {
2332 7 : pari_sp av = avma;
2333 7 : GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
2334 7 : if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
2335 7 : return gerepileupto(av, d);
2336 : }
2337 :
2338 : GEN
2339 2821 : nfX_disc(GEN nf, GEN x)
2340 : {
2341 2821 : pari_sp av = avma;
2342 2821 : GEN c, D, T = nf_get_pol(nf);
2343 : ulong bound;
2344 2821 : long d = degpol(x), v = varn(T);
2345 2821 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2346 2674 : x = Q_primitive_part(x, &c);
2347 2674 : bound = ZXQX_resultant_bound(nf, x, RgX_deriv(x));
2348 2674 : D = ZXQX_disc_all(x, T, bound);
2349 2674 : if (c) D = gmul(D, gpowgs(c, 2*d - 2));
2350 2674 : return gerepileupto(av, D);
2351 : }
2352 :
2353 : GEN
2354 478032 : QXQ_mul(GEN x, GEN y, GEN T)
2355 : {
2356 478032 : GEN dx, nx = Q_primitive_part(x, &dx);
2357 478032 : GEN dy, ny = Q_primitive_part(y, &dy);
2358 478031 : GEN z = ZXQ_mul(nx, ny, T);
2359 478032 : if (dx || dy)
2360 : {
2361 475401 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
2362 475401 : if (!gequal1(d)) z = ZX_Q_mul(z, d);
2363 : }
2364 478033 : return z;
2365 : }
2366 :
2367 : GEN
2368 76855 : QXQ_sqr(GEN x, GEN T)
2369 : {
2370 76855 : GEN dx, nx = Q_primitive_part(x, &dx);
2371 76855 : GEN z = ZXQ_sqr(nx, T);
2372 76855 : if (dx)
2373 75504 : z = ZX_Q_mul(z, gsqr(dx));
2374 76855 : return z;
2375 : }
2376 :
2377 : static GEN
2378 60564 : QXQ_inv_slice(GEN A, GEN B, GEN P, GEN *mod)
2379 : {
2380 60564 : pari_sp av = avma;
2381 60564 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2382 : GEN H, T;
2383 60564 : if (n == 1)
2384 : {
2385 42391 : ulong p = uel(P,1);
2386 42391 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2387 42391 : GEN U = Flxq_invsafe(a, b, p);
2388 42391 : if (!U)
2389 : {
2390 24 : set_avma(av);
2391 24 : *mod = gen_1; return pol_0(v);
2392 : }
2393 42367 : H = gerepilecopy(av, Flx_to_ZX(U));
2394 42367 : *mod = utoipos(p); return H;
2395 : }
2396 18173 : T = ZV_producttree(P);
2397 18171 : A = ZX_nv_mod_tree(A, P, T);
2398 18172 : B = ZX_nv_mod_tree(B, P, T);
2399 18172 : H = cgetg(n+1, t_VEC);
2400 89785 : for(i=1; i <= n; i++)
2401 : {
2402 71613 : ulong p = P[i];
2403 71613 : GEN a = gel(A,i), b = gel(B,i);
2404 71613 : GEN U = Flxq_invsafe(a, b, p);
2405 71600 : if (!U)
2406 : {
2407 601 : gel(H,i) = pol_0(v);
2408 601 : P[i] = 1; redo = 1;
2409 : }
2410 : else
2411 70999 : gel(H,i) = U;
2412 : }
2413 18172 : if (redo) T = ZV_producttree(P);
2414 18172 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2415 18172 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2416 : }
2417 :
2418 : GEN
2419 60564 : QXQ_inv_worker(GEN P, GEN A, GEN B)
2420 : {
2421 60564 : GEN V = cgetg(3, t_VEC);
2422 60564 : gel(V,1) = QXQ_inv_slice(A, B, P, &gel(V,2));
2423 60564 : return V;
2424 : }
2425 :
2426 : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
2427 : GEN
2428 29836 : QXQ_inv(GEN A, GEN B)
2429 : {
2430 : GEN D, Ap, Bp;
2431 : ulong pp;
2432 29836 : pari_sp av2, av = avma;
2433 : forprime_t S;
2434 29836 : GEN worker, U, H = NULL, mod = gen_1;
2435 : pari_timer ti;
2436 : long k, dA, dB;
2437 29836 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2438 : /* A a QX, B a ZX */
2439 29836 : A = Q_primitive_part(A, &D);
2440 29836 : dA = degpol(A); dB= degpol(B);
2441 : /* A, B in Z[X] */
2442 29836 : init_modular_small(&S);
2443 : do {
2444 29836 : pp = u_forprime_next(&S);
2445 29836 : Ap = ZX_to_Flx(A, pp);
2446 29836 : Bp = ZX_to_Flx(B, pp);
2447 29836 : } while (degpol(Ap) != dA || degpol(Bp) != dB);
2448 29836 : if (degpol(Flx_gcd(Ap, Bp, pp)) != 0 && degpol(ZX_gcd(A,B))!=0)
2449 14 : pari_err_INV("QXQ_inv",mkpolmod(A,B));
2450 29822 : worker = snm_closure(is_entry("_QXQ_inv_worker"), mkvec2(A, B));
2451 29822 : av2 = avma;
2452 29822 : for (k = 1; ;k *= 2)
2453 22644 : {
2454 : GEN res, b, N, den;
2455 52466 : gen_inccrt_i("QXQ_inv", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2456 : nxV_chinese_center, FpX_center);
2457 52466 : gerepileall(av2, 2, &H, &mod);
2458 52466 : b = sqrti(shifti(mod,-1));
2459 52466 : if (DEBUGLEVEL>5) timer_start(&ti);
2460 52466 : U = FpX_ratlift(H, mod, b, b, NULL);
2461 52466 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: ratlift");
2462 55491 : if (!U) continue;
2463 32847 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2464 32847 : res = Flx_rem(Flx_Fl_sub(Flx_mul(Ap, ZX_to_Flx(N,pp), pp),
2465 : umodiu(den, pp), pp), Bp, pp);
2466 32847 : if (degpol(res) >= 0) continue;
2467 29822 : res = ZX_Z_sub(ZX_mul(A, N), den);
2468 29822 : res = ZX_is_monic(B) ? ZX_rem(res, B): RgX_pseudorem(res, B);
2469 29822 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: final check");
2470 29822 : if (degpol(res)<0)
2471 : {
2472 29822 : if (D) U = RgX_Rg_div(U, D);
2473 29822 : return gerepilecopy(av, U);
2474 : }
2475 : }
2476 : }
2477 :
2478 : static GEN
2479 115095 : QXQ_div_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
2480 : {
2481 115095 : pari_sp av = avma;
2482 115095 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2483 : GEN H, T;
2484 115095 : if (n == 1)
2485 : {
2486 42037 : ulong p = uel(P,1);
2487 42037 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), c = ZX_to_Flx(C, p);
2488 42037 : GEN bi = Flxq_invsafe(b, c, p), U;
2489 42037 : if (!bi)
2490 : {
2491 0 : set_avma(av);
2492 0 : *mod = gen_1; return pol_0(v);
2493 : }
2494 42037 : U = Flxq_mul(a, bi, c, p);
2495 42037 : H = gerepilecopy(av, Flx_to_ZX(U));
2496 42038 : *mod = utoipos(p); return H;
2497 : }
2498 73058 : T = ZV_producttree(P);
2499 73056 : A = ZX_nv_mod_tree(A, P, T);
2500 73052 : B = ZX_nv_mod_tree(B, P, T);
2501 73055 : C = ZX_nv_mod_tree(C, P, T);
2502 73056 : H = cgetg(n+1, t_VEC);
2503 322242 : for(i=1; i <= n; i++)
2504 : {
2505 249185 : ulong p = P[i];
2506 249185 : GEN a = gel(A,i), b = gel(B,i), c = gel(C, i);
2507 249185 : GEN bi = Flxq_invsafe(b, c, p);
2508 249193 : if (!bi)
2509 : {
2510 0 : gel(H,i) = pol_0(v);
2511 0 : P[i] = 1; redo = 1;
2512 : }
2513 : else
2514 249193 : gel(H,i) = Flxq_mul(a, bi, c, p);
2515 : }
2516 73057 : if (redo) T = ZV_producttree(P);
2517 73057 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2518 73059 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2519 : }
2520 :
2521 : GEN
2522 115095 : QXQ_div_worker(GEN P, GEN A, GEN B, GEN C)
2523 : {
2524 115095 : GEN V = cgetg(3, t_VEC);
2525 115095 : gel(V,1) = QXQ_div_slice(A, B, C, P, &gel(V,2));
2526 115095 : return V;
2527 : }
2528 :
2529 : /* lift(Mod(A/B, C)). C a ZX, A, B a scalar or a QX */
2530 : GEN
2531 31248 : QXQ_div(GEN A, GEN B, GEN C)
2532 : {
2533 : GEN DA, DB, Ap, Bp, Cp;
2534 : ulong pp;
2535 31248 : pari_sp av2, av = avma;
2536 : forprime_t S;
2537 31248 : GEN worker, U, H = NULL, mod = gen_1;
2538 : pari_timer ti;
2539 : long k, dA, dB, dC;
2540 31248 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2541 : /* A a QX, B a ZX */
2542 31248 : A = Q_primitive_part(A, &DA);
2543 31248 : B = Q_primitive_part(B, &DB);
2544 31248 : dA = degpol(A); dB = degpol(B); dC = degpol(C);
2545 : /* A, B in Z[X] */
2546 31248 : init_modular_small(&S);
2547 : do {
2548 31248 : pp = u_forprime_next(&S);
2549 31248 : Ap = ZX_to_Flx(A, pp);
2550 31248 : Bp = ZX_to_Flx(B, pp);
2551 31248 : Cp = ZX_to_Flx(C, pp);
2552 31248 : } while (degpol(Ap) != dA || degpol(Bp) != dB || degpol(Cp) != dC);
2553 31248 : if (degpol(Flx_gcd(Bp, Cp, pp)) != 0 && degpol(ZX_gcd(B,C))!=0)
2554 0 : pari_err_INV("QXQ_div",mkpolmod(B,C));
2555 31248 : worker = snm_closure(is_entry("_QXQ_div_worker"), mkvec3(A, B, C));
2556 31248 : av2 = avma;
2557 31248 : for (k = 1; ;k *= 2)
2558 44377 : {
2559 : GEN res, b, N, den;
2560 75625 : gen_inccrt_i("QXQ_div", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2561 : nxV_chinese_center, FpX_center);
2562 75624 : gerepileall(av2, 2, &H, &mod);
2563 75624 : b = sqrti(shifti(mod,-1));
2564 75624 : if (DEBUGLEVEL>5) timer_start(&ti);
2565 75624 : U = FpX_ratlift(H, mod, b, b, NULL);
2566 75625 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: ratlift");
2567 86077 : if (!U) continue;
2568 41700 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2569 41699 : res = Flx_rem(Flx_sub(Flx_mul(Bp, ZX_to_Flx(N,pp), pp),
2570 : Flx_Fl_mul(Ap, umodiu(den, pp), pp), pp), Cp, pp);
2571 41699 : if (degpol(res) >= 0) continue;
2572 31247 : res = ZX_sub(ZX_mul(B, N), ZX_Z_mul(A,den));
2573 31248 : res = ZX_is_monic(C) ? ZX_rem(res, C): RgX_pseudorem(res, C);
2574 31248 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: final check");
2575 31248 : if (degpol(res)<0)
2576 : {
2577 31248 : if (DA && DB) U = RgX_Rg_mul(U, gdiv(DA,DB));
2578 26656 : else if (DA) U = RgX_Rg_mul(U, DA);
2579 14931 : else if (DB) U = RgX_Rg_div(U, DB);
2580 31247 : return gerepilecopy(av, U);
2581 : }
2582 : }
2583 : }
2584 :
2585 : /************************************************************************
2586 : * *
2587 : * ZXQ_minpoly *
2588 : * *
2589 : ************************************************************************/
2590 :
2591 : static GEN
2592 3523 : ZXQ_minpoly_slice(GEN A, GEN B, long d, GEN P, GEN *mod)
2593 : {
2594 3523 : pari_sp av = avma;
2595 3523 : long i, n = lg(P)-1, v = evalvarn(varn(B));
2596 : GEN H, T;
2597 3523 : if (n == 1)
2598 : {
2599 716 : ulong p = uel(P,1);
2600 716 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2601 716 : GEN Hp = Flxq_minpoly(a, b, p);
2602 716 : if (degpol(Hp) != d) { p = 1; Hp = pol0_Flx(v); }
2603 716 : H = gerepileupto(av, Flx_to_ZX(Hp));
2604 716 : *mod = utoipos(p); return H;
2605 : }
2606 2807 : T = ZV_producttree(P);
2607 2807 : A = ZX_nv_mod_tree(A, P, T);
2608 2807 : B = ZX_nv_mod_tree(B, P, T);
2609 2807 : H = cgetg(n+1, t_VEC);
2610 16838 : for(i=1; i <= n; i++)
2611 : {
2612 14031 : ulong p = P[i];
2613 14031 : GEN a = gel(A,i), b = gel(B,i);
2614 14031 : GEN m = Flxq_minpoly(a, b, p);
2615 14031 : if (degpol(m) != d) { P[i] = 1; m = pol0_Flx(v); }
2616 14031 : gel(H, i) = m;
2617 : }
2618 2807 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2619 2807 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2620 : }
2621 :
2622 : GEN
2623 3523 : ZXQ_minpoly_worker(GEN P, GEN A, GEN B, long d)
2624 : {
2625 3523 : GEN V = cgetg(3, t_VEC);
2626 3523 : gel(V,1) = ZXQ_minpoly_slice(A, B, d, P, &gel(V,2));
2627 3523 : return V;
2628 : }
2629 :
2630 : GEN
2631 1701 : ZXQ_minpoly(GEN A, GEN B, long d, ulong bound)
2632 : {
2633 1701 : pari_sp av = avma;
2634 : GEN worker, H, dB;
2635 : forprime_t S;
2636 1701 : B = Q_remove_denom(B, &dB);
2637 1701 : worker = strtoclosure("_ZXQ_minpoly_worker", 3, A, B, stoi(d));
2638 1701 : init_modular_big(&S);
2639 1701 : H = gen_crt("ZXQ_minpoly", worker, &S, dB, bound, 0, NULL,
2640 : nxV_chinese_center, FpX_center_i);
2641 1701 : return gerepilecopy(av, H);
2642 : }
2643 :
2644 : /************************************************************************
2645 : * *
2646 : * ZX_ZXY_resultant *
2647 : * *
2648 : ************************************************************************/
2649 :
2650 : static GEN
2651 107244 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
2652 : long degA, long degB, long dres, long sX)
2653 : {
2654 107244 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2655 107244 : GEN Hp = Flx_FlxY_resultant_polint(a, b, p, dres, sX);
2656 107245 : if (dropa && dropb)
2657 0 : Hp = zero_Flx(sX);
2658 : else {
2659 107245 : if (dropa)
2660 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2661 0 : GEN c = gel(b,degB+2); /* lc(B) */
2662 0 : if (odd(degB)) c = Flx_neg(c, p);
2663 0 : if (!Flx_equal1(c)) {
2664 0 : c = Flx_powu(c, dropa, p);
2665 0 : if (!Flx_equal1(c)) Hp = Flx_mul(Hp, c, p);
2666 : }
2667 : }
2668 107245 : else if (dropb)
2669 : { /* multiply by lc(A)^(deg B - deg b) */
2670 0 : ulong c = uel(a, degA+2); /* lc(A) */
2671 0 : c = Fl_powu(c, dropb, p);
2672 0 : if (c != 1) Hp = Flx_Fl_mul(Hp, c, p);
2673 : }
2674 : }
2675 107245 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
2676 107245 : return Hp;
2677 : }
2678 :
2679 : static GEN
2680 47098 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
2681 : GEN P, GEN *mod, long sX, long vY)
2682 : {
2683 47098 : pari_sp av = avma;
2684 47098 : long i, n = lg(P)-1;
2685 : GEN H, T, D;
2686 47098 : if (n == 1)
2687 : {
2688 38209 : ulong p = uel(P,1);
2689 38209 : ulong dp = dB ? umodiu(dB, p): 1;
2690 38209 : GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
2691 38209 : GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2692 38208 : H = gerepileupto(av, Flx_to_ZX(Hp));
2693 38207 : *mod = utoipos(p); return H;
2694 : }
2695 8889 : T = ZV_producttree(P);
2696 8889 : A = ZX_nv_mod_tree(A, P, T);
2697 8889 : B = ZXX_nv_mod_tree(B, P, T, vY);
2698 8889 : D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
2699 8889 : H = cgetg(n+1, t_VEC);
2700 32891 : for(i=1; i <= n; i++)
2701 : {
2702 24002 : ulong p = P[i];
2703 24002 : GEN a = gel(A,i), b = gel(B,i);
2704 24002 : ulong dp = D ? uel(D, i): 1;
2705 24002 : gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2706 : }
2707 8889 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2708 8889 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2709 : }
2710 :
2711 : GEN
2712 47098 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
2713 : {
2714 47098 : GEN V = cgetg(3, t_VEC);
2715 47098 : if (isintzero(dB)) dB = NULL;
2716 47098 : gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
2717 47096 : return V;
2718 : }
2719 :
2720 : GEN
2721 44677 : ZX_ZXY_resultant(GEN A, GEN B)
2722 : {
2723 44677 : pari_sp av = avma;
2724 : forprime_t S;
2725 : ulong bound;
2726 44677 : long v = fetch_var_higher();
2727 44677 : long degA = degpol(A), degB, dres = degA * degpol(B);
2728 44677 : long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
2729 44677 : long sX = evalvarn(vX);
2730 : GEN worker, H, dB;
2731 44677 : B = Q_remove_denom(B, &dB);
2732 44677 : if (!dB) B = leafcopy(B);
2733 44677 : A = leafcopy(A); setvarn(A,v);
2734 44677 : B = swap_vars(B, vY); setvarn(B,v); degB = degpol(B);
2735 44677 : bound = ZX_ZXY_ResBound(A, B, dB);
2736 44676 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
2737 89352 : worker = snm_closure(is_entry("_ZX_ZXY_resultant_worker"),
2738 44676 : mkvec4(A, B, dB? dB: gen_0,
2739 : mkvecsmall5(degA, degB, dres, sX, vY)));
2740 44678 : init_modular_big(&S);
2741 44678 : H = gen_crt("ZX_ZXY_resultant_all", worker, &S, dB, bound, 0, NULL,
2742 : nxV_chinese_center, FpX_center_i);
2743 44678 : setvarn(H, vX); (void)delete_var();
2744 44678 : return gerepilecopy(av, H);
2745 : }
2746 :
2747 : static long
2748 40347 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
2749 : {
2750 40347 : pari_sp av = avma;
2751 40347 : long degA = degpol(A), degB, dres = degA*degpol(B0);
2752 40347 : long v = fetch_var_higher();
2753 40347 : long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
2754 40347 : long sX = evalvarn(vX);
2755 : GEN dB, B, a, b, Hp;
2756 : forprime_t S;
2757 :
2758 40347 : B0 = Q_remove_denom(B0, &dB);
2759 40347 : if (!dB) B0 = leafcopy(B0);
2760 40347 : A = leafcopy(A);
2761 40347 : B = B0;
2762 40347 : setvarn(A,v);
2763 45036 : INIT:
2764 45036 : if (lambda) B = RgX_translate(B0, monomial(stoi(lambda), 1, vY));
2765 45036 : B = swap_vars(B, vY); setvarn(B,v);
2766 : /* B0(lambda v + x, v) */
2767 45036 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
2768 :
2769 45036 : degB = degpol(B);
2770 45036 : init_modular_big(&S);
2771 : while (1)
2772 0 : {
2773 45036 : ulong p = u_forprime_next(&S);
2774 45036 : ulong dp = dB ? umodiu(dB, p): 1;
2775 45036 : if (!dp) continue;
2776 45036 : a = ZX_to_Flx(A, p);
2777 45036 : b = ZXX_to_FlxX(B, p, v);
2778 45035 : Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2779 45036 : if (degpol(Hp) != dres) continue;
2780 45036 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
2781 45036 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
2782 40346 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
2783 40346 : (void)delete_var(); return gc_long(av,lambda);
2784 : }
2785 : }
2786 :
2787 : GEN
2788 41192 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
2789 : {
2790 41192 : if (lambda)
2791 : {
2792 40347 : *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
2793 40346 : if (*lambda) B = RgX_translate(B, monomial(stoi(*lambda), 1, varn(A)));
2794 : }
2795 41191 : return ZX_ZXY_resultant(A,B);
2796 : }
2797 :
2798 : static GEN
2799 10368 : ZX_direct_compositum_slice(GEN A, GEN B, GEN P, GEN *mod)
2800 : {
2801 10368 : pari_sp av = avma;
2802 10368 : long i, n = lg(P)-1;
2803 : GEN H, T;
2804 10368 : if (n == 1)
2805 : {
2806 9866 : ulong p = uel(P,1);
2807 9866 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2808 9862 : GEN Hp = Flx_direct_compositum(a, b, p);
2809 9864 : H = gerepileupto(av, Flx_to_ZX(Hp));
2810 9870 : *mod = utoipos(p); return H;
2811 : }
2812 502 : T = ZV_producttree(P);
2813 502 : A = ZX_nv_mod_tree(A, P, T);
2814 502 : B = ZX_nv_mod_tree(B, P, T);
2815 502 : H = cgetg(n+1, t_VEC);
2816 4526 : for(i=1; i <= n; i++)
2817 : {
2818 4024 : ulong p = P[i];
2819 4024 : GEN a = gel(A,i), b = gel(B,i);
2820 4024 : gel(H,i) = Flx_direct_compositum(a, b, p);
2821 : }
2822 502 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2823 502 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2824 : }
2825 :
2826 : GEN
2827 10367 : ZX_direct_compositum_worker(GEN P, GEN A, GEN B)
2828 : {
2829 10367 : GEN V = cgetg(3, t_VEC);
2830 10368 : gel(V,1) = ZX_direct_compositum_slice(A, B, P, &gel(V,2));
2831 10370 : return V;
2832 : }
2833 :
2834 : /* Assume A,B irreducible (in particular squarefree) and define linearly
2835 : * disjoint extensions: no factorisation needed */
2836 : static GEN
2837 10102 : ZX_direct_compositum(GEN A, GEN B, GEN lead)
2838 : {
2839 10102 : pari_sp av = avma;
2840 : forprime_t S;
2841 : ulong bound;
2842 : GEN H, worker, mod;
2843 10102 : if (degpol(A) < degpol(B)) swap(A, B);
2844 10102 : bound = ZX_ZXY_ResBound_1(A, B);
2845 10108 : worker = snm_closure(is_entry("_ZX_direct_compositum_worker"), mkvec2(A,B));
2846 10107 : init_modular_big(&S);
2847 10106 : H = gen_crt("ZX_direct_compositum", worker, &S, lead, bound, 0, &mod,
2848 : nxV_chinese_center, FpX_center);
2849 10105 : return gerepileupto(av, H);
2850 : }
2851 :
2852 : static long
2853 9718 : ZX_compositum_lambda(GEN A, GEN B, GEN lead, long lambda)
2854 : {
2855 9718 : pari_sp av = avma;
2856 : forprime_t S;
2857 : ulong p;
2858 9718 : init_modular_big(&S);
2859 9721 : p = u_forprime_next(&S);
2860 : while (1)
2861 112 : {
2862 : GEN Hp, a;
2863 9833 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
2864 9833 : if (lead && dvdiu(lead,p)) { p = u_forprime_next(&S); continue; }
2865 9825 : a = ZX_to_Flx(ZX_rescale(A, stoi(-lambda)), p);
2866 9825 : Hp = Flx_direct_compositum(a, ZX_to_Flx(B, p), p);
2867 9822 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); continue; }
2868 9716 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
2869 9716 : return gc_long(av, lambda);
2870 : }
2871 : }
2872 :
2873 : GEN
2874 10105 : ZX_compositum(GEN A, GEN B, long *lambda)
2875 : {
2876 10105 : GEN lead = mulii(leading_coeff(A),leading_coeff(B));
2877 10104 : if (lambda)
2878 : {
2879 9719 : *lambda = ZX_compositum_lambda(A, B, lead, *lambda);
2880 9716 : A = ZX_rescale(A, stoi(-*lambda));
2881 : }
2882 10101 : return ZX_direct_compositum(A, B, lead);
2883 : }
2884 :
2885 : GEN
2886 385 : ZX_compositum_disjoint(GEN A, GEN B)
2887 385 : { return ZX_compositum(A, B, NULL); }
2888 :
2889 : static GEN
2890 287 : ZXQX_direct_compositum_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
2891 : {
2892 287 : pari_sp av = avma;
2893 287 : long i, n = lg(P)-1, dC = degpol(C), v = varn(C);
2894 : GEN H, T;
2895 287 : if (n == 1)
2896 : {
2897 174 : ulong p = uel(P,1);
2898 174 : GEN a = ZXX_to_FlxX(A, p, v), b = ZXX_to_FlxX(B, p, v);
2899 174 : GEN c = ZX_to_Flx(C, p);
2900 174 : GEN Hp = FlxX_to_Flm(FlxqX_direct_compositum(a, b, c, p), dC);
2901 174 : H = gerepileupto(av, Flm_to_ZM(Hp));
2902 174 : *mod = utoipos(p); return H;
2903 : }
2904 113 : T = ZV_producttree(P);
2905 113 : A = ZXX_nv_mod_tree(A, P, T, v);
2906 113 : B = ZXX_nv_mod_tree(B, P, T, v);
2907 113 : C = ZX_nv_mod_tree(C, P, T);
2908 113 : H = cgetg(n+1, t_VEC);
2909 414 : for(i=1; i <= n; i++)
2910 : {
2911 301 : ulong p = P[i];
2912 301 : GEN a = gel(A,i), b = gel(B,i), c = gel(C,i);
2913 301 : gel(H,i) = FlxX_to_Flm(FlxqX_direct_compositum(a, b, c, p), dC);
2914 : }
2915 113 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
2916 113 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2917 : }
2918 :
2919 : GEN
2920 287 : ZXQX_direct_compositum_worker(GEN P, GEN A, GEN B, GEN C)
2921 : {
2922 287 : GEN V = cgetg(3, t_VEC);
2923 287 : gel(V,1) = ZXQX_direct_compositum_slice(A, B, C, P, &gel(V,2));
2924 287 : return V;
2925 : }
2926 :
2927 : static GEN
2928 252 : ZXQX_direct_compositum(GEN A, GEN B, GEN T, ulong bound)
2929 : {
2930 252 : pari_sp av = avma;
2931 : forprime_t S;
2932 : GEN H, worker, mod;
2933 252 : GEN lead = mulii(Q_content(leading_coeff(A)), Q_content(leading_coeff(B)));
2934 252 : worker = snm_closure(is_entry("_ZXQX_direct_compositum_worker")
2935 : , mkvec3(A,B,T));
2936 252 : init_modular_big(&S);
2937 252 : H = gen_crt("ZXQX_direct_compositum", worker, &S, lead, bound, 0, &mod,
2938 : nmV_chinese_center, FpM_center);
2939 252 : if (DEBUGLEVEL > 4)
2940 0 : err_printf("nfcompositum: a priori bound: %lu, a posteriori: %lu\n",
2941 : bound, expi(gsupnorm(H, DEFAULTPREC)));
2942 252 : return gerepilecopy(av, RgM_to_RgXX(H, varn(A), varn(T)));
2943 : }
2944 :
2945 : static long
2946 252 : ZXQX_direct_compositum_bound(GEN nf, GEN A, GEN B)
2947 252 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound_1); }
2948 :
2949 : GEN
2950 252 : nf_direct_compositum(GEN nf, GEN A, GEN B)
2951 : {
2952 252 : ulong bnd = ZXQX_direct_compositum_bound(nf, A, B);
2953 252 : return ZXQX_direct_compositum(A, B, nf_get_pol(nf), bnd);
2954 : }
2955 :
2956 : /************************************************************************
2957 : * *
2958 : * IRREDUCIBLE POLYNOMIAL / Fp *
2959 : * *
2960 : ************************************************************************/
2961 :
2962 : /* irreducible (unitary) polynomial of degree n over Fp */
2963 : GEN
2964 0 : ffinit_rand(GEN p,long n)
2965 : {
2966 0 : for(;;) {
2967 0 : pari_sp av = avma;
2968 0 : GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
2969 0 : if (FpX_is_irred(pol, p)) return pol;
2970 0 : set_avma(av);
2971 : }
2972 : }
2973 :
2974 : /* return an extension of degree 2^l of F_2, assume l > 0
2975 : * Not stack clean. */
2976 : static GEN
2977 590 : ffinit_Artin_Schreier_2(long l)
2978 : {
2979 : GEN Q, T, S;
2980 : long i, v;
2981 :
2982 590 : if (l == 1) return mkvecsmall4(0,1,1,1); /*x^2 + x + 1*/
2983 555 : v = fetch_var_higher();
2984 556 : S = mkvecsmall5(0, 0, 0, 1, 1); /* y(y^2 + y) */
2985 555 : Q = mkpoln(3, pol1_Flx(0), pol1_Flx(0), S); /* x^2 + x + y(y^2+y) */
2986 556 : setvarn(Q, v);
2987 :
2988 : /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
2989 556 : T = mkvecsmalln(6,evalvarn(v),1UL,1UL,0UL,0UL,1UL);
2990 : /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
2991 : * ==> x^2 + x + a(y) b irred. over K for any root b of Q
2992 : * ==> x^2 + x + (b^2+b)b */
2993 3182 : for (i=2; i<l; i++) T = Flx_FlxY_resultant(T, Q, 2); /* minpoly(b) / F2*/
2994 556 : (void)delete_var(); T[1] = 0; return T;
2995 : }
2996 :
2997 : /* return an extension of degree p^l of F_p, assume l > 0
2998 : * Not stack clean. */
2999 : GEN
3000 787 : ffinit_Artin_Schreier(ulong p, long l)
3001 : {
3002 : long i, v;
3003 : GEN Q, R, S, T, xp;
3004 787 : if (p==2) return ffinit_Artin_Schreier_2(l);
3005 196 : xp = polxn_Flx(p,0); /* x^p */
3006 196 : T = Flx_sub(xp, mkvecsmall3(0,1,1),p); /* x^p - x - 1 */
3007 196 : if (l == 1) return T;
3008 :
3009 7 : v = evalvarn(fetch_var_higher());
3010 7 : xp[1] = v;
3011 7 : R = Flx_sub(polxn_Flx(2*p-1,0), polxn_Flx(p,0),p);
3012 7 : S = Flx_sub(xp, polx_Flx(0), p);
3013 7 : Q = FlxX_Flx_sub(Flx_to_FlxX(S, v), R, p); /* x^p - x - (y^(2p-1)-y^p) */
3014 14 : for (i = 2; i <= l; ++i) T = Flx_FlxY_resultant(T, Q, p);
3015 7 : (void)delete_var(); T[1] = 0; return T;
3016 : }
3017 :
3018 : static long
3019 20859 : flinit_check(ulong p, long n, long l)
3020 : {
3021 : ulong q;
3022 20859 : if (!uisprime(n)) return 0;
3023 12862 : q = p % n; if (!q) return 0;
3024 10783 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3025 : }
3026 :
3027 : static GEN
3028 5038 : flinit(ulong p, long l)
3029 : {
3030 5038 : ulong n = 1+l;
3031 13597 : while (!flinit_check(p,n,l)) n += l;
3032 5038 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3033 5038 : return ZX_to_Flx(polsubcyclo(n,l,0), p);
3034 : }
3035 :
3036 : static GEN
3037 5196 : ffinit_fact_Flx(ulong p, long n)
3038 : {
3039 5196 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3040 5197 : long i, l = lg(Fm);
3041 5197 : P = cgetg(l, t_VEC);
3042 11022 : for (i = 1; i < l; ++i)
3043 5825 : gel(P,i) = p==uel(Fp,i) ?
3044 787 : ffinit_Artin_Schreier(uel(Fp,i), Fe[i])
3045 5825 : : flinit(p, uel(Fm,i));
3046 5197 : return FlxV_direct_compositum(P, p);
3047 : }
3048 :
3049 : static GEN
3050 7261 : init_Flxq_i(ulong p, long n, long sv)
3051 : {
3052 : GEN P;
3053 7261 : if (n == 1) return polx_Flx(sv);
3054 7261 : if (flinit_check(p, n+1, n))
3055 : {
3056 2065 : P = const_vecsmall(n+2,1);
3057 2065 : P[1] = sv; return P;
3058 : }
3059 5196 : P = ffinit_fact_Flx(p,n);
3060 5197 : P[1] = sv; return P;
3061 : }
3062 :
3063 : GEN
3064 0 : init_Flxq(ulong p, long n, long v)
3065 : {
3066 0 : pari_sp av = avma;
3067 0 : return gerepileupto(av, init_Flxq_i(p, n, v));
3068 : }
3069 :
3070 : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
3071 : static long
3072 2188 : fpinit_check(GEN p, long n, long l)
3073 : {
3074 : ulong q;
3075 2188 : if (!uisprime(n)) return 0;
3076 1537 : q = umodiu(p,n); if (!q) return 0;
3077 1537 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3078 : }
3079 :
3080 : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
3081 : * Return an irreducible polynomial of degree l over F_p.
3082 : * Variant of Adleman and Lenstra "Finding irreducible polynomials over
3083 : * finite fields", ACM, 1986 (5) 350--355.
3084 : * Not stack clean */
3085 : static GEN
3086 527 : fpinit(GEN p, long l)
3087 : {
3088 527 : ulong n = 1+l;
3089 1600 : while (!fpinit_check(p,n,l)) n += l;
3090 527 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3091 527 : return FpX_red(polsubcyclo(n,l,0),p);
3092 : }
3093 :
3094 : static GEN
3095 497 : ffinit_fact(GEN p, long n)
3096 : {
3097 497 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3098 497 : long i, l = lg(Fm);
3099 497 : P = cgetg(l, t_VEC);
3100 1024 : for (i = 1; i < l; ++i)
3101 1054 : gel(P,i) = absequaliu(p, Fp[i]) ?
3102 0 : Flx_to_ZX(ffinit_Artin_Schreier(Fp[i], Fe[i]))
3103 527 : : fpinit(p, Fm[i]);
3104 497 : return FpXV_direct_compositum(P, p);
3105 : }
3106 :
3107 : static GEN
3108 8101 : init_Fq_i(GEN p, long n, long v)
3109 : {
3110 : GEN P;
3111 8101 : if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
3112 8101 : if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
3113 8101 : if (signe(p) <= 0) pari_err_PRIME("ffinit",p);
3114 8101 : if (v < 0) v = 0;
3115 8101 : if (n == 1) return pol_x(v);
3116 7849 : if (lgefint(p) == 3)
3117 7261 : return Flx_to_ZX(init_Flxq_i(p[2], n, evalvarn(v)));
3118 588 : if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
3119 497 : P = ffinit_fact(p,n);
3120 497 : setvarn(P, v); return P;
3121 : }
3122 : GEN
3123 7646 : init_Fq(GEN p, long n, long v)
3124 : {
3125 7646 : pari_sp av = avma;
3126 7646 : return gerepileupto(av, init_Fq_i(p, n, v));
3127 : }
3128 : GEN
3129 455 : ffinit(GEN p, long n, long v)
3130 : {
3131 455 : pari_sp av = avma;
3132 455 : return gerepileupto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
3133 : }
3134 :
3135 : GEN
3136 3178 : ffnbirred(GEN p, long n)
3137 : {
3138 3178 : pari_sp av = avma;
3139 3178 : GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
3140 3178 : long j, l = lg(D);
3141 6797 : for (j = 2; j < l; j++) /* skip d = 1 */
3142 : {
3143 3619 : long md = D[j]; /* mu(d) * d, d squarefree */
3144 3619 : GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
3145 3619 : s = md > 0? addii(s, pd): subii(s,pd);
3146 : }
3147 3178 : return gerepileuptoint(av, diviuexact(s, n));
3148 : }
3149 :
3150 : GEN
3151 616 : ffsumnbirred(GEN p, long n)
3152 : {
3153 616 : pari_sp av = avma, av2;
3154 616 : GEN q, t = p, v = vecfactoru_i(1, n);
3155 : long i;
3156 616 : q = cgetg(n+1,t_VEC); gel(q,1) = p;
3157 1764 : for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
3158 616 : av2 = avma;
3159 1764 : for (i=2; i<=n; i++)
3160 : {
3161 1148 : GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
3162 1148 : long j, l = lg(D);
3163 2534 : for (j = 2; j < l; j++) /* skip 1 */
3164 : {
3165 1386 : long md = D[j];
3166 1386 : GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
3167 1386 : s = md > 0? addii(s, pd): subii(s, pd);
3168 : }
3169 1148 : t = gerepileuptoint(av2, addii(t, diviuexact(s, i)));
3170 : }
3171 616 : return gerepileuptoint(av, t);
3172 : }
3173 :
3174 : GEN
3175 140 : ffnbirred0(GEN p, long n, long flag)
3176 : {
3177 140 : if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
3178 140 : if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
3179 140 : switch(flag)
3180 : {
3181 70 : case 0: return ffnbirred(p, n);
3182 70 : case 1: return ffsumnbirred(p, n);
3183 : }
3184 0 : pari_err_FLAG("ffnbirred");
3185 : return NULL; /* LCOV_EXCL_LINE */
3186 : }
3187 :
3188 : static void
3189 2254 : checkmap(GEN m, const char *s)
3190 : {
3191 2254 : if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
3192 0 : pari_err_TYPE(s,m);
3193 2254 : }
3194 :
3195 : GEN
3196 182 : ffembed(GEN a, GEN b)
3197 : {
3198 182 : pari_sp av = avma;
3199 182 : GEN p, Ta, Tb, g, r = NULL;
3200 182 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
3201 182 : if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
3202 182 : p = FF_p_i(a); g = FF_gen(a);
3203 182 : if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
3204 182 : Ta = FF_mod(a);
3205 182 : Tb = FF_mod(b);
3206 182 : if (degpol(Tb)%degpol(Ta)!=0)
3207 7 : pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
3208 175 : r = gel(FFX_roots(Ta, b), 1);
3209 175 : return gerepilecopy(av, mkvec2(g,r));
3210 : }
3211 :
3212 : GEN
3213 91 : ffextend(GEN a, GEN P, long v)
3214 : {
3215 91 : pari_sp av = avma;
3216 : long n;
3217 : GEN p, T, R, g, m;
3218 91 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
3219 91 : T = a; p = FF_p_i(a);
3220 91 : if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
3221 49 : if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
3222 49 : if (v < 0) v = varn(P);
3223 49 : n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
3224 49 : m = ffembed(a, g);
3225 49 : R = FFX_roots(ffmap(m, P),g);
3226 49 : return gerepilecopy(av, mkvec2(gel(R,1), m));
3227 : }
3228 :
3229 : GEN
3230 42 : fffrobenius(GEN a, long n)
3231 : {
3232 42 : if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
3233 42 : retmkvec2(FF_gen(a), FF_Frobenius(a, n));
3234 : }
3235 :
3236 : GEN
3237 133 : ffinvmap(GEN m)
3238 : {
3239 133 : pari_sp av = avma;
3240 : long i, l;
3241 133 : GEN T, F, a, g, r, f = NULL;
3242 133 : checkmap(m, "ffinvmap");
3243 133 : a = gel(m,1); r = gel(m,2);
3244 133 : if (typ(r) != t_FFELT)
3245 7 : pari_err_TYPE("ffinvmap", m);
3246 126 : g = FF_gen(a);
3247 126 : T = FF_mod(r);
3248 126 : F = gel(FFX_factor(T, a), 1);
3249 126 : l = lg(F);
3250 490 : for(i=1; i<l; i++)
3251 : {
3252 490 : GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
3253 490 : if (degpol(s)==0 && gequal(constant_coeff(s),g)) { f = gel(F, i); break; }
3254 : }
3255 126 : if (f==NULL) pari_err_TYPE("ffinvmap", m);
3256 126 : if (degpol(f)==1) f = FF_neg_i(gel(f,2));
3257 126 : return gerepilecopy(av, mkvec2(FF_gen(r),f));
3258 : }
3259 :
3260 : static GEN
3261 1260 : ffpartmapimage(const char *s, GEN r)
3262 : {
3263 1260 : GEN a = NULL, p = NULL;
3264 1260 : if (typ(r)==t_POL && degpol(r) >= 1
3265 1260 : && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
3266 0 : pari_err_TYPE(s, r);
3267 : return NULL; /* LCOV_EXCL_LINE */
3268 : }
3269 :
3270 : static GEN
3271 2702 : ffeltmap_i(GEN m, GEN x)
3272 : {
3273 2702 : GEN r = gel(m,2);
3274 2702 : if (!FF_samefield(x, gel(m,1)))
3275 84 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3276 2618 : if (typ(r)==t_FFELT)
3277 1652 : return FF_map(r, x);
3278 : else
3279 966 : return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
3280 : }
3281 :
3282 : static GEN
3283 4452 : ffmap_i(GEN m, GEN x)
3284 : {
3285 : GEN y;
3286 4452 : long i, lx, tx = typ(x);
3287 4452 : switch(tx)
3288 : {
3289 2534 : case t_FFELT:
3290 2534 : return ffeltmap_i(m, x);
3291 1267 : case t_POL: case t_RFRAC: case t_SER:
3292 : case t_VEC: case t_COL: case t_MAT:
3293 1267 : y = cgetg_copy(x, &lx);
3294 1988 : for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
3295 4564 : for (i=lontyp[tx]; i<lx; i++)
3296 : {
3297 3339 : GEN yi = ffmap_i(m, gel(x,i));
3298 3297 : if (!yi) return NULL;
3299 3297 : gel(y,i) = yi;
3300 : }
3301 1225 : return y;
3302 : }
3303 651 : return gcopy(x);
3304 : }
3305 :
3306 : GEN
3307 1029 : ffmap(GEN m, GEN x)
3308 : {
3309 1029 : pari_sp ltop = avma;
3310 : GEN y;
3311 1029 : checkmap(m, "ffmap");
3312 1029 : y = ffmap_i(m, x);
3313 1029 : if (y) return y;
3314 42 : set_avma(ltop); return cgetg(1,t_VEC);
3315 : }
3316 :
3317 : static GEN
3318 252 : ffeltmaprel_i(GEN m, GEN x)
3319 : {
3320 252 : GEN g = gel(m,1), r = gel(m,2);
3321 252 : if (!FF_samefield(x, g))
3322 0 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3323 252 : if (typ(r)==t_FFELT)
3324 84 : retmkpolmod(FF_map(r, x), pol_x(FF_var(g)));
3325 : else
3326 168 : retmkpolmod(FFX_preimagerel(x, r, ffpartmapimage("ffmap", r)), gcopy(r));
3327 : }
3328 :
3329 : static GEN
3330 252 : ffmaprel_i(GEN m, GEN x)
3331 : {
3332 : GEN y;
3333 252 : long i, lx, tx = typ(x);
3334 252 : switch(tx)
3335 : {
3336 252 : case t_FFELT:
3337 252 : return ffeltmaprel_i(m, x);
3338 0 : case t_POL: case t_RFRAC: case t_SER:
3339 : case t_VEC: case t_COL: case t_MAT:
3340 0 : y = cgetg_copy(x, &lx);
3341 0 : for (i=1; i<lontyp[tx]; i++) y[i] = x[1];
3342 0 : for (i=lontyp[tx]; i<lx; i++)
3343 0 : gel(y,i) = ffmaprel_i(m, gel(x,i));
3344 0 : return y;
3345 : }
3346 0 : return gcopy(x);
3347 : }
3348 :
3349 : GEN
3350 252 : ffmaprel(GEN m, GEN x)
3351 : {
3352 252 : checkmap(m, "ffmaprel");
3353 252 : return ffmaprel_i(m, x);
3354 : }
3355 :
3356 : static void
3357 84 : err_compo(GEN m, GEN n)
3358 84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
3359 :
3360 : GEN
3361 420 : ffcompomap(GEN m, GEN n)
3362 : {
3363 420 : pari_sp av = avma;
3364 420 : GEN g = gel(n,1), r, m2, n2;
3365 420 : checkmap(m, "ffcompomap");
3366 420 : checkmap(n, "ffcompomap");
3367 420 : m2 = gel(m,2); n2 = gel(n,2);
3368 420 : switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
3369 : {
3370 84 : case 0:
3371 84 : if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
3372 42 : r = FF_map(gel(m,2), n2);
3373 42 : break;
3374 84 : case 2:
3375 84 : r = ffmap_i(m, n2);
3376 42 : if (lg(r) == 1) err_compo(m,n);
3377 42 : break;
3378 168 : case 1:
3379 168 : r = ffeltmap_i(m, n2);
3380 126 : if (!r)
3381 : {
3382 : GEN a, A, R, M;
3383 : long dm, dn;
3384 42 : a = ffpartmapimage("ffcompomap",m2);
3385 42 : A = FF_to_FpXQ_i(FF_neg(n2));
3386 42 : setvarn(A, 1);
3387 42 : R = deg1pol(gen_1, A, 0);
3388 42 : setvarn(R, 0);
3389 42 : M = gcopy(m2);
3390 42 : setvarn(M, 1);
3391 42 : r = polresultant0(R, M, 1, 0);
3392 42 : dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
3393 42 : if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
3394 42 : setvarn(r, varn(FF_mod(g)));
3395 : }
3396 126 : break;
3397 84 : case 3:
3398 : {
3399 : GEN M, R, T, p, a;
3400 84 : a = ffpartmapimage("ffcompomap",n2);
3401 84 : if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
3402 42 : p = FF_p_i(gel(n,1));
3403 42 : T = FF_mod(gel(n,1));
3404 42 : setvarn(T, 1);
3405 42 : R = RgX_to_FpXQX(n2,T,p);
3406 42 : setvarn(R, 0);
3407 42 : M = gcopy(m2);
3408 42 : setvarn(M, 1);
3409 42 : r = polresultant0(R, M, 1, 0);
3410 42 : setvarn(r, varn(n2));
3411 : }
3412 : }
3413 252 : return gerepilecopy(av, mkvec2(g,r));
3414 : }
|
