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 6594 : char_update_int(struct charact *S, GEN n)
43 : {
44 6594 : if (S->isprime)
45 : {
46 7 : if (dvdii(n, S->q)) return;
47 7 : pari_err_MODULUS("characteristic", S->q, n);
48 : }
49 6587 : S->q = gcdii(S->q, n);
50 : }
51 : static void
52 176946 : charact(struct charact *S, GEN x)
53 : {
54 176946 : const long tx = typ(x);
55 : long i, l;
56 176946 : switch(tx)
57 : {
58 5145 : 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 25711 : 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 25711 : l = lg(x);
64 175021 : for (i=lontyp[tx]; i < l; i++) charact(S,gel(x,i));
65 25697 : break;
66 7 : case t_LIST:
67 7 : x = list_data(x);
68 7 : if (x) charact(S, x);
69 7 : break;
70 : }
71 176918 : }
72 : static void
73 4739 : charact_res(struct charact *S, GEN x)
74 : {
75 4739 : const long tx = typ(x);
76 : long i, l;
77 4739 : switch(tx)
78 : {
79 1449 : case t_INTMOD:char_update_int(S, gel(x,1)); break;
80 0 : case t_FFELT: char_update_prime(S, gel(x,4)); break;
81 91 : case t_PADIC: char_update_prime(S, padic_p(x)); break;
82 1722 : 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 1722 : l = lg(x);
86 6132 : for (i=lontyp[tx]; i < l; i++) charact_res(S,gel(x,i));
87 1722 : 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 4739 : }
94 : GEN
95 27622 : characteristic(GEN x)
96 : {
97 : struct charact S;
98 27622 : S.q = gen_0; S.isprime = 0;
99 27622 : charact(&S, x); return S.q;
100 : }
101 : GEN
102 329 : residual_characteristic(GEN x)
103 : {
104 : struct charact S;
105 329 : S.q = gen_0; S.isprime = 0;
106 329 : charact_res(&S, x); return S.q;
107 : }
108 :
109 : int
110 71065574 : Rg_is_Fp(GEN x, GEN *pp)
111 : {
112 : GEN mod;
113 71065574 : switch(typ(x))
114 : {
115 2482676 : case t_INTMOD:
116 2482676 : mod = gel(x,1);
117 2482676 : if (!*pp) *pp = mod;
118 2341976 : 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 2482676 : return 1;
124 57168718 : case t_INT:
125 57168718 : return 1;
126 11414180 : default: return 0;
127 : }
128 : }
129 :
130 : int
131 28194137 : RgX_is_FpX(GEN x, GEN *pp)
132 : {
133 28194137 : long i, lx = lg(x);
134 87819392 : for (i=2; i<lx; i++)
135 71039431 : if (!Rg_is_Fp(gel(x, i), pp))
136 11414182 : return 0;
137 16779961 : 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 60802 : Rg_is_FpXQ(GEN x, GEN *pT, GEN *pp)
160 : {
161 : GEN pol, mod, p;
162 60802 : switch(typ(x))
163 : {
164 26131 : case t_INTMOD:
165 26131 : 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 4578 : if (!RgX_is_FpX(pol, pp)) return 0;
186 : }
187 7 : 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 3381 : RgX_is_FpXQX(GEN x, GEN *pT, GEN *pp)
202 : {
203 3381 : long i, lx = lg(x);
204 63427 : for (i = 2; i < lx; i++)
205 60144 : if (!Rg_is_FpXQ(gel(x,i), pT, pp)) return 0;
206 3283 : return 1;
207 : }
208 :
209 : /************************************************************************
210 : ** **
211 : ** Ring conversion **
212 : ** **
213 : ************************************************************************/
214 :
215 : /* p > 0 a t_INT, return lift(x * Mod(1,p)).
216 : * If x is an INTMOD, assume modulus is a multiple of p. */
217 : GEN
218 52334333 : Rg_to_Fp(GEN x, GEN p)
219 : {
220 52334333 : if (lgefint(p) == 3) return utoi(Rg_to_Fl(x, uel(p,2)));
221 4555476 : switch(typ(x))
222 : {
223 289093 : case t_INT: return modii(x, p);
224 18790 : case t_FRAC: {
225 18790 : pari_sp av = avma;
226 18790 : GEN z = modii(gel(x,1), p);
227 18790 : if (z == gen_0) return gen_0;
228 18785 : return gc_INT(av, remii(mulii(z, Fp_inv(gel(x,2), p)), p));
229 : }
230 70 : case t_PADIC: return padic_to_Fp(x, p);
231 4247525 : case t_INTMOD: {
232 4247525 : GEN q = gel(x,1), a = gel(x,2);
233 4247525 : if (equalii(q, p)) return icopy(a);
234 14 : if (!dvdii(q,p)) pari_err_MODULUS("Rg_to_Fp", q, p);
235 0 : return remii(a, p);
236 : }
237 0 : default: pari_err_TYPE("Rg_to_Fp",x);
238 : return NULL; /* LCOV_EXCL_LINE */
239 : }
240 : }
241 : /* If x is a POLMOD, assume modulus is a multiple of T. */
242 : GEN
243 1291986 : Rg_to_FpXQ(GEN x, GEN T, GEN p)
244 : {
245 1291986 : long ta, tx = typ(x), v = get_FpX_var(T);
246 : GEN a, b;
247 1291986 : if (is_const_t(tx))
248 : {
249 59182 : if (tx == t_FFELT)
250 : {
251 17355 : GEN z = FF_to_FpXQ(x);
252 17355 : setvarn(z, v);
253 17355 : return z;
254 : }
255 41827 : return scalar_ZX(degpol(T)? Rg_to_Fp(x, p): gen_0, v);
256 : }
257 1232804 : switch(tx)
258 : {
259 1230697 : case t_POLMOD:
260 1230697 : b = gel(x,1);
261 1230697 : a = gel(x,2); ta = typ(a);
262 1230697 : if (is_const_t(ta))
263 3885 : return scalar_ZX(degpol(T)? Rg_to_Fp(a, p): gen_0, v);
264 1226812 : b = RgX_to_FpX(b, p); if (varn(b) != v) break;
265 1226812 : a = RgX_to_FpX(a, p);
266 1226812 : if (ZX_equal(b,get_FpX_mod(T)) || signe(FpX_rem(b,T,p))==0)
267 1226812 : return FpX_rem(a, T, p);
268 0 : break;
269 2107 : case t_POL:
270 2107 : if (varn(x) != v) break;
271 2100 : return FpX_rem(RgX_to_FpX(x,p), T, p);
272 0 : case t_RFRAC:
273 0 : a = Rg_to_FpXQ(gel(x,1), T,p);
274 0 : b = Rg_to_FpXQ(gel(x,2), T,p);
275 0 : return FpXQ_div(a,b, T,p);
276 : }
277 7 : pari_err_TYPE("Rg_to_FpXQ",x);
278 : return NULL; /* LCOV_EXCL_LINE */
279 : }
280 : GEN
281 3335720 : RgX_to_FpX(GEN x, GEN p)
282 : {
283 : long i, l;
284 3335720 : GEN z = cgetg_copy(x, &l); z[1] = x[1];
285 14765324 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
286 3335720 : return FpX_renormalize(z, l);
287 : }
288 :
289 : GEN
290 140 : RgV_to_FpV(GEN x, GEN p)
291 483 : { pari_APPLY_type(t_VEC, Rg_to_Fp(gel(x,i), p)) }
292 :
293 : GEN
294 1751090 : RgC_to_FpC(GEN x, GEN p)
295 28499485 : { pari_APPLY_type(t_COL, Rg_to_Fp(gel(x,i), p)) }
296 :
297 : GEN
298 222349 : RgM_to_FpM(GEN x, GEN p)
299 1973397 : { pari_APPLY_same(RgC_to_FpC(gel(x,i), p)) }
300 :
301 : GEN
302 369001 : RgV_to_Flv(GEN x, ulong p)
303 1676894 : { pari_APPLY_ulong(Rg_to_Fl(gel(x,i), p)) }
304 :
305 : GEN
306 118296 : RgM_to_Flm(GEN x, ulong p)
307 422998 : { pari_APPLY_same(RgV_to_Flv(gel(x,i), p)) }
308 :
309 : GEN
310 5105 : RgX_to_FpXQX(GEN x, GEN T, GEN p)
311 : {
312 5105 : long i, l = lg(x);
313 5105 : GEN z = cgetg(l, t_POL); z[1] = x[1];
314 43394 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T,p);
315 5105 : return FpXQX_renormalize(z, l);
316 : }
317 : GEN
318 49186 : RgX_to_FqX(GEN x, GEN T, GEN p)
319 : {
320 49186 : long i, l = lg(x);
321 49186 : GEN z = cgetg(l, t_POL); z[1] = x[1];
322 49186 : if (T)
323 10990 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
324 : else
325 791394 : for (i = 2; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
326 49186 : return FpXQX_renormalize(z, l);
327 : }
328 :
329 : GEN
330 219128 : RgC_to_FqC(GEN x, GEN T, GEN p)
331 : {
332 219128 : long i, l = lg(x);
333 219128 : GEN z = cgetg(l, t_COL);
334 219128 : if (T)
335 1430310 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_FpXQ(gel(x,i), T, p);
336 : else
337 0 : for (i = 1; i < l; i++) gel(z,i) = Rg_to_Fp(gel(x,i), p);
338 219128 : return z;
339 : }
340 :
341 : GEN
342 52325 : RgM_to_FqM(GEN x, GEN T, GEN p)
343 271418 : { pari_APPLY_same(RgC_to_FqC(gel(x, i), T, p)) }
344 :
345 : /* lg(V) > 1 */
346 : GEN
347 851487 : FpXV_FpC_mul(GEN V, GEN W, GEN p)
348 : {
349 851487 : pari_sp av = avma;
350 851487 : long i, l = lg(V);
351 851487 : GEN z = ZX_Z_mul(gel(V,1),gel(W,1));
352 4201029 : for(i=2; i<l; i++)
353 : {
354 3349542 : z = ZX_add(z, ZX_Z_mul(gel(V,i),gel(W,i)));
355 3349542 : if ((i & 7) == 0) z = gc_upto(av, z);
356 : }
357 851487 : return gc_upto(av, FpX_red(z,p));
358 : }
359 :
360 : GEN
361 55832 : FqX_Fq_add(GEN y, GEN x, GEN T, GEN p)
362 : {
363 55832 : long i, lz = lg(y);
364 : GEN z;
365 55832 : if (!T) return FpX_Fp_add(y, x, p);
366 8666 : if (lz == 2) return scalarpol(x, varn(y));
367 7952 : z = cgetg(lz,t_POL); z[1] = y[1];
368 7952 : gel(z,2) = Fq_add(gel(y,2),x, T, p);
369 7952 : if (lz == 3) z = FpXX_renormalize(z,lz);
370 : else
371 33145 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
372 7952 : return z;
373 : }
374 :
375 : GEN
376 1059 : FqX_Fq_sub(GEN y, GEN x, GEN T, GEN p)
377 : {
378 1059 : long i, lz = lg(y);
379 : GEN z;
380 1059 : if (!T) return FpX_Fp_sub(y, x, p);
381 1059 : if (lz == 2) return scalarpol(x, varn(y));
382 1059 : z = cgetg(lz,t_POL); z[1] = y[1];
383 1059 : gel(z,2) = Fq_sub(gel(y,2), x, T, p);
384 1059 : if (lz == 3) z = FpXX_renormalize(z,lz);
385 : else
386 10278 : for(i=3;i<lz;i++) gel(z,i) = gcopy(gel(y,i));
387 1059 : return z;
388 : }
389 :
390 : GEN
391 149052 : FqX_Fq_mul_to_monic(GEN P, GEN U, GEN T, GEN p)
392 : {
393 : long i, lP;
394 149052 : GEN res = cgetg_copy(P, &lP); res[1] = P[1];
395 918799 : for(i=2; i<lP-1; i++) gel(res,i) = Fq_mul(U,gel(P,i), T,p);
396 149052 : gel(res,lP-1) = gen_1; return res;
397 : }
398 :
399 : GEN
400 38188 : FpXQX_normalize(GEN z, GEN T, GEN p)
401 : {
402 : GEN lc;
403 38188 : if (lg(z) == 2) return z;
404 38174 : lc = leading_coeff(z);
405 38174 : if (typ(lc) == t_POL)
406 : {
407 2167 : if (lg(lc) > 3) /* nonconstant */
408 1902 : return FqX_Fq_mul_to_monic(z, Fq_inv(lc,T,p), T,p);
409 : /* constant */
410 265 : lc = gel(lc,2);
411 265 : z = shallowcopy(z);
412 265 : gel(z, lg(z)-1) = lc;
413 : }
414 : /* lc a t_INT */
415 36272 : if (equali1(lc)) return z;
416 87 : return FqX_Fq_mul_to_monic(z, Fp_inv(lc,p), T,p);
417 : }
418 :
419 : GEN
420 390934 : FqX_eval(GEN x, GEN y, GEN T, GEN p)
421 : {
422 : pari_sp av;
423 : GEN p1, r;
424 390934 : long j, i=lg(x)-1;
425 390934 : if (i<=2)
426 45107 : return (i==2)? Fq_red(gel(x,2), T, p): gen_0;
427 345827 : av=avma; p1=gel(x,i);
428 : /* specific attention to sparse polynomials (see poleval)*/
429 : /*You've guessed it! It's a copy-paste(tm)*/
430 1150444 : for (i--; i>=2; i=j-1)
431 : {
432 805302 : for (j=i; !signe(gel(x,j)); j--)
433 686 : if (j==2)
434 : {
435 490 : if (i!=j) y = Fq_pow(y,utoipos(i-j+1), T, p);
436 490 : return gc_upto(av, Fq_mul(p1,y, T, p));
437 : }
438 804616 : r = (i==j)? y: Fq_pow(y, utoipos(i-j+1), T, p);
439 804616 : p1 = Fq_add(Fq_mul(p1,r,T,p), gel(x,j), T, p);
440 : }
441 345338 : return gc_upto(av, p1);
442 : }
443 :
444 : GEN
445 97321 : FqXY_evalx(GEN Q, GEN x, GEN T, GEN p)
446 : {
447 97321 : long i, lb = lg(Q);
448 : GEN z;
449 97321 : if (!T) return FpXY_evalx(Q, x, p);
450 86961 : z = cgetg(lb, t_POL); z[1] = Q[1];
451 462945 : for (i=2; i<lb; i++)
452 : {
453 375984 : GEN q = gel(Q,i);
454 375984 : gel(z,i) = typ(q) == t_INT? modii(q,p): FqX_eval(q, x, T, p);
455 : }
456 86961 : return FpXQX_renormalize(z, lb);
457 : }
458 :
459 : /* Q an FpXY, evaluate at (X,Y) = (x,y) */
460 : GEN
461 14623 : FqXY_eval(GEN Q, GEN y, GEN x, GEN T, GEN p)
462 : {
463 14623 : pari_sp av = avma;
464 14623 : if (!T) return FpXY_eval(Q, y, x, p);
465 966 : return gc_upto(av, FqX_eval(FqXY_evalx(Q, x, T, p), y, T, p));
466 : }
467 :
468 : /* a X^d */
469 : GEN
470 13640839 : monomial(GEN a, long d, long v)
471 : {
472 : long i, n;
473 : GEN P;
474 13640839 : if (d < 0) {
475 14 : if (isrationalzero(a)) return pol_0(v);
476 14 : retmkrfrac(a, pol_xn(-d, v));
477 : }
478 13640825 : if (gequal0(a))
479 : {
480 93989 : if (isexactzero(a)) return scalarpol_shallow(a,v);
481 0 : n = d+2; P = cgetg(n+1, t_POL);
482 0 : P[1] = evalsigne(0) | evalvarn(v);
483 : }
484 : else
485 : {
486 13546835 : n = d+2; P = cgetg(n+1, t_POL);
487 13546831 : P[1] = evalsigne(1) | evalvarn(v);
488 : }
489 32924018 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
490 13546831 : gel(P,i) = a; return P;
491 : }
492 : GEN
493 157969 : monomialcopy(GEN a, long d, long v)
494 : {
495 : long i, n;
496 : GEN P;
497 157969 : if (d < 0) {
498 14 : if (isrationalzero(a)) return pol_0(v);
499 14 : retmkrfrac(gcopy(a), pol_xn(-d, v));
500 : }
501 157955 : if (gequal0(a))
502 : {
503 14 : if (isexactzero(a)) return scalarpol(a,v);
504 7 : n = d+2; P = cgetg(n+1, t_POL);
505 7 : P[1] = evalsigne(0) | evalvarn(v);
506 : }
507 : else
508 : {
509 157941 : n = d+2; P = cgetg(n+1, t_POL);
510 157941 : P[1] = evalsigne(1) | evalvarn(v);
511 : }
512 314678 : for (i = 2; i < n; i++) gel(P,i) = gen_0;
513 157948 : gel(P,i) = gcopy(a); return P;
514 : }
515 : GEN
516 325951 : pol_x_powers(long N, long v)
517 : {
518 325951 : GEN L = cgetg(N+1,t_VEC);
519 : long i;
520 789057 : for (i=1; i<=N; i++) gel(L,i) = pol_xn(i-1, v);
521 325955 : return L;
522 : }
523 :
524 : GEN
525 0 : FqXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
526 : {
527 0 : return T ? FpXQXQ_powers(x, l, S, T, p): FpXQ_powers(x, l, S, p);
528 : }
529 :
530 : GEN
531 0 : FqXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
532 : {
533 0 : return T ? FpXQXQ_matrix_pow(y, n, m, S, T, p): FpXQ_matrix_pow(y, n, m, S, p);
534 : }
535 :
536 : /*******************************************************************/
537 : /* */
538 : /* Fq */
539 : /* */
540 : /*******************************************************************/
541 :
542 : GEN
543 11592898 : Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
544 : {
545 : (void)T;
546 11592898 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
547 : {
548 1143687 : case 0: return Fp_add(x,y,p);
549 748136 : case 1: return FpX_Fp_add(x,y,p);
550 91991 : case 2: return FpX_Fp_add(y,x,p);
551 9609084 : case 3: return FpX_add(x,y,p);
552 : }
553 : return NULL;/*LCOV_EXCL_LINE*/
554 : }
555 :
556 : GEN
557 8563061 : Fq_sub(GEN x, GEN y, GEN T/*unused*/, GEN p)
558 : {
559 : (void)T;
560 8563061 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
561 : {
562 256368 : case 0: return Fp_sub(x,y,p);
563 24540 : case 1: return FpX_Fp_sub(x,y,p);
564 23896 : case 2: return Fp_FpX_sub(x,y,p);
565 8258257 : case 3: return FpX_sub(x,y,p);
566 : }
567 : return NULL;/*LCOV_EXCL_LINE*/
568 : }
569 :
570 : GEN
571 1080519 : Fq_neg(GEN x, GEN T/*unused*/, GEN p)
572 : {
573 : (void)T;
574 1080519 : return (typ(x)==t_POL)? FpX_neg(x,p): Fp_neg(x,p);
575 : }
576 :
577 : GEN
578 81354 : Fq_halve(GEN x, GEN T/*unused*/, GEN p)
579 : {
580 : (void)T;
581 81354 : return (typ(x)==t_POL)? FpX_halve(x,p): Fp_halve(x,p);
582 : }
583 :
584 : /* If T==NULL do not reduce*/
585 : GEN
586 8608652 : Fq_mul(GEN x, GEN y, GEN T, GEN p)
587 : {
588 8608652 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
589 : {
590 1037917 : case 0: return Fp_mul(x,y,p);
591 128565 : case 1: return FpX_Fp_mul(x,y,p);
592 394685 : case 2: return FpX_Fp_mul(y,x,p);
593 7047488 : case 3: if (T) return FpXQ_mul(x,y,T,p);
594 4478770 : else return FpX_mul(x,y,p);
595 : }
596 : return NULL;/*LCOV_EXCL_LINE*/
597 : }
598 :
599 : /* If T==NULL do not reduce*/
600 : GEN
601 490588 : Fq_mulu(GEN x, ulong y, /*unused*/GEN T, GEN p)
602 : {
603 : (void) T;
604 490588 : return typ(x)==t_POL ? FpX_Fp_mul(x,utoi(y),p): Fp_mulu(x, y, p);
605 : }
606 :
607 : /* y t_INT */
608 : GEN
609 44041 : Fq_Fp_mul(GEN x, GEN y, GEN T/*unused*/, GEN p)
610 : {
611 : (void)T;
612 6822 : return (typ(x) == t_POL)? FpX_Fp_mul(x,y,p)
613 50863 : : Fp_mul(x,y,p);
614 : }
615 : /* If T==NULL do not reduce*/
616 : GEN
617 611417 : Fq_sqr(GEN x, GEN T, GEN p)
618 : {
619 611417 : if (typ(x) == t_POL)
620 : {
621 70585 : if (T) return FpXQ_sqr(x,T,p);
622 0 : else return FpX_sqr(x,p);
623 : }
624 : else
625 540832 : return Fp_sqr(x,p);
626 : }
627 :
628 : GEN
629 0 : Fq_neg_inv(GEN x, GEN T, GEN p)
630 : {
631 0 : if (typ(x) == t_INT) return Fp_inv(Fp_neg(x,p),p);
632 0 : return FpXQ_inv(FpX_neg(x,p),T,p);
633 : }
634 :
635 : GEN
636 0 : Fq_invsafe(GEN x, GEN pol, GEN p)
637 : {
638 0 : if (typ(x) == t_INT) return Fp_invsafe(x,p);
639 0 : return FpXQ_invsafe(x,pol,p);
640 : }
641 :
642 : GEN
643 89312 : Fq_inv(GEN x, GEN pol, GEN p)
644 : {
645 89312 : if (typ(x) == t_INT) return Fp_inv(x,p);
646 81546 : return FpXQ_inv(x,pol,p);
647 : }
648 :
649 : GEN
650 343791 : Fq_div(GEN x, GEN y, GEN pol, GEN p)
651 : {
652 343791 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
653 : {
654 318402 : case 0: return Fp_div(x,y,p);
655 16702 : case 1: return FpX_Fp_div(x,y,p);
656 406 : case 2: return FpX_Fp_mul(FpXQ_inv(y,pol,p),x,p);
657 8281 : case 3: return FpXQ_div(x,y,pol,p);
658 : }
659 : return NULL;/*LCOV_EXCL_LINE*/
660 : }
661 :
662 : GEN
663 1088145 : Fq_pow(GEN x, GEN n, GEN pol, GEN p)
664 : {
665 1088145 : if (typ(x) == t_INT) return Fp_pow(x,n,p);
666 137448 : return FpXQ_pow(x,n,pol,p);
667 : }
668 :
669 : GEN
670 15050 : Fq_powu(GEN x, ulong n, GEN pol, GEN p)
671 : {
672 15050 : if (typ(x) == t_INT) return Fp_powu(x,n,p);
673 1267 : return FpXQ_powu(x,n,pol,p);
674 : }
675 :
676 : GEN
677 1892926 : Fq_sqrt(GEN x, GEN T, GEN p)
678 : {
679 1892926 : if (typ(x) == t_INT)
680 : {
681 1825064 : if (!T || odd(get_FpX_degree(T))) return Fp_sqrt(x,p);
682 9621 : x = scalarpol_shallow(x, get_FpX_var(T));
683 : }
684 77483 : return FpXQ_sqrt(x,T,p);
685 : }
686 : GEN
687 170786 : Fq_sqrtn(GEN x, GEN n, GEN T, GEN p, GEN *zeta)
688 : {
689 170786 : if (typ(x) == t_INT)
690 : {
691 : long d;
692 170415 : if (!T) return Fp_sqrtn(x,n,p,zeta);
693 126 : d = get_FpX_degree(T);
694 126 : if (ugcdiu(n,d) == 1)
695 : {
696 56 : if (!zeta) return Fp_sqrtn(x,n,p,NULL);
697 : /* gcd(n,p-1)=gcd(n,q-1): same number of solutions in Fp and F_q */
698 21 : if (equalii(gcdii(subiu(p,1),n), gcdii(subiu(Fp_powu(p,d,n), 1), n)))
699 14 : return Fp_sqrtn(x,n,p,zeta);
700 : }
701 77 : x = scalarpol(x, get_FpX_var(T)); /* left on stack */
702 : }
703 448 : return FpXQ_sqrtn(x,n,T,p,zeta);
704 : }
705 :
706 : struct _Fq_field
707 : {
708 : GEN T, p;
709 : };
710 :
711 : static GEN
712 303247 : _Fq_red(void *E, GEN x)
713 303247 : { struct _Fq_field *s = (struct _Fq_field *)E;
714 303247 : return Fq_red(x, s->T, s->p);
715 : }
716 :
717 : static GEN
718 362523 : _Fq_add(void *E, GEN x, GEN y)
719 : {
720 : (void) E;
721 362523 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
722 : {
723 14 : case 0: return addii(x,y);
724 0 : case 1: return ZX_Z_add(x,y);
725 15918 : case 2: return ZX_Z_add(y,x);
726 346591 : default: return ZX_add(x,y);
727 : }
728 : }
729 :
730 : static GEN
731 315028 : _Fq_neg(void *E, GEN x) { (void) E; return typ(x)==t_POL?ZX_neg(x):negi(x); }
732 :
733 : static GEN
734 243341 : _Fq_mul(void *E, GEN x, GEN y)
735 : {
736 : (void) E;
737 243341 : switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
738 : {
739 679 : case 0: return mulii(x,y);
740 1085 : case 1: return ZX_Z_mul(x,y);
741 1043 : case 2: return ZX_Z_mul(y,x);
742 240534 : default: return ZX_mul(x,y);
743 : }
744 : }
745 :
746 : static GEN
747 9331 : _Fq_inv(void *E, GEN x)
748 9331 : { struct _Fq_field *s = (struct _Fq_field *)E;
749 9331 : return Fq_inv(x,s->T,s->p);
750 : }
751 :
752 : static int
753 38388 : _Fq_equal0(GEN x) { return signe(x)==0; }
754 :
755 : static GEN
756 4151 : _Fq_s(void *E, long x) { (void) E; return stoi(x); }
757 :
758 : static const struct bb_field Fq_field={_Fq_red,_Fq_add,_Fq_mul,_Fq_neg,
759 : _Fq_inv,_Fq_equal0,_Fq_s};
760 :
761 4725 : const struct bb_field *get_Fq_field(void **E, GEN T, GEN p)
762 : {
763 4725 : if (!T)
764 0 : return get_Fp_field(E, p);
765 : else
766 : {
767 4725 : GEN z = new_chunk(sizeof(struct _Fq_field));
768 4725 : struct _Fq_field *e = (struct _Fq_field *) z;
769 4725 : e->T = T; e->p = p; *E = (void*)e;
770 4725 : return &Fq_field;
771 : }
772 : }
773 :
774 : /*******************************************************************/
775 : /* */
776 : /* Fq[X] */
777 : /* */
778 : /*******************************************************************/
779 : /* FpV_red(vecbinomial(n), p) */
780 : static GEN
781 0 : vecbinomial_Fp(long n, GEN p)
782 : {
783 0 : GEN C = vecbinomial(n) + 1;
784 0 : long k, d = (n + 1) >> 1;
785 0 : for (k = d; k >= 1; k--)
786 : {
787 0 : GEN a = gel(C,k);
788 0 : if (cmpii(a, p) < 0) break;
789 0 : gel(C,k) = gel(C, n-k) = remii(a, p);
790 : }
791 0 : return C - 1;
792 : }
793 : /* (X+u)^n */
794 : static GEN
795 434 : Fp_XpN_powu(GEN u, ulong n, GEN p, long v)
796 : {
797 : pari_sp av;
798 : ulong k;
799 : GEN B, C, V;
800 434 : if (!n) return pol_1(v);
801 434 : if (is_pm1(u))
802 434 : return Xpm1_powu(n, signe(u), v);
803 0 : av = avma;
804 0 : V = Fp_powers(u, n, p);
805 0 : B = vecbinomial_Fp(n, p);
806 0 : C = cgetg(n+3, t_POL);
807 0 : C[1] = evalsigne(1)| evalvarn(v);
808 0 : for (k=1; k <= n+1; k++)
809 0 : gel(C,k+1) = Fp_mul(gel(V,n+2-k), gel(B,k), p);
810 0 : return gc_upto(av, C);
811 : }
812 :
813 : static GEN
814 700 : FpX_Fp_translate_basecase(GEN P, GEN c, GEN p)
815 : {
816 700 : pari_sp av = avma;
817 : GEN Q;
818 : long i, k, n;
819 :
820 700 : if (!signe(P) || !signe(c)) return ZX_copy(P);
821 560 : Q = leafcopy(P); n = degpol(P);
822 1316 : for (i=1; i<=n; i++)
823 : {
824 2016 : for (k=n-i; k<n; k++)
825 1260 : gel(Q,2+k) = Fp_add(gel(Q,2+k), Fp_mul(c, gel(Q,2+k+1), p), p);
826 :
827 756 : if (gc_needed(av,2))
828 : {
829 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FpX_Fp_translate, i = %ld/%ld", i,n);
830 0 : Q = gc_GEN(av, Q);
831 : }
832 : }
833 560 : return gc_GEN(av, FpX_renormalize(Q, lg(Q)));
834 : }
835 :
836 : GEN
837 1134 : FpX_Fp_translate(GEN P, GEN c, GEN p)
838 : {
839 1134 : pari_sp av = avma;
840 1134 : long n = degpol(P);
841 1134 : if (n<=3 || 25*(n-3) < expi(p))
842 700 : return FpX_Fp_translate_basecase(P, c, p);
843 : else
844 : {
845 434 : long d = n >> 1;
846 434 : GEN Q = FpX_Fp_translate(RgX_shift_shallow(P, -d), c, p);
847 434 : GEN R = FpX_Fp_translate(RgXn_red_shallow(P, d), c, p);
848 434 : GEN S = Fp_XpN_powu(c, d, p, varn(P));
849 434 : return gc_upto(av, FpX_add(FpX_mul(Q, S, p), R, p));
850 : }
851 : }
852 :
853 : /* (X+u)^n mod (T,p) */
854 : static GEN
855 0 : FpXQX_XpN_powu(GEN u, ulong n, GEN T, GEN p, long v)
856 : {
857 : pari_sp av;
858 : ulong k;
859 : GEN B, C, V;
860 0 : if (!n) return pol_1(v);
861 0 : av = avma;
862 0 : V = FpXQ_powers(u, n, T, p);
863 0 : B = vecbinomial_Fp(n, p);
864 0 : C = cgetg(n+3, t_POL);
865 0 : C[1] = evalsigne(1)| evalvarn(v);
866 0 : for (k=1; k <= n+1; k++)
867 0 : gel(C,k+1) = Fq_mul(gel(V,n+2-k), gel(B,k), T, p);
868 0 : return gc_upto(av, C);
869 : }
870 :
871 : static GEN
872 33887 : FpXQX_FpXQ_translate_basecase(GEN P, GEN c, GEN T, GEN p)
873 : {
874 33887 : pari_sp av = avma;
875 : GEN Q;
876 : long i, k, n;
877 :
878 : /* signe works for t_(INT|POL) */
879 33887 : if (!signe(P) || !signe(c)) return RgX_copy(P);
880 33887 : Q = leafcopy(P); n = degpol(P);
881 150080 : for (i=1; i<=n; i++)
882 : {
883 433594 : for (k=n-i; k<n; k++)
884 317401 : gel(Q,2+k) = Fq_add(gel(Q,2+k), Fq_mul(c, gel(Q,2+k+1), T, p), T, p);
885 :
886 116193 : if (gc_needed(av,2))
887 : {
888 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FqX_Fq_translate, i = %ld/%ld", i,n);
889 0 : Q = gc_GEN(av, Q);
890 : }
891 : }
892 33887 : return gc_GEN(av, FpXQX_renormalize(Q, lg(Q)));
893 : }
894 :
895 : GEN
896 33887 : FpXQX_FpXQ_translate(GEN P, GEN c, GEN T, GEN p)
897 : {
898 33887 : pari_sp av = avma;
899 33887 : long n = degpol(P);
900 33887 : if (n < 220)
901 33887 : return FpXQX_FpXQ_translate_basecase(P, c, T, p);
902 : else
903 : {
904 0 : long d = n >> 1;
905 0 : GEN Q = FpXQX_FpXQ_translate(RgX_shift_shallow(P, -d), c, T, p);
906 0 : GEN R = FpXQX_FpXQ_translate(RgXn_red_shallow(P, d), c, T, p);
907 0 : GEN S = FpXQX_XpN_powu(c, d, T, p, varn(P));
908 0 : return gc_upto(av, FpXX_add(FpXQX_mul(Q, S, T, p), R, p));
909 : }
910 : }
911 :
912 : GEN
913 33880 : FqX_Fq_translate(GEN x, GEN y, GEN T, GEN p)
914 : {
915 33880 : if (!T) return FpX_Fp_translate(x,y,p);
916 33880 : if (typ(y)==t_INT)
917 : {
918 0 : pari_sp av = avma;
919 0 : y = scalarpol_shallow(y,varn(T));
920 0 : return gc_upto(av, FpXQX_FpXQ_translate(x,y,T,p));
921 : }
922 33880 : return FpXQX_FpXQ_translate(x,y,T,p);
923 : }
924 :
925 : GEN
926 40453 : FqV_roots_to_pol(GEN V, GEN T, GEN p, long v)
927 : {
928 40453 : pari_sp ltop = avma;
929 : long k;
930 : GEN W;
931 40453 : if (lgefint(p) == 3)
932 : {
933 31742 : ulong pp = p[2];
934 31742 : GEN Tl = ZX_to_Flx(T, pp);
935 31744 : GEN Vl = FqC_to_FlxqC(V, Tl, pp);
936 31746 : Tl = FlxqV_roots_to_pol(Vl, Tl, pp, v);
937 31746 : return gc_upto(ltop, FlxX_to_ZXX(Tl));
938 : }
939 8711 : W = cgetg(lg(V),t_VEC);
940 78117 : for(k=1; k < lg(V); k++)
941 69406 : gel(W,k) = deg1pol_shallow(gen_1,Fq_neg(gel(V,k),T,p),v);
942 8711 : return gc_upto(ltop, FpXQXV_prod(W, T, p));
943 : }
944 :
945 : GEN
946 109458 : FqV_red(GEN x, GEN T, GEN p)
947 778056 : { pari_APPLY_type(t_VEC, Fq_red(gel(x,i), T, p)) }
948 :
949 : GEN
950 24058 : FqC_red(GEN x, GEN T, GEN p)
951 163384 : { pari_APPLY_type(t_COL, Fq_red(gel(x,i), T, p)) }
952 :
953 : GEN
954 1701 : FqM_red(GEN x, GEN T, GEN p)
955 5411 : { pari_APPLY_same(FqC_red(gel(x,i), T, p)) }
956 :
957 : GEN
958 0 : FqC_add(GEN x, GEN y, GEN T, GEN p)
959 : {
960 0 : if (!T) return FpC_add(x, y, p);
961 0 : pari_APPLY_type(t_COL, Fq_add(gel(x,i), gel(y,i), T, p))
962 : }
963 :
964 : GEN
965 0 : FqC_sub(GEN x, GEN y, GEN T, GEN p)
966 : {
967 0 : if (!T) return FpC_sub(x, y, p);
968 0 : pari_APPLY_type(t_COL, Fq_sub(gel(x,i), gel(y,i), T, p))
969 : }
970 :
971 : GEN
972 0 : FqC_Fq_mul(GEN x, GEN y, GEN T, GEN p)
973 : {
974 0 : if (!T) return FpC_Fp_mul(x, y, p);
975 0 : pari_APPLY_type(t_COL, Fq_mul(gel(x,i),y,T,p))
976 : }
977 :
978 : GEN
979 105 : FqC_FqV_mul(GEN x, GEN y, GEN T, GEN p)
980 : {
981 105 : long i,j, lx=lg(x), ly=lg(y);
982 : GEN z;
983 105 : if (ly==1) return cgetg(1,t_MAT);
984 105 : z = cgetg(ly,t_MAT);
985 819 : for (j=1; j < ly; j++)
986 : {
987 714 : GEN zj = cgetg(lx,t_COL);
988 4200 : for (i=1; i<lx; i++) gel(zj,i) = Fq_mul(gel(x,i),gel(y,j), T, p);
989 714 : gel(z, j) = zj;
990 : }
991 105 : return z;
992 : }
993 :
994 : GEN
995 5467 : FpXC_center(GEN x, GEN p, GEN pov2)
996 41476 : { pari_APPLY_type(t_COL, FpX_center(gel(x,i), p, pov2)) }
997 :
998 : GEN
999 109021 : FqC_to_FlxqC(GEN x, GEN T, ulong p)
1000 109021 : { long sv = get_Flx_var(T);
1001 4834758 : pari_APPLY_type(t_COL,typ(gel(x,i))==t_INT ?
1002 : Z_to_Flx(gel(x,i), p, sv): ZX_to_Flx(gel(x,i), p)) }
1003 :
1004 : GEN
1005 8708 : FqM_to_FlxqM(GEN x, GEN T, ulong p)
1006 85985 : { pari_APPLY_same(FqC_to_FlxqC(gel(x,i), T, p)) }
1007 :
1008 : GEN
1009 1800 : FpXM_center(GEN x, GEN p, GEN pov2)
1010 7267 : { pari_APPLY_same(FpXC_center(gel(x,i), p, pov2)) }
1011 :
1012 : /*******************************************************************/
1013 : /* */
1014 : /* GENERIC CRT */
1015 : /* */
1016 : /*******************************************************************/
1017 : static GEN
1018 8398608 : primelist(forprime_t *S, long n, GEN dB)
1019 : {
1020 8398608 : GEN P = cgetg(n+1, t_VECSMALL);
1021 8398593 : long i = 1;
1022 : ulong p;
1023 20251399 : while (i <= n && (p = u_forprime_next(S)))
1024 11852806 : if (!dB || umodiu(dB, p)) P[i++] = p;
1025 8398597 : return P;
1026 : }
1027 :
1028 : void
1029 7815775 : gen_inccrt_i(const char *str, GEN worker, GEN dB, long n, long mmin,
1030 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
1031 : GEN center(GEN, GEN, GEN))
1032 : {
1033 7815775 : long m = mmin? minss(mmin, n): usqrt(n);
1034 : GEN H, P, mod;
1035 : pari_timer ti;
1036 7815772 : if (DEBUGLEVEL > 4)
1037 : {
1038 0 : timer_start(&ti);
1039 0 : err_printf("%s: nb primes: %ld\n",str, n);
1040 : }
1041 7815770 : if (m == 1)
1042 : {
1043 7504545 : GEN P = primelist(S, n, dB);
1044 7504516 : GEN done = closure_callgen1(worker, P);
1045 7504502 : H = gel(done,1);
1046 7504502 : mod = gel(done,2);
1047 7504502 : if (!*pH && center) H = center(H, mod, shifti(mod,-1));
1048 7504441 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
1049 : }
1050 : else
1051 : {
1052 311225 : long i, s = (n+m-1)/m, r = m - (m*s-n), di = 0;
1053 : struct pari_mt pt;
1054 311225 : long pending = 0;
1055 311225 : H = cgetg(m+1, t_VEC); P = cgetg(m+1, t_VEC);
1056 311225 : mt_queue_start_lim(&pt, worker, m);
1057 1270909 : for (i=1; i<=m || pending; i++)
1058 : {
1059 : GEN done;
1060 959684 : GEN pr = i <= m ? mkvec(primelist(S, i<=r ? s: s-1, dB)): NULL;
1061 959685 : mt_queue_submit(&pt, i, pr);
1062 959684 : done = mt_queue_get(&pt, NULL, &pending);
1063 959684 : if (done)
1064 : {
1065 894066 : di++;
1066 894066 : gel(H, di) = gel(done,1);
1067 894066 : gel(P, di) = gel(done,2);
1068 894066 : if (DEBUGLEVEL>5) err_printf("%ld%% ",100*di/m);
1069 : }
1070 : }
1071 311225 : mt_queue_end(&pt);
1072 311225 : if (DEBUGLEVEL>5) err_printf("\n");
1073 311225 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: modular", str);
1074 311225 : H = crt(H, P, &mod);
1075 311225 : if (DEBUGLEVEL>4) timer_printf(&ti,"%s: chinese", str);
1076 : }
1077 7815666 : if (*pH) H = crt(mkvec2(*pH, H), mkvec2(*pmod, mod), &mod);
1078 7815666 : *pH = H; *pmod = mod;
1079 7815666 : }
1080 : void
1081 2062487 : gen_inccrt(const char *str, GEN worker, GEN dB, long n, long mmin,
1082 : forprime_t *S, GEN *pH, GEN *pmod, GEN crt(GEN, GEN, GEN*),
1083 : GEN center(GEN, GEN, GEN))
1084 : {
1085 2062487 : pari_sp av = avma;
1086 2062487 : gen_inccrt_i(str, worker, dB, n, mmin, S, pH, pmod, crt, center);
1087 2062423 : (void)gc_all(av, 2, pH, pmod);
1088 2062569 : }
1089 :
1090 : GEN
1091 1274670 : gen_crt(const char *str, GEN worker, forprime_t *S, GEN dB, ulong bound, long mmin, GEN *pmod,
1092 : GEN crt(GEN, GEN, GEN*), GEN center(GEN, GEN, GEN))
1093 : {
1094 1274670 : GEN mod = gen_1, H = NULL;
1095 : ulong e;
1096 :
1097 1274670 : bound++;
1098 2549410 : while (bound > (e = expi(mod)))
1099 : {
1100 1274625 : long n = (bound - e) / expu(S->p) + 1;
1101 1274654 : gen_inccrt(str, worker, dB, n, mmin, S, &H, &mod, crt, center);
1102 : }
1103 1274719 : if (pmod) *pmod = mod;
1104 1274719 : return H;
1105 : }
1106 :
1107 : /*******************************************************************/
1108 : /* */
1109 : /* MODULAR GCD */
1110 : /* */
1111 : /*******************************************************************/
1112 : /* return z = a mod q, b mod p (p,q) = 1; qinv = 1/q mod p; a in ]-q,q] */
1113 : static GEN
1114 5162228 : Fl_chinese_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq, GEN pq2)
1115 : {
1116 5162228 : ulong d, amod = umodiu(a, p);
1117 5162306 : pari_sp av = avma;
1118 : GEN ax;
1119 :
1120 5162306 : if (b == amod) return NULL;
1121 2126307 : d = Fl_mul(Fl_sub(b, amod, p), qinv, p); /* != 0 */
1122 2127031 : if (d >= 1 + (p>>1))
1123 1037809 : ax = subii(a, mului(p-d, q));
1124 : else
1125 : {
1126 1089222 : ax = addii(a, mului(d, q)); /* in ]0, pq[ assuming a in ]-q,q[ */
1127 1088634 : if (cmpii(ax,pq2) > 0) ax = subii(ax,pq);
1128 : }
1129 2125880 : return gc_INT(av, ax);
1130 : }
1131 : GEN
1132 406 : Z_init_CRT(ulong Hp, ulong p) { return stoi(Fl_center(Hp, p, p>>1)); }
1133 : GEN
1134 32487 : ZX_init_CRT(GEN Hp, ulong p, long v)
1135 : {
1136 32487 : long i, l = lg(Hp), lim = (long)(p>>1);
1137 32487 : GEN H = cgetg(l, t_POL);
1138 32487 : H[1] = evalsigne(1) | evalvarn(v);
1139 801143 : for (i=2; i<l; i++)
1140 768656 : gel(H,i) = stoi(Fl_center(Hp[i], p, lim));
1141 32487 : return ZX_renormalize(H,l);
1142 : }
1143 :
1144 : GEN
1145 5978 : ZM_init_CRT(GEN Hp, ulong p)
1146 : {
1147 5978 : long i,j, m, l = lg(Hp), lim = (long)(p>>1);
1148 5978 : GEN c, cp, H = cgetg(l, t_MAT);
1149 5978 : if (l==1) return H;
1150 5978 : m = lgcols(Hp);
1151 20223 : for (j=1; j<l; j++)
1152 : {
1153 14245 : cp = gel(Hp,j);
1154 14245 : c = cgetg(m, t_COL);
1155 14245 : gel(H,j) = c;
1156 169316 : for (i=1; i<m; i++) gel(c,i) = stoi(Fl_center(cp[i],p, lim));
1157 : }
1158 5978 : return H;
1159 : }
1160 :
1161 : int
1162 7742 : Z_incremental_CRT(GEN *H, ulong Hp, GEN *ptq, ulong p)
1163 : {
1164 7742 : GEN h, q = *ptq, qp = muliu(q,p);
1165 7742 : ulong qinv = Fl_inv(umodiu(q,p), p);
1166 7742 : int stable = 1;
1167 7742 : h = Fl_chinese_coprime(*H,Hp,q,p,qinv,qp,shifti(qp,-1));
1168 7742 : if (h) { *H = h; stable = 0; }
1169 7742 : *ptq = qp; return stable;
1170 : }
1171 :
1172 : static int
1173 148378 : ZX_incremental_CRT_raw(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p)
1174 : {
1175 148378 : GEN H = *ptH, h, qp2 = shifti(qp,-1);
1176 148373 : ulong qinv = Fl_inv(umodiu(q,p), p);
1177 148377 : long i, l = lg(H), lp = lg(Hp);
1178 148377 : int stable = 1;
1179 :
1180 148377 : if (l < lp)
1181 : { /* degree increases */
1182 0 : GEN x = cgetg(lp, t_POL);
1183 0 : for (i=1; i<l; i++) x[i] = H[i];
1184 0 : for ( ; i<lp; i++) gel(x,i) = gen_0;
1185 0 : *ptH = H = x;
1186 0 : stable = 0;
1187 148377 : } else if (l > lp)
1188 : { /* degree decreases */
1189 0 : GEN x = cgetg(l, t_VECSMALL);
1190 0 : for (i=1; i<lp; i++) x[i] = Hp[i];
1191 0 : for ( ; i<l; i++) x[i] = 0;
1192 0 : Hp = x; lp = l;
1193 : }
1194 4939092 : for (i=2; i<lp; i++)
1195 : {
1196 4790838 : h = Fl_chinese_coprime(gel(H,i),Hp[i],q,p,qinv,qp,qp2);
1197 4790715 : if (h) { gel(H,i) = h; stable = 0; }
1198 : }
1199 148254 : (void)ZX_renormalize(H,lp);
1200 148378 : return stable;
1201 : }
1202 :
1203 : int
1204 0 : ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN *ptq, ulong p)
1205 : {
1206 0 : GEN q = *ptq, qp = muliu(q,p);
1207 0 : int stable = ZX_incremental_CRT_raw(ptH, Hp, q, qp, p);
1208 0 : *ptq = qp; return stable;
1209 : }
1210 :
1211 : int
1212 7787 : ZM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1213 : {
1214 7787 : GEN h, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1215 7787 : ulong qinv = Fl_inv(umodiu(q,p), p);
1216 7787 : long i,j, l = lg(H), m = lgcols(H);
1217 7787 : int stable = 1;
1218 27451 : for (j=1; j<l; j++)
1219 205516 : for (i=1; i<m; i++)
1220 : {
1221 185852 : h = Fl_chinese_coprime(gcoeff(H,i,j), coeff(Hp,i,j),q,p,qinv,qp,qp2);
1222 185852 : if (h) { gcoeff(H,i,j) = h; stable = 0; }
1223 : }
1224 7787 : *ptq = qp; return stable;
1225 : }
1226 :
1227 : GEN
1228 679 : ZXM_init_CRT(GEN Hp, long deg, ulong p)
1229 : {
1230 : long i, j, k;
1231 : GEN H;
1232 679 : long m, l = lg(Hp), lim = (long)(p>>1), n;
1233 679 : H = cgetg(l, t_MAT);
1234 679 : if (l==1) return H;
1235 679 : m = lgcols(Hp);
1236 679 : n = deg + 3;
1237 2268 : for (j=1; j<l; j++)
1238 : {
1239 1589 : GEN cp = gel(Hp,j);
1240 1589 : GEN c = cgetg(m, t_COL);
1241 1589 : gel(H,j) = c;
1242 24465 : for (i=1; i<m; i++)
1243 : {
1244 22876 : GEN dp = gel(cp, i);
1245 22876 : long l = lg(dp);
1246 22876 : GEN d = cgetg(n, t_POL);
1247 22876 : gel(c, i) = d;
1248 22876 : d[1] = dp[1] | evalsigne(1);
1249 46459 : for (k=2; k<l; k++)
1250 23583 : gel(d,k) = stoi(Fl_center(dp[k], p, lim));
1251 45493 : for ( ; k<n; k++)
1252 22617 : gel(d,k) = gen_0;
1253 : }
1254 : }
1255 679 : return H;
1256 : }
1257 :
1258 : int
1259 653 : ZXM_incremental_CRT(GEN *pH, GEN Hp, GEN *ptq, ulong p)
1260 : {
1261 653 : GEN v, H = *pH, q = *ptq, qp = muliu(q, p), qp2 = shifti(qp,-1);
1262 653 : ulong qinv = Fl_inv(umodiu(q,p), p);
1263 653 : long i,j,k, l = lg(H), m = lgcols(H), n = lg(gmael(H,1,1));
1264 653 : int stable = 1;
1265 2225 : for (j=1; j<l; j++)
1266 90418 : for (i=1; i<m; i++)
1267 : {
1268 88846 : GEN h = gmael(H,j,i), hp = gmael(Hp,j,i);
1269 88846 : long lh = lg(hp);
1270 246641 : for (k=2; k<lh; k++)
1271 : {
1272 157795 : v = Fl_chinese_coprime(gel(h,k),uel(hp,k),q,p,qinv,qp,qp2);
1273 157795 : if (v) { gel(h,k) = v; stable = 0; }
1274 : }
1275 108763 : for (; k<n; k++)
1276 : {
1277 19917 : v = Fl_chinese_coprime(gel(h,k),0,q,p,qinv,qp,qp2);
1278 19917 : if (v) { gel(h,k) = v; stable = 0; }
1279 : }
1280 : }
1281 653 : *ptq = qp; return stable;
1282 : }
1283 :
1284 : /* record the degrees of Euclidean remainders (make them as large as
1285 : * possible : smaller values correspond to a degenerate sequence) */
1286 : static void
1287 23741 : Flx_resultant_set_dglist(GEN a, GEN b, GEN dglist, ulong p)
1288 : {
1289 : long da,db,dc, ind;
1290 23741 : pari_sp av = avma;
1291 :
1292 23741 : if (lgpol(a)==0 || lgpol(b)==0) return;
1293 22474 : da = degpol(a);
1294 22474 : db = degpol(b);
1295 22474 : if (db > da)
1296 0 : { swapspec(a,b, da,db); }
1297 22474 : else if (!da) return;
1298 22474 : ind = 0;
1299 145766 : while (db)
1300 : {
1301 123295 : GEN c = Flx_rem(a,b, p);
1302 123292 : a = b; b = c; dc = degpol(c);
1303 123292 : if (dc < 0) break;
1304 :
1305 123292 : ind++;
1306 123292 : if (dc > dglist[ind]) dglist[ind] = dc;
1307 123292 : if (gc_needed(av,2))
1308 : {
1309 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1310 0 : (void)gc_all(av, 2, &a,&b);
1311 : }
1312 123292 : db = dc; /* = degpol(b) */
1313 : }
1314 22471 : if (ind+1 > lg(dglist)) setlg(dglist,ind+1);
1315 22471 : set_avma(av);
1316 : }
1317 : /* assuming the PRS finishes on a degree 1 polynomial C0 + C1X, with
1318 : * "generic" degree sequence as given by dglist, set *Ci and return
1319 : * resultant(a,b). Modular version of Collins's subresultant */
1320 : static ulong
1321 2089707 : Flx_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p)
1322 : {
1323 : long da,db,dc, ind;
1324 2089707 : ulong lb, res, g = 1UL, h = 1UL, ca = 1UL, cb = 1UL;
1325 2089707 : int s = 1;
1326 2089707 : pari_sp av = avma;
1327 :
1328 2089707 : *C0 = 1; *C1 = 0;
1329 2089707 : if (lgpol(a)==0 || lgpol(b)==0) return 0;
1330 2080262 : da = degpol(a);
1331 2080281 : db = degpol(b);
1332 2080260 : if (db > da)
1333 : {
1334 0 : swapspec(a,b, da,db);
1335 0 : if (both_odd(da,db)) s = -s;
1336 : }
1337 2080260 : else if (!da) return 1; /* = a[2] ^ db, since 0 <= db <= da = 0 */
1338 2080260 : ind = 0;
1339 19812520 : while (db)
1340 : { /* sub-resultant algo., applied to ca * a and cb * b, ca,cb scalars,
1341 : * da = deg a, db = deg b */
1342 17737012 : GEN c = Flx_rem(a,b, p);
1343 17608050 : long delta = da - db;
1344 :
1345 17608050 : if (both_odd(da,db)) s = -s;
1346 17604971 : lb = Fl_mul(b[db+2], cb, p);
1347 17615179 : a = b; b = c; dc = degpol(c);
1348 17618752 : ind++;
1349 17618752 : if (dc != dglist[ind]) return gc_ulong(av,0); /* degenerates */
1350 17613776 : if (g == h)
1351 : { /* frequent */
1352 17553813 : ulong cc = Fl_mul(ca, Fl_powu(Fl_div(lb,g,p), delta+1, p), p);
1353 17672621 : ca = cb;
1354 17672621 : cb = cc;
1355 : }
1356 : else
1357 : {
1358 59963 : ulong cc = Fl_mul(ca, Fl_powu(lb, delta+1, p), p);
1359 59966 : ulong ghdelta = Fl_mul(g, Fl_powu(h, delta, p), p);
1360 59967 : ca = cb;
1361 59967 : cb = Fl_div(cc, ghdelta, p);
1362 : }
1363 17732971 : da = db; /* = degpol(a) */
1364 17732971 : db = dc; /* = degpol(b) */
1365 :
1366 17732971 : g = lb;
1367 17732971 : if (delta == 1)
1368 17632297 : h = g; /* frequent */
1369 : else
1370 100674 : h = Fl_mul(h, Fl_powu(Fl_div(g,h,p), delta, p), p);
1371 :
1372 17731156 : if (gc_needed(av,2))
1373 : {
1374 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_resultant_all");
1375 0 : (void)gc_all(av, 2, &a,&b);
1376 : }
1377 : }
1378 2075508 : if (da > 1) return 0; /* Failure */
1379 : /* last nonconstant polynomial has degree 1 */
1380 2075508 : *C0 = Fl_mul(ca, a[2], p);
1381 2075464 : *C1 = Fl_mul(ca, a[3], p);
1382 2075483 : res = Fl_mul(cb, b[2], p);
1383 2075474 : if (s == -1) res = p - res;
1384 2075474 : return gc_ulong(av,res);
1385 : }
1386 :
1387 : /* Q a vector of polynomials representing B in Fp[X][Y], evaluate at X = x,
1388 : * Return 0 in case of degree drop. */
1389 : static GEN
1390 2113585 : FlxY_evalx_drop(GEN Q, ulong x, ulong p)
1391 : {
1392 : GEN z;
1393 2113585 : long i, lb = lg(Q);
1394 2113585 : ulong leadz = Flx_eval(leading_coeff(Q), x, p);
1395 2113287 : long vs=mael(Q,2,1);
1396 2113287 : if (!leadz) return zero_Flx(vs);
1397 :
1398 2102627 : z = cgetg(lb, t_VECSMALL); z[1] = vs;
1399 20090686 : for (i=2; i<lb-1; i++) z[i] = Flx_eval(gel(Q,i), x, p);
1400 2101339 : z[i] = leadz; return z;
1401 : }
1402 :
1403 : GEN
1404 2072 : FpXY_FpXQ_evaly(GEN Q, GEN y, GEN T, GEN p, long vx)
1405 : {
1406 2072 : pari_sp av = avma;
1407 2072 : long i, lb = lg(Q);
1408 : GEN z;
1409 2072 : if (lb == 2) return pol_0(vx);
1410 2072 : z = gel(Q, lb-1);
1411 2072 : if (lb == 3 || !signe(y)) return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1412 :
1413 2072 : if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1414 48636 : for (i=lb-2; i>=2; i--)
1415 : {
1416 46564 : GEN c = gel(Q,i);
1417 46564 : z = FqX_Fq_mul(z, y, T, p);
1418 46564 : z = typ(c) == t_INT? FqX_Fq_add(z,c,T,p): FqX_add(z,c,T,p);
1419 : }
1420 2072 : return gc_upto(av, z);
1421 : }
1422 :
1423 : static GEN
1424 292120 : ZX_norml1(GEN x)
1425 : {
1426 292120 : long i, l = lg(x);
1427 : GEN s;
1428 :
1429 292120 : if (l == 2) return gen_0;
1430 199566 : s = gel(x, l-1); /* != 0 */
1431 698309 : for (i = l-2; i > 1; i--) {
1432 498748 : GEN xi = gel(x,i);
1433 498748 : if (!signe(xi)) continue;
1434 259796 : s = addii_sign(s,1, xi,1);
1435 : }
1436 199561 : return s;
1437 : }
1438 : /* x >= 0, y != 0, return x + |y| */
1439 : static GEN
1440 25552 : addii_abs(GEN x, GEN y)
1441 : {
1442 25552 : if (!signe(x)) return absi_shallow(y);
1443 16041 : return addii_sign(x,1, y,1);
1444 : }
1445 :
1446 : /* x a ZX, return sum_{i >= k} |x[i]| binomial(i, k) */
1447 : static GEN
1448 31644 : ZX_norml1_1(GEN x, long k)
1449 : {
1450 31644 : long i, d = degpol(x);
1451 : GEN s, C; /* = binomial(i, k) */
1452 :
1453 31643 : if (!d || k > d) return gen_0;
1454 31643 : s = absi_shallow(gel(x, k+2)); /* may be 0 */
1455 31644 : C = gen_1;
1456 68050 : for (i = k+1; i <= d; i++) {
1457 36404 : GEN xi = gel(x,i+2);
1458 36404 : if (k) C = diviuexact(muliu(C, i), i-k);
1459 36407 : if (signe(xi)) s = addii_abs(s, mulii(C, xi));
1460 : }
1461 31646 : return s;
1462 : }
1463 : /* x has non-negative real coefficients */
1464 : static GEN
1465 3283 : RgX_norml1_1(GEN x, long k)
1466 : {
1467 3283 : long i, d = degpol(x);
1468 : GEN s, C; /* = binomial(i, k) */
1469 :
1470 3283 : if (!d || k > d) return gen_0;
1471 3283 : s = gel(x, k+2); /* may be 0 */
1472 3283 : C = gen_1;
1473 9198 : for (i = k+1; i <= d; i++) {
1474 5915 : GEN xi = gel(x,i+2);
1475 5915 : if (k) C = diviuexact(muliu(C, i), i-k);
1476 5915 : if (!gequal0(xi)) s = gadd(s, gmul(C, xi));
1477 : }
1478 3283 : return s;
1479 : }
1480 :
1481 : /* N_2(A)^2 */
1482 : static GEN
1483 9019 : sqrN2(GEN A, long prec)
1484 : {
1485 9019 : pari_sp av = avma;
1486 9019 : long i, l = lg(A);
1487 9019 : GEN a = gen_0;
1488 43727 : for (i = 2; i < l; i++)
1489 : {
1490 34708 : a = gadd(a, gabs(gnorm(gel(A,i)), prec));
1491 34708 : if (gc_needed(av,1))
1492 : {
1493 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1494 0 : a = gc_upto(av, a);
1495 : }
1496 : }
1497 9019 : return a;
1498 : }
1499 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1500 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1501 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1502 : * Return e such that Res(A, B) < 2^e */
1503 : static GEN
1504 8165 : RgX_RgXY_ResBound(GEN A, GEN B, long prec)
1505 : {
1506 8165 : pari_sp av = avma;
1507 8165 : GEN b = gen_0, bnd;
1508 8165 : long i, lB = lg(B);
1509 31747 : for (i=2; i<lB; i++)
1510 : {
1511 23582 : GEN t = gel(B,i);
1512 23582 : if (typ(t) == t_POL) t = gnorml1(t, prec);
1513 23582 : b = gadd(b, gabs(gsqr(t), prec));
1514 23582 : if (gc_needed(av,1))
1515 : {
1516 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1517 0 : b = gc_upto(av, b);
1518 : }
1519 : }
1520 8165 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1521 : gpowgs(b, degpol(A))), prec);
1522 8165 : return gc_upto(av, bnd);
1523 : }
1524 : /* A,B in C[X] return RgX_RgXY_ResBound(A, B(x+y)) */
1525 : static GEN
1526 854 : RgX_RgXY_ResBound_1(GEN A, GEN B, long prec)
1527 : {
1528 854 : pari_sp av = avma, av2;
1529 854 : GEN b = gen_0, bnd;
1530 854 : long i, lB = lg(B);
1531 854 : B = shallowcopy(B);
1532 4137 : for (i=2; i<lB; i++) gel(B,i) = gabs(gel(B,i), prec);
1533 854 : av2 = avma;
1534 4137 : for (i=2; i<lB; i++)
1535 : {
1536 3283 : b = gadd(b, gsqr(RgX_norml1_1(B, i-2)));
1537 3283 : if (gc_needed(av2,1))
1538 : {
1539 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_RgXY_ResBound i = %ld",i);
1540 0 : b = gc_upto(av2, b);
1541 : }
1542 : }
1543 854 : bnd = gsqrt(gmul(gpowgs(sqrN2(A,prec), degpol(B)),
1544 : gpowgs(b, degpol(A))), prec);
1545 854 : return gc_upto(av, bnd);
1546 : }
1547 :
1548 : /* log2 N_2(A)^2 */
1549 : static double
1550 177023 : log2N2(GEN A)
1551 : {
1552 177023 : pari_sp av = avma;
1553 177023 : long i, l = lg(A);
1554 177023 : GEN a = gen_0;
1555 1337216 : for (i=2; i < l; i++)
1556 : {
1557 1160196 : a = addii(a, sqri(gel(A,i)));
1558 1160191 : if (gc_needed(av,1))
1559 : {
1560 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1561 0 : a = gc_upto(av, a);
1562 : }
1563 : }
1564 177020 : return gc_double(av, dbllog2(a));
1565 : }
1566 : /* Interpolate at roots of 1 and use Hadamard bound for univariate resultant:
1567 : * bound = N_2(A)^degpol B N_2(B)^degpol(A), where
1568 : * N_2(A) = sqrt(sum (N_1(Ai))^2)
1569 : * Return e such that Res(A, B) < 2^e */
1570 : ulong
1571 166943 : ZX_ZXY_ResBound(GEN A, GEN B, GEN dB)
1572 : {
1573 166943 : pari_sp av = avma;
1574 166943 : GEN b = gen_0;
1575 166943 : long i, lB = lg(B);
1576 : double logb;
1577 1262247 : for (i=2; i<lB; i++)
1578 : {
1579 1095322 : GEN t = gel(B,i);
1580 1095322 : if (typ(t) == t_POL) t = ZX_norml1(t);
1581 1095324 : b = addii(b, sqri(t));
1582 1095304 : if (gc_needed(av,1))
1583 : {
1584 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1585 0 : b = gc_upto(av, b);
1586 : }
1587 : }
1588 166925 : logb = dbllog2(b); if (dB) logb -= 2 * dbllog2(dB);
1589 166936 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * logb) / 2);
1590 166944 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1591 : }
1592 : /* A,B ZX. Return ZX_ZXY_ResBound(A(x), B(x+y)) */
1593 : static ulong
1594 10086 : ZX_ZXY_ResBound_1(GEN A, GEN B)
1595 : {
1596 10086 : pari_sp av = avma;
1597 10086 : GEN b = gen_0;
1598 10086 : long i, lB = lg(B);
1599 41728 : for (i=2; i<lB; i++)
1600 : {
1601 31644 : b = addii(b, sqri(ZX_norml1_1(B, i-2)));
1602 31642 : if (gc_needed(av,1))
1603 : {
1604 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_ResBound i = %ld",i);
1605 0 : b = gc_upto(av, b);
1606 : }
1607 : }
1608 10084 : i = (long)((degpol(B) * log2N2(A) + degpol(A) * dbllog2(b)) / 2);
1609 10085 : return gc_ulong(av, (i <= 0)? 1: 1 + (ulong)i);
1610 : }
1611 : /* special case B = A' */
1612 : static ulong
1613 1134923 : ZX_discbound(GEN A)
1614 : {
1615 1134923 : pari_sp av = avma;
1616 1134923 : GEN a = gen_0, b = gen_0;
1617 1134923 : long i , lA = lg(A), dA = degpol(A);
1618 : double loga, logb;
1619 6771960 : for (i = 2; i < lA; i++)
1620 : {
1621 5637247 : GEN c = sqri(gel(A,i));
1622 5637011 : a = addii(a, c);
1623 5637142 : if (i > 2) b = addii(b, mulii(c, sqru(i-2)));
1624 5637085 : if (gc_needed(av,1))
1625 : {
1626 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZX_discbound i = %ld",i);
1627 0 : (void)gc_all(av, 2, &a, &b);
1628 : }
1629 : }
1630 1134713 : loga = dbllog2(a);
1631 1134827 : logb = dbllog2(b); set_avma(av);
1632 1134855 : i = (long)(((dA-1) * loga + dA * logb) / 2);
1633 1134855 : return (i <= 0)? 1: 1 + (ulong)i;
1634 : }
1635 :
1636 : /* return Res(a(Y), b(n,Y)) over Fp. la = leading_coeff(a) [for efficiency] */
1637 : static ulong
1638 5539434 : Flx_FlxY_eval_resultant(GEN a, GEN b, ulong n, ulong p, ulong pi, ulong la)
1639 : {
1640 5539434 : GEN ev = FlxY_evalx_pre(b, n, p, pi);
1641 5540221 : long drop = lg(b) - lg(ev);
1642 5540221 : ulong r = Flx_resultant_pre(a, ev, p, pi);
1643 5539321 : if (drop && la != 1) r = Fl_mul(r, Fl_powu_pre(la, drop, p, pi), p);
1644 5539332 : return r;
1645 : }
1646 : static GEN
1647 284 : FpX_FpXY_eval_resultant(GEN a, GEN b, GEN n, GEN p, GEN la, long db, long vX)
1648 : {
1649 284 : GEN ev = FpXY_evaly(b, n, p, vX);
1650 284 : long drop = db-degpol(ev);
1651 284 : GEN r = FpX_resultant(a, ev, p);
1652 284 : if (drop && !gequal1(la)) r = Fp_mul(r, Fp_powu(la, drop,p),p);
1653 284 : return r;
1654 : }
1655 :
1656 : /* assume dres := deg(Res_X(a,b), Y) <= deg(a,X) * deg(b,Y) < p */
1657 : /* Return a Fly */
1658 : static GEN
1659 368517 : Flx_FlxY_resultant_polint(GEN a, GEN b, ulong p, ulong pi, long dres, long sx)
1660 : {
1661 : long i;
1662 368517 : ulong n, la = Flx_lead(a);
1663 368517 : GEN x = cgetg(dres+2, t_VECSMALL);
1664 368517 : GEN y = cgetg(dres+2, t_VECSMALL);
1665 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1666 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1667 2957527 : for (i=0,n = 1; i < dres; n++)
1668 : {
1669 2589015 : x[++i] = n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1670 2588907 : x[++i] = p-n; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1671 : }
1672 368512 : if (i == dres)
1673 : {
1674 363006 : x[++i] = 0; y[i] = Flx_FlxY_eval_resultant(a,b, x[i], p,pi,la);
1675 : }
1676 368513 : return Flv_polint(x,y, p, sx);
1677 : }
1678 :
1679 : static GEN
1680 7522 : FlxX_pseudorem(GEN x, GEN y, ulong p, ulong pi)
1681 : {
1682 7522 : long vx = varn(x), dx, dy, dz, i, lx, dp;
1683 7522 : pari_sp av = avma, av2;
1684 :
1685 7522 : if (!signe(y)) pari_err_INV("FlxX_pseudorem",y);
1686 7522 : (void)new_chunk(2);
1687 7525 : dx=degpol(x); x = RgX_recip_i(x)+2;
1688 7529 : dy=degpol(y); y = RgX_recip_i(y)+2; dz=dx-dy; dp = dz+1;
1689 7529 : av2 = avma;
1690 : for (;;)
1691 : {
1692 61652 : gel(x,0) = Flx_neg(gel(x,0), p); dp--;
1693 230846 : for (i=1; i<=dy; i++)
1694 168881 : gel(x,i) = Flx_add( Flx_mul_pre(gel(y,0), gel(x,i), p, pi),
1695 169158 : Flx_mul_pre(gel(x,0), gel(y,i), p, pi), p );
1696 1111810 : for ( ; i<=dx; i++)
1697 1051030 : gel(x,i) = Flx_mul_pre(gel(y,0), gel(x,i), p, pi);
1698 65546 : do { x++; dx--; } while (dx >= 0 && lg(gel(x,0))==2);
1699 60780 : if (dx < dy) break;
1700 53249 : if (gc_needed(av2,1))
1701 : {
1702 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_pseudorem dx = %ld >= %ld",dx,dy);
1703 0 : gc_slice(av2,x,dx+1);
1704 : }
1705 : }
1706 7531 : if (dx < 0) return zero_Flx(0);
1707 7531 : lx = dx+3; x -= 2;
1708 7531 : x[0]=evaltyp(t_POL) | _evallg(lx);
1709 7531 : x[1]=evalsigne(1) | evalvarn(vx);
1710 7531 : x = RgX_recip_i(x);
1711 7530 : if (dp)
1712 : { /* multiply by y[0]^dp [beware dummy vars from FpX_FpXY_resultant] */
1713 1967 : GEN t = Flx_powu_pre(gel(y,0), dp, p, pi);
1714 7867 : for (i=2; i<lx; i++) gel(x,i) = Flx_mul_pre(gel(x,i), t, p, pi);
1715 : }
1716 7528 : return gc_GEN(av, x);
1717 : }
1718 :
1719 : /* return a Flx */
1720 : GEN
1721 2516 : FlxX_resultant(GEN u, GEN v, ulong p, long sx)
1722 : {
1723 2516 : pari_sp av = avma, av2;
1724 : long degq, dx, dy, du, dv, dr, signh;
1725 : ulong pi;
1726 : GEN z, g, h, r, p1;
1727 :
1728 2516 : dx = degpol(u); dy = degpol(v); signh = 1;
1729 2518 : if (dx < dy)
1730 : {
1731 7 : swap(u,v); lswap(dx,dy);
1732 7 : if (both_odd(dx, dy)) signh = -signh;
1733 : }
1734 2518 : if (dy < 0) return zero_Flx(sx);
1735 2518 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1736 2518 : if (dy==0) return gc_upto(av, Flx_powu_pre(gel(v,2),dx,p,pi));
1737 :
1738 2518 : g = h = pol1_Flx(sx); av2 = avma;
1739 : for(;;)
1740 : {
1741 7523 : r = FlxX_pseudorem(u,v,p,pi); dr = lg(r);
1742 7531 : if (dr == 2) { set_avma(av); return zero_Flx(sx); }
1743 7531 : du = degpol(u); dv = degpol(v); degq = du-dv;
1744 7529 : u = v; p1 = g; g = leading_coeff(u);
1745 7529 : switch(degq)
1746 : {
1747 0 : case 0: break;
1748 5549 : case 1:
1749 5549 : p1 = Flx_mul_pre(h,p1, p, pi); h = g; break;
1750 1980 : default:
1751 1980 : p1 = Flx_mul_pre(Flx_powu_pre(h,degq,p,pi), p1, p, pi);
1752 1977 : h = Flx_div_pre(Flx_powu_pre(g,degq,p,pi),
1753 1979 : Flx_powu_pre(h,degq-1,p,pi), p, pi);
1754 : }
1755 7521 : if (both_odd(du,dv)) signh = -signh;
1756 7521 : v = FlxY_Flx_div(r, p1, p);
1757 7523 : if (dr==3) break;
1758 5004 : if (gc_needed(av2,1))
1759 : {
1760 0 : if(DEBUGMEM>1) pari_warn(warnmem,"FlxX_resultant, dr = %ld",dr);
1761 0 : (void)gc_all(av2,4, &u, &v, &g, &h);
1762 : }
1763 : }
1764 2519 : z = gel(v,2);
1765 2519 : if (dv > 1) z = Flx_div_pre(Flx_powu_pre(z,dv,p,pi),
1766 0 : Flx_powu_pre(h,dv-1,p,pi), p, pi);
1767 2519 : if (signh < 0) z = Flx_neg(z,p);
1768 2519 : return gc_upto(av, z);
1769 : }
1770 :
1771 : /* Warning:
1772 : * This function switches between valid and invalid variable ordering*/
1773 :
1774 : static GEN
1775 6130 : FlxY_to_FlyX(GEN b, long sv)
1776 : {
1777 6130 : long i, n=-1;
1778 6130 : long sw = b[1]&VARNBITS;
1779 20914 : for(i=2;i<lg(b);i++) n = maxss(n,lgpol(gel(b,i)));
1780 6129 : return Flm_to_FlxX(Flm_transpose(FlxX_to_Flm(b,n)),sv,sw);
1781 : }
1782 :
1783 : /* Return a Fly*/
1784 : GEN
1785 6130 : Flx_FlxY_resultant(GEN a, GEN b, ulong p)
1786 : {
1787 6130 : pari_sp ltop=avma;
1788 6130 : long dres = degpol(a)*degpol(b);
1789 6130 : long sx=a[1], sy=b[1]&VARNBITS;
1790 : GEN z;
1791 6130 : b = FlxY_to_FlyX(b,sx);
1792 6129 : if ((ulong)dres >= p)
1793 2519 : z = FlxX_resultant(Fly_to_FlxY(a, sy), b, p, sx);
1794 : else
1795 : {
1796 3610 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1797 3610 : z = Flx_FlxY_resultant_polint(a, b, p, pi, (ulong)dres, sy);
1798 : }
1799 6130 : return gc_upto(ltop,z);
1800 : }
1801 :
1802 : /* Return a t_POL in variable vc whose coeffs are the coeffs of b in
1803 : * variable v; vc must have higher priority than all variables occuring in b. */
1804 : GEN
1805 146886 : swap_vars(GEN b, long v, long vc)
1806 : {
1807 146886 : long i, n = RgX_degree(b, v);
1808 : GEN c, x;
1809 146885 : if (n < 0) return pol_0(vc);
1810 146885 : c = cgetg(n+3, t_POL); x = c + 2;
1811 146885 : c[1] = evalsigne(1) | evalvarn(vc);
1812 971257 : for (i = 0; i <= n; i++) gel(x,i) = polcoef_i(b, i, v);
1813 146882 : return c;
1814 : }
1815 :
1816 : /* assume varn(b) << varn(a) */
1817 : /* return a FpY*/
1818 : GEN
1819 15 : FpX_FpXY_resultant(GEN a, GEN b, GEN p)
1820 : {
1821 15 : long i,n,dres, db, vY = varn(b), vX = varn(a);
1822 : GEN la,x,y;
1823 :
1824 15 : if (lgefint(p) == 3)
1825 : {
1826 0 : ulong pp = uel(p,2);
1827 0 : b = ZXX_to_FlxX(b, pp, vX);
1828 0 : a = ZX_to_Flx(a, pp);
1829 0 : x = Flx_FlxY_resultant(a, b, pp);
1830 0 : return Flx_to_ZX(x);
1831 : }
1832 15 : db = RgXY_degreex(b);
1833 15 : dres = degpol(a)*degpol(b);
1834 15 : la = leading_coeff(a);
1835 15 : x = cgetg(dres+2, t_VEC);
1836 15 : y = cgetg(dres+2, t_VEC);
1837 : /* Evaluate at dres+ 1 points: 0 (if dres even) and +/- n, so that P_n(X) =
1838 : * P_{-n}(-X), where P_i is Lagrange polynomial: P_i(j) = delta_{i,j} */
1839 157 : for (i=0,n = 1; i < dres; n++)
1840 : {
1841 142 : gel(x,++i) = utoipos(n);
1842 142 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1843 142 : gel(x,++i) = subiu(p,n);
1844 142 : gel(y,i) = FpX_FpXY_eval_resultant(a,b,gel(x,i),p,la,db,vY);
1845 : }
1846 15 : if (i == dres)
1847 : {
1848 0 : gel(x,++i) = gen_0;
1849 0 : gel(y,i) = FpX_FpXY_eval_resultant(a,b, gel(x,i), p,la,db,vY);
1850 : }
1851 15 : return FpV_polint(x,y, p, vY);
1852 : }
1853 :
1854 : GEN
1855 191 : FpX_composedsum(GEN P, GEN Q, GEN p)
1856 : {
1857 191 : pari_sp av = avma;
1858 191 : if (lgefint(p)==3)
1859 : {
1860 0 : ulong pp = p[2];
1861 0 : GEN z = Flx_composedsum(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
1862 0 : return gc_upto(av, Flx_to_ZX(z));
1863 : }
1864 : else
1865 : {
1866 191 : long n = 1+ degpol(P)*degpol(Q);
1867 191 : GEN Pl = FpX_invLaplace(FpX_Newton(P,n,p), p);
1868 191 : GEN Ql = FpX_invLaplace(FpX_Newton(Q,n,p), p);
1869 191 : GEN L = FpX_Laplace(FpXn_mul(Pl, Ql, n, p), p);
1870 191 : GEN lead = Fp_mul(Fp_powu(leading_coeff(P),degpol(Q), p),
1871 191 : Fp_powu(leading_coeff(Q),degpol(P), p), p);
1872 191 : GEN R = FpX_fromNewton(L, p);
1873 191 : return gc_upto(av, FpX_Fp_mul(R, lead, p));
1874 : }
1875 : }
1876 :
1877 : GEN
1878 0 : FpX_composedprod(GEN P, GEN Q, GEN p)
1879 : {
1880 0 : pari_sp av = avma;
1881 0 : if (lgefint(p)==3)
1882 : {
1883 0 : ulong pp = p[2];
1884 0 : GEN z = Flx_composedprod(ZX_to_Flx(P, pp), ZX_to_Flx(Q, pp), pp);
1885 0 : return gc_upto(av, Flx_to_ZX(z));
1886 : }
1887 : else
1888 : {
1889 0 : long n = 1+ degpol(P)*degpol(Q);
1890 0 : GEN L = FpX_convol(FpX_Newton(P,n,p), FpX_Newton(Q,n,p), p);
1891 0 : return gc_upto(av,FpX_fromNewton(L, p));
1892 : }
1893 : }
1894 :
1895 : static GEN
1896 191 : _FpX_composedsum(void *E, GEN a, GEN b)
1897 191 : { return FpX_composedsum(a,b, (GEN)E); }
1898 :
1899 : GEN
1900 1637 : FpXV_composedsum(GEN V, GEN p)
1901 : {
1902 1637 : if (lgefint(p)==3)
1903 : {
1904 0 : ulong pp = p[2];
1905 0 : return Flx_to_ZX(FlxV_composedsum(ZXV_to_FlxV(V, pp), pp));
1906 : }
1907 1637 : return gen_product(V, (void *)p, &_FpX_composedsum);
1908 : }
1909 :
1910 : /* 0, 1, -1, 2, -2, ... */
1911 : #define next_lambda(a) (a>0 ? -a : 1-a)
1912 :
1913 : /* Assume A in Z[Y], B in Q[Y][X], B squarefree in (Q[Y]/(A))[X] and
1914 : * Res_Y(A, B) in Z[X]. Find a small lambda (start from *lambda, use
1915 : * next_lambda successively) such that C(X) = Res_Y(A(Y), B(X + lambda Y))
1916 : * is squarefree, reset *lambda to the chosen value and return C. Set LERS to
1917 : * the Last nonconstant polynomial in the Euclidean Remainder Sequence */
1918 : static GEN
1919 22421 : ZX_ZXY_resultant_LERS(GEN A, GEN B0, long *plambda, GEN *LERS)
1920 : {
1921 : ulong bound, dp;
1922 22421 : pari_sp av = avma, av2 = 0;
1923 22421 : long lambda = *plambda, degA = degpol(A), dres = degA*degpol(B0);
1924 : long stable, checksqfree, i,n, cnt, degB;
1925 22421 : long v, vX = varn(B0), vY = varn(A); /* vY < vX */
1926 : GEN x, y, dglist, B, q, a, b, ev, H, H0, H1, Hp, H0p, H1p, C0, C1;
1927 : forprime_t S;
1928 :
1929 22421 : if (degA == 1)
1930 : {
1931 1260 : GEN a1 = gel(A,3), a0 = gel(A,2);
1932 1260 : B = lambda? RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY)): B0;
1933 1260 : H = gsubst(B, vY, gdiv(gneg(a0),a1));
1934 1260 : if (!equali1(a1)) H = RgX_Rg_mul(H, powiu(a1, poldegree(B,vY)));
1935 1260 : *LERS = mkvec2(scalarpol_shallow(a0,vX), scalarpol_shallow(a1,vX));
1936 1260 : return gc_all(av, 2, &H, LERS);
1937 : }
1938 :
1939 21161 : dglist = Hp = H0p = H1p = C0 = C1 = NULL; /* gcc -Wall */
1940 21161 : C0 = cgetg(dres+2, t_VECSMALL);
1941 21161 : C1 = cgetg(dres+2, t_VECSMALL);
1942 21161 : dglist = cgetg(dres+1, t_VECSMALL);
1943 21161 : x = cgetg(dres+2, t_VECSMALL);
1944 21161 : y = cgetg(dres+2, t_VECSMALL);
1945 21161 : B0 = leafcopy(B0);
1946 21161 : A = leafcopy(A);
1947 21161 : B = B0;
1948 21161 : v = fetch_var_higher(); setvarn(A,v);
1949 : /* make sure p large enough */
1950 22320 : INIT:
1951 : /* always except the first time */
1952 22320 : if (av2) { set_avma(av2); lambda = next_lambda(lambda); }
1953 22320 : if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
1954 22320 : B = swap_vars(B, vY, v);
1955 : /* B0(lambda v + x, v) */
1956 22320 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
1957 22320 : av2 = avma;
1958 :
1959 22320 : if (degA <= 3)
1960 : { /* sub-resultant faster for small degrees */
1961 11319 : H = RgX_resultant_all(A,B,&q);
1962 11319 : if (typ(q) != t_POL || degpol(q)!=1) goto INIT;
1963 10374 : H0 = gel(q,2);
1964 10374 : if (typ(H0) == t_POL) setvarn(H0,vX); else H0 = scalarpol(H0,vX);
1965 10374 : H1 = gel(q,3);
1966 10374 : if (typ(H1) == t_POL) setvarn(H1,vX); else H1 = scalarpol(H1,vX);
1967 10374 : if (!ZX_is_squarefree(H)) goto INIT;
1968 10332 : goto END;
1969 : }
1970 :
1971 11001 : H = H0 = H1 = NULL;
1972 11001 : degB = degpol(B);
1973 11001 : bound = ZX_ZXY_ResBound(A, B, NULL);
1974 11001 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
1975 11001 : dp = 1;
1976 11001 : init_modular_big(&S);
1977 11001 : for(cnt = 0, checksqfree = 1;;)
1978 49459 : {
1979 60460 : ulong p = u_forprime_next(&S);
1980 : GEN Hi;
1981 60460 : a = ZX_to_Flx(A, p);
1982 60461 : b = ZXX_to_FlxX(B, p, varn(A));
1983 60461 : if (degpol(a) < degA || degpol(b) < degB) continue; /* p | lc(A)lc(B) */
1984 60460 : if (checksqfree)
1985 : { /* find degree list for generic Euclidean Remainder Sequence */
1986 11001 : long goal = minss(degpol(a), degpol(b)); /* longest possible */
1987 74267 : for (n=1; n <= goal; n++) dglist[n] = 0;
1988 11001 : setlg(dglist, 1);
1989 24154 : for (n=0; n <= dres; n++)
1990 : {
1991 23741 : ev = FlxY_evalx_drop(b, n, p);
1992 23741 : Flx_resultant_set_dglist(a, ev, dglist, p);
1993 23740 : if (lg(dglist)-1 == goal) break;
1994 : }
1995 : /* last pol in ERS has degree > 1 ? */
1996 11000 : goal = lg(dglist)-1;
1997 11000 : if (degpol(B) == 1) { if (!goal) goto INIT; }
1998 : else
1999 : {
2000 10937 : if (goal <= 1) goto INIT;
2001 10853 : if (dglist[goal] != 0 || dglist[goal-1] != 1) goto INIT;
2002 : }
2003 10916 : if (DEBUGLEVEL>4)
2004 0 : err_printf("Degree list for ERS (trials: %ld) = %Ps\n",n+1,dglist);
2005 : }
2006 :
2007 2150223 : for (i=0,n = 0; i <= dres; n++)
2008 : {
2009 2089853 : ev = FlxY_evalx_drop(b, n, p);
2010 2089721 : x[++i] = n; y[i] = Flx_resultant_all(a, ev, C0+i, C1+i, dglist, p);
2011 2089848 : if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */
2012 : }
2013 60370 : Hi = Flv_Flm_polint(x, mkvec3(y,C0,C1), p, 0);
2014 60376 : Hp = gel(Hi,1); H0p = gel(Hi,2); H1p = gel(Hi,3);
2015 60376 : if (!H && degpol(Hp) != dres) continue;
2016 60376 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
2017 60376 : if (checksqfree) {
2018 10917 : if (!Flx_is_squarefree(Hp, p)) goto INIT;
2019 10829 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
2020 10829 : checksqfree = 0;
2021 : }
2022 :
2023 60288 : if (!H)
2024 : { /* initialize */
2025 10829 : q = utoipos(p); stable = 0;
2026 10829 : H = ZX_init_CRT(Hp, p,vX);
2027 10829 : H0= ZX_init_CRT(H0p, p,vX);
2028 10829 : H1= ZX_init_CRT(H1p, p,vX);
2029 : }
2030 : else
2031 : {
2032 49459 : GEN qp = muliu(q,p);
2033 49458 : stable = ZX_incremental_CRT_raw(&H, Hp, q,qp, p)
2034 49460 : & ZX_incremental_CRT_raw(&H0,H0p, q,qp, p)
2035 49460 : & ZX_incremental_CRT_raw(&H1,H1p, q,qp, p);
2036 49459 : q = qp;
2037 : }
2038 : /* could make it probabilistic for H ? [e.g if stable twice, etc]
2039 : * Probabilistic anyway for H0, H1 */
2040 60288 : if (DEBUGLEVEL>5 && (stable || ++cnt==100))
2041 0 : { cnt=0; err_printf("%ld%%%s ",100*expi(q)/bound,stable?"s":""); }
2042 60288 : if (stable && (ulong)expi(q) >= bound) break; /* DONE */
2043 49459 : if (gc_needed(av,2))
2044 : {
2045 0 : if (DEBUGMEM>1) pari_warn(warnmem,"ZX_ZXY_rnfequation");
2046 0 : (void)gc_all(av2, 4, &H, &q, &H0, &H1);
2047 : }
2048 : }
2049 21161 : END:
2050 21161 : if (DEBUGLEVEL>5) err_printf(" done\n");
2051 21161 : setvarn(H, vX); (void)delete_var();
2052 21161 : *LERS = mkvec2(H0,H1);
2053 21161 : *plambda = lambda; return gc_all(av, 2, &H, LERS);
2054 : }
2055 :
2056 : GEN
2057 60213 : ZX_ZXY_resultant_all(GEN A, GEN B, long *plambda, GEN *LERS)
2058 : {
2059 60213 : if (LERS)
2060 : {
2061 22421 : if (!plambda)
2062 0 : pari_err_BUG("ZX_ZXY_resultant_all [LERS != NULL needs lambda]");
2063 22421 : return ZX_ZXY_resultant_LERS(A, B, plambda, LERS);
2064 : }
2065 37792 : return ZX_ZXY_rnfequation(A, B, plambda);
2066 : }
2067 :
2068 : /* If lambda = NULL, return caract(Mod(A, T)), T,A in Z[X].
2069 : * Otherwise find a small lambda such that caract (Mod(A + lambda X, T)) is
2070 : * squarefree */
2071 : GEN
2072 22595 : ZXQ_charpoly_sqf(GEN A, GEN T, long *lambda, long v)
2073 : {
2074 22595 : pari_sp av = avma;
2075 : GEN R, a;
2076 : long dA;
2077 : int delvar;
2078 :
2079 22595 : if (v < 0) v = 0;
2080 22595 : switch (typ(A))
2081 : {
2082 22595 : case t_POL: dA = degpol(A); if (dA > 0) break;
2083 0 : A = constant_coeff(A);
2084 0 : default:
2085 0 : if (lambda) { A = scalar_ZX_shallow(A,varn(T)); dA = 0; break;}
2086 0 : return gc_upto(av, gpowgs(gsub(pol_x(v), A), degpol(T)));
2087 : }
2088 22595 : delvar = 0;
2089 22595 : if (varncmp(varn(T), 0) <= 0)
2090 : {
2091 3681 : long v0 = fetch_var(); delvar = 1;
2092 3681 : T = leafcopy(T); setvarn(T,v0);
2093 3681 : A = leafcopy(A); setvarn(A,v0);
2094 : }
2095 22595 : R = ZX_ZXY_rnfequation(T, deg1pol_shallow(gen_1, gneg_i(A), 0), lambda);
2096 22595 : if (delvar) (void)delete_var();
2097 22595 : setvarn(R, v); a = leading_coeff(T);
2098 22595 : if (!gequal1(a)) R = gdiv(R, powiu(a, dA));
2099 22595 : return gc_upto(av, R);
2100 : }
2101 :
2102 : /* charpoly(Mod(A,T)), A may be in Q[X], but assume T and result are integral */
2103 : GEN
2104 1247393 : ZXQ_charpoly(GEN A, GEN T, long v)
2105 : {
2106 1247393 : return (degpol(T) < 16) ? RgXQ_charpoly_i(A,T,v): ZXQ_charpoly_sqf(A,T, NULL, v);
2107 : }
2108 :
2109 : GEN
2110 9772 : QXQ_charpoly(GEN A, GEN T, long v)
2111 : {
2112 9772 : pari_sp av = avma;
2113 9772 : GEN den, B = Q_remove_denom(A, &den);
2114 9772 : GEN P = ZXQ_charpoly(B, T, v);
2115 9772 : return gc_GEN(av, den ? RgX_rescale(P, ginv(den)): P);
2116 : }
2117 :
2118 : static ulong
2119 3864927 : ZX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, ulong p)
2120 : {
2121 3864927 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2122 : ulong H, dp;
2123 3864772 : if (dropa && dropb) return 0; /* p | lc(A), p | lc(B) */
2124 3864772 : H = Flx_resultant(a, b, p);
2125 3864425 : if (dropa)
2126 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2127 0 : ulong c = b[degB+2]; /* lc(B) */
2128 0 : if (odd(degB)) c = p - c;
2129 0 : c = Fl_powu(c, dropa, p);
2130 0 : if (c != 1) H = Fl_mul(H, c, p);
2131 : }
2132 3864425 : else if (dropb)
2133 : { /* multiply by lc(A)^(deg B - deg b) */
2134 0 : ulong c = a[degA+2]; /* lc(A) */
2135 0 : c = Fl_powu(c, dropb, p);
2136 0 : if (c != 1) H = Fl_mul(H, c, p);
2137 : }
2138 3864428 : dp = dB ? umodiu(dB, p): 1;
2139 3864428 : if (dp != 1) H = Fl_mul(H, Fl_powu(Fl_inv(dp,p), degA, p), p);
2140 3864426 : return H;
2141 : }
2142 :
2143 : /* If B=NULL, assume B=A' */
2144 : static GEN
2145 1495367 : ZX_resultant_slice(GEN A, GEN B, GEN dB, GEN P, GEN *mod)
2146 : {
2147 1495367 : pari_sp av = avma, av2;
2148 1495367 : long degA, degB, i, n = lg(P)-1;
2149 : GEN H, T;
2150 :
2151 1495367 : degA = degpol(A);
2152 1495364 : degB = B? degpol(B): degA - 1;
2153 1495364 : if (n == 1)
2154 : {
2155 811362 : ulong Hp, p = uel(P,1);
2156 811362 : GEN a = ZX_to_Flx(A, p), b = B? ZX_to_Flx(B, p): Flx_deriv(a, p);
2157 811372 : Hp = ZX_resultant_prime(a, b, dB, degA, degB, p);
2158 811352 : set_avma(av); *mod = utoipos(p); return utoi(Hp);
2159 : }
2160 684002 : T = ZV_producttree(P);
2161 684002 : A = ZX_nv_mod_tree(A, P, T);
2162 684000 : if (B) B = ZX_nv_mod_tree(B, P, T);
2163 684000 : H = cgetg(n+1, t_VECSMALL); av2 = avma;
2164 3737174 : for(i=1; i <= n; i++, set_avma(av2))
2165 : {
2166 3053181 : ulong p = P[i];
2167 3053181 : GEN a = gel(A,i), b = B? gel(B,i): Flx_deriv(a, p);
2168 3053562 : H[i] = ZX_resultant_prime(a, b, dB, degA, degB, p);
2169 : }
2170 683993 : H = ZV_chinese_tree(H, P, T, ZV_chinesetree(P,T));
2171 684002 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2172 : }
2173 :
2174 : GEN
2175 1495380 : ZX_resultant_worker(GEN P, GEN A, GEN B, GEN dB)
2176 : {
2177 1495380 : GEN V = cgetg(3, t_VEC);
2178 1495368 : if (typ(B) == t_INT) B = NULL;
2179 1495368 : if (!signe(dB)) dB = NULL;
2180 1495368 : gel(V,1) = ZX_resultant_slice(A, B, dB, P, &gel(V,2));
2181 1495369 : return V;
2182 : }
2183 :
2184 : /* Compute Res(A, B/dB) in Z, assuming A,B in Z[X], dB in Z or NULL (= 1)
2185 : * If B=NULL, take B = A' and assume deg A > 1 and 'bound' is set */
2186 : GEN
2187 1350924 : ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound)
2188 : {
2189 1350924 : pari_sp av = avma;
2190 : forprime_t S;
2191 : GEN H, worker;
2192 1350924 : if (!B && degpol(A)==2)
2193 : {
2194 114097 : GEN a = gel(A,4), b = gel(A,3), c = gel(A,2);
2195 114097 : H = mulii(a, subii(shifti(mulii(a, c), 2), sqri(b)));
2196 114088 : if (dB) H = diviiexact(H, sqri(dB));
2197 114088 : return gc_INT(av, H);
2198 : }
2199 1236827 : if (B)
2200 : {
2201 154024 : long a = degpol(A), b = degpol(B);
2202 154024 : if (a < 0 || b < 0) return gen_0;
2203 153994 : if (!a) return powiu(gel(A,2), b);
2204 153994 : if (!b) return powiu(gel(B,2), a);
2205 153369 : if (minss(a, b) <= 1)
2206 : {
2207 76623 : H = RgX_resultant_all(A, B, NULL);
2208 76623 : if (dB) H = diviiexact(H, powiu(dB, a));
2209 76623 : return gc_INT(av, H);
2210 : }
2211 76746 : if (!bound) bound = ZX_ZXY_ResBound(A, B, dB);
2212 : }
2213 1159557 : worker = snm_closure(is_entry("_ZX_resultant_worker"),
2214 : mkvec3(A, B? B: gen_0, dB? dB: gen_0));
2215 1159616 : init_modular_big(&S);
2216 1159571 : H = gen_crt("ZX_resultant_all", worker, &S, dB, bound, 0, NULL,
2217 : ZV_chinese_center, Fp_center);
2218 1159603 : return gc_INT(av, H);
2219 : }
2220 :
2221 : /* A0 and B0 in Q[X] */
2222 : GEN
2223 56 : QX_resultant(GEN A0, GEN B0)
2224 : {
2225 : GEN s, a, b, A, B;
2226 56 : pari_sp av = avma;
2227 :
2228 56 : A = Q_primitive_part(A0, &a);
2229 56 : B = Q_primitive_part(B0, &b);
2230 56 : s = ZX_resultant(A, B);
2231 56 : if (!signe(s)) { set_avma(av); return gen_0; }
2232 56 : if (a) s = gmul(s, gpowgs(a,degpol(B)));
2233 56 : if (b) s = gmul(s, gpowgs(b,degpol(A)));
2234 56 : return gc_upto(av, s);
2235 : }
2236 :
2237 : GEN
2238 56077 : ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); }
2239 :
2240 : GEN
2241 0 : QXQ_intnorm(GEN A, GEN B)
2242 : {
2243 : GEN c, n, R, lB;
2244 0 : long dA = degpol(A), dB = degpol(B);
2245 0 : pari_sp av = avma;
2246 0 : if (dA < 0) return gen_0;
2247 0 : A = Q_primitive_part(A, &c);
2248 0 : if (!c || typ(c) == t_INT) {
2249 0 : n = c;
2250 0 : R = ZX_resultant(B, A);
2251 : } else {
2252 0 : n = gel(c,1);
2253 0 : R = ZX_resultant_all(B, A, gel(c,2), 0);
2254 : }
2255 0 : if (n && !equali1(n)) R = mulii(R, powiu(n, dB));
2256 0 : lB = leading_coeff(B);
2257 0 : if (!equali1(lB)) R = diviiexact(R, powiu(lB, dA));
2258 0 : return gc_INT(av, R);
2259 : }
2260 :
2261 : GEN
2262 19418 : QXQ_norm(GEN A, GEN B)
2263 : {
2264 : GEN c, R, lB;
2265 19418 : long dA = degpol(A), dB = degpol(B);
2266 19418 : pari_sp av = avma;
2267 19418 : if (dA < 0) return gen_0;
2268 19418 : A = Q_primitive_part(A, &c);
2269 19418 : R = ZX_resultant(B, A);
2270 19418 : if (c) R = gmul(R, gpowgs(c, dB));
2271 19418 : lB = leading_coeff(B);
2272 19418 : if (!equali1(lB)) R = gdiv(R, gpowgs(lB, dA));
2273 19418 : return gc_upto(av, R);
2274 : }
2275 :
2276 : /* assume x has integral coefficients */
2277 : GEN
2278 1200237 : ZX_disc_all(GEN x, ulong bound)
2279 : {
2280 1200237 : pari_sp av = avma;
2281 1200237 : long s, d = degpol(x);
2282 : GEN l, R;
2283 :
2284 1200237 : if (d <= 1) return d == 1? gen_1: gen_0;
2285 1196944 : s = (d & 2) ? -1: 1;
2286 1196944 : l = leading_coeff(x);
2287 1196943 : if (!bound) bound = ZX_discbound(x);
2288 1196873 : R = ZX_resultant_all(x, NULL, NULL, bound);
2289 1196911 : if (is_pm1(l))
2290 1017877 : { if (signe(l) < 0) s = -s; }
2291 : else
2292 179037 : R = diviiexact(R,l);
2293 1196914 : if (s == -1) togglesign_safe(&R);
2294 1196911 : return gc_INT(av,R);
2295 : }
2296 :
2297 : GEN
2298 1138168 : ZX_disc(GEN x) { return ZX_disc_all(x,0); }
2299 :
2300 : static GEN
2301 11008 : ZXQX_resultant_prime(GEN a, GEN b, GEN dB, long degA, long degB, GEN T, ulong p)
2302 : {
2303 11008 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2304 : GEN H, dp;
2305 11008 : if (dropa && dropb) return pol0_Flx(T[1]); /* p | lc(A), p | lc(B) */
2306 11008 : H = FlxqX_saferesultant(a, b, T, p);
2307 11009 : if (!H) return NULL;
2308 11009 : if (dropa)
2309 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2310 0 : GEN c = gel(b,degB+2); /* lc(B) */
2311 0 : if (odd(degB)) c = Flx_neg(c, p);
2312 0 : c = Flxq_powu(c, dropa, T, p);
2313 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2314 : }
2315 11009 : else if (dropb)
2316 : { /* multiply by lc(A)^(deg B - deg b) */
2317 0 : GEN c = gel(a,degA+2); /* lc(A) */
2318 0 : c = Flxq_powu(c, dropb, T, p);
2319 0 : if (!Flx_equal1(c)) H = Flxq_mul(H, c, T, p);
2320 : }
2321 11009 : dp = dB ? ZX_to_Flx(dB, p): pol1_Flx(T[1]);
2322 11009 : if (!Flx_equal1(dp))
2323 : {
2324 0 : GEN idp = Flxq_invsafe(dp, T, p);
2325 0 : if (!idp) return NULL;
2326 0 : H = Flxq_mul(H, Flxq_powu(idp, degA, T, p), T, p);
2327 : }
2328 11009 : return H;
2329 : }
2330 :
2331 : /* If B=NULL, assume B=A' */
2332 : static GEN
2333 4911 : ZXQX_resultant_slice(GEN A, GEN B, GEN U, GEN dB, GEN P, GEN *mod)
2334 : {
2335 4911 : pari_sp av = avma;
2336 4911 : long degA, degB, i, n = lg(P)-1;
2337 : GEN H, T;
2338 4911 : long v = varn(U), redo = 0;
2339 :
2340 4911 : degA = degpol(A);
2341 4911 : degB = B? degpol(B): degA - 1;
2342 4911 : if (n == 1)
2343 : {
2344 3177 : ulong p = uel(P,1);
2345 3177 : GEN a = ZXX_to_FlxX(A, p, v), b = B? ZXX_to_FlxX(B, p, v): FlxX_deriv(a, p);
2346 3177 : GEN u = ZX_to_Flx(U, p);
2347 3177 : GEN Hp = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2348 3177 : if (!Hp) { set_avma(av); *mod = gen_1; return pol_0(v); }
2349 3177 : Hp = gc_upto(av, Flx_to_ZX(Hp)); *mod = utoipos(p); return Hp;
2350 : }
2351 1734 : T = ZV_producttree(P);
2352 1734 : A = ZXX_nv_mod_tree(A, P, T, v);
2353 1734 : if (B) B = ZXX_nv_mod_tree(B, P, T, v);
2354 1734 : U = ZX_nv_mod_tree(U, P, T);
2355 1734 : H = cgetg(n+1, t_VEC);
2356 9567 : for(i=1; i <= n; i++)
2357 : {
2358 7833 : ulong p = P[i];
2359 7833 : GEN a = gel(A,i), b = B? gel(B,i): FlxX_deriv(a, p), u = gel(U, i);
2360 7831 : GEN h = ZXQX_resultant_prime(a, b, dB, degA, degB, u, p);
2361 7832 : if (!h)
2362 : {
2363 0 : gel(H,i) = pol_0(v);
2364 0 : P[i] = 1; redo = 1;
2365 : }
2366 : else
2367 7832 : gel(H,i) = h;
2368 : }
2369 1734 : if (redo) T = ZV_producttree(P);
2370 1734 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2371 1734 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2372 : }
2373 :
2374 : GEN
2375 4911 : ZXQX_resultant_worker(GEN P, GEN A, GEN B, GEN T, GEN dB)
2376 : {
2377 4911 : GEN V = cgetg(3, t_VEC);
2378 4911 : if (isintzero(B)) B = NULL;
2379 4911 : if (!signe(dB)) dB = NULL;
2380 4911 : gel(V,1) = ZXQX_resultant_slice(A, B, T, dB, P, &gel(V,2));
2381 4911 : return V;
2382 : }
2383 :
2384 : static ulong
2385 4315 : ZXQX_resultant_bound_i(GEN nf, GEN A, GEN B, GEN (*f)(GEN,GEN,long))
2386 : {
2387 4315 : pari_sp av = avma;
2388 4315 : GEN r, M = nf_L2_bound(nf, NULL, &r);
2389 4315 : long v = nf_get_varn(nf), i, l = lg(r);
2390 4315 : GEN a = cgetg(l, t_COL);
2391 13334 : for (i = 1; i < l; i++)
2392 9019 : gel(a, i) = f(gsubst(A, v, gel(r,i)), gsubst(B, v, gel(r,i)), DEFAULTPREC);
2393 4315 : return gc_ulong(av, (ulong) dbllog2(gmul(M,RgC_fpnorml2(a, DEFAULTPREC))));
2394 : }
2395 : static ulong
2396 4000 : ZXQX_resultant_bound(GEN nf, GEN A, GEN B)
2397 4000 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound); }
2398 :
2399 : static GEN
2400 56 : _ZXQ_powu(GEN x, ulong u, GEN T)
2401 56 : { return typ(x) == t_INT? powiu(x, u): ZXQ_powu(x, u, T); }
2402 :
2403 : /* Compute Res(A, B/dB) in Z[X]/T, assuming A,B in Z[X,Y], dB in Z or NULL (= 1)
2404 : * If B=NULL, take B = A' and assume deg A > 1 */
2405 : static GEN
2406 3997 : ZXQX_resultant_all(GEN A, GEN B, GEN T, GEN dB, ulong bound)
2407 : {
2408 3997 : pari_sp av = avma;
2409 : forprime_t S;
2410 : GEN H, worker;
2411 3997 : if (B)
2412 : {
2413 63 : long a = degpol(A), b = degpol(B);
2414 63 : if (a < 0 || b < 0) return gen_0;
2415 63 : if (!a) return _ZXQ_powu(gel(A,2), b, T);
2416 63 : if (!b) return _ZXQ_powu(gel(B,2), a, T);
2417 : } else
2418 3934 : if (!bound) B = RgX_deriv(A);
2419 3997 : if (!bound) bound = ZXQX_resultant_bound(nfinit(T, DEFAULTPREC), A, B);
2420 3997 : worker = snm_closure(is_entry("_ZXQX_resultant_worker"),
2421 : mkvec4(A, B? B: gen_0, T, dB? dB: gen_0));
2422 3997 : init_modular_big(&S);
2423 3997 : H = gen_crt("ZXQX_resultant_all", worker, &S, dB, bound, 0, NULL,
2424 : nxV_chinese_center, FpX_center);
2425 3997 : if (DEBUGLEVEL)
2426 0 : err_printf("ZXQX_resultant_all: a priori bound: %lu, a posteriori: %lu\n",
2427 : bound, expi(gsupnorm(H, DEFAULTPREC)));
2428 3997 : return gc_upto(av, H);
2429 : }
2430 :
2431 : GEN
2432 119 : nfX_resultant(GEN nf, GEN x, GEN y)
2433 : {
2434 119 : pari_sp av = avma;
2435 119 : GEN cx, cy, D, T = nf_get_pol(nf);
2436 119 : long dx = degpol(x), dy = degpol(y);
2437 119 : if (dx < 0 || dy < 0) return gen_0;
2438 119 : x = Q_primitive_part(x, &cx); if (cx) cx = gpowgs(cx, dy);
2439 119 : y = Q_primitive_part(y, &cy); if (cy) cy = gpowgs(cy, dx);
2440 119 : if (!dx) D = _ZXQ_powu(gel(x,2), dy, T);
2441 119 : else if (!dy) D = _ZXQ_powu(gel(y,2), dx, T);
2442 : else
2443 : {
2444 63 : ulong bound = ZXQX_resultant_bound(nf, x, y);
2445 63 : D = ZXQX_resultant_all(x, y, T, NULL, bound);
2446 : }
2447 119 : cx = mul_content(cx, cy); if (cx) D = gmul(D, cx);
2448 119 : return gc_upto(av, D);
2449 : }
2450 :
2451 : static GEN
2452 252 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2453 :
2454 : static GEN
2455 3934 : ZXQX_disc_all(GEN x, GEN T, ulong bound)
2456 : {
2457 3934 : pari_sp av = avma;
2458 3934 : long s, d = degpol(x), v = varn(T);
2459 : GEN l, R;
2460 :
2461 3934 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2462 3934 : s = (d & 2) ? -1: 1;
2463 3934 : l = leading_coeff(x);
2464 3934 : R = ZXQX_resultant_all(x, NULL, T, NULL, bound);
2465 3934 : if (!gequal1(l)) R = QXQ_div(R, to_ZX(l,v), T);
2466 3934 : if (s == -1) R = RgX_neg(R);
2467 3934 : return gc_upto(av, R);
2468 : }
2469 :
2470 : GEN
2471 7 : QX_disc(GEN x)
2472 : {
2473 7 : pari_sp av = avma;
2474 7 : GEN c, d = ZX_disc( Q_primitive_part(x, &c) );
2475 7 : if (c) d = gmul(d, gpowgs(c, 2*degpol(x) - 2));
2476 7 : return gc_upto(av, d);
2477 : }
2478 :
2479 : GEN
2480 4165 : nfX_disc(GEN nf, GEN x)
2481 : {
2482 4165 : pari_sp av = avma;
2483 4165 : GEN c, D, T = nf_get_pol(nf);
2484 : ulong bound;
2485 4165 : long d = degpol(x), v = varn(T);
2486 4165 : if (d <= 1) return d == 1? pol_1(v): pol_0(v);
2487 3934 : x = Q_primitive_part(x, &c);
2488 3934 : bound = ZXQX_resultant_bound(nf, x, RgX_deriv(x));
2489 3934 : D = ZXQX_disc_all(x, T, bound);
2490 3934 : if (c) D = gmul(D, gpowgs(c, 2*d - 2));
2491 3934 : return gc_upto(av, D);
2492 : }
2493 :
2494 : GEN
2495 846037 : QXQ_mul(GEN x, GEN y, GEN T)
2496 : {
2497 846037 : GEN dx, nx = Q_primitive_part(x, &dx);
2498 846038 : GEN dy, ny = Q_primitive_part(y, &dy);
2499 846040 : GEN z = ZXQ_mul(nx, ny, T);
2500 846038 : if (dx || dy)
2501 : {
2502 843240 : GEN d = dx ? dy ? gmul(dx, dy): dx : dy;
2503 843240 : if (!gequal1(d)) z = ZX_Q_mul(z, d);
2504 : }
2505 846040 : return z;
2506 : }
2507 :
2508 : GEN
2509 407040 : QXQ_sqr(GEN x, GEN T)
2510 : {
2511 407040 : GEN dx, nx = Q_primitive_part(x, &dx);
2512 407040 : GEN z = ZXQ_sqr(nx, T);
2513 407040 : if (dx)
2514 405304 : z = ZX_Q_mul(z, gsqr(dx));
2515 407040 : return z;
2516 : }
2517 :
2518 : static GEN
2519 212913 : QXQ_inv_slice(GEN A, GEN B, GEN P, GEN *mod)
2520 : {
2521 212913 : pari_sp av = avma;
2522 212913 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2523 : GEN H, T;
2524 212913 : if (n == 1)
2525 : {
2526 165825 : ulong p = uel(P,1);
2527 165825 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2528 165825 : GEN U = Flxq_invsafe(a, b, p);
2529 165825 : if (!U)
2530 : {
2531 24 : set_avma(av);
2532 24 : *mod = gen_1; return pol_0(v);
2533 : }
2534 165801 : H = gc_GEN(av, Flx_to_ZX(U));
2535 165801 : *mod = utoipos(p); return H;
2536 : }
2537 47088 : T = ZV_producttree(P);
2538 47088 : A = ZX_nv_mod_tree(A, P, T);
2539 47087 : B = ZX_nv_mod_tree(B, P, T);
2540 47087 : H = cgetg(n+1, t_VEC);
2541 238080 : for(i=1; i <= n; i++)
2542 : {
2543 190993 : ulong p = P[i];
2544 190993 : GEN a = gel(A,i), b = gel(B,i);
2545 190993 : GEN U = Flxq_invsafe(a, b, p);
2546 190992 : if (!U)
2547 : {
2548 601 : gel(H,i) = pol_0(v);
2549 601 : P[i] = 1; redo = 1;
2550 : }
2551 : else
2552 190391 : gel(H,i) = U;
2553 : }
2554 47087 : if (redo) T = ZV_producttree(P);
2555 47087 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2556 47088 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2557 : }
2558 :
2559 : GEN
2560 212913 : QXQ_inv_worker(GEN P, GEN A, GEN B)
2561 : {
2562 212913 : GEN V = cgetg(3, t_VEC);
2563 212913 : gel(V,1) = QXQ_inv_slice(A, B, P, &gel(V,2));
2564 212913 : return V;
2565 : }
2566 :
2567 : /* lift(1 / Mod(A,B)). B a ZX, A a scalar or a QX */
2568 : GEN
2569 146112 : QXQ_inv(GEN A, GEN B)
2570 : {
2571 : GEN D, Ap, Bp;
2572 : ulong pp;
2573 146112 : pari_sp av2, av = avma;
2574 : forprime_t S;
2575 146112 : GEN worker, U, H = NULL, mod = gen_1;
2576 : pari_timer ti;
2577 : long k, dA, dB;
2578 146112 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2579 : /* A a QX, B a ZX */
2580 146112 : A = Q_primitive_part(A, &D);
2581 146112 : dA = degpol(A); dB= degpol(B);
2582 : /* A, B in Z[X] */
2583 146112 : init_modular_small(&S);
2584 : do {
2585 146112 : pp = u_forprime_next(&S);
2586 146112 : Ap = ZX_to_Flx(A, pp);
2587 146112 : Bp = ZX_to_Flx(B, pp);
2588 146112 : } while (degpol(Ap) != dA || degpol(Bp) != dB);
2589 146112 : if (degpol(Flx_gcd(Ap, Bp, pp)) != 0 && degpol(ZX_gcd(A,B))!=0)
2590 14 : pari_err_INV("QXQ_inv",mkpolmod(A,B));
2591 146098 : worker = snm_closure(is_entry("_QXQ_inv_worker"), mkvec2(A, B));
2592 146098 : av2 = avma;
2593 146098 : for (k = 1; ;k *= 2)
2594 42591 : {
2595 : GEN res, b, N, den;
2596 188689 : gen_inccrt_i("QXQ_inv", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2597 : nxV_chinese_center, FpX_center);
2598 188689 : (void)gc_all(av2, 2, &H, &mod);
2599 188689 : b = sqrti(shifti(mod,-1));
2600 188689 : if (DEBUGLEVEL>5) timer_start(&ti);
2601 188689 : U = FpX_ratlift(H, mod, b, b, NULL);
2602 188689 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: ratlift");
2603 194435 : if (!U) continue;
2604 151844 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2605 151844 : res = Flx_rem(Flx_Fl_sub(Flx_mul(Ap, ZX_to_Flx(N,pp), pp),
2606 : umodiu(den, pp), pp), Bp, pp);
2607 151844 : if (degpol(res) >= 0) continue;
2608 146098 : res = ZX_Z_sub(ZX_mul(A, N), den);
2609 146098 : res = ZX_is_monic(B) ? ZX_rem(res, B): RgX_pseudorem(res, B);
2610 146098 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_inv: final check");
2611 146098 : if (degpol(res)<0)
2612 : {
2613 146098 : if (D) U = RgX_Rg_div(U, D);
2614 146098 : return gc_GEN(av, U);
2615 : }
2616 : }
2617 : }
2618 :
2619 : static GEN
2620 121107 : QXQ_div_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
2621 : {
2622 121107 : pari_sp av = avma;
2623 121107 : long i, n = lg(P)-1, v = varn(A), redo = 0;
2624 : GEN H, T;
2625 121107 : if (n == 1)
2626 : {
2627 44701 : ulong p = uel(P,1);
2628 44701 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p), c = ZX_to_Flx(C, p);
2629 44701 : GEN bi = Flxq_invsafe(b, c, p), U;
2630 44701 : if (!bi)
2631 : {
2632 0 : set_avma(av);
2633 0 : *mod = gen_1; return pol_0(v);
2634 : }
2635 44701 : U = Flxq_mul(a, bi, c, p);
2636 44701 : H = gc_GEN(av, Flx_to_ZX(U));
2637 44701 : *mod = utoipos(p); return H;
2638 : }
2639 76406 : T = ZV_producttree(P);
2640 76406 : A = ZX_nv_mod_tree(A, P, T);
2641 76406 : B = ZX_nv_mod_tree(B, P, T);
2642 76406 : C = ZX_nv_mod_tree(C, P, T);
2643 76406 : H = cgetg(n+1, t_VEC);
2644 337441 : for(i=1; i <= n; i++)
2645 : {
2646 261035 : ulong p = P[i];
2647 261035 : GEN a = gel(A,i), b = gel(B,i), c = gel(C, i);
2648 261035 : GEN bi = Flxq_invsafe(b, c, p);
2649 261041 : if (!bi)
2650 : {
2651 4 : gel(H,i) = pol_0(v);
2652 4 : P[i] = 1; redo = 1;
2653 : }
2654 : else
2655 261037 : gel(H,i) = Flxq_mul(a, bi, c, p);
2656 : }
2657 76406 : if (redo) T = ZV_producttree(P);
2658 76406 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2659 76406 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2660 : }
2661 :
2662 : GEN
2663 121107 : QXQ_div_worker(GEN P, GEN A, GEN B, GEN C)
2664 : {
2665 121107 : GEN V = cgetg(3, t_VEC);
2666 121107 : gel(V,1) = QXQ_div_slice(A, B, C, P, &gel(V,2));
2667 121107 : return V;
2668 : }
2669 :
2670 : /* lift(Mod(A/B, C)). C a ZX, A, B a scalar or a QX */
2671 : GEN
2672 33207 : QXQ_div(GEN A, GEN B, GEN C)
2673 : {
2674 : GEN DA, DB, Ap, Bp, Cp;
2675 : ulong pp;
2676 33207 : pari_sp av2, av = avma;
2677 : forprime_t S;
2678 33207 : GEN worker, U, H = NULL, mod = gen_1;
2679 : pari_timer ti;
2680 : long k, dA, dB, dC;
2681 33207 : if (is_scalar_t(typ(A))) return scalarpol(ginv(A), varn(B));
2682 : /* A a QX, B a ZX */
2683 33207 : A = Q_primitive_part(A, &DA);
2684 33207 : B = Q_primitive_part(B, &DB);
2685 33207 : dA = degpol(A); dB = degpol(B); dC = degpol(C);
2686 : /* A, B in Z[X] */
2687 33207 : init_modular_small(&S);
2688 : do {
2689 33207 : pp = u_forprime_next(&S);
2690 33207 : Ap = ZX_to_Flx(A, pp);
2691 33207 : Bp = ZX_to_Flx(B, pp);
2692 33207 : Cp = ZX_to_Flx(C, pp);
2693 33207 : } while (degpol(Ap) != dA || degpol(Bp) != dB || degpol(Cp) != dC);
2694 33207 : if (degpol(Flx_gcd(Bp, Cp, pp)) != 0 && degpol(ZX_gcd(B,C))!=0)
2695 0 : pari_err_INV("QXQ_div",mkpolmod(B,C));
2696 33207 : worker = snm_closure(is_entry("_QXQ_div_worker"), mkvec3(A, B, C));
2697 33207 : av2 = avma;
2698 33207 : for (k = 1; ;k *= 2)
2699 46720 : {
2700 : GEN res, b, N, den;
2701 79927 : gen_inccrt_i("QXQ_div", worker, NULL, (k+1)>>1, 0, &S, &H, &mod,
2702 : nxV_chinese_center, FpX_center);
2703 79927 : (void)gc_all(av2, 2, &H, &mod);
2704 79927 : b = sqrti(shifti(mod,-1));
2705 79927 : if (DEBUGLEVEL>5) timer_start(&ti);
2706 79927 : U = FpX_ratlift(H, mod, b, b, NULL);
2707 79927 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: ratlift");
2708 90556 : if (!U) continue;
2709 43836 : N = Q_remove_denom(U, &den); if (!den) den = gen_1;
2710 43836 : res = Flx_rem(Flx_sub(Flx_mul(Bp, ZX_to_Flx(N,pp), pp),
2711 : Flx_Fl_mul(Ap, umodiu(den, pp), pp), pp), Cp, pp);
2712 43836 : if (degpol(res) >= 0) continue;
2713 33207 : res = ZX_sub(ZX_mul(B, N), ZX_Z_mul(A,den));
2714 33207 : res = ZX_is_monic(C) ? ZX_rem(res, C): RgX_pseudorem(res, C);
2715 33207 : if (DEBUGLEVEL>5) timer_printf(&ti,"QXQ_div: final check");
2716 33207 : if (degpol(res)<0)
2717 : {
2718 33207 : if (DA && DB) U = RgX_Rg_mul(U, gdiv(DA,DB));
2719 28069 : else if (DA) U = RgX_Rg_mul(U, DA);
2720 15981 : else if (DB) U = RgX_Rg_div(U, DB);
2721 33207 : return gc_GEN(av, U);
2722 : }
2723 : }
2724 : }
2725 :
2726 : /************************************************************************
2727 : * *
2728 : * ZXQ_minpoly *
2729 : * *
2730 : ************************************************************************/
2731 :
2732 : static GEN
2733 3523 : ZXQ_minpoly_slice(GEN A, GEN B, long d, GEN P, GEN *mod)
2734 : {
2735 3523 : pari_sp av = avma;
2736 3523 : long i, n = lg(P)-1, v = evalvarn(varn(B));
2737 : GEN H, T;
2738 3523 : if (n == 1)
2739 : {
2740 716 : ulong p = uel(P,1);
2741 716 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2742 716 : GEN Hp = Flxq_minpoly(a, b, p);
2743 716 : if (degpol(Hp) != d) { p = 1; Hp = pol0_Flx(v); }
2744 716 : H = gc_upto(av, Flx_to_ZX(Hp));
2745 716 : *mod = utoipos(p); return H;
2746 : }
2747 2807 : T = ZV_producttree(P);
2748 2807 : A = ZX_nv_mod_tree(A, P, T);
2749 2807 : B = ZX_nv_mod_tree(B, P, T);
2750 2807 : H = cgetg(n+1, t_VEC);
2751 16838 : for(i=1; i <= n; i++)
2752 : {
2753 14031 : ulong p = P[i];
2754 14031 : GEN a = gel(A,i), b = gel(B,i);
2755 14031 : GEN m = Flxq_minpoly(a, b, p);
2756 14031 : if (degpol(m) != d) { P[i] = 1; m = pol0_Flx(v); }
2757 14031 : gel(H, i) = m;
2758 : }
2759 2807 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2760 2807 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2761 : }
2762 :
2763 : GEN
2764 3523 : ZXQ_minpoly_worker(GEN P, GEN A, GEN B, long d)
2765 : {
2766 3523 : GEN V = cgetg(3, t_VEC);
2767 3523 : gel(V,1) = ZXQ_minpoly_slice(A, B, d, P, &gel(V,2));
2768 3523 : return V;
2769 : }
2770 :
2771 : GEN
2772 1701 : ZXQ_minpoly(GEN A, GEN B, long d, ulong bound)
2773 : {
2774 1701 : pari_sp av = avma;
2775 : GEN worker, H, dB;
2776 : forprime_t S;
2777 1701 : B = Q_remove_denom(B, &dB);
2778 1701 : worker = strtoclosure("_ZXQ_minpoly_worker", 3, A, B, stoi(d));
2779 1701 : init_modular_big(&S);
2780 1701 : H = gen_crt("ZXQ_minpoly", worker, &S, dB, bound, 0, NULL,
2781 : nxV_chinese_center, FpX_center_i);
2782 1701 : return gc_GEN(av, H);
2783 : }
2784 :
2785 : /************************************************************************
2786 : * *
2787 : * ZX_ZXY_resultant *
2788 : * *
2789 : ************************************************************************/
2790 :
2791 : static GEN
2792 364908 : ZX_ZXY_resultant_prime(GEN a, GEN b, ulong dp, ulong p,
2793 : long degA, long degB, long dres, long sX)
2794 : {
2795 364908 : long dropa = degA - degpol(a), dropb = degB - degpol(b);
2796 364908 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2797 364907 : GEN Hp = Flx_FlxY_resultant_polint(a, b, p, pi, dres, sX);
2798 364906 : if (dropa && dropb)
2799 0 : Hp = zero_Flx(sX);
2800 : else {
2801 364906 : if (dropa)
2802 : { /* multiply by ((-1)^deg B lc(B))^(deg A - deg a) */
2803 0 : GEN c = gel(b,degB+2); /* lc(B) */
2804 0 : if (odd(degB)) c = Flx_neg(c, p);
2805 0 : if (!Flx_equal1(c)) {
2806 0 : c = Flx_powu_pre(c, dropa, p, pi);
2807 0 : if (!Flx_equal1(c)) Hp = Flx_mul_pre(Hp, c, p, pi);
2808 : }
2809 : }
2810 364906 : else if (dropb)
2811 : { /* multiply by lc(A)^(deg B - deg b) */
2812 0 : ulong c = uel(a, degA+2); /* lc(A) */
2813 0 : c = Fl_powu(c, dropb, p);
2814 0 : if (c != 1) Hp = Flx_Fl_mul_pre(Hp, c, p, pi);
2815 : }
2816 : }
2817 364906 : if (dp != 1) Hp = Flx_Fl_mul_pre(Hp, Fl_powu_pre(Fl_inv(dp,p), degA, p, pi), p, pi);
2818 364905 : return Hp;
2819 : }
2820 :
2821 : static GEN
2822 124963 : ZX_ZXY_resultant_slice(GEN A, GEN B, GEN dB, long degA, long degB, long dres,
2823 : GEN P, GEN *mod, long sX, long vY)
2824 : {
2825 124963 : pari_sp av = avma;
2826 124963 : long i, n = lg(P)-1;
2827 : GEN H, T, D;
2828 124963 : if (n == 1)
2829 : {
2830 40164 : ulong p = uel(P,1);
2831 40164 : ulong dp = dB ? umodiu(dB, p): 1;
2832 40164 : GEN a = ZX_to_Flx(A, p), b = ZXX_to_FlxX(B, p, vY);
2833 40164 : GEN Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2834 40163 : H = gc_upto(av, Flx_to_ZX(Hp));
2835 40164 : *mod = utoipos(p); return H;
2836 : }
2837 84799 : T = ZV_producttree(P);
2838 84799 : A = ZX_nv_mod_tree(A, P, T);
2839 84799 : B = ZXX_nv_mod_tree(B, P, T, vY);
2840 84799 : D = dB ? Z_ZV_mod_tree(dB, P, T): NULL;
2841 84799 : H = cgetg(n+1, t_VEC);
2842 364207 : for(i=1; i <= n; i++)
2843 : {
2844 279411 : ulong p = P[i];
2845 279411 : GEN a = gel(A,i), b = gel(B,i);
2846 279411 : ulong dp = D ? uel(D, i): 1;
2847 279411 : gel(H,i) = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2848 : }
2849 84796 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2850 84798 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2851 : }
2852 :
2853 : GEN
2854 124963 : ZX_ZXY_resultant_worker(GEN P, GEN A, GEN B, GEN dB, GEN v)
2855 : {
2856 124963 : GEN V = cgetg(3, t_VEC);
2857 124963 : if (isintzero(dB)) dB = NULL;
2858 124963 : gel(V,1) = ZX_ZXY_resultant_slice(A, B, dB, v[1], v[2], v[3], P, &gel(V,2), v[4], v[5]);
2859 124964 : return V;
2860 : }
2861 :
2862 : GEN
2863 79197 : ZX_ZXY_resultant(GEN A, GEN B)
2864 : {
2865 79197 : pari_sp av = avma;
2866 : forprime_t S;
2867 : ulong bound;
2868 79197 : long v = fetch_var_higher();
2869 79197 : long degA = degpol(A), degB, dres = degA * degpol(B);
2870 79197 : long vX = varn(B), vY = varn(A); /* assume vY has lower priority */
2871 79197 : long sX = evalvarn(vX);
2872 : GEN worker, H, dB;
2873 79197 : B = Q_remove_denom(B, &dB);
2874 79197 : if (!dB) B = leafcopy(B);
2875 79197 : A = leafcopy(A); setvarn(A,v);
2876 79197 : B = swap_vars(B, vY, v); degB = degpol(B);
2877 79196 : bound = ZX_ZXY_ResBound(A, B, dB);
2878 79197 : if (DEBUGLEVEL>4) err_printf("bound for resultant coeffs: 2^%ld\n",bound);
2879 158394 : worker = snm_closure(is_entry("_ZX_ZXY_resultant_worker"),
2880 79197 : mkvec4(A, B, dB? dB: gen_0,
2881 : mkvecsmall5(degA, degB, dres, sX, vY)));
2882 79197 : init_modular_big(&S);
2883 79197 : H = gen_crt("ZX_ZXY_resultant_all", worker, &S, dB, bound, 0, NULL,
2884 : nxV_chinese_center, FpX_center_i);
2885 79197 : setvarn(H, vX); (void)delete_var();
2886 79197 : return gc_GEN(av, H);
2887 : }
2888 :
2889 : static long
2890 40522 : ZX_ZXY_rnfequation_lambda(GEN A, GEN B0, long lambda)
2891 : {
2892 40522 : pari_sp av = avma;
2893 40522 : long degA = degpol(A), degB, dres = degA*degpol(B0);
2894 40522 : long v = fetch_var_higher();
2895 40522 : long vX = varn(B0), vY = varn(A); /* assume vY has lower priority */
2896 40522 : long sX = evalvarn(vX);
2897 : GEN dB, B, a, b, Hp;
2898 : forprime_t S;
2899 :
2900 40522 : B0 = Q_remove_denom(B0, &dB);
2901 40522 : if (!dB) B0 = leafcopy(B0);
2902 40522 : A = leafcopy(A);
2903 40522 : B = B0;
2904 40522 : setvarn(A,v);
2905 45334 : INIT:
2906 45334 : if (lambda) B = RgX_Rg_translate(B0, monomial(stoi(lambda), 1, vY));
2907 45334 : B = swap_vars(B, vY, v);
2908 : /* B0(lambda v + x, v) */
2909 45334 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
2910 :
2911 45334 : degB = degpol(B);
2912 45334 : init_modular_big(&S);
2913 : while (1)
2914 0 : {
2915 45334 : ulong p = u_forprime_next(&S);
2916 45334 : ulong dp = dB ? umodiu(dB, p): 1;
2917 45334 : if (!dp) continue;
2918 45334 : a = ZX_to_Flx(A, p);
2919 45334 : b = ZXX_to_FlxX(B, p, v);
2920 45334 : Hp = ZX_ZXY_resultant_prime(a, b, dp, p, degA, degB, dres, sX);
2921 45334 : if (degpol(Hp) != dres) continue;
2922 45334 : if (dp != 1) Hp = Flx_Fl_mul(Hp, Fl_powu(Fl_inv(dp,p), degA, p), p);
2923 45334 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); goto INIT; }
2924 40523 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
2925 40523 : (void)delete_var(); return gc_long(av,lambda);
2926 : }
2927 : }
2928 :
2929 : GEN
2930 60562 : ZX_ZXY_rnfequation(GEN A, GEN B, long *lambda)
2931 : {
2932 60562 : if (lambda)
2933 : {
2934 40522 : *lambda = ZX_ZXY_rnfequation_lambda(A, B, *lambda);
2935 40523 : if (*lambda) B = RgX_Rg_translate(B, monomial(stoi(*lambda), 1, varn(A)));
2936 : }
2937 60563 : return ZX_ZXY_resultant(A,B);
2938 : }
2939 :
2940 : static GEN
2941 10347 : ZX_composedsum_slice(GEN A, GEN B, GEN P, GEN *mod)
2942 : {
2943 10347 : pari_sp av = avma;
2944 10347 : long i, n = lg(P)-1;
2945 : GEN H, T;
2946 10347 : if (n == 1)
2947 : {
2948 9845 : ulong p = uel(P,1);
2949 9845 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
2950 9846 : GEN Hp = Flx_composedsum(a, b, p);
2951 9846 : H = gc_upto(av, Flx_to_ZX(Hp));
2952 9848 : *mod = utoipos(p); return H;
2953 : }
2954 502 : T = ZV_producttree(P);
2955 502 : A = ZX_nv_mod_tree(A, P, T);
2956 502 : B = ZX_nv_mod_tree(B, P, T);
2957 502 : H = cgetg(n+1, t_VEC);
2958 4526 : for(i=1; i <= n; i++)
2959 : {
2960 4024 : ulong p = P[i];
2961 4024 : GEN a = gel(A,i), b = gel(B,i);
2962 4024 : gel(H,i) = Flx_composedsum(a, b, p);
2963 : }
2964 502 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
2965 502 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
2966 : }
2967 :
2968 : GEN
2969 10348 : ZX_composedsum_worker(GEN P, GEN A, GEN B)
2970 : {
2971 10348 : GEN V = cgetg(3, t_VEC);
2972 10347 : gel(V,1) = ZX_composedsum_slice(A, B, P, &gel(V,2));
2973 10350 : return V;
2974 : }
2975 :
2976 : static GEN
2977 10086 : ZX_composedsum_i(GEN A, GEN B, GEN lead)
2978 : {
2979 10086 : pari_sp av = avma;
2980 : forprime_t S;
2981 : ulong bound;
2982 : GEN H, worker, mod;
2983 10086 : if (degpol(A) < degpol(B)) swap(A, B);
2984 10086 : if (!lead) lead = mulii(leading_coeff(A),leading_coeff(B));
2985 10086 : bound = ZX_ZXY_ResBound_1(A, B);
2986 10085 : worker = snm_closure(is_entry("_ZX_composedsum_worker"), mkvec2(A,B));
2987 10087 : init_modular_big(&S);
2988 10084 : H = gen_crt("ZX_composedsum", worker, &S, lead, bound, 0, &mod,
2989 : nxV_chinese_center, FpX_center);
2990 10087 : return gc_upto(av, H);
2991 : }
2992 :
2993 : static long
2994 9700 : ZX_compositum_lambda(GEN A, GEN B, GEN lead, long lambda)
2995 : {
2996 9700 : pari_sp av = avma;
2997 : forprime_t S;
2998 : ulong p;
2999 9700 : init_modular_big(&S);
3000 9700 : p = u_forprime_next(&S);
3001 : while (1)
3002 112 : {
3003 : GEN Hp, a;
3004 9812 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
3005 9812 : if (lead && dvdiu(lead,p)) { p = u_forprime_next(&S); continue; }
3006 9805 : a = ZX_to_Flx(ZX_rescale(A, stoi(-lambda)), p);
3007 9807 : Hp = Flx_composedsum(a, ZX_to_Flx(B, p), p);
3008 9799 : if (!Flx_is_squarefree(Hp, p)) { lambda = next_lambda(lambda); continue; }
3009 9692 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
3010 9692 : return gc_long(av, lambda);
3011 : }
3012 : }
3013 :
3014 : GEN
3015 9700 : ZX_compositum(GEN A, GEN B, long *lambda)
3016 : {
3017 9700 : GEN lead = mulii(leading_coeff(A),leading_coeff(B));
3018 9700 : if (lambda)
3019 : {
3020 9700 : *lambda = ZX_compositum_lambda(A, B, lead, *lambda);
3021 9692 : A = ZX_rescale(A, stoi(-*lambda));
3022 : }
3023 9701 : return ZX_composedsum_i(A, B, lead);
3024 : }
3025 :
3026 : GEN
3027 385 : ZX_composedsum(GEN A, GEN B)
3028 385 : { return ZX_composedsum_i(A, B, NULL); }
3029 :
3030 : static GEN
3031 359 : ZXQX_composedsum_slice(GEN A, GEN B, GEN C, GEN P, GEN *mod)
3032 : {
3033 359 : pari_sp av = avma;
3034 359 : long i, n = lg(P)-1, dC = degpol(C), v = varn(C);
3035 : GEN H, T;
3036 359 : if (n == 1)
3037 : {
3038 181 : ulong p = uel(P,1);
3039 181 : GEN a = ZXX_to_FlxX(A, p, v), b = ZXX_to_FlxX(B, p, v);
3040 181 : GEN c = ZX_to_Flx(C, p);
3041 181 : GEN Hp = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
3042 181 : H = gc_upto(av, Flm_to_ZM(Hp));
3043 181 : *mod = utoipos(p); return H;
3044 : }
3045 178 : T = ZV_producttree(P);
3046 178 : A = ZXX_nv_mod_tree(A, P, T, v);
3047 178 : B = ZXX_nv_mod_tree(B, P, T, v);
3048 178 : C = ZX_nv_mod_tree(C, P, T);
3049 178 : H = cgetg(n+1, t_VEC);
3050 660 : for(i=1; i <= n; i++)
3051 : {
3052 482 : ulong p = P[i];
3053 482 : GEN a = gel(A,i), b = gel(B,i), c = gel(C,i);
3054 482 : gel(H,i) = FlxX_to_Flm(FlxqX_composedsum(a, b, c, p), dC);
3055 : }
3056 178 : H = nmV_chinese_center_tree_seq(H, P, T, ZV_chinesetree(P, T));
3057 178 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
3058 : }
3059 :
3060 : GEN
3061 359 : ZXQX_composedsum_worker(GEN P, GEN A, GEN B, GEN C)
3062 : {
3063 359 : GEN V = cgetg(3, t_VEC);
3064 359 : gel(V,1) = ZXQX_composedsum_slice(A, B, C, P, &gel(V,2));
3065 359 : return V;
3066 : }
3067 :
3068 : static GEN
3069 315 : ZXQX_composedsum(GEN A, GEN B, GEN T, ulong bound)
3070 : {
3071 315 : pari_sp av = avma;
3072 : forprime_t S;
3073 : GEN H, worker, mod;
3074 315 : GEN lead = mulii(Q_content(leading_coeff(A)), Q_content(leading_coeff(B)));
3075 315 : worker = snm_closure(is_entry("_ZXQX_composedsum_worker")
3076 : , mkvec3(A,B,T));
3077 315 : init_modular_big(&S);
3078 315 : H = gen_crt("ZXQX_composedsum", worker, &S, lead, bound, 0, &mod,
3079 : nmV_chinese_center, FpM_center);
3080 315 : if (DEBUGLEVEL > 4)
3081 0 : err_printf("nfcompositum: a priori bound: %lu, a posteriori: %lu\n",
3082 : bound, expi(gsupnorm(H, DEFAULTPREC)));
3083 315 : return gc_GEN(av, RgM_to_RgXX(H, varn(A), varn(T)));
3084 : }
3085 :
3086 : static long
3087 315 : ZXQX_composedsum_bound(GEN nf, GEN A, GEN B)
3088 315 : { return ZXQX_resultant_bound_i(nf, A, B, &RgX_RgXY_ResBound_1); }
3089 :
3090 : GEN
3091 315 : nf_direct_compositum(GEN nf, GEN A, GEN B)
3092 : {
3093 315 : ulong bnd = ZXQX_composedsum_bound(nf, A, B);
3094 315 : return ZXQX_composedsum(A, B, nf_get_pol(nf), bnd);
3095 : }
3096 :
3097 : /************************************************************************
3098 : * *
3099 : * IRREDUCIBLE POLYNOMIAL / Fp *
3100 : * *
3101 : ************************************************************************/
3102 :
3103 : /* irreducible (unitary) polynomial of degree n over Fp */
3104 : GEN
3105 0 : ffinit_rand(GEN p,long n)
3106 : {
3107 0 : for(;;) {
3108 0 : pari_sp av = avma;
3109 0 : GEN pol = ZX_add(pol_xn(n, 0), random_FpX(n-1,0, p));
3110 0 : if (FpX_is_irred(pol, p)) return pol;
3111 0 : set_avma(av);
3112 : }
3113 : }
3114 :
3115 : /* return an extension of degree 2^l of F_2, assume l > 0
3116 : * Not stack clean. */
3117 : static GEN
3118 602 : ffinit_Artin_Schreier_2(long l)
3119 : {
3120 : GEN Q, T, S;
3121 : long i, v;
3122 :
3123 602 : if (l == 1) return mkvecsmall4(0,1,1,1); /*x^2 + x + 1*/
3124 553 : v = fetch_var_higher();
3125 553 : S = mkvecsmall5(0, 0, 0, 1, 1); /* y(y^2 + y) */
3126 553 : Q = mkpoln(3, pol1_Flx(0), pol1_Flx(0), S); /* x^2 + x + y(y^2+y) */
3127 553 : setvarn(Q, v);
3128 :
3129 : /* x^4+x+1, irred over F_2, minimal polynomial of a root of Q */
3130 553 : T = mkvecsmalln(6,evalvarn(v),1UL,1UL,0UL,0UL,1UL);
3131 : /* Q = x^2 + x + a(y) irred. over K = F2[y] / (T(y))
3132 : * ==> x^2 + x + a(y) b irred. over K for any root b of Q
3133 : * ==> x^2 + x + (b^2+b)b */
3134 3052 : for (i=2; i<l; i++) T = Flx_FlxY_resultant(T, Q, 2); /* minpoly(b) / F2*/
3135 553 : (void)delete_var(); T[1] = 0; return T;
3136 : }
3137 :
3138 : /* return an extension of degree p^l of F_p, assume l > 0
3139 : * Not stack clean. */
3140 : GEN
3141 959 : ffinit_Artin_Schreier(ulong p, long l)
3142 : {
3143 : long i, v;
3144 : GEN Q, R, S, T, xp;
3145 959 : if (p==2) return ffinit_Artin_Schreier_2(l);
3146 357 : xp = polxn_Flx(p,0); /* x^p */
3147 357 : T = Flx_sub(xp, mkvecsmall3(0,1,1),p); /* x^p - x - 1 */
3148 357 : if (l == 1) return T;
3149 :
3150 7 : v = evalvarn(fetch_var_higher());
3151 7 : xp[1] = v;
3152 7 : R = Flx_sub(polxn_Flx(2*p-1,0), polxn_Flx(p,0),p);
3153 7 : S = Flx_sub(xp, polx_Flx(0), p);
3154 7 : Q = FlxX_Flx_sub(Flx_to_FlxX(S, v), R, p); /* x^p - x - (y^(2p-1)-y^p) */
3155 14 : for (i = 2; i <= l; ++i) T = Flx_FlxY_resultant(T, Q, p);
3156 7 : (void)delete_var(); T[1] = 0; return T;
3157 : }
3158 :
3159 : static long
3160 149713 : flinit_check(ulong p, long n, long l)
3161 : {
3162 : ulong q;
3163 149713 : if (!uisprime(n)) return 0;
3164 102496 : q = p % n; if (!q) return 0;
3165 99885 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3166 : }
3167 :
3168 : static GEN
3169 31944 : flinit(ulong p, long l)
3170 : {
3171 31944 : ulong n = 1+l;
3172 96782 : while (!flinit_check(p,n,l)) n += l;
3173 31944 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3174 31944 : return ZX_to_Flx(polsubcyclo(n,l,0), p);
3175 : }
3176 :
3177 : static GEN
3178 28993 : ffinit_fact_Flx(ulong p, long n)
3179 : {
3180 28993 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3181 28993 : long i, l = lg(Fm);
3182 28993 : P = cgetg(l, t_VEC);
3183 61896 : for (i = 1; i < l; i++)
3184 32903 : gel(P,i) = p==uel(Fp,i) ? ffinit_Artin_Schreier(p, Fe[i])
3185 32903 : : flinit(p, uel(Fm,i));
3186 28993 : return FlxV_composedsum(P, p);
3187 : }
3188 :
3189 : static GEN
3190 52938 : init_Flxq_i(ulong p, long n, long sv)
3191 : {
3192 : GEN P;
3193 52938 : if (!odd(p) && p != 2) pari_err_PRIME("ffinit", utoi(p));
3194 52931 : if (n == 1) return polx_Flx(sv);
3195 52931 : if (flinit_check(p, n+1, n))
3196 : {
3197 23938 : P = const_vecsmall(n+2,1);
3198 23938 : P[1] = sv; return P;
3199 : }
3200 28993 : P = ffinit_fact_Flx(p,n);
3201 28993 : P[1] = sv; return P;
3202 : }
3203 :
3204 : GEN
3205 0 : init_Flxq(ulong p, long n, long v)
3206 : {
3207 0 : pari_sp av = avma;
3208 0 : return gc_upto(av, init_Flxq_i(p, n, v));
3209 : }
3210 :
3211 : /* check if polsubcyclo(n,l,0) is irreducible modulo p */
3212 : static long
3213 8207 : fpinit_check(GEN p, long n, long l)
3214 : {
3215 : ulong q;
3216 8207 : if (!uisprime(n)) return 0;
3217 4842 : q = umodiu(p,n); if (!q) return 0;
3218 4842 : return ugcd((n-1)/Fl_order(q, n-1, n), l) == 1;
3219 : }
3220 :
3221 : /* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l.
3222 : * Return an irreducible polynomial of degree l over F_p.
3223 : * Variant of Adleman and Lenstra "Finding irreducible polynomials over
3224 : * finite fields", ACM, 1986 (5) 350--355.
3225 : * Not stack clean */
3226 : static GEN
3227 1828 : fpinit(GEN p, long l)
3228 : {
3229 1828 : ulong n = 1+l;
3230 6168 : while (!fpinit_check(p,n,l)) n += l;
3231 1828 : if (DEBUGLEVEL>=4) err_printf("FFInit: using polsubcyclo(%ld, %ld)\n",n,l);
3232 1828 : return FpX_red(polsubcyclo(n,l,0),p);
3233 : }
3234 :
3235 : static GEN
3236 1637 : ffinit_fact(GEN p, long n)
3237 : {
3238 1637 : GEN P, F = factoru_pow(n), Fp = gel(F,1), Fe = gel(F,2), Fm = gel(F,3);
3239 1637 : long i, l = lg(Fm);
3240 1637 : P = cgetg(l, t_VEC);
3241 3465 : for (i = 1; i < l; ++i)
3242 3656 : gel(P,i) = absequaliu(p, Fp[i]) ?
3243 0 : Flx_to_ZX(ffinit_Artin_Schreier(Fp[i], Fe[i]))
3244 1828 : : fpinit(p, Fm[i]);
3245 1637 : return FpXV_composedsum(P, p);
3246 : }
3247 :
3248 : static GEN
3249 55243 : init_Fq_i(GEN p, long n, long v)
3250 : {
3251 : GEN P;
3252 55243 : if (n <= 0) pari_err_DOMAIN("ffinit", "degree", "<=", gen_0, stoi(n));
3253 55243 : if (typ(p) != t_INT) pari_err_TYPE("ffinit",p);
3254 55243 : if (cmpiu(p, 2) < 0) pari_err_PRIME("ffinit",p);
3255 55236 : if (v < 0) v = 0;
3256 55236 : if (n == 1) return pol_x(v);
3257 54984 : if (lgefint(p) == 3)
3258 52938 : return Flx_to_ZX(init_Flxq_i(p[2], n, evalvarn(v)));
3259 2046 : if (!mpodd(p)) pari_err_PRIME("ffinit", p);
3260 2039 : if (fpinit_check(p, n+1, n)) return polcyclo(n+1, v);
3261 1637 : P = ffinit_fact(p,n);
3262 1637 : setvarn(P, v); return P;
3263 : }
3264 : GEN
3265 54676 : init_Fq(GEN p, long n, long v)
3266 : {
3267 54676 : pari_sp av = avma;
3268 54676 : return gc_upto(av, init_Fq_i(p, n, v));
3269 : }
3270 : GEN
3271 567 : ffinit(GEN p, long n, long v)
3272 : {
3273 567 : pari_sp av = avma;
3274 567 : return gc_upto(av, FpX_to_mod(init_Fq_i(p, n, v), p));
3275 : }
3276 :
3277 : GEN
3278 3178 : ffnbirred(GEN p, long n)
3279 : {
3280 3178 : pari_sp av = avma;
3281 3178 : GEN s = powiu(p,n), F = factoru(n), D = divisorsu_moebius(gel(F, 1));
3282 3178 : long j, l = lg(D);
3283 6797 : for (j = 2; j < l; j++) /* skip d = 1 */
3284 : {
3285 3619 : long md = D[j]; /* mu(d) * d, d squarefree */
3286 3619 : GEN pd = powiu(p, n / labs(md)); /* p^{n/d} */
3287 3619 : s = md > 0? addii(s, pd): subii(s,pd);
3288 : }
3289 3178 : return gc_INT(av, diviuexact(s, n));
3290 : }
3291 :
3292 : GEN
3293 616 : ffsumnbirred(GEN p, long n)
3294 : {
3295 616 : pari_sp av = avma, av2;
3296 616 : GEN q, t = p, v = vecfactoru_i(1, n);
3297 : long i;
3298 616 : q = cgetg(n+1,t_VEC); gel(q,1) = p;
3299 1764 : for (i=2; i<=n; i++) gel(q,i) = mulii(gel(q,i-1), p);
3300 616 : av2 = avma;
3301 1764 : for (i=2; i<=n; i++)
3302 : {
3303 1148 : GEN s = gel(q,i), F = gel(v,i), D = divisorsu_moebius(gel(F,1));
3304 1148 : long j, l = lg(D);
3305 2534 : for (j = 2; j < l; j++) /* skip 1 */
3306 : {
3307 1386 : long md = D[j];
3308 1386 : GEN pd = gel(q, i / labs(md)); /* p^{i/d} */
3309 1386 : s = md > 0? addii(s, pd): subii(s, pd);
3310 : }
3311 1148 : t = gc_INT(av2, addii(t, diviuexact(s, i)));
3312 : }
3313 616 : return gc_INT(av, t);
3314 : }
3315 :
3316 : GEN
3317 140 : ffnbirred0(GEN p, long n, long flag)
3318 : {
3319 140 : if (typ(p) != t_INT) pari_err_TYPE("ffnbirred", p);
3320 140 : if (n <= 0) pari_err_DOMAIN("ffnbirred", "degree", "<=", gen_0, stoi(n));
3321 140 : switch(flag)
3322 : {
3323 70 : case 0: return ffnbirred(p, n);
3324 70 : case 1: return ffsumnbirred(p, n);
3325 : }
3326 0 : pari_err_FLAG("ffnbirred");
3327 : return NULL; /* LCOV_EXCL_LINE */
3328 : }
3329 :
3330 : static void
3331 2261 : checkmap(GEN m, const char *s)
3332 : {
3333 2261 : if (typ(m)!=t_VEC || lg(m)!=3 || typ(gel(m,1))!=t_FFELT)
3334 0 : pari_err_TYPE(s,m);
3335 2261 : }
3336 :
3337 : GEN
3338 189 : ffembed(GEN a, GEN b)
3339 : {
3340 189 : pari_sp av = avma;
3341 189 : GEN p, Ta, Tb, g, r = NULL;
3342 189 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffembed",a);
3343 189 : if (typ(b)!=t_FFELT) pari_err_TYPE("ffembed",b);
3344 189 : p = FF_p_i(a); g = FF_gen(a);
3345 189 : if (!equalii(p, FF_p_i(b))) pari_err_MODULUS("ffembed",a,b);
3346 189 : Ta = FF_mod(a);
3347 189 : Tb = FF_mod(b);
3348 189 : if (degpol(Tb)%degpol(Ta)!=0)
3349 7 : pari_err_DOMAIN("ffembed",GENtostr_raw(a),"is not a subfield of",b,a);
3350 182 : r = gel(FFX_roots(Ta, b), 1);
3351 182 : return gc_GEN(av, mkvec2(g,r));
3352 : }
3353 :
3354 : GEN
3355 91 : ffextend(GEN a, GEN P, long v)
3356 : {
3357 91 : pari_sp av = avma;
3358 : long n;
3359 : GEN p, T, R, g, m;
3360 91 : if (typ(a)!=t_FFELT) pari_err_TYPE("ffextend",a);
3361 91 : T = a; p = FF_p_i(a);
3362 91 : if (typ(P)!=t_POL || !RgX_is_FpXQX(P,&T,&p)) pari_err_TYPE("ffextend", P);
3363 49 : if (!FF_samefield(a, T)) pari_err_MODULUS("ffextend",a,T);
3364 49 : if (v < 0) v = varn(P);
3365 49 : n = FF_f(T) * degpol(P); R = ffinit(p, n, v); g = ffgen(R, v);
3366 49 : m = ffembed(a, g);
3367 49 : R = FFX_roots(ffmap(m, P),g);
3368 49 : return gc_GEN(av, mkvec2(gel(R,1), m));
3369 : }
3370 :
3371 : GEN
3372 42 : fffrobenius(GEN a, long n)
3373 : {
3374 42 : if (typ(a)!=t_FFELT) pari_err_TYPE("fffrobenius",a);
3375 42 : retmkvec2(FF_gen(a), FF_Frobenius(a, n));
3376 : }
3377 :
3378 : GEN
3379 133 : ffinvmap(GEN m)
3380 : {
3381 133 : pari_sp av = avma;
3382 : long i, l;
3383 133 : GEN T, F, a, g, r, f = NULL;
3384 133 : checkmap(m, "ffinvmap");
3385 133 : a = gel(m,1); r = gel(m,2);
3386 133 : if (typ(r) != t_FFELT)
3387 7 : pari_err_TYPE("ffinvmap", m);
3388 126 : g = FF_gen(a);
3389 126 : T = FF_mod(r);
3390 126 : F = gel(FFX_factor(T, a), 1);
3391 126 : l = lg(F);
3392 490 : for(i=1; i<l; i++)
3393 : {
3394 490 : GEN s = FFX_rem(FF_to_FpXQ_i(r), gel(F, i), a);
3395 490 : if (degpol(s)==0 && gequal(constant_coeff(s),g)) { f = gel(F, i); break; }
3396 : }
3397 126 : if (f==NULL) pari_err_TYPE("ffinvmap", m);
3398 126 : if (degpol(f)==1) f = FF_neg_i(gel(f,2));
3399 126 : return gc_GEN(av, mkvec2(FF_gen(r),f));
3400 : }
3401 :
3402 : static GEN
3403 1260 : ffpartmapimage(const char *s, GEN r)
3404 : {
3405 1260 : GEN a = NULL, p = NULL;
3406 1260 : if (typ(r)==t_POL && degpol(r) >= 1
3407 1260 : && RgX_is_FpXQX(r,&a,&p) && a && typ(a)==t_FFELT) return a;
3408 0 : pari_err_TYPE(s, r);
3409 : return NULL; /* LCOV_EXCL_LINE */
3410 : }
3411 :
3412 : static GEN
3413 2709 : ffeltmap_i(GEN m, GEN x)
3414 : {
3415 2709 : GEN r = gel(m,2);
3416 2709 : if (!FF_samefield(x, gel(m,1)))
3417 84 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3418 2625 : if (typ(r)==t_FFELT)
3419 1659 : return FF_map(r, x);
3420 : else
3421 966 : return FFX_preimage(x, r, ffpartmapimage("ffmap", r));
3422 : }
3423 :
3424 : static GEN
3425 4459 : ffmap_i(GEN m, GEN x)
3426 : {
3427 : GEN y;
3428 4459 : long i, lx, tx = typ(x);
3429 4459 : switch(tx)
3430 : {
3431 2541 : case t_FFELT:
3432 2541 : return ffeltmap_i(m, x);
3433 1267 : case t_POL: case t_RFRAC: case t_SER:
3434 : case t_VEC: case t_COL: case t_MAT:
3435 1267 : y = cgetg_copy(x, &lx);
3436 1988 : for (i = 1; i < lontyp[tx]; i++) y[i] = x[i];
3437 4564 : for (; i < lx; i++)
3438 : {
3439 3339 : GEN yi = ffmap_i(m, gel(x,i));
3440 3297 : if (!yi) return NULL;
3441 3297 : gel(y,i) = yi;
3442 : }
3443 1225 : return y;
3444 : }
3445 651 : return gcopy(x);
3446 : }
3447 :
3448 : GEN
3449 1036 : ffmap(GEN m, GEN x)
3450 : {
3451 1036 : pari_sp ltop = avma;
3452 : GEN y;
3453 1036 : checkmap(m, "ffmap");
3454 1036 : y = ffmap_i(m, x);
3455 1036 : if (y) return y;
3456 42 : retgc_const(ltop, cgetg(1, t_VEC));
3457 : }
3458 :
3459 : static GEN
3460 252 : ffeltmaprel_i(GEN m, GEN x)
3461 : {
3462 252 : GEN g = gel(m,1), r = gel(m,2);
3463 252 : if (!FF_samefield(x, g))
3464 0 : pari_err_DOMAIN("ffmap","m","domain does not contain", x, r);
3465 252 : if (typ(r)==t_FFELT)
3466 84 : retmkpolmod(FF_map(r, x), pol_x(FF_var(g)));
3467 : else
3468 168 : retmkpolmod(FFX_preimagerel(x, r, ffpartmapimage("ffmap", r)), gcopy(r));
3469 : }
3470 :
3471 : static GEN
3472 252 : ffmaprel_i(GEN m, GEN x)
3473 : {
3474 252 : switch(typ(x))
3475 : {
3476 252 : case t_FFELT:
3477 252 : return ffeltmaprel_i(m, x);
3478 0 : case t_POL: pari_APPLY_pol_normalized(ffmaprel_i(m, gel(x,i)));
3479 0 : case t_SER: pari_APPLY_ser_normalized(ffmaprel_i(m, gel(x,i)));
3480 0 : case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
3481 0 : pari_APPLY_same(ffmaprel_i(m, gel(x,i)));
3482 : }
3483 0 : return gcopy(x);
3484 : }
3485 : GEN
3486 252 : ffmaprel(GEN m, GEN x) { checkmap(m, "ffmaprel"); return ffmaprel_i(m, x); }
3487 :
3488 : static void
3489 84 : err_compo(GEN m, GEN n)
3490 84 : { pari_err_DOMAIN("ffcompomap","m","domain does not contain codomain of",n,m); }
3491 :
3492 : GEN
3493 420 : ffcompomap(GEN m, GEN n)
3494 : {
3495 420 : pari_sp av = avma;
3496 420 : GEN g = gel(n,1), r, m2, n2;
3497 420 : checkmap(m, "ffcompomap");
3498 420 : checkmap(n, "ffcompomap");
3499 420 : m2 = gel(m,2); n2 = gel(n,2);
3500 420 : switch((typ(m2)==t_POL)|((typ(n2)==t_POL)<<1))
3501 : {
3502 84 : case 0:
3503 84 : if (!FF_samefield(gel(m,1),n2)) err_compo(m,n);
3504 42 : r = FF_map(gel(m,2), n2);
3505 42 : break;
3506 84 : case 2:
3507 84 : r = ffmap_i(m, n2);
3508 42 : if (lg(r) == 1) err_compo(m,n);
3509 42 : break;
3510 168 : case 1:
3511 168 : r = ffeltmap_i(m, n2);
3512 126 : if (!r)
3513 : {
3514 : GEN a, A, R, M;
3515 : long dm, dn;
3516 42 : a = ffpartmapimage("ffcompomap",m2);
3517 42 : A = FF_to_FpXQ_i(FF_neg(n2));
3518 42 : setvarn(A, 1);
3519 42 : R = deg1pol(gen_1, A, 0);
3520 42 : setvarn(R, 0);
3521 42 : M = gcopy(m2);
3522 42 : setvarn(M, 1);
3523 42 : r = polresultant0(R, M, 1, 0);
3524 42 : dm = FF_f(gel(m,1)); dn = FF_f(gel(n,1));
3525 42 : if (dm % dn || !FFX_ispower(r, dm/dn, a, &r)) err_compo(m,n);
3526 42 : setvarn(r, varn(FF_mod(g)));
3527 : }
3528 126 : break;
3529 84 : case 3:
3530 : {
3531 : GEN M, R, T, p, a;
3532 84 : a = ffpartmapimage("ffcompomap",n2);
3533 84 : if (!FF_samefield(a, gel(m,1))) err_compo(m,n);
3534 42 : p = FF_p_i(gel(n,1));
3535 42 : T = FF_mod(gel(n,1));
3536 42 : setvarn(T, 1);
3537 42 : R = RgX_to_FpXQX(n2,T,p);
3538 42 : setvarn(R, 0);
3539 42 : M = gcopy(m2);
3540 42 : setvarn(M, 1);
3541 42 : r = polresultant0(R, M, 1, 0);
3542 42 : setvarn(r, varn(n2));
3543 : }
3544 : }
3545 252 : return gc_GEN(av, mkvec2(g,r));
3546 : }
|
