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 1225 : char_update_prime(struct charact *S, GEN p)
37 : {
38 1225 : if (!S->isprime) { S->isprime = 1; S->q = p; }
39 1225 : if (!equalii(p, S->q)) pari_err_MODULUS("characteristic", S->q, p);
40 1218 : }
41 : static void
42 5873 : char_update_int(struct charact *S, GEN n)
43 : {
44 5873 : if (S->isprime)
45 : {
46 7 : if (dvdii(n, S->q)) return;
47 7 : pari_err_MODULUS("characteristic", S->q, n);
48 : }
49 5866 : S->q = gcdii(S->q, n);
50 : }
51 : static void
52 175105 : charact(struct charact *S, GEN x)
53 : {
54 175105 : const long tx = typ(x);
55 : long i, l;
56 175105 : switch(tx)
57 : {
58 4396 : case t_INTMOD:char_update_int(S, gel(x,1)); break;
59 1134 : case t_FFELT: char_update_prime(S, gel(x,4)); break;
60 23765 : 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 23765 : l = lg(x);
64 170849 : for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
65 23751 : break;
66 7 : case t_LIST:
67 7 : x = list_data(x);
68 7 : if (x) charact(S, x);
69 7 : break;
70 : }
71 175077 : }
72 : static void
73 4137 : charact_res(struct charact *S, GEN x)
74 : {
75 4137 : const long tx = typ(x);
76 : long i, l;
77 4137 : switch(tx)
78 : {
79 1477 : 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, padic_p(x)); break;
82 1029 : 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 1029 : l = lg(x);
86 4830 : for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
87 1029 : 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 4137 : }
94 : GEN
95 28007 : characteristic(GEN x)
96 : {
97 : struct charact S;
98 28007 : S.q = gen_0; S.isprime = 0;
99 28007 : charact(&S, x); return S.q;
100 : }
101 : GEN
102 336 : residual_characteristic(GEN x)
103 : {
104 : struct charact S;
105 336 : S.q = gen_0; S.isprime = 0;
106 336 : charact_res(&S, x); return S.q;
107 : }
108 :
109 : int
110 73378310 : Rg_is_Fp(GEN x, GEN *pp)
111 : {
112 : GEN mod;
113 73378310 : switch(typ(x))
114 : {
115 2569539 : case t_INTMOD:
116 2569539 : mod = gel(x,1);
117 2569539 : if (!*pp) *pp = mod;
118 2411927 : 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 2569539 : return 1;
124 59500242 : case t_INT:
125 59500242 : return 1;
126 11308529 : default: return 0;
127 : }
128 : }
129 :
130 : int
131 28216249 : RgX_is_FpX(GEN x, GEN *pp)
132 : {
133 28216249 : long i, lx = lg(x);
134 90127100 : for (i=2; i<lx; i++)
135 73219383 : if (!Rg_is_Fp(gel(x, i), pp))
136 11308524 : return 0;
137 16907717 : 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 193571 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
160 : {
161 : GEN pol, mod, p;
162 193571 : switch(typ(x))
163 : {
164 158900 : case t_INTMOD:
165 158900 : return Rg_is_Fp(x, pp);
166 8561 : case t_INT:
167 8561 : 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 4585 : case t_POLMOD:
181 4585 : mod = gel(x,1); pol = gel(x, 2);
182 4585 : if (!RgX_is_FpX(mod, pp)) return 0;
183 4585 : if (typ(pol)==t_POL)
184 : {
185 4557 : if (!RgX_is_FpX(pol, pp)) return 0;
186 : }
187 28 : else if (!Rg_is_Fp(pol, pp)) return 0;
188 4585 : 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 4585 : return 1;
195 :
196 154 : default: return 0;
197 : }
198 : }
199 :
200 : int
201 2021370 : RgX_is_ZXX(GEN x, long *v)
202 : {
203 : long i;
204 8485379 : for (i = lg(x)-1; i > 1; i--)
205 : {
206 6464219 : GEN xi = gel(x,i);
207 6464219 : long t = typ(xi), vi;
208 6464219 : if (t==t_INT) continue;
209 1012939 : if (t!=t_POL || !RgX_is_ZX(xi)) return 0;
210 1012785 : vi = varn(xi);
211 1012785 : if (*v < 0) *v = vi;
212 2107 : else if (vi!=*v) return 0;
213 : }
214 2021160 : return 1;
215 : }
216 :
217 : int
218 22974 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
219 : {
220 22974 : long i, lx = lg(x);
221 215789 : for (i = 2; i < lx; i++)
222 192913 : if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
223 22876 : return 1;
224 : }
225 :
226 : /************************************************************************
227 : ** **
228 : ** Ring conversion **
229 : ** **
230 : ************************************************************************/
231 :
232 : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
233 : * If x is an INTMOD, assume modulus is a multiple of p. */
234 : GEN
235 52389296 : Rg_to_Fp(GEN x, GEN p)
236 : {
237 52389296 : if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
238 4535758 : switch(typ(x))
239 : {
240 288513 : case t_INT: return modii(x, p);
241 18790 : case t_FRAC: {
242 18790 : pari_sp av = avma;
243 18790 : GEN z = modii(gel(x,1), p);
244 18790 : if (z == gen_0) return gen_0;
245 18785 : return gc_INT(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
246 : }
247 70 : case t_PADIC: return padic_to_Fp(x, p);
248 4228397 : case t_INTMOD: {
249 4228397 : GEN q = gel(x,1), a = gel(x,2);
250 4228397 : if (equalii(q, p)) return icopy(a);
251 14 : if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
252 0 : return remii(a, p);
253 : }
254 0 : default: pari_err_TYPE("Rg_to_Fp",x);
255 : return NULL; /* LCOV_EXCL_LINE */
256 : }
257 : }
258 : /* If x is a POLMOD, assume modulus is a multiple of T. */
259 : GEN
260 1288276 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
261 : {
262 1288276 : long ta, tx = typ(x), v = get_FpX_var(T);
263 : GEN a, b;
264 1288276 : if (is_const_t(tx))
265 : {
266 57432 : if (tx == t_FFELT)
267 : {
268 17355 : GEN z = FF_to_FpXQ(x);
269 17355 : setvarn(z, v);
270 17355 : return z;
271 : }
272 40077 : return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
273 : }
274 1230844 : switch(tx)
275 : {
276 1229122 : case t_POLMOD:
277 1229122 : b = gel(x,1);
278 1229122 : a = gel(x,2); ta = typ(a);
279 1229122 : if (is_const_t(ta))
280 5726 : return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
281 1223396 : b = RgX_to_FpX(b, p); if (varn(b) != v) break;
282 1223396 : a = RgX_to_FpX(a, p);
283 1223396 : if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
284 1223396 : return FpX_rem(a, T, p);
285 0 : break;
286 1722 : case t_POL:
287 1722 : if (varn(x) != v) break;
288 1715 : return FpX_rem(RgX_to_FpX(x,p), T, p);
289 0 : case t_RFRAC:
290 0 : a = Rg_to_FpXQ(gel(x,1), T,p);
291 0 : b = Rg_to_FpXQ(gel(x,2), T,p);
292 0 : return FpXQ_div(a,b, T,p);
293 : }
294 7 : pari_err_TYPE("Rg_to_FpXQ",x);
295 : return NULL; /* LCOV_EXCL_LINE */
296 : }
297 : GEN
298 3334075 : RgX_to_FpX(GEN x, GEN p)
299 : {
300 : long i, l;
301 3334075 : GEN z = cgetg_copy(x, &l); z[1] = x[1];
302 14813903 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
303 3334075 : return FpX_renormalize(z, l);
304 : }
305 :
306 : GEN
307 140 : RgV_to_FpV(GEN x, GEN p)
308 483 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
309 :
310 : GEN
311 1751090 : RgC_to_FpC(GEN x, GEN p)
312 28499485 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
313 :
314 : GEN
315 222349 : RgM_to_FpM(GEN x, GEN p)
316 1973397 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
317 :
318 : GEN
319 369001 : RgV_to_Flv(GEN x, ulong p)
320 1676894 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
321 :
322 : GEN
323 118296 : RgM_to_Flm(GEN x, ulong p)
324 422998 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
325 :
326 : GEN
327 3887 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
328 : {
329 3887 : long i, l = lg(x);
330 3887 : GEN z = cgetg(l, t_POL); z[1] = x[1];
331 38466 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
332 3887 : return FpXQX_renormalize(z, l);
333 : }
334 : GEN
335 49200 : RgX_to_FqX(GEN x, GEN T, GEN p)
336 : {
337 49200 : long i, l = lg(x);
338 49200 : GEN z = cgetg(l, t_POL); z[1] = x[1];
339 49200 : if (T)
340 10990 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
341 : else
342 791436 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
343 49200 : return FpXQX_renormalize(z, l);
344 : }
345 :
346 : GEN
347 219128 : RgC_to_FqC(GEN x, GEN T, GEN p)
348 : {
349 219128 : long i, l = lg(x);
350 219128 : GEN z = cgetg(l, t_COL);
351 219128 : if (T)
352 1430310 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
353 : else
354 0 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
355 219128 : return z;
356 : }
357 :
358 : GEN
359 52325 : RgM_to_FqM(GEN x, GEN T, GEN p)
360 271418 : { pari_APPLY_same(RgC_to_FqC(gel(x, i), T, p)) }
361 :
362 : /* lg(V) > 1 */
363 : GEN
364 851487 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
365 : {
366 851487 : pari_sp av = avma;
367 851487 : long i, l = lg(V);
368 851487 : GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
369 4201029 : for(i=2; i<l; i++)
370 : {
371 3349542 : z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
372 3349542 : if ((i & 7) == 0) z = gc_upto(av, z);
373 : }
374 851487 : return gc_upto(av, FpX_red(z,p));
375 : }
376 :
377 : GEN
378 55832 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
379 : {
380 55832 : long i, lz = lg(y);
381 : GEN z;
382 55832 : if (!T) return FpX_Fp_add(y, x, p);
383 8666 : if (lz == 2) return scalarpol(x, varn(y));
384 7952 : z = cgetg(lz,t_POL); z[1] = y[1];
385 7952 : gel(z,2) = Fq_add(gel(y,2),x, T, p);
386 7952 : if (lz == 3) z = FpXX_renormalize(z,lz);
387 : else
388 33145 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
389 7952 : return z;
390 : }
391 :
392 : GEN
393 1059 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
394 : {
395 1059 : long i, lz = lg(y);
396 : GEN z;
397 1059 : if (!T) return FpX_Fp_sub(y, x, p);
398 1059 : if (lz == 2) return scalarpol(x, varn(y));
399 1059 : z = cgetg(lz,t_POL); z[1] = y[1];
400 1059 : gel(z,2) = Fq_sub(gel(y,2), x, T, p);
401 1059 : if (lz == 3) z = FpXX_renormalize(z,lz);
402 : else
403 10278 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
404 1059 : return z;
405 : }
406 :
407 : GEN
408 149052 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
409 : {
410 : long i, lP;
411 149052 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
412 918799 : for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
413 149052 : gel(res,lP-1) = gen_1; return res;
414 : }
415 :
416 : GEN
417 38189 : FpXQX_normalize(GEN z, GEN T, GEN p)
418 : {
419 : GEN lc;
420 38189 : if (lg(z) == 2) return z;
421 38175 : lc = leading_coeff(z);
422 38175 : if (typ(lc) == t_POL)
423 : {
424 2167 : if (lg(lc) > 3) /* nonconstant */
425 1902 : return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
426 : /* constant */
427 265 : lc = gel(lc,2);
428 265 : z = shallowcopy(z);
429 265 : gel(z, lg(z)-1) = lc;
430 : }
431 : /* lc a t_INT */
432 36273 : if (equali1(lc)) return z;
433 87 : return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
434 : }
435 :
436 : GEN
437 391145 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
438 : {
439 : pari_sp av;
440 : GEN p1, r;
441 391145 : long j, i=lg(x)-1;
442 391145 : if (i<=2)
443 45107 : return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
444 346038 : av=avma; p1=gel(x,i);
445 : /* specific attention to sparse polynomials (see poleval)*/
446 : /*You've guessed it! It's a copy-paste(tm)*/
447 1151145 : for (i--; i>=2; i=j-1)
448 : {
449 805794 : for (j=i; !signe(gel(x,j)); j--)
450 686 : if (j==2)
451 : {
452 490 : if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
453 490 : return gc_upto(av, Fq_mul(p1,y, T, p));
454 : }
455 805108 : r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
456 805108 : p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
457 : }
458 345547 : return gc_upto(av, p1);
459 : }
460 :
461 : GEN
462 97391 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
463 : {
464 97391 : long i, lb = lg(Q);
465 : GEN z;
466 97391 : if (!T) return FpXY_evalx(Q, x, p);
467 87031 : z = cgetg(lb, t_POL); z[1] = Q[1];
468 463295 : for (i=2; i<lb; i++)
469 : {
470 376264 : GEN q = gel(Q,i);
471 376264 : gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
472 : }
473 87031 : return FpXQX_renormalize(z, lb);
474 : }
475 :
476 : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
477 : GEN
478 14623 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
479 : {
480 14623 : pari_sp av = avma;
481 14623 : if (!T) return FpXY_eval(Q, y, x, p);
482 966 : return gc_upto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
483 : }
484 :
485 : /* a X^d */
486 : GEN
487 12542309 : monomial(GEN a, long d, long v)
488 : {
489 : long i, n;
490 : GEN P;
491 12542309 : if (d < 0) {
492 14 : if (isrationalzero(a)) return pol_0(v);
493 14 : retmkrfrac(a, pol_xn(-d, v));
494 : }
495 12542295 : if (gequal0(a))
496 : {
497 89369 : if (isexactzero(a)) return scalarpol_shallow(a,v);
498 0 : n = d+2; P = cgetg(n+1, t_POL);
499 0 : P[1] = evalsigne(0) | evalvarn(v);
500 : }
501 : else
502 : {
503 12452925 : n = d+2; P = cgetg(n+1, t_POL);
504 12452925 : P[1] = evalsigne(1) | evalvarn(v);
505 : }
506 31853419 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
507 12452925 : gel(P,i) = a; return P;
508 : }
509 : GEN
510 156275 : monomialcopy(GEN a, long d, long v)
511 : {
512 : long i, n;
513 : GEN P;
514 156275 : if (d < 0) {
515 14 : if (isrationalzero(a)) return pol_0(v);
516 14 : retmkrfrac(gcopy(a), pol_xn(-d, v));
517 : }
518 156261 : if (gequal0(a))
519 : {
520 14 : if (isexactzero(a)) return scalarpol(a,v);
521 7 : n = d+2; P = cgetg(n+1, t_POL);
522 7 : P[1] = evalsigne(0) | evalvarn(v);
523 : }
524 : else
525 : {
526 156247 : n = d+2; P = cgetg(n+1, t_POL);
527 156247 : P[1] = evalsigne(1) | evalvarn(v);
528 : }
529 313103 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
530 156254 : gel(P,i) = gcopy(a); return P;
531 : }
532 : GEN
533 325893 : pol_x_powers(long N, long v)
534 : {
535 325893 : GEN L = cgetg(N+1,t_VEC);
536 : long i;
537 788813 : for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
538 325895 : return L;
539 : }
540 :
541 : GEN
542 0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
543 : {
544 0 : return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
545 : }
546 :
547 : GEN
548 0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
549 : {
550 0 : return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
551 : }
552 :
553 : /*******************************************************************/
554 : /* */
555 : /* Fq */
556 : /* */
557 : /*******************************************************************/
558 :
559 : GEN
560 11593382 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
561 : {
562 : (void)T;
563 11593382 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
564 : {
565 1143736 : case 0: return Fp_add(x,y,p);
566 748745 : case 1: return FpX_Fp_add(x,y,p);
567 95461 : case 2: return FpX_Fp_add(y,x,p);
568 9605440 : case 3: return FpX_add(x,y,p);
569 : }
570 : return NULL;/*LCOV_EXCL_LINE*/
571 : }
572 :
573 : GEN
574 8563005 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
575 : {
576 : (void)T;
577 8563005 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
578 : {
579 256936 : case 0: return Fp_sub(x,y,p);
580 24561 : case 1: return FpX_Fp_sub(x,y,p);
581 23924 : case 2: return Fp_FpX_sub(x,y,p);
582 8257584 : case 3: return FpX_sub(x,y,p);
583 : }
584 : return NULL;/*LCOV_EXCL_LINE*/
585 : }
586 :
587 : GEN
588 1080336 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
589 : {
590 : (void)T;
591 1080336 : return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
592 : }
593 :
594 : GEN
595 81417 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
596 : {
597 : (void)T;
598 81417 : return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
599 : }
600 :
601 : /* If T==NULL do not reduce*/
602 : GEN
603 8609154 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
604 : {
605 8609154 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
606 : {
607 1038716 : case 0: return Fp_mul(x,y,p);
608 170153 : case 1: return FpX_Fp_mul(x,y,p);
609 432345 : case 2: return FpX_Fp_mul(y,x,p);
610 6967942 : case 3: if (T) return FpXQ_mul(x,y,T,p);
611 4399046 : else return FpX_mul(x,y,p);
612 : }
613 : return NULL;/*LCOV_EXCL_LINE*/
614 : }
615 :
616 : /* If T==NULL do not reduce*/
617 : GEN
618 490240 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
619 : {
620 : (void) T;
621 490240 : return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
622 : }
623 :
624 : /* y t_INT */
625 : GEN
626 43985 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
627 : {
628 : (void)T;
629 6822 : return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
630 50807 : : Fp_mul(x,y,p);
631 : }
632 : /* If T==NULL do not reduce*/
633 : GEN
634 611236 : Fq_sqr(GEN x, GEN T, GEN p)
635 : {
636 611236 : if (typ(x) == t_POL)
637 : {
638 70647 : if (T) return FpXQ_sqr(x,T,p);
639 0 : else return FpX_sqr(x,p);
640 : }
641 : else
642 540589 : return Fp_sqr(x,p);
643 : }
644 :
645 : GEN
646 0 : Fq_neg_inv(GEN x, GEN T, GEN p)
647 : {
648 0 : if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
649 0 : return FpXQ_inv(FpX_neg(x,p),T,p);
650 : }
651 :
652 : GEN
653 0 : Fq_invsafe(GEN x, GEN pol, GEN p)
654 : {
655 0 : if (typ(x) == t_INT) return Fp_invsafe(x,p);
656 0 : return FpXQ_invsafe(x,pol,p);
657 : }
658 :
659 : GEN
660 89298 : Fq_inv(GEN x, GEN pol, GEN p)
661 : {
662 89298 : if (typ(x) == t_INT) return Fp_inv(x,p);
663 81504 : return FpXQ_inv(x,pol,p);
664 : }
665 :
666 : GEN
667 343791 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
668 : {
669 343791 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
670 : {
671 318402 : case 0: return Fp_div(x,y,p);
672 16702 : case 1: return FpX_Fp_div(x,y,p);
673 406 : case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
674 8281 : case 3: return FpXQ_div(x,y,pol,p);
675 : }
676 : return NULL;/*LCOV_EXCL_LINE*/
677 : }
678 :
679 : GEN
680 1088232 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
681 : {
682 1088232 : if (typ(x) == t_INT) return Fp_pow(x,n,p);
683 137466 : return FpXQ_pow(x,n,pol,p);
684 : }
685 :
686 : GEN
687 15050 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
688 : {
689 15050 : if (typ(x) == t_INT) return Fp_powu(x,n,p);
690 1267 : return FpXQ_powu(x,n,pol,p);
691 : }
692 :
693 : GEN
694 1893010 : Fq_sqrt(GEN x, GEN T, GEN p)
695 : {
696 1893010 : if (typ(x) == t_INT)
697 : {
698 1825085 : if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
699 9642 : x = scalarpol_shallow(x, get_FpX_var(T));
700 : }
701 77567 : return FpXQ_sqrt(x,T,p);
702 : }
703 : GEN
704 170786 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
705 : {
706 170786 : if (typ(x) == t_INT)
707 : {
708 : long d;
709 170415 : if (!T) return Fp_sqrtn(x,n,p,zeta);
710 126 : d = get_FpX_degree(T);
711 126 : if (ugcdiu(n,d) == 1)
712 : {
713 56 : if (!zeta) return Fp_sqrtn(x,n,p,NULL);
714 : /* gcd(n,p-1)=gcd(n,q-1): same number of solutions in Fp and F_q */
715 21 : if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
716 14 : return Fp_sqrtn(x,n,p,zeta);
717 : }
718 77 : x = scalarpol(x, get_FpX_var(T)); /* left on stack */
719 : }
720 448 : return FpXQ_sqrtn(x,n,T,p,zeta);
721 : }
722 :
723 : struct _Fq_field
724 : {
725 : GEN T, p;
726 : };
727 :
728 : static GEN
729 303247 : _Fq_red(void *E, GEN x)
730 303247 : { struct _Fq_field *s = (struct _Fq_field *)E;
731 303247 : return Fq_red(x, s->T, s->p);
732 : }
733 :
734 : static GEN
735 362523 : _Fq_add(void *E, GEN x, GEN y)
736 : {
737 : (void) E;
738 362523 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
739 : {
740 14 : case 0: return addii(x,y);
741 0 : case 1: return ZX_Z_add(x,y);
742 15918 : case 2: return ZX_Z_add(y,x);
743 346591 : default: return ZX_add(x,y);
744 : }
745 : }
746 :
747 : static GEN
748 315028 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
749 :
750 : static GEN
751 243341 : _Fq_mul(void *E, GEN x, GEN y)
752 : {
753 : (void) E;
754 243341 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
755 : {
756 679 : case 0: return mulii(x,y);
757 1085 : case 1: return ZX_Z_mul(x,y);
758 1043 : case 2: return ZX_Z_mul(y,x);
759 240534 : default: return ZX_mul(x,y);
760 : }
761 : }
762 :
763 : static GEN
764 9331 : _Fq_inv(void *E, GEN x)
765 9331 : { struct _Fq_field *s = (struct _Fq_field *)E;
766 9331 : return Fq_inv(x,s->T,s->p);
767 : }
768 :
769 : static int
770 38388 : _Fq_equal0(GEN x) { return signe(x)==0; }
771 :
772 : static GEN
773 4151 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
774 :
775 : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
776 : _Fq_inv,_Fq_equal0,_Fq_s};
777 :
778 4725 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
779 : {
780 4725 : if (!T)
781 0 : return get_Fp_field(E, p);
782 : else
783 : {
784 4725 : GEN z = new_chunk(sizeof(struct _Fq_field));
785 4725 : struct _Fq_field *e = (struct _Fq_field *) z;
786 4725 : e->T = T; e->p = p; *E = (void*)e;
787 4725 : return &Fq_field;
788 : }
789 : }
790 :
791 : /*******************************************************************/
792 : /* */
793 : /* Fq[X] */
794 : /* */
795 : /*******************************************************************/
796 : /* FpV_red(vecbinomial(n), p) */
797 : static GEN
798 0 : vecbinomial_Fp(long n, GEN p)
799 : {
800 0 : GEN C = vecbinomial(n) + 1;
801 0 : long k, d = (n + 1) >> 1;
802 0 : for (k = d; k >= 1; k--)
803 : {
804 0 : GEN a = gel(C,k);
805 0 : if (cmpii(a, p) < 0) break;
806 0 : gel(C,k) = gel(C, n-k) = remii(a, p);
807 : }
808 0 : return C - 1;
809 : }
810 : /* (X+u)^n */
811 : static GEN
812 434 : Fp_XpN_powu(GEN u, ulong n, GEN p, long v)
813 : {
814 : pari_sp av;
815 : ulong k;
816 : GEN B, C, V;
817 434 : if (!n) return pol_1(v);
818 434 : if (is_pm1(u))
819 434 : return Xpm1_powu(n, signe(u), v);
820 0 : av = avma;
821 0 : V = Fp_powers(u, n, p);
822 0 : B = vecbinomial_Fp(n, p);
823 0 : C = cgetg(n+3, t_POL);
824 0 : C[1] = evalsigne(1)| evalvarn(v);
825 0 : for (k=1; k <= n+1; k++)
826 0 : gel(C,k+1) = Fp_mul(gel(V,n+2-k), gel(B,k), p);
827 0 : return gc_upto(av, C);
828 : }
829 :
830 : static GEN
831 700 : FpX_Fp_translate_basecase(GEN P, GEN c, GEN p)
832 : {
833 700 : pari_sp av = avma;
834 : GEN Q;
835 : long i, k, n;
836 :
837 700 : if (!signe(P) || !signe(c)) return ZX_copy(P);
838 560 : Q = leafcopy(P); n = degpol(P);
839 1316 : for (i=1; i<=n; i++)
840 : {
841 2016 : for (k=n-i; k<n; k++)
842 1260 : gel(Q,2+k) = Fp_add(gel(Q,2+k), Fp_mul(c, gel(Q,2+k+1), p), p);
843 :
844 756 : if (gc_needed(av,2))
845 : {
846 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FpX_Fp_translate, i = %ld/%ld", i,n);
847 0 : Q = gc_GEN(av, Q);
848 : }
849 : }
850 560 : return gc_GEN(av, FpX_renormalize(Q, lg(Q)));
851 : }
852 :
853 : GEN
854 1134 : FpX_Fp_translate(GEN P, GEN c, GEN p)
855 : {
856 1134 : pari_sp av = avma;
857 1134 : long n = degpol(P);
858 1134 : if (n<=3 || 25*(n-3) < expi(p))
859 700 : return FpX_Fp_translate_basecase(P, c, p);
860 : else
861 : {
862 434 : long d = n >> 1;
863 434 : GEN Q = FpX_Fp_translate(RgX_shift_shallow(P, -d), c, p);
864 434 : GEN R = FpX_Fp_translate(RgXn_red_shallow(P, d), c, p);
865 434 : GEN S = Fp_XpN_powu(c, d, p, varn(P));
866 434 : return gc_upto(av, FpX_add(FpX_mul(Q, S, p), R, p));
867 : }
868 : }
869 :
870 : /* (X+u)^n mod (T,p) */
871 : static GEN
872 0 : FpXQX_XpN_powu(GEN u, ulong n, GEN T, GEN p, long v)
873 : {
874 : pari_sp av;
875 : ulong k;
876 : GEN B, C, V;
877 0 : if (!n) return pol_1(v);
878 0 : av = avma;
879 0 : V = FpXQ_powers(u, n, T, p);
880 0 : B = vecbinomial_Fp(n, p);
881 0 : C = cgetg(n+3, t_POL);
882 0 : C[1] = evalsigne(1)| evalvarn(v);
883 0 : for (k=1; k <= n+1; k++)
884 0 : gel(C,k+1) = Fq_mul(gel(V,n+2-k), gel(B,k), T, p);
885 0 : return gc_upto(av, C);
886 : }
887 :
888 : static GEN
889 33887 : FpXQX_FpXQ_translate_basecase(GEN P, GEN c, GEN T, GEN p)
890 : {
891 33887 : pari_sp av = avma;
892 : GEN Q;
893 : long i, k, n;
894 :
895 : /* signe works for t_(INT|POL) */
896 33887 : if (!signe(P) || !signe(c)) return RgX_copy(P);
897 33887 : Q = leafcopy(P); n = degpol(P);
898 150080 : for (i=1; i<=n; i++)
899 : {
900 433594 : for (k=n-i; k<n; k++)
901 317401 : gel(Q,2+k) = Fq_add(gel(Q,2+k), Fq_mul(c, gel(Q,2+k+1), T, p), T, p);
902 :
903 116193 : if (gc_needed(av,2))
904 : {
905 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FqX_Fq_translate, i = %ld/%ld", i,n);
906 0 : Q = gc_GEN(av, Q);
907 : }
908 : }
909 33887 : return gc_GEN(av, FpXQX_renormalize(Q, lg(Q)));
910 : }
911 :
912 : GEN
913 33887 : FpXQX_FpXQ_translate(GEN P, GEN c, GEN T, GEN p)
914 : {
915 33887 : pari_sp av = avma;
916 33887 : long n = degpol(P);
917 33887 : if (n < 220)
918 33887 : return FpXQX_FpXQ_translate_basecase(P, c, T, p);
919 : else
920 : {
921 0 : long d = n >> 1;
922 0 : GEN Q = FpXQX_FpXQ_translate(RgX_shift_shallow(P, -d), c, T, p);
923 0 : GEN R = FpXQX_FpXQ_translate(RgXn_red_shallow(P, d), c, T, p);
924 0 : GEN S = FpXQX_XpN_powu(c, d, T, p, varn(P));
925 0 : return gc_upto(av, FpXX_add(FpXQX_mul(Q, S, T, p), R, p));
926 : }
927 : }
928 :
929 : GEN
930 33880 : FqX_Fq_translate(GEN x, GEN y, GEN T, GEN p)
931 : {
932 33880 : if (!T) return FpX_Fp_translate(x,y,p);
933 33880 : if (typ(y)==t_INT)
934 : {
935 0 : pari_sp av = avma;
936 0 : y = scalarpol_shallow(y,varn(T));
937 0 : return gc_upto(av, FpXQX_FpXQ_translate(x,y,T,p));
938 : }
939 33880 : return FpXQX_FpXQ_translate(x,y,T,p);
940 : }
941 :
942 : GEN
943 40456 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
944 : {
945 40456 : pari_sp ltop = avma;
946 : long k;
947 : GEN W;
948 40456 : if (lgefint(p) == 3)
949 : {
950 31743 : ulong pp = p[2];
951 31743 : GEN Tl = ZX_to_Flx(T, pp);
952 31743 : GEN Vl = FqC_to_FlxqC(V, Tl, pp);
953 31743 : Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
954 31741 : return gc_upto(ltop, FlxX_to_ZXX(Tl));
955 : }
956 8713 : W = cgetg(lg(V),t_VEC);
957 78138 : for(k=1; k < lg(V); k++)
958 69424 : gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
959 8714 : return gc_upto(ltop, FpXQXV_prod(W, T, p));
960 : }
961 :
962 : GEN
963 109464 : FqV_red(GEN x, GEN T, GEN p)
964 778530 : { pari_APPLY_type(t_VEC, Fq_red(gel(x,i), T, p)) }
965 :
966 : GEN
967 24058 : FqC_red(GEN x, GEN T, GEN p)
968 163384 : { pari_APPLY_type(t_COL, Fq_red(gel(x,i), T, p)) }
969 :
970 : GEN
971 1701 : FqM_red(GEN x, GEN T, GEN p)
972 5411 : { pari_APPLY_same(FqC_red(gel(x,i), T, p)) }
973 :
974 : GEN
975 0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
976 : {
977 0 : if (!T) return FpC_add(x, y, p);
978 0 : pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
979 : }
980 :
981 : GEN
982 0 : FqM_add(GEN x, GEN y, GEN T, GEN p)
983 : {
984 0 : if (!T) return FpM_add(x, y, p);
985 0 : pari_APPLY_same(FqC_add(gel(x,i), gel(y,i), T, p))
986 : }
987 :
988 : GEN
989 0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
990 : {
991 0 : if (!T) return FpC_sub(x, y, p);
992 0 : pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
993 : }
994 :
995 : GEN
996 0 : FqM_sub(GEN x, GEN y, GEN T, GEN p)
997 : {
998 0 : if (!T) return FpM_sub(x, y, p);
999 0 : pari_APPLY_same(FqC_sub(gel(x,i), gel(y,i), T, p))
1000 : }
1001 :
1002 : GEN
1003 0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
1004 : {
1005 0 : if (!T) return FpC_Fp_mul(x, y, p);
1006 0 : pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
1007 : }
1008 :
1009 : GEN
1010 105 : FqC_FqV_mul(GEN x, GEN y, GEN T, GEN p)
1011 : {
1012 105 : long i,j, lx=lg(x), ly=lg(y);
1013 : GEN z;
1014 105 : if (ly==1) return cgetg(1,t_MAT);
1015 105 : z = cgetg(ly,t_MAT);
1016 819 : for (j=1; j < ly; j++)
1017 : {
1018 714 : GEN zj = cgetg(lx,t_COL);
1019 4200 : for (i=1; i<lx; i++) gel(zj,i) = Fq_mul(gel(x,i),gel(y,j), T, p);
1020 714 : gel(z, j) = zj;
1021 : }
1022 105 : return z;
1023 : }
1024 :
1025 : GEN
1026 5467 : FpXC_center(GEN x, GEN p, GEN pov2)
1027 41476 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
1028 :
1029 : GEN
1030 109020 : FqC_to_FlxqC(GEN x, GEN T, ulong p)
1031 109020 : { long sv = get_Flx_var(T);
1032 4834743 : pari_APPLY_type(t_COL,typ(gel(x,i))==t_INT ?
1033 : Z_to_Flx(gel(x,i), p, sv): ZX_to_Flx(gel(x,i), p)) }
1034 :
1035 : GEN
1036 8708 : FqM_to_FlxqM(GEN x, GEN T, ulong p)
1037 85985 : { pari_APPLY_same(FqC_to_FlxqC(gel(x,i), T, p)) }
1038 :
1039 : GEN
1040 1800 : FpXM_center(GEN x, GEN p, GEN pov2)
1041 7267 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
1042 :
1043 : /*******************************************************************/
1044 : /* */
1045 : /* GENERIC CRT */
1046 : /* */
1047 : /*******************************************************************/
1048 : static GEN
1049 9376931 : primelist(forprime_t *S, long n, GEN dB)
1050 : {
1051 9376931 : GEN P = cgetg(n+1, t_VECSMALL);
1052 9376914 : long i = 1;
1053 : ulong p;
1054 22290945 : while (i <= n && (p = u_forprime_next(S)))
1055 12914031 : if (!dB || umodiu(dB, p)) P[i++] = p;
1056 9376910 : return P;
1057 : }
1058 :
1059 : void
1060 8788835 : gen_inccrt_i(const char *str, GEN worker, GEN dB, long n, long mmin,
1061 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
1062 : GEN center(GEN, GEN, GEN))
1063 : {
1064 8788835 : long m = mmin? minss(mmin, n): usqrt(n);
1065 : GEN H, P, mod;
1066 : pari_timer ti;
1067 8788835 : if (DEBUGLEVEL > 4)
1068 : {
1069 0 : timer_start(&ti);
1070 0 : err_printf("%s: nb primes: %ld\n",str, n);
1071 : }
1072 8788824 : if (m == 1)
1073 : {
1074 8473150 : GEN P = primelist(S, n, dB);
1075 8473124 : GEN done = closure_callgen1(worker, P);
1076 8473093 : H = gel(done,1);
1077 8473093 : mod = gel(done,2);
1078 8473093 : if (!*pH && center) H = center(H, mod, shifti(mod,-1));
1079 8473053 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
1080 : }
1081 : else
1082 : {
1083 315674 : long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
1084 : struct pari_mt pt;
1085 315674 : long pending = 0;
1086 315674 : H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
1087 315674 : mt_queue_start_lim(&pt, worker, m);
1088 1285470 : for (i=1; i<=m || pending; i++)
1089 : {
1090 : GEN done;
1091 969796 : GEN pr = i <= m ? mkvec(primelist(S, i<=r ? s: s-1, dB)): NULL;
1092 969797 : mt_queue_submit(&pt, i, pr);
1093 969795 : done = mt_queue_get(&pt, NULL, &pending);
1094 969795 : if (done)
1095 : {
1096 903784 : di++;
1097 903784 : gel(H, di) = gel(done,1);
1098 903784 : gel(P, di) = gel(done,2);
1099 903784 : if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
1100 : }
1101 : }
1102 315674 : mt_queue_end(&pt);
1103 315674 : if (DEBUGLEVEL>5) err_printf("\n");
1104 315674 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
1105 315674 : H = crt(H, P, &mod);
1106 315674 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
1107 : }
1108 8788727 : if (*pH) H = crt(mkvec2(*pH, H), mkvec2(*pmod, mod), &mod);
1109 8788727 : *pH = H; *pmod = mod;
1110 8788727 : }
1111 : void
1112 3075715 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
1113 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
1114 : GEN center(GEN, GEN, GEN))
1115 : {
1116 3075715 : pari_sp av = avma;
1117 3075715 : gen_inccrt_i(str, worker, dB, n, mmin, S, pH, pmod, crt, center);
1118 3075657 : (void)gc_all(av, 2, pH, pmod);
1119 3075756 : }
1120 :
1121 : GEN
1122 2288080 : gen_crt(const char *str, GEN worker, forprime_t *S, GEN dB, ulong bound, long mmin, GEN *pmod,
1123 : GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
1124 : {
1125 2288080 : GEN mod = gen_1, H = NULL;
1126 : ulong e;
1127 :
1128 2288080 : bound++;
1129 4576202 : while (bound > (e = expi(mod)))
1130 : {
1131 2288038 : long n = (bound - e) / expu(S->p) + 1;
1132 2288073 : gen_inccrt(str, worker, dB, n, mmin, S, &H, &mod, crt, center);
1133 : }
1134 2288108 : if (pmod) *pmod = mod;
1135 2288108 : return H;
1136 : }
1137 :
1138 : /*******************************************************************/
1139 : /* */
1140 : /* MODULAR GCD */
1141 : /* */
1142 : /*******************************************************************/
1143 : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
1144 : static GEN
1145 5163125 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
1146 : {
1147 5163125 : ulong d, amod = umodiu(a, p);
1148 5163127 : pari_sp av = avma;
1149 : GEN ax;
1150 :
1151 5163127 : if (b == amod) return NULL;
1152 2126500 : d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
1153 2126921 : if (d >= 1 + (p>>1))
1154 1037807 : ax = subii(a, mului(p-d, q));
1155 : else
1156 : {
1157 1089114 : ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
1158 1088915 : if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
1159 : }
1160 2126502 : return gc_INT(av, ax);
1161 : }
1162 : GEN
1163 357 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
1164 : GEN
1165 32466 : ZX_init_CRT(GEN Hp, ulong p, long v)
1166 : {
1167 32466 : long i, l = lg(Hp), lim = (long)(p>>1);
1168 32466 : GEN H = cgetg(l, t_POL);
1169 32466 : H[1] = evalsigne(1) | evalvarn(v);
1170 801773 : for (i=2; i<l; i++)
1171 769307 : gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
1172 32466 : return ZX_renormalize(H,l);
1173 : }
1174 :
1175 : GEN
1176 5985 : ZM_init_CRT(GEN Hp, ulong p)
1177 : {
1178 5985 : long i,j, m, l = lg(Hp), lim = (long)(p>>1);
1179 5985 : GEN c, cp, H = cgetg(l, t_MAT);
1180 5985 : if (l==1) return H;
1181 5985 : m = lgcols(Hp);
1182 20244 : for (j=1; j<l; j++)
1183 : {
1184 14259 : cp = gel(Hp,j);
1185 14259 : c = cgetg(m, t_COL);
1186 14259 : gel(H,j) = c;
1187 169400 : for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
1188 : }
1189 5985 : return H;
1190 : }
1191 :
1192 : int
1193 7525 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
1194 : {
1195 7525 : GEN h, q = *ptq, qp = muliu(q,p);
1196 7525 : ulong qinv = Fl_inv(umodiu(q,p), p);
1197 7525 : int stable = 1;
1198 7525 : h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
1199 7525 : if (h) { *H = h; stable = 0; }
1200 7525 : *ptq = qp; return stable;
1201 : }
1202 :
1203 : static int
1204 148356 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
1205 : {
1206 148356 : GEN H = *ptH, h, qp2 = shifti(qp,-1);
1207 148354 : ulong qinv = Fl_inv(umodiu(q,p), p);
1208 148355 : long i, l = lg(H), lp = lg(Hp);
1209 148355 : int stable = 1;
1210 :
1211 148355 : if (l < lp)
1212 : { /* degree increases */
1213 0 : GEN x = cgetg(lp, t_POL);
1214 0 : for (i=1; i<l; i++) x[i] = H[i];
1215 0 : for ( ; i<lp; i++) gel(x,i) = gen_0;
1216 0 : *ptH = H = x;
1217 0 : stable = 0;
1218 148355 : } else if (l > lp)
1219 : { /* degree decreases */
1220 0 : GEN x = cgetg(l, t_VECSMALL);
1221 0 : for (i=1; i<lp; i++) x[i] = Hp[i];
1222 0 : for ( ; i<l; i++) x[i] = 0;
1223 0 : Hp = x; lp = l;
1224 : }
1225 4940231 : for (i=2; i<lp; i++)
1226 : {
1227 4791963 : h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
1228 4791876 : if (h) { gel(H,i) = h; stable = 0; }
1229 : }
1230 148268 : (void)ZX_renormalize(H,lp);
1231 148356 : return stable;
1232 : }
1233 :
1234 : int
1235 0 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
1236 : {
1237 0 : GEN q = *ptq, qp = muliu(q,p);
1238 0 : int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
1239 0 : *ptq = qp; return stable;
1240 : }
1241 :
1242 : int
1243 7787 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1244 : {
1245 7787 : GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1246 7787 : ulong qinv = Fl_inv(umodiu(q,p), p);
1247 7787 : long i,j, l = lg(H), m = lgcols(H);
1248 7787 : int stable = 1;
1249 27451 : for (j=1; j<l; j++)
1250 205516 : for (i=1; i<m; i++)
1251 : {
1252 185852 : h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
1253 185852 : if (h) { gcoeff(H,i,j) = h; stable = 0; }
1254 : }
1255 7787 : *ptq = qp; return stable;
1256 : }
1257 :
1258 : GEN
1259 679 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
1260 : {
1261 : long i, j, k;
1262 : GEN H;
1263 679 : long m, l = lg(Hp), lim = (long)(p>>1), n;
1264 679 : H = cgetg(l, t_MAT);
1265 679 : if (l==1) return H;
1266 679 : m = lgcols(Hp);
1267 679 : n = deg + 3;
1268 2268 : for (j=1; j<l; j++)
1269 : {
1270 1589 : GEN cp = gel(Hp,j);
1271 1589 : GEN c = cgetg(m, t_COL);
1272 1589 : gel(H,j) = c;
1273 24465 : for (i=1; i<m; i++)
1274 : {
1275 22876 : GEN dp = gel(cp, i);
1276 22876 : long l = lg(dp);
1277 22876 : GEN d = cgetg(n, t_POL);
1278 22876 : gel(c, i) = d;
1279 22876 : d[1] = dp[1] | evalsigne(1);
1280 46459 : for (k=2; k<l; k++)
1281 23583 : gel(d,k) = stoi(Fl_center(dp[k], p, lim));
1282 45493 : for ( ; k<n; k++)
1283 22617 : gel(d,k) = gen_0;
1284 : }
1285 : }
1286 679 : return H;
1287 : }
1288 :
1289 : int
1290 653 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1291 : {
1292 653 : GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1293 653 : ulong qinv = Fl_inv(umodiu(q,p), p);
1294 653 : long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
1295 653 : int stable = 1;
1296 2225 : for (j=1; j<l; j++)
1297 90418 : for (i=1; i<m; i++)
1298 : {
1299 88846 : GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
1300 88846 : long lh = lg(hp);
1301 246641 : for (k=2; k<lh; k++)
1302 : {
1303 157795 : v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
1304 157795 : if (v) { gel(h,k) = v; stable = 0; }
1305 : }
1306 108763 : for (; k<n; k++)
1307 : {
1308 19917 : v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
1309 19917 : if (v) { gel(h,k) = v; stable = 0; }
1310 : }
1311 : }
1312 653 : *ptq = qp; return stable;
1313 : }
1314 :
1315 : /* record the degrees of Euclidean remainders (make them as large as
1316 : * possible : smaller values correspond to a degenerate sequence) */
1317 : static void
1318 23720 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
1319 : {
1320 : long da,db,dc, ind;
1321 23720 : pari_sp av = avma;
1322 :
1323 23720 : if (lgpol(a)==0 || lgpol(b)==0) return;
1324 22453 : da = degpol(a);
1325 22453 : db = degpol(b);
1326 22453 : if (db > da)
1327 0 : { swapspec(a,b, da,db); }
1328 22453 : else if (!da) return;
1329 22453 : ind = 0;
1330 145746 : while (db)
1331 : {
1332 123294 : GEN c = Flx_rem(a,b, p);
1333 123291 : a = b; b = c; dc = degpol(c);
1334 123292 : if (dc < 0) break;
1335 :
1336 123292 : ind++;
1337 123292 : if (dc > dglist[ind]) dglist[ind] = dc;
1338 123292 : if (gc_needed(av,2))
1339 : {
1340 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1341 0 : (void)gc_all(av, 2, &a,&b);
1342 : }
1343 123293 : db = dc; /* = degpol(b) */
1344 : }
1345 22452 : if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
1346 22452 : set_avma(av);
1347 : }
1348 : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
1349 : * "generic" degree sequence as given by dglist, set *Ci and return
1350 : * resultant(a,b). Modular version of Collins's subresultant */
1351 : static ulong
1352 2090589 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
1353 : {
1354 : long da,db,dc, ind;
1355 2090589 : ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
1356 2090589 : int s = 1;
1357 2090589 : pari_sp av = avma;
1358 :
1359 2090589 : *C0 = 1; *C1 = 0;
1360 2090589 : if (lgpol(a)==0 || lgpol(b)==0) return 0;
1361 2081125 : da = degpol(a);
1362 2081125 : db = degpol(b);
1363 2081129 : if (db > da)
1364 : {
1365 0 : swapspec(a,b, da,db);
1366 0 : if (both_odd(da,db)) s = -s;
1367 : }
1368 2081129 : else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
1369 2081129 : ind = 0;
1370 19829603 : while (db)
1371 : { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
1372 : * da = deg a, db = deg b */
1373 17753285 : GEN c = Flx_rem(a,b, p);
1374 17704262 : long delta = da - db;
1375 :
1376 17704262 : if (both_odd(da,db)) s = -s;
1377 17704612 : lb = Fl_mul(b[db+2], cb, p);
1378 17724606 : a = b; b = c; dc = degpol(c);
1379 17726421 : ind++;
1380 17726421 : if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
1381 17721489 : if (g == h)
1382 : { /* frequent */
1383 17661519 : ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
1384 17689927 : ca = cb;
1385 17689927 : cb = cc;
1386 : }
1387 : else
1388 : {
1389 59970 : ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
1390 59970 : ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
1391 59970 : ca = cb;
1392 59970 : cb = Fl_div(cc, ghdelta, p);
1393 : }
1394 17749805 : da = db; /* = degpol(a) */
1395 17749805 : db = dc; /* = degpol(b) */
1396 :
1397 17749805 : g = lb;
1398 17749805 : if (delta == 1)
1399 17649400 : h = g; /* frequent */
1400 : else
1401 100405 : h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
1402 :
1403 17749732 : if (gc_needed(av,2))
1404 : {
1405 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1406 0 : (void)gc_all(av, 2, &a,&b);
1407 : }
1408 : }
1409 2076318 : if (da > 1) return 0; /* Failure */
1410 : /* last nonconstant polynomial has degree 1 */
1411 2076318 : *C0 = Fl_mul(ca, a[2], p);
1412 2076312 : *C1 = Fl_mul(ca, a[3], p);
1413 2076310 : res = Fl_mul(cb, b[2], p);
1414 2076304 : if (s == -1) res = p - res;
1415 2076304 : return gc_ulong(av,res);
1416 : }
1417 :
1418 : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
1419 : * Return 0 in case of degree drop. */
1420 : static GEN
1421 2114344 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
1422 : {
1423 : GEN z;
1424 2114344 : long i, lb = lg(Q);
1425 2114344 : ulong leadz = Flx_eval(leading_coeff(Q), x, p);
1426 2114210 : long vs=mael(Q,2,1);
1427 2114210 : if (!leadz) return zero_Flx(vs);
1428 :
1429 2103550 : z = cgetg(lb, t_VECSMALL); z[1] = vs;
1430 20122646 : for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
1431 2102526 : z[i] = leadz; return z;
1432 : }
1433 :
1434 : GEN
1435 2072 : FpXY_FpXQ_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
1436 : {
1437 2072 : pari_sp av = avma;
1438 2072 : long i, lb = lg(Q);
1439 : GEN z;
1440 2072 : if (lb == 2) return pol_0(vx);
1441 2072 : z = gel(Q, lb-1);
1442 2072 : if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1443 :
1444 2072 : if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1445 48636 : for (i=lb-2; i>=2; i--)
1446 : {
1447 46564 : GEN c = gel(Q,i);
1448 46564 : z = FqX_Fq_mul(z, y, T, p);
1449 46564 : z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
1450 : }
1451 2072 : return gc_upto(av, z);
1452 : }
1453 :
1454 : static GEN
1455 1296719 : ZX_norml1(GEN x)
1456 : {
1457 1296719 : long i, l = lg(x);
1458 : GEN s;
1459 :
1460 1296719 : if (l == 2) return gen_0;
1461 1208771 : s = gel(x, l-1); /* != 0 */
1462 2712574 : for (i = l-2; i > 1; i--) {
1463 1503811 : GEN xi = gel(x,i);
1464 1503811 : if (!signe(xi)) continue;
1465 1224265 : s = addii_sign(s,1, xi,1);
1466 : }
1467 1208763 : return s;
1468 : }
1469 : /* x >= 0, y != 0, return x + |y| */
1470 : static GEN
1471 25974 : addii_abs(GEN x, GEN y)
1472 : {
1473 25974 : if (!signe(x)) return absi_shallow(y);
1474 16443 : return addii_sign(x,1, y,1);
1475 : }
1476 :
1477 : /* x a ZX, return sum_{i >= k} |x[i]| binomial(i, k) */
1478 : static GEN
1479 32067 : ZX_norml1_1(GEN x, long k)
1480 : {
1481 32067 : long i, d = degpol(x);
1482 : GEN s, C; /* = binomial(i, k) */
1483 :
1484 32068 : if (!d || k > d) return gen_0;
1485 32068 : s = absi_shallow(gel(x, k+2)); /* may be 0 */
1486 32070 : C = gen_1;
1487 68921 : for (i = k+1; i <= d; i++) {
1488 36852 : GEN xi = gel(x,i+2);
1489 36852 : if (k) C = diviuexact(muliu(C, i), i-k);
1490 36853 : if (signe(xi)) s = addii_abs(s, mulii(C, xi));
1491 : }
1492 32069 : return s;
1493 : }
1494 : /* x has non-negative real coefficients */
1495 : static GEN
1496 3283 : RgX_norml1_1(GEN x, long k)
1497 : {
1498 3283 : long i, d = degpol(x);
1499 : GEN s, C; /* = binomial(i, k) */
1500 :
1501 3283 : if (!d || k > d) return gen_0;
1502 3283 : s = gel(x, k+2); /* may be 0 */
1503 3283 : C = gen_1;
1504 9198 : for (i = k+1; i <= d; i++) {
1505 5915 : GEN xi = gel(x,i+2);
1506 5915 : if (k) C = diviuexact(muliu(C, i), i-k);
1507 5915 : if (!gequal0(xi)) s = gadd(s, gmul(C, xi));
1508 : }
1509 3283 : return s;
1510 : }
1511 :
1512 : /* N_2(A)^2 */
1513 : static GEN
1514 9061 : sqrN2(GEN A, long prec)
1515 : {
1516 9061 : pari_sp av = avma;
1517 9061 : long i, l = lg(A);
1518 9061 : GEN a = gen_0;
1519 44189 : for (i = 2; i < l; i++)
1520 : {
1521 35128 : a = gadd(a, gabs(gnorm(gel(A,i)), prec));
1522 35128 : if (gc_needed(av,1))
1523 : {
1524 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1525 0 : a = gc_upto(av, a);
1526 : }
1527 : }
1528 9061 : return a;
1529 : }
1530 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1531 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1532 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1533 : * Return e such that Res(A, B) < 2^e */
1534 : static GEN
1535 8207 : RgX_RgXY_ResBound(GEN A, GEN B, long prec)
1536 : {
1537 8207 : pari_sp av = avma;
1538 8207 : GEN b = gen_0, bnd;
1539 8207 : long i, lB = lg(B);
1540 32167 : for (i=2; i<lB; i++)
1541 : {
1542 23960 : GEN t = gel(B,i);
1543 23960 : if (typ(t) == t_POL) t = gnorml1(t, prec);
1544 23960 : b = gadd(b, gabs(gsqr(t), prec));
1545 23960 : if (gc_needed(av,1))
1546 : {
1547 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1548 0 : b = gc_upto(av, b);
1549 : }
1550 : }
1551 8207 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1552 : gpowgs(b, degpol(A))), prec);
1553 8207 : return gc_upto(av, bnd);
1554 : }
1555 : /* A,B in C[X] return RgX_RgXY_ResBound(A, B(x+y)) */
1556 : static GEN
1557 854 : RgX_RgXY_ResBound_1(GEN A, GEN B, long prec)
1558 : {
1559 854 : pari_sp av = avma, av2;
1560 854 : GEN b = gen_0, bnd;
1561 854 : long i, lB = lg(B);
1562 854 : B = shallowcopy(B);
1563 4137 : for (i=2; i<lB; i++) gel(B,i) = gabs(gel(B,i), prec);
1564 854 : av2 = avma;
1565 4137 : for (i=2; i<lB; i++)
1566 : {
1567 3283 : b = gadd(b, gsqr(RgX_norml1_1(B, i-2)));
1568 3283 : if (gc_needed(av2,1))
1569 : {
1570 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1571 0 : b = gc_upto(av2, b);
1572 : }
1573 : }
1574 854 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1575 : gpowgs(b, degpol(A))), prec);
1576 854 : return gc_upto(av, bnd);
1577 : }
1578 :
1579 : /* log2 N_2(A)^2 */
1580 : static double
1581 176224 : log2N2(GEN A)
1582 : {
1583 176224 : pari_sp av = avma;
1584 176224 : long i, l = lg(A);
1585 176224 : GEN a = gen_0;
1586 1334632 : for (i=2; i < l; i++)
1587 : {
1588 1158417 : a = addii(a, sqri(gel(A,i)));
1589 1158408 : if (gc_needed(av,1))
1590 : {
1591 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1592 0 : a = gc_upto(av, a);
1593 : }
1594 : }
1595 176215 : return gc_double(av, dbllog2(a));
1596 : }
1597 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1598 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1599 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1600 : * Return e such that Res(A, B) < 2^e */
1601 :
1602 : static double
1603 2187045 : resbound(GEN B)
1604 : {
1605 2187045 : pari_sp av = avma;
1606 2187045 : long i, lB = lg(B);
1607 2187045 : GEN b = gen_0;
1608 9735430 : for (i=2; i<lB; i++)
1609 : {
1610 7548517 : GEN t = gel(B,i);
1611 7548517 : if (typ(t) == t_POL) t = ZX_norml1(t);
1612 7548517 : b = addii(b, sqri(t));
1613 7548387 : if (gc_needed(av,1))
1614 : {
1615 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1616 0 : b = gc_upto(av, b);
1617 : }
1618 : }
1619 2186913 : return gc_double(av, dbllog2(b));
1620 : }
1621 :
1622 : ulong
1623 166003 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
1624 : {
1625 166003 : pari_sp av = avma;
1626 166003 : double logb = resbound(B);
1627 : long i;
1628 166001 : if (dB) logb -= 2 * dbllog2(dB);
1629 166001 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * logb) / 2);
1630 166001 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1631 : }
1632 : static ulong
1633 1010526 : ZXX_ResBound(GEN A, GEN B)
1634 : {
1635 1010526 : pari_sp av = avma;
1636 1010526 : double loga = resbound(A), logb = resbound(B);
1637 1010525 : long i = (long)((degpol(B) * loga + degpol(A) * logb) / 2);
1638 1010527 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1639 : }
1640 :
1641 : /* A,B ZX. Return ZX_ZXY_ResBound(A(x), B(x+y)) */
1642 : static ulong
1643 10217 : ZX_ZXY_ResBound_1(GEN A, GEN B)
1644 : {
1645 10217 : pari_sp av = avma;
1646 10217 : GEN b = gen_0;
1647 10217 : long i, lB = lg(B);
1648 42282 : for (i=2; i<lB; i++)
1649 : {
1650 32067 : b = addii(b, sqri(ZX_norml1_1(B, i-2)));
1651 32065 : if (gc_needed(av,1))
1652 : {
1653 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1654 0 : b = gc_upto(av, b);
1655 : }
1656 : }
1657 10215 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * dbllog2(b)) / 2);
1658 10218 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1659 : }
1660 : /* special case B = A' */
1661 : static ulong
1662 1140651 : ZX_discbound(GEN A)
1663 : {
1664 1140651 : pari_sp av = avma;
1665 1140651 : GEN a = gen_0, b = gen_0;
1666 1140651 : long i , lA = lg(A), dA = degpol(A);
1667 : double loga, logb;
1668 6812392 : for (i = 2; i < lA; i++)
1669 : {
1670 5671848 : GEN c = sqri(gel(A,i));
1671 5671669 : a = addii(a, c);
1672 5671754 : if (i > 2) b = addii(b, mulii(c, sqru(i-2)));
1673 5671745 : if (gc_needed(av,1))
1674 : {
1675 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_discbound i = %ld",i);
1676 0 : (void)gc_all(av, 2, &a, &b);
1677 : }
1678 : }
1679 1140544 : loga = dbllog2(a);
1680 1140601 : logb = dbllog2(b); set_avma(av);
1681 1140613 : i = (long)(((dA-1) * loga + dA * logb) / 2);
1682 1140613 : return (i <= 0)? 1: 1 + (ulong)i;
1683 : }
1684 :
1685 : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
1686 : static ulong
1687 5519333 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong pi, ulong la)
1688 : {
1689 5519333 : GEN ev = FlxY_evalx_pre(b, n, p, pi);
1690 5519626 : long drop = lg(b) - lg(ev);
1691 5519626 : ulong r = Flx_resultant_pre(a, ev, p, pi);
1692 5519150 : if (drop && la != 1) r = Fl_mul(r, Fl_powu_pre(la, drop, p, pi), p);
1693 5519179 : return r;
1694 : }
1695 : /* return Res_Y(a(x,Y), b(x,Y)) | x=n over Fp;
1696 : * we write A = a(Y,n), B = b(Y,n) */
1697 : static ulong
1698 4348104 : FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong pi)
1699 : {
1700 4348104 : long da = degpol(a), db = degpol(b), dA, dB, dropa, dropb;
1701 : ulong r, lA, lB;
1702 : GEN A, B;
1703 4348092 : if (!da && !db) return 1UL;
1704 4348085 : if (da < 0 || db < 0) return 0UL;
1705 : /* a and b are non-zero */
1706 4348085 : A = FlxY_evalx_pre(a, n, p, pi); dA = degpol(A); lA = uel(A,dA+2);
1707 4348114 : B = FlxY_evalx_pre(b, n, p, pi); dB = degpol(B); lB = uel(B,dB+2);
1708 4348117 : dropa = da - dA;
1709 4348117 : dropb = db - dB;
1710 4348117 : if (dropa && dropb) return 0UL;
1711 : /* A = 0, implies dropa && !dropb */
1712 4348110 : if (dA < 0) return db? 0UL: Fl_powu_pre(lB, da, p, pi);
1713 : /* B = 0, implies dropb && !dropa */
1714 4348096 : if (dB < 0) return da? 0UL: Fl_powu_pre(lA, db, p, pi);
1715 4348089 : r = Flx_resultant_pre(A, B, p, pi);
1716 4348060 : if (dropa)
1717 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
1718 7 : ulong c = Fl_powu_pre(odd(db)? p - lB: lB, dropa, p, pi);
1719 7 : if (c != 1UL) r = Fl_mul_pre(r, c, p, pi);
1720 : }
1721 4348053 : else if (dropb)
1722 : { /* multiply by lc(A)^(deg B - deg b) */
1723 0 : ulong c = Fl_powu_pre(lA, dropb, p, pi);
1724 0 : if (c != 1UL) r = Fl_mul_pre(r, c, p, pi);
1725 : }
1726 4348059 : return r;
1727 : }
1728 :
1729 : static GEN
1730 284 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
1731 : {
1732 284 : GEN ev = FpXY_evaly(b, n, p, vX);
1733 284 : long drop = db-degpol(ev);
1734 284 : GEN r = FpX_resultant(a, ev, p);
1735 284 : if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
1736 284 : return r;
1737 : }
1738 :
1739 : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
1740 : /* Return a Fly */
1741 : static GEN
1742 367281 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, ulong pi, long dres, long sx)
1743 : {
1744 : long i;
1745 367281 : ulong n, la = Flx_lead(a);
1746 367280 : GEN x = cgetg(dres+2, t_VECSMALL);
1747 367280 : GEN y = cgetg(dres+2, t_VECSMALL);
1748 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1749 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1750 2946690 : for (i=0,n = 1; i < dres; n++)
1751 : {
1752 2579413 : x[++i] = n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1753 2579307 : x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1754 : }
1755 367277 : if (i == dres)
1756 : {
1757 361893 : x[++i] = 0; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1758 : }
1759 367276 : return Flv_polint(x,y, p, sx);
1760 : }
1761 :
1762 : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
1763 : /* Return a Fly */
1764 : static GEN
1765 1118088 : FlxX_resultant_polint(GEN a, GEN b, ulong p, ulong pi, long dres, long sx)
1766 : {
1767 : long i;
1768 : ulong n;
1769 1118088 : GEN x = cgetg(dres+2, t_VECSMALL);
1770 1118086 : GEN y = cgetg(dres+2, t_VECSMALL);
1771 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1772 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1773 2998451 : for (i=0,n = 1; i < dres; n++)
1774 : {
1775 1880362 : x[++i] = n; y[i] = FlxY_eval_resultant(a,b, x[i], p,pi);
1776 1880350 : x[++i] = p-n; y[i] = FlxY_eval_resultant(a,b, x[i], p,pi);
1777 : }
1778 1118089 : if (i == dres)
1779 : {
1780 587427 : x[++i] = 0; y[i] = FlxY_eval_resultant(a,b, x[i], p,pi);
1781 : }
1782 1118089 : return Flv_polint(x,y, p, sx);
1783 : }
1784 :
1785 : static GEN
1786 7568 : FlxX_pseudorem(GEN x, GEN y, ulong p, ulong pi)
1787 : {
1788 7568 : long vx = varn(x), dx, dy, dz, i, lx, dp;
1789 7568 : pari_sp av = avma, av2;
1790 :
1791 7568 : if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
1792 7568 : (void)new_chunk(2);
1793 7567 : dx=degpol(x); x = RgX_recip_i(x)+2;
1794 7568 : dy=degpol(y); y = RgX_recip_i(y)+2; dz=dx-dy; dp = dz+1;
1795 7568 : av2 = avma;
1796 : for (;;)
1797 : {
1798 62068 : gel(x,0) = Flx_neg(gel(x,0), p); dp--;
1799 232545 : for (i=1; i<=dy; i++)
1800 170246 : gel(x,i) = Flx_add( Flx_mul_pre(gel(y,0), gel(x,i), p, pi),
1801 170443 : Flx_mul_pre(gel(x,0), gel(y,i), p, pi), p );
1802 1131393 : for ( ; i<=dx; i++)
1803 1070049 : gel(x,i) = Flx_mul_pre(gel(y,0), gel(x,i), p, pi);
1804 66004 : do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
1805 61344 : if (dx < dy) break;
1806 53775 : if (gc_needed(av2,1))
1807 : {
1808 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
1809 0 : gc_slice(av2,x,dx+1);
1810 : }
1811 : }
1812 7569 : if (dx < 0) return zero_Flx(0);
1813 7569 : lx = dx+3; x -= 2;
1814 7569 : x[0]=evaltyp(t_POL) | _evallg(lx);
1815 7569 : x[1]=evalsigne(1) | evalvarn(vx);
1816 7569 : x = RgX_recip_i(x);
1817 7570 : if (dp)
1818 : { /* multiply by y[0]^dp [beware dummy vars from FpX_FpXY_resultant] */
1819 1977 : GEN t = Flx_powu_pre(gel(y,0), dp, p, pi);
1820 7911 : for (i=2; i<lx; i++) gel(x,i) = Flx_mul_pre(gel(x,i), t, p, pi);
1821 : }
1822 7574 : return gc_GEN(av, x);
1823 : }
1824 :
1825 : static GEN
1826 1991 : Flx_Lazard(GEN x, GEN y, long n, ulong p, ulong pi)
1827 : {
1828 : long a;
1829 : GEN c;
1830 1991 : if (n == 1) return x;
1831 1991 : a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
1832 1991 : c=x; n-=a;
1833 1991 : y = Flx_get_red_pre(y, p, pi);
1834 8917 : while (a>1)
1835 : {
1836 6925 : a>>=1; c = Flx_div_pre(Flx_sqr_pre(c, p, pi), y, p, pi);
1837 6924 : if (n>=a) { c = Flx_div_pre(Flx_mul_pre(c,x,p,pi),y, p,pi); n -= a; }
1838 : }
1839 1992 : return c;
1840 : }
1841 :
1842 : /* return a Flx */
1843 : static GEN
1844 2533 : FlxX_resultant_subres(GEN u, GEN v, ulong p, ulong pi, long sx)
1845 : {
1846 2533 : pari_sp av = avma, av2;
1847 : long degq, dx, dy, du, dv, dr, signh;
1848 : GEN z, g, h, r, p1;
1849 :
1850 2533 : dx = degpol(u); dy = degpol(v); signh = 1;
1851 2534 : if (dx < dy)
1852 : {
1853 7 : swap(u,v); lswap(dx,dy);
1854 7 : if (both_odd(dx, dy)) signh = -signh;
1855 : }
1856 2534 : if (dy < 0) return zero_Flx(sx);
1857 2534 : if (dy==0) return gc_upto(av, Flx_powu_pre(gel(v,2),dx,p,pi));
1858 :
1859 2534 : g = h = pol1_Flx(sx); av2 = avma;
1860 : for(;;)
1861 : {
1862 7568 : r = FlxX_pseudorem(u,v,p,pi); dr = lg(r);
1863 7576 : if (dr == 2) { set_avma(av); return zero_Flx(sx); }
1864 7576 : du = degpol(u); dv = degpol(v); degq = du-dv;
1865 7575 : u = v; p1 = g; g = leading_coeff(u);
1866 7575 : switch(degq)
1867 : {
1868 0 : case 0: break;
1869 5582 : case 1:
1870 5582 : p1 = Flx_mul_pre(h,p1, p, pi); h = g; break;
1871 1993 : default:
1872 1993 : p1 = Flx_mul_pre(Flx_powu_pre(h,degq,p,pi), p1, p, pi);
1873 1991 : h = Flx_Lazard(g,h,degq, p, pi);
1874 : }
1875 7566 : if (both_odd(du,dv)) signh = -signh;
1876 7565 : v = FlxY_Flx_div_pre(r, p1, p, pi);
1877 7564 : if (dr==3) break;
1878 5030 : if (gc_needed(av2,1))
1879 : {
1880 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
1881 0 : (void)gc_all(av2,4, &u, &v, &g, &h);
1882 : }
1883 : }
1884 2534 : z = gel(v,2);
1885 2534 : if (dv > 1) z = Flx_Lazard(z,h,dv,p,pi);
1886 2534 : if (signh < 0) z = Flx_neg(z,p);
1887 2534 : return gc_upto(av, z);
1888 : }
1889 :
1890 : /* Return a Flx*/
1891 : GEN
1892 0 : FlxX_resultant_pre(GEN a, GEN b, ulong p, ulong pi, long sx)
1893 : {
1894 0 : pari_sp ltop = avma;
1895 0 : long da = degpol(a), db = degpol(b);
1896 : ulong dres;
1897 : GEN z;
1898 0 : if (da<0 || db<0) return pol0_Flx(sx);
1899 0 : dres = FlxY_degreex(a)*db+FlxY_degreex(b)*da;
1900 0 : z = dres >= p ? FlxX_resultant_subres(a, b, p, pi, sx)
1901 0 : : FlxX_resultant_polint(a, b, p, pi, dres, sx);
1902 0 : return gc_upto(ltop, z);
1903 : }
1904 :
1905 : GEN
1906 0 : FlxX_resultant(GEN a, GEN b, ulong p, long sx)
1907 0 : { return FlxX_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), sx); }
1908 :
1909 : /* Warning:
1910 : * This function switches between valid and invalid variable ordering*/
1911 : static GEN
1912 6143 : FlxY_to_FlyX(GEN b, long sv)
1913 : {
1914 6143 : long i, n=-1;
1915 6143 : long sw = b[1]&VARNBITS;
1916 20968 : for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
1917 6143 : return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
1918 : }
1919 :
1920 : /* Return a Fly*/
1921 : GEN
1922 6144 : Flx_FlxY_resultant(GEN a, GEN b, ulong p)
1923 : {
1924 6144 : pari_sp ltop=avma;
1925 6144 : long dres = degpol(a)*degpol(b);
1926 6143 : long sx=a[1], sy=b[1]&VARNBITS;
1927 6143 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1928 : GEN z;
1929 6143 : b = FlxY_to_FlyX(b,sx);
1930 6143 : if ((ulong)dres >= p)
1931 2532 : z = FlxX_resultant_subres(Fly_to_FlxY(a, sy), b, p, pi, sy);
1932 : else
1933 3611 : z = Flx_FlxY_resultant_polint(a, b, p, pi, (ulong)dres, sy);
1934 6145 : return gc_upto(ltop,z);
1935 : }
1936 :
1937 : /* Return a t_POL in variable vc whose coeffs are the coeffs of b in
1938 : * variable v; vc must have higher priority than all variables occuring in b. */
1939 : GEN
1940 146150 : swap_vars(GEN b, long v, long vc)
1941 : {
1942 146150 : long i, n = RgX_degree(b, v);
1943 : GEN c, x;
1944 146147 : if (n < 0) return pol_0(vc);
1945 146147 : c = cgetg(n+3, t_POL); x = c + 2;
1946 146146 : c[1] = evalsigne(1) | evalvarn(vc);
1947 964237 : for (i = 0; i <= n; i++) gel(x,i) = polcoef_i(b, i, v);
1948 146145 : return c;
1949 : }
1950 :
1951 : /* assume varn(b) << varn(a) */
1952 : /* return a FpY*/
1953 : GEN
1954 15 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
1955 : {
1956 15 : long i,n,dres, db, vY = varn(b), vX = varn(a);
1957 : GEN la,x,y;
1958 :
1959 15 : if (lgefint(p) == 3)
1960 : {
1961 0 : ulong pp = uel(p,2);
1962 0 : b = ZXX_to_FlxX(b, pp, vX);
1963 0 : a = ZX_to_Flx(a, pp);
1964 0 : x = Flx_FlxY_resultant(a, b, pp);
1965 0 : return Flx_to_ZX(x);
1966 : }
1967 15 : db = RgXY_degreex(b);
1968 15 : dres = degpol(a)*degpol(b);
1969 15 : la = leading_coeff(a);
1970 15 : x = cgetg(dres+2, t_VEC);
1971 15 : y = cgetg(dres+2, t_VEC);
1972 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1973 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1974 157 : for (i=0,n = 1; i < dres; n++)
1975 : {
1976 142 : gel(x,++i) = utoipos(n);
1977 142 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1978 142 : gel(x,++i) = subiu(p,n);
1979 142 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1980 : }
1981 15 : if (i == dres)
1982 : {
1983 0 : gel(x,++i) = gen_0;
1984 0 : gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
1985 : }
1986 15 : return FpV_polint(x,y, p, vY);
1987 : }
1988 :
1989 : GEN
1990 93 : FpX_composedsum(GEN P, GEN Q, GEN p)
1991 : {
1992 93 : pari_sp av = avma;
1993 93 : if (lgefint(p)==3)
1994 : {
1995 0 : ulong pp = p[2];
1996 0 : GEN z = Flx_composedsum(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
1997 0 : return gc_upto(av, Flx_to_ZX(z));
1998 : }
1999 : else
2000 : {
2001 93 : long n = 1+ degpol(P)*degpol(Q);
2002 93 : GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
2003 93 : GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
2004 93 : GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
2005 93 : GEN lead = Fp_mul(Fp_powu(leading_coeff(P),degpol(Q), p),
2006 93 : Fp_powu(leading_coeff(Q),degpol(P), p), p);
2007 93 : GEN R = FpX_fromNewton(L, p);
2008 93 : return gc_upto(av, FpX_Fp_mul(R, lead, p));
2009 : }
2010 : }
2011 :
2012 : GEN
2013 0 : FpX_composedprod(GEN P, GEN Q, GEN p)
2014 : {
2015 0 : pari_sp av = avma;
2016 0 : if (lgefint(p)==3)
2017 : {
2018 0 : ulong pp = p[2];
2019 0 : GEN z = Flx_composedprod(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
2020 0 : return gc_upto(av, Flx_to_ZX(z));
2021 : }
2022 : else
2023 : {
2024 0 : long n = 1+ degpol(P)*degpol(Q);
2025 0 : GEN L = FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
2026 0 : return gc_upto(av,FpX_fromNewton(L, p));
2027 : }
2028 : }
2029 :
2030 : static GEN
2031 93 : _FpX_composedsum(void *E, GEN a, GEN b)
2032 93 : { return FpX_composedsum(a,b, (GEN)E); }
2033 :
2034 : GEN
2035 1588 : FpXV_composedsum(GEN V, GEN p)
2036 : {
2037 1588 : if (lgefint(p)==3)
2038 : {
2039 0 : ulong pp = p[2];
2040 0 : return Flx_to_ZX(FlxV_composedsum(ZXV_to_FlxV(V, pp), pp));
2041 : }
2042 1588 : return gen_product(V, (void *)p, &_FpX_composedsum);
2043 : }
2044 :
2045 : /* 0, 1, -1, 2, -2, ... */
2046 : #define next_lambda(a) (a>0 ? -a : 1-a)
2047 :
2048 : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
2049 : * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
2050 : * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
2051 : * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
2052 : * the Last nonconstant polynomial in the Euclidean Remainder Sequence */
2053 : static GEN
2054 22337 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
2055 : {
2056 : ulong bound, dp;
2057 22337 : pari_sp av = avma, av2 = 0;
2058 22337 : long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
2059 : long stable, checksqfree, i,n, cnt, degB;
2060 22337 : long v, vX = varn(B0), vY = varn(A); /* vY < vX */
2061 : GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
2062 : forprime_t S;
2063 :
2064 22337 : if (degA == 1)
2065 : {
2066 1260 : GEN a1 = gel(A,3), a0 = gel(A,2);
2067 1260 : B = lambda? RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
2068 1260 : H = gsubst(B, vY, gdiv(gneg(a0),a1));
2069 1260 : if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
2070 1260 : *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
2071 1260 : return gc_all(av, 2, &H, LERS);
2072 : }
2073 :
2074 21077 : dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
2075 21077 : C0 = cgetg(dres+2, t_VECSMALL);
2076 21077 : C1 = cgetg(dres+2, t_VECSMALL);
2077 21077 : dglist = cgetg(dres+1, t_VECSMALL);
2078 21077 : x = cgetg(dres+2, t_VECSMALL);
2079 21077 : y = cgetg(dres+2, t_VECSMALL);
2080 21077 : B0 = leafcopy(B0);
2081 21077 : A = leafcopy(A);
2082 21077 : B = B0;
2083 21077 : v = fetch_var_higher(); setvarn(A,v);
2084 : /* make sure p large enough */
2085 22236 : INIT:
2086 : /* always except the first time */
2087 22236 : if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
2088 22236 : if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
2089 22236 : B = swap_vars(B, vY, v);
2090 : /* B0(lambda v + x, v) */
2091 22236 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
2092 22236 : av2 = avma;
2093 :
2094 22236 : if (degA <= 3)
2095 : { /* sub-resultant faster for small degrees */
2096 11242 : H = RgX_resultant_all(A,B,&q);
2097 11242 : if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
2098 10297 : H0 = gel(q,2);
2099 10297 : if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
2100 10297 : H1 = gel(q,3);
2101 10297 : if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
2102 10297 : if (!ZX_is_squarefree(H)) goto INIT;
2103 10255 : goto END;
2104 : }
2105 :
2106 10994 : H = H0 = H1 = NULL;
2107 10994 : degB = degpol(B);
2108 10994 : bound = ZX_ZXY_ResBound(A, B, NULL);
2109 10994 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
2110 10994 : dp = 1;
2111 10994 : init_modular_big(&S);
2112 10994 : for(cnt = 0, checksqfree = 1;;)
2113 49452 : {
2114 60446 : ulong p = u_forprime_next(&S);
2115 : GEN Hi;
2116 60446 : a = ZX_to_Flx(A, p);
2117 60446 : b = ZXX_to_FlxX(B, p, varn(A));
2118 60446 : if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
2119 60446 : if (checksqfree)
2120 : { /* find degree list for generic Euclidean Remainder Sequence */
2121 10994 : long goal = minss(degpol(a), degpol(b)); /* longest possible */
2122 74267 : for (n=1; n <= goal; n++) dglist[n] = 0;
2123 10994 : setlg(dglist, 1);
2124 24132 : for (n=0; n <= dres; n++)
2125 : {
2126 23719 : ev = FlxY_evalx_drop(b, n, p);
2127 23720 : Flx_resultant_set_dglist(a, ev, dglist, p);
2128 23719 : if (lg(dglist)-1 == goal) break;
2129 : }
2130 : /* last pol in ERS has degree > 1 ? */
2131 10994 : goal = lg(dglist)-1;
2132 10994 : if (degpol(B) == 1) { if (!goal) goto INIT; }
2133 : else
2134 : {
2135 10931 : if (goal <= 1) goto INIT;
2136 10847 : if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
2137 : }
2138 10910 : if (DEBUGLEVEL>4)
2139 0 : err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
2140 : }
2141 :
2142 2150992 : for (i=0,n = 0; i <= dres; n++)
2143 : {
2144 2090633 : ev = FlxY_evalx_drop(b, n, p);
2145 2090576 : x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
2146 2090630 : if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
2147 : }
2148 60359 : Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
2149 60362 : Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
2150 60362 : if (!H && degpol(Hp) != dres) continue;
2151 60362 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
2152 60362 : if (checksqfree) {
2153 10910 : if (!Flx_is_squarefree(Hp, p)) goto INIT;
2154 10822 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
2155 10822 : checksqfree = 0;
2156 : }
2157 :
2158 60274 : if (!H)
2159 : { /* initialize */
2160 10822 : q = utoipos(p); stable = 0;
2161 10822 : H = ZX_init_CRT(Hp, p,vX);
2162 10822 : H0= ZX_init_CRT(H0p, p,vX);
2163 10822 : H1= ZX_init_CRT(H1p, p,vX);
2164 : }
2165 : else
2166 : {
2167 49452 : GEN qp = muliu(q,p);
2168 49452 : stable = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
2169 49452 : & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
2170 49452 : & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
2171 49452 : q = qp;
2172 : }
2173 : /* could make it probabilistic for H ? [e.g if stable twice, etc]
2174 : * Probabilistic anyway for H0, H1 */
2175 60274 : if (DEBUGLEVEL>5 && (stable || ++cnt==100))
2176 0 : { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
2177 60274 : if (stable && (ulong)expi(q) >= bound) break; /* DONE */
2178 49452 : if (gc_needed(av,2))
2179 : {
2180 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
2181 0 : (void)gc_all(av2, 4, &H, &q, &H0, &H1);
2182 : }
2183 : }
2184 21077 : END:
2185 21077 : if (DEBUGLEVEL>5) err_printf(" done\n");
2186 21077 : setvarn(H, vX); (void)delete_var();
2187 21077 : *LERS = mkvec2(H0,H1);
2188 21077 : *plambda = lambda; return gc_all(av, 2, &H, LERS);
2189 : }
2190 :
2191 : GEN
2192 60114 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
2193 : {
2194 60114 : if (LERS)
2195 : {
2196 22337 : if (!plambda)
2197 0 : pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
2198 22337 : return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
2199 : }
2200 37777 : return ZX_ZXY_rnfequation(A, B, plambda);
2201 : }
2202 :
2203 : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
2204 : * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
2205 : * squarefree */
2206 : GEN
2207 22600 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
2208 : {
2209 22600 : pari_sp av = avma;
2210 : GEN R, a;
2211 : long dA;
2212 : int delvar;
2213 :
2214 22600 : if (v < 0) v = 0;
2215 22600 : switch (typ(A))
2216 : {
2217 22600 : case t_POL: dA = degpol(A); if (dA > 0) break;
2218 0 : A = constant_coeff(A);
2219 0 : default:
2220 0 : if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
2221 0 : return gc_upto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
2222 : }
2223 22600 : delvar = 0;
2224 22600 : if (varncmp(varn(T), 0) <= 0)
2225 : {
2226 3681 : long v0 = fetch_var(); delvar = 1;
2227 3681 : T = leafcopy(T); setvarn(T,v0);
2228 3681 : A = leafcopy(A); setvarn(A,v0);
2229 : }
2230 22600 : R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
2231 22600 : if (delvar) (void)delete_var();
2232 22600 : setvarn(R, v); a = leading_coeff(T);
2233 22600 : if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
2234 22600 : return gc_upto(av, R);
2235 : }
2236 :
2237 : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
2238 : GEN
2239 1247304 : ZXQ_charpoly(GEN A, GEN T, long v)
2240 : {
2241 1247304 : return (degpol(T) < 16) ? RgXQ_charpoly_i(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
2242 : }
2243 :
2244 : GEN
2245 9793 : QXQ_charpoly(GEN A, GEN T, long v)
2246 : {
2247 9793 : pari_sp av = avma;
2248 9793 : GEN den, B = Q_remove_denom(A, &den);
2249 9793 : GEN P = ZXQ_charpoly(B, T, v);
2250 9793 : return gc_GEN(av, den ? RgX_rescale(P, ginv(den)): P);
2251 : }
2252 :
2253 : static ulong
2254 3870678 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
2255 : {
2256 3870678 : pari_sp av = avma;
2257 3870678 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2258 : ulong H, dp;
2259 3870626 : if (dropa && dropb) return 0; /* p | lc(A), p | lc(B) */
2260 3870626 : H = Flx_resultant(a, b, p);
2261 3870358 : if (dropa)
2262 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2263 0 : ulong c = b[degB+2]; /* lc(B) */
2264 0 : if (odd(degB)) c = p - c;
2265 0 : c = Fl_powu(c, dropa, p);
2266 0 : if (c != 1) H = Fl_mul(H, c, p);
2267 : }
2268 3870358 : else if (dropb)
2269 : { /* multiply by lc(A)^(deg B - deg b) */
2270 0 : ulong c = a[degA+2]; /* lc(A) */
2271 0 : c = Fl_powu(c, dropb, p);
2272 0 : if (c != 1) H = Fl_mul(H, c, p);
2273 : }
2274 3870362 : dp = dB ? umodiu(dB, p): 1;
2275 3870362 : if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
2276 3870364 : return gc_ulong(av, H);
2277 : }
2278 :
2279 : /* If B=NULL, assume B=A' */
2280 : static GEN
2281 1498547 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
2282 : {
2283 1498547 : pari_sp av = avma, av2;
2284 1498547 : long degA, degB, i, n = lg(P)-1;
2285 : GEN H, T;
2286 :
2287 1498547 : degA = degpol(A);
2288 1498545 : degB = B? degpol(B): degA - 1;
2289 1498544 : if (n == 1)
2290 : {
2291 812715 : ulong Hp, p = uel(P,1);
2292 812715 : GEN a = ZX_to_Flx(A, p), b = B? ZX_to_Flx(B, p): Flx_deriv(a, p);
2293 812694 : Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
2294 812705 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
2295 : }
2296 685829 : T = ZV_producttree(P);
2297 685824 : A = ZX_nv_mod_tree(A, P, T);
2298 685806 : if (B) B = ZX_nv_mod_tree(B, P, T);
2299 685806 : H = cgetg(n+1, t_VECSMALL); av2 = avma;
2300 3743509 : for(i=1; i <= n; i++, set_avma(av2))
2301 : {
2302 3057683 : ulong p = P[i];
2303 3057683 : GEN a = gel(A,i), b = B? gel(B,i): Flx_deriv(a, p);
2304 3057979 : H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
2305 : }
2306 685826 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
2307 685831 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2308 : }
2309 :
2310 : GEN
2311 1498566 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
2312 : {
2313 1498566 : GEN V = cgetg(3, t_VEC);
2314 1498552 : if (typ(B) == t_INT) B = NULL;
2315 1498552 : if (!signe(dB)) dB = NULL;
2316 1498552 : gel(V,1) = ZX_resultant_slice(A, B, dB, P, &gel(V,2));
2317 1498532 : return V;
2318 : }
2319 :
2320 : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
2321 : * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
2322 : GEN
2323 1352872 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
2324 : {
2325 1352872 : pari_sp av = avma;
2326 : forprime_t S;
2327 : GEN H, worker;
2328 1352872 : if (!B && degpol(A)==2)
2329 : {
2330 116190 : GEN a = gel(A,4), b = gel(A,3), c = gel(A,2);
2331 116190 : H = mulii(a, subii(shifti(mulii(a, c), 2), sqri(b)));
2332 116181 : if (dB) H = diviiexact(H, sqri(dB));
2333 116181 : return gc_INT(av, H);
2334 : }
2335 1236685 : if (B)
2336 : {
2337 150227 : long a = degpol(A), b = degpol(B);
2338 150227 : if (a < 0 || b < 0) return gen_0;
2339 150197 : if (!a) return powiu(gel(A,2), b);
2340 150197 : if (!b) return powiu(gel(B,2), a);
2341 149572 : if (minss(a, b) <= 1)
2342 : {
2343 73120 : H = RgX_resultant_all(A, B, NULL);
2344 73120 : if (dB) H = diviiexact(H, powiu(dB, a));
2345 73120 : return gc_INT(av, H);
2346 : }
2347 76452 : if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
2348 : }
2349 1162921 : worker = snm_closure(is_entry("_ZX_resultant_worker"),
2350 : mkvec3(A, B? B: gen_0, dB? dB: gen_0));
2351 1162975 : init_modular_big(&S);
2352 1162947 : H = gen_crt("ZX_resultant_all", worker, &S, dB, bound, 0, NULL,
2353 : ZV_chinese_center, Fp_center);
2354 1162969 : return gc_INT(av, H);
2355 : }
2356 :
2357 : /* A0 and B0 in Q[X] */
2358 : GEN
2359 56 : QX_resultant(GEN A0, GEN B0)
2360 : {
2361 : GEN s, a, b, A, B;
2362 56 : pari_sp av = avma;
2363 :
2364 56 : A = Q_primitive_part(A0, &a);
2365 56 : B = Q_primitive_part(B0, &b);
2366 56 : s = ZX_resultant(A, B);
2367 56 : if (!signe(s)) { set_avma(av); return gen_0; }
2368 56 : if (a) s = gmul(s, gpowgs(a,degpol(B)));
2369 56 : if (b) s = gmul(s, gpowgs(b,degpol(A)));
2370 56 : return gc_upto(av, s);
2371 : }
2372 :
2373 : GEN
2374 55867 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
2375 :
2376 : GEN
2377 0 : QXQ_intnorm(GEN A, GEN B)
2378 : {
2379 : GEN c, n, R, lB;
2380 0 : long dA = degpol(A), dB = degpol(B);
2381 0 : pari_sp av = avma;
2382 0 : if (dA < 0) return gen_0;
2383 0 : A = Q_primitive_part(A, &c);
2384 0 : if (!c || typ(c) == t_INT) {
2385 0 : n = c;
2386 0 : R = ZX_resultant(B, A);
2387 : } else {
2388 0 : n = gel(c,1);
2389 0 : R = ZX_resultant_all(B, A, gel(c,2), 0);
2390 : }
2391 0 : if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
2392 0 : lB = leading_coeff(B);
2393 0 : if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
2394 0 : return gc_INT(av, R);
2395 : }
2396 :
2397 : GEN
2398 19418 : QXQ_norm(GEN A, GEN B)
2399 : {
2400 : GEN c, R, lB;
2401 19418 : long dA = degpol(A), dB = degpol(B);
2402 19418 : pari_sp av = avma;
2403 19418 : if (dA < 0) return gen_0;
2404 19418 : A = Q_primitive_part(A, &c);
2405 19418 : R = ZX_resultant(B, A);
2406 19418 : if (c) R = gmul(R, gpowgs(c, dB));
2407 19418 : lB = leading_coeff(B);
2408 19418 : if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
2409 19418 : return gc_upto(av, R);
2410 : }
2411 :
2412 : /* assume x has integral coefficients */
2413 : GEN
2414 1205967 : ZX_disc_all(GEN x, ulong bound)
2415 : {
2416 1205967 : pari_sp av = avma;
2417 1205967 : long s, d = degpol(x);
2418 : GEN l, R;
2419 :
2420 1205970 : if (d <= 1) return d == 1? gen_1: gen_0;
2421 1202677 : s = (d & 2) ? -1: 1;
2422 1202677 : l = leading_coeff(x);
2423 1202679 : if (!bound) bound = ZX_discbound(x);
2424 1202636 : R = ZX_resultant_all(x, NULL, NULL, bound);
2425 1202680 : if (is_pm1(l))
2426 1021105 : { if (signe(l) < 0) s = -s; }
2427 : else
2428 181578 : R = diviiexact(R,l);
2429 1202683 : if (s == -1) togglesign_safe(&R);
2430 1202684 : return gc_INT(av,R);
2431 : }
2432 :
2433 : GEN
2434 1143894 : ZX_disc(GEN x) { return ZX_disc_all(x,0); }
2435 :
2436 : static GEN
2437 11040 : ZXQX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, GEN T, ulong p)
2438 : {
2439 11040 : pari_sp av = avma;
2440 11040 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2441 : GEN H, dp;
2442 11041 : if (dropa && dropb) return pol0_Flx(T[1]); /* p | lc(A), p | lc(B) */
2443 11041 : H = FlxqX_saferesultant(a, b, T, p);
2444 11036 : if (!H) return NULL;
2445 11036 : if (dropa)
2446 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2447 0 : GEN c = gel(b,degB+2); /* lc(B) */
2448 0 : if (odd(degB)) c = Flx_neg(c, p);
2449 0 : c = Flxq_powu(c, dropa, T, p);
2450 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2451 : }
2452 11036 : else if (dropb)
2453 : { /* multiply by lc(A)^(deg B - deg b) */
2454 0 : GEN c = gel(a,degA+2); /* lc(A) */
2455 0 : c = Flxq_powu(c, dropb, T, p);
2456 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2457 : }
2458 11036 : dp = dB ? ZX_to_Flx(dB, p): pol1_Flx(T[1]);
2459 11038 : if (!Flx_equal1(dp))
2460 : {
2461 0 : GEN idp = Flxq_invsafe(dp, T, p);
2462 0 : if (!idp) return NULL;
2463 0 : H = Flxq_mul(H, Flxq_powu(idp, degA, T, p), T, p);
2464 : }
2465 11038 : return gc_leaf(av, H);
2466 : }
2467 :
2468 : /* If B=NULL, assume B=A' */
2469 : static GEN
2470 4925 : ZXQX_resultant_slice(GEN A, GEN B, GEN U, GEN dB, GEN P, GEN *mod)
2471 : {
2472 4925 : pari_sp av = avma;
2473 4925 : long degA, degB, i, n = lg(P)-1;
2474 : GEN H, T;
2475 4925 : long v = varn(U), redo = 0;
2476 :
2477 4925 : degA = degpol(A);
2478 4925 : degB = B? degpol(B): degA - 1;
2479 4925 : if (n == 1)
2480 : {
2481 3177 : ulong p = uel(P,1);
2482 3177 : GEN a = ZXX_to_FlxX(A, p, v), b = B? ZXX_to_FlxX(B, p, v): FlxX_deriv(a, p);
2483 3177 : GEN u = ZX_to_Flx(U, p);
2484 3177 : GEN Hp = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2485 3177 : if (!Hp) { set_avma(av); *mod = gen_1; return pol_0(v); }
2486 3177 : Hp = gc_upto(av, Flx_to_ZX(Hp)); *mod = utoipos(p); return Hp;
2487 : }
2488 1748 : T = ZV_producttree(P);
2489 1748 : A = ZXX_nv_mod_tree(A, P, T, v);
2490 1748 : if (B) B = ZXX_nv_mod_tree(B, P, T, v);
2491 1748 : U = ZX_nv_mod_tree(U, P, T);
2492 1748 : H = cgetg(n+1, t_VEC);
2493 9610 : for(i=1; i <= n; i++)
2494 : {
2495 7862 : ulong p = P[i];
2496 7862 : GEN a = gel(A,i), b = B? gel(B,i): FlxX_deriv(a, p), u = gel(U, i);
2497 7863 : GEN h = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2498 7861 : if (!h)
2499 : {
2500 0 : gel(H,i) = pol_0(v);
2501 0 : P[i] = 1; redo = 1;
2502 : }
2503 : else
2504 7861 : gel(H,i) = h;
2505 : }
2506 1748 : if (redo) T = ZV_producttree(P);
2507 1748 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2508 1748 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2509 : }
2510 :
2511 : GEN
2512 4925 : ZXQX_resultant_worker(GEN P, GEN A, GEN B, GEN T, GEN dB)
2513 : {
2514 4925 : GEN V = cgetg(3, t_VEC);
2515 4925 : if (isintzero(B)) B = NULL;
2516 4925 : if (!signe(dB)) dB = NULL;
2517 4925 : gel(V,1) = ZXQX_resultant_slice(A, B, T, dB, P, &gel(V,2));
2518 4925 : return V;
2519 : }
2520 :
2521 : static ulong
2522 4322 : ZXQX_resultant_bound_i(GEN nf, GEN A, GEN B, GEN (*f)(GEN,GEN,long))
2523 : {
2524 4322 : pari_sp av = avma;
2525 4322 : GEN r, M = nf_L2_bound(nf, NULL, &r);
2526 4322 : long v = nf_get_varn(nf), i, l = lg(r);
2527 4322 : GEN a = cgetg(l, t_COL);
2528 13383 : for (i = 1; i < l; i++)
2529 9061 : gel(a, i) = f(gsubst(A, v, gel(r,i)), gsubst(B, v, gel(r,i)), DEFAULTPREC);
2530 4322 : return gc_ulong(av, (ulong) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
2531 : }
2532 : static ulong
2533 4007 : ZXQX_resultant_bound(GEN nf, GEN A, GEN B)
2534 4007 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound); }
2535 :
2536 : static GEN
2537 56 : _ZXQ_powu(GEN x, ulong u, GEN T)
2538 56 : { return typ(x) == t_INT? powiu(x, u): ZXQ_powu(x, u, T); }
2539 :
2540 : /* Compute Res(A, B/dB) in Z[X]/T, assuming A,B in Z[X,Y], dB in Z or NULL (= 1)
2541 : * If B=NULL, take B = A' and assume deg A > 1 */
2542 : static GEN
2543 4004 : ZXQX_resultant_all(GEN A, GEN B, GEN T, GEN dB, ulong bound)
2544 : {
2545 4004 : pari_sp av = avma;
2546 : forprime_t S;
2547 : GEN H, worker;
2548 4004 : if (B)
2549 : {
2550 63 : long a = degpol(A), b = degpol(B);
2551 63 : if (a < 0 || b < 0) return gen_0;
2552 63 : if (!a) return _ZXQ_powu(gel(A,2), b, T);
2553 63 : if (!b) return _ZXQ_powu(gel(B,2), a, T);
2554 : } else
2555 3941 : if (!bound) B = RgX_deriv(A);
2556 4004 : if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, B);
2557 4004 : worker = snm_closure(is_entry("_ZXQX_resultant_worker"),
2558 : mkvec4(A, B? B: gen_0, T, dB? dB: gen_0));
2559 4004 : init_modular_big(&S);
2560 4004 : H = gen_crt("ZXQX_resultant_all", worker, &S, dB, bound, 0, NULL,
2561 : nxV_chinese_center, FpX_center);
2562 4004 : if (DEBUGLEVEL)
2563 0 : err_printf("ZXQX_resultant_all: a priori bound: %lu, a posteriori: %lu\n",
2564 : bound, expi(gsupnorm(H, DEFAULTPREC)));
2565 4004 : return gc_upto(av, H);
2566 : }
2567 :
2568 : GEN
2569 119 : nfX_resultant(GEN nf, GEN x, GEN y)
2570 : {
2571 119 : pari_sp av = avma;
2572 119 : GEN cx, cy, D, T = nf_get_pol(nf);
2573 119 : long dx = degpol(x), dy = degpol(y);
2574 119 : if (dx < 0 || dy < 0) return gen_0;
2575 119 : x = Q_primitive_part(x, &cx); if (cx) cx = gpowgs(cx, dy);
2576 119 : y = Q_primitive_part(y, &cy); if (cy) cy = gpowgs(cy, dx);
2577 119 : if (!dx) D = _ZXQ_powu(gel(x,2), dy, T);
2578 119 : else if (!dy) D = _ZXQ_powu(gel(y,2), dx, T);
2579 : else
2580 : {
2581 63 : ulong bound = ZXQX_resultant_bound(nf, x, y);
2582 63 : D = ZXQX_resultant_all(x, y, T, NULL, bound);
2583 : }
2584 119 : cx = mul_content(cx, cy); if (cx) D = gmul(D, cx);
2585 119 : return gc_upto(av, D);
2586 : }
2587 :
2588 : static GEN
2589 252 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2590 :
2591 : static GEN
2592 3941 : ZXQX_disc_all(GEN x, GEN T, ulong bound)
2593 : {
2594 3941 : pari_sp av = avma;
2595 3941 : long s, d = degpol(x), v = varn(T);
2596 : GEN l, R;
2597 :
2598 3941 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2599 3941 : s = (d & 2) ? -1: 1;
2600 3941 : l = leading_coeff(x);
2601 3941 : R = ZXQX_resultant_all(x, NULL, T, NULL, bound);
2602 3941 : if (!gequal1(l)) R = QXQ_div(R, to_ZX(l,v), T);
2603 3941 : if (s == -1) R = RgX_neg(R);
2604 3941 : return gc_upto(av, R);
2605 : }
2606 :
2607 : GEN
2608 7 : QX_disc(GEN x)
2609 : {
2610 7 : pari_sp av = avma;
2611 7 : GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
2612 7 : if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
2613 7 : return gc_upto(av, d);
2614 : }
2615 :
2616 : GEN
2617 4172 : nfX_disc(GEN nf, GEN x)
2618 : {
2619 4172 : pari_sp av = avma;
2620 4172 : GEN c, D, T = nf_get_pol(nf);
2621 : ulong bound;
2622 4172 : long d = degpol(x), v = varn(T);
2623 4172 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2624 3941 : x = Q_primitive_part(x, &c);
2625 3941 : bound = ZXQX_resultant_bound(nf, x, RgX_deriv(x));
2626 3941 : D = ZXQX_disc_all(x, T, bound);
2627 3941 : if (c) D = gmul(D, gpowgs(c, 2*d - 2));
2628 3941 : return gc_upto(av, D);
2629 : }
2630 :
2631 : GEN
2632 845493 : QXQ_mul(GEN x, GEN y, GEN T)
2633 : {
2634 845493 : GEN dx, nx = Q_primitive_part(x, &dx);
2635 845494 : GEN dy, ny = Q_primitive_part(y, &dy);
2636 845494 : GEN z = ZXQ_mul(nx, ny, T);
2637 845493 : if (dx || dy)
2638 : {
2639 842693 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
2640 842693 : if (!gequal1(d)) z = ZX_Q_mul(z, d);
2641 : }
2642 845493 : return z;
2643 : }
2644 :
2645 : GEN
2646 407060 : QXQ_sqr(GEN x, GEN T)
2647 : {
2648 407060 : GEN dx, nx = Q_primitive_part(x, &dx);
2649 407060 : GEN z = ZXQ_sqr(nx, T);
2650 407060 : if (dx)
2651 405324 : z = ZX_Q_mul(z, gsqr(dx));
2652 407060 : return z;
2653 : }
2654 :
2655 : static GEN
2656 212723 : QXQ_inv_slice(GEN A, GEN B, GEN P, GEN *mod)
2657 : {
2658 212723 : pari_sp av = avma;
2659 212723 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2660 : GEN H, T;
2661 212723 : if (n == 1)
2662 : {
2663 165637 : ulong p = uel(P,1);
2664 165637 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2665 165637 : GEN U = Flxq_invsafe(a, b, p);
2666 165637 : if (!U)
2667 : {
2668 24 : set_avma(av);
2669 24 : *mod = gen_1; return pol_0(v);
2670 : }
2671 165613 : H = gc_GEN(av, Flx_to_ZX(U));
2672 165613 : *mod = utoipos(p); return H;
2673 : }
2674 47086 : T = ZV_producttree(P);
2675 47086 : A = ZX_nv_mod_tree(A, P, T);
2676 47084 : B = ZX_nv_mod_tree(B, P, T);
2677 47083 : H = cgetg(n+1, t_VEC);
2678 238055 : for(i=1; i <= n; i++)
2679 : {
2680 190969 : ulong p = P[i];
2681 190969 : GEN a = gel(A,i), b = gel(B,i);
2682 190969 : GEN U = Flxq_invsafe(a, b, p);
2683 190966 : if (!U)
2684 : {
2685 602 : gel(H,i) = pol_0(v);
2686 601 : P[i] = 1; redo = 1;
2687 : }
2688 : else
2689 190364 : gel(H,i) = U;
2690 : }
2691 47086 : if (redo) T = ZV_producttree(P);
2692 47086 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2693 47085 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2694 : }
2695 :
2696 : GEN
2697 212723 : QXQ_inv_worker(GEN P, GEN A, GEN B)
2698 : {
2699 212723 : GEN V = cgetg(3, t_VEC);
2700 212723 : gel(V,1) = QXQ_inv_slice(A, B, P, &gel(V,2));
2701 212723 : return V;
2702 : }
2703 :
2704 : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
2705 : GEN
2706 145999 : QXQ_inv(GEN A, GEN B)
2707 : {
2708 : GEN D, Ap, Bp;
2709 : ulong pp;
2710 145999 : pari_sp av2, av = avma;
2711 : forprime_t S;
2712 145999 : GEN worker, U, H = NULL, mod = gen_1;
2713 : pari_timer ti;
2714 : long k, dA, dB;
2715 145999 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2716 : /* A a QX, B a ZX */
2717 145999 : A = Q_primitive_part(A, &D);
2718 145999 : dA = degpol(A); dB= degpol(B);
2719 : /* A, B in Z[X] */
2720 145999 : init_modular_small(&S);
2721 : do {
2722 145999 : pp = u_forprime_next(&S);
2723 145999 : Ap = ZX_to_Flx(A, pp);
2724 145999 : Bp = ZX_to_Flx(B, pp);
2725 145999 : } while (degpol(Ap) != dA || degpol(Bp) != dB);
2726 145999 : if (degpol(Flx_gcd(Ap, Bp, pp)) != 0 && degpol(ZX_gcd(A,B))!=0)
2727 14 : pari_err_INV("QXQ_inv",mkpolmod(A,B));
2728 145985 : worker = snm_closure(is_entry("_QXQ_inv_worker"), mkvec2(A, B));
2729 145985 : av2 = avma;
2730 145985 : for (k = 1; ;k *= 2)
2731 42510 : {
2732 : GEN res, b, N, den;
2733 188495 : gen_inccrt_i("QXQ_inv", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2734 : nxV_chinese_center, FpX_center);
2735 188495 : (void)gc_all(av2, 2, &H, &mod);
2736 188495 : b = sqrti(shifti(mod,-1));
2737 188495 : if (DEBUGLEVEL>5) timer_start(&ti);
2738 188495 : U = FpX_ratlift(H, mod, b, b, NULL);
2739 188495 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: ratlift");
2740 194200 : if (!U) continue;
2741 151690 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2742 151690 : res = Flx_rem(Flx_Fl_sub(Flx_mul(Ap, ZX_to_Flx(N,pp), pp),
2743 : umodiu(den, pp), pp), Bp, pp);
2744 151690 : if (degpol(res) >= 0) continue;
2745 145985 : res = ZX_Z_sub(ZX_mul(A, N), den);
2746 145985 : res = ZX_is_monic(B) ? ZX_rem(res, B): RgX_pseudorem(res, B);
2747 145985 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: final check");
2748 145985 : if (degpol(res)<0)
2749 : {
2750 145985 : if (D) U = RgX_Rg_div(U, D);
2751 145985 : return gc_GEN(av, U);
2752 : }
2753 : }
2754 : }
2755 :
2756 : static GEN
2757 121144 : QXQ_div_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
2758 : {
2759 121144 : pari_sp av = avma;
2760 121144 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2761 : GEN H, T;
2762 121144 : if (n == 1)
2763 : {
2764 44714 : ulong p = uel(P,1);
2765 44714 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), c = ZX_to_Flx(C, p);
2766 44714 : GEN bi = Flxq_invsafe(b, c, p), U;
2767 44714 : if (!bi)
2768 : {
2769 0 : set_avma(av);
2770 0 : *mod = gen_1; return pol_0(v);
2771 : }
2772 44714 : U = Flxq_mul(a, bi, c, p);
2773 44714 : H = gc_GEN(av, Flx_to_ZX(U));
2774 44714 : *mod = utoipos(p); return H;
2775 : }
2776 76430 : T = ZV_producttree(P);
2777 76430 : A = ZX_nv_mod_tree(A, P, T);
2778 76430 : B = ZX_nv_mod_tree(B, P, T);
2779 76430 : C = ZX_nv_mod_tree(C, P, T);
2780 76430 : H = cgetg(n+1, t_VEC);
2781 337524 : for(i=1; i <= n; i++)
2782 : {
2783 261094 : ulong p = P[i];
2784 261094 : GEN a = gel(A,i), b = gel(B,i), c = gel(C, i);
2785 261094 : GEN bi = Flxq_invsafe(b, c, p);
2786 261093 : if (!bi)
2787 : {
2788 4 : gel(H,i) = pol_0(v);
2789 4 : P[i] = 1; redo = 1;
2790 : }
2791 : else
2792 261089 : gel(H,i) = Flxq_mul(a, bi, c, p);
2793 : }
2794 76430 : if (redo) T = ZV_producttree(P);
2795 76430 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2796 76430 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2797 : }
2798 :
2799 : GEN
2800 121144 : QXQ_div_worker(GEN P, GEN A, GEN B, GEN C)
2801 : {
2802 121144 : GEN V = cgetg(3, t_VEC);
2803 121144 : gel(V,1) = QXQ_div_slice(A, B, C, P, &gel(V,2));
2804 121144 : return V;
2805 : }
2806 :
2807 : /* lift(Mod(A/B, C)). C a ZX, A, B a scalar or a QX */
2808 : GEN
2809 33207 : QXQ_div(GEN A, GEN B, GEN C)
2810 : {
2811 : GEN DA, DB, Ap, Bp, Cp;
2812 : ulong pp;
2813 33207 : pari_sp av2, av = avma;
2814 : forprime_t S;
2815 33207 : GEN worker, U, H = NULL, mod = gen_1;
2816 : pari_timer ti;
2817 : long k, dA, dB, dC;
2818 33207 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2819 : /* A a QX, B a ZX */
2820 33207 : A = Q_primitive_part(A, &DA);
2821 33207 : B = Q_primitive_part(B, &DB);
2822 33207 : dA = degpol(A); dB = degpol(B); dC = degpol(C);
2823 : /* A, B in Z[X] */
2824 33207 : init_modular_small(&S);
2825 : do {
2826 33207 : pp = u_forprime_next(&S);
2827 33207 : Ap = ZX_to_Flx(A, pp);
2828 33207 : Bp = ZX_to_Flx(B, pp);
2829 33207 : Cp = ZX_to_Flx(C, pp);
2830 33207 : } while (degpol(Ap) != dA || degpol(Bp) != dB || degpol(Cp) != dC);
2831 33207 : if (degpol(Flx_gcd(Bp, Cp, pp)) != 0 && degpol(ZX_gcd(B,C))!=0)
2832 0 : pari_err_INV("QXQ_div",mkpolmod(B,C));
2833 33207 : worker = snm_closure(is_entry("_QXQ_div_worker"), mkvec3(A, B, C));
2834 33207 : av2 = avma;
2835 33207 : for (k = 1; ;k *= 2)
2836 46749 : {
2837 : GEN res, b, N, den;
2838 79956 : gen_inccrt_i("QXQ_div", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2839 : nxV_chinese_center, FpX_center);
2840 79956 : (void)gc_all(av2, 2, &H, &mod);
2841 79956 : b = sqrti(shifti(mod,-1));
2842 79956 : if (DEBUGLEVEL>5) timer_start(&ti);
2843 79956 : U = FpX_ratlift(H, mod, b, b, NULL);
2844 79956 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: ratlift");
2845 90586 : if (!U) continue;
2846 43837 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2847 43837 : res = Flx_rem(Flx_sub(Flx_mul(Bp, ZX_to_Flx(N,pp), pp),
2848 : Flx_Fl_mul(Ap, umodiu(den, pp), pp), pp), Cp, pp);
2849 43837 : if (degpol(res) >= 0) continue;
2850 33207 : res = ZX_sub(ZX_mul(B, N), ZX_Z_mul(A,den));
2851 33207 : res = ZX_is_monic(C) ? ZX_rem(res, C): RgX_pseudorem(res, C);
2852 33207 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: final check");
2853 33207 : if (degpol(res)<0)
2854 : {
2855 33207 : if (DA && DB) U = RgX_Rg_mul(U, gdiv(DA,DB));
2856 28069 : else if (DA) U = RgX_Rg_mul(U, DA);
2857 15981 : else if (DB) U = RgX_Rg_div(U, DB);
2858 33207 : return gc_GEN(av, U);
2859 : }
2860 : }
2861 : }
2862 :
2863 : /************************************************************************
2864 : * *
2865 : * ZXQ_minpoly *
2866 : * *
2867 : ************************************************************************/
2868 :
2869 : static GEN
2870 3523 : ZXQ_minpoly_slice(GEN A, GEN B, GEN dA, long d, GEN P, GEN *mod)
2871 : {
2872 3523 : pari_sp av = avma;
2873 3523 : long i, n = lg(P)-1, v = evalvarn(varn(B)), redo = 0;
2874 : GEN H, T;
2875 3523 : if (n == 1)
2876 : {
2877 716 : ulong p = uel(P,1);
2878 716 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), Hp;
2879 716 : if (dA) a = Flx_Fl_mul(a, Fl_inv(umodiu(dA,p), p), p);
2880 716 : Hp = Flxq_minpoly(a, b, p);
2881 716 : if (degpol(Hp) != d) { p = 1; Hp = pol0_Flx(v); }
2882 716 : H = gc_upto(av, Flx_to_ZX(Hp));
2883 716 : *mod = utoipos(p); return H;
2884 : }
2885 2807 : T = ZV_producttree(P);
2886 2807 : A = ZX_nv_mod_tree(A, P, T);
2887 2807 : B = ZX_nv_mod_tree(B, P, T);
2888 2807 : H = cgetg(n+1, t_VEC);
2889 16838 : for(i=1; i <= n; i++)
2890 : {
2891 14031 : ulong p = P[i];
2892 14031 : GEN a = gel(A,i), b = gel(B,i), Hp;
2893 14031 : if (dA) a = Flx_Fl_mul(a, Fl_inv(umodiu(dA,p), p), p);
2894 14031 : Hp = Flxq_minpoly(a, b, p);
2895 14031 : if (degpol(Hp) != d) { P[i] = 1; redo = 1; Hp = pol0_Flx(v); }
2896 14031 : gel(H, i) = Hp;
2897 : }
2898 2807 : if (redo) T = ZV_producttree(P);
2899 2807 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2900 2807 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2901 : }
2902 :
2903 : GEN
2904 3523 : ZXQ_minpoly_worker(GEN P, GEN A, GEN B, GEN dA, long d)
2905 : {
2906 3523 : GEN V = cgetg(3, t_VEC);
2907 3523 : gel(V,1) = ZXQ_minpoly_slice(A, B, isintzero(dA) ? NULL: dA, d, P, &gel(V,2));
2908 3523 : return V;
2909 : }
2910 :
2911 : GEN
2912 1701 : ZXQ_minpoly(GEN A, GEN B, long d, ulong bound)
2913 : {
2914 1701 : pari_sp av = avma;
2915 : GEN worker, H, dA;
2916 : forprime_t S;
2917 1701 : A = Q_remove_denom(A, &dA);
2918 1701 : worker = strtoclosure("_ZXQ_minpoly_worker", 4, A, B, dA ? dA : gen_0, stoi(d));
2919 1701 : init_modular_big(&S);
2920 1701 : H = gen_crt("ZXQ_minpoly", worker, &S, dA, bound, 0, NULL,
2921 : nxV_chinese_center, FpX_center_i);
2922 1701 : return gc_GEN(av, H);
2923 : }
2924 :
2925 : /************************************************************************
2926 : * *
2927 : * ZXQXQ_minpoly *
2928 : * *
2929 : ************************************************************************/
2930 :
2931 : static GEN
2932 0 : ZXQXQ_minpoly_slice(GEN A, GEN B, GEN C, GEN dA, long d, GEN P, GEN *mod)
2933 : {
2934 0 : pari_sp av = avma;
2935 0 : long i, n = lg(P)-1, w = varn(C), dC = degpol(C), redo = 0;
2936 : GEN H, T;
2937 0 : if (n == 1)
2938 : {
2939 0 : ulong p = uel(P,1);
2940 0 : GEN a = ZXX_to_FlxX(A, p, w), b = ZXX_to_FlxX(B, p, w);
2941 0 : GEN c = ZX_to_Flx(C, p), Hp;
2942 0 : if (dA) a = FlxX_Fl_mul(a, Fl_inv(umodiu(dA, p), p), p);
2943 0 : Hp = FlxqXQ_minpoly(a, b, c, p);
2944 0 : if (degpol(Hp) != d)
2945 0 : { p = 1; H = gc_upto(av, zeromat(dC, d+1)); }
2946 : else
2947 0 : H = gc_upto(av, Flm_to_ZM(FlxX_to_Flm(Hp, dC)));
2948 0 : *mod = utoipos(p); return H;
2949 : }
2950 0 : T = ZV_producttree(P);
2951 0 : A = ZXX_nv_mod_tree(A, P, T, w);
2952 0 : B = ZXX_nv_mod_tree(B, P, T, w);
2953 0 : C = ZX_nv_mod_tree(C, P, T);
2954 0 : H = cgetg(n+1, t_VEC);
2955 0 : for (i = 1; i <= n; i++)
2956 : {
2957 0 : pari_sp av2 = avma;
2958 0 : ulong p = P[i];
2959 0 : GEN a = gel(A,i), b = gel(B,i), c = gel(C,i), Hp;
2960 0 : if (dA) a = FlxX_Fl_mul(a, Fl_inv(umodiu(dA, p), p), p);
2961 0 : Hp = FlxqXQ_minpoly(a, b, c, p);
2962 0 : if (degpol(Hp) != d) { P[i] = 1; redo = 1; gel(H, i) = zero_Flm(dC, d+1); }
2963 0 : else gel(H, i) = gc_upto(av2, FlxX_to_Flm(Hp, dC));
2964 : }
2965 0 : if (redo) T = ZV_producttree(P);
2966 0 : H = nmV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2967 0 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2968 : }
2969 :
2970 : GEN
2971 0 : ZXQXQ_minpoly_worker(GEN P, GEN A, GEN B, GEN C, GEN dA, long d)
2972 : {
2973 0 : GEN V = cgetg(3, t_VEC);
2974 0 : gel(V,1) = ZXQXQ_minpoly_slice(A, B, C, isintzero(dA) ? NULL: dA, d, P, &gel(V,2));
2975 0 : return V;
2976 : }
2977 :
2978 : GEN
2979 0 : ZXQXQ_minpoly(GEN A, GEN B, GEN T, long d, ulong bound)
2980 : {
2981 0 : pari_sp av = avma;
2982 : GEN worker, H, dA, dB, dAB;
2983 : forprime_t S;
2984 0 : A = Q_remove_denom(A, &dA);
2985 0 : B = Q_remove_denom(B, &dB);
2986 0 : dAB = dA ? (dB ? mulii(dA, dB): dA) : dB;
2987 0 : worker = strtoclosure("_ZXQXQ_minpoly_worker", 5, A, B, T, dA ? dA : gen_0, stoi(d));
2988 0 : init_modular_big(&S);
2989 0 : H = gen_crt("ZXQXQ_minpoly", worker, &S, dAB, bound, 0, NULL,
2990 : nmV_chinese_center, FpM_center);
2991 0 : return gc_GEN(av, RgM_to_RgXX(H, varn(A), varn(T)));
2992 : }
2993 :
2994 : /************************************************************************
2995 : * *
2996 : * ZX_ZXY_resultant *
2997 : * *
2998 : ************************************************************************/
2999 :
3000 : static GEN
3001 363671 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
3002 : long degA, long degB, long dres, long sX)
3003 : {
3004 363671 : pari_sp av = avma;
3005 363671 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
3006 363670 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3007 363670 : GEN Hp = Flx_FlxY_resultant_polint(a, b, p, pi, dres, sX);
3008 363668 : if (dropa && dropb)
3009 0 : Hp = zero_Flx(sX);
3010 : else {
3011 363668 : if (dropa)
3012 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
3013 0 : GEN c = gel(b,degB+2); /* lc(B) */
3014 0 : if (odd(degB)) c = Flx_neg(c, p);
3015 0 : if (!Flx_equal1(c)) {
3016 0 : c = Flx_powu_pre(c, dropa, p, pi);
3017 0 : if (!Flx_equal1(c)) Hp = Flx_mul_pre(Hp, c, p, pi);
3018 : }
3019 : }
3020 363668 : else if (dropb)
3021 : { /* multiply by lc(A)^(deg B - deg b) */
3022 0 : ulong c = uel(a, degA+2); /* lc(A) */
3023 0 : c = Fl_powu(c, dropb, p);
3024 0 : if (c != 1) Hp = Flx_Fl_mul_pre(Hp, c, p, pi);
3025 : }
3026 : }
3027 363668 : if (dp != 1) Hp = Flx_Fl_mul_pre(Hp, Fl_powu_pre(Fl_inv(dp,p), degA, p, pi), p, pi);
3028 363669 : return gc_leaf(av, Hp);
3029 : }
3030 :
3031 : static GEN
3032 124317 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
3033 : GEN P, GEN *mod, long sX, long vY)
3034 : {
3035 124317 : pari_sp av = avma;
3036 124317 : long i, n = lg(P)-1;
3037 : GEN H, T, D;
3038 124317 : if (n == 1)
3039 : {
3040 39995 : ulong p = uel(P,1);
3041 39995 : ulong dp = dB ? umodiu(dB, p): 1;
3042 39995 : GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
3043 39995 : GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
3044 39996 : H = gc_upto(av, Flx_to_ZX(Hp));
3045 39996 : *mod = utoipos(p); return H;
3046 : }
3047 84322 : T = ZV_producttree(P);
3048 84322 : A = ZX_nv_mod_tree(A, P, T);
3049 84322 : B = ZXX_nv_mod_tree(B, P, T, vY);
3050 84322 : D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
3051 84322 : H = cgetg(n+1, t_VEC);
3052 362677 : for(i=1; i <= n; i++)
3053 : {
3054 278355 : ulong p = P[i];
3055 278355 : GEN a = gel(A,i), b = gel(B,i);
3056 278355 : ulong dp = D ? uel(D, i): 1;
3057 278355 : gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
3058 : }
3059 84322 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
3060 84322 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3061 : }
3062 :
3063 : GEN
3064 124318 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
3065 : {
3066 124318 : GEN V = cgetg(3, t_VEC);
3067 124317 : if (isintzero(dB)) dB = NULL;
3068 124317 : gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
3069 124318 : return V;
3070 : }
3071 :
3072 : GEN
3073 78557 : ZX_ZXY_resultant(GEN A, GEN B)
3074 : {
3075 78557 : pari_sp av = avma;
3076 : forprime_t S;
3077 : ulong bound;
3078 78557 : long v = fetch_var_higher();
3079 78557 : long degA = degpol(A), degB, dres = degA * degpol(B);
3080 78557 : long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
3081 78557 : long sX = evalvarn(vX);
3082 : GEN worker, H, dB;
3083 78557 : B = Q_remove_denom(B, &dB);
3084 78558 : if (!dB) B = leafcopy(B);
3085 78558 : A = leafcopy(A); setvarn(A,v);
3086 78558 : B = swap_vars(B, vY, v); degB = degpol(B);
3087 78557 : bound = ZX_ZXY_ResBound(A, B, dB);
3088 78557 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
3089 157114 : worker = snm_closure(is_entry("_ZX_ZXY_resultant_worker"),
3090 78557 : mkvec4(A, B, dB? dB: gen_0,
3091 : mkvecsmall5(degA, degB, dres, sX, vY)));
3092 78558 : init_modular_big(&S);
3093 78558 : H = gen_crt("ZX_ZXY_resultant_all", worker, &S, dB, bound, 0, NULL,
3094 : nxV_chinese_center, FpX_center_i);
3095 78558 : setvarn(H, vX); (void)delete_var();
3096 78558 : return gc_GEN(av, H);
3097 : }
3098 :
3099 : static GEN
3100 1118089 : ZXX_resultant_prime(GEN a, GEN b, ulong p, long degA, long degB, long dres, long sX)
3101 : {
3102 1118089 : pari_sp av = avma;
3103 1118089 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
3104 1118089 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3105 1118089 : GEN Hp = FlxX_resultant_polint(a, b, p, pi, dres, sX);
3106 1118086 : if (dropa && dropb)
3107 0 : Hp = zero_Flx(sX);
3108 : else {
3109 1118086 : if (dropa)
3110 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
3111 0 : GEN c = gel(b,degB+2); /* lc(B) */
3112 0 : if (odd(degB)) c = Flx_neg(c, p);
3113 0 : if (!Flx_equal1(c)) {
3114 0 : c = Flx_powu_pre(c, dropa, p, pi);
3115 0 : if (!Flx_equal1(c)) Hp = Flx_mul_pre(Hp, c, p, pi);
3116 : }
3117 : }
3118 1118086 : else if (dropb)
3119 : { /* multiply by lc(A)^(deg B - deg b) */
3120 0 : ulong c = uel(a, degA+2); /* lc(A) */
3121 0 : c = Fl_powu(c, dropb, p);
3122 0 : if (c != 1) Hp = Flx_Fl_mul_pre(Hp, c, p, pi);
3123 : }
3124 : }
3125 1118086 : return gc_leaf(av, Hp);
3126 : }
3127 :
3128 : static GEN
3129 1015832 : ZXX_resultant_slice(GEN A, GEN B, long degA, long degB, long dres,
3130 : GEN P, GEN *mod, long sX, long vY)
3131 : {
3132 1015832 : pari_sp av = avma;
3133 1015832 : long i, n = lg(P)-1;
3134 : GEN H, T;
3135 1015832 : if (n == 1)
3136 : {
3137 931768 : ulong p = uel(P,1);
3138 931768 : GEN a = ZXX_to_FlxX(A, p, vY), b = ZXX_to_FlxX(B, p, vY);
3139 931768 : GEN Hp = ZXX_resultant_prime(a, b, p, degA, degB, dres, sX);
3140 931767 : H = gc_upto(av, Flx_to_ZX(Hp));
3141 931766 : *mod = utoipos(p); return H;
3142 : }
3143 84064 : T = ZV_producttree(P);
3144 84064 : A = ZXX_nv_mod_tree(A, P, T, vY);
3145 84064 : B = ZXX_nv_mod_tree(B, P, T, vY);
3146 84064 : H = cgetg(n+1, t_VEC);
3147 270386 : for(i=1; i <= n; i++)
3148 : {
3149 186322 : ulong p = P[i];
3150 186322 : GEN a = gel(A,i), b = gel(B,i);
3151 186322 : gel(H,i) = ZXX_resultant_prime(a, b, p, degA, degB, dres, sX);
3152 : }
3153 84064 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
3154 84064 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3155 : }
3156 :
3157 : GEN
3158 1015832 : ZXX_resultant_worker(GEN P, GEN A, GEN B, GEN v)
3159 : {
3160 1015832 : GEN V = cgetg(3, t_VEC);
3161 1015832 : gel(V,1) = ZXX_resultant_slice(A, B, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
3162 1015830 : return V;
3163 : }
3164 :
3165 : static int
3166 2021106 : ZXX_is_sparse(GEN x, long r)
3167 : {
3168 2021106 : long c = 0, i, l = lg(x), ly = 2, thr;
3169 8483020 : for (i = 2; i < l; i++)
3170 : {
3171 6461914 : GEN xi = gel(x,i);
3172 6461914 : if (typ(xi)==t_INT)
3173 5450318 : c += !!signe(xi);
3174 : else
3175 : {
3176 1011596 : long j, li = lg(xi);
3177 1011596 : ly = maxss(ly, li);
3178 3072469 : for (j = 2; j < li; j++)
3179 2060873 : c += !!signe(gel(xi,j));
3180 : }
3181 : }
3182 2021106 : thr = (l-2)*(ly-2)/r;
3183 2021106 : if (DEBUGLEVEL >= 5)
3184 0 : err_printf("ZXX_is_sparse: lx=%ld ly=%ld r=%ld c=%ld < thr=%ld \n",l-2, ly-2,r,c,thr);
3185 2021108 : return c < thr;
3186 : }
3187 :
3188 :
3189 : static GEN
3190 1010528 : ZXX_resultant_interp(GEN A, GEN B, long vX)
3191 : {
3192 1010528 : pari_sp av = avma;
3193 : forprime_t S;
3194 : ulong bound;
3195 1010528 : long degA = degpol(A), degB = degpol(B), dres = degA * RgXY_degreex(B) + RgXY_degreex(A)*degB;
3196 1010526 : long vY = varn(A), sX = evalvarn(vX);
3197 : GEN worker, H;
3198 1010526 : bound = ZXX_ResBound(A, B);
3199 1010528 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
3200 1010528 : worker = snm_closure(is_entry("_ZXX_resultant_worker"),
3201 : mkvec3(A, B, mkvecsmall5(degA, degB, dres, sX, vY)));
3202 1010531 : init_modular_big(&S);
3203 1010531 : H = gen_crt("ZXX_resultant", worker, &S, NULL, bound, 0, NULL, nxV_chinese_center, FpX_center_i);
3204 1010529 : return gc_GEN(av, H);
3205 : }
3206 :
3207 : GEN
3208 1010560 : ZXX_resultant(GEN A, GEN B, long vX)
3209 : {
3210 1010560 : if (!ZXX_is_sparse(A,6) && !ZXX_is_sparse(B,6))
3211 1010528 : return ZXX_resultant_interp(A, B, vX);
3212 35 : return RgX_resultant_all(A, B, NULL);
3213 : }
3214 :
3215 : static long
3216 40507 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
3217 : {
3218 40507 : pari_sp av = avma;
3219 40507 : long degA = degpol(A), degB, dres = degA*degpol(B0);
3220 40507 : long v = fetch_var_higher();
3221 40507 : long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
3222 40507 : long sX = evalvarn(vX);
3223 : GEN dB, B, a, b, Hp;
3224 : forprime_t S;
3225 :
3226 40507 : B0 = Q_remove_denom(B0, &dB);
3227 40508 : if (!dB) B0 = leafcopy(B0);
3228 40508 : A = leafcopy(A);
3229 40509 : B = B0;
3230 40509 : setvarn(A,v);
3231 45322 : INIT:
3232 45322 : if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
3233 45322 : B = swap_vars(B, vY, v);
3234 : /* B0(lambda v + x, v) */
3235 45320 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
3236 :
3237 45320 : degB = degpol(B);
3238 45320 : init_modular_big(&S);
3239 : while (1)
3240 0 : {
3241 45320 : ulong p = u_forprime_next(&S);
3242 45320 : ulong dp = dB ? umodiu(dB, p): 1;
3243 45320 : if (!dp) continue;
3244 45320 : a = ZX_to_Flx(A, p);
3245 45321 : b = ZXX_to_FlxX(B, p, v);
3246 45322 : Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
3247 45321 : if (degpol(Hp) != dres) continue;
3248 45321 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
3249 45321 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
3250 40509 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
3251 40509 : (void)delete_var(); return gc_long(av,lambda);
3252 : }
3253 : }
3254 :
3255 : GEN
3256 60552 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
3257 : {
3258 60552 : if (lambda)
3259 : {
3260 40507 : *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
3261 40509 : if (*lambda) B = RgX_Rg_translate(B, monomial(stoi(*lambda), 1, varn(A)));
3262 : }
3263 60553 : return ZX_ZXY_resultant(A,B);
3264 : }
3265 :
3266 : static GEN
3267 10483 : ZX_composedsum_slice(GEN A, GEN B, GEN P, GEN *mod)
3268 : {
3269 10483 : pari_sp av = avma;
3270 10483 : long i, n = lg(P)-1;
3271 : GEN H, T;
3272 10483 : if (n == 1)
3273 : {
3274 9977 : ulong p = uel(P,1);
3275 9977 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
3276 9976 : GEN Hp = Flx_composedsum(a, b, p);
3277 9978 : H = gc_upto(av, Flx_to_ZX(Hp));
3278 9978 : *mod = utoipos(p); return H;
3279 : }
3280 506 : T = ZV_producttree(P);
3281 506 : A = ZX_nv_mod_tree(A, P, T);
3282 506 : B = ZX_nv_mod_tree(B, P, T);
3283 506 : H = cgetg(n+1, t_VEC);
3284 4538 : for(i=1; i <= n; i++)
3285 : {
3286 4032 : ulong p = P[i];
3287 4032 : GEN a = gel(A,i), b = gel(B,i);
3288 4032 : gel(H,i) = Flx_composedsum(a, b, p);
3289 : }
3290 506 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
3291 506 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3292 : }
3293 :
3294 : GEN
3295 10483 : ZX_composedsum_worker(GEN P, GEN A, GEN B)
3296 : {
3297 10483 : GEN V = cgetg(3, t_VEC);
3298 10483 : gel(V,1) = ZX_composedsum_slice(A, B, P, &gel(V,2));
3299 10483 : return V;
3300 : }
3301 :
3302 : static GEN
3303 10217 : ZX_composedsum_i(GEN A, GEN B, GEN lead)
3304 : {
3305 10217 : pari_sp av = avma;
3306 : forprime_t S;
3307 : ulong bound;
3308 : GEN H, worker, mod;
3309 10217 : if (degpol(A) < degpol(B)) swap(A, B);
3310 10217 : if (!lead) lead = mulii(leading_coeff(A),leading_coeff(B));
3311 10217 : bound = ZX_ZXY_ResBound_1(A, B);
3312 10218 : worker = snm_closure(is_entry("_ZX_composedsum_worker"), mkvec2(A,B));
3313 10220 : init_modular_big(&S);
3314 10220 : H = gen_crt("ZX_composedsum", worker, &S, lead, bound, 0, &mod,
3315 : nxV_chinese_center, FpX_center);
3316 10220 : return gc_upto(av, H);
3317 : }
3318 :
3319 : static long
3320 9708 : ZX_compositum_lambda(GEN A, GEN B, GEN lead, long lambda)
3321 : {
3322 9708 : pari_sp av = avma;
3323 : forprime_t S;
3324 : ulong p;
3325 9708 : init_modular_big(&S);
3326 9708 : p = u_forprime_next(&S);
3327 : while (1)
3328 112 : {
3329 : GEN Hp, a;
3330 9820 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
3331 9820 : if (lead && dvdiu(lead,p)) { p = u_forprime_next(&S); continue; }
3332 9813 : a = ZX_to_Flx(ZX_rescale(A, stoi(-lambda)), p);
3333 9812 : Hp = Flx_composedsum(a, ZX_to_Flx(B, p), p);
3334 9806 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); continue; }
3335 9705 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
3336 9705 : return gc_long(av, lambda);
3337 : }
3338 : }
3339 :
3340 : GEN
3341 9708 : ZX_compositum(GEN A, GEN B, long *lambda)
3342 : {
3343 9708 : GEN lead = mulii(leading_coeff(A),leading_coeff(B));
3344 9708 : if (lambda)
3345 : {
3346 9708 : *lambda = ZX_compositum_lambda(A, B, lead, *lambda);
3347 9704 : A = ZX_rescale(A, stoi(-*lambda));
3348 : }
3349 9706 : return ZX_composedsum_i(A, B, lead);
3350 : }
3351 :
3352 : GEN
3353 511 : ZX_composedsum(GEN A, GEN B)
3354 511 : { return ZX_composedsum_i(A, B, NULL); }
3355 :
3356 : static GEN
3357 359 : ZXQX_composedsum_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
3358 : {
3359 359 : pari_sp av = avma;
3360 359 : long i, n = lg(P)-1, dC = degpol(C), v = varn(C);
3361 : GEN H, T;
3362 359 : if (n == 1)
3363 : {
3364 181 : ulong p = uel(P,1);
3365 181 : GEN a = ZXX_to_FlxX(A, p, v), b = ZXX_to_FlxX(B, p, v);
3366 181 : GEN c = ZX_to_Flx(C, p);
3367 181 : GEN Hp = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
3368 181 : H = gc_upto(av, Flm_to_ZM(Hp));
3369 181 : *mod = utoipos(p); return H;
3370 : }
3371 178 : T = ZV_producttree(P);
3372 178 : A = ZXX_nv_mod_tree(A, P, T, v);
3373 178 : B = ZXX_nv_mod_tree(B, P, T, v);
3374 178 : C = ZX_nv_mod_tree(C, P, T);
3375 178 : H = cgetg(n+1, t_VEC);
3376 660 : for(i=1; i <= n; i++)
3377 : {
3378 482 : ulong p = P[i];
3379 482 : GEN a = gel(A,i), b = gel(B,i), c = gel(C,i);
3380 482 : gel(H,i) = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
3381 : }
3382 178 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
3383 178 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3384 : }
3385 :
3386 : GEN
3387 359 : ZXQX_composedsum_worker(GEN P, GEN A, GEN B, GEN C)
3388 : {
3389 359 : GEN V = cgetg(3, t_VEC);
3390 359 : gel(V,1) = ZXQX_composedsum_slice(A, B, C, P, &gel(V,2));
3391 359 : return V;
3392 : }
3393 :
3394 : static GEN
3395 315 : ZXQX_composedsum(GEN A, GEN B, GEN T, ulong bound)
3396 : {
3397 315 : pari_sp av = avma;
3398 : forprime_t S;
3399 : GEN H, worker, mod;
3400 315 : GEN lead = mulii(Q_content(leading_coeff(A)), Q_content(leading_coeff(B)));
3401 315 : worker = snm_closure(is_entry("_ZXQX_composedsum_worker")
3402 : , mkvec3(A,B,T));
3403 315 : init_modular_big(&S);
3404 315 : H = gen_crt("ZXQX_composedsum", worker, &S, lead, bound, 0, &mod,
3405 : nmV_chinese_center, FpM_center);
3406 315 : if (DEBUGLEVEL > 4)
3407 0 : err_printf("nfcompositum: a priori bound: %lu, a posteriori: %lu\n",
3408 : bound, expi(gsupnorm(H, DEFAULTPREC)));
3409 315 : return gc_GEN(av, RgM_to_RgXX(H, varn(A), varn(T)));
3410 : }
3411 :
3412 : static long
3413 315 : ZXQX_composedsum_bound(GEN nf, GEN A, GEN B)
3414 315 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound_1); }
3415 :
3416 : GEN
3417 315 : nf_direct_compositum(GEN nf, GEN A, GEN B)
3418 : {
3419 315 : ulong bnd = ZXQX_composedsum_bound(nf, A, B);
3420 315 : return ZXQX_composedsum(A, B, nf_get_pol(nf), bnd);
3421 : }
3422 :
3423 : /************************************************************************
3424 : * *
3425 : * IRREDUCIBLE POLYNOMIAL / Fp *
3426 : * *
3427 : ************************************************************************/
3428 :
3429 : /* irreducible (unitary) polynomial of degree n over Fp */
3430 : GEN
3431 0 : ffinit_rand(GEN p,long n)
3432 : {
3433 0 : for(;;) {
3434 0 : pari_sp av = avma;
3435 0 : GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
3436 0 : if (FpX_is_irred(pol, p)) return pol;
3437 0 : set_avma(av);
3438 : }
3439 : }
3440 :
3441 : /* return an extension of degree 2^l of F_2, assume l > 0
3442 : * Not stack clean. */
3443 : static GEN
3444 602 : ffinit_Artin_Schreier_2(long l)
3445 : {
3446 : GEN Q, T, S;
3447 : long i, v;
3448 :
3449 602 : if (l == 1) return mkvecsmall4(0,1,1,1); /*x^2 + x + 1*/
3450 553 : v = fetch_var_higher();
3451 554 : S = mkvecsmall5(0, 0, 0, 1, 1); /* y(y^2 + y) */
3452 553 : Q = mkpoln(3, pol1_Flx(0), pol1_Flx(0), S); /* x^2 + x + y(y^2+y) */
3453 555 : setvarn(Q, v);
3454 :
3455 : /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
3456 555 : T = mkvecsmalln(6,evalvarn(v),1UL,1UL,0UL,0UL,1UL);
3457 : /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
3458 : * ==> x^2 + x + a(y) b irred. over K for any root b of Q
3459 : * ==> x^2 + x + (b^2+b)b */
3460 3067 : for (i=2; i<l; i++) T = Flx_FlxY_resultant(T, Q, 2); /* minpoly(b) / F2*/
3461 556 : (void)delete_var(); T[1] = 0; return T;
3462 : }
3463 :
3464 : /* return an extension of degree p^l of F_p, assume l > 0
3465 : * Not stack clean. */
3466 : GEN
3467 959 : ffinit_Artin_Schreier(ulong p, long l)
3468 : {
3469 : long i, v;
3470 : GEN Q, R, S, T, xp;
3471 959 : if (p==2) return ffinit_Artin_Schreier_2(l);
3472 357 : xp = polxn_Flx(p,0); /* x^p */
3473 357 : T = Flx_sub(xp, mkvecsmall3(0,1,1),p); /* x^p - x - 1 */
3474 357 : if (l == 1) return T;
3475 :
3476 7 : v = evalvarn(fetch_var_higher());
3477 7 : xp[1] = v;
3478 7 : R = Flx_sub(polxn_Flx(2*p-1,0), polxn_Flx(p,0),p);
3479 7 : S = Flx_sub(xp, polx_Flx(0), p);
3480 7 : Q = FlxY_Flx_sub(Flx_to_FlxX(S, v), R, p); /* x^p - x - (y^(2p-1)-y^p) */
3481 14 : for (i = 2; i <= l; ++i) T = Flx_FlxY_resultant(T, Q, p);
3482 7 : (void)delete_var(); T[1] = 0; return T;
3483 : }
3484 :
3485 : static long
3486 149483 : flinit_check(ulong p, long n, long l)
3487 : {
3488 : ulong q;
3489 149483 : if (!uisprime(n)) return 0;
3490 102377 : q = p % n; if (!q) return 0;
3491 99794 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3492 : }
3493 :
3494 : static GEN
3495 31902 : flinit(ulong p, long l)
3496 : {
3497 31902 : ulong n = 1+l;
3498 96544 : while (!flinit_check(p,n,l)) n += l;
3499 31902 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3500 31902 : return ZX_to_Flx(polsubcyclo(n,l,0), p);
3501 : }
3502 :
3503 : static GEN
3504 28980 : ffinit_fact_Flx(ulong p, long n)
3505 : {
3506 28980 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3507 28980 : long i, l = lg(Fm);
3508 28980 : P = cgetg(l, t_VEC);
3509 61843 : for (i = 1; i < l; i++)
3510 32861 : gel(P,i) = p==uel(Fp,i) ? ffinit_Artin_Schreier(p, Fe[i])
3511 32861 : : flinit(p, uel(Fm,i));
3512 28982 : return FlxV_composedsum(P, p);
3513 : }
3514 :
3515 : static GEN
3516 52946 : init_Flxq_i(ulong p, long n, long sv)
3517 : {
3518 : GEN P;
3519 52946 : if (!odd(p) && p != 2) pari_err_PRIME("ffinit", utoi(p));
3520 52939 : if (n == 1) return polx_Flx(sv);
3521 52939 : if (flinit_check(p, n+1, n))
3522 : {
3523 23959 : P = const_vecsmall(n+2,1);
3524 23959 : P[1] = sv; return P;
3525 : }
3526 28980 : P = ffinit_fact_Flx(p,n);
3527 28982 : P[1] = sv; return P;
3528 : }
3529 :
3530 : GEN
3531 0 : init_Flxq(ulong p, long n, long v)
3532 : {
3533 0 : pari_sp av = avma;
3534 0 : return gc_upto(av, init_Flxq_i(p, n, v));
3535 : }
3536 :
3537 : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
3538 : static long
3539 7332 : fpinit_check(GEN p, long n, long l)
3540 : {
3541 : ulong q;
3542 7332 : if (!uisprime(n)) return 0;
3543 4506 : q = umodiu(p,n); if (!q) return 0;
3544 4506 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3545 : }
3546 :
3547 : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
3548 : * Return an irreducible polynomial of degree l over F_p.
3549 : * Variant of Adleman and Lenstra "Finding irreducible polynomials over
3550 : * finite fields", ACM, 1986 (5) 350--355.
3551 : * Not stack clean */
3552 : static GEN
3553 1681 : fpinit(GEN p, long l)
3554 : {
3555 1681 : ulong n = 1+l;
3556 5335 : while (!fpinit_check(p,n,l)) n += l;
3557 1681 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3558 1681 : return FpX_red(polsubcyclo(n,l,0),p);
3559 : }
3560 :
3561 : static GEN
3562 1588 : ffinit_fact(GEN p, long n)
3563 : {
3564 1588 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3565 1588 : long i, l = lg(Fm);
3566 1588 : P = cgetg(l, t_VEC);
3567 3269 : for (i = 1; i < l; ++i)
3568 3362 : gel(P,i) = absequaliu(p, Fp[i]) ?
3569 0 : Flx_to_ZX(ffinit_Artin_Schreier(Fp[i], Fe[i]))
3570 1681 : : fpinit(p, Fm[i]);
3571 1588 : return FpXV_composedsum(P, p);
3572 : }
3573 :
3574 : static GEN
3575 55223 : init_Fq_i(GEN p, long n, long v)
3576 : {
3577 : GEN P;
3578 55223 : if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
3579 55223 : if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
3580 55223 : if (cmpiu(p, 2) < 0) pari_err_PRIME("ffinit",p);
3581 55216 : if (v < 0) v = 0;
3582 55216 : if (n == 1) return pol_x(v);
3583 54950 : if (lgefint(p) == 3)
3584 52946 : return Flx_to_ZX(init_Flxq_i(p[2], n, evalvarn(v)));
3585 2004 : if (!mpodd(p)) pari_err_PRIME("ffinit", p);
3586 1997 : if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
3587 1588 : P = ffinit_fact(p,n);
3588 1588 : setvarn(P, v); return P;
3589 : }
3590 : GEN
3591 54656 : init_Fq(GEN p, long n, long v)
3592 : {
3593 54656 : pari_sp av = avma;
3594 54656 : return gc_upto(av, init_Fq_i(p, n, v));
3595 : }
3596 : GEN
3597 567 : ffinit(GEN p, long n, long v)
3598 : {
3599 567 : pari_sp av = avma;
3600 567 : return gc_upto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
3601 : }
3602 :
3603 : GEN
3604 3178 : ffnbirred(GEN p, long n)
3605 : {
3606 3178 : pari_sp av = avma;
3607 3178 : GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
3608 3178 : long j, l = lg(D);
3609 6797 : for (j = 2; j < l; j++) /* skip d = 1 */
3610 : {
3611 3619 : long md = D[j]; /* mu(d) * d, d squarefree */
3612 3619 : GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
3613 3619 : s = md > 0? addii(s, pd): subii(s,pd);
3614 : }
3615 3178 : return gc_INT(av, diviuexact(s, n));
3616 : }
3617 :
3618 : GEN
3619 616 : ffsumnbirred(GEN p, long n)
3620 : {
3621 616 : pari_sp av = avma, av2;
3622 616 : GEN q, t = p, v = vecfactoru_i(1, n);
3623 : long i;
3624 616 : q = cgetg(n+1,t_VEC); gel(q,1) = p;
3625 1764 : for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
3626 616 : av2 = avma;
3627 1764 : for (i=2; i<=n; i++)
3628 : {
3629 1148 : GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
3630 1148 : long j, l = lg(D);
3631 2534 : for (j = 2; j < l; j++) /* skip 1 */
3632 : {
3633 1386 : long md = D[j];
3634 1386 : GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
3635 1386 : s = md > 0? addii(s, pd): subii(s, pd);
3636 : }
3637 1148 : t = gc_INT(av2, addii(t, diviuexact(s, i)));
3638 : }
3639 616 : return gc_INT(av, t);
3640 : }
3641 :
3642 : GEN
3643 140 : ffnbirred0(GEN p, long n, long flag)
3644 : {
3645 140 : if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
3646 140 : if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
3647 140 : switch(flag)
3648 : {
3649 70 : case 0: return ffnbirred(p, n);
3650 70 : case 1: return ffsumnbirred(p, n);
3651 : }
3652 0 : pari_err_FLAG("ffnbirred");
3653 : return NULL; /* LCOV_EXCL_LINE */
3654 : }
3655 :
3656 : static void
3657 2261 : checkmap(GEN m, const char *s)
3658 : {
3659 2261 : if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
3660 0 : pari_err_TYPE(s,m);
3661 2261 : }
3662 :
3663 : GEN
3664 189 : ffembed(GEN a, GEN b)
3665 : {
3666 189 : pari_sp av = avma;
3667 189 : GEN p, Ta, Tb, g, r = NULL;
3668 189 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
3669 189 : if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
3670 189 : p = FF_p_i(a); g = FF_gen(a);
3671 189 : if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
3672 189 : Ta = FF_mod(a);
3673 189 : Tb = FF_mod(b);
3674 189 : if (degpol(Tb)%degpol(Ta)!=0)
3675 7 : pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
3676 182 : r = gel(FFX_roots(Ta, b), 1);
3677 182 : return gc_GEN(av, mkvec2(g,r));
3678 : }
3679 :
3680 : GEN
3681 91 : ffextend(GEN a, GEN P, long v)
3682 : {
3683 91 : pari_sp av = avma;
3684 : long n;
3685 : GEN p, T, R, g, m;
3686 91 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
3687 91 : T = a; p = FF_p_i(a);
3688 91 : if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
3689 49 : if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
3690 49 : if (v < 0) v = varn(P);
3691 49 : n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
3692 49 : m = ffembed(a, g);
3693 49 : R = FFX_roots(ffmap(m, P),g);
3694 49 : return gc_GEN(av, mkvec2(gel(R,1), m));
3695 : }
3696 :
3697 : GEN
3698 42 : fffrobenius(GEN a, long n)
3699 : {
3700 42 : if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
3701 42 : retmkvec2(FF_gen(a), FF_Frobenius(a, n));
3702 : }
3703 :
3704 : GEN
3705 133 : ffinvmap(GEN m)
3706 : {
3707 133 : pari_sp av = avma;
3708 : long i, l;
3709 133 : GEN T, F, a, g, r, f = NULL;
3710 133 : checkmap(m, "ffinvmap");
3711 133 : a = gel(m,1); r = gel(m,2);
3712 133 : if (typ(r) != t_FFELT)
3713 7 : pari_err_TYPE("ffinvmap", m);
3714 126 : g = FF_gen(a);
3715 126 : T = FF_mod(r);
3716 126 : F = gel(FFX_factor(T, a), 1);
3717 126 : l = lg(F);
3718 490 : for(i=1; i<l; i++)
3719 : {
3720 490 : GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
3721 490 : if (degpol(s)==0 && gequal(constant_coeff(s),g)) { f = gel(F, i); break; }
3722 : }
3723 126 : if (f==NULL) pari_err_TYPE("ffinvmap", m);
3724 126 : if (degpol(f)==1) f = FF_neg_i(gel(f,2));
3725 126 : return gc_GEN(av, mkvec2(FF_gen(r),f));
3726 : }
3727 :
3728 : static GEN
3729 1260 : ffpartmapimage(const char *s, GEN r)
3730 : {
3731 1260 : GEN a = NULL, p = NULL;
3732 1260 : if (typ(r)==t_POL && degpol(r) >= 1
3733 1260 : && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
3734 0 : pari_err_TYPE(s, r);
3735 : return NULL; /* LCOV_EXCL_LINE */
3736 : }
3737 :
3738 : static GEN
3739 2709 : ffeltmap_i(GEN m, GEN x)
3740 : {
3741 2709 : GEN r = gel(m,2);
3742 2709 : if (!FF_samefield(x, gel(m,1)))
3743 84 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3744 2625 : if (typ(r)==t_FFELT)
3745 1659 : return FF_map(r, x);
3746 : else
3747 966 : return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
3748 : }
3749 :
3750 : static GEN
3751 4459 : ffmap_i(GEN m, GEN x)
3752 : {
3753 : GEN y;
3754 4459 : long i, lx, tx = typ(x);
3755 4459 : switch(tx)
3756 : {
3757 2541 : case t_FFELT:
3758 2541 : return ffeltmap_i(m, x);
3759 1267 : case t_POL: case t_RFRAC: case t_SER:
3760 : case t_VEC: case t_COL: case t_MAT:
3761 1267 : y = cgetg_copy(x, &lx);
3762 1988 : for (i = 1; i < lontyp[tx]; i++) y[i] = x[i];
3763 4564 : for (; i < lx; i++)
3764 : {
3765 3339 : GEN yi = ffmap_i(m, gel(x,i));
3766 3297 : if (!yi) return NULL;
3767 3297 : gel(y,i) = yi;
3768 : }
3769 1225 : return y;
3770 : }
3771 651 : return gcopy(x);
3772 : }
3773 :
3774 : GEN
3775 1036 : ffmap(GEN m, GEN x)
3776 : {
3777 1036 : pari_sp ltop = avma;
3778 : GEN y;
3779 1036 : checkmap(m, "ffmap");
3780 1036 : y = ffmap_i(m, x);
3781 1036 : if (y) return y;
3782 42 : retgc_const(ltop, cgetg(1, t_VEC));
3783 : }
3784 :
3785 : static GEN
3786 252 : ffeltmaprel_i(GEN m, GEN x)
3787 : {
3788 252 : GEN g = gel(m,1), r = gel(m,2);
3789 252 : if (!FF_samefield(x, g))
3790 0 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3791 252 : if (typ(r)==t_FFELT)
3792 84 : retmkpolmod(FF_map(r, x), pol_x(FF_var(g)));
3793 : else
3794 168 : retmkpolmod(FFX_preimagerel(x, r, ffpartmapimage("ffmap", r)), gcopy(r));
3795 : }
3796 :
3797 : static GEN
3798 252 : ffmaprel_i(GEN m, GEN x)
3799 : {
3800 252 : switch(typ(x))
3801 : {
3802 252 : case t_FFELT:
3803 252 : return ffeltmaprel_i(m, x);
3804 0 : case t_POL: pari_APPLY_pol_normalized(ffmaprel_i(m, gel(x,i)));
3805 0 : case t_SER: pari_APPLY_ser_normalized(ffmaprel_i(m, gel(x,i)));
3806 0 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
3807 0 : pari_APPLY_same(ffmaprel_i(m, gel(x,i)));
3808 : }
3809 0 : return gcopy(x);
3810 : }
3811 : GEN
3812 252 : ffmaprel(GEN m, GEN x) { checkmap(m, "ffmaprel"); return ffmaprel_i(m, x); }
3813 :
3814 : static void
3815 84 : err_compo(GEN m, GEN n)
3816 84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
3817 :
3818 : GEN
3819 420 : ffcompomap(GEN m, GEN n)
3820 : {
3821 420 : pari_sp av = avma;
3822 420 : GEN g = gel(n,1), r, m2, n2;
3823 420 : checkmap(m, "ffcompomap");
3824 420 : checkmap(n, "ffcompomap");
3825 420 : m2 = gel(m,2); n2 = gel(n,2);
3826 420 : switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
3827 : {
3828 84 : case 0:
3829 84 : if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
3830 42 : r = FF_map(gel(m,2), n2);
3831 42 : break;
3832 84 : case 2:
3833 84 : r = ffmap_i(m, n2);
3834 42 : if (lg(r) == 1) err_compo(m,n);
3835 42 : break;
3836 168 : case 1:
3837 168 : r = ffeltmap_i(m, n2);
3838 126 : if (!r)
3839 : {
3840 : GEN a, A, R, M;
3841 : long dm, dn;
3842 42 : a = ffpartmapimage("ffcompomap",m2);
3843 42 : A = FF_to_FpXQ_i(FF_neg(n2));
3844 42 : setvarn(A, 1);
3845 42 : R = deg1pol_shallow(gen_1, A, 0);
3846 42 : setvarn(R, 0);
3847 42 : M = gcopy(m2);
3848 42 : setvarn(M, 1);
3849 42 : r = polresultant0(R, M, 1, 0);
3850 42 : dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
3851 42 : if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
3852 42 : setvarn(r, varn(FF_mod(g)));
3853 : }
3854 126 : break;
3855 84 : case 3:
3856 : {
3857 : GEN M, R, T, p, a;
3858 84 : a = ffpartmapimage("ffcompomap",n2);
3859 84 : if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
3860 42 : p = FF_p_i(gel(n,1));
3861 42 : T = FF_mod(gel(n,1));
3862 42 : setvarn(T, 1);
3863 42 : R = RgX_to_FpXQX(n2,T,p);
3864 42 : setvarn(R, 0);
3865 42 : M = gcopy(m2);
3866 42 : setvarn(M, 1);
3867 42 : r = polresultant0(R, M, 1, 0);
3868 42 : setvarn(r, varn(n2));
3869 : }
3870 : }
3871 252 : return gc_GEN(av, mkvec2(g,r));
3872 : }
|
