1 /* Copyright (C) 2004 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 with small coefficients. */
19
20 static GEN
get_Flx_red(GEN T,GEN * B)21 get_Flx_red(GEN T, GEN *B)
22 {
23 if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 *B = gel(T,1); return gel(T,2);
25 }
26
27 /***********************************************************************/
28 /** **/
29 /** Flx **/
30 /** **/
31 /***********************************************************************/
32 /* Flx objects are defined as follows:
33 Let l an ulong. An Flx is a t_VECSMALL:
34 x[0] = codeword
35 x[1] = evalvarn(variable number) (signe is not stored).
36 x[2] = a_0 x[3] = a_1, etc.
37 With 0 <= a_i < l
38
39 signe(x) is not valid. Use degpol(x)>=0 instead.
40 */
41 /***********************************************************************/
42 /** **/
43 /** Conversion from Flx **/
44 /** **/
45 /***********************************************************************/
46
47 GEN
Flx_to_ZX(GEN z)48 Flx_to_ZX(GEN z)
49 {
50 long i, l = lg(z);
51 GEN x = cgetg(l,t_POL);
52 for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
53 x[1] = evalsigne(l-2!=0)| z[1]; return x;
54 }
55
56 GEN
Flx_to_FlxX(GEN z,long sv)57 Flx_to_FlxX(GEN z, long sv)
58 {
59 long i, l = lg(z);
60 GEN x = cgetg(l,t_POL);
61 for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
62 x[1] = evalsigne(l-2!=0)| z[1]; return x;
63 }
64
65 /* same as Flx_to_ZX, in place */
66 GEN
Flx_to_ZX_inplace(GEN z)67 Flx_to_ZX_inplace(GEN z)
68 {
69 long i, l = lg(z);
70 for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
71 settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
72 }
73
74 /*Flx_to_Flv=zx_to_zv*/
75 GEN
Flx_to_Flv(GEN x,long N)76 Flx_to_Flv(GEN x, long N)
77 {
78 long i, l;
79 GEN z = cgetg(N+1,t_VECSMALL);
80 if (typ(x) != t_VECSMALL) pari_err_TYPE("Flx_to_Flv",x);
81 l = lg(x)-1; x++;
82 for (i=1; i<l ; i++) z[i]=x[i];
83 for ( ; i<=N; i++) z[i]=0;
84 return z;
85 }
86
87 /*Flv_to_Flx=zv_to_zx*/
88 GEN
Flv_to_Flx(GEN x,long sv)89 Flv_to_Flx(GEN x, long sv)
90 {
91 long i, l=lg(x)+1;
92 GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
93 x--;
94 for (i=2; i<l ; i++) z[i]=x[i];
95 return Flx_renormalize(z,l);
96 }
97
98 /*Flm_to_FlxV=zm_to_zxV*/
99 GEN
Flm_to_FlxV(GEN x,long sv)100 Flm_to_FlxV(GEN x, long sv)
101 { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
102
103 /*FlxC_to_ZXC=zxC_to_ZXC*/
104 GEN
FlxC_to_ZXC(GEN x)105 FlxC_to_ZXC(GEN x)
106 { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
107
108 /*FlxC_to_ZXC=zxV_to_ZXV*/
109 GEN
FlxV_to_ZXV(GEN x)110 FlxV_to_ZXV(GEN x)
111 { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
112
113 void
FlxV_to_ZXV_inplace(GEN v)114 FlxV_to_ZXV_inplace(GEN v)
115 {
116 long i;
117 for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
118 }
119
120 /*FlxM_to_ZXM=zxM_to_ZXM*/
121 GEN
FlxM_to_ZXM(GEN x)122 FlxM_to_ZXM(GEN x)
123 { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
124
125 GEN
FlxV_to_FlxX(GEN x,long v)126 FlxV_to_FlxX(GEN x, long v)
127 {
128 long i, l = lg(x)+1;
129 GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
130 x--;
131 for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
132 return FlxX_renormalize(z,l);
133 }
134
135 GEN
FlxM_Flx_add_shallow(GEN x,GEN y,ulong p)136 FlxM_Flx_add_shallow(GEN x, GEN y, ulong p)
137 {
138 long l = lg(x), i, j;
139 GEN z = cgetg(l,t_MAT);
140
141 if (l==1) return z;
142 if (l != lgcols(x)) pari_err_OP( "+", x, y);
143 for (i=1; i<l; i++)
144 {
145 GEN zi = cgetg(l,t_COL), xi = gel(x,i);
146 gel(z,i) = zi;
147 for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
148 gel(zi,i) = Flx_add(gel(zi,i), y, p);
149 }
150 return z;
151 }
152
153 /***********************************************************************/
154 /** **/
155 /** Conversion to Flx **/
156 /** **/
157 /***********************************************************************/
158 /* Take an integer and return a scalar polynomial mod p, with evalvarn=vs */
159 GEN
Fl_to_Flx(ulong x,long sv)160 Fl_to_Flx(ulong x, long sv)
161 {
162 return x? mkvecsmall2(sv, x): pol0_Flx(sv);
163 }
164
165 /* a X^d */
166 GEN
monomial_Flx(ulong a,long d,long vs)167 monomial_Flx(ulong a, long d, long vs)
168 {
169 GEN P;
170 if (a==0) return pol0_Flx(vs);
171 P = const_vecsmall(d+2, 0);
172 P[1] = vs; P[d+2] = a;
173 return P;
174 }
175
176 GEN
Z_to_Flx(GEN x,ulong p,long sv)177 Z_to_Flx(GEN x, ulong p, long sv)
178 {
179 long u = umodiu(x,p);
180 return u? mkvecsmall2(sv, u): pol0_Flx(sv);
181 }
182
183 /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
184 GEN
ZX_to_Flx(GEN x,ulong p)185 ZX_to_Flx(GEN x, ulong p)
186 {
187 long i, lx = lg(x);
188 GEN a = cgetg(lx, t_VECSMALL);
189 a[1]=((ulong)x[1])&VARNBITS;
190 for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
191 return Flx_renormalize(a,lx);
192 }
193
194 /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
195 GEN
zx_to_Flx(GEN x,ulong p)196 zx_to_Flx(GEN x, ulong p)
197 {
198 long i, lx = lg(x);
199 GEN a = cgetg(lx, t_VECSMALL);
200 a[1] = x[1];
201 for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
202 return Flx_renormalize(a,lx);
203 }
204
205 ulong
Rg_to_Fl(GEN x,ulong p)206 Rg_to_Fl(GEN x, ulong p)
207 {
208 switch(typ(x))
209 {
210 case t_INT: return umodiu(x, p);
211 case t_FRAC: {
212 ulong z = umodiu(gel(x,1), p);
213 if (!z) return 0;
214 return Fl_div(z, umodiu(gel(x,2), p), p);
215 }
216 case t_PADIC: return padic_to_Fl(x, p);
217 case t_INTMOD: {
218 GEN q = gel(x,1), a = gel(x,2);
219 if (absequaliu(q, p)) return itou(a);
220 if (!dvdiu(q,p)) pari_err_MODULUS("Rg_to_Fl", q, utoi(p));
221 return umodiu(a, p);
222 }
223 default: pari_err_TYPE("Rg_to_Fl",x);
224 return 0; /* LCOV_EXCL_LINE */
225 }
226 }
227
228 ulong
Rg_to_F2(GEN x)229 Rg_to_F2(GEN x)
230 {
231 switch(typ(x))
232 {
233 case t_INT: return mpodd(x);
234 case t_FRAC:
235 if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
236 return mpodd(gel(x,1));
237 case t_PADIC:
238 if (!absequaliu(gel(x,2),2)) pari_err_OP("",x, mkintmodu(1,2));
239 if (valp(x) < 0) (void)Fl_inv(0,2);
240 return valp(x) & 1;
241 case t_INTMOD: {
242 GEN q = gel(x,1), a = gel(x,2);
243 if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
244 return mpodd(a);
245 }
246 default: pari_err_TYPE("Rg_to_F2",x);
247 return 0; /* LCOV_EXCL_LINE */
248 }
249 }
250
251 GEN
RgX_to_Flx(GEN x,ulong p)252 RgX_to_Flx(GEN x, ulong p)
253 {
254 long i, lx = lg(x);
255 GEN a = cgetg(lx, t_VECSMALL);
256 a[1]=((ulong)x[1])&VARNBITS;
257 for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
258 return Flx_renormalize(a,lx);
259 }
260
261 GEN
RgXV_to_FlxV(GEN x,ulong p)262 RgXV_to_FlxV(GEN x, ulong p)
263 { pari_APPLY_type(t_VEC, RgX_to_Flx(gel(x,i), p)) }
264
265 /* If x is a POLMOD, assume modulus is a multiple of T. */
266 GEN
Rg_to_Flxq(GEN x,GEN T,ulong p)267 Rg_to_Flxq(GEN x, GEN T, ulong p)
268 {
269 long ta, tx = typ(x), v = get_Flx_var(T);
270 GEN a, b;
271 if (is_const_t(tx))
272 {
273 if (tx == t_FFELT) return FF_to_Flxq(x);
274 return Fl_to_Flx(Rg_to_Fl(x, p), v);
275 }
276 switch(tx)
277 {
278 case t_POLMOD:
279 b = gel(x,1);
280 a = gel(x,2); ta = typ(a);
281 if (is_const_t(ta)) return Fl_to_Flx(Rg_to_Fl(a, p), v);
282 b = RgX_to_Flx(b, p); if (b[1] != v) break;
283 a = RgX_to_Flx(a, p); if (Flx_equal(b,T)) return a;
284 if (lgpol(Flx_rem(b,T,p))==0) return Flx_rem(a, T, p);
285 break;
286 case t_POL:
287 x = RgX_to_Flx(x,p);
288 if (x[1] != v) break;
289 return Flx_rem(x, T, p);
290 case t_RFRAC:
291 a = Rg_to_Flxq(gel(x,1), T,p);
292 b = Rg_to_Flxq(gel(x,2), T,p);
293 return Flxq_div(a,b, T,p);
294 }
295 pari_err_TYPE("Rg_to_Flxq",x);
296 return NULL; /* LCOV_EXCL_LINE */
297 }
298
299 /***********************************************************************/
300 /** **/
301 /** Basic operation on Flx **/
302 /** **/
303 /***********************************************************************/
304 /* = zx_renormalize. Similar to normalizepol, in place */
305 GEN
Flx_renormalize(GEN x,long lx)306 Flx_renormalize(GEN /*in place*/ x, long lx)
307 {
308 long i;
309 for (i = lx-1; i>1; i--)
310 if (x[i]) break;
311 stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
312 setlg(x, i+1); return x;
313 }
314
315 GEN
Flx_red(GEN z,ulong p)316 Flx_red(GEN z, ulong p)
317 {
318 long i, l = lg(z);
319 GEN x = cgetg(l, t_VECSMALL);
320 x[1] = z[1];
321 for (i=2; i<l; i++) x[i] = uel(z,i)%p;
322 return Flx_renormalize(x,l);
323 }
324
325 int
Flx_equal(GEN V,GEN W)326 Flx_equal(GEN V, GEN W)
327 {
328 long l = lg(V);
329 if (lg(W) != l) return 0;
330 while (--l > 1) /* do not compare variables, V[1] */
331 if (V[l] != W[l]) return 0;
332 return 1;
333 }
334
335 GEN
random_Flx(long d1,long vs,ulong p)336 random_Flx(long d1, long vs, ulong p)
337 {
338 long i, d = d1+2;
339 GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
340 for (i=2; i<d; i++) y[i] = random_Fl(p);
341 return Flx_renormalize(y,d);
342 }
343
344 static GEN
Flx_addspec(GEN x,GEN y,ulong p,long lx,long ly)345 Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
346 {
347 long i,lz;
348 GEN z;
349
350 if (ly>lx) swapspec(x,y, lx,ly);
351 lz = lx+2; z = cgetg(lz, t_VECSMALL);
352 for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
353 for ( ; i<lx; i++) z[i+2] = x[i];
354 z[1] = 0; return Flx_renormalize(z, lz);
355 }
356
357 GEN
Flx_add(GEN x,GEN y,ulong p)358 Flx_add(GEN x, GEN y, ulong p)
359 {
360 long i,lz;
361 GEN z;
362 long lx=lg(x);
363 long ly=lg(y);
364 if (ly>lx) swapspec(x,y, lx,ly);
365 lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
366 for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
367 for ( ; i<lx; i++) z[i] = x[i];
368 return Flx_renormalize(z, lz);
369 }
370
371 GEN
Flx_Fl_add(GEN y,ulong x,ulong p)372 Flx_Fl_add(GEN y, ulong x, ulong p)
373 {
374 GEN z;
375 long lz, i;
376 if (!lgpol(y))
377 return Fl_to_Flx(x,y[1]);
378 lz=lg(y);
379 z=cgetg(lz,t_VECSMALL);
380 z[1]=y[1];
381 z[2] = Fl_add(y[2],x,p);
382 for(i=3;i<lz;i++)
383 z[i] = y[i];
384 if (lz==3) z = Flx_renormalize(z,lz);
385 return z;
386 }
387
388 static GEN
Flx_subspec(GEN x,GEN y,ulong p,long lx,long ly)389 Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
390 {
391 long i,lz;
392 GEN z;
393
394 if (ly <= lx)
395 {
396 lz = lx+2; z = cgetg(lz, t_VECSMALL);
397 for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
398 for ( ; i<lx; i++) z[i+2] = x[i];
399 }
400 else
401 {
402 lz = ly+2; z = cgetg(lz, t_VECSMALL);
403 for (i=0; i<lx; i++) z[i+2] = Fl_sub(x[i],y[i],p);
404 for ( ; i<ly; i++) z[i+2] = Fl_neg(y[i],p);
405 }
406 z[1] = 0; return Flx_renormalize(z, lz);
407 }
408
409 GEN
Flx_sub(GEN x,GEN y,ulong p)410 Flx_sub(GEN x, GEN y, ulong p)
411 {
412 long i,lz,lx = lg(x), ly = lg(y);
413 GEN z;
414
415 if (ly <= lx)
416 {
417 lz = lx; z = cgetg(lz, t_VECSMALL);
418 for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
419 for ( ; i<lx; i++) z[i] = x[i];
420 }
421 else
422 {
423 lz = ly; z = cgetg(lz, t_VECSMALL);
424 for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
425 for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
426 }
427 z[1]=x[1]; return Flx_renormalize(z, lz);
428 }
429
430 GEN
Flx_Fl_sub(GEN y,ulong x,ulong p)431 Flx_Fl_sub(GEN y, ulong x, ulong p)
432 {
433 GEN z;
434 long lz = lg(y), i;
435 if (lz==2)
436 return Fl_to_Flx(Fl_neg(x, p),y[1]);
437 z = cgetg(lz, t_VECSMALL);
438 z[1] = y[1];
439 z[2] = Fl_sub(uel(y,2), x, p);
440 for(i=3; i<lz; i++)
441 z[i] = y[i];
442 if (lz==3) z = Flx_renormalize(z,lz);
443 return z;
444 }
445
446 static GEN
Flx_negspec(GEN x,ulong p,long l)447 Flx_negspec(GEN x, ulong p, long l)
448 {
449 long i;
450 GEN z = cgetg(l+2, t_VECSMALL) + 2;
451 for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
452 return z-2;
453 }
454
455 GEN
Flx_neg(GEN x,ulong p)456 Flx_neg(GEN x, ulong p)
457 {
458 GEN z = Flx_negspec(x+2, p, lgpol(x));
459 z[1] = x[1];
460 return z;
461 }
462
463 GEN
Flx_neg_inplace(GEN x,ulong p)464 Flx_neg_inplace(GEN x, ulong p)
465 {
466 long i, l = lg(x);
467 for (i=2; i<l; i++)
468 if (x[i]) x[i] = p - x[i];
469 return x;
470 }
471
472 GEN
Flx_double(GEN y,ulong p)473 Flx_double(GEN y, ulong p)
474 {
475 long i, l;
476 GEN z = cgetg_copy(y, &l); z[1] = y[1];
477 for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
478 return Flx_renormalize(z, l);
479 }
480 GEN
Flx_triple(GEN y,ulong p)481 Flx_triple(GEN y, ulong p)
482 {
483 long i, l;
484 GEN z = cgetg_copy(y, &l); z[1] = y[1];
485 for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
486 return Flx_renormalize(z, l);
487 }
488 GEN
Flx_Fl_mul(GEN y,ulong x,ulong p)489 Flx_Fl_mul(GEN y, ulong x, ulong p)
490 {
491 GEN z;
492 long i, l;
493 if (!x) return pol0_Flx(y[1]);
494 z = cgetg_copy(y, &l); z[1] = y[1];
495 if (HIGHWORD(x | p))
496 for(i=2; i<l; i++) z[i] = Fl_mul(y[i], x, p);
497 else
498 for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
499 return Flx_renormalize(z, l);
500 }
501 GEN
Flx_Fl_mul_to_monic(GEN y,ulong x,ulong p)502 Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
503 {
504 GEN z;
505 long i, l;
506 z = cgetg_copy(y, &l); z[1] = y[1];
507 if (HIGHWORD(x | p))
508 for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
509 else
510 for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
511 z[l-1] = 1; return z;
512 }
513
514 /* Return a*x^n if n>=0 and a\x^(-n) if n<0 */
515 GEN
Flx_shift(GEN a,long n)516 Flx_shift(GEN a, long n)
517 {
518 long i, l = lg(a);
519 GEN b;
520 if (l==2 || !n) return Flx_copy(a);
521 if (l+n<=2) return pol0_Flx(a[1]);
522 b = cgetg(l+n, t_VECSMALL);
523 b[1] = a[1];
524 if (n < 0)
525 for (i=2-n; i<l; i++) b[i+n] = a[i];
526 else
527 {
528 for (i=0; i<n; i++) b[2+i] = 0;
529 for (i=2; i<l; i++) b[i+n] = a[i];
530 }
531 return b;
532 }
533
534 GEN
Flx_normalize(GEN z,ulong p)535 Flx_normalize(GEN z, ulong p)
536 {
537 long l = lg(z)-1;
538 ulong p1 = z[l]; /* leading term */
539 if (p1 == 1) return z;
540 return Flx_Fl_mul_to_monic(z, Fl_inv(p1,p), p);
541 }
542
543 /* return (x * X^d) + y. Assume d > 0, shallow if x == 0*/
544 static GEN
Flx_addshift(GEN x,GEN y,ulong p,long d)545 Flx_addshift(GEN x, GEN y, ulong p, long d)
546 {
547 GEN xd,yd,zd = (GEN)avma;
548 long a,lz,ny = lgpol(y), nx = lgpol(x);
549 long vs = x[1];
550 if (nx == 0) return y;
551 x += 2; y += 2; a = ny-d;
552 if (a <= 0)
553 {
554 lz = (a>nx)? ny+2: nx+d+2;
555 (void)new_chunk(lz); xd = x+nx; yd = y+ny;
556 while (xd > x) *--zd = *--xd;
557 x = zd + a;
558 while (zd > x) *--zd = 0;
559 }
560 else
561 {
562 xd = new_chunk(d); yd = y+d;
563 x = Flx_addspec(x,yd,p, nx,a);
564 lz = (a>nx)? ny+2: lg(x)+d;
565 x += 2; while (xd > x) *--zd = *--xd;
566 }
567 while (yd > y) *--zd = *--yd;
568 *--zd = vs;
569 *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
570 }
571
572 /* shift polynomial + gerepile */
573 /* Do not set evalvarn*/
574 static GEN
Flx_shiftip(pari_sp av,GEN x,long v)575 Flx_shiftip(pari_sp av, GEN x, long v)
576 {
577 long i, lx = lg(x), ly;
578 GEN y;
579 if (!v || lx==2) return gerepileuptoleaf(av, x);
580 ly = lx + v; /* result length */
581 (void)new_chunk(ly); /* check that result fits */
582 x += lx; y = (GEN)av;
583 for (i = 2; i<lx; i++) *--y = *--x;
584 for (i = 0; i< v; i++) *--y = 0;
585 y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
586 return gc_const((pari_sp)y, y);
587 }
588
589 static long
get_Fl_threshold(ulong p,long mul,long mul2)590 get_Fl_threshold(ulong p, long mul, long mul2)
591 {
592 return SMALL_ULONG(p) ? mul: mul2;
593 }
594
595 #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
596 #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
597 #define LLQUARTWORD(x) ((x) & QUARTMASK)
598 #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
599 #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
600 #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
601 INLINE long
maxbitcoeffpol(ulong p,long n)602 maxbitcoeffpol(ulong p, long n)
603 {
604 GEN z = muliu(sqru(p - 1), n);
605 long b = expi(z) + 1;
606 /* only do expensive bit-packing if it saves at least 1 limb */
607 if (b <= BITS_IN_QUARTULONG)
608 {
609 if (nbits2nlong(n*b) == (n + 3)>>2)
610 b = BITS_IN_QUARTULONG;
611 }
612 else if (b <= BITS_IN_HALFULONG)
613 {
614 if (nbits2nlong(n*b) == (n + 1)>>1)
615 b = BITS_IN_HALFULONG;
616 }
617 else
618 {
619 long l = lgefint(z) - 2;
620 if (nbits2nlong(n*b) == n*l)
621 b = l*BITS_IN_LONG;
622 }
623 return b;
624 }
625
626 INLINE ulong
Flx_mullimb_ok(GEN x,GEN y,ulong p,long a,long b)627 Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
628 { /* Assume OK_ULONG*/
629 ulong p1 = 0;
630 long i;
631 for (i=a; i<b; i++)
632 if (y[i])
633 {
634 p1 += y[i] * x[-i];
635 if (p1 & HIGHBIT) p1 %= p;
636 }
637 return p1 % p;
638 }
639
640 INLINE ulong
Flx_mullimb(GEN x,GEN y,ulong p,ulong pi,long a,long b)641 Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
642 {
643 ulong p1 = 0;
644 long i;
645 for (i=a; i<b; i++)
646 if (y[i])
647 p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
648 return p1;
649 }
650
651 /* assume nx >= ny > 0 */
652 static GEN
Flx_mulspec_basecase(GEN x,GEN y,ulong p,long nx,long ny)653 Flx_mulspec_basecase(GEN x, GEN y, ulong p, long nx, long ny)
654 {
655 long i,lz,nz;
656 GEN z;
657
658 lz = nx+ny+1; nz = lz-2;
659 z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
660 if (SMALL_ULONG(p))
661 {
662 for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
663 for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
664 for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
665 }
666 else
667 {
668 ulong pi = get_Fl_red(p);
669 for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
670 for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
671 for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
672 }
673 z -= 2; return Flx_renormalize(z, lz);
674 }
675
676 static GEN
int_to_Flx(GEN z,ulong p)677 int_to_Flx(GEN z, ulong p)
678 {
679 long i, l = lgefint(z);
680 GEN x = cgetg(l, t_VECSMALL);
681 for (i=2; i<l; i++) x[i] = uel(z,i)%p;
682 return Flx_renormalize(x, l);
683 }
684
685 INLINE GEN
Flx_mulspec_mulii(GEN a,GEN b,ulong p,long na,long nb)686 Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
687 {
688 GEN z=muliispec(a,b,na,nb);
689 return int_to_Flx(z,p);
690 }
691
692 static GEN
Flx_to_int_halfspec(GEN a,long na)693 Flx_to_int_halfspec(GEN a, long na)
694 {
695 long j;
696 long n = (na+1)>>1UL;
697 GEN V = cgetipos(2+n);
698 GEN w;
699 for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
700 *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
701 if (j<na)
702 *w = a[j];
703 return V;
704 }
705
706 static GEN
int_to_Flx_half(GEN z,ulong p)707 int_to_Flx_half(GEN z, ulong p)
708 {
709 long i;
710 long lx = (lgefint(z)-2)*2+2;
711 GEN w, x = cgetg(lx, t_VECSMALL);
712 for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
713 {
714 x[i] = LOWWORD((ulong)*w)%p;
715 x[i+1] = HIGHWORD((ulong)*w)%p;
716 }
717 return Flx_renormalize(x, lx);
718 }
719
720 static GEN
Flx_mulspec_halfmulii(GEN a,GEN b,ulong p,long na,long nb)721 Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
722 {
723 GEN A = Flx_to_int_halfspec(a,na);
724 GEN B = Flx_to_int_halfspec(b,nb);
725 GEN z = mulii(A,B);
726 return int_to_Flx_half(z,p);
727 }
728
729 static GEN
Flx_to_int_quartspec(GEN a,long na)730 Flx_to_int_quartspec(GEN a, long na)
731 {
732 long j;
733 long n = (na+3)>>2UL;
734 GEN V = cgetipos(2+n);
735 GEN w;
736 for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
737 *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
738 switch (na-j)
739 {
740 case 3:
741 *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
742 break;
743 case 2:
744 *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
745 break;
746 case 1:
747 *w = a[j];
748 break;
749 case 0:
750 break;
751 }
752 return V;
753 }
754
755 static GEN
int_to_Flx_quart(GEN z,ulong p)756 int_to_Flx_quart(GEN z, ulong p)
757 {
758 long i;
759 long lx = (lgefint(z)-2)*4+2;
760 GEN w, x = cgetg(lx, t_VECSMALL);
761 for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
762 {
763 x[i] = LLQUARTWORD((ulong)*w)%p;
764 x[i+1] = HLQUARTWORD((ulong)*w)%p;
765 x[i+2] = LHQUARTWORD((ulong)*w)%p;
766 x[i+3] = HHQUARTWORD((ulong)*w)%p;
767 }
768 return Flx_renormalize(x, lx);
769 }
770
771 static GEN
Flx_mulspec_quartmulii(GEN a,GEN b,ulong p,long na,long nb)772 Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
773 {
774 GEN A = Flx_to_int_quartspec(a,na);
775 GEN B = Flx_to_int_quartspec(b,nb);
776 GEN z = mulii(A,B);
777 return int_to_Flx_quart(z,p);
778 }
779
780 /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
781 static GEN
Flx_eval2BILspec(GEN x,long k,long l)782 Flx_eval2BILspec(GEN x, long k, long l)
783 {
784 long i, lz = k*l, ki;
785 GEN pz = cgetipos(2+lz);
786 for (i=0; i < lz; i++)
787 *int_W(pz,i) = 0UL;
788 for (i=0, ki=0; i<l; i++, ki+=k)
789 *int_W(pz,ki) = x[i];
790 return int_normalize(pz,0);
791 }
792
793 static GEN
Z_mod2BIL_Flx_2(GEN x,long d,ulong p)794 Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
795 {
796 long i, offset, lm = lgefint(x)-2, l = d+3;
797 ulong pi = get_Fl_red(p);
798 GEN pol = cgetg(l, t_VECSMALL);
799 pol[1] = 0;
800 for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
801 pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
802 if (offset < lm)
803 pol[i+2] = (*int_W(x,offset)) % p;
804 return Flx_renormalize(pol,l);
805 }
806
807 static GEN
Z_mod2BIL_Flx_3(GEN x,long d,ulong p)808 Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
809 {
810 long i, offset, lm = lgefint(x)-2, l = d+3;
811 ulong pi = get_Fl_red(p);
812 GEN pol = cgetg(l, t_VECSMALL);
813 pol[1] = 0;
814 for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
815 pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
816 *int_W(x,offset), p, pi);
817 if (offset+1 < lm)
818 pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
819 else if (offset < lm)
820 pol[i+2] = (*int_W(x,offset)) % p;
821 return Flx_renormalize(pol,l);
822 }
823
824 static GEN
Z_mod2BIL_Flx(GEN x,long bs,long d,ulong p)825 Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
826 {
827 return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
828 }
829
830 static GEN
Flx_mulspec_mulii_inflate(GEN x,GEN y,long N,ulong p,long nx,long ny)831 Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
832 {
833 pari_sp av = avma;
834 GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
835 return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
836 }
837
838 static GEN
kron_pack_Flx_spec_bits(GEN x,long b,long l)839 kron_pack_Flx_spec_bits(GEN x, long b, long l) {
840 GEN y;
841 long i;
842 if (l == 0)
843 return gen_0;
844 y = cgetg(l + 1, t_VECSMALL);
845 for(i = 1; i <= l; i++)
846 y[i] = x[l - i];
847 return nv_fromdigits_2k(y, b);
848 }
849
850 /* assume b < BITS_IN_LONG */
851 static GEN
kron_unpack_Flx_bits_narrow(GEN z,long b,ulong p)852 kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
853 GEN v = binary_2k_nv(z, b), x;
854 long i, l = lg(v) + 1;
855 x = cgetg(l, t_VECSMALL);
856 for (i = 2; i < l; i++)
857 x[i] = v[l - i] % p;
858 return Flx_renormalize(x, l);
859 }
860
861 static GEN
kron_unpack_Flx_bits_wide(GEN z,long b,ulong p,ulong pi)862 kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
863 GEN v = binary_2k(z, b), x, y;
864 long i, l = lg(v) + 1, ly;
865 x = cgetg(l, t_VECSMALL);
866 for (i = 2; i < l; i++) {
867 y = gel(v, l - i);
868 ly = lgefint(y);
869 switch (ly) {
870 case 2: x[i] = 0; break;
871 case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
872 case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
873 case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
874 *int_W_lg(y, 0, ly), p, pi); break;
875 default: x[i] = umodiu(gel(v, l - i), p);
876 }
877 }
878 return Flx_renormalize(x, l);
879 }
880
881 static GEN
Flx_mulspec_Kronecker(GEN A,GEN B,long b,ulong p,long lA,long lB)882 Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
883 {
884 GEN C, D;
885 pari_sp av = avma;
886 A = kron_pack_Flx_spec_bits(A, b, lA);
887 B = kron_pack_Flx_spec_bits(B, b, lB);
888 C = gerepileuptoint(av, mulii(A, B));
889 if (b < BITS_IN_LONG)
890 D = kron_unpack_Flx_bits_narrow(C, b, p);
891 else
892 {
893 ulong pi = get_Fl_red(p);
894 D = kron_unpack_Flx_bits_wide(C, b, p, pi);
895 }
896 return D;
897 }
898
899 static GEN
Flx_sqrspec_Kronecker(GEN A,long b,ulong p,long lA)900 Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
901 {
902 GEN C, D;
903 A = kron_pack_Flx_spec_bits(A, b, lA);
904 C = sqri(A);
905 if (b < BITS_IN_LONG)
906 D = kron_unpack_Flx_bits_narrow(C, b, p);
907 else
908 {
909 ulong pi = get_Fl_red(p);
910 D = kron_unpack_Flx_bits_wide(C, b, p, pi);
911 }
912 return D;
913 }
914
915 /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
916 * b+2 were sent instead. na, nb = number of terms of a, b.
917 * Only c, c0, c1, c2 are genuine GEN.
918 */
919 static GEN
Flx_mulspec(GEN a,GEN b,ulong p,long na,long nb)920 Flx_mulspec(GEN a, GEN b, ulong p, long na, long nb)
921 {
922 GEN a0,c,c0;
923 long n0, n0a, i, v = 0;
924 pari_sp av;
925
926 while (na && !a[0]) { a++; na--; v++; }
927 while (nb && !b[0]) { b++; nb--; v++; }
928 if (na < nb) swapspec(a,b, na,nb);
929 if (!nb) return pol0_Flx(0);
930
931 av = avma;
932 if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
933 {
934 long m = maxbitcoeffpol(p,nb);
935 switch (m)
936 {
937 case BITS_IN_QUARTULONG:
938 return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
939 case BITS_IN_HALFULONG:
940 return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
941 case BITS_IN_LONG:
942 return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
943 case 2*BITS_IN_LONG:
944 return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
945 case 3*BITS_IN_LONG:
946 return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
947 default:
948 return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
949 }
950 }
951 if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
952 return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,na,nb), v);
953 i=(na>>1); n0=na-i; na=i;
954 a0=a+n0; n0a=n0;
955 while (n0a && !a[n0a-1]) n0a--;
956
957 if (nb > n0)
958 {
959 GEN b0,c1,c2;
960 long n0b;
961
962 nb -= n0; b0 = b+n0; n0b = n0;
963 while (n0b && !b[n0b-1]) n0b--;
964 c = Flx_mulspec(a,b,p,n0a,n0b);
965 c0 = Flx_mulspec(a0,b0,p,na,nb);
966
967 c2 = Flx_addspec(a0,a,p,na,n0a);
968 c1 = Flx_addspec(b0,b,p,nb,n0b);
969
970 c1 = Flx_mul(c1,c2,p);
971 c2 = Flx_add(c0,c,p);
972
973 c2 = Flx_neg_inplace(c2,p);
974 c2 = Flx_add(c1,c2,p);
975 c0 = Flx_addshift(c0,c2 ,p, n0);
976 }
977 else
978 {
979 c = Flx_mulspec(a,b,p,n0a,nb);
980 c0 = Flx_mulspec(a0,b,p,na,nb);
981 }
982 c0 = Flx_addshift(c0,c,p,n0);
983 return Flx_shiftip(av,c0, v);
984 }
985
986 GEN
Flx_mul(GEN x,GEN y,ulong p)987 Flx_mul(GEN x, GEN y, ulong p)
988 {
989 GEN z = Flx_mulspec(x+2,y+2,p, lgpol(x),lgpol(y));
990 z[1] = x[1]; return z;
991 }
992
993 static GEN
Flx_sqrspec_basecase(GEN x,ulong p,long nx)994 Flx_sqrspec_basecase(GEN x, ulong p, long nx)
995 {
996 long i, lz, nz;
997 ulong p1;
998 GEN z;
999
1000 if (!nx) return pol0_Flx(0);
1001 lz = (nx << 1) + 1, nz = lz-2;
1002 z = cgetg(lz, t_VECSMALL) + 2;
1003 if (SMALL_ULONG(p))
1004 {
1005 z[0] = x[0]*x[0]%p;
1006 for (i=1; i<nx; i++)
1007 {
1008 p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1009 p1 <<= 1;
1010 if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1011 z[i] = p1 % p;
1012 }
1013 for ( ; i<nz; i++)
1014 {
1015 p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1016 p1 <<= 1;
1017 if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1018 z[i] = p1 % p;
1019 }
1020 }
1021 else
1022 {
1023 ulong pi = get_Fl_red(p);
1024 z[0] = Fl_sqr_pre(x[0], p, pi);
1025 for (i=1; i<nx; i++)
1026 {
1027 p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1028 p1 = Fl_add(p1, p1, p);
1029 if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1030 z[i] = p1;
1031 }
1032 for ( ; i<nz; i++)
1033 {
1034 p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1035 p1 = Fl_add(p1, p1, p);
1036 if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1037 z[i] = p1;
1038 }
1039 }
1040 z -= 2; return Flx_renormalize(z, lz);
1041 }
1042
1043 static GEN
Flx_sqrspec_sqri(GEN a,ulong p,long na)1044 Flx_sqrspec_sqri(GEN a, ulong p, long na)
1045 {
1046 GEN z=sqrispec(a,na);
1047 return int_to_Flx(z,p);
1048 }
1049
1050 static GEN
Flx_sqrspec_halfsqri(GEN a,ulong p,long na)1051 Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
1052 {
1053 GEN z = sqri(Flx_to_int_halfspec(a,na));
1054 return int_to_Flx_half(z,p);
1055 }
1056
1057 static GEN
Flx_sqrspec_quartsqri(GEN a,ulong p,long na)1058 Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1059 {
1060 GEN z = sqri(Flx_to_int_quartspec(a,na));
1061 return int_to_Flx_quart(z,p);
1062 }
1063
1064 static GEN
Flx_sqrspec_sqri_inflate(GEN x,long N,ulong p,long nx)1065 Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1066 {
1067 pari_sp av = avma;
1068 GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1069 return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1070 }
1071
1072 static GEN
Flx_sqrspec(GEN a,ulong p,long na)1073 Flx_sqrspec(GEN a, ulong p, long na)
1074 {
1075 GEN a0, c, c0;
1076 long n0, n0a, i, v = 0, m;
1077 pari_sp av;
1078
1079 while (na && !a[0]) { a++; na--; v += 2; }
1080 if (!na) return pol0_Flx(0);
1081
1082 av = avma;
1083 if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1084 {
1085 m = maxbitcoeffpol(p,na);
1086 switch(m)
1087 {
1088 case BITS_IN_QUARTULONG:
1089 return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
1090 case BITS_IN_HALFULONG:
1091 return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
1092 case BITS_IN_LONG:
1093 return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1094 case 2*BITS_IN_LONG:
1095 return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
1096 case 3*BITS_IN_LONG:
1097 return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
1098 default:
1099 return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1100 }
1101 }
1102 if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1103 return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,na), v);
1104 i=(na>>1); n0=na-i; na=i;
1105 a0=a+n0; n0a=n0;
1106 while (n0a && !a[n0a-1]) n0a--;
1107
1108 c = Flx_sqrspec(a,p,n0a);
1109 c0= Flx_sqrspec(a0,p,na);
1110 if (p == 2) n0 *= 2;
1111 else
1112 {
1113 GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1114 t = Flx_sqr(t,p);
1115 c1= Flx_add(c0,c, p);
1116 c1= Flx_sub(t, c1, p);
1117 c0 = Flx_addshift(c0,c1,p,n0);
1118 }
1119 c0 = Flx_addshift(c0,c,p,n0);
1120 return Flx_shiftip(av,c0,v);
1121 }
1122
1123 GEN
Flx_sqr(GEN x,ulong p)1124 Flx_sqr(GEN x, ulong p)
1125 {
1126 GEN z = Flx_sqrspec(x+2,p, lgpol(x));
1127 z[1] = x[1]; return z;
1128 }
1129
1130 GEN
Flx_powu(GEN x,ulong n,ulong p)1131 Flx_powu(GEN x, ulong n, ulong p)
1132 {
1133 GEN y = pol1_Flx(x[1]), z;
1134 ulong m;
1135 if (n == 0) return y;
1136 m = n; z = x;
1137 for (;;)
1138 {
1139 if (m&1UL) y = Flx_mul(y,z, p);
1140 m >>= 1; if (!m) return y;
1141 z = Flx_sqr(z, p);
1142 }
1143 }
1144
1145 GEN
Flx_halve(GEN y,ulong p)1146 Flx_halve(GEN y, ulong p)
1147 {
1148 GEN z;
1149 long i, l;
1150 z = cgetg_copy(y, &l); z[1] = y[1];
1151 for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
1152 return z;
1153 }
1154
1155 static GEN
Flx_recipspec(GEN x,long l,long n)1156 Flx_recipspec(GEN x, long l, long n)
1157 {
1158 long i;
1159 GEN z=cgetg(n+2,t_VECSMALL)+2;
1160 for(i=0; i<l; i++)
1161 z[n-i-1] = x[i];
1162 for( ; i<n; i++)
1163 z[n-i-1] = 0;
1164 return Flx_renormalize(z-2,n+2);
1165 }
1166
1167 GEN
Flx_recip(GEN x)1168 Flx_recip(GEN x)
1169 {
1170 GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
1171 z[1]=x[1];
1172 return z;
1173 }
1174
1175 /* Return h^degpol(P) P(x / h) */
1176 GEN
Flx_rescale(GEN P,ulong h,ulong p)1177 Flx_rescale(GEN P, ulong h, ulong p)
1178 {
1179 long i, l = lg(P);
1180 GEN Q = cgetg(l,t_VECSMALL);
1181 ulong hi = h;
1182 Q[l-1] = P[l-1];
1183 for (i=l-2; i>=2; i--)
1184 {
1185 Q[i] = Fl_mul(P[i], hi, p);
1186 if (i == 2) break;
1187 hi = Fl_mul(hi,h, p);
1188 }
1189 Q[1] = P[1]; return Q;
1190 }
1191
1192 /*
1193 * x/polrecip(P)+O(x^n)
1194 */
1195 static GEN
Flx_invBarrett_basecase(GEN T,ulong p)1196 Flx_invBarrett_basecase(GEN T, ulong p)
1197 {
1198 long i, l=lg(T)-1, lr=l-1, k;
1199 GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1200 r[2] = 1;
1201 if (SMALL_ULONG(p))
1202 for (i=3;i<lr;i++)
1203 {
1204 ulong u = uel(T, l-i+2);
1205 for (k=3; k<i; k++)
1206 { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1207 r[i] = Fl_neg(u % p, p);
1208 }
1209 else
1210 for (i=3;i<lr;i++)
1211 {
1212 ulong u = Fl_neg(uel(T,l-i+2), p);
1213 for (k=3; k<i; k++)
1214 u = Fl_sub(u, Fl_mul(uel(T,l-i+k), uel(r,k), p), p);
1215 r[i] = u;
1216 }
1217 return Flx_renormalize(r,lr);
1218 }
1219
1220 /* Return new lgpol */
1221 static long
Flx_lgrenormalizespec(GEN x,long lx)1222 Flx_lgrenormalizespec(GEN x, long lx)
1223 {
1224 long i;
1225 for (i = lx-1; i>=0; i--)
1226 if (x[i]) break;
1227 return i+1;
1228 }
1229 static GEN
Flx_invBarrett_Newton(GEN T,ulong p)1230 Flx_invBarrett_Newton(GEN T, ulong p)
1231 {
1232 long nold, lx, lz, lq, l = degpol(T), lQ;
1233 GEN q, y, z, x = zero_zv(l+1) + 2;
1234 ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1235 pari_sp av;
1236
1237 y = T+2;
1238 q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1239 av = avma;
1240 /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1241
1242 /* initialize */
1243 x[0] = Fl_inv(q[0], p);
1244 if (lQ>1 && q[1])
1245 {
1246 ulong u = q[1];
1247 if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1248 x[1] = p - u; lx = 2;
1249 }
1250 else
1251 lx = 1;
1252 nold = 1;
1253 for (; mask > 1; set_avma(av))
1254 { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1255 long i, lnew, nnew = nold << 1;
1256
1257 if (mask & 1) nnew--;
1258 mask >>= 1;
1259
1260 lnew = nnew + 1;
1261 lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1262 z = Flx_mulspec(x, q, p, lx, lq); /* FIXME: high product */
1263 lz = lgpol(z); if (lz > lnew) lz = lnew;
1264 z += 2;
1265 /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1266 for (i = nold; i < lz; i++) if (z[i]) break;
1267 nold = nnew;
1268 if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1269
1270 /* z + i represents (x*q - 1) / t^i */
1271 lz = Flx_lgrenormalizespec (z+i, lz-i);
1272 z = Flx_mulspec(x, z+i, p, lx, lz); /* FIXME: low product */
1273 lz = lgpol(z); z += 2;
1274 if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1275
1276 lx = lz+ i;
1277 y = x + i; /* x -= z * t^i, in place */
1278 for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1279 }
1280 x -= 2; setlg(x, lx + 2); x[1] = T[1];
1281 return x;
1282 }
1283
1284 GEN
Flx_invBarrett(GEN T,ulong p)1285 Flx_invBarrett(GEN T, ulong p)
1286 {
1287 pari_sp ltop = avma;
1288 long l = lgpol(T);
1289 GEN r;
1290 if (l < 3) return pol0_Flx(T[1]);
1291 if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1292 {
1293 ulong c = T[l+1];
1294 if (c!=1)
1295 {
1296 ulong ci = Fl_inv(c,p);
1297 T=Flx_Fl_mul(T, ci, p);
1298 r=Flx_invBarrett_basecase(T,p);
1299 r=Flx_Fl_mul(r,ci,p);
1300 }
1301 else
1302 r=Flx_invBarrett_basecase(T,p);
1303 }
1304 else
1305 r = Flx_invBarrett_Newton(T,p);
1306 return gerepileuptoleaf(ltop, r);
1307 }
1308
1309 GEN
Flx_get_red(GEN T,ulong p)1310 Flx_get_red(GEN T, ulong p)
1311 {
1312 if (typ(T)!=t_VECSMALL
1313 || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1314 Flx_BARRETT2_LIMIT))
1315 return T;
1316 retmkvec2(Flx_invBarrett(T,p),T);
1317 }
1318
1319 /* separate from Flx_divrem for maximal speed. */
1320 static GEN
Flx_rem_basecase(GEN x,GEN y,ulong p)1321 Flx_rem_basecase(GEN x, GEN y, ulong p)
1322 {
1323 pari_sp av;
1324 GEN z, c;
1325 long dx,dy,dy1,dz,i,j;
1326 ulong p1,inv;
1327 long vs=x[1];
1328
1329 dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1330 dx = degpol(x);
1331 dz = dx-dy; if (dz < 0) return Flx_copy(x);
1332 x += 2; y += 2;
1333 inv = y[dy];
1334 if (inv != 1UL) inv = Fl_inv(inv,p);
1335 for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1336
1337 c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1338 z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1339
1340 if (SMALL_ULONG(p))
1341 {
1342 z[dz] = (inv*x[dx]) % p;
1343 for (i=dx-1; i>=dy; --i)
1344 {
1345 p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1346 for (j=i-dy1; j<=i && j<=dz; j++)
1347 {
1348 p1 += z[j]*y[i-j];
1349 if (p1 & HIGHBIT) p1 %= p;
1350 }
1351 p1 %= p;
1352 z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1353 }
1354 for (i=0; i<dy; i++)
1355 {
1356 p1 = z[0]*y[i];
1357 for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1358 {
1359 p1 += z[j]*y[i-j];
1360 if (p1 & HIGHBIT) p1 %= p;
1361 }
1362 c[i] = Fl_sub(x[i], p1%p, p);
1363 }
1364 }
1365 else
1366 {
1367 ulong pi = get_Fl_red(p);
1368 z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1369 for (i=dx-1; i>=dy; --i)
1370 {
1371 p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1372 for (j=i-dy1; j<=i && j<=dz; j++)
1373 p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1374 z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1375 }
1376 for (i=0; i<dy; i++)
1377 {
1378 p1 = Fl_mul_pre(z[0],y[i],p,pi);
1379 for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1380 p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1381 c[i] = Fl_sub(x[i], p1, p);
1382 }
1383 }
1384 i = dy-1; while (i>=0 && !c[i]) i--;
1385 set_avma(av); return Flx_renormalize(c-2, i+3);
1386 }
1387
1388 /* as FpX_divrem but working only on ulong types.
1389 * if relevant, *pr is the last object on stack */
1390 static GEN
Flx_divrem_basecase(GEN x,GEN y,ulong p,GEN * pr)1391 Flx_divrem_basecase(GEN x, GEN y, ulong p, GEN *pr)
1392 {
1393 GEN z,q,c;
1394 long dx,dy,dy1,dz,i,j;
1395 ulong p1,inv;
1396 long sv=x[1];
1397
1398 dy = degpol(y);
1399 if (dy<0) pari_err_INV("Flx_divrem",y);
1400 if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p);
1401 if (!dy)
1402 {
1403 if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1404 if (y[2] == 1UL) return Flx_copy(x);
1405 return Flx_Fl_mul(x, Fl_inv(y[2], p), p);
1406 }
1407 dx = degpol(x);
1408 dz = dx-dy;
1409 if (dz < 0)
1410 {
1411 q = pol0_Flx(sv);
1412 if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1413 return q;
1414 }
1415 x += 2;
1416 y += 2;
1417 z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1418 inv = uel(y, dy);
1419 if (inv != 1UL) inv = Fl_inv(inv,p);
1420 for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1421
1422 if (SMALL_ULONG(p))
1423 {
1424 z[dz] = (inv*x[dx]) % p;
1425 for (i=dx-1; i>=dy; --i)
1426 {
1427 p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1428 for (j=i-dy1; j<=i && j<=dz; j++)
1429 {
1430 p1 += z[j]*y[i-j];
1431 if (p1 & HIGHBIT) p1 %= p;
1432 }
1433 p1 %= p;
1434 z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1435 }
1436 }
1437 else
1438 {
1439 z[dz] = Fl_mul(inv, x[dx], p);
1440 for (i=dx-1; i>=dy; --i)
1441 { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1442 p1 = p - uel(x,i);
1443 for (j=i-dy1; j<=i && j<=dz; j++)
1444 p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1445 z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1446 }
1447 }
1448 q = Flx_renormalize(z-2, dz+3);
1449 if (!pr) return q;
1450
1451 c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1452 if (SMALL_ULONG(p))
1453 {
1454 for (i=0; i<dy; i++)
1455 {
1456 p1 = (ulong)z[0]*y[i];
1457 for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1458 {
1459 p1 += (ulong)z[j]*y[i-j];
1460 if (p1 & HIGHBIT) p1 %= p;
1461 }
1462 c[i] = Fl_sub(x[i], p1%p, p);
1463 }
1464 }
1465 else
1466 {
1467 for (i=0; i<dy; i++)
1468 {
1469 p1 = Fl_mul(z[0],y[i],p);
1470 for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1471 p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1472 c[i] = Fl_sub(x[i], p1, p);
1473 }
1474 }
1475 i=dy-1; while (i>=0 && !c[i]) i--;
1476 c = Flx_renormalize(c-2, i+3);
1477 if (pr == ONLY_DIVIDES)
1478 { if (lg(c) != 2) return NULL; }
1479 else
1480 *pr = c;
1481 return q;
1482 }
1483
1484 /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1485 * and mg is the Barrett inverse of T. */
1486 static GEN
Flx_divrem_Barrettspec(GEN x,long l,GEN mg,GEN T,ulong p,GEN * pr)1487 Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, GEN *pr)
1488 {
1489 GEN q, r;
1490 long lt = degpol(T); /*We discard the leading term*/
1491 long ld, lm, lT, lmg;
1492 ld = l-lt;
1493 lm = minss(ld, lgpol(mg));
1494 lT = Flx_lgrenormalizespec(T+2,lt);
1495 lmg = Flx_lgrenormalizespec(mg+2,lm);
1496 q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1497 q = Flx_mulspec(q+2,mg+2,p,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1498 q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1499 if (!pr) return q;
1500 r = Flx_mulspec(q+2,T+2,p,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1501 r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1502 if (pr == ONLY_REM) return r;
1503 *pr = r; return q;
1504 }
1505
1506 static GEN
Flx_divrem_Barrett(GEN x,GEN mg,GEN T,ulong p,GEN * pr)1507 Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, GEN *pr)
1508 {
1509 GEN q = NULL, r = Flx_copy(x);
1510 long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1511 long i;
1512 if (l <= lt)
1513 {
1514 if (pr == ONLY_REM) return Flx_copy(x);
1515 if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1516 if (pr) *pr = Flx_copy(x);
1517 return pol0_Flx(v);
1518 }
1519 if (lt <= 1)
1520 return Flx_divrem_basecase(x,T,p,pr);
1521 if (pr != ONLY_REM && l>lm)
1522 { q = zero_zv(l-lt+1); q[1] = T[1]; }
1523 while (l>lm)
1524 {
1525 GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,&zr);
1526 long lz = lgpol(zr);
1527 if (pr != ONLY_REM)
1528 {
1529 long lq = lgpol(zq);
1530 for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1531 }
1532 for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1533 l = l-lm+lz;
1534 }
1535 if (pr == ONLY_REM)
1536 {
1537 if (l > lt)
1538 r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,ONLY_REM);
1539 else
1540 r = Flx_renormalize(r, l+2);
1541 r[1] = v; return r;
1542 }
1543 if (l > lt)
1544 {
1545 GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p, pr ? &r: NULL);
1546 if (!q) q = zq;
1547 else
1548 {
1549 long lq = lgpol(zq);
1550 for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1551 }
1552 }
1553 else if (pr)
1554 r = Flx_renormalize(r, l+2);
1555 q[1] = v; q = Flx_renormalize(q, lg(q));
1556 if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1557 if (pr) { r[1] = v; *pr = r; }
1558 return q;
1559 }
1560
1561 GEN
Flx_divrem(GEN x,GEN T,ulong p,GEN * pr)1562 Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1563 {
1564 GEN B, y;
1565 long dy, dx, d;
1566 if (pr==ONLY_REM) return Flx_rem(x, T, p);
1567 y = get_Flx_red(T, &B);
1568 dy = degpol(y); dx = degpol(x); d = dx-dy;
1569 if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1570 return Flx_divrem_basecase(x,y,p,pr);
1571 else
1572 {
1573 pari_sp av = avma;
1574 GEN mg = B? B: Flx_invBarrett(y, p);
1575 GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pr);
1576 if (!q1) return gc_NULL(av);
1577 if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
1578 gerepileall(av,2,&q1,pr);
1579 return q1;
1580 }
1581 }
1582
1583 GEN
Flx_rem(GEN x,GEN T,ulong p)1584 Flx_rem(GEN x, GEN T, ulong p)
1585 {
1586 GEN B, y = get_Flx_red(T, &B);
1587 long dy = degpol(y), dx = degpol(x), d = dx-dy;
1588 if (d < 0) return Flx_copy(x);
1589 if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1590 return Flx_rem_basecase(x,y,p);
1591 else
1592 {
1593 pari_sp av=avma;
1594 GEN mg = B ? B: Flx_invBarrett(y, p);
1595 GEN r = Flx_divrem_Barrett(x, mg, y, p, ONLY_REM);
1596 return gerepileuptoleaf(av, r);
1597 }
1598 }
1599
1600 /* reduce T mod (X^n - 1, p). Shallow function */
1601 GEN
Flx_mod_Xnm1(GEN T,ulong n,ulong p)1602 Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1603 {
1604 long i, j, L = lg(T), l = n+2;
1605 GEN S;
1606 if (L <= l || n & ~LGBITS) return T;
1607 S = cgetg(l, t_VECSMALL);
1608 S[1] = T[1];
1609 for (i = 2; i < l; i++) S[i] = T[i];
1610 for (j = 2; i < L; i++) {
1611 S[j] = Fl_add(S[j], T[i], p);
1612 if (++j == l) j = 2;
1613 }
1614 return Flx_renormalize(S, l);
1615 }
1616 /* reduce T mod (X^n + 1, p). Shallow function */
1617 GEN
Flx_mod_Xn1(GEN T,ulong n,ulong p)1618 Flx_mod_Xn1(GEN T, ulong n, ulong p)
1619 {
1620 long i, j, L = lg(T), l = n+2;
1621 GEN S;
1622 if (L <= l || n & ~LGBITS) return T;
1623 S = cgetg(l, t_VECSMALL);
1624 S[1] = T[1];
1625 for (i = 2; i < l; i++) S[i] = T[i];
1626 for (j = 2; i < L; i++) {
1627 S[j] = Fl_sub(S[j], T[i], p);
1628 if (++j == l) j = 2;
1629 }
1630 return Flx_renormalize(S, l);
1631 }
1632
1633 struct _Flxq {
1634 GEN aut;
1635 GEN T;
1636 ulong p;
1637 };
1638
1639 static GEN
_Flx_divrem(void * E,GEN x,GEN y,GEN * r)1640 _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1641 {
1642 struct _Flxq *D = (struct _Flxq*) E;
1643 return Flx_divrem(x, y, D->p, r);
1644 }
1645 static GEN
_Flx_add(void * E,GEN x,GEN y)1646 _Flx_add(void * E, GEN x, GEN y) {
1647 struct _Flxq *D = (struct _Flxq*) E;
1648 return Flx_add(x, y, D->p);
1649 }
1650 static GEN
_Flx_mul(void * E,GEN x,GEN y)1651 _Flx_mul(void *E, GEN x, GEN y) {
1652 struct _Flxq *D = (struct _Flxq*) E;
1653 return Flx_mul(x, y, D->p);
1654 }
1655 static GEN
_Flx_sqr(void * E,GEN x)1656 _Flx_sqr(void *E, GEN x) {
1657 struct _Flxq *D = (struct _Flxq*) E;
1658 return Flx_sqr(x, D->p);
1659 }
1660
1661 static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1662
1663 GEN
Flx_digits(GEN x,GEN T,ulong p)1664 Flx_digits(GEN x, GEN T, ulong p)
1665 {
1666 pari_sp av = avma;
1667 struct _Flxq D;
1668 long d = degpol(T), n = (lgpol(x)+d-1)/d;
1669 GEN z;
1670 D.p = p;
1671 z = gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1672 return gerepileupto(av, z);
1673 }
1674
1675 GEN
FlxV_Flx_fromdigits(GEN x,GEN T,ulong p)1676 FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1677 {
1678 pari_sp av = avma;
1679 struct _Flxq D;
1680 GEN z;
1681 D.p = p;
1682 z = gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1683 return gerepileupto(av, z);
1684 }
1685
1686 long
Flx_val(GEN x)1687 Flx_val(GEN x)
1688 {
1689 long i, l=lg(x);
1690 if (l==2) return LONG_MAX;
1691 for (i=2; i<l && x[i]==0; i++) /*empty*/;
1692 return i-2;
1693 }
1694 long
Flx_valrem(GEN x,GEN * Z)1695 Flx_valrem(GEN x, GEN *Z)
1696 {
1697 long v, i, l=lg(x);
1698 GEN y;
1699 if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1700 for (i=2; i<l && x[i]==0; i++) /*empty*/;
1701 v = i-2;
1702 if (v == 0) { *Z = x; return 0; }
1703 l -= v;
1704 y = cgetg(l, t_VECSMALL); y[1] = x[1];
1705 for (i=2; i<l; i++) y[i] = x[i+v];
1706 *Z = y; return v;
1707 }
1708
1709 GEN
Flx_deriv(GEN z,ulong p)1710 Flx_deriv(GEN z, ulong p)
1711 {
1712 long i,l = lg(z)-1;
1713 GEN x;
1714 if (l < 2) l = 2;
1715 x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1716 if (HIGHWORD(l | p))
1717 for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1718 else
1719 for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1720 return Flx_renormalize(x,l);
1721 }
1722
1723 static GEN
Flx_integXn(GEN x,long n,ulong p)1724 Flx_integXn(GEN x, long n, ulong p)
1725 {
1726 long i, lx = lg(x);
1727 GEN y;
1728 if (lx == 2) return Flx_copy(x);
1729 y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1730 for (i=2; i<lx; i++)
1731 {
1732 ulong xi = uel(x,i);
1733 if (xi == 0)
1734 uel(y,i) = 0;
1735 else
1736 {
1737 ulong j = n+i-1;
1738 ulong d = ugcd(j, xi);
1739 if (d==1)
1740 uel(y,i) = Fl_div(xi, j, p);
1741 else
1742 uel(y,i) = Fl_div(xi/d, j/d, p);
1743 }
1744 }
1745 return Flx_renormalize(y, lx);;
1746 }
1747
1748 GEN
Flx_integ(GEN x,ulong p)1749 Flx_integ(GEN x, ulong p)
1750 {
1751 long i, lx = lg(x);
1752 GEN y;
1753 if (lx == 2) return Flx_copy(x);
1754 y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1755 uel(y,2) = 0;
1756 for (i=3; i<=lx; i++)
1757 uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1758 return Flx_renormalize(y, lx+1);;
1759 }
1760
1761 /* assume p prime */
1762 GEN
Flx_diff1(GEN P,ulong p)1763 Flx_diff1(GEN P, ulong p)
1764 {
1765 return Flx_sub(Flx_translate1(P, p), P, p);
1766 }
1767
1768 GEN
Flx_deflate(GEN x0,long d)1769 Flx_deflate(GEN x0, long d)
1770 {
1771 GEN z, y, x;
1772 long i,id, dy, dx = degpol(x0);
1773 if (d == 1 || dx <= 0) return Flx_copy(x0);
1774 dy = dx/d;
1775 y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1776 z = y + 2;
1777 x = x0+ 2;
1778 for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1779 return y;
1780 }
1781
1782 GEN
Flx_inflate(GEN x0,long d)1783 Flx_inflate(GEN x0, long d)
1784 {
1785 long i, id, dy, dx = degpol(x0);
1786 GEN x = x0 + 2, z, y;
1787 if (dx <= 0) return Flx_copy(x0);
1788 dy = dx*d;
1789 y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1790 z = y + 2;
1791 for (i=0; i<=dy; i++) z[i] = 0;
1792 for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1793 return y;
1794 }
1795
1796 /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1797 GEN
Flx_splitting(GEN p,long k)1798 Flx_splitting(GEN p, long k)
1799 {
1800 long n = degpol(p), v = p[1], m, i, j, l;
1801 GEN r;
1802
1803 m = n/k;
1804 r = cgetg(k+1,t_VEC);
1805 for(i=1; i<=k; i++)
1806 {
1807 gel(r,i) = cgetg(m+3, t_VECSMALL);
1808 mael(r,i,1) = v;
1809 }
1810 for (j=1, i=0, l=2; i<=n; i++)
1811 {
1812 mael(r,j,l) = p[2+i];
1813 if (j==k) { j=1; l++; } else j++;
1814 }
1815 for(i=1; i<=k; i++)
1816 gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1817 return r;
1818 }
1819 static GEN
Flx_halfgcd_basecase(GEN a,GEN b,ulong p)1820 Flx_halfgcd_basecase(GEN a, GEN b, ulong p)
1821 {
1822 pari_sp av=avma;
1823 GEN u,u1,v,v1;
1824 long vx = a[1];
1825 long n = lgpol(a)>>1;
1826 u1 = v = pol0_Flx(vx);
1827 u = v1 = pol1_Flx(vx);
1828 while (lgpol(b)>n)
1829 {
1830 GEN r, q = Flx_divrem(a,b,p, &r);
1831 a = b; b = r; swap(u,u1); swap(v,v1);
1832 u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
1833 v1 = Flx_sub(v1, Flx_mul(v, q ,p), p);
1834 if (gc_needed(av,2))
1835 {
1836 if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
1837 gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
1838 }
1839 }
1840 return gerepilecopy(av, mkmat2(mkcol2(u,u1), mkcol2(v,v1)));
1841 }
1842 /* ux + vy */
1843 static GEN
Flx_addmulmul(GEN u,GEN v,GEN x,GEN y,ulong p)1844 Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p)
1845 { return Flx_add(Flx_mul(u,x, p), Flx_mul(v,y, p), p); }
1846
1847 static GEN
FlxM_Flx_mul2(GEN M,GEN x,GEN y,ulong p)1848 FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p)
1849 {
1850 GEN res = cgetg(3, t_COL);
1851 gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p);
1852 gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p);
1853 return res;
1854 }
1855
1856 #if 0
1857 static GEN
1858 FlxM_mul2_old(GEN M, GEN N, ulong p)
1859 {
1860 GEN res = cgetg(3, t_MAT);
1861 gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1862 gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1863 return res;
1864 }
1865 #endif
1866 /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1867 static GEN
FlxM_mul2(GEN A,GEN B,ulong p)1868 FlxM_mul2(GEN A, GEN B, ulong p)
1869 {
1870 GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1871 GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1872 GEN M1 = Flx_mul(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p);
1873 GEN M2 = Flx_mul(Flx_add(A21,A22, p), B11, p);
1874 GEN M3 = Flx_mul(A11, Flx_sub(B12,B22, p), p);
1875 GEN M4 = Flx_mul(A22, Flx_sub(B21,B11, p), p);
1876 GEN M5 = Flx_mul(Flx_add(A11,A12, p), B22, p);
1877 GEN M6 = Flx_mul(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p);
1878 GEN M7 = Flx_mul(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p);
1879 GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1880 GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1881 retmkmat2(mkcol2(Flx_add(T1,T2, p), Flx_add(M2,M4, p)),
1882 mkcol2(Flx_add(M3,M5, p), Flx_add(T3,T4, p)));
1883 }
1884
1885 /* Return [0,1;1,-q]*M */
1886 static GEN
Flx_FlxM_qmul(GEN q,GEN M,ulong p)1887 Flx_FlxM_qmul(GEN q, GEN M, ulong p)
1888 {
1889 GEN u, v, res = cgetg(3, t_MAT);
1890 u = Flx_sub(gcoeff(M,1,1), Flx_mul(gcoeff(M,2,1), q, p), p);
1891 gel(res,1) = mkcol2(gcoeff(M,2,1), u);
1892 v = Flx_sub(gcoeff(M,1,2), Flx_mul(gcoeff(M,2,2), q, p), p);
1893 gel(res,2) = mkcol2(gcoeff(M,2,2), v);
1894 return res;
1895 }
1896
1897 static GEN
matid2_FlxM(long v)1898 matid2_FlxM(long v)
1899 {
1900 return mkmat2(mkcol2(pol1_Flx(v),pol0_Flx(v)),
1901 mkcol2(pol0_Flx(v),pol1_Flx(v)));
1902 }
1903
1904 static GEN
Flx_halfgcd_split(GEN x,GEN y,ulong p)1905 Flx_halfgcd_split(GEN x, GEN y, ulong p)
1906 {
1907 pari_sp av=avma;
1908 GEN R, S, V;
1909 GEN y1, r, q;
1910 long l = lgpol(x), n = l>>1, k;
1911 if (lgpol(y)<=n) return matid2_FlxM(x[1]);
1912 R = Flx_halfgcd(Flx_shift(x,-n),Flx_shift(y,-n),p);
1913 V = FlxM_Flx_mul2(R,x,y,p); y1 = gel(V,2);
1914 if (lgpol(y1)<=n) return gerepilecopy(av, R);
1915 q = Flx_divrem(gel(V,1), y1, p, &r);
1916 k = 2*n-degpol(y1);
1917 S = Flx_halfgcd(Flx_shift(y1,-k), Flx_shift(r,-k),p);
1918 return gerepileupto(av, FlxM_mul2(S,Flx_FlxM_qmul(q,R,p),p));
1919 }
1920
1921 /* Return M in GL_2(Fl[X]) such that:
1922 if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
1923 */
1924
1925 static GEN
Flx_halfgcd_i(GEN x,GEN y,ulong p)1926 Flx_halfgcd_i(GEN x, GEN y, ulong p)
1927 {
1928 if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
1929 return Flx_halfgcd_basecase(x,y,p);
1930 return Flx_halfgcd_split(x,y,p);
1931 }
1932
1933 GEN
Flx_halfgcd(GEN x,GEN y,ulong p)1934 Flx_halfgcd(GEN x, GEN y, ulong p)
1935 {
1936 pari_sp av;
1937 GEN M,q,r;
1938 long lx=lgpol(x), ly=lgpol(y);
1939 if (!lx)
1940 {
1941 long v = x[1];
1942 retmkmat2(mkcol2(pol0_Flx(v),pol1_Flx(v)),
1943 mkcol2(pol1_Flx(v),pol0_Flx(v)));
1944 }
1945 if (ly < lx) return Flx_halfgcd_i(x,y,p);
1946 av = avma;
1947 q = Flx_divrem(y,x,p,&r);
1948 M = Flx_halfgcd_i(x,r,p);
1949 gcoeff(M,1,1) = Flx_sub(gcoeff(M,1,1), Flx_mul(q, gcoeff(M,1,2), p), p);
1950 gcoeff(M,2,1) = Flx_sub(gcoeff(M,2,1), Flx_mul(q, gcoeff(M,2,2), p), p);
1951 return gerepilecopy(av, M);
1952 }
1953
1954 /*Do not garbage collect*/
1955 static GEN
Flx_gcd_basecase(GEN a,GEN b,ulong p)1956 Flx_gcd_basecase(GEN a, GEN b, ulong p)
1957 {
1958 pari_sp av = avma;
1959 ulong iter = 0;
1960 if (lg(b) > lg(a)) swap(a, b);
1961 while (lgpol(b))
1962 {
1963 GEN c = Flx_rem(a,b,p);
1964 iter++; a = b; b = c;
1965 if (gc_needed(av,2))
1966 {
1967 if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
1968 gerepileall(av,2, &a,&b);
1969 }
1970 }
1971 return iter < 2 ? Flx_copy(a) : a;
1972 }
1973
1974 GEN
Flx_gcd(GEN x,GEN y,ulong p)1975 Flx_gcd(GEN x, GEN y, ulong p)
1976 {
1977 pari_sp av = avma;
1978 long lim;
1979 if (!lgpol(x)) return Flx_copy(y);
1980 lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
1981 while (lgpol(y) >= lim)
1982 {
1983 GEN c;
1984 if (lgpol(y)<=(lgpol(x)>>1))
1985 {
1986 GEN r = Flx_rem(x, y, p);
1987 x = y; y = r;
1988 }
1989 c = FlxM_Flx_mul2(Flx_halfgcd(x,y, p), x, y, p);
1990 x = gel(c,1); y = gel(c,2);
1991 if (gc_needed(av,2))
1992 {
1993 if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
1994 gerepileall(av,2,&x,&y);
1995 }
1996 }
1997 return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p));
1998 }
1999
2000 int
Flx_is_squarefree(GEN z,ulong p)2001 Flx_is_squarefree(GEN z, ulong p)
2002 {
2003 pari_sp av = avma;
2004 GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
2005 return gc_bool(av, degpol(d) == 0);
2006 }
2007
2008 static long
Flx_is_smooth_squarefree(GEN f,long r,ulong p)2009 Flx_is_smooth_squarefree(GEN f, long r, ulong p)
2010 {
2011 pari_sp av = avma;
2012 long i;
2013 GEN sx = polx_Flx(f[1]), a = sx;
2014 for(i=1;;i++)
2015 {
2016 if (degpol(f)<=r) return gc_long(av,1);
2017 a = Flxq_powu(Flx_rem(a,f,p), p, f, p);
2018 if (Flx_equal(a, sx)) return gc_long(av,1);
2019 if (i==r) return gc_long(av,0);
2020 f = Flx_div(f, Flx_gcd(Flx_sub(a,sx,p),f,p),p);
2021 }
2022 }
2023
2024 static long
Flx_is_l_pow(GEN x,ulong p)2025 Flx_is_l_pow(GEN x, ulong p)
2026 {
2027 ulong i, lx = lgpol(x);
2028 for (i=1; i<lx; i++)
2029 if (x[i+2] && i%p) return 0;
2030 return 1;
2031 }
2032
2033 int
Flx_is_smooth(GEN g,long r,ulong p)2034 Flx_is_smooth(GEN g, long r, ulong p)
2035 {
2036 while (1)
2037 {
2038 GEN f = Flx_gcd(g, Flx_deriv(g, p), p);
2039 if (!Flx_is_smooth_squarefree(Flx_div(g, f, p), r, p))
2040 return 0;
2041 if (degpol(f)==0) return 1;
2042 g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2043 }
2044 }
2045
2046 static GEN
Flx_extgcd_basecase(GEN a,GEN b,ulong p,GEN * ptu,GEN * ptv)2047 Flx_extgcd_basecase(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv)
2048 {
2049 pari_sp av=avma;
2050 GEN u,v,d,d1,v1;
2051 long vx = a[1];
2052 d = a; d1 = b;
2053 v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2054 while (lgpol(d1))
2055 {
2056 GEN r, q = Flx_divrem(d,d1,p, &r);
2057 v = Flx_sub(v,Flx_mul(q,v1,p),p);
2058 u=v; v=v1; v1=u;
2059 u=r; d=d1; d1=u;
2060 if (gc_needed(av,2))
2061 {
2062 if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(d));
2063 gerepileall(av,5, &d,&d1,&u,&v,&v1);
2064 }
2065 }
2066 if (ptu) *ptu = Flx_div(Flx_sub(d, Flx_mul(b,v,p), p), a, p);
2067 *ptv = v; return d;
2068 }
2069
2070 static GEN
Flx_extgcd_halfgcd(GEN x,GEN y,ulong p,GEN * ptu,GEN * ptv)2071 Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2072 {
2073 pari_sp av=avma;
2074 GEN u,v,R = matid2_FlxM(x[1]);
2075 long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2076 while (lgpol(y) >= lim)
2077 {
2078 GEN M, c;
2079 if (lgpol(y)<=(lgpol(x)>>1))
2080 {
2081 GEN r, q = Flx_divrem(x, y, p, &r);
2082 x = y; y = r;
2083 R = Flx_FlxM_qmul(q, R, p);
2084 }
2085 M = Flx_halfgcd(x,y, p);
2086 c = FlxM_Flx_mul2(M, x,y, p);
2087 R = FlxM_mul2(M, R, p);
2088 x = gel(c,1); y = gel(c,2);
2089 gerepileall(av,3,&x,&y,&R);
2090 }
2091 y = Flx_extgcd_basecase(x,y,p,&u,&v);
2092 if (ptu) *ptu = Flx_addmulmul(u,v,gcoeff(R,1,1),gcoeff(R,2,1),p);
2093 *ptv = Flx_addmulmul(u,v,gcoeff(R,1,2),gcoeff(R,2,2),p);
2094 return y;
2095 }
2096
2097 /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2098 * ux + vy = gcd (mod p) */
2099 GEN
Flx_extgcd(GEN x,GEN y,ulong p,GEN * ptu,GEN * ptv)2100 Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2101 {
2102 GEN d;
2103 pari_sp ltop=avma;
2104 long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2105 if (lgpol(y) >= lim)
2106 d = Flx_extgcd_halfgcd(x, y, p, ptu, ptv);
2107 else
2108 d = Flx_extgcd_basecase(x, y, p, ptu, ptv);
2109 gerepileall(ltop,ptu?3:2,&d,ptv,ptu);
2110 return d;
2111 }
2112
2113 ulong
Flx_resultant(GEN a,GEN b,ulong p)2114 Flx_resultant(GEN a, GEN b, ulong p)
2115 {
2116 long da,db,dc,cnt;
2117 ulong lb, res = 1UL;
2118 pari_sp av;
2119 GEN c;
2120
2121 if (lgpol(a)==0 || lgpol(b)==0) return 0;
2122 da = degpol(a);
2123 db = degpol(b);
2124 if (db > da)
2125 {
2126 swapspec(a,b, da,db);
2127 if (both_odd(da,db)) res = p-res;
2128 }
2129 else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2130 cnt = 0; av = avma;
2131 while (db)
2132 {
2133 lb = b[db+2];
2134 c = Flx_rem(a,b, p);
2135 a = b; b = c; dc = degpol(c);
2136 if (dc < 0) return gc_long(av,0);
2137
2138 if (both_odd(da,db)) res = p - res;
2139 if (lb != 1) res = Fl_mul(res, Fl_powu(lb, da - dc, p), p);
2140 if (++cnt == 100) { cnt = 0; gerepileall(av, 2, &a, &b); }
2141 da = db; /* = degpol(a) */
2142 db = dc; /* = degpol(b) */
2143 }
2144 return gc_ulong(av, Fl_mul(res, Fl_powu(b[2], da, p), p));
2145 }
2146
2147 /* If resultant is 0, *ptU and *ptV are not set */
2148 ulong
Flx_extresultant(GEN a,GEN b,ulong p,GEN * ptU,GEN * ptV)2149 Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2150 {
2151 GEN z,q,u,v, x = a, y = b;
2152 ulong lb, res = 1UL;
2153 pari_sp av = avma;
2154 long dx, dy, dz;
2155 long vs=a[1];
2156
2157 dx = degpol(x);
2158 dy = degpol(y);
2159 if (dy > dx)
2160 {
2161 swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
2162 a = x; b = y;
2163 if (both_odd(dx,dy)) res = p-res;
2164 }
2165 /* dy <= dx */
2166 if (dy < 0) return 0;
2167
2168 u = pol0_Flx(vs);
2169 v = pol1_Flx(vs); /* v = 1 */
2170 while (dy)
2171 { /* b u = x (a), b v = y (a) */
2172 lb = y[dy+2];
2173 q = Flx_divrem(x,y, p, &z);
2174 x = y; y = z; /* (x,y) = (y, x - q y) */
2175 dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2176 z = Flx_sub(u, Flx_mul(q,v, p), p);
2177 u = v; v = z; /* (u,v) = (v, u - q v) */
2178
2179 if (both_odd(dx,dy)) res = p - res;
2180 if (lb != 1) res = Fl_mul(res, Fl_powu(lb, dx-dz, p), p);
2181 dx = dy; /* = degpol(x) */
2182 dy = dz; /* = degpol(y) */
2183 }
2184 res = Fl_mul(res, Fl_powu(y[2], dx, p), p);
2185 lb = Fl_mul(res, Fl_inv(y[2],p), p);
2186 v = gerepileuptoleaf(av, Flx_Fl_mul(v, lb, p));
2187 av = avma;
2188 u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul(b,v,p), p);
2189 u = gerepileuptoleaf(av, Flx_div(u,a,p)); /* = (res - b v) / a */
2190 *ptU = u;
2191 *ptV = v; return res;
2192 }
2193
2194 ulong
Flx_eval_powers_pre(GEN x,GEN y,ulong p,ulong pi)2195 Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2196 {
2197 ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2198 LOCAL_OVERFLOW;
2199 LOCAL_HIREMAINDER;
2200 x++;
2201
2202 if (lx == 1)
2203 return 0;
2204 l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2205 while (++i < lx) {
2206 l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2207 l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2208 }
2209 if (v1 == 0) return remll_pre(h1, l1, p, pi);
2210 else return remlll_pre(v1, h1, l1, p, pi);
2211 }
2212
2213 INLINE ulong
Flx_eval_pre_i(GEN x,ulong y,ulong p,ulong pi)2214 Flx_eval_pre_i(GEN x, ulong y, ulong p, ulong pi)
2215 {
2216 ulong p1;
2217 long i=lg(x)-1;
2218 if (i<=2)
2219 return (i==2)? x[2]: 0;
2220 p1 = x[i];
2221 for (i--; i>=2; i--)
2222 p1 = Fl_addmul_pre(uel(x, i), p1, y, p, pi);
2223 return p1;
2224 }
2225
2226 ulong
Flx_eval_pre(GEN x,ulong y,ulong p,ulong pi)2227 Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2228 {
2229 if (degpol(x) > 15)
2230 {
2231 pari_sp av = avma;
2232 GEN v = Fl_powers_pre(y, degpol(x), p, pi);
2233 ulong r = Flx_eval_powers_pre(x, v, p, pi);
2234 return gc_ulong(av,r);
2235 }
2236 else
2237 return Flx_eval_pre_i(x, y, p, pi);
2238 }
2239
2240 ulong
Flx_eval(GEN x,ulong y,ulong p)2241 Flx_eval(GEN x, ulong y, ulong p)
2242 {
2243 return Flx_eval_pre(x, y, p, get_Fl_red(p));
2244 }
2245
2246 ulong
Flv_prod_pre(GEN x,ulong p,ulong pi)2247 Flv_prod_pre(GEN x, ulong p, ulong pi)
2248 {
2249 pari_sp ltop = avma;
2250 GEN v;
2251 long i,k,lx = lg(x);
2252 if (lx == 1) return 1UL;
2253 if (lx == 2) return uel(x,1);
2254 v = cgetg(1+(lx << 1), t_VECSMALL);
2255 k = 1;
2256 for (i=1; i<lx-1; i+=2)
2257 uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2258 if (i < lx) uel(v,k++) = uel(x,i);
2259 while (k > 2)
2260 {
2261 lx = k; k = 1;
2262 for (i=1; i<lx-1; i+=2)
2263 uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2264 if (i < lx) uel(v,k++) = uel(v,i);
2265 }
2266 return gc_ulong(ltop, uel(v,1));
2267 }
2268
2269 ulong
Flv_prod(GEN v,ulong p)2270 Flv_prod(GEN v, ulong p)
2271 {
2272 return Flv_prod_pre(v, p, get_Fl_red(p));
2273 }
2274
2275 GEN
FlxV_prod(GEN V,ulong p)2276 FlxV_prod(GEN V, ulong p)
2277 {
2278 struct _Flxq D;
2279 D.T = NULL; D.aut = NULL; D.p = p;
2280 return gen_product(V, (void *)&D, &_Flx_mul);
2281 }
2282
2283 /* compute prod (x - a[i]) */
2284 GEN
Flv_roots_to_pol(GEN a,ulong p,long vs)2285 Flv_roots_to_pol(GEN a, ulong p, long vs)
2286 {
2287 struct _Flxq D;
2288 long i,k,lx = lg(a);
2289 GEN p1;
2290 if (lx == 1) return pol1_Flx(vs);
2291 p1 = cgetg(lx, t_VEC);
2292 for (k=1,i=1; i<lx-1; i+=2)
2293 gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2294 Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2295 if (i < lx)
2296 gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2297 D.T = NULL; D.aut = NULL; D.p = p;
2298 setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2299 }
2300
2301 /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2302 INLINE void
Flv_inv_pre_indir(GEN w,GEN v,ulong p,ulong pi)2303 Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2304 {
2305 pari_sp av = avma;
2306 long n = lg(w), i;
2307 ulong u;
2308 GEN c;
2309
2310 if (n == 1) return;
2311 c = cgetg(n, t_VECSMALL); c[1] = w[1];
2312 for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2313 i = n-1; u = Fl_inv(c[i], p);
2314 for ( ; i > 1; --i)
2315 {
2316 ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2317 u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2318 }
2319 v[1] = u; set_avma(av);
2320 }
2321
2322 void
Flv_inv_pre_inplace(GEN v,ulong p,ulong pi)2323 Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2324
2325 GEN
Flv_inv_pre(GEN w,ulong p,ulong pi)2326 Flv_inv_pre(GEN w, ulong p, ulong pi)
2327 { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2328
2329 /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2330 INLINE void
Flv_inv_indir(GEN w,GEN v,ulong p)2331 Flv_inv_indir(GEN w, GEN v, ulong p)
2332 {
2333 pari_sp av = avma;
2334 long n = lg(w), i;
2335 ulong u;
2336 GEN c;
2337
2338 if (n == 1) return;
2339 c = cgetg(n, t_VECSMALL); c[1] = w[1];
2340 for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2341 i = n-1; u = Fl_inv(c[i], p);
2342 for ( ; i > 1; --i)
2343 {
2344 ulong t = Fl_mul(u, c[i-1], p);
2345 u = Fl_mul(u, w[i], p); v[i] = t;
2346 }
2347 v[1] = u; set_avma(av);
2348 }
2349 static void
Flv_inv_i(GEN v,GEN w,ulong p)2350 Flv_inv_i(GEN v, GEN w, ulong p)
2351 {
2352 if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2353 else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2354 }
2355 void
Flv_inv_inplace(GEN v,ulong p)2356 Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2357 GEN
Flv_inv(GEN w,ulong p)2358 Flv_inv(GEN w, ulong p)
2359 { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2360
2361 GEN
Flx_div_by_X_x(GEN a,ulong x,ulong p,ulong * rem)2362 Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2363 {
2364 long l = lg(a), i;
2365 GEN a0, z0, z;
2366 if (l <= 3)
2367 {
2368 if (rem) *rem = l == 2? 0: a[2];
2369 return zero_Flx(a[1]);
2370 }
2371 z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2372 a0 = a + l-1;
2373 z0 = z + l-2; *z0 = *a0--;
2374 if (SMALL_ULONG(p))
2375 {
2376 for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2377 {
2378 ulong t = (*a0-- + x * *z0--) % p;
2379 *z0 = (long)t;
2380 }
2381 if (rem) *rem = (*a0 + x * *z0) % p;
2382 }
2383 else
2384 {
2385 for (i=l-3; i>1; i--)
2386 {
2387 ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2388 *z0 = (long)t;
2389 }
2390 if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2391 }
2392 return z;
2393 }
2394
2395 /* xa, ya = t_VECSMALL */
2396 static GEN
Flv_producttree(GEN xa,GEN s,ulong p,long vs)2397 Flv_producttree(GEN xa, GEN s, ulong p, long vs)
2398 {
2399 long n = lg(xa)-1;
2400 long m = n==1 ? 1: expu(n-1)+1;
2401 long i, j, k, ls = lg(s);
2402 GEN T = cgetg(m+1, t_VEC);
2403 GEN t = cgetg(ls, t_VEC);
2404 for (j=1, k=1; j<ls; k+=s[j++])
2405 gel(t, j) = s[j] == 1 ?
2406 mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2407 mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2408 Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2409 gel(T,1) = t;
2410 for (i=2; i<=m; i++)
2411 {
2412 GEN u = gel(T, i-1);
2413 long n = lg(u)-1;
2414 GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2415 for (j=1, k=1; k<n; j++, k+=2)
2416 gel(t, j) = Flx_mul(gel(u, k), gel(u, k+1), p);
2417 gel(T, i) = t;
2418 }
2419 return T;
2420 }
2421
2422 static GEN
Flx_Flv_multieval_tree(GEN P,GEN xa,GEN T,ulong p)2423 Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p)
2424 {
2425 long i,j,k;
2426 long m = lg(T)-1;
2427 GEN R = cgetg(lg(xa), t_VECSMALL);
2428 GEN Tp = cgetg(m+1, t_VEC), t;
2429 gel(Tp, m) = mkvec(P);
2430 for (i=m-1; i>=1; i--)
2431 {
2432 GEN u = gel(T, i), v = gel(Tp, i+1);
2433 long n = lg(u)-1;
2434 t = cgetg(n+1, t_VEC);
2435 for (j=1, k=1; k<n; j++, k+=2)
2436 {
2437 gel(t, k) = Flx_rem(gel(v, j), gel(u, k), p);
2438 gel(t, k+1) = Flx_rem(gel(v, j), gel(u, k+1), p);
2439 }
2440 gel(Tp, i) = t;
2441 }
2442 {
2443 GEN u = gel(T, i+1), v = gel(Tp, i+1);
2444 long n = lg(u)-1;
2445 for (j=1, k=1; j<=n; j++)
2446 {
2447 long c, d = degpol(gel(u,j));
2448 for (c=1; c<=d; c++, k++) R[k] = Flx_eval(gel(v, j), xa[k], p);
2449 }
2450 return gc_const((pari_sp)R, R);
2451 }
2452 }
2453
2454 static GEN
FlvV_polint_tree(GEN T,GEN R,GEN s,GEN xa,GEN ya,ulong p,long vs)2455 FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, long vs)
2456 {
2457 pari_sp av = avma;
2458 long m = lg(T)-1;
2459 long i, j, k, ls = lg(s);
2460 GEN Tp = cgetg(m+1, t_VEC);
2461 GEN t = cgetg(ls, t_VEC);
2462 for (j=1, k=1; j<ls; k+=s[j++])
2463 if (s[j]==2)
2464 {
2465 ulong a = Fl_mul(ya[k], R[k], p);
2466 ulong b = Fl_mul(ya[k+1], R[k+1], p);
2467 gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2468 Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2469 gel(t, j) = Flx_renormalize(gel(t, j), 4);
2470 }
2471 else
2472 gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2473 gel(Tp, 1) = t;
2474 for (i=2; i<=m; i++)
2475 {
2476 GEN u = gel(T, i-1);
2477 GEN t = cgetg(lg(gel(T,i)), t_VEC);
2478 GEN v = gel(Tp, i-1);
2479 long n = lg(v)-1;
2480 for (j=1, k=1; k<n; j++, k+=2)
2481 gel(t, j) = Flx_add(Flx_mul(gel(u, k), gel(v, k+1), p),
2482 Flx_mul(gel(u, k+1), gel(v, k), p), p);
2483 gel(Tp, i) = t;
2484 }
2485 return gerepileuptoleaf(av, gmael(Tp,m,1));
2486 }
2487
2488 GEN
Flx_Flv_multieval(GEN P,GEN xa,ulong p)2489 Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2490 {
2491 pari_sp av = avma;
2492 GEN s = producttree_scheme(lg(xa)-1);
2493 GEN T = Flv_producttree(xa, s, p, P[1]);
2494 return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p));
2495 }
2496
2497 static GEN
FlxV_Flv_multieval_tree(GEN x,GEN xa,GEN T,ulong p)2498 FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p)
2499 { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p)) }
2500
2501 GEN
FlxV_Flv_multieval(GEN P,GEN xa,ulong p)2502 FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2503 {
2504 pari_sp av = avma;
2505 GEN s = producttree_scheme(lg(xa)-1);
2506 GEN T = Flv_producttree(xa, s, p, P[1]);
2507 return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p));
2508 }
2509
2510 GEN
Flv_polint(GEN xa,GEN ya,ulong p,long vs)2511 Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2512 {
2513 pari_sp av = avma;
2514 GEN s = producttree_scheme(lg(xa)-1);
2515 GEN T = Flv_producttree(xa, s, p, vs);
2516 long m = lg(T)-1;
2517 GEN P = Flx_deriv(gmael(T, m, 1), p);
2518 GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2519 return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, vs));
2520 }
2521
2522 GEN
Flv_Flm_polint(GEN xa,GEN ya,ulong p,long vs)2523 Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2524 {
2525 pari_sp av = avma;
2526 GEN s = producttree_scheme(lg(xa)-1);
2527 GEN T = Flv_producttree(xa, s, p, vs);
2528 long i, m = lg(T)-1, l = lg(ya)-1;
2529 GEN P = Flx_deriv(gmael(T, m, 1), p);
2530 GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2531 GEN M = cgetg(l+1, t_VEC);
2532 for (i=1; i<=l; i++)
2533 gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, vs);
2534 return gerepileupto(av, M);
2535 }
2536
2537 GEN
Flv_invVandermonde(GEN L,ulong den,ulong p)2538 Flv_invVandermonde(GEN L, ulong den, ulong p)
2539 {
2540 pari_sp av = avma;
2541 long i, n = lg(L);
2542 GEN M, R;
2543 GEN s = producttree_scheme(n-1);
2544 GEN tree = Flv_producttree(L, s, p, 0);
2545 long m = lg(tree)-1;
2546 GEN T = gmael(tree, m, 1);
2547 R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p), p);
2548 if (den!=1) R = Flv_Fl_mul(R, den, p);
2549 M = cgetg(n, t_MAT);
2550 for (i = 1; i < n; i++)
2551 {
2552 GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2553 gel(M,i) = Flx_to_Flv(P, n-1);
2554 }
2555 return gerepilecopy(av, M);
2556 }
2557
2558 /***********************************************************************/
2559 /** **/
2560 /** Flxq **/
2561 /** **/
2562 /***********************************************************************/
2563 /* Flxq objects are defined as follows:
2564 They are Flx modulo another Flx called q.
2565 */
2566
2567 /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2568 GEN
Flxq_mul(GEN x,GEN y,GEN T,ulong p)2569 Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2570 {
2571 return Flx_rem(Flx_mul(x,y,p),T,p);
2572 }
2573
2574 /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2575 GEN
Flxq_sqr(GEN x,GEN T,ulong p)2576 Flxq_sqr(GEN x,GEN T,ulong p)
2577 {
2578 return Flx_rem(Flx_sqr(x,p),T,p);
2579 }
2580
2581 static GEN
_Flxq_red(void * E,GEN x)2582 _Flxq_red(void *E, GEN x)
2583 { struct _Flxq *s = (struct _Flxq *)E;
2584 return Flx_rem(x, s->T, s->p); }
2585 #if 0
2586 static GEN
2587 _Flx_sub(void *E, GEN x, GEN y)
2588 { struct _Flxq *s = (struct _Flxq *)E;
2589 return Flx_sub(x,y,s->p); }
2590 #endif
2591 static GEN
_Flxq_sqr(void * data,GEN x)2592 _Flxq_sqr(void *data, GEN x)
2593 {
2594 struct _Flxq *D = (struct _Flxq*)data;
2595 return Flxq_sqr(x, D->T, D->p);
2596 }
2597 static GEN
_Flxq_mul(void * data,GEN x,GEN y)2598 _Flxq_mul(void *data, GEN x, GEN y)
2599 {
2600 struct _Flxq *D = (struct _Flxq*)data;
2601 return Flxq_mul(x,y, D->T, D->p);
2602 }
2603 static GEN
_Flxq_one(void * data)2604 _Flxq_one(void *data)
2605 {
2606 struct _Flxq *D = (struct _Flxq*)data;
2607 return pol1_Flx(get_Flx_var(D->T));
2608 }
2609 #if 0
2610 static GEN
2611 _Flxq_zero(void *data)
2612 {
2613 struct _Flxq *D = (struct _Flxq*)data;
2614 return pol0_Flx(get_Flx_var(D->T));
2615 }
2616 static GEN
2617 _Flxq_cmul(void *data, GEN P, long a, GEN x)
2618 {
2619 struct _Flxq *D = (struct _Flxq*)data;
2620 return Flx_Fl_mul(x, P[a+2], D->p);
2621 }
2622 #endif
2623
2624 /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2625 GEN
Flxq_powu(GEN x,ulong n,GEN T,ulong p)2626 Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2627 {
2628 pari_sp av = avma;
2629 struct _Flxq D;
2630 GEN y;
2631 switch(n)
2632 {
2633 case 0: return pol1_Flx(get_Flx_var(T));
2634 case 1: return Flx_copy(x);
2635 case 2: return Flxq_sqr(x, T, p);
2636 }
2637 D.T = Flx_get_red(T, p); D.p = p;
2638 y = gen_powu_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2639 return gerepileuptoleaf(av, y);
2640 }
2641
2642 /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2643 GEN
Flxq_pow(GEN x,GEN n,GEN T,ulong p)2644 Flxq_pow(GEN x, GEN n, GEN T, ulong p)
2645 {
2646 pari_sp av = avma;
2647 struct _Flxq D;
2648 GEN y;
2649 long s = signe(n);
2650 if (!s) return pol1_Flx(get_Flx_var(T));
2651 if (s < 0)
2652 x = Flxq_inv(x,T,p);
2653 if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2654 D.T = Flx_get_red(T, p); D.p = p;
2655 y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2656 return gerepileuptoleaf(av, y);
2657 }
2658
2659 GEN
Flxq_pow_init(GEN x,GEN n,long k,GEN T,ulong p)2660 Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
2661 {
2662 struct _Flxq D;
2663 D.T = Flx_get_red(T, p); D.p = p;
2664 return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2665 }
2666
2667 GEN
Flxq_pow_table(GEN R,GEN n,GEN T,ulong p)2668 Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
2669 {
2670 struct _Flxq D;
2671 D.T = Flx_get_red(T, p); D.p = p;
2672 return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
2673 }
2674
2675 /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
2676 * not stack clean.
2677 */
2678 GEN
Flxq_invsafe(GEN x,GEN T,ulong p)2679 Flxq_invsafe(GEN x, GEN T, ulong p)
2680 {
2681 GEN V, z = Flx_extgcd(get_Flx_mod(T), x, p, NULL, &V);
2682 ulong iz;
2683 if (degpol(z)) return NULL;
2684 iz = Fl_inv (uel(z,2), p);
2685 return Flx_Fl_mul(V, iz, p);
2686 }
2687
2688 GEN
Flxq_inv(GEN x,GEN T,ulong p)2689 Flxq_inv(GEN x,GEN T,ulong p)
2690 {
2691 pari_sp av=avma;
2692 GEN U = Flxq_invsafe(x, T, p);
2693 if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
2694 return gerepileuptoleaf(av, U);
2695 }
2696
2697 GEN
Flxq_div(GEN x,GEN y,GEN T,ulong p)2698 Flxq_div(GEN x,GEN y,GEN T,ulong p)
2699 {
2700 pari_sp av = avma;
2701 return gerepileuptoleaf(av, Flxq_mul(x,Flxq_inv(y,T,p),T,p));
2702 }
2703
2704 GEN
Flxq_powers(GEN x,long l,GEN T,ulong p)2705 Flxq_powers(GEN x, long l, GEN T, ulong p)
2706 {
2707 struct _Flxq D;
2708 int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
2709 D.T = Flx_get_red(T, p); D.p = p;
2710 return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
2711 }
2712
2713 GEN
Flxq_matrix_pow(GEN y,long n,long m,GEN P,ulong l)2714 Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
2715 {
2716 return FlxV_to_Flm(Flxq_powers(y,m-1,P,l),n);
2717 }
2718
2719 GEN
Flx_Frobenius(GEN T,ulong p)2720 Flx_Frobenius(GEN T, ulong p)
2721 {
2722 return Flxq_powu(polx_Flx(get_Flx_var(T)), p, T, p);
2723 }
2724
2725 GEN
Flx_matFrobenius(GEN T,ulong p)2726 Flx_matFrobenius(GEN T, ulong p)
2727 {
2728 long n = get_Flx_degree(T);
2729 return Flxq_matrix_pow(Flx_Frobenius(T, p), n, n, T, p);
2730 }
2731
2732 static GEN
Flx_blocks_Flm(GEN P,long n,long m)2733 Flx_blocks_Flm(GEN P, long n, long m)
2734 {
2735 GEN z = cgetg(m+1,t_MAT);
2736 long i,j, k=2, l = lg(P);
2737 for(i=1; i<=m; i++)
2738 {
2739 GEN zi = cgetg(n+1,t_VECSMALL);
2740 gel(z,i) = zi;
2741 for(j=1; j<=n; j++)
2742 uel(zi, j) = k==l ? 0 : uel(P,k++);
2743 }
2744 return z;
2745 }
2746
2747 GEN
Flx_blocks(GEN P,long n,long m)2748 Flx_blocks(GEN P, long n, long m)
2749 {
2750 GEN z = cgetg(m+1,t_VEC);
2751 long i,j, k=2, l = lg(P);
2752 for(i=1; i<=m; i++)
2753 {
2754 GEN zi = cgetg(n+2,t_VECSMALL);
2755 zi[1] = P[1];
2756 gel(z,i) = zi;
2757 for(j=2; j<n+2; j++)
2758 uel(zi, j) = k==l ? 0 : uel(P,k++);
2759 zi = Flx_renormalize(zi, n+2);
2760 }
2761 return z;
2762 }
2763
2764 static GEN
FlxV_to_Flm_lg(GEN x,long m,long n)2765 FlxV_to_Flm_lg(GEN x, long m, long n)
2766 {
2767 long i;
2768 GEN y = cgetg(n+1, t_MAT);
2769 for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
2770 return y;
2771 }
2772
2773 GEN
Flx_FlxqV_eval(GEN Q,GEN x,GEN T,ulong p)2774 Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
2775 {
2776 pari_sp btop, av = avma;
2777 long sv = get_Flx_var(T), m = get_Flx_degree(T);
2778 long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
2779 GEN A, B, C, S, g;
2780 if (lQ == 0) return pol0_Flx(sv);
2781 if (lQ <= l)
2782 {
2783 n = l;
2784 d = 1;
2785 }
2786 else
2787 {
2788 n = l-1;
2789 d = (lQ+n-1)/n;
2790 }
2791 A = FlxV_to_Flm_lg(x, m, n);
2792 B = Flx_blocks_Flm(Q, n, d);
2793 C = gerepileupto(av, Flm_mul(A, B, p));
2794 g = gel(x, l);
2795 T = Flx_get_red(T, p);
2796 btop = avma;
2797 S = Flv_to_Flx(gel(C, d), sv);
2798 for (i = d-1; i>0; i--)
2799 {
2800 S = Flx_add(Flxq_mul(S, g, T, p), Flv_to_Flx(gel(C,i), sv), p);
2801 if (gc_needed(btop,1))
2802 S = gerepileuptoleaf(btop, S);
2803 }
2804 return gerepileuptoleaf(av, S);
2805 }
2806
2807 GEN
Flx_Flxq_eval(GEN Q,GEN x,GEN T,ulong p)2808 Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
2809 {
2810 pari_sp av = avma;
2811 GEN z, V;
2812 long d = degpol(Q), rtd;
2813 if (d < 0) return pol0_Flx(get_Flx_var(T));
2814 rtd = (long) sqrt((double)d);
2815 T = Flx_get_red(T, p);
2816 V = Flxq_powers(x, rtd, T, p);
2817 z = Flx_FlxqV_eval(Q, V, T, p);
2818 return gerepileupto(av, z);
2819 }
2820
2821 GEN
FlxC_FlxqV_eval(GEN x,GEN v,GEN T,ulong p)2822 FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
2823 { pari_APPLY_type(t_COL, Flx_FlxqV_eval(gel(x,i), v, T, p))
2824 }
2825
2826 GEN
FlxC_Flxq_eval(GEN x,GEN F,GEN T,ulong p)2827 FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
2828 {
2829 long d = brent_kung_optpow(degpol(T)-1,lg(x)-1,1);
2830 GEN Fp = Flxq_powers(F, d, T, p);
2831 return FlxC_FlxqV_eval(x, Fp, T, p);
2832 }
2833
2834 #if 0
2835 static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
2836 _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
2837 #endif
2838
2839 static GEN
Flxq_autpow_sqr(void * E,GEN x)2840 Flxq_autpow_sqr(void *E, GEN x)
2841 {
2842 struct _Flxq *D = (struct _Flxq*)E;
2843 return Flx_Flxq_eval(x, x, D->T, D->p);
2844 }
2845 static GEN
Flxq_autpow_msqr(void * E,GEN x)2846 Flxq_autpow_msqr(void *E, GEN x)
2847 {
2848 struct _Flxq *D = (struct _Flxq*)E;
2849 return Flx_FlxqV_eval(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p);
2850 }
2851
2852 GEN
Flxq_autpow(GEN x,ulong n,GEN T,ulong p)2853 Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
2854 {
2855 pari_sp av = avma;
2856 struct _Flxq D;
2857 long d;
2858 if (n==0) return Flx_rem(polx_Flx(x[1]), T, p);
2859 if (n==1) return Flx_rem(x, T, p);
2860 D.T = Flx_get_red(T, p); D.p = p;
2861 d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
2862 D.aut = Flxq_powers(x, d, T, p);
2863 x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
2864 return gerepilecopy(av, x);
2865 }
2866
2867 GEN
Flxq_autpowers(GEN x,ulong l,GEN T,ulong p)2868 Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
2869 {
2870 long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
2871 ulong i;
2872 pari_sp av = avma;
2873 GEN xp, V = cgetg(l+2,t_VEC);
2874 gel(V,1) = polx_Flx(vT); if (l==0) return V;
2875 gel(V,2) = gcopy(x); if (l==1) return V;
2876 T = Flx_get_red(T, p);
2877 d = brent_kung_optpow(dT-1, l-1, 1);
2878 xp = Flxq_powers(x, d, T, p);
2879 for(i = 3; i < l+2; i++)
2880 gel(V,i) = Flx_FlxqV_eval(gel(V,i-1), xp, T, p);
2881 return gerepilecopy(av, V);
2882 }
2883
2884 static GEN
Flxq_autsum_mul(void * E,GEN x,GEN y)2885 Flxq_autsum_mul(void *E, GEN x, GEN y)
2886 {
2887 struct _Flxq *D = (struct _Flxq*)E;
2888 GEN T = D->T;
2889 ulong p = D->p;
2890 GEN phi1 = gel(x,1), a1 = gel(x,2);
2891 GEN phi2 = gel(y,1), a2 = gel(y,2);
2892 ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
2893 GEN V2 = Flxq_powers(phi2, d, T, p);
2894 GEN phi3 = Flx_FlxqV_eval(phi1, V2, T, p);
2895 GEN aphi = Flx_FlxqV_eval(a1, V2, T, p);
2896 GEN a3 = Flxq_mul(aphi, a2, T, p);
2897 return mkvec2(phi3, a3);
2898 }
2899 static GEN
Flxq_autsum_sqr(void * E,GEN x)2900 Flxq_autsum_sqr(void *E, GEN x)
2901 { return Flxq_autsum_mul(E, x, x); }
2902
2903 GEN
Flxq_autsum(GEN x,ulong n,GEN T,ulong p)2904 Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
2905 {
2906 pari_sp av = avma;
2907 struct _Flxq D;
2908 D.T = Flx_get_red(T, p); D.p = p;
2909 x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
2910 return gerepilecopy(av, x);
2911 }
2912
2913 static GEN
Flxq_auttrace_mul(void * E,GEN x,GEN y)2914 Flxq_auttrace_mul(void *E, GEN x, GEN y)
2915 {
2916 struct _Flxq *D = (struct _Flxq*)E;
2917 GEN T = D->T;
2918 ulong p = D->p;
2919 GEN phi1 = gel(x,1), a1 = gel(x,2);
2920 GEN phi2 = gel(y,1), a2 = gel(y,2);
2921 ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
2922 GEN V1 = Flxq_powers(phi1, d, T, p);
2923 GEN phi3 = Flx_FlxqV_eval(phi2, V1, T, p);
2924 GEN aphi = Flx_FlxqV_eval(a2, V1, T, p);
2925 GEN a3 = Flx_add(a1, aphi, p);
2926 return mkvec2(phi3, a3);
2927 }
2928
2929 static GEN
Flxq_auttrace_sqr(void * E,GEN x)2930 Flxq_auttrace_sqr(void *E, GEN x)
2931 { return Flxq_auttrace_mul(E, x, x); }
2932
2933 GEN
Flxq_auttrace(GEN x,ulong n,GEN T,ulong p)2934 Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
2935 {
2936 pari_sp av = avma;
2937 struct _Flxq D;
2938 D.T = Flx_get_red(T, p); D.p = p;
2939 x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
2940 return gerepilecopy(av, x);
2941 }
2942
2943 static long
bounded_order(ulong p,GEN b,long k)2944 bounded_order(ulong p, GEN b, long k)
2945 {
2946 long i;
2947 GEN a=modii(utoi(p),b);
2948 for(i=1;i<k;i++)
2949 {
2950 if (equali1(a))
2951 return i;
2952 a = modii(muliu(a,p),b);
2953 }
2954 return 0;
2955 }
2956
2957 /*
2958 n = (p^d-a)\b
2959 b = bb*p^vb
2960 p^k = 1 [bb]
2961 d = m*k+r+vb
2962 u = (p^k-1)/bb;
2963 v = (p^(r+vb)-a)/b;
2964 w = (p^(m*k)-1)/(p^k-1)
2965 n = p^r*w*u+v
2966 w*u = p^vb*(p^(m*k)-1)/b
2967 n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b
2968 */
2969
2970 static GEN
Flxq_pow_Frobenius(GEN x,GEN n,GEN aut,GEN T,ulong p)2971 Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p)
2972 {
2973 pari_sp av=avma;
2974 long d = get_Flx_degree(T);
2975 GEN an = absi_shallow(n), z, q;
2976 if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow(x, n, T, p);
2977 q = powuu(p, d);
2978 if (dvdii(q, n))
2979 {
2980 long vn = logint(an,utoi(p));
2981 GEN autvn = vn==1 ? aut: Flxq_autpow(aut,vn,T,p);
2982 z = Flx_Flxq_eval(x,autvn,T,p);
2983 } else
2984 {
2985 GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
2986 GEN bb, u, v, autk;
2987 long vb = Z_lvalrem(b,p,&bb);
2988 long m, r, k = is_pm1(bb) ? 1 : bounded_order(p,bb,d);
2989 if (!k || d-vb<k) return Flxq_pow(x,n, T, p);
2990 m = (d-vb)/k; r = (d-vb)%k;
2991 u = diviiexact(subiu(powuu(p,k),1),bb);
2992 v = diviiexact(subii(powuu(p,r+vb),a),b);
2993 autk = k==1 ? aut: Flxq_autpow(aut,k,T,p);
2994 if (r)
2995 {
2996 GEN autr = r==1 ? aut: Flxq_autpow(aut,r,T,p);
2997 z = Flx_Flxq_eval(x,autr,T,p);
2998 } else z = x;
2999 if (m > 1) z = gel(Flxq_autsum(mkvec2(autk, z), m, T, p), 2);
3000 if (!is_pm1(u)) z = Flxq_pow(z, u, T, p);
3001 if (signe(v)) z = Flxq_mul(z, Flxq_pow(x, v, T, p), T, p);
3002 }
3003 return gerepileupto(av,signe(n)>0 ? z : Flxq_inv(z,T,p));
3004 }
3005
3006 static GEN
_Flxq_pow(void * data,GEN x,GEN n)3007 _Flxq_pow(void *data, GEN x, GEN n)
3008 {
3009 struct _Flxq *D = (struct _Flxq*)data;
3010 return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p);
3011 }
3012
3013 static GEN
_Flxq_rand(void * data)3014 _Flxq_rand(void *data)
3015 {
3016 pari_sp av=avma;
3017 struct _Flxq *D = (struct _Flxq*)data;
3018 GEN z;
3019 do
3020 {
3021 set_avma(av);
3022 z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3023 } while (lgpol(z)==0);
3024 return z;
3025 }
3026
3027 /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3028 static GEN
Fl_Flxq_log(ulong a,GEN g,GEN o,GEN T,ulong p)3029 Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3030 {
3031 pari_sp av = avma;
3032 GEN q,n_q,ord,ordp, op;
3033
3034 if (a == 1UL) return gen_0;
3035 /* p > 2 */
3036
3037 ordp = utoi(p - 1);
3038 ord = get_arith_Z(o);
3039 if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3040 if (a == p - 1) /* -1 */
3041 return gerepileuptoint(av, shifti(ord,-1));
3042 ordp = gcdii(ordp, ord);
3043 op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3044
3045 q = NULL;
3046 if (T)
3047 { /* we want < g > = Fp^* */
3048 if (!equalii(ord,ordp)) {
3049 q = diviiexact(ord,ordp);
3050 g = Flxq_pow(g,q,T,p);
3051 }
3052 }
3053 n_q = Fp_log(utoi(a), utoi(uel(g,2)), op, utoi(p));
3054 if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
3055 if (q) n_q = mulii(q, n_q);
3056 return gerepileuptoint(av, n_q);
3057 }
3058
3059 static GEN
Flxq_easylog(void * E,GEN a,GEN g,GEN ord)3060 Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3061 {
3062 struct _Flxq *f = (struct _Flxq *)E;
3063 GEN T = f->T;
3064 ulong p = f->p;
3065 long d = get_Flx_degree(T);
3066 if (Flx_equal1(a)) return gen_0;
3067 if (Flx_equal(a,g)) return gen_1;
3068 if (!degpol(a))
3069 return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3070 if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3071 return NULL;
3072 return Flxq_log_index(a, g, ord, T, p);
3073 }
3074
3075 static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3076
3077 const struct bb_group *
get_Flxq_star(void ** E,GEN T,ulong p)3078 get_Flxq_star(void **E, GEN T, ulong p)
3079 {
3080 struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3081 e->T = T; e->p = p; e->aut = Flx_Frobenius(T, p);
3082 *E = (void*)e; return &Flxq_star;
3083 }
3084
3085 GEN
Flxq_order(GEN a,GEN ord,GEN T,ulong p)3086 Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3087 {
3088 void *E;
3089 const struct bb_group *S = get_Flxq_star(&E,T,p);
3090 return gen_order(a,ord,E,S);
3091 }
3092
3093 GEN
Flxq_log(GEN a,GEN g,GEN ord,GEN T,ulong p)3094 Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3095 {
3096 void *E;
3097 pari_sp av = avma;
3098 const struct bb_group *S = get_Flxq_star(&E,T,p);
3099 GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3100 if (Flxq_log_use_index(gel(F,lg(F)-1), T, p))
3101 v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3102 return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
3103 }
3104
3105 GEN
Flxq_sqrtn(GEN a,GEN n,GEN T,ulong p,GEN * zeta)3106 Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3107 {
3108 if (!lgpol(a))
3109 {
3110 if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3111 if (zeta)
3112 *zeta=pol1_Flx(get_Flx_var(T));
3113 return pol0_Flx(get_Flx_var(T));
3114 }
3115 else
3116 {
3117 void *E;
3118 pari_sp av = avma;
3119 const struct bb_group *S = get_Flxq_star(&E,T,p);
3120 GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3121 GEN s = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
3122 if (s) gerepileall(av, zeta?2:1, &s, zeta);
3123 return s;
3124 }
3125 }
3126
3127 GEN
Flxq_sqrt(GEN a,GEN T,ulong p)3128 Flxq_sqrt(GEN a, GEN T, ulong p)
3129 {
3130 return Flxq_sqrtn(a, gen_2, T, p, NULL);
3131 }
3132
3133 /* assume T irreducible mod p */
3134 int
Flxq_issquare(GEN x,GEN T,ulong p)3135 Flxq_issquare(GEN x, GEN T, ulong p)
3136 {
3137 if (lgpol(x) == 0 || p == 2) return 1;
3138 return krouu(Flxq_norm(x,T,p), p) == 1;
3139 }
3140
3141 /* assume T irreducible mod p */
3142 int
Flxq_is2npower(GEN x,long n,GEN T,ulong p)3143 Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3144 {
3145 pari_sp av;
3146 GEN m;
3147 if (n==1) return Flxq_issquare(x, T, p);
3148 if (lgpol(x) == 0 || p == 2) return 1;
3149 av = avma;
3150 m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3151 return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3152 }
3153
3154 GEN
Flxq_lroot_fast(GEN a,GEN sqx,GEN T,long p)3155 Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
3156 {
3157 pari_sp av=avma;
3158 GEN A = Flx_splitting(a,p);
3159 return gerepileuptoleaf(av, FlxqV_dotproduct(A,sqx,T,p));
3160 }
3161
3162 GEN
Flxq_lroot(GEN a,GEN T,long p)3163 Flxq_lroot(GEN a, GEN T, long p)
3164 {
3165 pari_sp av=avma;
3166 long n = get_Flx_degree(T), d = degpol(a);
3167 GEN sqx, V;
3168 if (n==1) return leafcopy(a);
3169 if (n==2) return Flxq_powu(a, p, T, p);
3170 sqx = Flxq_autpow(Flx_Frobenius(T, p), n-1, T, p);
3171 if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
3172 if (d>=p)
3173 {
3174 V = Flxq_powers(sqx,p-1,T,p);
3175 return gerepileuptoleaf(av, Flxq_lroot_fast(a,V,T,p));
3176 } else
3177 return gerepileuptoleaf(av, Flx_Flxq_eval(a,sqx,T,p));
3178 }
3179
3180 ulong
Flxq_norm(GEN x,GEN TB,ulong p)3181 Flxq_norm(GEN x, GEN TB, ulong p)
3182 {
3183 GEN T = get_Flx_mod(TB);
3184 ulong y = Flx_resultant(T, x, p);
3185 ulong L = Flx_lead(T);
3186 if ( L==1 || lgpol(x)==0) return y;
3187 return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3188 }
3189
3190 ulong
Flxq_trace(GEN x,GEN TB,ulong p)3191 Flxq_trace(GEN x, GEN TB, ulong p)
3192 {
3193 pari_sp av = avma;
3194 ulong t;
3195 GEN T = get_Flx_mod(TB);
3196 long n = degpol(T)-1;
3197 GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3198 t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3199 return gc_ulong(av, t);
3200 }
3201
3202 /*x must be reduced*/
3203 GEN
Flxq_charpoly(GEN x,GEN TB,ulong p)3204 Flxq_charpoly(GEN x, GEN TB, ulong p)
3205 {
3206 pari_sp ltop=avma;
3207 GEN T = get_Flx_mod(TB);
3208 long vs = evalvarn(fetch_var());
3209 GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3210 GEN r = Flx_FlxY_resultant(T, xm1, p);
3211 r[1] = x[1];
3212 (void)delete_var(); return gerepileupto(ltop, r);
3213 }
3214
3215 /* Computing minimal polynomial : */
3216 /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3217 /* in Algebraic Extensions of Finite Fields' */
3218
3219 /* Let v a linear form, return the linear form z->v(tau*z)
3220 that is, v*(M_tau) */
3221
3222 static GEN
Flxq_transmul_init(GEN tau,GEN T,ulong p)3223 Flxq_transmul_init(GEN tau, GEN T, ulong p)
3224 {
3225 GEN bht;
3226 GEN h, Tp = get_Flx_red(T, &h);
3227 long n = degpol(Tp), vT = Tp[1];
3228 GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3229 GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3230 ft[1] = vT; bt[1] = vT;
3231 if (h)
3232 bht = Flxn_mul(bt, h, n-1, p);
3233 else
3234 {
3235 GEN bh = Flx_div(Flx_shift(tau, n-1), T, p);
3236 bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3237 bht[1] = vT;
3238 }
3239 return mkvec3(bt, bht, ft);
3240 }
3241
3242 static GEN
Flxq_transmul(GEN tau,GEN a,long n,ulong p)3243 Flxq_transmul(GEN tau, GEN a, long n, ulong p)
3244 {
3245 pari_sp ltop = avma;
3246 GEN t1, t2, t3, vec;
3247 GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3248 if (lgpol(a)==0) return pol0_Flx(a[1]);
3249 t2 = Flx_shift(Flx_mul(bt, a, p),1-n);
3250 if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
3251 t1 = Flx_shift(Flx_mul(ft, a, p),-n);
3252 t3 = Flxn_mul(t1, bht, n-1, p);
3253 vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3254 return gerepileuptoleaf(ltop, vec);
3255 }
3256
3257 GEN
Flxq_minpoly(GEN x,GEN T,ulong p)3258 Flxq_minpoly(GEN x, GEN T, ulong p)
3259 {
3260 pari_sp ltop = avma;
3261 long vT = get_Flx_var(T), n = get_Flx_degree(T);
3262 GEN v_x;
3263 GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3264 T = Flx_get_red(T, p);
3265 v_x = Flxq_powers(x, usqrt(2*n), T, p);
3266 while (lgpol(tau) != 0)
3267 {
3268 long i, j, m, k1;
3269 GEN M, v, tr;
3270 GEN g_prime, c;
3271 if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3272 v = random_Flx(n, vT, p);
3273 tr = Flxq_transmul_init(tau, T, p);
3274 v = Flxq_transmul(tr, v, n, p);
3275 m = 2*(n-degpol(g));
3276 k1 = usqrt(m);
3277 tr = Flxq_transmul_init(gel(v_x,k1+1), T, p);
3278 c = cgetg(m+2,t_VECSMALL);
3279 c[1] = vT;
3280 for (i=0; i<m; i+=k1)
3281 {
3282 long mj = minss(m-i, k1);
3283 for (j=0; j<mj; j++)
3284 uel(c,m+1-(i+j)) = Flx_dotproduct(v, gel(v_x,j+1), p);
3285 v = Flxq_transmul(tr, v, n, p);
3286 }
3287 c = Flx_renormalize(c, m+2);
3288 /* now c contains <v,x^i> , i = 0..m-1 */
3289 M = Flx_halfgcd(monomial_Flx(1, m, vT), c, p);
3290 g_prime = gmael(M, 2, 2);
3291 if (degpol(g_prime) < 1) continue;
3292 g = Flx_mul(g, g_prime, p);
3293 tau = Flxq_mul(tau, Flx_FlxqV_eval(g_prime, v_x, T, p), T, p);
3294 }
3295 g = Flx_normalize(g,p);
3296 return gerepileuptoleaf(ltop,g);
3297 }
3298
3299 GEN
Flxq_conjvec(GEN x,GEN T,ulong p)3300 Flxq_conjvec(GEN x, GEN T, ulong p)
3301 {
3302 long i, l = 1+get_Flx_degree(T);
3303 GEN z = cgetg(l,t_COL);
3304 T = Flx_get_red(T,p);
3305 gel(z,1) = Flx_copy(x);
3306 for (i=2; i<l; i++) gel(z,i) = Flxq_powu(gel(z,i-1), p, T, p);
3307 return z;
3308 }
3309
3310 GEN
gener_Flxq(GEN T,ulong p,GEN * po)3311 gener_Flxq(GEN T, ulong p, GEN *po)
3312 {
3313 long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
3314 ulong p_1;
3315 GEN g, L, L2, o, q, F;
3316 pari_sp av0, av;
3317
3318 if (f == 1) {
3319 GEN fa;
3320 o = utoipos(p-1);
3321 fa = Z_factor(o);
3322 L = gel(fa,1);
3323 L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
3324 g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
3325 if (po) *po = mkvec2(o, fa);
3326 return g;
3327 }
3328
3329 av0 = avma; p_1 = p - 1;
3330 q = diviuexact(subiu(powuu(p,f), 1), p_1);
3331
3332 L = cgetg(1, t_VECSMALL);
3333 if (p > 3)
3334 {
3335 ulong t = p_1 >> vals(p_1);
3336 GEN P = gel(factoru(t), 1);
3337 L = cgetg_copy(P, &i);
3338 while (--i) L[i] = p_1 / P[i];
3339 }
3340 o = factor_pn_1(utoipos(p),f);
3341 L2 = leafcopy( gel(o, 1) );
3342 for (i = j = 1; i < lg(L2); i++)
3343 {
3344 if (umodui(p_1, gel(L2,i)) == 0) continue;
3345 gel(L2,j++) = diviiexact(q, gel(L2,i));
3346 }
3347 setlg(L2, j);
3348 F = Flx_Frobenius(T, p);
3349 for (av = avma;; set_avma(av))
3350 {
3351 GEN tt;
3352 g = random_Flx(f, vT, p);
3353 if (degpol(g) < 1) continue;
3354 if (p == 2) tt = g;
3355 else
3356 {
3357 ulong t = Flxq_norm(g, T, p);
3358 if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
3359 tt = Flxq_powu(g, p_1>>1, T, p);
3360 }
3361 for (i = 1; i < j; i++)
3362 {
3363 GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p);
3364 if (!degpol(a) && uel(a,2) == p_1) break;
3365 }
3366 if (i == j) break;
3367 }
3368 if (!po)
3369 {
3370 set_avma((pari_sp)g);
3371 g = gerepileuptoleaf(av0, g);
3372 }
3373 else {
3374 *po = mkvec2(subiu(powuu(p,f), 1), o);
3375 gerepileall(av0, 2, &g, po);
3376 }
3377 return g;
3378 }
3379
3380 static GEN
_Flxq_neg(void * E,GEN x)3381 _Flxq_neg(void *E, GEN x)
3382 { struct _Flxq *s = (struct _Flxq *)E;
3383 return Flx_neg(x,s->p); }
3384
3385 static GEN
_Flxq_rmul(void * E,GEN x,GEN y)3386 _Flxq_rmul(void *E, GEN x, GEN y)
3387 { struct _Flxq *s = (struct _Flxq *)E;
3388 return Flx_mul(x,y,s->p); }
3389
3390 static GEN
_Flxq_inv(void * E,GEN x)3391 _Flxq_inv(void *E, GEN x)
3392 { struct _Flxq *s = (struct _Flxq *)E;
3393 return Flxq_inv(x,s->T,s->p); }
3394
3395 static int
_Flxq_equal0(GEN x)3396 _Flxq_equal0(GEN x) { return lgpol(x)==0; }
3397
3398 static GEN
_Flxq_s(void * E,long x)3399 _Flxq_s(void *E, long x)
3400 { struct _Flxq *s = (struct _Flxq *)E;
3401 ulong u = x<0 ? s->p+x: (ulong)x;
3402 return Fl_to_Flx(u, get_Flx_var(s->T));
3403 }
3404
3405 static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
3406 _Flxq_inv,_Flxq_equal0,_Flxq_s};
3407
get_Flxq_field(void ** E,GEN T,ulong p)3408 const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
3409 {
3410 GEN z = new_chunk(sizeof(struct _Flxq));
3411 struct _Flxq *e = (struct _Flxq *) z;
3412 e->T = Flx_get_red(T, p); e->p = p; *E = (void*)e;
3413 return &Flxq_field;
3414 }
3415
3416 /***********************************************************************/
3417 /** **/
3418 /** Flxn **/
3419 /** **/
3420 /***********************************************************************/
3421
3422 GEN
Flx_invLaplace(GEN x,ulong p)3423 Flx_invLaplace(GEN x, ulong p)
3424 {
3425 long i, d = degpol(x);
3426 ulong t;
3427 GEN y;
3428 if (d <= 1) return Flx_copy(x);
3429 t = Fl_inv(factorial_Fl(d, p), p);
3430 y = cgetg(d+3, t_VECSMALL);
3431 y[1] = x[1];
3432 for (i=d; i>=2; i--)
3433 {
3434 uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3435 t = Fl_mul(t, i, p);
3436 }
3437 uel(y,3) = uel(x,3);
3438 uel(y,2) = uel(x,2);
3439 return y;
3440 }
3441
3442 GEN
Flx_Laplace(GEN x,ulong p)3443 Flx_Laplace(GEN x, ulong p)
3444 {
3445 long i, d = degpol(x);
3446 ulong t = 1;
3447 GEN y;
3448 if (d <= 1) return Flx_copy(x);
3449 y = cgetg(d+3, t_VECSMALL);
3450 y[1] = x[1];
3451 uel(y,2) = uel(x,2);
3452 uel(y,3) = uel(x,3);
3453 for (i=2; i<=d; i++)
3454 {
3455 t = Fl_mul(t, i%p, p);
3456 uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3457 }
3458 return y;
3459 }
3460
3461 GEN
Flxn_red(GEN a,long n)3462 Flxn_red(GEN a, long n)
3463 {
3464 long i, L, l = lg(a);
3465 GEN b;
3466 if (l == 2 || !n) return zero_Flx(a[1]);
3467 L = n+2; if (L > l) L = l;
3468 b = cgetg(L, t_VECSMALL); b[1] = a[1];
3469 for (i=2; i<L; i++) b[i] = a[i];
3470 return Flx_renormalize(b,L);
3471 }
3472
3473 GEN
Flxn_mul(GEN a,GEN b,long n,ulong p)3474 Flxn_mul(GEN a, GEN b, long n, ulong p)
3475 { return Flxn_red(Flx_mul(a, b, p), n); }
3476
3477 GEN
Flxn_sqr(GEN a,long n,ulong p)3478 Flxn_sqr(GEN a, long n, ulong p)
3479 { return Flxn_red(Flx_sqr(a, p), n); }
3480
3481 /* (f*g) \/ x^n */
3482 static GEN
Flx_mulhigh_i(GEN f,GEN g,long n,ulong p)3483 Flx_mulhigh_i(GEN f, GEN g, long n, ulong p)
3484 {
3485 return Flx_shift(Flx_mul(f,g, p),-n);
3486 }
3487
3488 static GEN
Flxn_mulhigh(GEN f,GEN g,long n2,long n,ulong p)3489 Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p)
3490 {
3491 GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
3492 return Flx_add(Flx_mulhigh_i(fl, g, n2, p), Flxn_mul(fh, g, n - n2, p), p);
3493 }
3494
3495 GEN
Flxn_inv(GEN f,long e,ulong p)3496 Flxn_inv(GEN f, long e, ulong p)
3497 {
3498 pari_sp av = avma, av2;
3499 ulong mask;
3500 GEN W;
3501 long n=1;
3502 if (lg(f)==2) pari_err_INV("Flxn_inv",f);
3503 W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
3504 mask = quadratic_prec_mask(e);
3505 av2 = avma;
3506 for (;mask>1;)
3507 {
3508 GEN u, fr;
3509 long n2 = n;
3510 n<<=1; if (mask & 1) n--;
3511 mask >>= 1;
3512 fr = Flxn_red(f, n);
3513 u = Flxn_mul(W, Flxn_mulhigh(fr, W, n2, n, p), n-n2, p);
3514 W = Flx_sub(W, Flx_shift(u, n2), p);
3515 if (gc_needed(av2,2))
3516 {
3517 if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_inv, e = %ld", n);
3518 W = gerepileupto(av2, W);
3519 }
3520 }
3521 return gerepileupto(av, W);
3522 }
3523
3524 GEN
Flxn_expint(GEN h,long e,ulong p)3525 Flxn_expint(GEN h, long e, ulong p)
3526 {
3527 pari_sp av = avma, av2;
3528 long v = h[1], n=1;
3529 GEN f = pol1_Flx(v), g = pol1_Flx(v);
3530 ulong mask = quadratic_prec_mask(e);
3531 av2 = avma;
3532 for (;mask>1;)
3533 {
3534 GEN u, w;
3535 long n2 = n;
3536 n<<=1; if (mask & 1) n--;
3537 mask >>= 1;
3538 u = Flxn_mul(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p), n-n2, p);
3539 u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
3540 w = Flxn_mul(f, Flx_integXn(u, n2-1, p), n-n2, p);
3541 f = Flx_add(f, Flx_shift(w, n2), p);
3542 if (mask<=1) break;
3543 u = Flxn_mul(g, Flxn_mulhigh(f, g, n2, n, p), n-n2, p);
3544 g = Flx_sub(g, Flx_shift(u, n2), p);
3545 if (gc_needed(av2,2))
3546 {
3547 if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
3548 gerepileall(av2, 2, &f, &g);
3549 }
3550 }
3551 return gerepileupto(av, f);
3552 }
3553
3554 GEN
Flxn_exp(GEN h,long e,ulong p)3555 Flxn_exp(GEN h, long e, ulong p)
3556 {
3557 if (degpol(h)<1 || uel(h,2)!=0)
3558 pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
3559 return Flxn_expint(Flx_deriv(h, p), e, p);
3560 }
3561
3562 INLINE GEN
Flxn_recip(GEN x,long n)3563 Flxn_recip(GEN x, long n)
3564 {
3565 GEN z=Flx_recipspec(x+2,lgpol(x),n);
3566 z[1]=x[1];
3567 return z;
3568 }
3569
3570 GEN
Flx_Newton(GEN P,long n,ulong p)3571 Flx_Newton(GEN P, long n, ulong p)
3572 {
3573 pari_sp av = avma;
3574 long d = degpol(P);
3575 GEN dP = Flxn_recip(Flx_deriv(P, p), d);
3576 GEN Q = Flxn_mul(Flxn_inv(Flxn_recip(P, d+1), n, p), dP, n, p);
3577 return gerepileuptoleaf(av, Q);
3578 }
3579
3580 GEN
Flx_fromNewton(GEN P,ulong p)3581 Flx_fromNewton(GEN P, ulong p)
3582 {
3583 pari_sp av = avma;
3584 ulong n = Flx_constant(P)+1;
3585 GEN z = Flx_neg(Flx_shift(P, -1), p);
3586 GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
3587 return gerepileuptoleaf(av, Q);
3588 }
3589
3590 static long
newtonlogint(ulong n,ulong pp)3591 newtonlogint(ulong n, ulong pp)
3592 {
3593 long s = 0;
3594 while (n > pp)
3595 {
3596 s += ulogint(n, pp);
3597 n = (n+1)>>1;
3598 }
3599 return s;
3600 }
3601
3602 static void
init_invlaplace(long d,ulong p,GEN * pt_P,GEN * pt_V)3603 init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
3604 {
3605 long i;
3606 ulong e;
3607 GEN P = cgetg(d+1, t_VECSMALL);
3608 GEN V = cgetg(d+1, t_VECSMALL);
3609 for (i=1, e=1; i<=d; i++, e++)
3610 {
3611 if (e==p)
3612 {
3613 e = 0;
3614 V[i] = u_lvalrem(i, p, &uel(P,i));
3615 } else
3616 {
3617 V[i] = 0; uel(P,i) = i;
3618 }
3619 }
3620 *pt_P = P; *pt_V = V;
3621 }
3622
3623 /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
3624 * val large enough to compensate for the power of p in the factorials */
3625
3626 static GEN
ZpX_invLaplace_init(long n,GEN q,ulong p,long v,long var)3627 ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long var)
3628 {
3629 pari_sp av = avma;
3630 long i, d = n-1, w;
3631 GEN y, W, E, t;
3632 init_invlaplace(d, p, &E, &W);
3633 t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
3634 w = zv_sum(W);
3635 if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
3636 y = cgetg(d+3,t_POL);
3637 y[1] = evalsigne(1) | evalvarn(var);
3638 for (i=d; i>=1; i--)
3639 {
3640 gel(y,i+2) = t;
3641 t = Fp_mulu(t, uel(E,i), q);
3642 if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
3643 }
3644 gel(y,2) = t;
3645 return gerepilecopy(av, ZX_renormalize(y, d+3));
3646 }
3647
3648 static GEN
Flx_composedsum(GEN P,GEN Q,ulong p)3649 Flx_composedsum(GEN P, GEN Q, ulong p)
3650 {
3651 long n = 1 + degpol(P)*degpol(Q);
3652 ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
3653 Fl_powu(Flx_lead(Q), degpol(P), p), p);
3654 GEN R;
3655 if (p >= (ulong)n)
3656 {
3657 GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
3658 GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
3659 GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
3660 R = Flx_fromNewton(L, p);
3661 } else
3662 {
3663 long v = factorial_lval(n-1, p);
3664 long w = 1 + newtonlogint(n-1, p);
3665 GEN pv = powuu(p, v);
3666 GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
3667 GEN iL = ZpX_invLaplace_init(n, q, p, v, varn(P));
3668 GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
3669 GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
3670 GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
3671 GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
3672 R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
3673 }
3674 return Flx_Fl_mul(R, lead, p);
3675 }
3676
3677 GEN
Flx_direct_compositum(GEN a,GEN b,ulong p)3678 Flx_direct_compositum(GEN a, GEN b, ulong p)
3679 {
3680 return Flx_composedsum(a, b, p);
3681 }
3682
3683 static GEN
_Flx_direct_compositum(void * E,GEN a,GEN b)3684 _Flx_direct_compositum(void *E, GEN a, GEN b)
3685 { return Flx_direct_compositum(a, b, (ulong)E); }
3686
3687 GEN
FlxV_direct_compositum(GEN V,ulong p)3688 FlxV_direct_compositum(GEN V, ulong p)
3689 {
3690 return gen_product(V, (void *)p, &_Flx_direct_compositum);
3691 }
3692
3693 /* (x+1)^n mod p; assume 2 <= n < 2p prime */
3694 static GEN
Fl_Xp1_powu(ulong n,ulong p,long v)3695 Fl_Xp1_powu(ulong n, ulong p, long v)
3696 {
3697 ulong k, d = (n + 1) >> 1;
3698 GEN C, V = identity_zv(d);
3699
3700 Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
3701 C = cgetg(n+3, t_VECSMALL);
3702 C[1] = v;
3703 uel(C,2) = 1UL;
3704 uel(C,3) = n%p;
3705 uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
3706 /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
3707 if (SMALL_ULONG(p))
3708 for (k = 3; k <= d; k++)
3709 uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
3710 else
3711 {
3712 ulong pi = get_Fl_red(p);
3713 for (k = 3; k <= d; k++)
3714 uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
3715 }
3716 for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
3717 return C; /* normalized */
3718 }
3719
3720 /* p arbitrary */
3721 GEN
Flx_translate1_basecase(GEN P,ulong p)3722 Flx_translate1_basecase(GEN P, ulong p)
3723 {
3724 GEN R = Flx_copy(P);
3725 long i, k, n = degpol(P);
3726 for (i = 1; i <= n; i++)
3727 for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
3728 return R;
3729 }
3730
3731 static int
translate_basecase(long n,ulong p)3732 translate_basecase(long n, ulong p)
3733 {
3734 #ifdef LONG_IS_64BIT
3735 if (p <= 19) return n < 40;
3736 if (p < 1UL<<30) return n < 58;
3737 if (p < 1UL<<59) return n < 100;
3738 if (p < 1UL<<62) return n < 120;
3739 if (p < 1UL<<63) return n < 240;
3740 return n < 250;
3741 #else
3742 if (p <= 13) return n < 18;
3743 if (p <= 17) return n < 22;
3744 if (p <= 29) return n < 39;
3745 if (p <= 67) return n < 69;
3746 if (p < 1UL<< 15) return n < 80;
3747 if (p < 1UL<< 16) return n < 100;
3748 if (p < 1UL<< 28) return n < 300;
3749 return n < 650;
3750 #endif
3751 }
3752 /* assume p prime */
3753 GEN
Flx_translate1(GEN P,ulong p)3754 Flx_translate1(GEN P, ulong p)
3755 {
3756 long d, n = degpol(P);
3757 GEN R, Q, S;
3758 if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
3759 /* n > 0 */
3760 d = n >> 1;
3761 if ((ulong)n < p)
3762 {
3763 R = Flx_translate1(Flxn_red(P, d), p);
3764 Q = Flx_translate1(Flx_shift(P, -d), p);
3765 S = Fl_Xp1_powu(d, p, P[1]);
3766 return Flx_add(Flx_mul(Q, S, p), R, p);
3767 }
3768 else
3769 {
3770 ulong q;
3771 if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
3772 R = Flx_translate1(Flxn_red(P, q), p);
3773 Q = Flx_translate1(Flx_shift(P, -q), p);
3774 S = Flx_add(Flx_shift(Q, q), Q, p);
3775 return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
3776 }
3777 }
3778
3779 static GEN
zl_Xp1_powu(ulong n,ulong p,ulong q,long e,long vs)3780 zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
3781 {
3782 ulong k, d = n >> 1, c, v = 0;
3783 GEN C, V, W, U = upowers(p, e-1);
3784 init_invlaplace(d, p, &V, &W);
3785 Flv_inv_inplace(V, q);
3786 C = cgetg(n+3, t_VECSMALL);
3787 C[1] = vs;
3788 uel(C,2) = 1UL;
3789 uel(C,3) = n%q;
3790 v = u_lvalrem(n, p, &c);
3791 for (k = 2; k <= d; k++)
3792 {
3793 ulong w;
3794 v += u_lvalrem(n-k+1, p, &w) - W[k];
3795 c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
3796 uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
3797 }
3798 for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
3799 return C; /* normalized */
3800 }
3801
3802 GEN
zlx_translate1(GEN P,ulong p,long e)3803 zlx_translate1(GEN P, ulong p, long e)
3804 {
3805 ulong d, q = upowuu(p,e), n = degpol(P);
3806 GEN R, Q, S;
3807 if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
3808 /* n > 0 */
3809 d = n >> 1;
3810 R = zlx_translate1(Flxn_red(P, d), p, e);
3811 Q = zlx_translate1(Flx_shift(P, -d), p, e);
3812 S = zl_Xp1_powu(d, p, q, e, P[1]);
3813 return Flx_add(Flx_mul(Q, S, q), R, q);
3814 }
3815
3816 /***********************************************************************/
3817 /** **/
3818 /** Fl2 **/
3819 /** **/
3820 /***********************************************************************/
3821 /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
3822 * a nonsquare D. */
3823
3824 INLINE GEN
mkF2(ulong a,ulong b)3825 mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
3826
3827 GEN
Fl2_mul_pre(GEN x,GEN y,ulong D,ulong p,ulong pi)3828 Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
3829 {
3830 ulong xaya, xbyb, Db2, mid;
3831 ulong z1, z2;
3832 ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
3833 xaya = Fl_mul_pre(x1,y1,p,pi);
3834 if (x2==0 && y2==0) return mkF2(xaya,0);
3835 if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
3836 if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
3837 xbyb = Fl_mul_pre(x2,y2,p,pi);
3838 mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
3839 Db2 = Fl_mul_pre(D, xbyb, p,pi);
3840 z1 = Fl_add(xaya,Db2,p);
3841 z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
3842 return mkF2(z1,z2);
3843 }
3844
3845 GEN
Fl2_sqr_pre(GEN x,ulong D,ulong p,ulong pi)3846 Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
3847 {
3848 ulong a = x[1], b = x[2];
3849 ulong a2, Db2, ab;
3850 a2 = Fl_sqr_pre(a,p,pi);
3851 if (b==0) return mkF2(a2,0);
3852 Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
3853 ab = Fl_mul_pre(a,b,p,pi);
3854 return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
3855 }
3856
3857 ulong
Fl2_norm_pre(GEN x,ulong D,ulong p,ulong pi)3858 Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
3859 {
3860 ulong a2 = Fl_sqr_pre(x[1],p,pi);
3861 return x[2]? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(x[2], p,pi), p,pi), p): a2;
3862 }
3863
3864 GEN
Fl2_inv_pre(GEN x,ulong D,ulong p,ulong pi)3865 Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
3866 {
3867 ulong n, ni;
3868 if (x[2] == 0) return mkF2(Fl_inv(x[1],p),0);
3869 n = Fl_sub(Fl_sqr_pre(x[1], p,pi),
3870 Fl_mul_pre(D, Fl_sqr_pre(x[2], p,pi), p,pi), p);
3871 ni = Fl_inv(n,p);
3872 return mkF2(Fl_mul_pre(x[1], ni, p,pi),
3873 Fl_neg(Fl_mul_pre(x[2], ni, p,pi), p));
3874 }
3875
3876 int
Fl2_equal1(GEN x)3877 Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
3878
3879 struct _Fl2 {
3880 ulong p, pi, D;
3881 };
3882
3883 static GEN
_Fl2_sqr(void * data,GEN x)3884 _Fl2_sqr(void *data, GEN x)
3885 {
3886 struct _Fl2 *D = (struct _Fl2*)data;
3887 return Fl2_sqr_pre(x, D->D, D->p, D->pi);
3888 }
3889 static GEN
_Fl2_mul(void * data,GEN x,GEN y)3890 _Fl2_mul(void *data, GEN x, GEN y)
3891 {
3892 struct _Fl2 *D = (struct _Fl2*)data;
3893 return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
3894 }
3895
3896 /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
3897 GEN
Fl2_pow_pre(GEN x,GEN n,ulong D,ulong p,ulong pi)3898 Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
3899 {
3900 pari_sp av = avma;
3901 struct _Fl2 d;
3902 GEN y;
3903 long s = signe(n);
3904 if (!s) return mkF2(1,0);
3905 if (s < 0)
3906 x = Fl2_inv_pre(x,D,p,pi);
3907 if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
3908 d.p = p; d.pi = pi; d.D=D;
3909 y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
3910 return gerepileuptoleaf(av, y);
3911 }
3912
3913 static GEN
_Fl2_pow(void * data,GEN x,GEN n)3914 _Fl2_pow(void *data, GEN x, GEN n)
3915 {
3916 struct _Fl2 *D = (struct _Fl2*)data;
3917 return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
3918 }
3919
3920 static GEN
_Fl2_rand(void * data)3921 _Fl2_rand(void *data)
3922 {
3923 struct _Fl2 *D = (struct _Fl2*)data;
3924 ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
3925 return mkF2(a,b);
3926 }
3927
3928 static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
3929 hash_GEN, zv_equal, Fl2_equal1, NULL};
3930
3931 GEN
Fl2_sqrtn_pre(GEN a,GEN n,ulong D,ulong p,ulong pi,GEN * zeta)3932 Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
3933 {
3934 struct _Fl2 E;
3935 GEN o;
3936 if (a[1]==0 && a[2]==0)
3937 {
3938 if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3939 if (zeta) *zeta=mkF2(1,0);
3940 return zv_copy(a);
3941 }
3942 E.p=p; E.pi = pi; E.D = D;
3943 o = subiu(powuu(p,2), 1);
3944 return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
3945 }
3946
3947 GEN
Flx_Fl2_eval_pre(GEN x,GEN y,ulong D,ulong p,ulong pi)3948 Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
3949 {
3950 GEN p1;
3951 long i = lg(x)-1;
3952 if (i <= 2)
3953 return mkF2(i == 2? x[2]: 0, 0);
3954 p1 = mkF2(x[i], 0);
3955 for (i--; i>=2; i--)
3956 {
3957 p1 = Fl2_mul_pre(p1, y, D, p, pi);
3958 uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
3959 }
3960 return p1;
3961 }
3962
3963 /***********************************************************************/
3964 /** **/
3965 /** FlxV **/
3966 /** **/
3967 /***********************************************************************/
3968 /* FlxV are t_VEC with Flx coefficients. */
3969
3970 GEN
FlxV_Flc_mul(GEN V,GEN W,ulong p)3971 FlxV_Flc_mul(GEN V, GEN W, ulong p)
3972 {
3973 pari_sp ltop=avma;
3974 long i;
3975 GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
3976 for(i=2;i<lg(V);i++)
3977 z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
3978 return gerepileuptoleaf(ltop,z);
3979 }
3980
3981 GEN
ZXV_to_FlxV(GEN x,ulong p)3982 ZXV_to_FlxV(GEN x, ulong p)
3983 { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
3984
3985 GEN
ZXT_to_FlxT(GEN x,ulong p)3986 ZXT_to_FlxT(GEN x, ulong p)
3987 {
3988 if (typ(x) == t_POL)
3989 return ZX_to_Flx(x, p);
3990 else
3991 pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
3992 }
3993
3994 GEN
FlxV_to_Flm(GEN x,long n)3995 FlxV_to_Flm(GEN x, long n)
3996 { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
3997
3998 GEN
FlxV_red(GEN x,ulong p)3999 FlxV_red(GEN x, ulong p)
4000 { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4001
4002 GEN
FlxT_red(GEN x,ulong p)4003 FlxT_red(GEN x, ulong p)
4004 {
4005 if (typ(x) == t_VECSMALL)
4006 return Flx_red(x, p);
4007 else
4008 pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4009 }
4010
4011 GEN
FlxqV_dotproduct(GEN x,GEN y,GEN T,ulong p)4012 FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4013 {
4014 long i, lx = lg(x);
4015 pari_sp av;
4016 GEN c;
4017 if (lx == 1) return pol0_Flx(get_Flx_var(T));
4018 av = avma; c = Flx_mul(gel(x,1),gel(y,1), p);
4019 for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul(gel(x,i),gel(y,i), p), p);
4020 return gerepileuptoleaf(av, Flx_rem(c,T,p));
4021 }
4022
4023 GEN
FlxqX_dotproduct(GEN x,GEN y,GEN T,ulong p)4024 FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4025 {
4026 long i, l = minss(lg(x), lg(y));
4027 pari_sp av;
4028 GEN c;
4029 if (l == 2) return pol0_Flx(get_Flx_var(T));
4030 av = avma; c = Flx_mul(gel(x,2),gel(y,2), p);
4031 for (i=3; i<l; i++) c = Flx_add(c, Flx_mul(gel(x,i),gel(y,i), p), p);
4032 return gerepileuptoleaf(av, Flx_rem(c,T,p));
4033 }
4034
4035 GEN
FlxC_eval_powers_pre(GEN z,GEN x,ulong p,ulong pi)4036 FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4037 {
4038 long i, l = lg(z);
4039 GEN y = cgetg(l, t_VECSMALL);
4040 for (i=1; i<l; i++)
4041 uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4042 return y;
4043 }
4044
4045 /***********************************************************************/
4046 /** **/
4047 /** FlxM **/
4048 /** **/
4049 /***********************************************************************/
4050
4051 GEN
FlxM_eval_powers_pre(GEN z,GEN x,ulong p,ulong pi)4052 FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4053 {
4054 long i, l = lg(z);
4055 GEN y = cgetg(l, t_MAT);
4056 for (i=1; i<l; i++)
4057 gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4058 return y;
4059 }
4060
4061 GEN
zero_FlxC(long n,long sv)4062 zero_FlxC(long n, long sv)
4063 {
4064 long i;
4065 GEN x = cgetg(n + 1, t_COL);
4066 GEN z = zero_Flx(sv);
4067 for (i = 1; i <= n; i++)
4068 gel(x, i) = z;
4069 return x;
4070 }
4071
4072 GEN
FlxC_neg(GEN x,ulong p)4073 FlxC_neg(GEN x, ulong p)
4074 { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4075
4076 GEN
FlxC_sub(GEN x,GEN y,ulong p)4077 FlxC_sub(GEN x, GEN y, ulong p)
4078 { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4079
4080 GEN
zero_FlxM(long r,long c,long sv)4081 zero_FlxM(long r, long c, long sv)
4082 {
4083 long j;
4084 GEN x = cgetg(c + 1, t_MAT);
4085 GEN z = zero_FlxC(r, sv);
4086 for (j = 1; j <= c; j++)
4087 gel(x, j) = z;
4088 return x;
4089 }
4090
4091 GEN
FlxM_neg(GEN x,ulong p)4092 FlxM_neg(GEN x, ulong p)
4093 { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4094
4095 GEN
FlxM_sub(GEN x,GEN y,ulong p)4096 FlxM_sub(GEN x, GEN y, ulong p)
4097 { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4098
4099 GEN
FlxqC_Flxq_mul(GEN x,GEN y,GEN T,ulong p)4100 FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4101 { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4102
4103 GEN
FlxqM_Flxq_mul(GEN x,GEN y,GEN T,ulong p)4104 FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4105 { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4106
4107 static GEN
FlxM_pack_ZM(GEN M,GEN (* pack)(GEN,long))4108 FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4109 long i, j, l, lc;
4110 GEN N = cgetg_copy(M, &l), x;
4111 if (l == 1)
4112 return N;
4113 lc = lgcols(M);
4114 for (j = 1; j < l; j++) {
4115 gel(N, j) = cgetg(lc, t_COL);
4116 for (i = 1; i < lc; i++) {
4117 x = gcoeff(M, i, j);
4118 gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4119 }
4120 }
4121 return N;
4122 }
4123
4124 static GEN
kron_pack_Flx_spec_half(GEN x,long l)4125 kron_pack_Flx_spec_half(GEN x, long l) {
4126 if (l == 0)
4127 return gen_0;
4128 return Flx_to_int_halfspec(x, l);
4129 }
4130
4131 static GEN
kron_pack_Flx_spec(GEN x,long l)4132 kron_pack_Flx_spec(GEN x, long l) {
4133 long i;
4134 GEN w, y;
4135 if (l == 0)
4136 return gen_0;
4137 y = cgetipos(l + 2);
4138 for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4139 *w = x[i];
4140 return y;
4141 }
4142
4143 static GEN
kron_pack_Flx_spec_2(GEN x,long l)4144 kron_pack_Flx_spec_2(GEN x, long l) {
4145 return Flx_eval2BILspec(x, 2, l);
4146 }
4147
4148 static GEN
kron_pack_Flx_spec_3(GEN x,long l)4149 kron_pack_Flx_spec_3(GEN x, long l) {
4150 return Flx_eval2BILspec(x, 3, l);
4151 }
4152
4153 static GEN
kron_unpack_Flx(GEN z,ulong p)4154 kron_unpack_Flx(GEN z, ulong p)
4155 {
4156 long i, l = lgefint(z);
4157 GEN x = cgetg(l, t_VECSMALL), w;
4158 for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4159 x[i] = ((ulong) *w) % p;
4160 return Flx_renormalize(x, l);
4161 }
4162
4163 static GEN
kron_unpack_Flx_2(GEN x,ulong p)4164 kron_unpack_Flx_2(GEN x, ulong p) {
4165 long d = (lgefint(x)-1)/2 - 1;
4166 return Z_mod2BIL_Flx_2(x, d, p);
4167 }
4168
4169 static GEN
kron_unpack_Flx_3(GEN x,ulong p)4170 kron_unpack_Flx_3(GEN x, ulong p) {
4171 long d = lgefint(x)/3 - 1;
4172 return Z_mod2BIL_Flx_3(x, d, p);
4173 }
4174
4175 static GEN
FlxM_pack_ZM_bits(GEN M,long b)4176 FlxM_pack_ZM_bits(GEN M, long b)
4177 {
4178 long i, j, l, lc;
4179 GEN N = cgetg_copy(M, &l), x;
4180 if (l == 1)
4181 return N;
4182 lc = lgcols(M);
4183 for (j = 1; j < l; j++) {
4184 gel(N, j) = cgetg(lc, t_COL);
4185 for (i = 1; i < lc; i++) {
4186 x = gcoeff(M, i, j);
4187 gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
4188 }
4189 }
4190 return N;
4191 }
4192
4193 static GEN
ZM_unpack_FlxqM(GEN M,GEN T,ulong p,GEN (* unpack)(GEN,ulong))4194 ZM_unpack_FlxqM(GEN M, GEN T, ulong p, GEN (*unpack)(GEN, ulong))
4195 {
4196 long i, j, l, lc, sv = get_Flx_var(T);
4197 GEN N = cgetg_copy(M, &l), x;
4198 if (l == 1)
4199 return N;
4200 lc = lgcols(M);
4201 for (j = 1; j < l; j++) {
4202 gel(N, j) = cgetg(lc, t_COL);
4203 for (i = 1; i < lc; i++) {
4204 x = unpack(gcoeff(M, i, j), p);
4205 x[1] = sv;
4206 gcoeff(N, i, j) = Flx_rem(x, T, p);
4207 }
4208 }
4209 return N;
4210 }
4211
4212 static GEN
ZM_unpack_FlxqM_bits(GEN M,long b,GEN T,ulong p)4213 ZM_unpack_FlxqM_bits(GEN M, long b, GEN T, ulong p)
4214 {
4215 long i, j, l, lc, sv = get_Flx_var(T);
4216 GEN N = cgetg_copy(M, &l), x;
4217 if (l == 1)
4218 return N;
4219 lc = lgcols(M);
4220 if (b < BITS_IN_LONG) {
4221 for (j = 1; j < l; j++) {
4222 gel(N, j) = cgetg(lc, t_COL);
4223 for (i = 1; i < lc; i++) {
4224 x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
4225 x[1] = sv;
4226 gcoeff(N, i, j) = Flx_rem(x, T, p);
4227 }
4228 }
4229 } else {
4230 ulong pi = get_Fl_red(p);
4231 for (j = 1; j < l; j++) {
4232 gel(N, j) = cgetg(lc, t_COL);
4233 for (i = 1; i < lc; i++) {
4234 x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
4235 x[1] = sv;
4236 gcoeff(N, i, j) = Flx_rem(x, T, p);
4237 }
4238 }
4239 }
4240 return N;
4241 }
4242
4243 GEN
FlxqM_mul_Kronecker(GEN A,GEN B,GEN T,ulong p)4244 FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
4245 {
4246 pari_sp av = avma;
4247 long b, d = get_Flx_degree(T), n = lg(A) - 1;
4248 GEN C, D, z;
4249 GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
4250 int is_sqr = A==B;
4251
4252 z = muliu(muliu(sqru(p - 1), d), n);
4253 b = expi(z) + 1;
4254 /* only do expensive bit-packing if it saves at least 1 limb */
4255 if (b <= BITS_IN_HALFULONG) {
4256 if (nbits2nlong(d*b) == (d + 1)/2)
4257 b = BITS_IN_HALFULONG;
4258 } else {
4259 long l = lgefint(z) - 2;
4260 if (nbits2nlong(d*b) == d*l)
4261 b = l*BITS_IN_LONG;
4262 }
4263 set_avma(av);
4264
4265 switch (b) {
4266 case BITS_IN_HALFULONG:
4267 pack = kron_pack_Flx_spec_half;
4268 unpack = int_to_Flx_half;
4269 break;
4270 case BITS_IN_LONG:
4271 pack = kron_pack_Flx_spec;
4272 unpack = kron_unpack_Flx;
4273 break;
4274 case 2*BITS_IN_LONG:
4275 pack = kron_pack_Flx_spec_2;
4276 unpack = kron_unpack_Flx_2;
4277 break;
4278 case 3*BITS_IN_LONG:
4279 pack = kron_pack_Flx_spec_3;
4280 unpack = kron_unpack_Flx_3;
4281 break;
4282 default:
4283 A = FlxM_pack_ZM_bits(A, b);
4284 B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
4285 C = ZM_mul(A, B);
4286 D = ZM_unpack_FlxqM_bits(C, b, T, p);
4287 return gerepilecopy(av, D);
4288 }
4289 A = FlxM_pack_ZM(A, pack);
4290 B = is_sqr? A: FlxM_pack_ZM(B, pack);
4291 C = ZM_mul(A, B);
4292 D = ZM_unpack_FlxqM(C, T, p, unpack);
4293 return gerepilecopy(av, D);
4294 }
4295