1 /* Copyright (C) 2012 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 #include "pari.h"
16 #include "paripriv.h"
17
18 /* Not so fast arithmetic with polynomials over FpX */
19
20 /*******************************************************************/
21 /* */
22 /* FpXX */
23 /* */
24 /*******************************************************************/
25 /*Polynomials whose coefficients are either polynomials or integers*/
26
27 static ulong
to_FlxqX(GEN P,GEN Q,GEN T,GEN p,GEN * pt_P,GEN * pt_Q,GEN * pt_T)28 to_FlxqX(GEN P, GEN Q, GEN T, GEN p, GEN *pt_P, GEN *pt_Q, GEN *pt_T)
29 {
30 ulong pp = uel(p,2);
31 long v = get_FpX_var(T);
32 *pt_P = ZXX_to_FlxX(P, pp, v);
33 if (pt_Q) *pt_Q = ZXX_to_FlxX(Q, pp, v);
34 *pt_T = ZXT_to_FlxT(T, pp);
35 return pp;
36 }
37
38 static GEN
ZXX_copy(GEN a)39 ZXX_copy(GEN a) { return gcopy(a); }
40
41 GEN
FpXX_red(GEN z,GEN p)42 FpXX_red(GEN z, GEN p)
43 {
44 GEN res;
45 long i, l = lg(z);
46 res = cgetg(l,t_POL); res[1] = z[1];
47 for (i=2; i<l; i++)
48 {
49 GEN zi = gel(z,i), c;
50 if (typ(zi)==t_INT)
51 c = modii(zi,p);
52 else
53 {
54 pari_sp av = avma;
55 c = FpX_red(zi,p);
56 switch(lg(c)) {
57 case 2: set_avma(av); c = gen_0; break;
58 case 3: c = gerepilecopy(av, gel(c,2)); break;
59 }
60 }
61 gel(res,i) = c;
62 }
63 return FpXX_renormalize(res,lg(res));
64 }
65 GEN
FpXX_add(GEN x,GEN y,GEN p)66 FpXX_add(GEN x, GEN y, GEN p)
67 {
68 long i,lz;
69 GEN z;
70 long lx=lg(x);
71 long ly=lg(y);
72 if (ly>lx) swapspec(x,y, lx,ly);
73 lz = lx; z = cgetg(lz, t_POL); z[1]=x[1];
74 for (i=2; i<ly; i++) gel(z,i) = Fq_add(gel(x,i), gel(y,i), NULL, p);
75 for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
76 return FpXX_renormalize(z, lz);
77 }
78 GEN
FpXX_sub(GEN x,GEN y,GEN p)79 FpXX_sub(GEN x, GEN y, GEN p)
80 {
81 long i,lz;
82 GEN z;
83 long lx=lg(x);
84 long ly=lg(y);
85 if (ly <= lx)
86 {
87 lz = lx; z = cgetg(lz, t_POL); z[1]=x[1];
88 for (i=2; i<ly; i++) gel(z,i) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
89 for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
90 }
91 else
92 {
93 lz = ly; z = cgetg(lz, t_POL); z[1]=x[1];
94 for (i=2; i<lx; i++) gel(z,i) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
95 for ( ; i<ly; i++) gel(z,i) = Fq_neg(gel(y,i), NULL, p);
96 }
97 return FpXX_renormalize(z, lz);
98 }
99
100 static GEN
FpXX_subspec(GEN x,GEN y,GEN p,long nx,long ny)101 FpXX_subspec(GEN x, GEN y, GEN p, long nx, long ny)
102 {
103 long i,lz;
104 GEN z;
105 if (ny <= nx)
106 {
107 lz = nx+2; z = cgetg(lz, t_POL);
108 for (i=0; i<ny; i++) gel(z,i+2) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
109 for ( ; i<nx; i++) gel(z,i+2) = gcopy(gel(x,i));
110 }
111 else
112 {
113 lz = ny+2; z = cgetg(lz, t_POL);
114 for (i=0; i<nx; i++) gel(z,i+2) = Fq_sub(gel(x,i), gel(y,i), NULL, p);
115 for ( ; i<ny; i++) gel(z,i+2) = Fq_neg(gel(y,i), NULL, p);
116 }
117 z[1] = 0; return FpXX_renormalize(z, lz);
118 }
119
120 GEN
FpXX_neg(GEN x,GEN p)121 FpXX_neg(GEN x, GEN p)
122 {
123 long i, lx = lg(x);
124 GEN y = cgetg(lx,t_POL);
125 y[1] = x[1];
126 for(i=2; i<lx; i++) gel(y,i) = Fq_neg(gel(x,i), NULL, p);
127 return FpXX_renormalize(y, lx);
128 }
129
130 GEN
FpXX_Fp_mul(GEN P,GEN u,GEN p)131 FpXX_Fp_mul(GEN P, GEN u, GEN p)
132 {
133 long i, lP;
134 GEN res = cgetg_copy(P, &lP); res[1] = P[1];
135 for(i=2; i<lP; i++)
136 {
137 GEN x = gel(P,i);
138 gel(res,i) = typ(x)==t_INT? Fp_mul(x,u,p): FpX_Fp_mul(x,u,p);
139 }
140 return FpXX_renormalize(res,lP);
141 }
142
143 GEN
FpXX_mulu(GEN P,ulong u,GEN p)144 FpXX_mulu(GEN P, ulong u, GEN p)
145 {
146 long i, lP;
147 GEN res = cgetg_copy(P, &lP); res[1] = P[1];
148 for(i=2; i<lP; i++)
149 {
150 GEN x = gel(P,i);
151 gel(res,i) = typ(x)==t_INT? Fp_mulu(x,u,p): FpX_mulu(x,u,p);
152 }
153 return FpXX_renormalize(res,lP);
154 }
155
156 GEN
FpXX_halve(GEN P,GEN p)157 FpXX_halve(GEN P, GEN p)
158 {
159 long i, lP;
160 GEN res = cgetg_copy(P, &lP); res[1] = P[1];
161 for(i=2; i<lP; i++)
162 {
163 GEN x = gel(P,i);
164 gel(res,i) = typ(x)==t_INT? Fp_halve(x,p): FpX_halve(x,p);
165 }
166 return FpXX_renormalize(res,lP);
167 }
168
169 GEN
FpXX_deriv(GEN P,GEN p)170 FpXX_deriv(GEN P, GEN p)
171 {
172 long i, l = lg(P)-1;
173 GEN res;
174
175 if (l < 3) return pol_0(varn(P));
176 res = cgetg(l, t_POL);
177 res[1] = P[1];
178 for (i=2; i<l ; i++)
179 {
180 GEN x = gel(P,i+1);
181 gel(res,i) = typ(x)==t_INT? Fp_mulu(x,i-1,p): FpX_mulu(x,i-1,p);
182 }
183 return FpXX_renormalize(res, l);
184 }
185
186 GEN
FpXX_integ(GEN P,GEN p)187 FpXX_integ(GEN P, GEN p)
188 {
189 long i, l = lg(P);
190 GEN res;
191
192 if (l == 2) return pol_0(varn(P));
193 res = cgetg(l+1, t_POL);
194 res[1] = P[1];
195 gel(res,2) = gen_0;
196 for (i=3; i<=l ; i++)
197 {
198 GEN x = gel(P,i-1);
199 if (signe(x))
200 {
201 GEN i1 = Fp_inv(utoi(i-2), p);
202 gel(res,i) = typ(x)==t_INT? Fp_mul(x,i1,p): FpX_Fp_mul(x,i1,p);
203 } else
204 gel(res,i) = gen_0;
205 }
206 return FpXX_renormalize(res, l+1);
207 }
208
209 /*******************************************************************/
210 /* */
211 /* (Fp[X]/(Q))[Y] */
212 /* */
213 /*******************************************************************/
214
215 static GEN
get_FpXQX_red(GEN T,GEN * B)216 get_FpXQX_red(GEN T, GEN *B)
217 {
218 if (typ(T)!=t_VEC) { *B=NULL; return T; }
219 *B = gel(T,1); return gel(T,2);
220 }
221
222 GEN
random_FpXQX(long d1,long v,GEN T,GEN p)223 random_FpXQX(long d1, long v, GEN T, GEN p)
224 {
225 long dT = get_FpX_degree(T), vT = get_FpX_var(T);
226 long i, d = d1+2;
227 GEN y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);
228 for (i=2; i<d; i++) gel(y,i) = random_FpX(dT, vT, p);
229 return FpXQX_renormalize(y,d);
230 }
231
232 /*Not stack clean*/
233 GEN
Kronecker_to_FpXQX(GEN Z,GEN T,GEN p)234 Kronecker_to_FpXQX(GEN Z, GEN T, GEN p)
235 {
236 long i,j,lx,l, N = (get_FpX_degree(T)<<1) + 1;
237 GEN x, t = cgetg(N,t_POL), z = FpX_red(Z, p);
238 t[1] = evalvarn(get_FpX_var(T));
239 l = lg(z); lx = (l-2) / (N-2);
240 x = cgetg(lx+3,t_POL);
241 x[1] = z[1];
242 for (i=2; i<lx+2; i++)
243 {
244 for (j=2; j<N; j++) gel(t,j) = gel(z,j);
245 z += (N-2);
246 gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
247 }
248 N = (l-2) % (N-2) + 2;
249 for (j=2; j<N; j++) gel(t,j) = gel(z,j);
250 gel(x,i) = FpX_rem(FpX_renormalize(t,N), T,p);
251 return FpXQX_renormalize(x, i+1);
252 }
253
254 GEN
FpXQX_red(GEN z,GEN T,GEN p)255 FpXQX_red(GEN z, GEN T, GEN p)
256 {
257 long i, l = lg(z);
258 GEN res = cgetg(l,t_POL); res[1] = z[1];
259 for(i=2;i<l;i++)
260 if (typ(gel(z,i)) == t_INT)
261 gel(res,i) = modii(gel(z,i),p);
262 else
263 gel(res,i) = FpXQ_red(gel(z,i),T,p);
264 return FpXQX_renormalize(res,l);
265 }
266
267 static GEN
to_intmod(GEN x,GEN p)268 to_intmod(GEN x, GEN p) { retmkintmod(modii(x, p), p); }
269
270 GEN
FpXQX_to_mod(GEN z,GEN T,GEN p)271 FpXQX_to_mod(GEN z, GEN T, GEN p)
272 {
273 long i,l = lg(z);
274 GEN x = cgetg(l, t_POL);
275 x[1] = z[1];
276 if (l == 2) return x;
277 p = icopy(p);
278 T = FpX_to_mod_raw(T, p);
279 for (i=2; i<l; i++)
280 {
281 GEN zi = gel(z,i);
282 gel(x,i) = typ(zi) == t_POL? mkpolmod(FpX_to_mod_raw(zi, p), T)
283 : to_intmod(zi, p);
284 }
285 return normalizepol_lg(x,l);
286 }
287
288 static GEN
FpXQX_to_mod_raw(GEN z,GEN T,GEN p)289 FpXQX_to_mod_raw(GEN z, GEN T, GEN p)
290 {
291 long i,l = lg(z);
292 GEN x = cgetg(l, t_POL);
293 x[1] = z[1];
294 if (l == 2) return x;
295 for (i=2; i<l; i++)
296 {
297 GEN zi = gel(z,i);
298 gel(x,i) = typ(zi) == t_POL? mkpolmod(FpX_to_mod_raw(zi, p), T)
299 : to_intmod(zi, p);
300 }
301 return normalizepol_lg(x,l);
302 }
303
304 INLINE GEN
FqX_to_mod_raw(GEN f,GEN T,GEN p)305 FqX_to_mod_raw(GEN f, GEN T, GEN p)
306 { return T?FpXQX_to_mod_raw(f, T, p): FpX_to_mod_raw(f, p); }
307
308 static GEN
FqXC_to_mod_raw(GEN x,GEN T,GEN p)309 FqXC_to_mod_raw(GEN x, GEN T, GEN p)
310 { pari_APPLY_type(t_COL, FqX_to_mod_raw(gel(x,i), T, p)) }
311
312 GEN
FqXC_to_mod(GEN z,GEN T,GEN p)313 FqXC_to_mod(GEN z, GEN T, GEN p)
314 {
315 GEN x;
316 long i,l = lg(z);
317 if (!T) return FpXC_to_mod(z, p);
318 x = cgetg(l, t_COL);
319 if (l == 1) return x;
320 p = icopy(p);
321 T = FpX_to_mod_raw(T, p);
322 for (i=1; i<l; i++)
323 gel(x,i) = FqX_to_mod_raw(gel(z, i), T, p);
324 return x;
325 }
326
327 GEN
FqXM_to_mod(GEN z,GEN T,GEN p)328 FqXM_to_mod(GEN z, GEN T, GEN p)
329 {
330 GEN x;
331 long i,l = lg(z);
332 if (!T) return FpXM_to_mod(z, p);
333 x = cgetg(l, t_MAT);
334 if (l == 1) return x;
335 p = icopy(p);
336 T = FpX_to_mod_raw(T, p);
337 for (i=1; i<l; i++)
338 gel(x,i) = FqXC_to_mod_raw(gel(z, i), T, p);
339 return x;
340 }
341
342 static int
ZXX_is_ZX_spec(GEN a,long na)343 ZXX_is_ZX_spec(GEN a,long na)
344 {
345 long i;
346 for(i=0;i<na;i++)
347 if(typ(gel(a,i))!=t_INT) return 0;
348 return 1;
349 }
350
351 static int
ZXX_is_ZX(GEN a)352 ZXX_is_ZX(GEN a) { return ZXX_is_ZX_spec(a+2,lgpol(a)); }
353
354 static GEN
FpXX_FpX_mulspec(GEN P,GEN U,GEN p,long v,long lU)355 FpXX_FpX_mulspec(GEN P, GEN U, GEN p, long v, long lU)
356 {
357 long i, lP =lg(P);
358 GEN res;
359 res = cgetg(lP, t_POL); res[1] = P[1];
360 for(i=2; i<lP; i++)
361 {
362 GEN Pi = gel(P,i);
363 gel(res,i) = typ(Pi)==t_INT? FpX_Fp_mulspec(U, Pi, p, lU):
364 FpX_mulspec(U, Pi+2, p, lU, lgpol(Pi));
365 setvarn(gel(res,i),v);
366 }
367 return FpXQX_renormalize(res,lP);
368 }
369
370 GEN
FpXX_FpX_mul(GEN P,GEN U,GEN p)371 FpXX_FpX_mul(GEN P, GEN U, GEN p)
372 { return FpXX_FpX_mulspec(P,U+2,p,varn(U),lgpol(U)); }
373
374 static GEN
FpXY_FpY_mulspec(GEN x,GEN y,GEN T,GEN p,long lx,long ly)375 FpXY_FpY_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
376 {
377 pari_sp av = avma;
378 long v = fetch_var();
379 GEN z = RgXY_swapspec(x,get_FpX_degree(T)-1,v,lx);
380 z = FpXX_FpX_mulspec(z,y,p,v,ly);
381 z = RgXY_swapspec(z+2,lx+ly+3,get_FpX_var(T),lgpol(z));
382 (void)delete_var(); return gerepilecopy(av,z);
383 }
384
385 static GEN
FpXQX_mulspec(GEN x,GEN y,GEN T,GEN p,long lx,long ly)386 FpXQX_mulspec(GEN x, GEN y, GEN T, GEN p, long lx, long ly)
387 {
388 pari_sp av = avma;
389 GEN z, kx, ky;
390 long n;
391 if (ZXX_is_ZX_spec(y,ly))
392 {
393 if (ZXX_is_ZX_spec(x,lx))
394 return FpX_mulspec(x,y,p,lx,ly);
395 else
396 return FpXY_FpY_mulspec(x,y,T,p,lx,ly);
397 } else if (ZXX_is_ZX_spec(x,lx))
398 return FpXY_FpY_mulspec(y,x,T,p,ly,lx);
399 n = get_FpX_degree(T);
400 kx = ZXX_to_Kronecker_spec(x, lx, n);
401 ky = ZXX_to_Kronecker_spec(y, ly, n);
402 z = Kronecker_to_FpXQX(ZX_mul(ky,kx), T, p);
403 return gerepileupto(av, z);
404 }
405
406 GEN
FpXQX_mul(GEN x,GEN y,GEN T,GEN p)407 FpXQX_mul(GEN x, GEN y, GEN T, GEN p)
408 {
409 GEN z = FpXQX_mulspec(x+2,y+2,T,p,lgpol(x),lgpol(y));
410 setvarn(z,varn(x)); return z;
411 }
412
413 GEN
FpXQX_sqr(GEN x,GEN T,GEN p)414 FpXQX_sqr(GEN x, GEN T, GEN p)
415 {
416 pari_sp av = avma;
417 GEN z, kx;
418 if (ZXX_is_ZX(x)) return ZX_sqr(x);
419 kx= ZXX_to_Kronecker(x, get_FpX_degree(T));
420 z = Kronecker_to_FpXQX(ZX_sqr(kx), T, p);
421 return gerepileupto(av, z);
422 }
423
424 GEN
FpXQX_FpXQ_mul(GEN P,GEN U,GEN T,GEN p)425 FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p)
426 {
427 long i, lP;
428 GEN res;
429 res = cgetg_copy(P, &lP); res[1] = P[1];
430 for(i=2; i<lP; i++)
431 gel(res,i) = typ(gel(P,i))==t_INT? FpX_Fp_mul(U, gel(P,i), p):
432 FpXQ_mul(U, gel(P,i), T,p);
433 return FpXQX_renormalize(res,lP);
434 }
435
436 /* x and y in Z[Y][X]. Assume T irreducible mod p */
437 static GEN
FpXQX_divrem_basecase(GEN x,GEN y,GEN T,GEN p,GEN * pr)438 FpXQX_divrem_basecase(GEN x, GEN y, GEN T, GEN p, GEN *pr)
439 {
440 long vx, dx, dy, dy1, dz, i, j, sx, lr;
441 pari_sp av0, av, tetpil;
442 GEN z,p1,rem,lead;
443
444 if (!signe(y)) pari_err_INV("FpX_divrem",y);
445 vx=varn(x); dy=degpol(y); dx=degpol(x);
446 if (dx < dy)
447 {
448 if (pr)
449 {
450 av0 = avma; x = FpXQX_red(x, T, p);
451 if (pr == ONLY_DIVIDES) { set_avma(av0); return signe(x)? NULL: pol_0(vx); }
452 if (pr == ONLY_REM) return x;
453 *pr = x;
454 }
455 return pol_0(vx);
456 }
457 lead = leading_coeff(y);
458 if (!dy) /* y is constant */
459 {
460 if (pr && pr != ONLY_DIVIDES)
461 {
462 if (pr == ONLY_REM) return pol_0(vx);
463 *pr = pol_0(vx);
464 }
465 if (gequal1(lead)) return FpXQX_red(x,T,p);
466 av0 = avma; x = FqX_Fq_mul(x, Fq_inv(lead, T,p), T,p);
467 return gerepileupto(av0,x);
468 }
469 av0 = avma; dz = dx-dy;
470 lead = gequal1(lead)? NULL: gclone(Fq_inv(lead,T,p));
471 set_avma(av0);
472 z = cgetg(dz+3,t_POL); z[1] = x[1];
473 x += 2; y += 2; z += 2;
474 for (dy1=dy-1; dy1>=0 && !signe(gel(y, dy1)); dy1--);
475
476 p1 = gel(x,dx); av = avma;
477 gel(z,dz) = lead? gerepileupto(av, Fq_mul(p1,lead, T, p)): gcopy(p1);
478 for (i=dx-1; i>=dy; i--)
479 {
480 av=avma; p1=gel(x,i);
481 for (j=i-dy1; j<=i && j<=dz; j++)
482 p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
483 if (lead) p1 = Fq_mul(p1, lead, NULL,p);
484 tetpil=avma; gel(z,i-dy) = gerepile(av,tetpil,Fq_red(p1,T,p));
485 }
486 if (!pr) { guncloneNULL(lead); return z-2; }
487
488 rem = (GEN)avma; av = (pari_sp)new_chunk(dx+3);
489 for (sx=0; ; i--)
490 {
491 p1 = gel(x,i);
492 for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
493 p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j),NULL,p),NULL,p);
494 tetpil=avma; p1 = Fq_red(p1, T, p); if (signe(p1)) { sx = 1; break; }
495 if (!i) break;
496 set_avma(av);
497 }
498 if (pr == ONLY_DIVIDES)
499 {
500 guncloneNULL(lead);
501 if (sx) return gc_NULL(av0);
502 return gc_const((pari_sp)rem, z-2);
503 }
504 lr=i+3; rem -= lr;
505 rem[0] = evaltyp(t_POL) | evallg(lr);
506 rem[1] = z[-1];
507 p1 = gerepile((pari_sp)rem,tetpil,p1);
508 rem += 2; gel(rem,i) = p1;
509 for (i--; i>=0; i--)
510 {
511 av=avma; p1 = gel(x,i);
512 for (j=maxss(0,i-dy1); j<=i && j<=dz; j++)
513 p1 = Fq_sub(p1, Fq_mul(gel(z,j),gel(y,i-j), NULL,p), NULL,p);
514 tetpil=avma; gel(rem,i) = gerepile(av,tetpil, Fq_red(p1, T, p));
515 }
516 rem -= 2;
517 guncloneNULL(lead);
518 if (!sx) (void)FpXQX_renormalize(rem, lr);
519 if (pr == ONLY_REM) return gerepileupto(av0,rem);
520 *pr = rem; return z-2;
521 }
522
523 static GEN
FpXQX_halfgcd_basecase(GEN a,GEN b,GEN T,GEN p)524 FpXQX_halfgcd_basecase(GEN a, GEN b, GEN T, GEN p)
525 {
526 pari_sp av=avma;
527 GEN u,u1,v,v1;
528 long vx = varn(a);
529 long n = lgpol(a)>>1;
530 u1 = v = pol_0(vx);
531 u = v1 = pol_1(vx);
532 while (lgpol(b)>n)
533 {
534 GEN r, q = FpXQX_divrem(a,b, T, p, &r);
535 a = b; b = r; swap(u,u1); swap(v,v1);
536 u1 = FpXX_sub(u1, FpXQX_mul(u, q, T, p), p);
537 v1 = FpXX_sub(v1, FpXQX_mul(v, q ,T, p), p);
538 if (gc_needed(av,2))
539 {
540 if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_halfgcd (d = %ld)",degpol(b));
541 gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
542 }
543 }
544 return gerepilecopy(av, mkmat2(mkcol2(u,u1), mkcol2(v,v1)));
545 }
546 static GEN
FpXQX_addmulmul(GEN u,GEN v,GEN x,GEN y,GEN T,GEN p)547 FpXQX_addmulmul(GEN u, GEN v, GEN x, GEN y, GEN T, GEN p)
548 {
549 return FpXX_add(FpXQX_mul(u, x, T, p),FpXQX_mul(v, y, T, p), p);
550 }
551
552 static GEN
FpXQXM_FpXQX_mul2(GEN M,GEN x,GEN y,GEN T,GEN p)553 FpXQXM_FpXQX_mul2(GEN M, GEN x, GEN y, GEN T, GEN p)
554 {
555 GEN res = cgetg(3, t_COL);
556 gel(res, 1) = FpXQX_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, T, p);
557 gel(res, 2) = FpXQX_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, T, p);
558 return res;
559 }
560
561 static GEN
FpXQXM_mul2(GEN A,GEN B,GEN T,GEN p)562 FpXQXM_mul2(GEN A, GEN B, GEN T, GEN p)
563 {
564 GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
565 GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
566 GEN M1 = FpXQX_mul(FpXX_add(A11,A22, p), FpXX_add(B11,B22, p), T, p);
567 GEN M2 = FpXQX_mul(FpXX_add(A21,A22, p), B11, T, p);
568 GEN M3 = FpXQX_mul(A11, FpXX_sub(B12,B22, p), T, p);
569 GEN M4 = FpXQX_mul(A22, FpXX_sub(B21,B11, p), T, p);
570 GEN M5 = FpXQX_mul(FpXX_add(A11,A12, p), B22, T, p);
571 GEN M6 = FpXQX_mul(FpXX_sub(A21,A11, p), FpXX_add(B11,B12, p), T, p);
572 GEN M7 = FpXQX_mul(FpXX_sub(A12,A22, p), FpXX_add(B21,B22, p), T, p);
573 GEN T1 = FpXX_add(M1,M4, p), T2 = FpXX_sub(M7,M5, p);
574 GEN T3 = FpXX_sub(M1,M2, p), T4 = FpXX_add(M3,M6, p);
575 retmkmat2(mkcol2(FpXX_add(T1,T2, p), FpXX_add(M2,M4, p)),
576 mkcol2(FpXX_add(M3,M5, p), FpXX_add(T3,T4, p)));
577 }
578 /* Return [0,1;1,-q]*M */
579 static GEN
FpXQX_FpXQXM_qmul(GEN q,GEN M,GEN T,GEN p)580 FpXQX_FpXQXM_qmul(GEN q, GEN M, GEN T, GEN p)
581 {
582 GEN u, v, res = cgetg(3, t_MAT);
583 u = FpXX_sub(gcoeff(M,1,1), FpXQX_mul(gcoeff(M,2,1), q, T, p), p);
584 gel(res,1) = mkcol2(gcoeff(M,2,1), u);
585 v = FpXX_sub(gcoeff(M,1,2), FpXQX_mul(gcoeff(M,2,2), q, T, p), p);
586 gel(res,2) = mkcol2(gcoeff(M,2,2), v);
587 return res;
588 }
589
590 static GEN
matid2_FpXQXM(long v)591 matid2_FpXQXM(long v)
592 {
593 retmkmat2(mkcol2(pol_1(v),pol_0(v)),
594 mkcol2(pol_0(v),pol_1(v)));
595 }
596
597 static GEN
FpXX_shift(GEN a,long n)598 FpXX_shift(GEN a, long n) { return RgX_shift_shallow(a, n); }
599
600 static GEN
FpXQX_halfgcd_split(GEN x,GEN y,GEN T,GEN p)601 FpXQX_halfgcd_split(GEN x, GEN y, GEN T, GEN p)
602 {
603 pari_sp av=avma;
604 GEN R, S, V;
605 GEN y1, r, q;
606 long l = lgpol(x), n = l>>1, k;
607 if (lgpol(y)<=n) return matid2_FpXQXM(varn(x));
608 R = FpXQX_halfgcd(FpXX_shift(x,-n),FpXX_shift(y,-n), T, p);
609 V = FpXQXM_FpXQX_mul2(R,x,y, T, p); y1 = gel(V,2);
610 if (lgpol(y1)<=n) return gerepilecopy(av, R);
611 q = FpXQX_divrem(gel(V,1), y1, T, p, &r);
612 k = 2*n-degpol(y1);
613 S = FpXQX_halfgcd(FpXX_shift(y1,-k), FpXX_shift(r,-k), T, p);
614 return gerepileupto(av, FpXQXM_mul2(S,FpXQX_FpXQXM_qmul(q,R, T, p), T, p));
615 }
616
617 /* Return M in GL_2(Fp[X]) such that:
618 if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
619 */
620
621 static GEN
FpXQX_halfgcd_i(GEN x,GEN y,GEN T,GEN p)622 FpXQX_halfgcd_i(GEN x, GEN y, GEN T, GEN p)
623 {
624 if (lg(x)<=FpXQX_HALFGCD_LIMIT) return FpXQX_halfgcd_basecase(x, y, T, p);
625 return FpXQX_halfgcd_split(x, y, T, p);
626 }
627
628 GEN
FpXQX_halfgcd(GEN x,GEN y,GEN T,GEN p)629 FpXQX_halfgcd(GEN x, GEN y, GEN T, GEN p)
630 {
631 pari_sp av = avma;
632 GEN M,q,r;
633 if (lgefint(p)==3)
634 {
635 ulong pp = to_FlxqX(x, y, T, p, &x, &y, &T);
636 M = FlxXM_to_ZXXM(FlxqX_halfgcd(x, y, T, pp));
637 }
638 else
639 {
640 if (!signe(x))
641 {
642 long v = varn(x);
643 retmkmat2(mkcol2(pol_0(v),pol_1(v)),
644 mkcol2(pol_1(v),pol_0(v)));
645 }
646 if (degpol(y)<degpol(x)) return FpXQX_halfgcd_i(x, y, T, p);
647 q = FpXQX_divrem(y, x, T, p, &r);
648 M = FpXQX_halfgcd_i(x, r, T, p);
649 gcoeff(M,1,1) = FpXX_sub(gcoeff(M,1,1), FpXQX_mul(q, gcoeff(M,1,2), T, p), p);
650 gcoeff(M,2,1) = FpXX_sub(gcoeff(M,2,1), FpXQX_mul(q, gcoeff(M,2,2), T, p), p);
651 }
652 return gerepilecopy(av, M);
653 }
654
655 static GEN
FpXQX_gcd_basecase(GEN a,GEN b,GEN T,GEN p)656 FpXQX_gcd_basecase(GEN a, GEN b, GEN T, GEN p)
657 {
658 pari_sp av = avma, av0=avma;
659 while (signe(b))
660 {
661 GEN c;
662 if (gc_needed(av0,2))
663 {
664 if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_gcd (d = %ld)",degpol(b));
665 gerepileall(av0,2, &a,&b);
666 }
667 av = avma; c = FpXQX_rem(a, b, T, p); a=b; b=c;
668 }
669 return gc_const(av, a);
670 }
671
672 GEN
FpXQX_gcd(GEN x,GEN y,GEN T,GEN p)673 FpXQX_gcd(GEN x, GEN y, GEN T, GEN p)
674 {
675 pari_sp av = avma;
676 if (lgefint(p) == 3)
677 {
678 GEN Pl, Ql, Tl, U;
679 ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
680 U = FlxqX_gcd(Pl, Ql, Tl, pp);
681 return gerepileupto(av, FlxX_to_ZXX(U));
682 }
683 x = FpXQX_red(x, T, p);
684 y = FpXQX_red(y, T, p);
685 if (!signe(x)) return gerepileupto(av, y);
686 while (lg(y)>FpXQX_GCD_LIMIT)
687 {
688 GEN c;
689 if (lgpol(y)<=(lgpol(x)>>1))
690 {
691 GEN r = FpXQX_rem(x, y, T, p);
692 x = y; y = r;
693 }
694 c = FpXQXM_FpXQX_mul2(FpXQX_halfgcd(x,y, T, p), x, y, T, p);
695 x = gel(c,1); y = gel(c,2);
696 gerepileall(av,2,&x,&y);
697 }
698 return gerepileupto(av, FpXQX_gcd_basecase(x, y, T, p));
699 }
700
701 static GEN
FpXQX_extgcd_basecase(GEN a,GEN b,GEN T,GEN p,GEN * ptu,GEN * ptv)702 FpXQX_extgcd_basecase(GEN a, GEN b, GEN T, GEN p, GEN *ptu, GEN *ptv)
703 {
704 pari_sp av=avma;
705 GEN u,v,d,d1,v1;
706 long vx = varn(a);
707 d = a; d1 = b;
708 v = pol_0(vx); v1 = pol_1(vx);
709 while (signe(d1))
710 {
711 GEN r, q = FpXQX_divrem(d, d1, T, p, &r);
712 v = FpXX_sub(v,FpXQX_mul(q,v1,T, p),p);
713 u=v; v=v1; v1=u;
714 u=r; d=d1; d1=u;
715 if (gc_needed(av,2))
716 {
717 if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_extgcd (d = %ld)",degpol(d));
718 gerepileall(av,5, &d,&d1,&u,&v,&v1);
719 }
720 }
721 if (ptu) *ptu = FpXQX_div(FpXX_sub(d,FpXQX_mul(b,v, T, p), p), a, T, p);
722 *ptv = v; return d;
723 }
724
725 static GEN
FpXQX_extgcd_halfgcd(GEN x,GEN y,GEN T,GEN p,GEN * ptu,GEN * ptv)726 FpXQX_extgcd_halfgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
727 {
728 pari_sp av=avma;
729 GEN u,v,R = matid2_FpXQXM(varn(x));
730 while (lg(y)>FpXQX_EXTGCD_LIMIT)
731 {
732 GEN M, c;
733 if (lgpol(y)<=(lgpol(x)>>1))
734 {
735 GEN r, q = FpXQX_divrem(x, y, T, p, &r);
736 x = y; y = r;
737 R = FpXQX_FpXQXM_qmul(q, R, T, p);
738 }
739 M = FpXQX_halfgcd(x,y, T, p);
740 c = FpXQXM_FpXQX_mul2(M, x,y, T, p);
741 R = FpXQXM_mul2(M, R, T, p);
742 x = gel(c,1); y = gel(c,2);
743 gerepileall(av,3,&x,&y,&R);
744 }
745 y = FpXQX_extgcd_basecase(x,y, T, p, &u,&v);
746 if (ptu) *ptu = FpXQX_addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1), T, p);
747 *ptv = FpXQX_addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2), T, p);
748 return y;
749 }
750
751 /* x and y in Z[Y][X], return lift(gcd(x mod T,p, y mod T,p)). Set u and v st
752 * ux + vy = gcd (mod T,p) */
753 GEN
FpXQX_extgcd(GEN x,GEN y,GEN T,GEN p,GEN * ptu,GEN * ptv)754 FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv)
755 {
756 GEN d;
757 pari_sp ltop=avma;
758 if (lgefint(p) == 3)
759 {
760 GEN Pl, Ql, Tl, Dl;
761 ulong pp = to_FlxqX(x, y, T, p, &Pl, &Ql, &Tl);
762 Dl = FlxqX_extgcd(Pl, Ql, Tl, pp, ptu, ptv);
763 if (ptu) *ptu = FlxX_to_ZXX(*ptu);
764 *ptv = FlxX_to_ZXX(*ptv);
765 d = FlxX_to_ZXX(Dl);
766 }
767 else
768 {
769 x = FpXQX_red(x, T, p);
770 y = FpXQX_red(y, T, p);
771 if (lg(y)>FpXQX_EXTGCD_LIMIT)
772 d = FpXQX_extgcd_halfgcd(x, y, T, p, ptu, ptv);
773 else
774 d = FpXQX_extgcd_basecase(x, y, T, p, ptu, ptv);
775 }
776 gerepileall(ltop,ptu?3:2,&d,ptv,ptu);
777 return d;
778 }
779
780 GEN
FpXQX_dotproduct(GEN x,GEN y,GEN T,GEN p)781 FpXQX_dotproduct(GEN x, GEN y, GEN T, GEN p)
782 {
783 long i, l = minss(lg(x), lg(y));
784 pari_sp av;
785 GEN c;
786 if (l == 2) return gen_0;
787 av = avma; c = gmul(gel(x,2),gel(y,2));
788 for (i=3; i<l; i++) c = gadd(c, gmul(gel(x,i),gel(y,i)));
789 return gerepileupto(av, Fq_red(c,T,p));
790 }
791
792 /* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab), with R=A%B, a=deg(A) ...*/
793 GEN
FpXQX_resultant(GEN a,GEN b,GEN T,GEN p)794 FpXQX_resultant(GEN a, GEN b, GEN T, GEN p)
795 {
796 long da,db,dc;
797 pari_sp av;
798 long vT = get_FpX_var(T);
799 GEN c,lb, res = pol_1(vT);
800
801 if (!signe(a) || !signe(b)) return pol_0(vT);
802 if (lgefint(p) == 3)
803 {
804 pari_sp av = avma;
805 GEN Pl, Ql, Tl, R;
806 ulong pp = to_FlxqX(a, b, T, p, &Pl, &Ql, &Tl);
807 R = FlxqX_resultant(Pl, Ql, Tl, pp);
808 return gerepileupto(av, Flx_to_ZX(R));
809 }
810
811 da = degpol(a);
812 db = degpol(b);
813 if (db > da)
814 {
815 swapspec(a,b, da,db);
816 if (both_odd(da,db)) res = FpX_neg(res, p);
817 }
818 if (!da) return pol_1(vT); /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
819 av = avma;
820 while (db)
821 {
822 lb = gel(b,db+2);
823 c = FpXQX_rem(a,b, T,p);
824 a = b; b = c; dc = degpol(c);
825 if (dc < 0) { set_avma(av); return pol_0(vT); }
826
827 if (both_odd(da,db)) res = FpX_neg(res, p);
828 if (!equali1(lb)) res = FpXQ_mul(res, FpXQ_powu(lb, da - dc, T, p), T, p);
829 if (gc_needed(av,2))
830 {
831 if (DEBUGMEM>1) pari_warn(warnmem,"FpXQX_resultant (da = %ld)",da);
832 gerepileall(av,3, &a,&b,&res);
833 }
834 da = db; /* = degpol(a) */
835 db = dc; /* = degpol(b) */
836 }
837 res = FpXQ_mul(res, FpXQ_powu(gel(b,2), da, T, p), T, p);
838 return gerepileupto(av, res);
839 }
840
841 /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
842 GEN
FpXQX_disc(GEN P,GEN T,GEN p)843 FpXQX_disc(GEN P, GEN T, GEN p)
844 {
845 pari_sp av = avma;
846 GEN L, dP = FpXX_deriv(P, p), D = FpXQX_resultant(P, dP, T, p);
847 long dd;
848 if (!signe(D)) return pol_0(get_FpX_var(T));
849 dd = degpol(P) - 2 - degpol(dP); /* >= -1; > -1 iff p | deg(P) */
850 L = leading_coeff(P);
851 if (dd && !gequal1(L))
852 D = (dd == -1)? FpXQ_div(D,L,T,p): FpXQ_mul(D, FpXQ_powu(L, dd, T, p), T, p);
853 if (degpol(P) & 2) D = FpX_neg(D, p);
854 return gerepileupto(av, D);
855 }
856
857 /***********************************************************************/
858 /** **/
859 /** Barrett reduction **/
860 /** **/
861 /***********************************************************************/
862
863 /* Return new lgpol */
864 static long
ZXX_lgrenormalizespec(GEN x,long lx)865 ZXX_lgrenormalizespec(GEN x, long lx)
866 {
867 long i;
868 for (i = lx-1; i>=0; i--)
869 if (signe(gel(x,i))) break;
870 return i+1;
871 }
872
873 static GEN
FpXQX_invBarrett_basecase(GEN S,GEN T,GEN p)874 FpXQX_invBarrett_basecase(GEN S, GEN T, GEN p)
875 {
876 long i, l=lg(S)-1, lr = l-1, k;
877 GEN r=cgetg(lr, t_POL); r[1]=S[1];
878 gel(r,2) = gen_1;
879 for (i=3; i<lr; i++)
880 {
881 pari_sp av = avma;
882 GEN u = gel(S,l-i+2);
883 for (k=3; k<i; k++)
884 u = Fq_add(u, Fq_mul(gel(S,l-i+k), gel(r,k), NULL, p), NULL, p);
885 gel(r,i) = gerepileupto(av, Fq_red(Fq_neg(u, NULL, p), T, p));
886 }
887 return FpXQX_renormalize(r,lr);
888 }
889
890 INLINE GEN
FpXX_recipspec(GEN x,long l,long n)891 FpXX_recipspec(GEN x, long l, long n)
892 {
893 return RgX_recipspec_shallow(x, l, n);
894 }
895
896 static GEN
FpXQX_invBarrett_Newton(GEN S,GEN T,GEN p)897 FpXQX_invBarrett_Newton(GEN S, GEN T, GEN p)
898 {
899 pari_sp av = avma;
900 long nold, lx, lz, lq, l = degpol(S), i, lQ;
901 GEN q, y, z, x = cgetg(l+2, t_POL) + 2;
902 ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
903 for (i=0;i<l;i++) gel(x,i) = gen_0;
904 q = RgX_recipspec_shallow(S+2,l+1,l+1); lQ = lgpol(q); q+=2;
905 /* We work on _spec_ FpX's, all the l[xzq] below are lgpol's */
906
907 /* initialize */
908 gel(x,0) = Fq_inv(gel(q,0), T, p);
909 if (lQ>1) gel(q,1) = Fq_red(gel(q,1), T, p);
910 if (lQ>1 && signe(gel(q,1)))
911 {
912 GEN u = gel(q, 1);
913 if (!gequal1(gel(x,0))) u = Fq_mul(u, Fq_sqr(gel(x,0), T, p), T, p);
914 gel(x,1) = Fq_neg(u, T, p); lx = 2;
915 }
916 else
917 lx = 1;
918 nold = 1;
919 for (; mask > 1; )
920 { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
921 long i, lnew, nnew = nold << 1;
922
923 if (mask & 1) nnew--;
924 mask >>= 1;
925
926 lnew = nnew + 1;
927 lq = ZXX_lgrenormalizespec(q, minss(lQ,lnew));
928 z = FpXQX_mulspec(x, q, T, p, lx, lq); /* FIXME: high product */
929 lz = lgpol(z); if (lz > lnew) lz = lnew;
930 z += 2;
931 /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
932 for (i = nold; i < lz; i++) if (signe(gel(z,i))) break;
933 nold = nnew;
934 if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
935
936 /* z + i represents (x*q - 1) / t^i */
937 lz = ZXX_lgrenormalizespec (z+i, lz-i);
938 z = FpXQX_mulspec(x, z+i, T, p, lx, lz); /* FIXME: low product */
939 lz = lgpol(z); z += 2;
940 if (lz > lnew-i) lz = ZXX_lgrenormalizespec(z, lnew-i);
941
942 lx = lz+ i;
943 y = x + i; /* x -= z * t^i, in place */
944 for (i = 0; i < lz; i++) gel(y,i) = Fq_neg(gel(z,i), T, p);
945 }
946 x -= 2; setlg(x, lx + 2); x[1] = S[1];
947 return gerepilecopy(av, x);
948 }
949
950 GEN
FpXQX_invBarrett(GEN S,GEN T,GEN p)951 FpXQX_invBarrett(GEN S, GEN T, GEN p)
952 {
953 pari_sp ltop = avma;
954 long l = lg(S);
955 GEN r;
956 if (l<5) return pol_0(varn(S));
957 if (l<=FpXQX_INVBARRETT_LIMIT)
958 {
959 GEN c = gel(S,l-1), ci=gen_1;
960 if (!gequal1(c))
961 {
962 ci = Fq_inv(c, T, p);
963 S = FqX_Fq_mul(S, ci, T, p);
964 r = FpXQX_invBarrett_basecase(S, T, p);
965 r = FqX_Fq_mul(r, ci, T, p);
966 } else
967 r = FpXQX_invBarrett_basecase(S, T, p);
968 }
969 else
970 r = FpXQX_invBarrett_Newton(S, T, p);
971 return gerepileupto(ltop, r);
972 }
973
974 GEN
FpXQX_get_red(GEN S,GEN T,GEN p)975 FpXQX_get_red(GEN S, GEN T, GEN p)
976 {
977 if (typ(S)==t_POL && lg(S)>FpXQX_BARRETT_LIMIT)
978 retmkvec2(FpXQX_invBarrett(S,T,p),S);
979 return S;
980 }
981
982 /* Compute x mod S where 2 <= degpol(S) <= l+1 <= 2*(degpol(S)-1)
983 * and mg is the Barrett inverse of S. */
984 static GEN
FpXQX_divrem_Barrettspec(GEN x,long l,GEN mg,GEN S,GEN T,GEN p,GEN * pr)985 FpXQX_divrem_Barrettspec(GEN x, long l, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
986 {
987 GEN q, r;
988 long lt = degpol(S); /*We discard the leading term*/
989 long ld, lm, lT, lmg;
990 ld = l-lt;
991 lm = minss(ld, lgpol(mg));
992 lT = ZXX_lgrenormalizespec(S+2,lt);
993 lmg = ZXX_lgrenormalizespec(mg+2,lm);
994 q = FpXX_recipspec(x+lt,ld,ld); /* q = rec(x) lq<=ld*/
995 q = FpXQX_mulspec(q+2,mg+2,T,p,lgpol(q),lmg); /* q = rec(x) * mg lq<=ld+lm*/
996 q = FpXX_recipspec(q+2,minss(ld,lgpol(q)),ld); /* q = rec (rec(x) * mg) lq<=ld*/
997 if (!pr) return q;
998 r = FpXQX_mulspec(q+2,S+2,T,p,lgpol(q),lT); /* r = q*pol lr<=ld+lt*/
999 r = FpXX_subspec(x,r+2,p,lt,minss(lt,lgpol(r))); /* r = x - r lr<=lt */
1000 if (pr == ONLY_REM) return r;
1001 *pr = r; return q;
1002 }
1003
1004 static GEN
FpXQX_divrem_Barrett(GEN x,GEN mg,GEN S,GEN T,GEN p,GEN * pr)1005 FpXQX_divrem_Barrett(GEN x, GEN mg, GEN S, GEN T, GEN p, GEN *pr)
1006 {
1007 GEN q = NULL, r = FpXQX_red(x, T, p);
1008 long l = lgpol(r), lt = degpol(S), lm = 2*lt-1, v = varn(S);
1009 long i;
1010 if (l <= lt)
1011 {
1012 if (pr == ONLY_REM) return r;
1013 if (pr == ONLY_DIVIDES) return signe(r)? NULL: pol_0(v);
1014 if (pr) *pr = r;
1015 return pol_0(v);
1016 }
1017 if (lt <= 1)
1018 return FpXQX_divrem_basecase(r,S,T,p,pr);
1019 if (pr != ONLY_REM && l>lm)
1020 {
1021 q = cgetg(l-lt+2, t_POL); q[1] = S[1];
1022 for (i=0;i<l-lt;i++) gel(q+2,i) = gen_0;
1023 }
1024 while (l>lm)
1025 {
1026 GEN zr, zq = FpXQX_divrem_Barrettspec(r+2+l-lm,lm,mg,S,T,p,&zr);
1027 long lz = lgpol(zr);
1028 if (pr != ONLY_REM)
1029 {
1030 long lq = lgpol(zq);
1031 for(i=0; i<lq; i++) gel(q+2+l-lm,i) = gel(zq,2+i);
1032 }
1033 for(i=0; i<lz; i++) gel(r+2+l-lm,i) = gel(zr,2+i);
1034 l = l-lm+lz;
1035 }
1036 if (pr == ONLY_REM)
1037 {
1038 if (l > lt)
1039 r = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,ONLY_REM);
1040 else
1041 r = FpXQX_renormalize(r, l+2);
1042 setvarn(r, v); return r;
1043 }
1044 if (l > lt)
1045 {
1046 GEN zq = FpXQX_divrem_Barrettspec(r+2,l,mg,S,T,p,pr ? &r: NULL);
1047 if (!q) q = zq;
1048 else
1049 {
1050 long lq = lgpol(zq);
1051 for(i=0; i<lq; i++) gel(q+2,i) = gel(zq,2+i);
1052 }
1053 }
1054 else if (pr)
1055 r = FpX_renormalize(r, l+2);
1056 setvarn(q, v); q = FpXQX_renormalize(q, lg(q));
1057 if (pr == ONLY_DIVIDES) return signe(r)? NULL: q;
1058 if (pr) { setvarn(r, v); *pr = r; }
1059 return q;
1060 }
1061
1062 GEN
FpXQX_divrem(GEN x,GEN S,GEN T,GEN p,GEN * pr)1063 FpXQX_divrem(GEN x, GEN S, GEN T, GEN p, GEN *pr)
1064 {
1065 GEN B, y;
1066 long dy, dx, d;
1067 if (pr==ONLY_REM) return FpXQX_rem(x, S, T, p);
1068 y = get_FpXQX_red(S, &B);
1069 dy = degpol(y); dx = degpol(x); d = dx-dy;
1070 if (lgefint(p) == 3)
1071 {
1072 GEN a, b, t, z;
1073 pari_sp tetpil, av = avma;
1074 ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1075 z = FlxqX_divrem(a, b, t, pp, pr);
1076 if (pr == ONLY_DIVIDES && !z) return gc_NULL(av);
1077 tetpil=avma;
1078 z = FlxX_to_ZXX(z);
1079 if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM)
1080 *pr = FlxX_to_ZXX(*pr);
1081 else return gerepile(av, tetpil, z);
1082 gerepileallsp(av,tetpil,2, pr, &z);
1083 return z;
1084 }
1085 if (!B && d+3 < FpXQX_DIVREM_BARRETT_LIMIT)
1086 return FpXQX_divrem_basecase(x,y,T,p,pr);
1087 else
1088 {
1089 pari_sp av=avma;
1090 GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1091 GEN q = FpXQX_divrem_Barrett(x,mg,y,T,p,pr);
1092 if (!q) return gc_NULL(av);
1093 if (!pr || pr==ONLY_DIVIDES) return gerepilecopy(av, q);
1094 gerepileall(av,2,&q,pr);
1095 return q;
1096 }
1097 }
1098
1099 GEN
FpXQX_rem(GEN x,GEN S,GEN T,GEN p)1100 FpXQX_rem(GEN x, GEN S, GEN T, GEN p)
1101 {
1102 GEN B, y = get_FpXQX_red(S, &B);
1103 long dy = degpol(y), dx = degpol(x), d = dx-dy;
1104 if (d < 0) return FpXQX_red(x, T, p);
1105 if (lgefint(p) == 3)
1106 {
1107 pari_sp av = avma;
1108 GEN a, b, t, z;
1109 ulong pp = to_FlxqX(x, y, T, p, &a, &b, &t);
1110 z = FlxqX_rem(a, b, t, pp);
1111 z = FlxX_to_ZXX(z);
1112 return gerepileupto(av, z);
1113 }
1114 if (!B && d+3 < FpXQX_REM_BARRETT_LIMIT)
1115 return FpXQX_divrem_basecase(x,y, T, p, ONLY_REM);
1116 else
1117 {
1118 pari_sp av=avma;
1119 GEN mg = B? B: FpXQX_invBarrett(y, T, p);
1120 GEN r = FpXQX_divrem_Barrett(x, mg, y, T, p, ONLY_REM);
1121 return gerepileupto(av, r);
1122 }
1123 }
1124
1125 /* x + y*z mod p */
1126 INLINE GEN
Fq_addmul(GEN x,GEN y,GEN z,GEN T,GEN p)1127 Fq_addmul(GEN x, GEN y, GEN z, GEN T, GEN p)
1128 {
1129 pari_sp av;
1130 if (!signe(y) || !signe(z)) return Fq_red(x, T, p);
1131 if (!signe(x)) return Fq_mul(z,y, T, p);
1132 av = avma;
1133 return gerepileupto(av, Fq_add(x, Fq_mul(y, z, T, p), T, p));
1134 }
1135
1136 GEN
FpXQX_div_by_X_x(GEN a,GEN x,GEN T,GEN p,GEN * r)1137 FpXQX_div_by_X_x(GEN a, GEN x, GEN T, GEN p, GEN *r)
1138 {
1139 long l = lg(a), i;
1140 GEN z;
1141 if (l <= 3)
1142 {
1143 if (r) *r = l == 2? gen_0: gcopy(gel(a,2));
1144 return pol_0(0);
1145 }
1146 l--; z = cgetg(l, t_POL); z[1] = evalsigne(1) | evalvarn(0);
1147 gel(z, l-1) = gel(a,l);
1148 for (i=l-2; i>1; i--) /* z[i] = a[i+1] + x*z[i+1] */
1149 gel(z, i) = Fq_addmul(gel(a,i+1), x, gel(z,i+1), T, p);
1150 if (r) *r = Fq_addmul(gel(a,2), x, gel(z,2), T, p);
1151 return z;
1152 }
1153
1154 struct _FpXQXQ {
1155 GEN T, S;
1156 GEN p;
1157 };
1158
_FpXQX_mul(void * data,GEN a,GEN b)1159 static GEN _FpXQX_mul(void *data, GEN a,GEN b)
1160 {
1161 struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1162 return FpXQX_mul(a,b,d->T,d->p);
1163 }
1164
_FpXQX_sqr(void * data,GEN a)1165 static GEN _FpXQX_sqr(void *data, GEN a)
1166 {
1167 struct _FpXQXQ *d=(struct _FpXQXQ*)data;
1168 return FpXQX_sqr(a, d->T, d->p);
1169 }
1170
1171 GEN
FpXQX_powu(GEN x,ulong n,GEN T,GEN p)1172 FpXQX_powu(GEN x, ulong n, GEN T, GEN p)
1173 {
1174 struct _FpXQXQ D;
1175 if (n==0) return pol_1(varn(x));
1176 D.T = T; D.p = p;
1177 return gen_powu(x, n, (void *)&D, _FpXQX_sqr, _FpXQX_mul);
1178 }
1179
1180 GEN
FpXQXV_prod(GEN V,GEN T,GEN p)1181 FpXQXV_prod(GEN V, GEN T, GEN p)
1182 {
1183 if (lgefint(p) == 3)
1184 {
1185 pari_sp av = avma;
1186 ulong pp = p[2];
1187 GEN Tl = ZXT_to_FlxT(T, pp);
1188 GEN Vl = ZXXV_to_FlxXV(V, pp, get_FpX_var(T));
1189 Tl = FlxqXV_prod(Vl, Tl, pp);
1190 return gerepileupto(av, FlxX_to_ZXX(Tl));
1191 }
1192 else
1193 {
1194 struct _FpXQXQ d;
1195 d.T=T; d.p=p;
1196 return gen_product(V, (void*)&d, &_FpXQX_mul);
1197 }
1198 }
1199
1200 static GEN
_FpXQX_divrem(void * E,GEN x,GEN y,GEN * r)1201 _FpXQX_divrem(void * E, GEN x, GEN y, GEN *r)
1202 {
1203 struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1204 return FpXQX_divrem(x, y, d->T, d->p, r);
1205 }
1206
1207 static GEN
_FpXQX_add(void * E,GEN x,GEN y)1208 _FpXQX_add(void * E, GEN x, GEN y)
1209 {
1210 struct _FpXQXQ *d = (struct _FpXQXQ *) E;
1211 return FpXX_add(x, y, d->p);
1212 }
1213
1214 static GEN
_FpXQX_sub(void * E,GEN x,GEN y)1215 _FpXQX_sub(void * E, GEN x, GEN y) {
1216 struct _FpXQXQ *d = (struct _FpXQXQ*) E;
1217 return FpXX_sub(x,y, d->p);
1218 }
1219
1220 static struct bb_ring FpXQX_ring = { _FpXQX_add, _FpXQX_mul, _FpXQX_sqr };
1221
1222 GEN
FpXQX_digits(GEN x,GEN B,GEN T,GEN p)1223 FpXQX_digits(GEN x, GEN B, GEN T, GEN p)
1224 {
1225 pari_sp av = avma;
1226 long d = degpol(B), n = (lgpol(x)+d-1)/d;
1227 GEN z;
1228 struct _FpXQXQ D;
1229 D.T = T; D.p = p;
1230 z = gen_digits(x, B, n, (void *)&D, &FpXQX_ring, _FpXQX_divrem);
1231 return gerepileupto(av, z);
1232 }
1233
1234 GEN
FpXQXV_FpXQX_fromdigits(GEN x,GEN B,GEN T,GEN p)1235 FpXQXV_FpXQX_fromdigits(GEN x, GEN B, GEN T, GEN p)
1236 {
1237 pari_sp av = avma;
1238 struct _FpXQXQ D;
1239 GEN z;
1240 D.T = T; D.p = p;
1241 z = gen_fromdigits(x,B,(void *)&D, &FpXQX_ring);
1242 return gerepileupto(av, z);
1243 }
1244
1245 /* Q an FpXY (t_POL with FpX coeffs), evaluate at X = x */
1246 GEN
FpXY_evalx(GEN Q,GEN x,GEN p)1247 FpXY_evalx(GEN Q, GEN x, GEN p)
1248 {
1249 long i, lb = lg(Q);
1250 GEN z;
1251 z = cgetg(lb, t_POL); z[1] = Q[1];
1252 for (i=2; i<lb; i++)
1253 {
1254 GEN q = gel(Q,i);
1255 gel(z,i) = typ(q) == t_INT? modii(q,p): FpX_eval(q, x, p);
1256 }
1257 return FpX_renormalize(z, lb);
1258 }
1259 /* Q an FpXY, evaluate at Y = y */
1260 GEN
FpXY_evaly(GEN Q,GEN y,GEN p,long vx)1261 FpXY_evaly(GEN Q, GEN y, GEN p, long vx)
1262 {
1263 pari_sp av = avma;
1264 long i, lb = lg(Q);
1265 GEN z;
1266 if (!signe(Q)) return pol_0(vx);
1267 if (lb == 3 || !signe(y)) {
1268 z = gel(Q, 2);
1269 return typ(z)==t_INT? scalar_ZX(z, vx): ZX_copy(z);
1270 }
1271 z = gel(Q, lb-1);
1272 if (typ(z) == t_INT) z = scalar_ZX_shallow(z, vx);
1273 for (i=lb-2; i>=2; i--) z = Fq_add(gel(Q,i), FpX_Fp_mul(z, y, p), NULL, p);
1274 return gerepileupto(av, z);
1275 }
1276 /* Q an FpXY, evaluate at (X,Y) = (x,y) */
1277 GEN
FpXY_eval(GEN Q,GEN y,GEN x,GEN p)1278 FpXY_eval(GEN Q, GEN y, GEN x, GEN p)
1279 {
1280 pari_sp av = avma;
1281 return gerepileuptoint(av, FpX_eval(FpXY_evalx(Q, x, p), y, p));
1282 }
1283
1284 GEN
FpXY_FpXQV_evalx(GEN P,GEN x,GEN T,GEN p)1285 FpXY_FpXQV_evalx(GEN P, GEN x, GEN T, GEN p)
1286 {
1287 long i, lP = lg(P);
1288 GEN res = cgetg(lP,t_POL);
1289 res[1] = P[1];
1290 for(i=2; i<lP; i++)
1291 gel(res,i) = typ(gel(P,i))==t_INT? icopy(gel(P,i)):
1292 FpX_FpXQV_eval(gel(P,i), x, T, p);
1293 return FlxX_renormalize(res, lP);
1294 }
1295
1296 GEN
FpXY_FpXQ_evalx(GEN P,GEN x,GEN T,GEN p)1297 FpXY_FpXQ_evalx(GEN P, GEN x, GEN T, GEN p)
1298 {
1299 pari_sp av = avma;
1300 long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(P),1);
1301 GEN xp = FpXQ_powers(x, n, T, p);
1302 return gerepileupto(av, FpXY_FpXQV_evalx(P, xp, T, p));
1303 }
1304
1305 /*******************************************************************/
1306 /* */
1307 /* (Fp[X]/T(X))[Y] / S(Y) */
1308 /* */
1309 /*******************************************************************/
1310
1311 /*Preliminary implementation to speed up FpX_ffisom*/
1312 typedef struct {
1313 GEN S, T, p;
1314 } FpXYQQ_muldata;
1315
1316 /* reduce x in Fp[X, Y] in the algebra Fp[X,Y]/ (S(X),T(Y)) */
1317 static GEN
FpXYQQ_redswap(GEN x,GEN S,GEN T,GEN p)1318 FpXYQQ_redswap(GEN x, GEN S, GEN T, GEN p)
1319 {
1320 pari_sp ltop=avma;
1321 long n = get_FpX_degree(S);
1322 long m = get_FpX_degree(T);
1323 long v = get_FpX_var(T);
1324 GEN V = RgXY_swap(x,m,v);
1325 V = FpXQX_red(V,S,p);
1326 V = RgXY_swap(V,n,v);
1327 return gerepilecopy(ltop,V);
1328 }
1329 static GEN
FpXYQQ_sqr(void * data,GEN x)1330 FpXYQQ_sqr(void *data, GEN x)
1331 {
1332 FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1333 return FpXYQQ_redswap(FpXQX_sqr(x, D->T, D->p),D->S,D->T,D->p);
1334
1335 }
1336 static GEN
FpXYQQ_mul(void * data,GEN x,GEN y)1337 FpXYQQ_mul(void *data, GEN x, GEN y)
1338 {
1339 FpXYQQ_muldata *D = (FpXYQQ_muldata*)data;
1340 return FpXYQQ_redswap(FpXQX_mul(x,y, D->T, D->p),D->S,D->T,D->p);
1341 }
1342
1343 /* x in Z[X,Y], S in Z[X] over Fq = Z[Y]/(p,T); compute lift(x^n mod (S,T,p)) */
1344 GEN
FpXYQQ_pow(GEN x,GEN n,GEN S,GEN T,GEN p)1345 FpXYQQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1346 {
1347 pari_sp av = avma;
1348 FpXYQQ_muldata D;
1349 GEN y;
1350 if (lgefint(p) == 3)
1351 {
1352 ulong pp = to_FlxqX(x, NULL, T, p, &x, NULL, &T);
1353 S = ZX_to_Flx(S, pp);
1354 y = FlxX_to_ZXX( FlxYqq_pow(x, n, S, T, pp) );
1355 y = gerepileupto(av, y);
1356 }
1357 else
1358 {
1359 D.S = S;
1360 D.T = T;
1361 D.p = p;
1362 y = gen_pow(x, n, (void*)&D, &FpXYQQ_sqr, &FpXYQQ_mul);
1363 }
1364 return y;
1365 }
1366
1367 GEN
FpXQXQ_mul(GEN x,GEN y,GEN S,GEN T,GEN p)1368 FpXQXQ_mul(GEN x, GEN y, GEN S, GEN T, GEN p) {
1369 return FpXQX_rem(FpXQX_mul(x, y, T, p), S, T, p);
1370 }
1371
1372 GEN
FpXQXQ_sqr(GEN x,GEN S,GEN T,GEN p)1373 FpXQXQ_sqr(GEN x, GEN S, GEN T, GEN p) {
1374 return FpXQX_rem(FpXQX_sqr(x, T, p), S, T, p);
1375 }
1376
1377 /* Inverse of x in Z/pZ[X]/(pol) or NULL if inverse doesn't exist
1378 * return lift(1 / (x mod (p,pol))) */
1379 GEN
FpXQXQ_invsafe(GEN x,GEN S,GEN T,GEN p)1380 FpXQXQ_invsafe(GEN x, GEN S, GEN T, GEN p)
1381 {
1382 GEN V, z = FpXQX_extgcd(get_FpXQX_mod(S), x, T, p, NULL, &V);
1383 if (degpol(z)) return NULL;
1384 z = gel(z,2);
1385 z = typ(z)==t_INT ? Fp_invsafe(z,p) : FpXQ_invsafe(z,T,p);
1386 if (!z) return NULL;
1387 return typ(z)==t_INT ? FpXX_Fp_mul(V, z, p): FpXQX_FpXQ_mul(V, z, T, p);
1388 }
1389
1390 GEN
FpXQXQ_inv(GEN x,GEN S,GEN T,GEN p)1391 FpXQXQ_inv(GEN x, GEN S, GEN T,GEN p)
1392 {
1393 pari_sp av = avma;
1394 GEN U = FpXQXQ_invsafe(x, S, T, p);
1395 if (!U) pari_err_INV("FpXQXQ_inv",x);
1396 return gerepileupto(av, U);
1397 }
1398
1399 GEN
FpXQXQ_div(GEN x,GEN y,GEN S,GEN T,GEN p)1400 FpXQXQ_div(GEN x,GEN y,GEN S, GEN T,GEN p)
1401 {
1402 pari_sp av = avma;
1403 return gerepileupto(av, FpXQXQ_mul(x, FpXQXQ_inv(y,S,T,p),S,T,p));
1404 }
1405
1406 static GEN
_FpXQXQ_cmul(void * data,GEN P,long a,GEN x)1407 _FpXQXQ_cmul(void *data, GEN P, long a, GEN x) {
1408 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1409 GEN y = gel(P,a+2);
1410 return typ(y)==t_INT ? FpXX_Fp_mul(x,y, d->p):
1411 FpXX_FpX_mul(x,y,d->p);
1412 }
1413 static GEN
_FpXQXQ_red(void * data,GEN x)1414 _FpXQXQ_red(void *data, GEN x) {
1415 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1416 return FpXQX_red(x, d->T, d->p);
1417 }
1418 static GEN
_FpXQXQ_mul(void * data,GEN x,GEN y)1419 _FpXQXQ_mul(void *data, GEN x, GEN y) {
1420 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1421 return FpXQXQ_mul(x,y, d->S,d->T, d->p);
1422 }
1423 static GEN
_FpXQXQ_sqr(void * data,GEN x)1424 _FpXQXQ_sqr(void *data, GEN x) {
1425 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1426 return FpXQXQ_sqr(x, d->S,d->T, d->p);
1427 }
1428
1429 static GEN
_FpXQXQ_one(void * data)1430 _FpXQXQ_one(void *data) {
1431 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1432 return pol_1(get_FpXQX_var(d->S));
1433 }
1434
1435 static GEN
_FpXQXQ_zero(void * data)1436 _FpXQXQ_zero(void *data) {
1437 struct _FpXQXQ *d = (struct _FpXQXQ*) data;
1438 return pol_0(get_FpXQX_var(d->S));
1439 }
1440
1441 static struct bb_algebra FpXQXQ_algebra = { _FpXQXQ_red, _FpXQX_add,
1442 _FpXQX_sub, _FpXQXQ_mul, _FpXQXQ_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1443
1444 const struct bb_algebra *
get_FpXQXQ_algebra(void ** E,GEN S,GEN T,GEN p)1445 get_FpXQXQ_algebra(void **E, GEN S, GEN T, GEN p)
1446 {
1447 GEN z = new_chunk(sizeof(struct _FpXQXQ));
1448 struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1449 e->T = FpX_get_red(T, p);
1450 e->S = FpXQX_get_red(S, e->T, p);
1451 e->p = p; *E = (void*)e;
1452 return &FpXQXQ_algebra;
1453 }
1454
1455 static struct bb_algebra FpXQX_algebra = { _FpXQXQ_red, _FpXQX_add,
1456 _FpXQX_sub, _FpXQX_mul, _FpXQX_sqr, _FpXQXQ_one, _FpXQXQ_zero };
1457
1458 const struct bb_algebra *
get_FpXQX_algebra(void ** E,GEN T,GEN p,long v)1459 get_FpXQX_algebra(void **E, GEN T, GEN p, long v)
1460 {
1461 GEN z = new_chunk(sizeof(struct _FpXQXQ));
1462 struct _FpXQXQ *e = (struct _FpXQXQ *) z;
1463 e->T = FpX_get_red(T, p);
1464 e->S = pol_x(v);
1465 e->p = p; *E = (void*)e;
1466 return &FpXQX_algebra;
1467 }
1468
1469 /* x over Fq, return lift(x^n) mod S */
1470 GEN
FpXQXQ_pow(GEN x,GEN n,GEN S,GEN T,GEN p)1471 FpXQXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p)
1472 {
1473 pari_sp ltop = avma;
1474 GEN y;
1475 struct _FpXQXQ D;
1476 long s = signe(n);
1477 if (!s) return pol_1(varn(x));
1478 if (is_pm1(n)) /* +/- 1 */
1479 return (s < 0)? FpXQXQ_inv(x,S,T,p): ZXX_copy(x);
1480 if (lgefint(p) == 3)
1481 {
1482 ulong pp = to_FlxqX(x, S, T, p, &x, &S, &T);
1483 GEN z = FlxqXQ_pow(x, n, S, T, pp);
1484 y = FlxX_to_ZXX(z);
1485 return gerepileupto(ltop, y);
1486 }
1487 else
1488 {
1489 T = FpX_get_red(T, p);
1490 S = FpXQX_get_red(S, T, p);
1491 D.S = S; D.T = T; D.p = p;
1492 if (s < 0) x = FpXQXQ_inv(x,S,T,p);
1493 y = gen_pow_i(x, n, (void*)&D,&_FpXQXQ_sqr,&_FpXQXQ_mul);
1494 return gerepilecopy(ltop, y);
1495 }
1496 }
1497
1498 /* generates the list of powers of x of degree 0,1,2,...,l*/
1499 GEN
FpXQXQ_powers(GEN x,long l,GEN S,GEN T,GEN p)1500 FpXQXQ_powers(GEN x, long l, GEN S, GEN T, GEN p)
1501 {
1502 struct _FpXQXQ D;
1503 int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1504 T = FpX_get_red(T, p);
1505 S = FpXQX_get_red(S, T, p);
1506 D.S = S; D.T = T; D.p = p;
1507 return gen_powers(x, l, use_sqr, (void*)&D, &_FpXQXQ_sqr, &_FpXQXQ_mul,&_FpXQXQ_one);
1508 }
1509
1510 /* Let v a linear form, return the linear form z->v(tau*z)
1511 that is, v*(M_tau) */
1512
1513 INLINE GEN
FpXQX_recipspec(GEN x,long l,long n)1514 FpXQX_recipspec(GEN x, long l, long n)
1515 {
1516 return RgX_recipspec_shallow(x, l, n);
1517 }
1518
1519 static GEN
FpXQXQ_transmul_init(GEN tau,GEN S,GEN T,GEN p)1520 FpXQXQ_transmul_init(GEN tau, GEN S, GEN T, GEN p)
1521 {
1522 GEN bht;
1523 GEN h, Sp = get_FpXQX_red(S, &h);
1524 long n = degpol(Sp), vT = varn(Sp);
1525 GEN ft = FpXQX_recipspec(Sp+2, n+1, n+1);
1526 GEN bt = FpXQX_recipspec(tau+2, lgpol(tau), n);
1527 setvarn(ft, vT); setvarn(bt, vT);
1528 if (h)
1529 bht = FpXQXn_mul(bt, h, n-1, T, p);
1530 else
1531 {
1532 GEN bh = FpXQX_div(FpXX_shift(tau, n-1), S, T, p);
1533 bht = FpXQX_recipspec(bh+2, lgpol(bh), n-1);
1534 setvarn(bht, vT);
1535 }
1536 return mkvec3(bt, bht, ft);
1537 }
1538
1539 static GEN
FpXQXQ_transmul(GEN tau,GEN a,long n,GEN T,GEN p)1540 FpXQXQ_transmul(GEN tau, GEN a, long n, GEN T, GEN p)
1541 {
1542 pari_sp ltop = avma;
1543 GEN t1, t2, t3, vec;
1544 GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
1545 if (signe(a)==0) return pol_0(varn(a));
1546 t2 = FpXX_shift(FpXQX_mul(bt, a, T, p),1-n);
1547 if (signe(bht)==0) return gerepilecopy(ltop, t2);
1548 t1 = FpXX_shift(FpXQX_mul(ft, a, T, p),-n);
1549 t3 = FpXQXn_mul(t1, bht, n-1, T, p);
1550 vec = FpXX_sub(t2, FpXX_shift(t3, 1), p);
1551 return gerepileupto(ltop, vec);
1552 }
1553
1554 static GEN
polxn_FpXX(long n,long v,long vT)1555 polxn_FpXX(long n, long v, long vT)
1556 {
1557 long i, a = n+2;
1558 GEN p = cgetg(a+1, t_POL);
1559 p[1] = evalsigne(1)|evalvarn(v);
1560 for (i = 2; i < a; i++) gel(p,i) = pol_0(vT);
1561 gel(p,a) = pol_1(vT); return p;
1562 }
1563
1564 GEN
FpXQXQ_minpoly(GEN x,GEN S,GEN T,GEN p)1565 FpXQXQ_minpoly(GEN x, GEN S, GEN T, GEN p)
1566 {
1567 pari_sp ltop = avma;
1568 long vS, vT, n;
1569 GEN v_x, g, tau;
1570 vS = get_FpXQX_var(S);
1571 vT = get_FpX_var(T);
1572 n = get_FpXQX_degree(S);
1573 g = pol_1(vS);
1574 tau = pol_1(vS);
1575 S = FpXQX_get_red(S, T, p);
1576 v_x = FpXQXQ_powers(x, usqrt(2*n), S, T, p);
1577 while(signe(tau) != 0)
1578 {
1579 long i, j, m, k1;
1580 GEN M, v, tr;
1581 GEN g_prime, c;
1582 if (degpol(g) == n) { tau = pol_1(vS); g = pol_1(vS); }
1583 v = random_FpXQX(n, vS, T, p);
1584 tr = FpXQXQ_transmul_init(tau, S, T, p);
1585 v = FpXQXQ_transmul(tr, v, n, T, p);
1586 m = 2*(n-degpol(g));
1587 k1 = usqrt(m);
1588 tr = FpXQXQ_transmul_init(gel(v_x,k1+1), S, T, p);
1589 c = cgetg(m+2,t_POL);
1590 c[1] = evalsigne(1)|evalvarn(vS);
1591 for (i=0; i<m; i+=k1)
1592 {
1593 long mj = minss(m-i, k1);
1594 for (j=0; j<mj; j++)
1595 gel(c,m+1-(i+j)) = FpXQX_dotproduct(v, gel(v_x,j+1), T, p);
1596 v = FpXQXQ_transmul(tr, v, n, T, p);
1597 }
1598 c = FpXX_renormalize(c, m+2);
1599 /* now c contains <v,x^i> , i = 0..m-1 */
1600 M = FpXQX_halfgcd(polxn_FpXX(m, vS, vT), c, T, p);
1601 g_prime = gmael(M, 2, 2);
1602 if (degpol(g_prime) < 1) continue;
1603 g = FpXQX_mul(g, g_prime, T, p);
1604 tau = FpXQXQ_mul(tau, FpXQX_FpXQXQV_eval(g_prime, v_x, S, T, p), S, T, p);
1605 }
1606 g = FpXQX_normalize(g,T, p);
1607 return gerepilecopy(ltop,g);
1608 }
1609
1610 GEN
FpXQXQ_matrix_pow(GEN y,long n,long m,GEN S,GEN T,GEN p)1611 FpXQXQ_matrix_pow(GEN y, long n, long m, GEN S, GEN T, GEN p)
1612 {
1613 return RgXV_to_RgM(FpXQXQ_powers(y,m-1,S,T,p),n);
1614 }
1615
1616 GEN
FpXQX_FpXQXQV_eval(GEN P,GEN V,GEN S,GEN T,GEN p)1617 FpXQX_FpXQXQV_eval(GEN P, GEN V, GEN S, GEN T, GEN p)
1618 {
1619 struct _FpXQXQ D;
1620 T = FpX_get_red(T, p);
1621 S = FpXQX_get_red(S, T, p);
1622 D.S=S; D.T=T; D.p=p;
1623 return gen_bkeval_powers(P, degpol(P), V, (void*)&D, &FpXQXQ_algebra,
1624 _FpXQXQ_cmul);
1625 }
1626
1627 GEN
FpXQX_FpXQXQ_eval(GEN Q,GEN x,GEN S,GEN T,GEN p)1628 FpXQX_FpXQXQ_eval(GEN Q, GEN x, GEN S, GEN T, GEN p)
1629 {
1630 struct _FpXQXQ D;
1631 int use_sqr = 2*degpol(x) >= get_FpXQX_degree(S);
1632 T = FpX_get_red(T, p);
1633 S = FpXQX_get_red(S, T, p);
1634 D.S=S; D.T=T; D.p=p;
1635 return gen_bkeval(Q, degpol(Q), x, use_sqr, (void*)&D, &FpXQXQ_algebra,
1636 _FpXQXQ_cmul);
1637 }
1638
1639 static GEN
FpXQXQ_autpow_sqr(void * E,GEN x)1640 FpXQXQ_autpow_sqr(void * E, GEN x)
1641 {
1642 struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1643 GEN S = D->S, T = D->T, p = D->p;
1644 GEN phi = gel(x,1), S1 = gel(x,2);
1645 long n = brent_kung_optpow(get_FpX_degree(T)-1,lgpol(S1)+1,1);
1646 GEN V = FpXQ_powers(phi, n, T, p);
1647 GEN phi2 = FpX_FpXQV_eval(phi, V, T, p);
1648 GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1649 GEN S2 = FpXQX_FpXQXQ_eval(Sphi, S1, S, T, p);
1650 return mkvec2(phi2, S2);
1651 }
1652
1653 static GEN
FpXQXQ_autpow_mul(void * E,GEN x,GEN y)1654 FpXQXQ_autpow_mul(void * E, GEN x, GEN y)
1655 {
1656 struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1657 GEN S = D->S, T = D->T, p = D->p;
1658 GEN phi1 = gel(x,1), S1 = gel(x,2);
1659 GEN phi2 = gel(y,1), S2 = gel(y,2);
1660 long n = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+1, 1);
1661 GEN V = FpXQ_powers(phi2, n, T, p);
1662 GEN phi3 = FpX_FpXQV_eval(phi1, V, T, p);
1663 GEN Sphi = FpXY_FpXQV_evalx(S1, V, T, p);
1664 GEN S3 = FpXQX_FpXQXQ_eval(Sphi, S2, S, T, p);
1665 return mkvec2(phi3, S3);
1666 }
1667
1668 GEN
FpXQXQ_autpow(GEN aut,long n,GEN S,GEN T,GEN p)1669 FpXQXQ_autpow(GEN aut, long n, GEN S, GEN T, GEN p)
1670 {
1671 pari_sp av = avma;
1672 struct _FpXQXQ D;
1673 T = FpX_get_red(T, p);
1674 S = FpXQX_get_red(S, T, p);
1675 D.S=S; D.T=T; D.p=p;
1676 aut = gen_powu_i(aut,n,&D,FpXQXQ_autpow_sqr,FpXQXQ_autpow_mul);
1677 return gerepilecopy(av, aut);
1678 }
1679
1680 static GEN
FpXQXQ_auttrace_mul(void * E,GEN x,GEN y)1681 FpXQXQ_auttrace_mul(void *E, GEN x, GEN y)
1682 {
1683 struct _FpXQXQ *D = (struct _FpXQXQ *)E;
1684 GEN S = D->S, T = D->T;
1685 GEN p = D->p;
1686 GEN S1 = gel(x,1), a1 = gel(x,2);
1687 GEN S2 = gel(y,1), a2 = gel(y,2);
1688 long n = brent_kung_optpow(maxss(degpol(S1),degpol(a1)),2,1);
1689 GEN V = FpXQXQ_powers(S2, n, S, T, p);
1690 GEN S3 = FpXQX_FpXQXQV_eval(S1, V, S, T, p);
1691 GEN aS = FpXQX_FpXQXQV_eval(a1, V, S, T, p);
1692 GEN a3 = FpXX_add(aS, a2, p);
1693 return mkvec2(S3, a3);
1694 }
1695
1696 static GEN
FpXQXQ_auttrace_sqr(void * E,GEN x)1697 FpXQXQ_auttrace_sqr(void *E, GEN x)
1698 { return FpXQXQ_auttrace_mul(E, x, x); }
1699
1700 GEN
FpXQXQ_auttrace(GEN aut,long n,GEN S,GEN T,GEN p)1701 FpXQXQ_auttrace(GEN aut, long n, GEN S, GEN T, GEN p)
1702 {
1703 pari_sp av = avma;
1704 struct _FpXQXQ D;
1705 T = FpX_get_red(T, p);
1706 S = FpXQX_get_red(S, T, p);
1707 D.S=S; D.T=T; D.p=p;
1708 aut = gen_powu_i(aut,n,&D,FpXQXQ_auttrace_sqr,FpXQXQ_auttrace_mul);
1709 return gerepilecopy(av, aut);
1710 }
1711
1712 static GEN
FpXQXQ_autsum_mul(void * E,GEN x,GEN y)1713 FpXQXQ_autsum_mul(void *E, GEN x, GEN y)
1714 {
1715 struct _FpXQXQ *D = (struct _FpXQXQ *) E;
1716 GEN S = D->S, T = D->T, p = D->p;
1717 GEN phi1 = gel(x,1), S1 = gel(x,2), a1 = gel(x,3);
1718 GEN phi2 = gel(y,1), S2 = gel(y,2), a2 = gel(y,3);
1719 long n2 = brent_kung_optpow(get_FpX_degree(T)-1, lgpol(S1)+lgpol(a1)+1, 1);
1720 GEN V2 = FpXQ_powers(phi2, n2, T, p);
1721 GEN phi3 = FpX_FpXQV_eval(phi1, V2, T, p);
1722 GEN Sphi = FpXY_FpXQV_evalx(S1, V2, T, p);
1723 GEN aphi = FpXY_FpXQV_evalx(a1, V2, T, p);
1724 long n = brent_kung_optpow(maxss(degpol(Sphi),degpol(aphi)),2,1);
1725 GEN V = FpXQXQ_powers(S2, n, S, T, p);
1726 GEN S3 = FpXQX_FpXQXQV_eval(Sphi, V, S, T, p);
1727 GEN aS = FpXQX_FpXQXQV_eval(aphi, V, S, T, p);
1728 GEN a3 = FpXQXQ_mul(aS, a2, S, T, p);
1729 return mkvec3(phi3, S3, a3);
1730 }
1731
1732 static GEN
FpXQXQ_autsum_sqr(void * T,GEN x)1733 FpXQXQ_autsum_sqr(void * T, GEN x)
1734 { return FpXQXQ_autsum_mul(T,x,x); }
1735
1736 GEN
FpXQXQ_autsum(GEN aut,long n,GEN S,GEN T,GEN p)1737 FpXQXQ_autsum(GEN aut, long n, GEN S, GEN T, GEN p)
1738 {
1739 pari_sp av = avma;
1740 struct _FpXQXQ D;
1741 T = FpX_get_red(T, p);
1742 S = FpXQX_get_red(S, T, p);
1743 D.S=S; D.T=T; D.p=p;
1744 aut = gen_powu_i(aut,n,&D,FpXQXQ_autsum_sqr,FpXQXQ_autsum_mul);
1745 return gerepilecopy(av, aut);
1746 }
1747
1748 GEN
FpXQXn_mul(GEN x,GEN y,long n,GEN T,GEN p)1749 FpXQXn_mul(GEN x, GEN y, long n, GEN T, GEN p)
1750 {
1751 pari_sp av = avma;
1752 GEN z, kx, ky;
1753 long d;
1754 if (ZXX_is_ZX(y) && ZXX_is_ZX(x))
1755 return FpXn_mul(x,y,n,p);
1756 d = get_FpX_degree(T);
1757 kx = ZXX_to_Kronecker(x, d);
1758 ky = ZXX_to_Kronecker(y, d);
1759 z = Kronecker_to_FpXQX(ZXn_mul(ky,kx,(2*d-1)*n), T, p);
1760 return gerepileupto(av, z);
1761 }
1762
1763 GEN
FpXQXn_sqr(GEN x,long n,GEN T,GEN p)1764 FpXQXn_sqr(GEN x, long n, GEN T, GEN p)
1765 {
1766 pari_sp av = avma;
1767 GEN z, kx;
1768 long d;
1769 if (ZXX_is_ZX(x)) return ZXn_sqr(x, n);
1770 d = get_FpX_degree(T);
1771 kx= ZXX_to_Kronecker(x, d);
1772 z = Kronecker_to_FpXQX(ZXn_sqr(kx, (2*d-1)*n), T, p);
1773 return gerepileupto(av, z);
1774 }
1775
1776 INLINE GEN
FpXXn_red(GEN a,long n)1777 FpXXn_red(GEN a, long n)
1778 { return RgXn_red_shallow(a, n); }
1779
1780 /* (f*g) \/ x^n */
1781 static GEN
FpXQX_mulhigh_i(GEN f,GEN g,long n,GEN T,GEN p)1782 FpXQX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
1783 {
1784 return FpXX_shift(FpXQX_mul(f,g,T, p),-n);
1785 }
1786
1787 static GEN
FpXQXn_mulhigh(GEN f,GEN g,long n2,long n,GEN T,GEN p)1788 FpXQXn_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
1789 {
1790 GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
1791 return FpXX_add(FpXQX_mulhigh_i(fl, g, n2, T, p), FpXQXn_mul(fh, g, n - n2, T, p), p);
1792 }
1793
1794 /* Compute intformal(x^n*S)/x^(n+1) */
1795 static GEN
FpXX_integXn(GEN x,long n,GEN p)1796 FpXX_integXn(GEN x, long n, GEN p)
1797 {
1798 long i, lx = lg(x);
1799 GEN y;
1800 if (lx == 2) return ZXX_copy(x);
1801 y = cgetg(lx, t_POL); y[1] = x[1];
1802 for (i=2; i<lx; i++)
1803 {
1804 ulong j = n+i-1;
1805 GEN xi = gel(x,i);
1806 if (!signe(xi))
1807 gel(y,i) = gen_0;
1808 else
1809 gel(y,i) = typ(xi)==t_INT ? Fp_divu(xi, j, p)
1810 : FpX_divu(xi, j, p);
1811 }
1812 return ZXX_renormalize(y, lx);;
1813 }
1814
1815 /* Compute intformal(x^n*S)/x^(n+1) */
1816 static GEN
ZlXX_integXn(GEN x,long n,GEN p,ulong pp)1817 ZlXX_integXn(GEN x, long n, GEN p, ulong pp)
1818 {
1819 long i, lx = lg(x);
1820 GEN y;
1821 if (lx == 2) return ZXX_copy(x);
1822 if (!pp) return FpXX_integXn(x, n, p);
1823 y = cgetg(lx, t_POL); y[1] = x[1];
1824 for (i=2; i<lx; i++)
1825 {
1826 GEN xi = gel(x,i);
1827 if (!signe(xi))
1828 gel(y,i) = gen_0;
1829 else
1830 {
1831 ulong j;
1832 long v = u_lvalrem(n+i-1, pp, &j);
1833 if (typ(xi)==t_INT)
1834 {
1835 if (v==0)
1836 gel(y,i) = Fp_divu(xi, j, p);
1837 else
1838 gel(y,i) = Fp_divu(diviuexact(xi, upowuu(pp, v)), j, p);
1839 } else
1840 {
1841 if (v==0)
1842 gel(y,i) = FpX_divu(xi, j, p);
1843 else
1844 gel(y,i) = FpX_divu(ZX_divuexact(xi, upowuu(pp, v)), j, p);
1845 }
1846 }
1847 }
1848 return ZXX_renormalize(y, lx);;
1849 }
1850
1851 GEN
ZlXQXn_expint(GEN h,long e,GEN T,GEN p,ulong pp)1852 ZlXQXn_expint(GEN h, long e, GEN T, GEN p, ulong pp)
1853 {
1854 pari_sp av = avma, av2;
1855 long v = varn(h), n=1;
1856 GEN f = pol_1(v), g = pol_1(v);
1857 ulong mask = quadratic_prec_mask(e);
1858 av2 = avma;
1859 for (;mask>1;)
1860 {
1861 GEN u, w;
1862 long n2 = n;
1863 n<<=1; if (mask & 1) n--;
1864 mask >>= 1;
1865 u = FpXQXn_mul(g, FpXQX_mulhigh_i(f, FpXXn_red(h, n2-1), n2-1, T, p), n-n2, T, p);
1866 u = FpXX_add(u, FpXX_shift(FpXXn_red(h, n-1), 1-n2), p);
1867 w = FpXQXn_mul(f, ZlXX_integXn(u, n2-1, p, pp), n-n2, T, p);
1868 f = FpXX_add(f, FpXX_shift(w, n2), p);
1869 if (mask<=1) break;
1870 u = FpXQXn_mul(g, FpXQXn_mulhigh(f, g, n2, n, T, p), n-n2, T, p);
1871 g = FpXX_sub(g, FpXX_shift(u, n2), p);
1872 if (gc_needed(av2,2))
1873 {
1874 if (DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_exp, e = %ld", n);
1875 gerepileall(av2, 2, &f, &g);
1876 }
1877 }
1878 return gerepileupto(av, f);
1879 }
1880
1881 GEN
FpXQXn_expint(GEN h,long e,GEN T,GEN p)1882 FpXQXn_expint(GEN h, long e, GEN T, GEN p)
1883 { return ZlXQXn_expint(h, e, T, p, 0); }
1884
1885 GEN
FpXQXn_exp(GEN h,long e,GEN T,GEN p)1886 FpXQXn_exp(GEN h, long e, GEN T, GEN p)
1887 {
1888 if (signe(h)==0 || degpol(h)<1 || !gequal0(gel(h,2)))
1889 pari_err_DOMAIN("FpXQXn_exp","valuation", "<", gen_1, h);
1890 return FpXQXn_expint(FpXX_deriv(h, p), e, T, p);
1891 }
1892
1893 GEN
FpXQXn_inv(GEN f,long e,GEN T,GEN p)1894 FpXQXn_inv(GEN f, long e, GEN T, GEN p)
1895 {
1896 pari_sp av = avma, av2;
1897 ulong mask;
1898 GEN W, a;
1899 long v = varn(f), n = 1;
1900
1901 if (!signe(f)) pari_err_INV("FpXXn_inv",f);
1902 a = Fq_inv(gel(f,2), T, p);
1903 if (e == 1) return scalarpol(a, v);
1904 else if (e == 2)
1905 {
1906 GEN b;
1907 if (degpol(f) <= 0) return scalarpol(a, v);
1908 b = Fq_neg(gel(f,3),T,p);
1909 if (signe(b)==0) return scalarpol(a, v);
1910 if (!is_pm1(a)) b = Fq_mul(b, Fq_sqr(a, T, p), T, p);
1911 W = deg1pol_shallow(b, a, v);
1912 return gerepilecopy(av, W);
1913 }
1914 W = scalarpol_shallow(Fq_inv(gel(f,2), T, p),v);
1915 mask = quadratic_prec_mask(e);
1916 av2 = avma;
1917 for (;mask>1;)
1918 {
1919 GEN u, fr;
1920 long n2 = n;
1921 n<<=1; if (mask & 1) n--;
1922 mask >>= 1;
1923 fr = FpXXn_red(f, n);
1924 u = FpXQXn_mul(W, FpXQXn_mulhigh(fr, W, n2, n, T, p), n-n2, T, p);
1925 W = FpXX_sub(W, FpXX_shift(u, n2), p);
1926 if (gc_needed(av2,2))
1927 {
1928 if(DEBUGMEM>1) pari_warn(warnmem,"FpXQXn_inv, e = %ld", n);
1929 W = gerepileupto(av2, W);
1930 }
1931 }
1932 return gerepileupto(av, W);
1933 }
1934