Line data Source code
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
21 886789313 : get_Flx_red(GEN T, GEN *B)
22 : {
23 886789313 : if (typ(T)!=t_VEC) { *B=NULL; return T; }
24 586202 : *B = gel(T,1); return gel(T,2);
25 : }
26 :
27 : /***********************************************************************/
28 : /** Flx **/
29 : /***********************************************************************/
30 : /* Flx objects are defined as follows:
31 : * Let l an ulong. An Flx is a t_VECSMALL:
32 : * x[0] = codeword
33 : * x[1] = evalvarn(variable number) (signe is not stored).
34 : * x[2] = a_0 x[3] = a_1, etc. with 0 <= a_i < l
35 : *
36 : * signe(x) is not valid. Use degpol(x)>0 instead. */
37 : /***********************************************************************/
38 : /** Conversion from Flx **/
39 : /***********************************************************************/
40 :
41 : GEN
42 35611914 : Flx_to_ZX(GEN z)
43 : {
44 35611914 : long i, l = lg(z);
45 35611914 : GEN x = cgetg(l,t_POL);
46 236380012 : for (i=2; i<l; i++) gel(x,i) = utoi(z[i]);
47 35575944 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
48 : }
49 :
50 : GEN
51 68638 : Flx_to_FlxX(GEN z, long sv)
52 : {
53 68638 : long i, l = lg(z);
54 68638 : GEN x = cgetg(l,t_POL);
55 269873 : for (i=2; i<l; i++) gel(x,i) = Fl_to_Flx(z[i], sv);
56 68638 : x[1] = evalsigne(l-2!=0)| z[1]; return x;
57 : }
58 :
59 : /* same as Flx_to_ZX, in place */
60 : GEN
61 36120660 : Flx_to_ZX_inplace(GEN z)
62 : {
63 36120660 : long i, l = lg(z);
64 228057490 : for (i=2; i<l; i++) gel(z,i) = utoi(z[i]);
65 36071709 : settyp(z, t_POL); z[1]=evalsigne(l-2!=0)|z[1]; return z;
66 : }
67 :
68 : /*Flx_to_Flv=zx_to_zv*/
69 : GEN
70 61853797 : Flx_to_Flv(GEN x, long N)
71 : {
72 61853797 : GEN z = cgetg(N+1,t_VECSMALL);
73 61840402 : long i, l = lg(x)-1;
74 61840402 : x++;
75 698958221 : for (i=1; i<l ; i++) z[i]=x[i];
76 324510218 : for ( ; i<=N; i++) z[i]=0;
77 61840402 : return z;
78 : }
79 :
80 : /*Flv_to_Flx=zv_to_zx*/
81 : GEN
82 26166059 : Flv_to_Flx(GEN x, long sv)
83 : {
84 26166059 : long i, l=lg(x)+1;
85 26166059 : GEN z = cgetg(l,t_VECSMALL); z[1]=sv;
86 26154444 : x--;
87 291539492 : for (i=2; i<l ; i++) z[i]=x[i];
88 26154444 : return Flx_renormalize(z,l);
89 : }
90 :
91 : /*Flm_to_FlxV=zm_to_zxV*/
92 : GEN
93 2268 : Flm_to_FlxV(GEN x, long sv)
94 6153 : { pari_APPLY_type(t_VEC, Flv_to_Flx(gel(x,i), sv)) }
95 :
96 : /*FlxC_to_ZXC=zxC_to_ZXC*/
97 : GEN
98 106341 : FlxC_to_ZXC(GEN x)
99 555431 : { pari_APPLY_type(t_COL, Flx_to_ZX(gel(x,i))) }
100 :
101 : /*FlxC_to_ZXC=zxV_to_ZXV*/
102 : GEN
103 591763 : FlxV_to_ZXV(GEN x)
104 2399268 : { pari_APPLY_type(t_VEC, Flx_to_ZX(gel(x,i))) }
105 :
106 : void
107 2765537 : FlxV_to_ZXV_inplace(GEN v)
108 : {
109 : long i;
110 7390833 : for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
111 2765467 : }
112 :
113 : /*FlxM_to_ZXM=zxM_to_ZXM*/
114 : GEN
115 4515 : FlxM_to_ZXM(GEN x)
116 15734 : { pari_APPLY_same(FlxC_to_ZXC(gel(x,i))) }
117 :
118 : GEN
119 378439 : FlxV_to_FlxX(GEN x, long v)
120 : {
121 378439 : long i, l = lg(x)+1;
122 378439 : GEN z = cgetg(l,t_POL); z[1] = evalvarn(v);
123 378439 : x--;
124 4911174 : for (i=2; i<l ; i++) gel(z,i) = gel(x,i);
125 378439 : return FlxX_renormalize(z,l);
126 : }
127 :
128 : GEN
129 0 : FlxM_to_FlxXV(GEN x, long v)
130 0 : { pari_APPLY_type(t_COL, FlxV_to_FlxX(gel(x,i), v)) }
131 :
132 : GEN
133 0 : FlxM_Flx_add_shallow(GEN x, GEN y, ulong p)
134 : {
135 0 : long l = lg(x), i, j;
136 0 : GEN z = cgetg(l,t_MAT);
137 :
138 0 : if (l==1) return z;
139 0 : if (l != lgcols(x)) pari_err_OP( "+", x, y);
140 0 : for (i=1; i<l; i++)
141 : {
142 0 : GEN zi = cgetg(l,t_COL), xi = gel(x,i);
143 0 : gel(z,i) = zi;
144 0 : for (j=1; j<l; j++) gel(zi,j) = gel(xi,j);
145 0 : gel(zi,i) = Flx_add(gel(zi,i), y, p);
146 : }
147 0 : return z;
148 : }
149 :
150 : /***********************************************************************/
151 : /** Conversion to Flx **/
152 : /***********************************************************************/
153 : /* Take an integer and return a scalar polynomial mod p, with evalvarn=vs */
154 : GEN
155 15941153 : Fl_to_Flx(ulong x, long sv) { return x? mkvecsmall2(sv, x): pol0_Flx(sv); }
156 :
157 : /* a X^d */
158 : GEN
159 785808 : monomial_Flx(ulong a, long d, long vs)
160 : {
161 : GEN P;
162 785808 : if (a==0) return pol0_Flx(vs);
163 785808 : P = const_vecsmall(d+2, 0);
164 785823 : P[1] = vs; P[d+2] = a; return P;
165 : }
166 :
167 : GEN
168 2576749 : Z_to_Flx(GEN x, ulong p, long sv)
169 : {
170 2576749 : long u = umodiu(x,p);
171 2576880 : return u? mkvecsmall2(sv, u): pol0_Flx(sv);
172 : }
173 :
174 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
175 : GEN
176 156500466 : ZX_to_Flx(GEN x, ulong p)
177 : {
178 156500466 : long i, lx = lg(x);
179 156500466 : GEN a = cgetg(lx, t_VECSMALL);
180 156278154 : a[1]=((ulong)x[1])&VARNBITS;
181 1061363943 : for (i=2; i<lx; i++) a[i] = umodiu(gel(x,i), p);
182 156489476 : return Flx_renormalize(a,lx);
183 : }
184 :
185 : /* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL*/
186 : GEN
187 5234563 : zx_to_Flx(GEN x, ulong p)
188 : {
189 5234563 : long i, lx = lg(x);
190 5234563 : GEN a = cgetg(lx, t_VECSMALL);
191 5225521 : a[1] = x[1];
192 16168029 : for (i=2; i<lx; i++) uel(a,i) = umodsu(x[i], p);
193 5223626 : return Flx_renormalize(a,lx);
194 : }
195 :
196 : ulong
197 57436189 : Rg_to_Fl(GEN x, ulong p)
198 : {
199 57436189 : switch(typ(x))
200 : {
201 28830589 : case t_INT: return umodiu(x, p);
202 330451 : case t_FRAC: {
203 330451 : ulong z = umodiu(gel(x,1), p);
204 330453 : if (!z) return 0;
205 326152 : return Fl_div(z, umodiu(gel(x,2), p), p);
206 : }
207 205957 : case t_PADIC: return padic_to_Fl(x, p);
208 28069209 : case t_INTMOD: {
209 28069209 : GEN q = gel(x,1), a = gel(x,2);
210 28069209 : if (absequaliu(q, p)) return itou(a);
211 0 : if (!dvdiu(q,p)) pari_err_MODULUS("Rg_to_Fl", q, utoipos(p));
212 0 : return umodiu(a, p);
213 : }
214 0 : default: pari_err_TYPE("Rg_to_Fl",x);
215 : return 0; /* LCOV_EXCL_LINE */
216 : }
217 : }
218 :
219 : ulong
220 1665346 : Rg_to_F2(GEN x)
221 : {
222 1665346 : switch(typ(x))
223 : {
224 262330 : case t_INT: return mpodd(x);
225 0 : case t_FRAC:
226 0 : if (!mpodd(gel(x,2))) (void)Fl_inv(0,2); /* error */
227 0 : return mpodd(gel(x,1));
228 0 : case t_PADIC:
229 0 : if (!absequaliu(gel(x,2),2)) pari_err_OP("",x, mkintmodu(1,2));
230 0 : if (valp(x) < 0) (void)Fl_inv(0,2);
231 0 : return valp(x) & 1;
232 1403017 : case t_INTMOD: {
233 1403017 : GEN q = gel(x,1), a = gel(x,2);
234 1403017 : if (mpodd(q)) pari_err_MODULUS("Rg_to_F2", q, gen_2);
235 1403017 : return mpodd(a);
236 : }
237 0 : default: pari_err_TYPE("Rg_to_F2",x);
238 : return 0; /* LCOV_EXCL_LINE */
239 : }
240 : }
241 :
242 : GEN
243 2700471 : RgX_to_Flx(GEN x, ulong p)
244 : {
245 2700471 : long i, lx = lg(x);
246 2700471 : GEN a = cgetg(lx, t_VECSMALL);
247 2700471 : a[1]=((ulong)x[1])&VARNBITS;
248 23696100 : for (i=2; i<lx; i++) a[i] = Rg_to_Fl(gel(x,i), p);
249 2700471 : return Flx_renormalize(a,lx);
250 : }
251 :
252 : GEN
253 7 : RgXV_to_FlxV(GEN x, ulong p)
254 175 : { pari_APPLY_type(t_VEC, RgX_to_Flx(gel(x,i), p)) }
255 :
256 : /* If x is a POLMOD, assume modulus is a multiple of T. */
257 : GEN
258 3328545 : Rg_to_Flxq(GEN x, GEN T, ulong p)
259 : {
260 3328545 : long ta, tx = typ(x), v = get_Flx_var(T);
261 : ulong pi;
262 : GEN a, b;
263 3328548 : if (is_const_t(tx))
264 : {
265 3168827 : if (tx == t_FFELT) return FF_to_Flxq(x);
266 2437917 : return Fl_to_Flx(Rg_to_Fl(x, p), v);
267 : }
268 159717 : switch(tx)
269 : {
270 8576 : case t_POLMOD:
271 8576 : b = gel(x,1);
272 8576 : a = gel(x,2); ta = typ(a);
273 8576 : if (is_const_t(ta)) return Fl_to_Flx(Rg_to_Fl(a, p), v);
274 8422 : b = RgX_to_Flx(b, p); if (b[1] != v) break;
275 8422 : a = RgX_to_Flx(a, p); if (Flx_equal(b,T)) return a;
276 0 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
277 0 : if (lgpol(Flx_rem_pre(b,T,p,pi))==0) return Flx_rem_pre(a, T, p, pi);
278 0 : break;
279 151141 : case t_POL:
280 151141 : x = RgX_to_Flx(x,p);
281 151141 : if (x[1] != v) break;
282 151141 : return Flx_rem(x, T, p);
283 0 : case t_RFRAC:
284 0 : a = Rg_to_Flxq(gel(x,1), T,p);
285 0 : b = Rg_to_Flxq(gel(x,2), T,p);
286 0 : return Flxq_div(a,b, T,p);
287 : }
288 0 : pari_err_TYPE("Rg_to_Flxq",x);
289 : return NULL; /* LCOV_EXCL_LINE */
290 : }
291 :
292 : /***********************************************************************/
293 : /** Basic operation on Flx **/
294 : /***********************************************************************/
295 : /* = zx_renormalize. Similar to normalizepol, in place */
296 : GEN
297 1971353147 : Flx_renormalize(GEN /*in place*/ x, long lx)
298 : {
299 : long i;
300 2220776163 : for (i = lx-1; i>1; i--)
301 2129333321 : if (x[i]) break;
302 1971353147 : stackdummy((pari_sp)(x + lg(x)), (pari_sp)(x + i+1));
303 1969137634 : setlg(x, i+1); return x;
304 : }
305 :
306 : GEN
307 1884426 : Flx_red(GEN z, ulong p)
308 : {
309 1884426 : long i, l = lg(z);
310 1884426 : GEN x = cgetg(l, t_VECSMALL);
311 1884324 : x[1] = z[1];
312 34749873 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
313 1884324 : return Flx_renormalize(x,l);
314 : }
315 :
316 : int
317 27972407 : Flx_equal(GEN V, GEN W)
318 : {
319 27972407 : long l = lg(V);
320 27972407 : if (lg(W) != l) return 0;
321 28993161 : while (--l > 1) /* do not compare variables, V[1] */
322 27736117 : if (V[l] != W[l]) return 0;
323 1257044 : return 1;
324 : }
325 :
326 : GEN
327 2447597 : random_Flx(long d1, long vs, ulong p)
328 : {
329 2447597 : long i, d = d1+2;
330 2447597 : GEN y = cgetg(d,t_VECSMALL); y[1] = vs;
331 17303144 : for (i=2; i<d; i++) y[i] = random_Fl(p);
332 2447944 : return Flx_renormalize(y,d);
333 : }
334 :
335 : static GEN
336 5986746 : Flx_addspec(GEN x, GEN y, ulong p, long lx, long ly)
337 : {
338 : long i,lz;
339 : GEN z;
340 :
341 5986746 : if (ly>lx) swapspec(x,y, lx,ly);
342 5986746 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
343 101048052 : for (i=0; i<ly; i++) z[i+2] = Fl_add(x[i], y[i], p);
344 85389080 : for ( ; i<lx; i++) z[i+2] = x[i];
345 5986746 : z[1] = 0; return Flx_renormalize(z, lz);
346 : }
347 :
348 : GEN
349 56227355 : Flx_add(GEN x, GEN y, ulong p)
350 : {
351 : long i,lz;
352 : GEN z;
353 56227355 : long lx=lg(x);
354 56227355 : long ly=lg(y);
355 56227355 : if (ly>lx) swapspec(x,y, lx,ly);
356 56227355 : lz = lx; z = cgetg(lz, t_VECSMALL); z[1]=x[1];
357 530873383 : for (i=2; i<ly; i++) z[i] = Fl_add(x[i], y[i], p);
358 116250428 : for ( ; i<lx; i++) z[i] = x[i];
359 56175971 : return Flx_renormalize(z, lz);
360 : }
361 :
362 : GEN
363 9259468 : Flx_Fl_add(GEN y, ulong x, ulong p)
364 : {
365 : GEN z;
366 : long lz, i;
367 9259468 : if (!lgpol(y))
368 225951 : return Fl_to_Flx(x,y[1]);
369 9034424 : lz=lg(y);
370 9034424 : z=cgetg(lz,t_VECSMALL);
371 9032911 : z[1]=y[1];
372 9032911 : z[2] = Fl_add(y[2],x,p);
373 45057838 : for(i=3;i<lz;i++)
374 36024996 : z[i] = y[i];
375 9032842 : if (lz==3) z = Flx_renormalize(z,lz);
376 9032881 : return z;
377 : }
378 :
379 : static GEN
380 897414 : Flx_subspec(GEN x, GEN y, ulong p, long lx, long ly)
381 : {
382 : long i,lz;
383 : GEN z;
384 :
385 897414 : if (ly <= lx)
386 : {
387 897459 : lz = lx+2; z = cgetg(lz, t_VECSMALL);
388 54194716 : for (i=0; i<ly; i++) z[i+2] = Fl_sub(x[i],y[i],p);
389 1442850 : for ( ; i<lx; i++) z[i+2] = x[i];
390 : }
391 : else
392 : {
393 0 : lz = ly+2; z = cgetg(lz, t_VECSMALL);
394 0 : for (i=0; i<lx; i++) z[i+2] = Fl_sub(x[i],y[i],p);
395 0 : for ( ; i<ly; i++) z[i+2] = Fl_neg(y[i],p);
396 : }
397 896613 : z[1] = 0; return Flx_renormalize(z, lz);
398 : }
399 :
400 : GEN
401 135083766 : Flx_sub(GEN x, GEN y, ulong p)
402 : {
403 135083766 : long i,lz,lx = lg(x), ly = lg(y);
404 : GEN z;
405 :
406 135083766 : if (ly <= lx)
407 : {
408 87108240 : lz = lx; z = cgetg(lz, t_VECSMALL);
409 466934465 : for (i=2; i<ly; i++) z[i] = Fl_sub(x[i],y[i],p);
410 174577712 : for ( ; i<lx; i++) z[i] = x[i];
411 : }
412 : else
413 : {
414 47975526 : lz = ly; z = cgetg(lz, t_VECSMALL);
415 328012067 : for (i=2; i<lx; i++) z[i] = Fl_sub(x[i],y[i],p);
416 236062286 : for ( ; i<ly; i++) z[i] = y[i]? (long)(p - y[i]): y[i];
417 : }
418 135012468 : z[1]=x[1]; return Flx_renormalize(z, lz);
419 : }
420 :
421 : GEN
422 32918 : Flx_Fl_sub(GEN y, ulong x, ulong p)
423 : {
424 : GEN z;
425 32918 : long lz = lg(y), i;
426 32918 : if (lz==2)
427 506 : return Fl_to_Flx(Fl_neg(x, p),y[1]);
428 32412 : z = cgetg(lz, t_VECSMALL);
429 32412 : z[1] = y[1];
430 32412 : z[2] = Fl_sub(uel(y,2), x, p);
431 208798 : for(i=3; i<lz; i++)
432 176386 : z[i] = y[i];
433 32412 : if (lz==3) z = Flx_renormalize(z,lz);
434 32412 : return z;
435 : }
436 :
437 : static GEN
438 3331002 : Flx_negspec(GEN x, ulong p, long l)
439 : {
440 : long i;
441 3331002 : GEN z = cgetg(l+2, t_VECSMALL) + 2;
442 20901101 : for (i=0; i<l; i++) z[i] = Fl_neg(x[i], p);
443 3331025 : return z-2;
444 : }
445 :
446 : GEN
447 3331030 : Flx_neg(GEN x, ulong p)
448 : {
449 3331030 : GEN z = Flx_negspec(x+2, p, lgpol(x));
450 3331332 : z[1] = x[1];
451 3331332 : return z;
452 : }
453 :
454 : GEN
455 1463215 : Flx_neg_inplace(GEN x, ulong p)
456 : {
457 1463215 : long i, l = lg(x);
458 49859636 : for (i=2; i<l; i++)
459 48396421 : if (x[i]) x[i] = p - x[i];
460 1463215 : return x;
461 : }
462 :
463 : GEN
464 2163730 : Flx_double(GEN y, ulong p)
465 : {
466 : long i, l;
467 2163730 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
468 20016322 : for(i=2; i<l; i++) z[i] = Fl_double(y[i], p);
469 2163730 : return Flx_renormalize(z, l);
470 : }
471 : GEN
472 901988 : Flx_triple(GEN y, ulong p)
473 : {
474 : long i, l;
475 901988 : GEN z = cgetg_copy(y, &l); z[1] = y[1];
476 7972336 : for(i=2; i<l; i++) z[i] = Fl_triple(y[i], p);
477 901988 : return Flx_renormalize(z, l);
478 : }
479 : GEN
480 17424015 : Flx_Fl_mul(GEN y, ulong x, ulong p)
481 : {
482 : GEN z;
483 : long i, l;
484 17424015 : if (!x) return pol0_Flx(y[1]);
485 16737504 : z = cgetg_copy(y, &l); z[1] = y[1];
486 16736515 : if (HIGHWORD(x | p))
487 14346508 : for(i=2; i<l; i++) z[i] = Fl_mul(y[i], x, p);
488 : else
489 89294583 : for(i=2; i<l; i++) z[i] = (y[i] * x) % p;
490 16736838 : return Flx_renormalize(z, l);
491 : }
492 : GEN
493 11492990 : Flx_Fl_mul_to_monic(GEN y, ulong x, ulong p)
494 : {
495 : GEN z;
496 : long i, l;
497 11492990 : z = cgetg_copy(y, &l); z[1] = y[1];
498 11487933 : if (HIGHWORD(x | p))
499 5308324 : for(i=2; i<l-1; i++) z[i] = Fl_mul(y[i], x, p);
500 : else
501 25551584 : for(i=2; i<l-1; i++) z[i] = (y[i] * x) % p;
502 11487918 : z[l-1] = 1; return z;
503 : }
504 :
505 : /* Return a*x^n if n>=0 and a\x^(-n) if n<0 */
506 : GEN
507 24388681 : Flx_shift(GEN a, long n)
508 : {
509 24388681 : long i, l = lg(a);
510 : GEN b;
511 24388681 : if (l==2 || !n) return Flx_copy(a);
512 24047531 : if (l+n<=2) return pol0_Flx(a[1]);
513 23834336 : b = cgetg(l+n, t_VECSMALL);
514 23830749 : b[1] = a[1];
515 23830749 : if (n < 0)
516 65695219 : for (i=2-n; i<l; i++) b[i+n] = a[i];
517 : else
518 : {
519 43063311 : for (i=0; i<n; i++) b[2+i] = 0;
520 140350290 : for (i=2; i<l; i++) b[i+n] = a[i];
521 : }
522 23830749 : return b;
523 : }
524 :
525 : GEN
526 59254171 : Flx_normalize(GEN z, ulong p)
527 : {
528 59254171 : long l = lg(z)-1;
529 59254171 : ulong p1 = z[l]; /* leading term */
530 59254171 : if (p1 == 1) return z;
531 11471017 : return Flx_Fl_mul_to_monic(z, Fl_inv(p1,p), p);
532 : }
533 :
534 : /* return (x * X^d) + y. Assume d > 0, shallow if x == 0*/
535 : static GEN
536 3090830 : Flx_addshift(GEN x, GEN y, ulong p, long d)
537 : {
538 3090830 : GEN xd,yd,zd = (GEN)avma;
539 3090830 : long a,lz,ny = lgpol(y), nx = lgpol(x);
540 3090830 : long vs = x[1];
541 3090830 : if (nx == 0) return y;
542 3088978 : x += 2; y += 2; a = ny-d;
543 3088978 : if (a <= 0)
544 : {
545 85040 : lz = (a>nx)? ny+2: nx+d+2;
546 85040 : (void)new_chunk(lz); xd = x+nx; yd = y+ny;
547 1729752 : while (xd > x) *--zd = *--xd;
548 85040 : x = zd + a;
549 166127 : while (zd > x) *--zd = 0;
550 : }
551 : else
552 : {
553 3003938 : xd = new_chunk(d); yd = y+d;
554 3003938 : x = Flx_addspec(x,yd,p, nx,a);
555 3003938 : lz = (a>nx)? ny+2: lg(x)+d;
556 126606614 : x += 2; while (xd > x) *--zd = *--xd;
557 : }
558 57245727 : while (yd > y) *--zd = *--yd;
559 3088978 : *--zd = vs;
560 3088978 : *--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd;
561 : }
562 :
563 : /* shift polynomial + gerepile */
564 : /* Do not set evalvarn*/
565 : static GEN
566 601363555 : Flx_shiftip(pari_sp av, GEN x, long v)
567 : {
568 601363555 : long i, lx = lg(x), ly;
569 : GEN y;
570 601363555 : if (!v || lx==2) return gerepileuptoleaf(av, x);
571 169926031 : ly = lx + v; /* result length */
572 169926031 : (void)new_chunk(ly); /* check that result fits */
573 169772153 : x += lx; y = (GEN)av;
574 1240717159 : for (i = 2; i<lx; i++) *--y = *--x;
575 693453562 : for (i = 0; i< v; i++) *--y = 0;
576 169772153 : y -= 2; y[0] = evaltyp(t_VECSMALL) | evallg(ly);
577 170154654 : return gc_const((pari_sp)y, y);
578 : }
579 :
580 : static long
581 2145269889 : get_Fl_threshold(ulong p, long mul, long mul2)
582 : {
583 2145269889 : return SMALL_ULONG(p) ? mul: mul2;
584 : }
585 :
586 : #define BITS_IN_QUARTULONG (BITS_IN_HALFULONG >> 1)
587 : #define QUARTMASK ((1UL<<BITS_IN_QUARTULONG)-1UL)
588 : #define LLQUARTWORD(x) ((x) & QUARTMASK)
589 : #define HLQUARTWORD(x) (((x) >> BITS_IN_QUARTULONG) & QUARTMASK)
590 : #define LHQUARTWORD(x) (((x) >> (2*BITS_IN_QUARTULONG)) & QUARTMASK)
591 : #define HHQUARTWORD(x) (((x) >> (3*BITS_IN_QUARTULONG)) & QUARTMASK)
592 : INLINE long
593 7563037 : maxbitcoeffpol(ulong p, long n)
594 : {
595 7563037 : GEN z = muliu(sqru(p - 1), n);
596 7557992 : long b = expi(z) + 1;
597 : /* only do expensive bit-packing if it saves at least 1 limb */
598 7560789 : if (b <= BITS_IN_QUARTULONG)
599 : {
600 887044 : if (nbits2nlong(n*b) == (n + 3)>>2)
601 106928 : b = BITS_IN_QUARTULONG;
602 : }
603 6673745 : else if (b <= BITS_IN_HALFULONG)
604 : {
605 1566245 : if (nbits2nlong(n*b) == (n + 1)>>1)
606 5584 : b = BITS_IN_HALFULONG;
607 : }
608 : else
609 : {
610 5107500 : long l = lgefint(z) - 2;
611 5107500 : if (nbits2nlong(n*b) == n*l)
612 306027 : b = l*BITS_IN_LONG;
613 : }
614 7560796 : return b;
615 : }
616 :
617 : INLINE ulong
618 3465328292 : Flx_mullimb_ok(GEN x, GEN y, ulong p, long a, long b)
619 : { /* Assume OK_ULONG*/
620 3465328292 : ulong p1 = 0;
621 : long i;
622 16428188279 : for (i=a; i<b; i++)
623 12962859987 : if (y[i])
624 : {
625 10882074586 : p1 += y[i] * x[-i];
626 10882074586 : if (p1 & HIGHBIT) p1 %= p;
627 : }
628 3465328292 : return p1 % p;
629 : }
630 :
631 : INLINE ulong
632 1063714445 : Flx_mullimb(GEN x, GEN y, ulong p, ulong pi, long a, long b)
633 : {
634 1063714445 : ulong p1 = 0;
635 : long i;
636 3322670461 : for (i=a; i<b; i++)
637 2257317100 : if (y[i])
638 2215896772 : p1 = Fl_addmul_pre(p1, y[i], x[-i], p, pi);
639 1065353361 : return p1;
640 : }
641 :
642 : /* assume nx >= ny > 0 */
643 : static GEN
644 322243247 : Flx_mulspec_basecase(GEN x, GEN y, ulong p, ulong pi, long nx, long ny)
645 : {
646 : long i,lz,nz;
647 : GEN z;
648 :
649 322243247 : lz = nx+ny+1; nz = lz-2;
650 322243247 : z = cgetg(lz, t_VECSMALL) + 2; /* x:y:z [i] = term of degree i */
651 321745151 : if (!pi)
652 : {
653 1150402166 : for (i=0; i<ny; i++)z[i] = Flx_mullimb_ok(x+i,y,p,0,i+1);
654 772723656 : for ( ; i<nx; i++) z[i] = Flx_mullimb_ok(x+i,y,p,0,ny);
655 903498133 : for ( ; i<nz; i++) z[i] = Flx_mullimb_ok(x+i,y,p,i-nx+1,ny);
656 : }
657 : else
658 : {
659 263792675 : for (i=0; i<ny; i++)z[i] = Flx_mullimb(x+i,y,p,pi,0,i+1);
660 197340356 : for ( ; i<nx; i++) z[i] = Flx_mullimb(x+i,y,p,pi,0,ny);
661 191257131 : for ( ; i<nz; i++) z[i] = Flx_mullimb(x+i,y,p,pi,i-nx+1,ny);
662 : }
663 322121431 : z -= 2; return Flx_renormalize(z, lz);
664 : }
665 :
666 : static GEN
667 12567 : int_to_Flx(GEN z, ulong p)
668 : {
669 12567 : long i, l = lgefint(z);
670 12567 : GEN x = cgetg(l, t_VECSMALL);
671 1083877 : for (i=2; i<l; i++) x[i] = uel(z,i)%p;
672 12562 : return Flx_renormalize(x, l);
673 : }
674 :
675 : INLINE GEN
676 10298 : Flx_mulspec_mulii(GEN a, GEN b, ulong p, long na, long nb)
677 : {
678 10298 : GEN z=muliispec(a,b,na,nb);
679 10302 : return int_to_Flx(z,p);
680 : }
681 :
682 : static GEN
683 491319 : Flx_to_int_halfspec(GEN a, long na)
684 : {
685 : long j;
686 491319 : long n = (na+1)>>1UL;
687 491319 : GEN V = cgetipos(2+n);
688 : GEN w;
689 1402700 : for (w = int_LSW(V), j=0; j+1<na; j+=2, w=int_nextW(w))
690 911381 : *w = a[j]|(a[j+1]<<BITS_IN_HALFULONG);
691 491319 : if (j<na)
692 337962 : *w = a[j];
693 491319 : return V;
694 : }
695 :
696 : static GEN
697 589294 : int_to_Flx_half(GEN z, ulong p)
698 : {
699 : long i;
700 589294 : long lx = (lgefint(z)-2)*2+2;
701 589294 : GEN w, x = cgetg(lx, t_VECSMALL);
702 2004691 : for (w = int_LSW(z), i=2; i<lx; i+=2, w=int_nextW(w))
703 : {
704 1415397 : x[i] = LOWWORD((ulong)*w)%p;
705 1415397 : x[i+1] = HIGHWORD((ulong)*w)%p;
706 : }
707 589294 : return Flx_renormalize(x, lx);
708 : }
709 :
710 : static GEN
711 5448 : Flx_mulspec_halfmulii(GEN a, GEN b, ulong p, long na, long nb)
712 : {
713 5448 : GEN A = Flx_to_int_halfspec(a,na);
714 5448 : GEN B = Flx_to_int_halfspec(b,nb);
715 5448 : GEN z = mulii(A,B);
716 5448 : return int_to_Flx_half(z,p);
717 : }
718 :
719 : static GEN
720 203661 : Flx_to_int_quartspec(GEN a, long na)
721 : {
722 : long j;
723 203661 : long n = (na+3)>>2UL;
724 203661 : GEN V = cgetipos(2+n);
725 : GEN w;
726 4348964 : for (w = int_LSW(V), j=0; j+3<na; j+=4, w=int_nextW(w))
727 4145303 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG))|(a[j+3]<<(3*BITS_IN_QUARTULONG));
728 203661 : switch (na-j)
729 : {
730 116094 : case 3:
731 116094 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG)|(a[j+2]<<(2*BITS_IN_QUARTULONG));
732 116094 : break;
733 34240 : case 2:
734 34240 : *w = a[j]|(a[j+1]<<BITS_IN_QUARTULONG);
735 34240 : break;
736 26973 : case 1:
737 26973 : *w = a[j];
738 26973 : break;
739 26354 : case 0:
740 26354 : break;
741 : }
742 203661 : return V;
743 : }
744 :
745 : static GEN
746 106930 : int_to_Flx_quart(GEN z, ulong p)
747 : {
748 : long i;
749 106930 : long lx = (lgefint(z)-2)*4+2;
750 106930 : GEN w, x = cgetg(lx, t_VECSMALL);
751 4844610 : for (w = int_LSW(z), i=2; i<lx; i+=4, w=int_nextW(w))
752 : {
753 4737680 : x[i] = LLQUARTWORD((ulong)*w)%p;
754 4737680 : x[i+1] = HLQUARTWORD((ulong)*w)%p;
755 4737680 : x[i+2] = LHQUARTWORD((ulong)*w)%p;
756 4737680 : x[i+3] = HHQUARTWORD((ulong)*w)%p;
757 : }
758 106930 : return Flx_renormalize(x, lx);
759 : }
760 :
761 : static GEN
762 96731 : Flx_mulspec_quartmulii(GEN a, GEN b, ulong p, long na, long nb)
763 : {
764 96731 : GEN A = Flx_to_int_quartspec(a,na);
765 96732 : GEN B = Flx_to_int_quartspec(b,nb);
766 96732 : GEN z = mulii(A,B);
767 96733 : return int_to_Flx_quart(z,p);
768 : }
769 :
770 : /*Eval x in 2^(k*BIL) in linear time, k==2 or 3*/
771 : static GEN
772 578822 : Flx_eval2BILspec(GEN x, long k, long l)
773 : {
774 578822 : long i, lz = k*l, ki;
775 578822 : GEN pz = cgetipos(2+lz);
776 16299316 : for (i=0; i < lz; i++)
777 15720494 : *int_W(pz,i) = 0UL;
778 8439069 : for (i=0, ki=0; i<l; i++, ki+=k)
779 7860247 : *int_W(pz,ki) = x[i];
780 578822 : return int_normalize(pz,0);
781 : }
782 :
783 : static GEN
784 296393 : Z_mod2BIL_Flx_2(GEN x, long d, ulong p)
785 : {
786 296393 : long i, offset, lm = lgefint(x)-2, l = d+3;
787 296393 : ulong pi = get_Fl_red(p);
788 296393 : GEN pol = cgetg(l, t_VECSMALL);
789 296393 : pol[1] = 0;
790 7976864 : for (i=0, offset=0; offset+1 < lm; i++, offset += 2)
791 7680471 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
792 296393 : if (offset < lm)
793 224805 : pol[i+2] = (*int_W(x,offset)) % p;
794 296393 : return Flx_renormalize(pol,l);
795 : }
796 :
797 : static GEN
798 0 : Z_mod2BIL_Flx_3(GEN x, long d, ulong p)
799 : {
800 0 : long i, offset, lm = lgefint(x)-2, l = d+3;
801 0 : ulong pi = get_Fl_red(p);
802 0 : GEN pol = cgetg(l, t_VECSMALL);
803 0 : pol[1] = 0;
804 0 : for (i=0, offset=0; offset+2 < lm; i++, offset += 3)
805 0 : pol[i+2] = remlll_pre(*int_W(x,offset+2), *int_W(x,offset+1),
806 0 : *int_W(x,offset), p, pi);
807 0 : if (offset+1 < lm)
808 0 : pol[i+2] = remll_pre(*int_W(x,offset+1), *int_W(x,offset), p, pi);
809 0 : else if (offset < lm)
810 0 : pol[i+2] = (*int_W(x,offset)) % p;
811 0 : return Flx_renormalize(pol,l);
812 : }
813 :
814 : static GEN
815 293463 : Z_mod2BIL_Flx(GEN x, long bs, long d, ulong p)
816 : {
817 293463 : return bs==2 ? Z_mod2BIL_Flx_2(x, d, p): Z_mod2BIL_Flx_3(x, d, p);
818 : }
819 :
820 : static GEN
821 281970 : Flx_mulspec_mulii_inflate(GEN x, GEN y, long N, ulong p, long nx, long ny)
822 : {
823 281970 : pari_sp av = avma;
824 281970 : GEN z = mulii(Flx_eval2BILspec(x,N,nx), Flx_eval2BILspec(y,N,ny));
825 281970 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, nx+ny-2, p));
826 : }
827 :
828 : static GEN
829 19183411 : kron_pack_Flx_spec_bits(GEN x, long b, long l) {
830 : GEN y;
831 : long i;
832 19183411 : if (l == 0)
833 3418657 : return gen_0;
834 15764754 : y = cgetg(l + 1, t_VECSMALL);
835 789199241 : for(i = 1; i <= l; i++)
836 773448700 : y[i] = x[l - i];
837 15750541 : return nv_fromdigits_2k(y, b);
838 : }
839 :
840 : /* assume b < BITS_IN_LONG */
841 : static GEN
842 5679989 : kron_unpack_Flx_bits_narrow(GEN z, long b, ulong p) {
843 5679989 : GEN v = binary_2k_nv(z, b), x;
844 5679964 : long i, l = lg(v) + 1;
845 5679964 : x = cgetg(l, t_VECSMALL);
846 614350215 : for (i = 2; i < l; i++)
847 608669935 : x[i] = v[l - i] % p;
848 5680280 : return Flx_renormalize(x, l);
849 : }
850 :
851 : static GEN
852 4739273 : kron_unpack_Flx_bits_wide(GEN z, long b, ulong p, ulong pi) {
853 4739273 : GEN v = binary_2k(z, b), x, y;
854 4741339 : long i, l = lg(v) + 1, ly;
855 4741339 : x = cgetg(l, t_VECSMALL);
856 219181675 : for (i = 2; i < l; i++) {
857 214439916 : y = gel(v, l - i);
858 214439916 : ly = lgefint(y);
859 214439916 : switch (ly) {
860 8851513 : case 2: x[i] = 0; break;
861 29093460 : case 3: x[i] = *int_W_lg(y, 0, ly) % p; break;
862 161193556 : case 4: x[i] = remll_pre(*int_W_lg(y, 1, ly), *int_W_lg(y, 0, ly), p, pi); break;
863 30602904 : case 5: x[i] = remlll_pre(*int_W_lg(y, 2, ly), *int_W_lg(y, 1, ly),
864 15301387 : *int_W_lg(y, 0, ly), p, pi); break;
865 0 : default: x[i] = umodiu(gel(v, l - i), p);
866 : }
867 : }
868 4741759 : return Flx_renormalize(x, l);
869 : }
870 :
871 : static GEN
872 6439560 : Flx_mulspec_Kronecker(GEN A, GEN B, long b, ulong p, long lA, long lB)
873 : {
874 : GEN C, D;
875 6439560 : pari_sp av = avma;
876 6439560 : A = kron_pack_Flx_spec_bits(A, b, lA);
877 6443957 : B = kron_pack_Flx_spec_bits(B, b, lB);
878 6444150 : C = gerepileuptoint(av, mulii(A, B));
879 6441392 : if (b < BITS_IN_LONG)
880 2074007 : D = kron_unpack_Flx_bits_narrow(C, b, p);
881 : else
882 : {
883 4367385 : ulong pi = get_Fl_red(p);
884 4365456 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
885 : }
886 6440096 : return D;
887 : }
888 :
889 : static GEN
890 702538 : Flx_sqrspec_Kronecker(GEN A, long b, ulong p, long lA)
891 : {
892 : GEN C, D;
893 702538 : A = kron_pack_Flx_spec_bits(A, b, lA);
894 702608 : C = sqri(A);
895 702625 : if (b < BITS_IN_LONG)
896 496084 : D = kron_unpack_Flx_bits_narrow(C, b, p);
897 : else
898 : {
899 206541 : ulong pi = get_Fl_red(p);
900 206533 : D = kron_unpack_Flx_bits_wide(C, b, p, pi);
901 : }
902 702594 : return D;
903 : }
904 :
905 : /* fast product (Karatsuba) of polynomials a,b. These are not real GENs, a+2,
906 : * b+2 were sent instead. na, nb = number of terms of a, b.
907 : * Only c, c0, c1, c2 are genuine GEN.
908 : */
909 : static GEN
910 359175045 : Flx_mulspec(GEN a, GEN b, ulong p, ulong pi, long na, long nb)
911 : {
912 : GEN a0,c,c0;
913 359175045 : long n0, n0a, i, v = 0;
914 : pari_sp av;
915 :
916 464398636 : while (na && !a[0]) { a++; na--; v++; }
917 545003221 : while (nb && !b[0]) { b++; nb--; v++; }
918 359175045 : if (na < nb) swapspec(a,b, na,nb);
919 359175045 : if (!nb) return pol0_Flx(0);
920 :
921 330299245 : av = avma;
922 330299245 : if (nb >= get_Fl_threshold(p, Flx_MUL_MULII_LIMIT, Flx_MUL2_MULII_LIMIT))
923 : {
924 6836597 : long m = maxbitcoeffpol(p,nb);
925 6833623 : switch (m)
926 : {
927 96731 : case BITS_IN_QUARTULONG:
928 96731 : return Flx_shiftip(av,Flx_mulspec_quartmulii(a,b,p,na,nb), v);
929 5448 : case BITS_IN_HALFULONG:
930 5448 : return Flx_shiftip(av,Flx_mulspec_halfmulii(a,b,p,na,nb), v);
931 10299 : case BITS_IN_LONG:
932 10299 : return Flx_shiftip(av,Flx_mulspec_mulii(a,b,p,na,nb), v);
933 281970 : case 2*BITS_IN_LONG:
934 281970 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,2,p,na,nb), v);
935 0 : case 3*BITS_IN_LONG:
936 0 : return Flx_shiftip(av,Flx_mulspec_mulii_inflate(a,b,3,p,na,nb), v);
937 6439175 : default:
938 6439175 : return Flx_shiftip(av,Flx_mulspec_Kronecker(a,b,m,p,na,nb), v);
939 : }
940 : }
941 323719343 : if (nb < get_Fl_threshold(p, Flx_MUL_KARATSUBA_LIMIT, Flx_MUL2_KARATSUBA_LIMIT))
942 322221386 : return Flx_shiftip(av,Flx_mulspec_basecase(a,b,p,pi,na,nb), v);
943 1514840 : i=(na>>1); n0=na-i; na=i;
944 1514840 : a0=a+n0; n0a=n0;
945 2292713 : while (n0a && !a[n0a-1]) n0a--;
946 :
947 1514840 : if (nb > n0)
948 : {
949 : GEN b0,c1,c2;
950 : long n0b;
951 :
952 1463215 : nb -= n0; b0 = b+n0; n0b = n0;
953 2549660 : while (n0b && !b[n0b-1]) n0b--;
954 1463215 : c = Flx_mulspec(a,b,p,pi,n0a,n0b);
955 1463215 : c0 = Flx_mulspec(a0,b0,p,pi,na,nb);
956 :
957 1463215 : c2 = Flx_addspec(a0,a,p,na,n0a);
958 1463215 : c1 = Flx_addspec(b0,b,p,nb,n0b);
959 :
960 1463215 : c1 = Flx_mul_pre(c1,c2,p,pi);
961 1463215 : c2 = Flx_add(c0,c,p);
962 :
963 1463215 : c2 = Flx_neg_inplace(c2,p);
964 1463215 : c2 = Flx_add(c1,c2,p);
965 1463215 : c0 = Flx_addshift(c0,c2 ,p, n0);
966 : }
967 : else
968 : {
969 51625 : c = Flx_mulspec(a,b,p,pi,n0a,nb);
970 51625 : c0 = Flx_mulspec(a0,b,p,pi,na,nb);
971 : }
972 1514840 : c0 = Flx_addshift(c0,c,p,n0);
973 1514840 : return Flx_shiftip(av,c0, v);
974 : }
975 :
976 : GEN
977 354008031 : Flx_mul_pre(GEN x, GEN y, ulong p, ulong pi)
978 : {
979 354008031 : GEN z = Flx_mulspec(x+2,y+2,p, pi, lgpol(x),lgpol(y));
980 354202051 : z[1] = x[1]; return z;
981 : }
982 : GEN
983 20045605 : Flx_mul(GEN x, GEN y, ulong p)
984 20045605 : { return Flx_mul_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
985 :
986 : static GEN
987 272242196 : Flx_sqrspec_basecase(GEN x, ulong p, ulong pi, long nx)
988 : {
989 : long i, lz, nz;
990 : ulong p1;
991 : GEN z;
992 :
993 272242196 : if (!nx) return pol0_Flx(0);
994 272242196 : lz = (nx << 1) + 1, nz = lz-2;
995 272242196 : z = cgetg(lz, t_VECSMALL) + 2;
996 271346649 : if (!pi)
997 : {
998 210020607 : z[0] = x[0]*x[0]%p;
999 915895748 : for (i=1; i<nx; i++)
1000 : {
1001 705546196 : p1 = Flx_mullimb_ok(x+i,x,p,0, (i+1)>>1);
1002 705875141 : p1 <<= 1;
1003 705875141 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1004 705875141 : z[i] = p1 % p;
1005 : }
1006 921323222 : for ( ; i<nz; i++)
1007 : {
1008 710434003 : p1 = Flx_mullimb_ok(x+i,x,p,i-nx+1, (i+1)>>1);
1009 710973670 : p1 <<= 1;
1010 710973670 : if ((i&1) == 0) p1 += x[i>>1] * x[i>>1];
1011 710973670 : z[i] = p1 % p;
1012 : }
1013 : }
1014 : else
1015 : {
1016 61326042 : z[0] = Fl_sqr_pre(x[0], p, pi);
1017 382813975 : for (i=1; i<nx; i++)
1018 : {
1019 321472087 : p1 = Flx_mullimb(x+i,x,p,pi,0, (i+1)>>1);
1020 322099918 : p1 = Fl_add(p1, p1, p);
1021 321662540 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1022 321347656 : z[i] = p1;
1023 : }
1024 383062026 : for ( ; i<nz; i++)
1025 : {
1026 321580473 : p1 = Flx_mullimb(x+i,x,p,pi,i-nx+1, (i+1)>>1);
1027 322499584 : p1 = Fl_add(p1, p1, p);
1028 322106838 : if ((i&1) == 0) p1 = Fl_add(p1, Fl_sqr_pre(x[i>>1], p, pi), p);
1029 321720138 : z[i] = p1;
1030 : }
1031 : }
1032 272370772 : z -= 2; return Flx_renormalize(z, lz);
1033 : }
1034 :
1035 : static GEN
1036 2265 : Flx_sqrspec_sqri(GEN a, ulong p, long na)
1037 : {
1038 2265 : GEN z=sqrispec(a,na);
1039 2265 : return int_to_Flx(z,p);
1040 : }
1041 :
1042 : static GEN
1043 136 : Flx_sqrspec_halfsqri(GEN a, ulong p, long na)
1044 : {
1045 136 : GEN z = sqri(Flx_to_int_halfspec(a,na));
1046 136 : return int_to_Flx_half(z,p);
1047 : }
1048 :
1049 : static GEN
1050 10197 : Flx_sqrspec_quartsqri(GEN a, ulong p, long na)
1051 : {
1052 10197 : GEN z = sqri(Flx_to_int_quartspec(a,na));
1053 10197 : return int_to_Flx_quart(z,p);
1054 : }
1055 :
1056 : static GEN
1057 11493 : Flx_sqrspec_sqri_inflate(GEN x, long N, ulong p, long nx)
1058 : {
1059 11493 : pari_sp av = avma;
1060 11493 : GEN z = sqri(Flx_eval2BILspec(x,N,nx));
1061 11493 : return gerepileupto(av, Z_mod2BIL_Flx(z, N, (nx-1)*2, p));
1062 : }
1063 :
1064 : static GEN
1065 272615334 : Flx_sqrspec(GEN a, ulong p, ulong pi, long na)
1066 : {
1067 : GEN a0, c, c0;
1068 272615334 : long n0, n0a, i, v = 0, m;
1069 : pari_sp av;
1070 :
1071 390789324 : while (na && !a[0]) { a++; na--; v += 2; }
1072 272615334 : if (!na) return pol0_Flx(0);
1073 :
1074 272379310 : av = avma;
1075 272379310 : if (na >= get_Fl_threshold(p, Flx_SQR_SQRI_LIMIT, Flx_SQR2_SQRI_LIMIT))
1076 : {
1077 726628 : m = maxbitcoeffpol(p,na);
1078 726626 : switch(m)
1079 : {
1080 10197 : case BITS_IN_QUARTULONG:
1081 10197 : return Flx_shiftip(av, Flx_sqrspec_quartsqri(a,p,na), v);
1082 136 : case BITS_IN_HALFULONG:
1083 136 : return Flx_shiftip(av, Flx_sqrspec_halfsqri(a,p,na), v);
1084 2265 : case BITS_IN_LONG:
1085 2265 : return Flx_shiftip(av, Flx_sqrspec_sqri(a,p,na), v);
1086 11493 : case 2*BITS_IN_LONG:
1087 11493 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,2,p,na), v);
1088 0 : case 3*BITS_IN_LONG:
1089 0 : return Flx_shiftip(av, Flx_sqrspec_sqri_inflate(a,3,p,na), v);
1090 702535 : default:
1091 702535 : return Flx_shiftip(av, Flx_sqrspec_Kronecker(a,m,p,na), v);
1092 : }
1093 : }
1094 272001792 : if (na < get_Fl_threshold(p, Flx_SQR_KARATSUBA_LIMIT, Flx_SQR2_KARATSUBA_LIMIT))
1095 271941299 : return Flx_shiftip(av, Flx_sqrspec_basecase(a,p,pi,na), v);
1096 56397 : i=(na>>1); n0=na-i; na=i;
1097 56397 : a0=a+n0; n0a=n0;
1098 71162 : while (n0a && !a[n0a-1]) n0a--;
1099 :
1100 56397 : c = Flx_sqrspec(a,p,pi,n0a);
1101 56397 : c0= Flx_sqrspec(a0,p,pi,na);
1102 56397 : if (p == 2) n0 *= 2;
1103 : else
1104 : {
1105 56378 : GEN c1, t = Flx_addspec(a0,a,p,na,n0a);
1106 56378 : t = Flx_sqr_pre(t,p,pi);
1107 56378 : c1= Flx_add(c0,c, p);
1108 56378 : c1= Flx_sub(t, c1, p);
1109 56378 : c0 = Flx_addshift(c0,c1,p,n0);
1110 : }
1111 56397 : c0 = Flx_addshift(c0,c,p,n0);
1112 56397 : return Flx_shiftip(av,c0,v);
1113 : }
1114 :
1115 : GEN
1116 272330763 : Flx_sqr_pre(GEN x, ulong p, ulong pi)
1117 : {
1118 272330763 : GEN z = Flx_sqrspec(x+2,p, pi, lgpol(x));
1119 273460660 : z[1] = x[1]; return z;
1120 : }
1121 : GEN
1122 282064 : Flx_sqr(GEN x, ulong p)
1123 282064 : { return Flx_sqr_pre(x, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1124 :
1125 : GEN
1126 8316 : Flx_powu_pre(GEN x, ulong n, ulong p, ulong pi)
1127 : {
1128 8316 : GEN y = pol1_Flx(x[1]), z;
1129 : ulong m;
1130 8307 : if (n == 0) return y;
1131 8307 : m = n; z = x;
1132 : for (;;)
1133 : {
1134 32042 : if (m&1UL) y = Flx_mul_pre(y,z, p, pi);
1135 32019 : m >>= 1; if (!m) return y;
1136 23717 : z = Flx_sqr_pre(z, p, pi);
1137 : }
1138 : }
1139 : GEN
1140 0 : Flx_powu(GEN x, ulong n, ulong p)
1141 : {
1142 0 : if (n == 0) return pol1_Flx(x[1]);
1143 0 : return Flx_powu_pre(x, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
1144 : }
1145 :
1146 : GEN
1147 13633 : Flx_halve(GEN y, ulong p)
1148 : {
1149 : GEN z;
1150 : long i, l;
1151 13633 : z = cgetg_copy(y, &l); z[1] = y[1];
1152 58446 : for(i=2; i<l; i++) uel(z,i) = Fl_halve(uel(y,i), p);
1153 13633 : return z;
1154 : }
1155 :
1156 : static GEN
1157 6375097 : Flx_recipspec(GEN x, long l, long n)
1158 : {
1159 : long i;
1160 6375097 : GEN z=cgetg(n+2,t_VECSMALL)+2;
1161 113487622 : for(i=0; i<l; i++)
1162 107114122 : z[n-i-1] = x[i];
1163 14039552 : for( ; i<n; i++)
1164 7666052 : z[n-i-1] = 0;
1165 6373500 : return Flx_renormalize(z-2,n+2);
1166 : }
1167 :
1168 : GEN
1169 0 : Flx_recip(GEN x)
1170 : {
1171 0 : GEN z=Flx_recipspec(x+2,lgpol(x),lgpol(x));
1172 0 : z[1]=x[1];
1173 0 : return z;
1174 : }
1175 :
1176 : /* Return h^degpol(P) P(x / h) */
1177 : GEN
1178 1117 : Flx_rescale(GEN P, ulong h, ulong p)
1179 : {
1180 1117 : long i, l = lg(P);
1181 1117 : GEN Q = cgetg(l,t_VECSMALL);
1182 1117 : ulong hi = h;
1183 1117 : Q[l-1] = P[l-1];
1184 12538 : for (i=l-2; i>=2; i--)
1185 : {
1186 12538 : Q[i] = Fl_mul(P[i], hi, p);
1187 12538 : if (i == 2) break;
1188 11421 : hi = Fl_mul(hi,h, p);
1189 : }
1190 1117 : Q[1] = P[1]; return Q;
1191 : }
1192 :
1193 : /* x/polrecip(P)+O(x^n); allow pi = 0 */
1194 : static GEN
1195 134080 : Flx_invBarrett_basecase(GEN T, ulong p, ulong pi)
1196 : {
1197 134080 : long i, l=lg(T)-1, lr=l-1, k;
1198 134080 : GEN r=cgetg(lr,t_VECSMALL); r[1] = T[1];
1199 134080 : r[2] = 1;
1200 134080 : if (!pi)
1201 761670 : for (i=3;i<lr;i++)
1202 : {
1203 754724 : ulong u = uel(T, l-i+2);
1204 45286854 : for (k=3; k<i; k++)
1205 44532130 : { u += uel(T,l-i+k) * uel(r, k); if (u & HIGHBIT) u %= p; }
1206 754724 : r[i] = Fl_neg(u % p, p);
1207 : }
1208 : else
1209 2108040 : for (i=3;i<lr;i++)
1210 : {
1211 1980907 : ulong u = Fl_neg(uel(T,l-i+2), p);
1212 59495455 : for (k=3; k<i; k++)
1213 : {
1214 57514549 : ulong t = Fl_neg(uel(T,l-i+k), p);
1215 57514550 : u = Fl_addmul_pre(u, t, uel(r,k), p, pi);
1216 : }
1217 1980906 : r[i] = u;
1218 : }
1219 134079 : return Flx_renormalize(r,lr);
1220 : }
1221 :
1222 : /* Return new lgpol */
1223 : static long
1224 2131064 : Flx_lgrenormalizespec(GEN x, long lx)
1225 : {
1226 : long i;
1227 7401017 : for (i = lx-1; i>=0; i--)
1228 7400200 : if (x[i]) break;
1229 2131064 : return i+1;
1230 : }
1231 : /* allow pi = 0 */
1232 : static GEN
1233 23054 : Flx_invBarrett_Newton(GEN T, ulong p, ulong pi)
1234 : {
1235 23054 : long nold, lx, lz, lq, l = degpol(T), lQ;
1236 23054 : GEN q, y, z, x = zero_zv(l+1) + 2;
1237 23056 : ulong mask = quadratic_prec_mask(l-2); /* assume l > 2 */
1238 : pari_sp av;
1239 :
1240 23056 : y = T+2;
1241 23056 : q = Flx_recipspec(y,l+1,l+1); lQ = lgpol(q); q+=2;
1242 23054 : av = avma;
1243 : /* We work on _spec_ Flx's, all the l[xzq12] below are lgpol's */
1244 :
1245 : /* initialize */
1246 23054 : x[0] = Fl_inv(q[0], p);
1247 23055 : if (lQ>1 && q[1])
1248 5106 : {
1249 5106 : ulong u = q[1];
1250 5106 : if (x[0] != 1) u = Fl_mul(u, Fl_sqr(x[0],p), p);
1251 5106 : x[1] = p - u; lx = 2;
1252 : }
1253 : else
1254 17949 : lx = 1;
1255 23055 : nold = 1;
1256 158279 : for (; mask > 1; set_avma(av))
1257 : { /* set x -= x(x*q - 1) + O(t^(nnew + 1)), knowing x*q = 1 + O(t^(nold+1)) */
1258 135228 : long i, lnew, nnew = nold << 1;
1259 :
1260 135228 : if (mask & 1) nnew--;
1261 135228 : mask >>= 1;
1262 :
1263 135228 : lnew = nnew + 1;
1264 135228 : lq = Flx_lgrenormalizespec(q, minss(lQ, lnew));
1265 135239 : z = Flx_mulspec(x, q, p, pi, lx, lq); /* FIXME: high product */
1266 135226 : lz = lgpol(z); if (lz > lnew) lz = lnew;
1267 135225 : z += 2;
1268 : /* subtract 1 [=>first nold words are 0]: renormalize so that z(0) != 0 */
1269 290427 : for (i = nold; i < lz; i++) if (z[i]) break;
1270 135225 : nold = nnew;
1271 135225 : if (i >= lz) continue; /* z-1 = 0(t^(nnew + 1)) */
1272 :
1273 : /* z + i represents (x*q - 1) / t^i */
1274 100825 : lz = Flx_lgrenormalizespec (z+i, lz-i);
1275 100829 : z = Flx_mulspec(x, z+i, p, pi, lx, lz); /* FIXME: low product */
1276 100825 : lz = lgpol(z); z += 2;
1277 100825 : if (lz > lnew-i) lz = Flx_lgrenormalizespec(z, lnew-i);
1278 :
1279 100826 : lx = lz+ i;
1280 100826 : y = x + i; /* x -= z * t^i, in place */
1281 913408 : for (i = 0; i < lz; i++) y[i] = Fl_neg(z[i], p);
1282 : }
1283 23056 : x -= 2; setlg(x, lx + 2); x[1] = T[1];
1284 23056 : return x;
1285 : }
1286 :
1287 : /* allow pi = 0 */
1288 : static GEN
1289 158500 : Flx_invBarrett_pre(GEN T, ulong p, ulong pi)
1290 : {
1291 158500 : pari_sp ltop = avma;
1292 158500 : long l = lgpol(T);
1293 : GEN r;
1294 158500 : if (l < 3) return pol0_Flx(T[1]);
1295 157134 : if (l < get_Fl_threshold(p, Flx_INVBARRETT_LIMIT, Flx_INVBARRETT2_LIMIT))
1296 : {
1297 134080 : ulong c = T[l+1];
1298 134080 : if (c != 1)
1299 : {
1300 98101 : ulong ci = Fl_inv(c,p);
1301 98101 : T = Flx_Fl_mul(T, ci, p);
1302 98101 : r = Flx_invBarrett_basecase(T, p, pi);
1303 98100 : r = Flx_Fl_mul(r, ci, p);
1304 : }
1305 : else
1306 35979 : r = Flx_invBarrett_basecase(T, p, pi);
1307 : }
1308 : else
1309 23054 : r = Flx_invBarrett_Newton(T, p, pi);
1310 157136 : return gerepileuptoleaf(ltop, r);
1311 : }
1312 : GEN
1313 0 : Flx_invBarrett(GEN T, ulong p)
1314 0 : { return Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1315 :
1316 : /* allow pi = 0 */
1317 : GEN
1318 96618957 : Flx_get_red_pre(GEN T, ulong p, ulong pi)
1319 : {
1320 96618957 : if (typ(T)!=t_VECSMALL
1321 96588326 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1322 : Flx_BARRETT2_LIMIT))
1323 96592730 : return T;
1324 7629 : retmkvec2(Flx_invBarrett_pre(T, p, pi),T);
1325 : }
1326 : GEN
1327 13563433 : Flx_get_red(GEN T, ulong p)
1328 : {
1329 13563433 : if (typ(T)!=t_VECSMALL
1330 13563460 : || lgpol(T) < get_Fl_threshold(p, Flx_BARRETT_LIMIT,
1331 : Flx_BARRETT2_LIMIT))
1332 13556518 : return T;
1333 5220 : retmkvec2(Flx_invBarrett_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)),T);
1334 : }
1335 :
1336 : /* separate from Flx_divrem for maximal speed. */
1337 : static GEN
1338 704210591 : Flx_rem_basecase(GEN x, GEN y, ulong p, ulong pi)
1339 : {
1340 : pari_sp av;
1341 : GEN z, c;
1342 : long dx,dy,dy1,dz,i,j;
1343 : ulong p1,inv;
1344 704210591 : long vs=x[1];
1345 :
1346 704210591 : dy = degpol(y); if (!dy) return pol0_Flx(x[1]);
1347 673526082 : dx = degpol(x);
1348 673739140 : dz = dx-dy; if (dz < 0) return Flx_copy(x);
1349 673739140 : x += 2; y += 2;
1350 673739140 : inv = y[dy];
1351 673739140 : if (inv != 1UL) inv = Fl_inv(inv,p);
1352 813885202 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1353 :
1354 675116064 : c = cgetg(dy+3, t_VECSMALL); c[1]=vs; c += 2; av=avma;
1355 672953718 : z = cgetg(dz+3, t_VECSMALL); z[1]=vs; z += 2;
1356 :
1357 671465617 : if (!pi)
1358 : {
1359 471744183 : z[dz] = (inv*x[dx]) % p;
1360 1816678069 : for (i=dx-1; i>=dy; --i)
1361 : {
1362 1344933886 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1363 10730206315 : for (j=i-dy1; j<=i && j<=dz; j++)
1364 : {
1365 9385272429 : p1 += z[j]*y[i-j];
1366 9385272429 : if (p1 & HIGHBIT) p1 %= p;
1367 : }
1368 1344933886 : p1 %= p;
1369 1344933886 : z[i-dy] = p1? ((p - p1)*inv) % p: 0;
1370 : }
1371 3299479142 : for (i=0; i<dy; i++)
1372 : {
1373 2827137742 : p1 = z[0]*y[i];
1374 14805807026 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1375 : {
1376 11978669284 : p1 += z[j]*y[i-j];
1377 11978669284 : if (p1 & HIGHBIT) p1 %= p;
1378 : }
1379 2827513792 : c[i] = Fl_sub(x[i], p1%p, p);
1380 : }
1381 : }
1382 : else
1383 : {
1384 199721434 : z[dz] = Fl_mul_pre(inv, x[dx], p, pi);
1385 663534219 : for (i=dx-1; i>=dy; --i)
1386 : {
1387 463803774 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1388 1994558547 : for (j=i-dy1; j<=i && j<=dz; j++)
1389 1530883956 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1390 463674591 : z[i-dy] = p1? Fl_mul_pre(p - p1, inv, p, pi): 0;
1391 : }
1392 1465287400 : for (i=0; i<dy; i++)
1393 : {
1394 1266028474 : p1 = Fl_mul_pre(z[0],y[i],p,pi);
1395 3639711105 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1396 2366682489 : p1 = Fl_addmul_pre(p1, z[j], y[i - j], p, pi);
1397 1255172978 : c[i] = Fl_sub(x[i], p1, p);
1398 : }
1399 : }
1400 835896431 : i = dy-1; while (i>=0 && !c[i]) i--;
1401 671600326 : set_avma(av); return Flx_renormalize(c-2, i+3);
1402 : }
1403 :
1404 : /* as FpX_divrem but working only on ulong types.
1405 : * if relevant, *pr is the last object on stack */
1406 : static GEN
1407 60530125 : Flx_divrem_basecase(GEN x, GEN y, ulong p, ulong pi, GEN *pr)
1408 : {
1409 : GEN z,q,c;
1410 : long dx,dy,dy1,dz,i,j;
1411 : ulong p1,inv;
1412 60530125 : long sv=x[1];
1413 :
1414 60530125 : dy = degpol(y);
1415 60528458 : if (dy<0) pari_err_INV("Flx_divrem",y);
1416 60528894 : if (pr == ONLY_REM) return Flx_rem_basecase(x, y, p, pi);
1417 60528448 : if (!dy)
1418 : {
1419 7166693 : if (pr && pr != ONLY_DIVIDES) *pr = pol0_Flx(sv);
1420 7166561 : if (y[2] == 1UL) return Flx_copy(x);
1421 4979630 : return Flx_Fl_mul(x, Fl_inv(y[2], p), p);
1422 : }
1423 53361755 : dx = degpol(x);
1424 53367576 : dz = dx-dy;
1425 53367576 : if (dz < 0)
1426 : {
1427 928838 : q = pol0_Flx(sv);
1428 928824 : if (pr && pr != ONLY_DIVIDES) *pr = Flx_copy(x);
1429 928824 : return q;
1430 : }
1431 52438738 : x += 2;
1432 52438738 : y += 2;
1433 52438738 : z = cgetg(dz + 3, t_VECSMALL); z[1] = sv; z += 2;
1434 52424732 : inv = uel(y, dy);
1435 52424732 : if (inv != 1UL) inv = Fl_inv(inv,p);
1436 77392965 : for (dy1=dy-1; dy1>=0 && !y[dy1]; dy1--);
1437 :
1438 52430882 : if (SMALL_ULONG(p))
1439 : {
1440 50570611 : z[dz] = (inv*x[dx]) % p;
1441 129893024 : for (i=dx-1; i>=dy; --i)
1442 : {
1443 79322413 : p1 = p - x[i]; /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1444 259964268 : for (j=i-dy1; j<=i && j<=dz; j++)
1445 : {
1446 180641855 : p1 += z[j]*y[i-j];
1447 180641855 : if (p1 & HIGHBIT) p1 %= p;
1448 : }
1449 79322413 : p1 %= p;
1450 79322413 : z[i-dy] = p1? (long) ((p - p1)*inv) % p: 0;
1451 : }
1452 : }
1453 : else
1454 : {
1455 1860271 : z[dz] = Fl_mul(inv, x[dx], p);
1456 9191718 : for (i=dx-1; i>=dy; --i)
1457 : { /* compute -p1 instead of p1 (pb with ulongs otherwise) */
1458 7333721 : p1 = p - uel(x,i);
1459 26278588 : for (j=i-dy1; j<=i && j<=dz; j++)
1460 18944870 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1461 7333718 : z[i-dy] = p1? Fl_mul(p - p1, inv, p): 0;
1462 : }
1463 : }
1464 52428608 : q = Flx_renormalize(z-2, dz+3);
1465 52435448 : if (!pr) return q;
1466 :
1467 26352230 : c = cgetg(dy + 3, t_VECSMALL); c[1] = sv; c += 2;
1468 26355857 : if (SMALL_ULONG(p))
1469 : {
1470 233205488 : for (i=0; i<dy; i++)
1471 : {
1472 208473359 : p1 = (ulong)z[0]*y[i];
1473 488860740 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1474 : {
1475 280387381 : p1 += (ulong)z[j]*y[i-j];
1476 280387381 : if (p1 & HIGHBIT) p1 %= p;
1477 : }
1478 208473200 : c[i] = Fl_sub(x[i], p1%p, p);
1479 : }
1480 : }
1481 : else
1482 : {
1483 17168531 : for (i=0; i<dy; i++)
1484 : {
1485 15544301 : p1 = Fl_mul(z[0],y[i],p);
1486 52509007 : for (j=maxss(1,i-dy1); j<=i && j<=dz; j++)
1487 36964704 : p1 = Fl_add(p1, Fl_mul(z[j],y[i-j],p), p);
1488 15544303 : c[i] = Fl_sub(x[i], p1, p);
1489 : }
1490 : }
1491 35452429 : i=dy-1; while (i>=0 && !c[i]) i--;
1492 26356359 : c = Flx_renormalize(c-2, i+3);
1493 26356112 : if (pr == ONLY_DIVIDES)
1494 420 : { if (lg(c) != 2) return NULL; }
1495 : else
1496 26355692 : *pr = c;
1497 26355972 : return q;
1498 : }
1499 :
1500 : /* Compute x mod T where 2 <= degpol(T) <= l+1 <= 2*(degpol(T)-1)
1501 : * and mg is the Barrett inverse of T. */
1502 : static GEN
1503 905077 : Flx_divrem_Barrettspec(GEN x, long l, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1504 : {
1505 : GEN q, r;
1506 905077 : long lt = degpol(T); /*We discard the leading term*/
1507 : long ld, lm, lT, lmg;
1508 905016 : ld = l-lt;
1509 905016 : lm = minss(ld, lgpol(mg));
1510 905173 : lT = Flx_lgrenormalizespec(T+2,lt);
1511 905275 : lmg = Flx_lgrenormalizespec(mg+2,lm);
1512 905206 : q = Flx_recipspec(x+lt,ld,ld); /* q = rec(x) lz<=ld*/
1513 904640 : q = Flx_mulspec(q+2,mg+2,p,pi,lgpol(q),lmg); /* q = rec(x) * mg lz<=ld+lm*/
1514 905165 : q = Flx_recipspec(q+2,minss(ld,lgpol(q)),ld);/* q = rec (rec(x) * mg) lz<=ld*/
1515 904553 : if (!pr) return q;
1516 896899 : r = Flx_mulspec(q+2,T+2,p,pi,lgpol(q),lT); /* r = q*pol lz<=ld+lt*/
1517 897558 : r = Flx_subspec(x,r+2,p,lt,minss(lt,lgpol(r)));/* r = x - q*pol lz<=lt */
1518 897346 : if (pr == ONLY_REM) return r;
1519 424116 : *pr = r; return q;
1520 : }
1521 :
1522 : static GEN
1523 608341 : Flx_divrem_Barrett(GEN x, GEN mg, GEN T, ulong p, ulong pi, GEN *pr)
1524 : {
1525 608341 : GEN q = NULL, r = Flx_copy(x);
1526 608373 : long l = lgpol(x), lt = degpol(T), lm = 2*lt-1, v = T[1];
1527 : long i;
1528 608367 : if (l <= lt)
1529 : {
1530 0 : if (pr == ONLY_REM) return Flx_copy(x);
1531 0 : if (pr == ONLY_DIVIDES) return lgpol(x)? NULL: pol0_Flx(v);
1532 0 : if (pr) *pr = Flx_copy(x);
1533 0 : return pol0_Flx(v);
1534 : }
1535 608367 : if (lt <= 1)
1536 1366 : return Flx_divrem_basecase(x,T,p,pi,pr);
1537 607001 : if (pr != ONLY_REM && l>lm)
1538 28877 : { q = zero_zv(l-lt+1); q[1] = T[1]; }
1539 906578 : while (l>lm)
1540 : {
1541 299727 : GEN zr, zq = Flx_divrem_Barrettspec(r+2+l-lm,lm,mg,T,p,pi,&zr);
1542 299623 : long lz = lgpol(zr);
1543 299577 : if (pr != ONLY_REM)
1544 : {
1545 57046 : long lq = lgpol(zq);
1546 851892 : for(i=0; i<lq; i++) q[2+l-lm+i] = zq[2+i];
1547 : }
1548 4348608 : for(i=0; i<lz; i++) r[2+l-lm+i] = zr[2+i];
1549 299577 : l = l-lm+lz;
1550 : }
1551 606851 : if (pr == ONLY_REM)
1552 : {
1553 473300 : if (l > lt)
1554 473258 : r = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi,ONLY_REM);
1555 : else
1556 42 : r = Flx_renormalize(r, l+2);
1557 473272 : r[1] = v; return r;
1558 : }
1559 133551 : if (l > lt)
1560 : {
1561 132123 : GEN zq = Flx_divrem_Barrettspec(r+2,l,mg,T,p,pi, pr ? &r: NULL);
1562 132123 : if (!q) q = zq;
1563 : else
1564 : {
1565 27300 : long lq = lgpol(zq);
1566 158521 : for(i=0; i<lq; i++) q[2+i] = zq[2+i];
1567 : }
1568 : }
1569 1428 : else if (pr)
1570 1535 : r = Flx_renormalize(r, l+2);
1571 133551 : q[1] = v; q = Flx_renormalize(q, lg(q));
1572 133700 : if (pr == ONLY_DIVIDES) return lgpol(r)? NULL: q;
1573 133700 : if (pr) { r[1] = v; *pr = r; }
1574 133700 : return q;
1575 : }
1576 :
1577 : /* allow pi = 0 (SMALL_ULONG) */
1578 : GEN
1579 77678001 : Flx_divrem_pre(GEN x, GEN T, ulong p, ulong pi, GEN *pr)
1580 : {
1581 : GEN B, y;
1582 : long dy, dx, d;
1583 77678001 : if (pr==ONLY_REM) return Flx_rem_pre(x, T, p, pi);
1584 60654648 : y = get_Flx_red(T, &B);
1585 60674652 : dy = degpol(y); dx = degpol(x); d = dx-dy;
1586 60663241 : if (!B && d+3 < get_Fl_threshold(p, Flx_DIVREM_BARRETT_LIMIT,Flx_DIVREM2_BARRETT_LIMIT))
1587 60526993 : return Flx_divrem_basecase(x,y,p,pi,pr);
1588 : else
1589 : {
1590 134620 : pari_sp av = avma;
1591 134620 : GEN mg = B? B: Flx_invBarrett_pre(y, p, pi);
1592 134620 : GEN q1 = Flx_divrem_Barrett(x,mg,y,p,pi,pr);
1593 134620 : if (!q1) return gc_NULL(av);
1594 134620 : if (!pr || pr==ONLY_DIVIDES) return gerepileuptoleaf(av, q1);
1595 126369 : return gc_all(av, 2, &q1, pr);
1596 : }
1597 : }
1598 : GEN
1599 29669924 : Flx_divrem(GEN x, GEN T, ulong p, GEN *pr)
1600 29669924 : { return Flx_divrem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p), pr); }
1601 :
1602 : GEN
1603 825158653 : Flx_rem_pre(GEN x, GEN T, ulong p, ulong pi)
1604 : {
1605 825158653 : GEN B, y = get_Flx_red(T, &B);
1606 824879197 : long d = degpol(x) - degpol(y);
1607 824671730 : if (d < 0) return Flx_copy(x);
1608 704864591 : if (!B && d+3 < get_Fl_threshold(p, Flx_REM_BARRETT_LIMIT,Flx_REM2_BARRETT_LIMIT))
1609 704104073 : return Flx_rem_basecase(x,y,p, pi);
1610 : else
1611 : {
1612 473720 : pari_sp av=avma;
1613 473720 : GEN mg = B ? B: Flx_invBarrett_pre(y, p, pi);
1614 473722 : GEN r = Flx_divrem_Barrett(x, mg, y, p, pi, ONLY_REM);
1615 473716 : return gerepileuptoleaf(av, r);
1616 : }
1617 : }
1618 : GEN
1619 40812235 : Flx_rem(GEN x, GEN T, ulong p)
1620 40812235 : { return Flx_rem_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1621 :
1622 : /* reduce T mod (X^n - 1, p). Shallow function */
1623 : GEN
1624 5023973 : Flx_mod_Xnm1(GEN T, ulong n, ulong p)
1625 : {
1626 5023973 : long i, j, L = lg(T), l = n+2;
1627 : GEN S;
1628 5023973 : if (L <= l || n & ~LGBITS) return T;
1629 3470 : S = cgetg(l, t_VECSMALL);
1630 3470 : S[1] = T[1];
1631 13838 : for (i = 2; i < l; i++) S[i] = T[i];
1632 9594 : for (j = 2; i < L; i++) {
1633 6124 : S[j] = Fl_add(S[j], T[i], p);
1634 6124 : if (++j == l) j = 2;
1635 : }
1636 3470 : return Flx_renormalize(S, l);
1637 : }
1638 : /* reduce T mod (X^n + 1, p). Shallow function */
1639 : GEN
1640 27539 : Flx_mod_Xn1(GEN T, ulong n, ulong p)
1641 : {
1642 27539 : long i, j, L = lg(T), l = n+2;
1643 : GEN S;
1644 27539 : if (L <= l || n & ~LGBITS) return T;
1645 2721 : S = cgetg(l, t_VECSMALL);
1646 2721 : S[1] = T[1];
1647 11227 : for (i = 2; i < l; i++) S[i] = T[i];
1648 7214 : for (j = 2; i < L; i++) {
1649 4493 : S[j] = Fl_sub(S[j], T[i], p);
1650 4493 : if (++j == l) j = 2;
1651 : }
1652 2721 : return Flx_renormalize(S, l);
1653 : }
1654 :
1655 : struct _Flxq {
1656 : GEN aut, T;
1657 : ulong p, pi;
1658 : };
1659 : /* allow pi = 0 */
1660 : static void
1661 70366224 : set_Flxq_pre(struct _Flxq *D, GEN T, ulong p, ulong pi)
1662 : {
1663 70366224 : D->p = p;
1664 70366224 : D->pi = pi;
1665 70366224 : D->T = Flx_get_red_pre(T, p, pi);
1666 70375884 : }
1667 : static void
1668 67942 : set_Flxq(struct _Flxq *D, GEN T, ulong p)
1669 67942 : { set_Flxq_pre(D, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
1670 :
1671 : static GEN
1672 0 : _Flx_divrem(void * E, GEN x, GEN y, GEN *r)
1673 : {
1674 0 : struct _Flxq *D = (struct _Flxq*) E;
1675 0 : return Flx_divrem_pre(x, y, D->p, D->pi, r);
1676 : }
1677 : static GEN
1678 582784 : _Flx_add(void * E, GEN x, GEN y) {
1679 582784 : struct _Flxq *D = (struct _Flxq*) E;
1680 582784 : return Flx_add(x, y, D->p);
1681 : }
1682 : static GEN
1683 10175980 : _Flx_mul(void *E, GEN x, GEN y) {
1684 10175980 : struct _Flxq *D = (struct _Flxq*) E;
1685 10175980 : return Flx_mul_pre(x, y, D->p, D->pi);
1686 : }
1687 : static GEN
1688 0 : _Flx_sqr(void *E, GEN x) {
1689 0 : struct _Flxq *D = (struct _Flxq*) E;
1690 0 : return Flx_sqr_pre(x, D->p, D->pi);
1691 : }
1692 :
1693 : static struct bb_ring Flx_ring = { _Flx_add,_Flx_mul,_Flx_sqr };
1694 :
1695 : GEN
1696 0 : Flx_digits(GEN x, GEN T, ulong p)
1697 : {
1698 : struct _Flxq D;
1699 0 : long d = degpol(T), n = (lgpol(x)+d-1)/d;
1700 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1701 0 : return gen_digits(x,T,n,(void *)&D, &Flx_ring, _Flx_divrem);
1702 : }
1703 :
1704 : GEN
1705 0 : FlxV_Flx_fromdigits(GEN x, GEN T, ulong p)
1706 : {
1707 : struct _Flxq D;
1708 0 : D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
1709 0 : return gen_fromdigits(x,T,(void *)&D, &Flx_ring);
1710 : }
1711 :
1712 : long
1713 3657626 : Flx_val(GEN x)
1714 : {
1715 3657626 : long i, l=lg(x);
1716 3657626 : if (l==2) return LONG_MAX;
1717 3666025 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1718 3657626 : return i-2;
1719 : }
1720 : long
1721 25504004 : Flx_valrem(GEN x, GEN *Z)
1722 : {
1723 25504004 : long v, i, l=lg(x);
1724 : GEN y;
1725 25504004 : if (l==2) { *Z = Flx_copy(x); return LONG_MAX; }
1726 27673064 : for (i=2; i<l && x[i]==0; i++) /*empty*/;
1727 25504004 : v = i-2;
1728 25504004 : if (v == 0) { *Z = x; return 0; }
1729 1022115 : l -= v;
1730 1022115 : y = cgetg(l, t_VECSMALL); y[1] = x[1];
1731 2612663 : for (i=2; i<l; i++) y[i] = x[i+v];
1732 1017952 : *Z = y; return v;
1733 : }
1734 :
1735 : GEN
1736 18268873 : Flx_deriv(GEN z, ulong p)
1737 : {
1738 18268873 : long i,l = lg(z)-1;
1739 : GEN x;
1740 18268873 : if (l < 2) l = 2;
1741 18268873 : x = cgetg(l, t_VECSMALL); x[1] = z[1]; z++;
1742 18260294 : if (HIGHWORD(l | p))
1743 52616903 : for (i=2; i<l; i++) x[i] = Fl_mul((ulong)i-1, z[i], p);
1744 : else
1745 74656874 : for (i=2; i<l; i++) x[i] = ((i-1) * z[i]) % p;
1746 18261216 : return Flx_renormalize(x,l);
1747 : }
1748 :
1749 : static GEN
1750 422556 : Flx_integXn(GEN x, long n, ulong p)
1751 : {
1752 422556 : long i, lx = lg(x);
1753 : GEN y;
1754 422556 : if (lx == 2) return Flx_copy(x);
1755 412742 : y = cgetg(lx, t_VECSMALL); y[1] = x[1];
1756 2087976 : for (i=2; i<lx; i++)
1757 : {
1758 1674662 : ulong xi = uel(x,i);
1759 1674662 : if (xi == 0)
1760 12505 : uel(y,i) = 0;
1761 : else
1762 : {
1763 1662157 : ulong j = n+i-1;
1764 1662157 : ulong d = ugcd(j, xi);
1765 1662084 : if (d==1)
1766 1014521 : uel(y,i) = Fl_div(xi, j, p);
1767 : else
1768 647563 : uel(y,i) = Fl_div(xi/d, j/d, p);
1769 : }
1770 : }
1771 413314 : return Flx_renormalize(y, lx);;
1772 : }
1773 :
1774 : GEN
1775 0 : Flx_integ(GEN x, ulong p)
1776 : {
1777 0 : long i, lx = lg(x);
1778 : GEN y;
1779 0 : if (lx == 2) return Flx_copy(x);
1780 0 : y = cgetg(lx+1, t_VECSMALL); y[1] = x[1];
1781 0 : uel(y,2) = 0;
1782 0 : for (i=3; i<=lx; i++)
1783 0 : uel(y,i) = uel(x,i-1) ? Fl_div(uel(x,i-1), (i-2)%p, p): 0UL;
1784 0 : return Flx_renormalize(y, lx+1);;
1785 : }
1786 :
1787 : /* assume p prime */
1788 : GEN
1789 13447 : Flx_diff1(GEN P, ulong p)
1790 : {
1791 13447 : return Flx_sub(Flx_translate1(P, p), P, p);
1792 : }
1793 :
1794 : GEN
1795 344365 : Flx_deflate(GEN x0, long d)
1796 : {
1797 : GEN z, y, x;
1798 344365 : long i,id, dy, dx = degpol(x0);
1799 344365 : if (d == 1 || dx <= 0) return Flx_copy(x0);
1800 281155 : dy = dx/d;
1801 281155 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1802 281147 : z = y + 2;
1803 281147 : x = x0+ 2;
1804 929564 : for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1805 281147 : return y;
1806 : }
1807 :
1808 : GEN
1809 163784 : Flx_inflate(GEN x0, long d)
1810 : {
1811 163784 : long i, id, dy, dx = degpol(x0);
1812 163776 : GEN x = x0 + 2, z, y;
1813 163776 : if (dx <= 0) return Flx_copy(x0);
1814 162651 : dy = dx*d;
1815 162651 : y = cgetg(dy+3, t_VECSMALL); y[1] = x0[1];
1816 162624 : z = y + 2;
1817 9256178 : for (i=0; i<=dy; i++) z[i] = 0;
1818 4510545 : for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1819 162624 : return y;
1820 : }
1821 :
1822 : /* write p(X) = a_0(X^k) + X*a_1(X^k) + ... + X^(k-1)*a_{k-1}(X^k) */
1823 : GEN
1824 149030 : Flx_splitting(GEN p, long k)
1825 : {
1826 149030 : long n = degpol(p), v = p[1], m, i, j, l;
1827 : GEN r;
1828 :
1829 149024 : m = n/k;
1830 149024 : r = cgetg(k+1,t_VEC);
1831 684329 : for(i=1; i<=k; i++)
1832 : {
1833 535313 : gel(r,i) = cgetg(m+3, t_VECSMALL);
1834 535315 : mael(r,i,1) = v;
1835 : }
1836 4673565 : for (j=1, i=0, l=2; i<=n; i++)
1837 : {
1838 4524549 : mael(r,j,l) = p[2+i];
1839 4524549 : if (j==k) { j=1; l++; } else j++;
1840 : }
1841 684370 : for(i=1; i<=k; i++)
1842 535375 : gel(r,i) = Flx_renormalize(gel(r,i),i<j?l+1:l);
1843 148995 : return r;
1844 : }
1845 :
1846 : /* ux + vy */
1847 : static GEN
1848 23076 : Flx_addmulmul(GEN u, GEN v, GEN x, GEN y, ulong p, ulong pi)
1849 23076 : { return Flx_add(Flx_mul_pre(u,x, p,pi), Flx_mul_pre(v,y, p,pi), p); }
1850 :
1851 : static GEN
1852 11532 : FlxM_Flx_mul2(GEN M, GEN x, GEN y, ulong p, ulong pi)
1853 : {
1854 11532 : GEN res = cgetg(3, t_COL);
1855 11532 : gel(res, 1) = Flx_addmulmul(gcoeff(M,1,1), gcoeff(M,1,2), x, y, p, pi);
1856 11532 : gel(res, 2) = Flx_addmulmul(gcoeff(M,2,1), gcoeff(M,2,2), x, y, p, pi);
1857 11532 : return res;
1858 : }
1859 :
1860 : #if 0
1861 : static GEN
1862 : FlxM_mul2_old(GEN M, GEN N, ulong p)
1863 : {
1864 : GEN res = cgetg(3, t_MAT);
1865 : gel(res, 1) = FlxM_Flx_mul2(M,gcoeff(N,1,1),gcoeff(N,2,1),p);
1866 : gel(res, 2) = FlxM_Flx_mul2(M,gcoeff(N,1,2),gcoeff(N,2,2),p);
1867 : return res;
1868 : }
1869 : #endif
1870 : /* A,B are 2x2 matrices, Flx entries. Return A x B using Strassen 7M formula */
1871 : static GEN
1872 2109 : FlxM_mul2(GEN A, GEN B, ulong p, ulong pi)
1873 : {
1874 2109 : GEN A11=gcoeff(A,1,1),A12=gcoeff(A,1,2), B11=gcoeff(B,1,1),B12=gcoeff(B,1,2);
1875 2109 : GEN A21=gcoeff(A,2,1),A22=gcoeff(A,2,2), B21=gcoeff(B,2,1),B22=gcoeff(B,2,2);
1876 2109 : GEN M1 = Flx_mul_pre(Flx_add(A11,A22, p), Flx_add(B11,B22, p), p, pi);
1877 2109 : GEN M2 = Flx_mul_pre(Flx_add(A21,A22, p), B11, p, pi);
1878 2109 : GEN M3 = Flx_mul_pre(A11, Flx_sub(B12,B22, p), p, pi);
1879 2109 : GEN M4 = Flx_mul_pre(A22, Flx_sub(B21,B11, p), p, pi);
1880 2109 : GEN M5 = Flx_mul_pre(Flx_add(A11,A12, p), B22, p, pi);
1881 2109 : GEN M6 = Flx_mul_pre(Flx_sub(A21,A11, p), Flx_add(B11,B12, p), p, pi);
1882 2109 : GEN M7 = Flx_mul_pre(Flx_sub(A12,A22, p), Flx_add(B21,B22, p), p, pi);
1883 2109 : GEN T1 = Flx_add(M1,M4, p), T2 = Flx_sub(M7,M5, p);
1884 2109 : GEN T3 = Flx_sub(M1,M2, p), T4 = Flx_add(M3,M6, p);
1885 2109 : retmkmat2(mkcol2(Flx_add(T1,T2, p), Flx_add(M2,M4, p)),
1886 : mkcol2(Flx_add(M3,M5, p), Flx_add(T3,T4, p)));
1887 : }
1888 :
1889 : /* Return [0,1;1,-q]*M */
1890 : static GEN
1891 2103 : Flx_FlxM_qmul(GEN q, GEN M, ulong p, ulong pi)
1892 : {
1893 2103 : GEN u, v, res = cgetg(3, t_MAT);
1894 2103 : u = Flx_sub(gcoeff(M,1,1), Flx_mul_pre(gcoeff(M,2,1), q, p, pi), p);
1895 2103 : gel(res,1) = mkcol2(gcoeff(M,2,1), u);
1896 2103 : v = Flx_sub(gcoeff(M,1,2), Flx_mul_pre(gcoeff(M,2,2), q, p, pi), p);
1897 2103 : gel(res,2) = mkcol2(gcoeff(M,2,2), v);
1898 2103 : return res;
1899 : }
1900 :
1901 : static GEN
1902 851 : matid2_FlxM(long v)
1903 : {
1904 851 : return mkmat2(mkcol2(pol1_Flx(v),pol0_Flx(v)),
1905 : mkcol2(pol0_Flx(v),pol1_Flx(v)));
1906 : }
1907 :
1908 : struct Flx_res
1909 : {
1910 : ulong res, lc;
1911 : long deg0, deg1, off;
1912 : };
1913 :
1914 : INLINE void
1915 1696 : Flx_halfres_update_pre(long da, long db, long dr, ulong p, ulong pi, struct Flx_res *res)
1916 : {
1917 1696 : if (dr >= 0)
1918 : {
1919 1696 : if (res->lc != 1)
1920 : {
1921 3 : if (pi)
1922 : {
1923 0 : res->lc = Fl_powu_pre(res->lc, da - dr, p, pi);
1924 0 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1925 : } else
1926 : {
1927 3 : res->lc = Fl_powu(res->lc, da - dr, p);
1928 3 : res->res = Fl_mul(res->res, res->lc, p);
1929 : }
1930 : }
1931 1696 : if (both_odd(da + res->off, db + res->off))
1932 0 : res->res = Fl_neg(res->res, p);
1933 : } else
1934 : {
1935 0 : if (db == 0)
1936 : {
1937 0 : if (res->lc != 1)
1938 : {
1939 0 : if (pi)
1940 : {
1941 0 : res->lc = Fl_powu_pre(res->lc, da, p, pi);
1942 0 : res->res = Fl_mul_pre(res->res, res->lc, p, pi);
1943 : } else
1944 : {
1945 0 : res->lc = Fl_powu(res->lc, da, p);
1946 0 : res->res = Fl_mul(res->res, res->lc, p);
1947 : }
1948 : }
1949 : } else
1950 0 : res->res = 0;
1951 : }
1952 1696 : }
1953 :
1954 : static GEN
1955 728067 : Flx_halfres_basecase(GEN a, GEN b, ulong p, ulong pi, struct Flx_res *res)
1956 : {
1957 728067 : pari_sp av=avma;
1958 : GEN u,u1,v,v1;
1959 728067 : long vx = a[1];
1960 728067 : long n = lgpol(a)>>1;
1961 728061 : u1 = v = pol0_Flx(vx);
1962 728054 : u = v1 = pol1_Flx(vx);
1963 2792716 : while (lgpol(b)>n)
1964 : {
1965 : GEN r, q;
1966 2064634 : q = Flx_divrem_pre(a,b,p,pi, &r);
1967 2064675 : if (res)
1968 : {
1969 851 : long da = degpol(a), db=degpol(b), dr = degpol(r);
1970 851 : res->lc = Flx_lead(b);
1971 851 : if (dr >= n)
1972 3 : Flx_halfres_update_pre(da, db, dr, p, pi, res);
1973 : else
1974 : {
1975 848 : res->deg0 = da;
1976 848 : res->deg1 = db;
1977 : }
1978 : }
1979 2064675 : a = b; b = r; swap(u,u1); swap(v,v1);
1980 2064675 : u1 = Flx_sub(u1, Flx_mul(u, q, p), p);
1981 2064633 : v1 = Flx_sub(v1, Flx_mul(v, q, p), p);
1982 2064673 : if (gc_needed(av,2))
1983 : {
1984 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_halfgcd (d = %ld)",degpol(b));
1985 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
1986 : }
1987 : }
1988 727898 : return gerepilecopy(av, mkmat2(mkcol2(u,u1), mkcol2(v,v1)));
1989 : }
1990 :
1991 : static GEN Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi, struct Flx_res *res);
1992 :
1993 : static GEN
1994 11514 : Flx_halfres_split(GEN x, GEN y, ulong p, ulong pi, struct Flx_res *res)
1995 : {
1996 11514 : pari_sp av = avma;
1997 : GEN R, S, V;
1998 : GEN x1, y1, r, q;
1999 11514 : long l = lgpol(x), n = l>>1, k;
2000 11514 : if (lgpol(y) <= n) return matid2_FlxM(x[1]);
2001 10669 : if (res)
2002 : {
2003 3054 : res->lc = Flx_lead(y);
2004 3054 : res->deg0 -= n;
2005 3054 : res->deg1 -= n;
2006 3054 : res->off += n;
2007 : }
2008 10669 : R = Flx_halfres_i(Flx_shift(x,-n),Flx_shift(y,-n),p,pi,res);
2009 10669 : if (res)
2010 : {
2011 3054 : res->off -= n;
2012 3054 : res->deg0 += n;
2013 3054 : res->deg1 += n;
2014 : }
2015 10669 : V = FlxM_Flx_mul2(R,x,y,p,pi); x1 = gel(V,1); y1 = gel(V,2);
2016 10669 : if (lgpol(y1) <= n) return gerepilecopy(av, R);
2017 2103 : k = 2*n-degpol(y1);
2018 2103 : q = Flx_divrem_pre(x1, y1, p, pi, &r);
2019 2103 : if (res)
2020 : {
2021 845 : long dx1 = degpol(x1), dy1 = degpol(y1), dr = degpol(r);
2022 845 : if (dy1 < degpol(y))
2023 0 : Flx_halfres_update_pre(res->deg0, res->deg1, dy1, p, pi, res);
2024 845 : res->lc = Flx_lead(y1);
2025 845 : res->deg0 = dx1;
2026 845 : res->deg1 = dy1;
2027 845 : if (dr >= n)
2028 : {
2029 845 : Flx_halfres_update_pre(dx1, dy1, dr, p, pi, res);
2030 845 : res->deg0 = dy1;
2031 845 : res->deg1 = dr;
2032 : }
2033 845 : res->deg0 -= k;
2034 845 : res->deg1 -= k;
2035 845 : res->off += k;
2036 : }
2037 2103 : S = Flx_halfres_i(Flx_shift(y1,-k), Flx_shift(r,-k),p,pi,res);
2038 2103 : if (res)
2039 : {
2040 845 : res->deg0 += k;
2041 845 : res->deg1 += k;
2042 845 : res->off -= k;
2043 : }
2044 2103 : return gerepileupto(av, FlxM_mul2(S,Flx_FlxM_qmul(q,R,p,pi),p,pi));
2045 : }
2046 :
2047 : static GEN
2048 739586 : Flx_halfres_i(GEN x, GEN y, ulong p, ulong pi,struct Flx_res *res)
2049 : {
2050 739586 : if (lgpol(x) < get_Fl_threshold(p, Flx_HALFGCD_LIMIT, Flx_HALFGCD2_LIMIT))
2051 728067 : return Flx_halfres_basecase(x, y, p, pi, res);
2052 11514 : return Flx_halfres_split(x, y, p, pi, res);
2053 : }
2054 :
2055 : static GEN
2056 725966 : Flx_halfgcd_i(GEN x, GEN y, ulong p, ulong pi)
2057 725966 : { return Flx_halfres_i(x, y, p, pi, NULL); }
2058 :
2059 : /* Return M in GL_2(Fl[X]) such that:
2060 : if [a',b']~=M*[a,b]~ then degpol(a')>= (lgpol(a)>>1) >degpol(b')
2061 : */
2062 :
2063 : GEN
2064 725970 : Flx_halfgcd_pre(GEN x, GEN y, ulong p, ulong pi)
2065 : {
2066 : pari_sp av;
2067 : GEN M,q,r;
2068 725970 : long lx=lgpol(x), ly=lgpol(y);
2069 725966 : if (!lx)
2070 : {
2071 0 : long v = x[1];
2072 0 : retmkmat2(mkcol2(pol0_Flx(v),pol1_Flx(v)),
2073 : mkcol2(pol1_Flx(v),pol0_Flx(v)));
2074 : }
2075 725966 : if (ly < lx) return Flx_halfgcd_i(x,y,p,pi);
2076 4477 : av = avma;
2077 4477 : q = Flx_divrem(y,x,p,&r);
2078 4477 : M = Flx_halfgcd_i(x,r,p,pi);
2079 4477 : gcoeff(M,1,1) = Flx_sub(gcoeff(M,1,1), Flx_mul_pre(q,gcoeff(M,1,2), p,pi), p);
2080 4477 : gcoeff(M,2,1) = Flx_sub(gcoeff(M,2,1), Flx_mul_pre(q,gcoeff(M,2,2), p,pi), p);
2081 4477 : return gerepilecopy(av, M);
2082 : }
2083 :
2084 : GEN
2085 112 : Flx_halfgcd(GEN x, GEN y, ulong p)
2086 112 : { return Flx_halfgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2087 :
2088 : /*Do not garbage collect*/
2089 : static GEN
2090 75375533 : Flx_gcd_basecase(GEN a, GEN b, ulong p, ulong pi)
2091 : {
2092 75375533 : pari_sp av = avma;
2093 75375533 : ulong iter = 0;
2094 75375533 : if (lg(b) > lg(a)) swap(a, b);
2095 263441229 : while (lgpol(b))
2096 : {
2097 187448839 : GEN c = Flx_rem_pre(a,b,p,pi);
2098 188065696 : iter++; a = b; b = c;
2099 188065696 : if (gc_needed(av,2))
2100 : {
2101 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (d = %ld)",degpol(c));
2102 0 : gerepileall(av,2, &a,&b);
2103 : }
2104 : }
2105 75291072 : return iter < 2 ? Flx_copy(a) : a;
2106 : }
2107 :
2108 : GEN
2109 76919595 : Flx_gcd_pre(GEN x, GEN y, ulong p, ulong pi)
2110 : {
2111 76919595 : pari_sp av = avma;
2112 : long lim;
2113 76919595 : if (!lgpol(x)) return Flx_copy(y);
2114 75370364 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2115 75387951 : while (lgpol(y) >= lim)
2116 : {
2117 : GEN c;
2118 9 : if (lgpol(y)<=(lgpol(x)>>1))
2119 : {
2120 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2121 0 : x = y; y = r;
2122 : }
2123 9 : c = FlxM_Flx_mul2(Flx_halfgcd_pre(x, y, p, pi), x, y, p, pi);
2124 9 : x = gel(c,1); y = gel(c,2);
2125 9 : if (gc_needed(av,2))
2126 : {
2127 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_gcd (y = %ld)",degpol(y));
2128 0 : gerepileall(av,2,&x,&y);
2129 : }
2130 : }
2131 75376353 : return gerepileuptoleaf(av, Flx_gcd_basecase(x,y,p,pi));
2132 : }
2133 : GEN
2134 27319961 : Flx_gcd(GEN x, GEN y, ulong p)
2135 27319961 : { return Flx_gcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2136 :
2137 : int
2138 6231368 : Flx_is_squarefree(GEN z, ulong p)
2139 : {
2140 6231368 : pari_sp av = avma;
2141 6231368 : GEN d = Flx_gcd(z, Flx_deriv(z,p) , p);
2142 6231347 : return gc_bool(av, degpol(d) == 0);
2143 : }
2144 :
2145 : static long
2146 139421 : Flx_is_smooth_squarefree(GEN f, long r, ulong p, ulong pi)
2147 : {
2148 139421 : pari_sp av = avma;
2149 : long i;
2150 139421 : GEN sx = polx_Flx(f[1]), a = sx;
2151 593027 : for(i=1;;i++)
2152 : {
2153 593027 : if (degpol(f)<=r) return gc_long(av,1);
2154 567226 : a = Flxq_powu_pre(Flx_rem_pre(a,f,p,pi), p, f, p, pi);
2155 571528 : if (Flx_equal(a, sx)) return gc_long(av,1);
2156 567335 : if (i==r) return gc_long(av,0);
2157 455691 : f = Flx_div_pre(f, Flx_gcd_pre(Flx_sub(a,sx,p),f,p,pi),p,pi);
2158 : }
2159 : }
2160 :
2161 : static long
2162 9090 : Flx_is_l_pow(GEN x, ulong p)
2163 : {
2164 9090 : ulong i, lx = lgpol(x);
2165 18284 : for (i=1; i<lx; i++)
2166 16337 : if (x[i+2] && i%p) return 0;
2167 1947 : return 1;
2168 : }
2169 :
2170 : int
2171 139393 : Flx_is_smooth_pre(GEN g, long r, ulong p, ulong pi)
2172 : {
2173 : while (1)
2174 9091 : {
2175 139393 : GEN f = Flx_gcd_pre(g, Flx_deriv(g, p), p, pi);
2176 139347 : if (!Flx_is_smooth_squarefree(Flx_div_pre(g, f, p, pi), r, p, pi))
2177 111613 : return 0;
2178 27918 : if (degpol(f)==0) return 1;
2179 9078 : g = Flx_is_l_pow(f,p) ? Flx_deflate(f, p): f;
2180 : }
2181 : }
2182 : int
2183 74256 : Flx_is_smooth(GEN g, long r, ulong p)
2184 74256 : { return Flx_is_smooth_pre(g, r, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2185 :
2186 : static GEN
2187 6281148 : Flx_extgcd_basecase(GEN a, GEN b, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2188 : {
2189 6281148 : pari_sp av=avma;
2190 : GEN u,v,d,d1,v1;
2191 6281148 : long vx = a[1];
2192 6281148 : d = a; d1 = b;
2193 6281148 : v = pol0_Flx(vx); v1 = pol1_Flx(vx);
2194 31636638 : while (lgpol(d1))
2195 : {
2196 25355242 : GEN r, q = Flx_divrem_pre(d,d1,p,pi, &r);
2197 25355504 : v = Flx_sub(v,Flx_mul_pre(q,v1,p,pi),p);
2198 25355815 : u=v; v=v1; v1=u;
2199 25355815 : u=r; d=d1; d1=u;
2200 25355815 : if (gc_needed(av,2))
2201 : {
2202 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_extgcd (d = %ld)",degpol(d));
2203 0 : gerepileall(av,5, &d,&d1,&u,&v,&v1);
2204 : }
2205 : }
2206 6279445 : if (ptu) *ptu = Flx_div_pre(Flx_sub(d, Flx_mul_pre(b,v,p,pi), p), a, p, pi);
2207 6281003 : *ptv = v; return d;
2208 : }
2209 :
2210 : static GEN
2211 6 : Flx_extgcd_halfgcd(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2212 : {
2213 6 : pari_sp av=avma;
2214 6 : GEN u,v,R = matid2_FlxM(x[1]);
2215 6 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2216 12 : while (lgpol(y) >= lim)
2217 : {
2218 : GEN M, c;
2219 6 : if (lgpol(y)<=(lgpol(x)>>1))
2220 : {
2221 0 : GEN r, q = Flx_divrem_pre(x, y, p, pi, &r);
2222 0 : x = y; y = r;
2223 0 : R = Flx_FlxM_qmul(q, R, p, pi);
2224 : }
2225 6 : M = Flx_halfgcd_pre(x,y, p, pi);
2226 6 : c = FlxM_Flx_mul2(M, x,y, p, pi);
2227 6 : R = FlxM_mul2(M, R, p, pi);
2228 6 : x = gel(c,1); y = gel(c,2);
2229 6 : gerepileall(av,3,&x,&y,&R);
2230 : }
2231 6 : y = Flx_extgcd_basecase(x,y,p,pi,&u,&v);
2232 6 : if (ptu) *ptu = Flx_addmulmul(u, v, gcoeff(R,1,1), gcoeff(R,2,1), p, pi);
2233 6 : *ptv = Flx_addmulmul(u, v, gcoeff(R,1,2), gcoeff(R,2,2), p, pi);
2234 6 : return y;
2235 : }
2236 :
2237 : /* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st
2238 : * ux + vy = gcd (mod p) */
2239 : GEN
2240 6281140 : Flx_extgcd_pre(GEN x, GEN y, ulong p, ulong pi, GEN *ptu, GEN *ptv)
2241 : {
2242 6281140 : pari_sp av = avma;
2243 : GEN d;
2244 6281140 : long lim = get_Fl_threshold(p, Flx_EXTGCD_LIMIT, Flx_EXTGCD2_LIMIT);
2245 6281168 : if (lgpol(y) >= lim)
2246 6 : d = Flx_extgcd_halfgcd(x, y, p, pi, ptu, ptv);
2247 : else
2248 6281144 : d = Flx_extgcd_basecase(x, y, p, pi, ptu, ptv);
2249 6281001 : return gc_all(av, ptu?3:2, &d, ptv, ptu);
2250 : }
2251 : GEN
2252 836840 : Flx_extgcd(GEN x, GEN y, ulong p, GEN *ptu, GEN *ptv)
2253 836840 : { return Flx_extgcd_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptu, ptv); }
2254 :
2255 : static GEN
2256 848 : Flx_halfres_pre(GEN a, GEN b, ulong p, ulong pi, ulong *r)
2257 : {
2258 : struct Flx_res res;
2259 : GEN M, V, A, B;
2260 : long dB;
2261 :
2262 848 : res.res = *r;
2263 848 : res.lc = Flx_lead(b);
2264 848 : res.deg0 = degpol(a);
2265 848 : res.deg1 = degpol(b);
2266 848 : res.off = 0;
2267 848 : M = Flx_halfres_i(a, b, p, pi, &res);
2268 848 : V = FlxM_Flx_mul2(M, a, b, p, pi);
2269 848 : A = gel(V,1); B = gel(V,2); dB = degpol(B);
2270 :
2271 848 : if (dB < degpol(b))
2272 848 : Flx_halfres_update_pre(res.deg0, res.deg1, dB, p, pi, &res);
2273 848 : *r = res.res;
2274 848 : return mkvec3(M,A,B);
2275 : }
2276 :
2277 : static ulong
2278 5940648 : Flx_resultant_basecase_pre(GEN a, GEN b, ulong p, ulong pi)
2279 : {
2280 : long da,db,dc,cnt;
2281 5940648 : ulong lb, res = 1UL;
2282 : pari_sp av;
2283 : GEN c;
2284 :
2285 5940648 : if (lgpol(a)==0 || lgpol(b)==0) return 0;
2286 5940628 : da = degpol(a);
2287 5940654 : db = degpol(b);
2288 5940652 : if (db > da)
2289 : {
2290 0 : swapspec(a,b, da,db);
2291 0 : if (both_odd(da,db)) res = p-res;
2292 : }
2293 5940652 : else if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */
2294 5940652 : cnt = 0; av = avma;
2295 50124072 : while (db)
2296 : {
2297 44195455 : lb = b[db+2];
2298 44195455 : c = Flx_rem_pre(a,b, p,pi);
2299 44012927 : a = b; b = c; dc = degpol(c);
2300 44006699 : if (dc < 0) return gc_long(av,0);
2301 :
2302 43999194 : if (both_odd(da,db)) res = p - res;
2303 44005487 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, da - dc, p, pi), p);
2304 44183115 : if (++cnt == 100) { cnt = 0; gerepileall(av, 2, &a, &b); }
2305 44183420 : da = db; /* = degpol(a) */
2306 44183420 : db = dc; /* = degpol(b) */
2307 : }
2308 5928617 : return gc_ulong(av, Fl_mul(res, Fl_powu_pre(b[2], da, p, pi), p));
2309 : }
2310 :
2311 : ulong
2312 5942158 : Flx_resultant_pre(GEN x, GEN y, ulong p, ulong pi)
2313 : {
2314 5942158 : pari_sp av = avma;
2315 : long lim;
2316 5942158 : ulong res = 1;
2317 5942158 : long dx = degpol(x), dy = degpol(y);
2318 5942138 : if (dx < 0 || dy < 0) return 0;
2319 5940696 : if (dx < dy)
2320 : {
2321 448537 : swap(x,y);
2322 448537 : if (both_odd(dx, dy))
2323 1906 : res = Fl_neg(res, p);
2324 : }
2325 5940697 : lim = get_Fl_threshold(p, Flx_GCD_LIMIT, Flx_GCD2_LIMIT);
2326 5941536 : while (lgpol(y) >= lim)
2327 : {
2328 : GEN c;
2329 848 : if (lgpol(y)<=(lgpol(x)>>1))
2330 : {
2331 0 : GEN r = Flx_rem_pre(x, y, p, pi);
2332 0 : long dx = degpol(x), dy = degpol(y), dr = degpol(r);
2333 0 : ulong ly = y[dy+2];
2334 0 : if (ly != 1) res = Fl_mul(res, Fl_powu_pre(ly, dx - dr, p, pi), p);
2335 0 : if (both_odd(dx, dy))
2336 0 : res = Fl_neg(res, p);
2337 0 : x = y; y = r;
2338 : }
2339 848 : c = Flx_halfres_pre(x, y, p, pi, &res);
2340 848 : x = gel(c,2); y = gel(c,3);
2341 848 : if (gc_needed(av,2))
2342 : {
2343 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flx_res (y = %ld)",degpol(y));
2344 0 : gerepileall(av,2,&x,&y);
2345 : }
2346 : }
2347 5940657 : return gc_ulong(av, Fl_mul(res, Flx_resultant_basecase_pre(x, y, p, pi), p));
2348 : }
2349 :
2350 : ulong
2351 4731912 : Flx_resultant(GEN a, GEN b, ulong p)
2352 4731912 : { return Flx_resultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2353 :
2354 : /* If resultant is 0, *ptU and *ptV are not set */
2355 : ulong
2356 0 : Flx_extresultant_pre(GEN a, GEN b, ulong p, ulong pi, GEN *ptU, GEN *ptV)
2357 : {
2358 0 : GEN z,q,u,v, x = a, y = b;
2359 0 : ulong lb, res = 1UL;
2360 0 : pari_sp av = avma;
2361 : long dx, dy, dz;
2362 0 : long vs=a[1];
2363 :
2364 0 : dx = degpol(x);
2365 0 : dy = degpol(y);
2366 0 : if (dy > dx)
2367 : {
2368 0 : swap(x,y); lswap(dx,dy); pswap(ptU, ptV);
2369 0 : a = x; b = y;
2370 0 : if (both_odd(dx,dy)) res = p-res;
2371 : }
2372 : /* dy <= dx */
2373 0 : if (dy < 0) return 0;
2374 :
2375 0 : u = pol0_Flx(vs);
2376 0 : v = pol1_Flx(vs); /* v = 1 */
2377 0 : while (dy)
2378 : { /* b u = x (a), b v = y (a) */
2379 0 : lb = y[dy+2];
2380 0 : q = Flx_divrem_pre(x,y, p, pi, &z);
2381 0 : x = y; y = z; /* (x,y) = (y, x - q y) */
2382 0 : dz = degpol(z); if (dz < 0) return gc_ulong(av,0);
2383 0 : z = Flx_sub(u, Flx_mul_pre(q,v, p, pi), p);
2384 0 : u = v; v = z; /* (u,v) = (v, u - q v) */
2385 :
2386 0 : if (both_odd(dx,dy)) res = p - res;
2387 0 : if (lb != 1) res = Fl_mul(res, Fl_powu_pre(lb, dx-dz, p, pi), p);
2388 0 : dx = dy; /* = degpol(x) */
2389 0 : dy = dz; /* = degpol(y) */
2390 : }
2391 0 : res = Fl_mul(res, Fl_powu_pre(y[2], dx, p, pi), p);
2392 0 : lb = Fl_mul(res, Fl_inv(y[2],p), p);
2393 0 : v = gerepileuptoleaf(av, Flx_Fl_mul(v, lb, p));
2394 0 : av = avma;
2395 0 : u = Flx_sub(Fl_to_Flx(res,vs), Flx_mul_pre(b,v,p,pi), p);
2396 0 : u = gerepileuptoleaf(av, Flx_div_pre(u,a,p,pi)); /* = (res - b v) / a */
2397 0 : *ptU = u;
2398 0 : *ptV = v; return res;
2399 : }
2400 :
2401 : ulong
2402 0 : Flx_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV)
2403 0 : { return Flx_extresultant_pre(a, b, p, SMALL_ULONG(p)? 0: get_Fl_red(p), ptU, ptV); }
2404 :
2405 : /* allow pi = 0 (SMALL_ULONG) */
2406 : ulong
2407 41618163 : Flx_eval_powers_pre(GEN x, GEN y, ulong p, ulong pi)
2408 : {
2409 41618163 : ulong l0, l1, h0, h1, v1, i = 1, lx = lg(x)-1;
2410 :
2411 41618163 : if (lx == 1) return 0;
2412 38940504 : x++;
2413 38940504 : if (pi)
2414 : {
2415 : LOCAL_OVERFLOW;
2416 : LOCAL_HIREMAINDER;
2417 38880274 : l1 = mulll(uel(x,i), uel(y,i)); h1 = hiremainder; v1 = 0;
2418 93579534 : while (++i < lx)
2419 : {
2420 54699260 : l0 = mulll(uel(x,i), uel(y,i)); h0 = hiremainder;
2421 54699260 : l1 = addll(l0, l1); h1 = addllx(h0, h1); v1 += overflow;
2422 : }
2423 80859 : return v1? remlll_pre(v1, h1, l1, p, pi)
2424 38961133 : : remll_pre(h1, l1, p, pi);
2425 : }
2426 : else
2427 : {
2428 60230 : l1 = x[i] * y[i];
2429 30921419 : while (++i < lx) { l1 += x[i] * y[i]; if (l1 & HIGHBIT) l1 %= p; }
2430 60230 : return l1 % p;
2431 : }
2432 : }
2433 :
2434 : /* allow pi = 0 (SMALL_ULONG) */
2435 : ulong
2436 34438692 : Flx_eval_pre(GEN x, ulong y, ulong p, ulong pi)
2437 : {
2438 34438692 : long i, n = degpol(x);
2439 : ulong t;
2440 34437564 : if (n <= 0) return n? 0: x[2];
2441 27058668 : if (n > 15)
2442 : {
2443 180418 : pari_sp av = avma;
2444 180418 : GEN v = Fl_powers_pre(y, n, p, pi);
2445 180464 : return gc_ulong(av, Flx_eval_powers_pre(x, v, p, pi));
2446 : }
2447 26878250 : i = n+2; t = x[i];
2448 26878250 : if (pi)
2449 : {
2450 103520058 : for (i--; i>=2; i--) t = Fl_addmul_pre(uel(x, i), t, y, p, pi);
2451 25776869 : return t;
2452 : }
2453 2633565 : for (i--; i>=2; i--) t = (t * y + x[i]) % p;
2454 1100018 : return t %= p;
2455 : }
2456 : ulong
2457 20352905 : Flx_eval(GEN x, ulong y, ulong p)
2458 20352905 : { return Flx_eval_pre(x, y, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2459 :
2460 : ulong
2461 3255 : Flv_prod_pre(GEN x, ulong p, ulong pi)
2462 : {
2463 3255 : pari_sp ltop = avma;
2464 : GEN v;
2465 3255 : long i,k,lx = lg(x);
2466 3255 : if (lx == 1) return 1UL;
2467 3255 : if (lx == 2) return uel(x,1);
2468 3003 : v = cgetg(1+(lx << 1), t_VECSMALL);
2469 3003 : k = 1;
2470 28593 : for (i=1; i<lx-1; i+=2)
2471 25590 : uel(v,k++) = Fl_mul_pre(uel(x,i), uel(x,i+1), p, pi);
2472 3003 : if (i < lx) uel(v,k++) = uel(x,i);
2473 13529 : while (k > 2)
2474 : {
2475 10526 : lx = k; k = 1;
2476 36116 : for (i=1; i<lx-1; i+=2)
2477 25590 : uel(v,k++) = Fl_mul_pre(uel(v,i), uel(v,i+1), p, pi);
2478 10526 : if (i < lx) uel(v,k++) = uel(v,i);
2479 : }
2480 3003 : return gc_ulong(ltop, uel(v,1));
2481 : }
2482 :
2483 : ulong
2484 0 : Flv_prod(GEN v, ulong p)
2485 : {
2486 0 : return Flv_prod_pre(v, p, get_Fl_red(p));
2487 : }
2488 :
2489 : GEN
2490 0 : FlxV_prod(GEN V, ulong p)
2491 : {
2492 : struct _Flxq D;
2493 0 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2494 0 : return gen_product(V, (void *)&D, &_Flx_mul);
2495 : }
2496 :
2497 : /* compute prod (x - a[i]) */
2498 : GEN
2499 671756 : Flv_roots_to_pol(GEN a, ulong p, long vs)
2500 : {
2501 : struct _Flxq D;
2502 671756 : long i,k,lx = lg(a);
2503 : GEN p1;
2504 671756 : if (lx == 1) return pol1_Flx(vs);
2505 671756 : p1 = cgetg(lx, t_VEC);
2506 11443796 : for (k=1,i=1; i<lx-1; i+=2)
2507 10772269 : gel(p1,k++) = mkvecsmall4(vs, Fl_mul(a[i], a[i+1], p),
2508 10772540 : Fl_neg(Fl_add(a[i],a[i+1],p),p), 1);
2509 671256 : if (i < lx)
2510 54497 : gel(p1,k++) = mkvecsmall3(vs, Fl_neg(a[i],p), 1);
2511 671251 : D.T = NULL; D.aut = NULL; D.p = p; D.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2512 671247 : setlg(p1, k); return gen_product(p1, (void *)&D, _Flx_mul);
2513 : }
2514 :
2515 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for "large" p */
2516 : INLINE void
2517 16301371 : Flv_inv_pre_indir(GEN w, GEN v, ulong p, ulong pi)
2518 : {
2519 16301371 : pari_sp av = avma;
2520 16301371 : long n = lg(w), i;
2521 : ulong u;
2522 : GEN c;
2523 :
2524 16301371 : if (n == 1) return;
2525 16301371 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2526 72038075 : for (i = 2; i < n; ++i) c[i] = Fl_mul_pre(w[i], c[i-1], p, pi);
2527 16459404 : i = n-1; u = Fl_inv(c[i], p);
2528 72409707 : for ( ; i > 1; --i)
2529 : {
2530 55979952 : ulong t = Fl_mul_pre(u, c[i-1], p, pi);
2531 55966412 : u = Fl_mul_pre(u, w[i], p, pi); v[i] = t;
2532 : }
2533 16429755 : v[1] = u; set_avma(av);
2534 : }
2535 :
2536 : void
2537 15955121 : Flv_inv_pre_inplace(GEN v, ulong p, ulong pi) { Flv_inv_pre_indir(v,v, p, pi); }
2538 :
2539 : GEN
2540 10864 : Flv_inv_pre(GEN w, ulong p, ulong pi)
2541 10864 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_pre_indir(w, v, p, pi); return v; }
2542 :
2543 : /* set v[i] = w[i]^{-1}; may be called with w = v, suitable for SMALL_ULONG p */
2544 : INLINE void
2545 45739 : Flv_inv_indir(GEN w, GEN v, ulong p)
2546 : {
2547 45739 : pari_sp av = avma;
2548 45739 : long n = lg(w), i;
2549 : ulong u;
2550 : GEN c;
2551 :
2552 45739 : if (n == 1) return;
2553 45739 : c = cgetg(n, t_VECSMALL); c[1] = w[1];
2554 1673407 : for (i = 2; i < n; ++i) c[i] = Fl_mul(w[i], c[i-1], p);
2555 45741 : i = n-1; u = Fl_inv(c[i], p);
2556 1673415 : for ( ; i > 1; --i)
2557 : {
2558 1627675 : ulong t = Fl_mul(u, c[i-1], p);
2559 1627667 : u = Fl_mul(u, w[i], p); v[i] = t;
2560 : }
2561 45740 : v[1] = u; set_avma(av);
2562 : }
2563 : static void
2564 373201 : Flv_inv_i(GEN v, GEN w, ulong p)
2565 : {
2566 373201 : if (SMALL_ULONG(p)) Flv_inv_indir(w, v, p);
2567 327463 : else Flv_inv_pre_indir(w, v, p, get_Fl_red(p));
2568 373203 : }
2569 : void
2570 12017 : Flv_inv_inplace(GEN v, ulong p) { Flv_inv_i(v, v, p); }
2571 : GEN
2572 361186 : Flv_inv(GEN w, ulong p)
2573 361186 : { GEN v = cgetg(lg(w), t_VECSMALL); Flv_inv_i(v, w, p); return v; }
2574 :
2575 : GEN
2576 31964122 : Flx_div_by_X_x(GEN a, ulong x, ulong p, ulong *rem)
2577 : {
2578 31964122 : long l = lg(a), i;
2579 : GEN a0, z0, z;
2580 31964122 : if (l <= 3)
2581 : {
2582 0 : if (rem) *rem = l == 2? 0: a[2];
2583 0 : return zero_Flx(a[1]);
2584 : }
2585 31964122 : z = cgetg(l-1,t_VECSMALL); z[1] = a[1];
2586 31785621 : a0 = a + l-1;
2587 31785621 : z0 = z + l-2; *z0 = *a0--;
2588 31785621 : if (SMALL_ULONG(p))
2589 : {
2590 76814707 : for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */
2591 : {
2592 57101910 : ulong t = (*a0-- + x * *z0--) % p;
2593 57101910 : *z0 = (long)t;
2594 : }
2595 19712797 : if (rem) *rem = (*a0 + x * *z0) % p;
2596 : }
2597 : else
2598 : {
2599 47814392 : for (i=l-3; i>1; i--)
2600 : {
2601 35672956 : ulong t = Fl_add((ulong)*a0--, Fl_mul(x, *z0--, p), p);
2602 35741568 : *z0 = (long)t;
2603 : }
2604 12141436 : if (rem) *rem = Fl_add((ulong)*a0, Fl_mul(x, *z0, p), p);
2605 : }
2606 31894031 : return z;
2607 : }
2608 :
2609 : /* xa, ya = t_VECSMALL */
2610 : static GEN
2611 362428 : Flv_producttree(GEN xa, GEN s, ulong p, long vs)
2612 : {
2613 362428 : long n = lg(xa)-1;
2614 362428 : long m = n==1 ? 1: expu(n-1)+1;
2615 362428 : long i, j, k, ls = lg(s);
2616 : ulong pi;
2617 362428 : GEN T = cgetg(m+1, t_VEC);
2618 362422 : GEN t = cgetg(ls, t_VEC);
2619 3639532 : for (j=1, k=1; j<ls; k+=s[j++])
2620 3276922 : gel(t, j) = s[j] == 1 ?
2621 3277105 : mkvecsmall3(vs, Fl_neg(xa[k], p), 1):
2622 1064626 : mkvecsmall4(vs, Fl_mul(xa[k], xa[k+1], p),
2623 1064630 : Fl_neg(Fl_add(xa[k],xa[k+1],p),p), 1);
2624 362427 : gel(T,1) = t; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2625 1098808 : for (i=2; i<=m; i++)
2626 : {
2627 736460 : GEN u = gel(T, i-1);
2628 736460 : long n = lg(u)-1;
2629 736460 : GEN t = cgetg(((n+1)>>1)+1, t_VEC);
2630 3650610 : for (j=1, k=1; k<n; j++, k+=2)
2631 2914229 : gel(t, j) = Flx_mul_pre(gel(u, k), gel(u, k+1), p, pi);
2632 736381 : gel(T, i) = t;
2633 : }
2634 362348 : return T;
2635 : }
2636 :
2637 : static GEN
2638 402732 : Flx_Flv_multieval_tree(GEN P, GEN xa, GEN T, ulong p)
2639 : {
2640 : long i,j,k;
2641 402732 : long m = lg(T)-1;
2642 402732 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2643 402731 : GEN R = cgetg(lg(xa), t_VECSMALL);
2644 402725 : GEN Tp = cgetg(m+1, t_VEC), t;
2645 402718 : gel(Tp, m) = mkvec(P);
2646 1324512 : for (i=m-1; i>=1; i--)
2647 : {
2648 921789 : GEN u = gel(T, i), v = gel(Tp, i+1);
2649 921789 : long n = lg(u)-1;
2650 921789 : t = cgetg(n+1, t_VEC);
2651 4865838 : for (j=1, k=1; k<n; j++, k+=2)
2652 : {
2653 3944065 : gel(t, k) = Flx_rem_pre(gel(v, j), gel(u, k), p, pi);
2654 3944057 : gel(t, k+1) = Flx_rem_pre(gel(v, j), gel(u, k+1), p, pi);
2655 : }
2656 921773 : gel(Tp, i) = t;
2657 : }
2658 : {
2659 402723 : GEN u = gel(T, i+1), v = gel(Tp, i+1);
2660 402723 : long n = lg(u)-1;
2661 4751927 : for (j=1, k=1; j<=n; j++)
2662 : {
2663 4349174 : long c, d = degpol(gel(u,j));
2664 10012110 : for (c=1; c<=d; c++, k++) R[k] = Flx_eval_pre(gel(v, j), xa[k], p, pi);
2665 : }
2666 402753 : return gc_const((pari_sp)R, R);
2667 : }
2668 : }
2669 :
2670 : static GEN
2671 1077421 : FlvV_polint_tree(GEN T, GEN R, GEN s, GEN xa, GEN ya, ulong p, long vs)
2672 : {
2673 1077421 : pari_sp av = avma;
2674 1077421 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
2675 1077406 : long m = lg(T)-1;
2676 1077406 : long i, j, k, ls = lg(s);
2677 1077406 : GEN Tp = cgetg(m+1, t_VEC);
2678 1076820 : GEN t = cgetg(ls, t_VEC);
2679 20240608 : for (j=1, k=1; j<ls; k+=s[j++])
2680 19164035 : if (s[j]==2)
2681 : {
2682 6320037 : ulong a = Fl_mul(ya[k], R[k], p);
2683 6322081 : ulong b = Fl_mul(ya[k+1], R[k+1], p);
2684 6332082 : gel(t, j) = mkvecsmall3(vs, Fl_neg(Fl_add(Fl_mul(xa[k], b, p ),
2685 6327035 : Fl_mul(xa[k+1], a, p), p), p), Fl_add(a, b, p));
2686 6323841 : gel(t, j) = Flx_renormalize(gel(t, j), 4);
2687 : }
2688 : else
2689 12843998 : gel(t, j) = Fl_to_Flx(Fl_mul(ya[k], R[k], p), vs);
2690 1076573 : gel(Tp, 1) = t;
2691 4933919 : for (i=2; i<=m; i++)
2692 : {
2693 3857407 : GEN u = gel(T, i-1);
2694 3857407 : GEN t = cgetg(lg(gel(T,i)), t_VEC);
2695 3854080 : GEN v = gel(Tp, i-1);
2696 3854080 : long n = lg(v)-1;
2697 21900926 : for (j=1, k=1; k<n; j++, k+=2)
2698 18056287 : gel(t, j) = Flx_add(Flx_mul_pre(gel(u, k), gel(v, k+1), p, pi),
2699 18043580 : Flx_mul_pre(gel(u, k+1), gel(v, k), p, pi), p);
2700 3857346 : gel(Tp, i) = t;
2701 : }
2702 1076512 : return gerepileuptoleaf(av, gmael(Tp,m,1));
2703 : }
2704 :
2705 : GEN
2706 0 : Flx_Flv_multieval(GEN P, GEN xa, ulong p)
2707 : {
2708 0 : pari_sp av = avma;
2709 0 : GEN s = producttree_scheme(lg(xa)-1);
2710 0 : GEN T = Flv_producttree(xa, s, p, P[1]);
2711 0 : return gerepileuptoleaf(av, Flx_Flv_multieval_tree(P, xa, T, p));
2712 : }
2713 :
2714 : static GEN
2715 2471 : FlxV_Flv_multieval_tree(GEN x, GEN xa, GEN T, ulong p)
2716 45248 : { pari_APPLY_same(Flx_Flv_multieval_tree(gel(x,i), xa, T, p)) }
2717 :
2718 : GEN
2719 2471 : FlxV_Flv_multieval(GEN P, GEN xa, ulong p)
2720 : {
2721 2471 : pari_sp av = avma;
2722 2471 : GEN s = producttree_scheme(lg(xa)-1);
2723 2471 : GEN T = Flv_producttree(xa, s, p, P[1]);
2724 2471 : return gerepileupto(av, FlxV_Flv_multieval_tree(P, xa, T, p));
2725 : }
2726 :
2727 : GEN
2728 110642 : Flv_polint(GEN xa, GEN ya, ulong p, long vs)
2729 : {
2730 110642 : pari_sp av = avma;
2731 110642 : GEN s = producttree_scheme(lg(xa)-1);
2732 110643 : GEN T = Flv_producttree(xa, s, p, vs);
2733 110642 : long m = lg(T)-1;
2734 110642 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2735 110643 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2736 110644 : return gerepileuptoleaf(av, FlvV_polint_tree(T, R, s, xa, ya, p, vs));
2737 : }
2738 :
2739 : GEN
2740 96387 : Flv_Flm_polint(GEN xa, GEN ya, ulong p, long vs)
2741 : {
2742 96387 : pari_sp av = avma;
2743 96387 : GEN s = producttree_scheme(lg(xa)-1);
2744 96386 : GEN T = Flv_producttree(xa, s, p, vs);
2745 96384 : long i, m = lg(T)-1, l = lg(ya)-1;
2746 96384 : GEN P = Flx_deriv(gmael(T, m, 1), p);
2747 96383 : GEN R = Flv_inv(Flx_Flv_multieval_tree(P, xa, T, p), p);
2748 96383 : GEN M = cgetg(l+1, t_VEC);
2749 1062996 : for (i=1; i<=l; i++)
2750 966619 : gel(M,i) = FlvV_polint_tree(T, R, s, xa, gel(ya,i), p, vs);
2751 96377 : return gerepileupto(av, M);
2752 : }
2753 :
2754 : GEN
2755 152929 : Flv_invVandermonde(GEN L, ulong den, ulong p)
2756 : {
2757 152929 : pari_sp av = avma;
2758 152929 : long i, n = lg(L);
2759 : GEN M, R;
2760 152929 : GEN s = producttree_scheme(n-1);
2761 152929 : GEN tree = Flv_producttree(L, s, p, 0);
2762 152929 : long m = lg(tree)-1;
2763 152929 : GEN T = gmael(tree, m, 1);
2764 152929 : R = Flv_inv(Flx_Flv_multieval_tree(Flx_deriv(T, p), L, tree, p), p);
2765 152929 : if (den!=1) R = Flv_Fl_mul(R, den, p);
2766 152929 : M = cgetg(n, t_MAT);
2767 599895 : for (i = 1; i < n; i++)
2768 : {
2769 446966 : GEN P = Flx_Fl_mul(Flx_div_by_X_x(T, uel(L,i), p, NULL), uel(R,i), p);
2770 446966 : gel(M,i) = Flx_to_Flv(P, n-1);
2771 : }
2772 152929 : return gerepilecopy(av, M);
2773 : }
2774 :
2775 : /***********************************************************************/
2776 : /** Flxq **/
2777 : /***********************************************************************/
2778 : /* Flxq objects are Flx modulo another Flx called q. */
2779 :
2780 : /* Product of y and x in Z/pZ[X]/(T), as t_VECSMALL. */
2781 : GEN
2782 187735999 : Flxq_mul_pre(GEN x,GEN y,GEN T,ulong p,ulong pi)
2783 187735999 : { return Flx_rem_pre(Flx_mul_pre(x,y,p,pi),T,p,pi); }
2784 : GEN
2785 13142050 : Flxq_mul(GEN x,GEN y,GEN T,ulong p)
2786 13142050 : { return Flxq_mul_pre(x,y,T,p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2787 :
2788 : GEN
2789 271380270 : Flxq_sqr_pre(GEN x,GEN T,ulong p,ulong pi)
2790 271380270 : { return Flx_rem_pre(Flx_sqr_pre(x, p,pi), T, p,pi); }
2791 : /* Square of y in Z/pZ[X]/(T), as t_VECSMALL. */
2792 : GEN
2793 2712612 : Flxq_sqr(GEN x,GEN T,ulong p)
2794 2712612 : { return Flxq_sqr_pre(x,T,p,SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2795 :
2796 : static GEN
2797 1724600 : _Flxq_red(void *E, GEN x)
2798 1724600 : { struct _Flxq *s = (struct _Flxq *)E;
2799 1724600 : return Flx_rem_pre(x, s->T, s->p, s->pi); }
2800 : #if 0
2801 : static GEN
2802 : _Flx_sub(void *E, GEN x, GEN y)
2803 : { struct _Flxq *s = (struct _Flxq *)E;
2804 : return Flx_sub(x,y,s->p); }
2805 : #endif
2806 : static GEN
2807 264468613 : _Flxq_sqr(void *data, GEN x)
2808 : {
2809 264468613 : struct _Flxq *D = (struct _Flxq*)data;
2810 264468613 : return Flxq_sqr_pre(x, D->T, D->p, D->pi);
2811 : }
2812 : static GEN
2813 147140055 : _Flxq_mul(void *data, GEN x, GEN y)
2814 : {
2815 147140055 : struct _Flxq *D = (struct _Flxq*)data;
2816 147140055 : return Flxq_mul_pre(x,y, D->T, D->p, D->pi);
2817 : }
2818 : static GEN
2819 21903598 : _Flxq_one(void *data)
2820 : {
2821 21903598 : struct _Flxq *D = (struct _Flxq*)data;
2822 21903598 : return pol1_Flx(get_Flx_var(D->T));
2823 : }
2824 :
2825 : static GEN
2826 21631286 : _Flxq_powu_i(struct _Flxq *D, GEN x, ulong n)
2827 21631286 : { return gen_powu_i(x, n, (void*)D, &_Flxq_sqr, &_Flxq_mul); }
2828 : static GEN
2829 68 : _Flxq_powu(struct _Flxq *D, GEN x, ulong n)
2830 68 : { pari_sp av = avma; return gerepileuptoleaf(av, _Flxq_powu_i(D, x, n)); }
2831 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2832 : GEN
2833 22345790 : Flxq_powu_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
2834 : {
2835 : pari_sp av;
2836 : struct _Flxq D;
2837 22345790 : switch(n)
2838 : {
2839 0 : case 0: return pol1_Flx(get_Flx_var(T));
2840 135411 : case 1: return Flx_copy(x);
2841 587454 : case 2: return Flxq_sqr_pre(x, T, p, pi);
2842 : }
2843 21622925 : av = avma; set_Flxq_pre(&D, T, p, pi);
2844 21631250 : return gerepileuptoleaf(av, _Flxq_powu_i(&D, x, n));
2845 : }
2846 : GEN
2847 489847 : Flxq_powu(GEN x, ulong n, GEN T, ulong p)
2848 489847 : { return Flxq_powu_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2849 :
2850 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL. */
2851 : GEN
2852 26859752 : Flxq_pow_pre(GEN x, GEN n, GEN T, ulong p, ulong pi)
2853 : {
2854 26859752 : pari_sp av = avma;
2855 : struct _Flxq D;
2856 : GEN y;
2857 26859752 : long s = signe(n);
2858 26859752 : if (!s) return pol1_Flx(get_Flx_var(T));
2859 26598894 : if (s < 0) x = Flxq_inv_pre(x,T,p,pi);
2860 26598888 : if (is_pm1(n)) return s < 0 ? x : Flx_copy(x);
2861 25354382 : set_Flxq_pre(&D, T, p, pi);
2862 25354434 : y = gen_pow_i(x, n, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2863 25354290 : return gerepileuptoleaf(av, y);
2864 : }
2865 : GEN
2866 1014963 : Flxq_pow(GEN x, GEN n, GEN T, ulong p)
2867 1014963 : { return Flxq_pow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2868 :
2869 : GEN
2870 28 : Flxq_pow_init_pre(GEN x, GEN n, long k, GEN T, ulong p, ulong pi)
2871 : {
2872 28 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2873 28 : return gen_pow_init(x, n, k, (void*)&D, &_Flxq_sqr, &_Flxq_mul);
2874 : }
2875 : GEN
2876 0 : Flxq_pow_init(GEN x, GEN n, long k, GEN T, ulong p)
2877 0 : { return Flxq_pow_init_pre(x, n, k, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2878 :
2879 : GEN
2880 4400 : Flxq_pow_table_pre(GEN R, GEN n, GEN T, ulong p, ulong pi)
2881 : {
2882 4400 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2883 4400 : return gen_pow_table(R, n, (void*)&D, &_Flxq_one, &_Flxq_mul);
2884 : }
2885 : GEN
2886 0 : Flxq_pow_table(GEN R, GEN n, GEN T, ulong p)
2887 0 : { return Flxq_pow_table_pre(R, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2888 :
2889 : /* Inverse of x in Z/lZ[X]/(T) or NULL if inverse doesn't exist
2890 : * not stack clean. */
2891 : GEN
2892 5444336 : Flxq_invsafe_pre(GEN x, GEN T, ulong p, ulong pi)
2893 : {
2894 5444336 : GEN V, z = Flx_extgcd_pre(get_Flx_mod(T), x, p, pi, NULL, &V);
2895 : ulong iz;
2896 5444426 : if (degpol(z)) return NULL;
2897 5443735 : iz = Fl_inv(uel(z,2), p);
2898 5443757 : return Flx_Fl_mul(V, iz, p);
2899 : }
2900 : GEN
2901 418532 : Flxq_invsafe(GEN x, GEN T, ulong p)
2902 418532 : { return Flxq_invsafe_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2903 :
2904 : GEN
2905 4624441 : Flxq_inv_pre(GEN x, GEN T, ulong p, ulong pi)
2906 : {
2907 4624441 : pari_sp av=avma;
2908 4624441 : GEN U = Flxq_invsafe_pre(x, T, p, pi);
2909 4624393 : if (!U) pari_err_INV("Flxq_inv",Flx_to_ZX(x));
2910 4624365 : return gerepileuptoleaf(av, U);
2911 : }
2912 : GEN
2913 345415 : Flxq_inv(GEN x, GEN T, ulong p)
2914 345415 : { return Flxq_inv_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2915 :
2916 : GEN
2917 2183427 : Flxq_div_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
2918 : {
2919 2183427 : pari_sp av = avma;
2920 2183427 : return gerepileuptoleaf(av, Flxq_mul_pre(x,Flxq_inv_pre(y,T,p,pi),T,p,pi));
2921 : }
2922 : GEN
2923 220847 : Flxq_div(GEN x, GEN y, GEN T, ulong p)
2924 220847 : { return Flxq_div_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2925 :
2926 : GEN
2927 21905262 : Flxq_powers_pre(GEN x, long l, GEN T, ulong p, ulong pi)
2928 : {
2929 21905262 : int use_sqr = 2*degpol(x) >= get_Flx_degree(T);
2930 21894066 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
2931 21892433 : return gen_powers(x, l, use_sqr, (void*)&D, &_Flxq_sqr, &_Flxq_mul, &_Flxq_one);
2932 : }
2933 : GEN
2934 231804 : Flxq_powers(GEN x, long l, GEN T, ulong p)
2935 231804 : { return Flxq_powers_pre(x, l, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2936 :
2937 : GEN
2938 168497 : Flxq_matrix_pow_pre(GEN y, long n, long m, GEN P, ulong l, ulong li)
2939 168497 : { return FlxV_to_Flm(Flxq_powers_pre(y,m-1,P,l,li),n); }
2940 : GEN
2941 399 : Flxq_matrix_pow(GEN y, long n, long m, GEN P, ulong l)
2942 399 : { return Flxq_matrix_pow_pre(y, n, m, P, l, SMALL_ULONG(l)? 0: get_Fl_red(l)); }
2943 :
2944 : GEN
2945 12476697 : Flx_Frobenius_pre(GEN T, ulong p, ulong pi)
2946 12476697 : { return Flxq_powu_pre(polx_Flx(get_Flx_var(T)), p, T, p, pi); }
2947 : GEN
2948 82296 : Flx_Frobenius(GEN T, ulong p)
2949 82296 : { return Flx_Frobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2950 :
2951 : GEN
2952 85501 : Flx_matFrobenius_pre(GEN T, ulong p, ulong pi)
2953 : {
2954 85501 : long n = get_Flx_degree(T);
2955 85501 : return Flxq_matrix_pow_pre(Flx_Frobenius_pre(T, p, pi), n, n, T, p, pi);
2956 : }
2957 : GEN
2958 0 : Flx_matFrobenius(GEN T, ulong p)
2959 0 : { return Flx_matFrobenius_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
2960 :
2961 : static GEN
2962 12436160 : Flx_blocks_Flm(GEN P, long n, long m)
2963 : {
2964 12436160 : GEN z = cgetg(m+1,t_MAT);
2965 12433525 : long i,j, k=2, l = lg(P);
2966 37234185 : for(i=1; i<=m; i++)
2967 : {
2968 24804938 : GEN zi = cgetg(n+1,t_VECSMALL);
2969 24800660 : gel(z,i) = zi;
2970 111560655 : for(j=1; j<=n; j++)
2971 86759995 : uel(zi, j) = k==l ? 0 : uel(P,k++);
2972 : }
2973 12429247 : return z;
2974 : }
2975 :
2976 : GEN
2977 516470 : Flx_blocks(GEN P, long n, long m)
2978 : {
2979 516470 : GEN z = cgetg(m+1,t_VEC);
2980 516255 : long i,j, k=2, l = lg(P);
2981 1546681 : for(i=1; i<=m; i++)
2982 : {
2983 1030702 : GEN zi = cgetg(n+2,t_VECSMALL);
2984 1029713 : zi[1] = P[1];
2985 1029713 : gel(z,i) = zi;
2986 6454393 : for(j=2; j<n+2; j++)
2987 5424680 : uel(zi, j) = k==l ? 0 : uel(P,k++);
2988 1029713 : zi = Flx_renormalize(zi, n+2);
2989 : }
2990 515979 : return z;
2991 : }
2992 :
2993 : static GEN
2994 12436348 : FlxV_to_Flm_lg(GEN x, long m, long n)
2995 : {
2996 : long i;
2997 12436348 : GEN y = cgetg(n+1, t_MAT);
2998 57264571 : for (i=1; i<=n; i++) gel(y,i) = Flx_to_Flv(gel(x,i), m);
2999 12432002 : return y;
3000 : }
3001 :
3002 : /* allow pi = 0 (SMALL_ULONG) */
3003 : GEN
3004 12633603 : Flx_FlxqV_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3005 : {
3006 12633603 : pari_sp btop, av = avma;
3007 12633603 : long sv = get_Flx_var(T), m = get_Flx_degree(T);
3008 12633333 : long i, l = lg(x)-1, lQ = lgpol(Q), n, d;
3009 : GEN A, B, C, S, g;
3010 12634444 : if (lQ == 0) return pol0_Flx(sv);
3011 12437640 : if (lQ <= l)
3012 : {
3013 5344714 : n = l;
3014 5344714 : d = 1;
3015 : }
3016 : else
3017 : {
3018 7092926 : n = l-1;
3019 7092926 : d = (lQ+n-1)/n;
3020 : }
3021 12437640 : A = FlxV_to_Flm_lg(x, m, n);
3022 12435113 : B = Flx_blocks_Flm(Q, n, d);
3023 12432871 : C = gerepileupto(av, Flm_mul(A, B, p));
3024 12438519 : g = gel(x, l);
3025 12438519 : if (pi && SMALL_ULONG(p)) pi = 0;
3026 12438519 : T = Flx_get_red_pre(T, p, pi);
3027 12437192 : btop = avma;
3028 12437192 : S = Flv_to_Flx(gel(C, d), sv);
3029 24813945 : for (i = d-1; i>0; i--)
3030 : {
3031 12379769 : S = Flx_add(Flxq_mul_pre(S, g, T, p, pi), Flv_to_Flx(gel(C,i), sv), p);
3032 12378942 : if (gc_needed(btop,1))
3033 0 : S = gerepileuptoleaf(btop, S);
3034 : }
3035 12434176 : return gerepileuptoleaf(av, S);
3036 : }
3037 : GEN
3038 5103 : Flx_FlxqV_eval(GEN Q, GEN x, GEN T, ulong p)
3039 5103 : { return Flx_FlxqV_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3040 :
3041 : /* allow pi = 0 (SMALL_ULONG) */
3042 : GEN
3043 2484122 : Flx_Flxq_eval_pre(GEN Q, GEN x, GEN T, ulong p, ulong pi)
3044 : {
3045 2484122 : pari_sp av = avma;
3046 : GEN z, V;
3047 2484122 : long d = degpol(Q), rtd;
3048 2484120 : if (d < 0) return pol0_Flx(get_Flx_var(T));
3049 2484029 : rtd = (long) sqrt((double)d);
3050 2484029 : T = Flx_get_red_pre(T, p, pi);
3051 2484040 : V = Flxq_powers_pre(x, rtd, T, p, pi);
3052 2484079 : z = Flx_FlxqV_eval_pre(Q, V, T, p, pi);
3053 2484070 : return gerepileupto(av, z);
3054 : }
3055 : GEN
3056 787559 : Flx_Flxq_eval(GEN Q, GEN x, GEN T, ulong p)
3057 787559 : { return Flx_Flxq_eval_pre(Q, x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3058 :
3059 : /* allow pi = 0 (SMALL_ULONG) */
3060 : GEN
3061 0 : FlxC_FlxqV_eval_pre(GEN x, GEN v, GEN T, ulong p, ulong pi)
3062 0 : { pari_APPLY_type(t_COL, Flx_FlxqV_eval_pre(gel(x,i), v, T, p, pi)) }
3063 : GEN
3064 0 : FlxC_FlxqV_eval(GEN x, GEN v, GEN T, ulong p)
3065 0 : { return FlxC_FlxqV_eval_pre(x, v, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3066 :
3067 : /* allow pi = 0 (SMALL_ULONG) */
3068 : GEN
3069 0 : FlxC_Flxq_eval_pre(GEN x, GEN F, GEN T, ulong p, ulong pi)
3070 : {
3071 0 : long d = brent_kung_optpow(get_Flx_degree(T)-1,lg(x)-1,1);
3072 0 : GEN Fp = Flxq_powers_pre(F, d, T, p, pi);
3073 0 : return FlxC_FlxqV_eval_pre(x, Fp, T, p, pi);
3074 : }
3075 : GEN
3076 0 : FlxC_Flxq_eval(GEN x, GEN F, GEN T, ulong p)
3077 0 : { return FlxC_Flxq_eval_pre(x, F, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3078 :
3079 : #if 0
3080 : static struct bb_algebra Flxq_algebra = { _Flxq_red, _Flx_add, _Flx_sub,
3081 : _Flxq_mul, _Flxq_sqr, _Flxq_one, _Flxq_zero};
3082 : #endif
3083 :
3084 : static GEN
3085 389109 : Flxq_autpow_sqr(void *E, GEN x)
3086 : {
3087 389109 : struct _Flxq *D = (struct _Flxq*)E;
3088 389109 : return Flx_Flxq_eval_pre(x, x, D->T, D->p, D->pi);
3089 : }
3090 : static GEN
3091 20734 : Flxq_autpow_msqr(void *E, GEN x)
3092 : {
3093 20734 : struct _Flxq *D = (struct _Flxq*)E;
3094 20734 : return Flx_FlxqV_eval_pre(Flxq_autpow_sqr(E, x), D->aut, D->T, D->p, D->pi);
3095 : }
3096 :
3097 : GEN
3098 307189 : Flxq_autpow_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3099 : {
3100 307189 : pari_sp av = avma;
3101 : struct _Flxq D;
3102 : long d;
3103 307189 : if (n==0) return Flx_rem_pre(polx_Flx(x[1]), T, p, pi);
3104 307182 : if (n==1) return Flx_rem_pre(x, T, p, pi);
3105 306692 : set_Flxq_pre(&D, T, p, pi);
3106 306692 : d = brent_kung_optpow(get_Flx_degree(T), hammingl(n)-1, 1);
3107 306692 : D.aut = Flxq_powers_pre(x, d, T, p, D.pi);
3108 306692 : x = gen_powu_fold_i(x,n,(void*)&D,Flxq_autpow_sqr,Flxq_autpow_msqr);
3109 306692 : return gerepilecopy(av, x);
3110 : }
3111 : GEN
3112 7 : Flxq_autpow(GEN x, ulong n, GEN T, ulong p)
3113 7 : { return Flxq_autpow_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3114 :
3115 : GEN
3116 1705 : Flxq_autpowers(GEN x, ulong l, GEN T, ulong p)
3117 : {
3118 1705 : long d, vT = get_Flx_var(T), dT = get_Flx_degree(T);
3119 : ulong i, pi;
3120 1705 : pari_sp av = avma;
3121 1705 : GEN xp, V = cgetg(l+2,t_VEC);
3122 1705 : gel(V,1) = polx_Flx(vT); if (l==0) return V;
3123 1705 : gel(V,2) = gcopy(x); if (l==1) return V;
3124 1705 : pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3125 1705 : T = Flx_get_red_pre(T, p, pi);
3126 1705 : d = brent_kung_optpow(dT-1, l-1, 1);
3127 1705 : xp = Flxq_powers_pre(x, d, T, p, pi);
3128 7188 : for(i = 3; i < l+2; i++)
3129 5483 : gel(V,i) = Flx_FlxqV_eval_pre(gel(V,i-1), xp, T, p, pi);
3130 1705 : return gerepilecopy(av, V);
3131 : }
3132 :
3133 : static GEN
3134 591678 : Flxq_autsum_mul(void *E, GEN x, GEN y)
3135 : {
3136 591678 : struct _Flxq *D = (struct _Flxq*)E;
3137 591678 : GEN T = D->T;
3138 591678 : ulong p = D->p, pi = D->pi;
3139 591678 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3140 591678 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3141 591678 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3142 591678 : GEN V2 = Flxq_powers_pre(phi2, d, T, p, pi);
3143 591678 : GEN phi3 = Flx_FlxqV_eval_pre(phi1, V2, T, p, pi);
3144 591678 : GEN aphi = Flx_FlxqV_eval_pre(a1, V2, T, p, pi);
3145 591678 : GEN a3 = Flxq_mul_pre(aphi, a2, T, p, pi);
3146 591678 : return mkvec2(phi3, a3);
3147 : }
3148 : static GEN
3149 381076 : Flxq_autsum_sqr(void *E, GEN x)
3150 381076 : { return Flxq_autsum_mul(E, x, x); }
3151 :
3152 : static GEN
3153 317636 : Flxq_autsum_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3154 : {
3155 317636 : pari_sp av = avma;
3156 317636 : struct _Flxq D; set_Flxq_pre(&D, T, p, pi);
3157 317636 : x = gen_powu_i(x,n,(void*)&D,Flxq_autsum_sqr,Flxq_autsum_mul);
3158 317636 : return gerepilecopy(av, x);
3159 : }
3160 : GEN
3161 0 : Flxq_autsum(GEN x, ulong n, GEN T, ulong p)
3162 0 : { return Flxq_autsum_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3163 :
3164 : static GEN
3165 674336 : Flxq_auttrace_mul(void *E, GEN x, GEN y)
3166 : {
3167 674336 : struct _Flxq *D = (struct _Flxq*)E;
3168 674336 : GEN T = D->T;
3169 674336 : ulong p = D->p, pi = D->pi;
3170 674336 : GEN phi1 = gel(x,1), a1 = gel(x,2);
3171 674336 : GEN phi2 = gel(y,1), a2 = gel(y,2);
3172 674336 : ulong d = brent_kung_optpow(maxss(degpol(phi1),degpol(a1)),2,1);
3173 674350 : GEN V1 = Flxq_powers_pre(phi1, d, T, p, pi);
3174 674296 : GEN phi3 = Flx_FlxqV_eval_pre(phi2, V1, T, p, pi);
3175 674332 : GEN aphi = Flx_FlxqV_eval_pre(a2, V1, T, p, pi);
3176 674333 : GEN a3 = Flx_add(a1, aphi, p);
3177 674332 : return mkvec2(phi3, a3);
3178 : }
3179 :
3180 : static GEN
3181 548646 : Flxq_auttrace_sqr(void *E, GEN x)
3182 548646 : { return Flxq_auttrace_mul(E, x, x); }
3183 :
3184 : GEN
3185 819970 : Flxq_auttrace_pre(GEN x, ulong n, GEN T, ulong p, ulong pi)
3186 : {
3187 819970 : pari_sp av = avma;
3188 : struct _Flxq D;
3189 819970 : set_Flxq_pre(&D, T, p, pi);
3190 819974 : x = gen_powu_i(x,n,(void*)&D,Flxq_auttrace_sqr,Flxq_auttrace_mul);
3191 819957 : return gerepilecopy(av, x);
3192 : }
3193 : GEN
3194 0 : Flxq_auttrace(GEN x, ulong n, GEN T, ulong p)
3195 0 : { return Flxq_auttrace_pre(x, n, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3196 :
3197 : static long
3198 862204 : bounded_order(ulong p, GEN b, long k)
3199 : {
3200 862204 : GEN a = modii(utoipos(p), b);
3201 : long i;
3202 2089981 : for(i = 1; i < k; i++)
3203 : {
3204 1614933 : if (equali1(a)) return i;
3205 1227780 : a = modii(muliu(a,p),b);
3206 : }
3207 475048 : return 0;
3208 : }
3209 :
3210 : /*
3211 : n = (p^d-a)\b
3212 : b = bb*p^vb
3213 : p^k = 1 [bb]
3214 : d = m*k+r+vb
3215 : u = (p^k-1)/bb;
3216 : v = (p^(r+vb)-a)/b;
3217 : w = (p^(m*k)-1)/(p^k-1)
3218 : n = p^r*w*u+v
3219 : w*u = p^vb*(p^(m*k)-1)/b
3220 : n = p^(r+vb)*(p^(m*k)-1)/b+(p^(r+vb)-a)/b
3221 : */
3222 :
3223 : static GEN
3224 26012361 : Flxq_pow_Frobenius(GEN x, GEN n, GEN aut, GEN T, ulong p, ulong pi)
3225 : {
3226 26012361 : pari_sp av=avma;
3227 26012361 : long d = get_Flx_degree(T);
3228 26012382 : GEN an = absi_shallow(n), z, q;
3229 26012358 : if (abscmpiu(an,p)<0 || cmpis(an,d)<=0) return Flxq_pow_pre(x, n, T, p, pi);
3230 862580 : q = powuu(p, d);
3231 862579 : if (dvdii(q, n))
3232 : {
3233 353 : long vn = logint(an, utoipos(p));
3234 353 : GEN autvn = vn==1 ? aut: Flxq_autpow_pre(aut,vn,T,p,pi);
3235 353 : z = Flx_Flxq_eval_pre(x,autvn,T,p,pi);
3236 : } else
3237 : {
3238 862227 : GEN b = diviiround(q, an), a = subii(q, mulii(an,b));
3239 : GEN bb, u, v, autk;
3240 862226 : long vb = Z_lvalrem(b,p,&bb);
3241 862227 : long m, r, k = is_pm1(bb)? 1: bounded_order(p,bb,d);
3242 862224 : if (!k || d-vb < k) return Flxq_pow_pre(x,n, T,p,pi);
3243 387169 : m = (d-vb)/k; r = (d-vb)%k;
3244 387169 : u = diviiexact(subiu(powuu(p,k),1),bb);
3245 387169 : v = diviiexact(subii(powuu(p,r+vb),a),b);
3246 387169 : autk = k==1 ? aut: Flxq_autpow_pre(aut,k,T,p,pi);
3247 387169 : if (r)
3248 : {
3249 95471 : GEN autr = r==1 ? aut: Flxq_autpow_pre(aut,r,T,p,pi);
3250 95471 : z = Flx_Flxq_eval_pre(x,autr,T,p,pi);
3251 291698 : } else z = x;
3252 387169 : if (m > 1) z = gel(Flxq_autsum_pre(mkvec2(autk, z), m, T, p, pi), 2);
3253 387169 : if (!is_pm1(u)) z = Flxq_pow_pre(z, u, T, p, pi);
3254 387169 : if (signe(v)) z = Flxq_mul_pre(z, Flxq_pow_pre(x, v, T, p, pi), T, p, pi);
3255 : }
3256 387522 : return gerepileupto(av,signe(n)>0 ? z : Flxq_inv_pre(z,T,p,pi));
3257 : }
3258 :
3259 : static GEN
3260 26004872 : _Flxq_pow(void *data, GEN x, GEN n)
3261 : {
3262 26004872 : struct _Flxq *D = (struct _Flxq*)data;
3263 26004872 : return Flxq_pow_Frobenius(x, n, D->aut, D->T, D->p, D->pi);
3264 : }
3265 :
3266 : static GEN
3267 363815 : _Flxq_rand(void *data)
3268 : {
3269 363815 : pari_sp av=avma;
3270 363815 : struct _Flxq *D = (struct _Flxq*)data;
3271 : GEN z;
3272 : do
3273 : {
3274 367214 : set_avma(av);
3275 367214 : z = random_Flx(get_Flx_degree(D->T),get_Flx_var(D->T),D->p);
3276 367215 : } while (lgpol(z)==0);
3277 363816 : return z;
3278 : }
3279 :
3280 : /* discrete log in FpXQ for a in Fp^*, g in FpXQ^* of order ord */
3281 : static GEN
3282 28970 : Fl_Flxq_log(ulong a, GEN g, GEN o, GEN T, ulong p)
3283 : {
3284 28970 : pari_sp av = avma;
3285 : GEN q,n_q,ord,ordp, op;
3286 :
3287 28970 : if (a == 1UL) return gen_0;
3288 : /* p > 2 */
3289 :
3290 28970 : ordp = utoi(p - 1);
3291 28970 : ord = get_arith_Z(o);
3292 28970 : if (!ord) ord = T? subiu(powuu(p, get_FpX_degree(T)), 1): ordp;
3293 28970 : if (a == p - 1) /* -1 */
3294 7501 : return gerepileuptoint(av, shifti(ord,-1));
3295 21469 : ordp = gcdii(ordp, ord);
3296 21469 : op = typ(o)==t_MAT ? famat_Z_gcd(o, ordp) : ordp;
3297 :
3298 21469 : q = NULL;
3299 21469 : if (T)
3300 : { /* we want < g > = Fp^* */
3301 21469 : if (!equalii(ord,ordp)) {
3302 10566 : q = diviiexact(ord,ordp);
3303 10566 : g = Flxq_pow(g,q,T,p);
3304 : }
3305 : }
3306 21469 : n_q = Fp_log(utoi(a), utoipos(uel(g,2)), op, utoipos(p));
3307 21469 : if (lg(n_q)==1) return gerepileuptoleaf(av, n_q);
3308 21469 : if (q) n_q = mulii(q, n_q);
3309 21469 : return gerepileuptoint(av, n_q);
3310 : }
3311 :
3312 : static GEN
3313 692488 : Flxq_easylog(void* E, GEN a, GEN g, GEN ord)
3314 : {
3315 692488 : struct _Flxq *f = (struct _Flxq *)E;
3316 692488 : GEN T = f->T;
3317 692488 : ulong p = f->p;
3318 692488 : long d = get_Flx_degree(T);
3319 692488 : if (Flx_equal1(a)) return gen_0;
3320 534252 : if (Flx_equal(a,g)) return gen_1;
3321 164852 : if (!degpol(a))
3322 28970 : return Fl_Flxq_log(uel(a,2), g, ord, T, p);
3323 135882 : if (typ(ord)!=t_INT || d <= 4 || d == 6 || abscmpiu(ord,1UL<<27)<0)
3324 135854 : return NULL;
3325 28 : return Flxq_log_index(a, g, ord, T, p);
3326 : }
3327 :
3328 : static const struct bb_group Flxq_star={_Flxq_mul,_Flxq_pow,_Flxq_rand,hash_GEN,Flx_equal,Flx_equal1,Flxq_easylog};
3329 :
3330 : const struct bb_group *
3331 441087 : get_Flxq_star(void **E, GEN T, ulong p)
3332 : {
3333 441087 : struct _Flxq *e = (struct _Flxq *) stack_malloc(sizeof(struct _Flxq));
3334 441086 : e->T = T; e->p = p; e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3335 441086 : e->aut = Flx_Frobenius_pre(T, p, e->pi);
3336 441088 : *E = (void*)e; return &Flxq_star;
3337 : }
3338 :
3339 : GEN
3340 95497 : Flxq_order(GEN a, GEN ord, GEN T, ulong p)
3341 : {
3342 : void *E;
3343 95497 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3344 95497 : return gen_order(a,ord,E,S);
3345 : }
3346 :
3347 : GEN
3348 158384 : Flxq_log(GEN a, GEN g, GEN ord, GEN T, ulong p)
3349 : {
3350 : void *E;
3351 158384 : pari_sp av = avma;
3352 158384 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3353 158384 : GEN v = get_arith_ZZM(ord), F = gmael(v,2,1);
3354 158384 : if (Flxq_log_use_index(gel(F,lg(F)-1), T, p))
3355 24590 : v = mkvec2(gel(v, 1), ZM_famat_limit(gel(v, 2), int2n(27)));
3356 158384 : return gerepileuptoleaf(av, gen_PH_log(a, g, v, E, S));
3357 : }
3358 :
3359 : GEN
3360 190273 : Flxq_sqrtn(GEN a, GEN n, GEN T, ulong p, GEN *zeta)
3361 : {
3362 190273 : if (!lgpol(a))
3363 : {
3364 3066 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
3365 3059 : if (zeta)
3366 0 : *zeta=pol1_Flx(get_Flx_var(T));
3367 3059 : return pol0_Flx(get_Flx_var(T));
3368 : }
3369 : else
3370 : {
3371 : void *E;
3372 187207 : pari_sp av = avma;
3373 187207 : const struct bb_group *S = get_Flxq_star(&E,T,p);
3374 187207 : GEN o = subiu(powuu(p,get_Flx_degree(T)), 1);
3375 187202 : GEN s = gen_Shanks_sqrtn(a,n,o,zeta,E,S);
3376 187207 : if (!s) return gc_NULL(av);
3377 184365 : return gc_all(av, zeta?2:1, &s, zeta);
3378 : }
3379 : }
3380 :
3381 : GEN
3382 164388 : Flxq_sqrt(GEN a, GEN T, ulong p)
3383 164388 : { return Flxq_sqrtn(a, gen_2, T, p, NULL); }
3384 :
3385 : /* assume T irreducible mod p */
3386 : int
3387 364126 : Flxq_issquare(GEN x, GEN T, ulong p)
3388 : {
3389 364126 : if (lgpol(x) == 0 || p == 2) return 1;
3390 361011 : return krouu(Flxq_norm(x,T,p), p) == 1;
3391 : }
3392 :
3393 : /* assume T irreducible mod p */
3394 : int
3395 0 : Flxq_is2npower(GEN x, long n, GEN T, ulong p)
3396 : {
3397 : pari_sp av;
3398 : GEN m;
3399 0 : if (n==1) return Flxq_issquare(x, T, p);
3400 0 : if (lgpol(x) == 0 || p == 2) return 1;
3401 0 : av = avma;
3402 0 : m = shifti(subiu(powuu(p, get_Flx_degree(T)), 1), -n);
3403 0 : return gc_bool(av, Flx_equal1(Flxq_pow(x, m, T, p)));
3404 : }
3405 :
3406 : GEN
3407 113589 : Flxq_lroot_fast_pre(GEN a, GEN sqx, GEN T, long p, ulong pi)
3408 : {
3409 113589 : pari_sp av=avma;
3410 113589 : GEN A = Flx_splitting(a,p);
3411 113589 : return gerepileuptoleaf(av, FlxqV_dotproduct_pre(A,sqx,T,p,pi));
3412 : }
3413 : GEN
3414 0 : Flxq_lroot_fast(GEN a, GEN sqx, GEN T, long p)
3415 0 : { return Flxq_lroot_fast_pre(a, sqx, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3416 :
3417 : GEN
3418 25053 : Flxq_lroot_pre(GEN a, GEN T, long p, ulong pi)
3419 : {
3420 25053 : pari_sp av=avma;
3421 25053 : long n = get_Flx_degree(T), d = degpol(a);
3422 : GEN sqx, V;
3423 25053 : if (n==1) return leafcopy(a);
3424 25053 : if (n==2) return Flxq_powu_pre(a, p, T, p, pi);
3425 25053 : sqx = Flxq_autpow_pre(Flx_Frobenius_pre(T, p, pi), n-1, T, p, pi);
3426 25053 : if (d==1 && a[2]==0 && a[3]==1) return gerepileuptoleaf(av, sqx);
3427 0 : if (d>=p)
3428 : {
3429 0 : V = Flxq_powers_pre(sqx,p-1,T,p,pi);
3430 0 : return gerepileuptoleaf(av, Flxq_lroot_fast_pre(a,V,T,p,pi));
3431 : } else
3432 0 : return gerepileuptoleaf(av, Flx_Flxq_eval_pre(a,sqx,T,p,pi));
3433 : }
3434 : GEN
3435 0 : Flxq_lroot(GEN a, GEN T, long p)
3436 0 : { return Flxq_lroot_pre(a, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3437 :
3438 : ulong
3439 397773 : Flxq_norm(GEN x, GEN TB, ulong p)
3440 : {
3441 397773 : GEN T = get_Flx_mod(TB);
3442 397773 : ulong y = Flx_resultant(T, x, p), L = Flx_lead(T);
3443 397773 : if (L==1 || lgpol(x)==0) return y;
3444 0 : return Fl_div(y, Fl_powu(L, (ulong)degpol(x), p), p);
3445 : }
3446 :
3447 : ulong
3448 3359 : Flxq_trace(GEN x, GEN TB, ulong p)
3449 : {
3450 3359 : pari_sp av = avma;
3451 : ulong t;
3452 3359 : GEN T = get_Flx_mod(TB);
3453 3359 : long n = degpol(T)-1;
3454 3359 : GEN z = Flxq_mul(x, Flx_deriv(T, p), TB, p);
3455 3359 : t = degpol(z)<n ? 0 : Fl_div(z[2+n],T[3+n],p);
3456 3359 : return gc_ulong(av, t);
3457 : }
3458 :
3459 : /*x must be reduced*/
3460 : GEN
3461 3626 : Flxq_charpoly(GEN x, GEN TB, ulong p)
3462 : {
3463 3626 : pari_sp ltop=avma;
3464 3626 : GEN T = get_Flx_mod(TB);
3465 3626 : long vs = evalvarn(fetch_var());
3466 3626 : GEN xm1 = deg1pol_shallow(pol1_Flx(x[1]),Flx_neg(x,p),vs);
3467 3626 : GEN r = Flx_FlxY_resultant(T, xm1, p);
3468 3626 : r[1] = x[1];
3469 3626 : (void)delete_var(); return gerepileupto(ltop, r);
3470 : }
3471 :
3472 : /* Computing minimal polynomial : */
3473 : /* cf Shoup 'Efficient Computation of Minimal Polynomials */
3474 : /* in Algebraic Extensions of Finite Fields' */
3475 :
3476 : /* Let v a linear form, return the linear form z->v(tau*z)
3477 : that is, v*(M_tau) */
3478 :
3479 : static GEN
3480 1442909 : Flxq_transmul_init(GEN tau, GEN T, ulong p, ulong pi)
3481 : {
3482 : GEN bht;
3483 1442909 : GEN h, Tp = get_Flx_red(T, &h);
3484 1442896 : long n = degpol(Tp), vT = Tp[1];
3485 1442896 : GEN ft = Flx_recipspec(Tp+2, n+1, n+1);
3486 1442879 : GEN bt = Flx_recipspec(tau+2, lgpol(tau), n);
3487 1442877 : ft[1] = vT; bt[1] = vT;
3488 1442877 : if (h)
3489 2688 : bht = Flxn_mul_pre(bt, h, n-1, p, pi);
3490 : else
3491 : {
3492 1440189 : GEN bh = Flx_div_pre(Flx_shift(tau, n-1), T, p, pi);
3493 1440183 : bht = Flx_recipspec(bh+2, lgpol(bh), n-1);
3494 1440185 : bht[1] = vT;
3495 : }
3496 1442873 : return mkvec3(bt, bht, ft);
3497 : }
3498 :
3499 : static GEN
3500 3489614 : Flxq_transmul(GEN tau, GEN a, long n, ulong p, ulong pi)
3501 : {
3502 3489614 : pari_sp ltop = avma;
3503 : GEN t1, t2, t3, vec;
3504 3489614 : GEN bt = gel(tau, 1), bht = gel(tau, 2), ft = gel(tau, 3);
3505 3489614 : if (lgpol(a)==0) return pol0_Flx(a[1]);
3506 3461924 : t2 = Flx_shift(Flx_mul_pre(bt, a, p, pi),1-n);
3507 3461640 : if (lgpol(bht)==0) return gerepileuptoleaf(ltop, t2);
3508 2617872 : t1 = Flx_shift(Flx_mul_pre(ft, a, p, pi),-n);
3509 2617916 : t3 = Flxn_mul_pre(t1, bht, n-1, p, pi);
3510 2617867 : vec = Flx_sub(t2, Flx_shift(t3, 1), p);
3511 2617869 : return gerepileuptoleaf(ltop, vec);
3512 : }
3513 :
3514 : GEN
3515 668359 : Flxq_minpoly_pre(GEN x, GEN T, ulong p, ulong pi)
3516 : {
3517 668359 : pari_sp ltop = avma;
3518 668359 : long vT = get_Flx_var(T), n = get_Flx_degree(T);
3519 : GEN v_x;
3520 668351 : GEN g = pol1_Flx(vT), tau = pol1_Flx(vT);
3521 668326 : T = Flx_get_red_pre(T, p, pi);
3522 668325 : v_x = Flxq_powers_pre(x, usqrt(2*n), T, p, pi);
3523 1389772 : while (lgpol(tau) != 0)
3524 : {
3525 : long i, j, m, k1;
3526 : GEN M, v, tr, g_prime, c;
3527 721434 : if (degpol(g) == n) { tau = pol1_Flx(vT); g = pol1_Flx(vT); }
3528 721438 : v = random_Flx(n, vT, p);
3529 721470 : tr = Flxq_transmul_init(tau, T, p, pi);
3530 721441 : v = Flxq_transmul(tr, v, n, p, pi);
3531 721465 : m = 2*(n-degpol(g));
3532 721465 : k1 = usqrt(m);
3533 721466 : tr = Flxq_transmul_init(gel(v_x,k1+1), T, p, pi);
3534 721441 : c = cgetg(m+2,t_VECSMALL);
3535 721367 : c[1] = vT;
3536 3489468 : for (i=0; i<m; i+=k1)
3537 : {
3538 2767997 : long mj = minss(m-i, k1);
3539 11229952 : for (j=0; j<mj; j++)
3540 8461478 : uel(c,m+1-(i+j)) = Flx_dotproduct_pre(v, gel(v_x,j+1), p, pi);
3541 2768474 : v = Flxq_transmul(tr, v, n, p, pi);
3542 : }
3543 721471 : c = Flx_renormalize(c, m+2);
3544 : /* now c contains <v,x^i> , i = 0..m-1 */
3545 721465 : M = Flx_halfgcd_pre(monomial_Flx(1, m, vT), c, p, pi);
3546 721478 : g_prime = gmael(M, 2, 2);
3547 721478 : if (degpol(g_prime) < 1) continue;
3548 710780 : g = Flx_mul_pre(g, g_prime, p, pi);
3549 710764 : tau = Flxq_mul_pre(tau, Flx_FlxqV_eval_pre(g_prime, v_x, T,p,pi), T,p,pi);
3550 : }
3551 668305 : g = Flx_normalize(g,p);
3552 668354 : return gerepileuptoleaf(ltop,g);
3553 : }
3554 : GEN
3555 40590 : Flxq_minpoly(GEN x, GEN T, ulong p)
3556 40590 : { return Flxq_minpoly_pre(x, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3557 :
3558 : GEN
3559 20 : Flxq_conjvec(GEN x, GEN T, ulong p)
3560 : {
3561 20 : long i, l = 1+get_Flx_degree(T);
3562 20 : GEN z = cgetg(l,t_COL);
3563 20 : struct _Flxq D; set_Flxq(&D, T, p);
3564 20 : gel(z,1) = Flx_copy(x);
3565 88 : for (i=2; i<l; i++) gel(z,i) = _Flxq_powu(&D, gel(z,i-1), p);
3566 20 : return z;
3567 : }
3568 :
3569 : GEN
3570 7166 : gener_Flxq(GEN T, ulong p, GEN *po)
3571 : {
3572 7166 : long i, j, vT = get_Flx_var(T), f = get_Flx_degree(T);
3573 : ulong p_1, pi;
3574 : GEN g, L, L2, o, q, F;
3575 : pari_sp av0, av;
3576 :
3577 7166 : if (f == 1) {
3578 : GEN fa;
3579 28 : o = utoipos(p-1);
3580 28 : fa = Z_factor(o);
3581 28 : L = gel(fa,1);
3582 28 : L = vecslice(L, 2, lg(L)-1); /* remove 2 for efficiency */
3583 28 : g = Fl_to_Flx(pgener_Fl_local(p, vec_to_vecsmall(L)), vT);
3584 28 : if (po) *po = mkvec2(o, fa);
3585 28 : return g;
3586 : }
3587 :
3588 7138 : av0 = avma; p_1 = p - 1;
3589 7138 : q = diviuexact(subiu(powuu(p,f), 1), p_1);
3590 :
3591 7138 : L = cgetg(1, t_VECSMALL);
3592 7138 : if (p > 3)
3593 : {
3594 2343 : ulong t = p_1 >> vals(p_1);
3595 2343 : GEN P = gel(factoru(t), 1);
3596 2343 : L = cgetg_copy(P, &i);
3597 3738 : while (--i) L[i] = p_1 / P[i];
3598 : }
3599 7138 : o = factor_pn_1(utoipos(p),f);
3600 7138 : L2 = leafcopy( gel(o, 1) );
3601 19093 : for (i = j = 1; i < lg(L2); i++)
3602 : {
3603 11955 : if (umodui(p_1, gel(L2,i)) == 0) continue;
3604 6460 : gel(L2,j++) = diviiexact(q, gel(L2,i));
3605 : }
3606 7138 : setlg(L2, j); pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3607 7138 : F = Flx_Frobenius_pre(T, p, pi);
3608 17737 : for (av = avma;; set_avma(av))
3609 10599 : {
3610 : GEN tt;
3611 17737 : g = random_Flx(f, vT, p);
3612 17737 : if (degpol(g) < 1) continue;
3613 12089 : if (p == 2) tt = g;
3614 : else
3615 : {
3616 8890 : ulong t = Flxq_norm(g, T, p);
3617 8890 : if (t == 1 || !is_gener_Fl(t, p, p_1, L)) continue;
3618 4782 : tt = Flxq_powu_pre(g, p_1>>1, T, p, pi);
3619 : }
3620 14576 : for (i = 1; i < j; i++)
3621 : {
3622 7438 : GEN a = Flxq_pow_Frobenius(tt, gel(L2,i), F, T, p, pi);
3623 7438 : if (!degpol(a) && uel(a,2) == p_1) break;
3624 : }
3625 7981 : if (i == j) break;
3626 : }
3627 7138 : if (!po)
3628 : {
3629 180 : set_avma((pari_sp)g);
3630 180 : g = gerepileuptoleaf(av0, g);
3631 : }
3632 : else {
3633 6958 : *po = mkvec2(subiu(powuu(p,f), 1), o);
3634 6958 : gerepileall(av0, 2, &g, po);
3635 : }
3636 7138 : return g;
3637 : }
3638 :
3639 : static GEN
3640 557727 : _Flxq_neg(void *E, GEN x)
3641 557727 : { struct _Flxq *s = (struct _Flxq *)E;
3642 557727 : return Flx_neg(x,s->p); }
3643 :
3644 : static GEN
3645 1528908 : _Flxq_rmul(void *E, GEN x, GEN y)
3646 1528908 : { struct _Flxq *s = (struct _Flxq *)E;
3647 1528908 : return Flx_mul_pre(x,y,s->p,s->pi); }
3648 :
3649 : static GEN
3650 19035 : _Flxq_inv(void *E, GEN x)
3651 19035 : { struct _Flxq *s = (struct _Flxq *)E;
3652 19035 : return Flxq_inv(x,s->T,s->p); }
3653 :
3654 : static int
3655 164940 : _Flxq_equal0(GEN x) { return lgpol(x)==0; }
3656 :
3657 : static GEN
3658 24793 : _Flxq_s(void *E, long x)
3659 24793 : { struct _Flxq *s = (struct _Flxq *)E;
3660 24793 : ulong u = x<0 ? s->p+x: (ulong)x;
3661 24793 : return Fl_to_Flx(u, get_Flx_var(s->T));
3662 : }
3663 :
3664 : static const struct bb_field Flxq_field={_Flxq_red,_Flx_add,_Flxq_rmul,_Flxq_neg,
3665 : _Flxq_inv,_Flxq_equal0,_Flxq_s};
3666 :
3667 67922 : const struct bb_field *get_Flxq_field(void **E, GEN T, ulong p)
3668 : {
3669 67922 : GEN z = new_chunk(sizeof(struct _Flxq));
3670 67922 : set_Flxq((struct _Flxq *)z, T, p); *E = (void*)z; return &Flxq_field;
3671 : }
3672 :
3673 : /***********************************************************************/
3674 : /** Flxn **/
3675 : /***********************************************************************/
3676 :
3677 : GEN
3678 54342 : Flx_invLaplace(GEN x, ulong p)
3679 : {
3680 54342 : long i, d = degpol(x);
3681 : ulong t;
3682 : GEN y;
3683 54333 : if (d <= 1) return Flx_copy(x);
3684 54333 : t = Fl_inv(factorial_Fl(d, p), p);
3685 54384 : y = cgetg(d+3, t_VECSMALL);
3686 54326 : y[1] = x[1];
3687 1328834 : for (i=d; i>=2; i--)
3688 : {
3689 1274464 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3690 1274492 : t = Fl_mul(t, i, p);
3691 : }
3692 54370 : uel(y,3) = uel(x,3);
3693 54370 : uel(y,2) = uel(x,2);
3694 54370 : return y;
3695 : }
3696 :
3697 : GEN
3698 27333 : Flx_Laplace(GEN x, ulong p)
3699 : {
3700 27333 : long i, d = degpol(x);
3701 27330 : ulong t = 1;
3702 : GEN y;
3703 27330 : if (d <= 1) return Flx_copy(x);
3704 27330 : y = cgetg(d+3, t_VECSMALL);
3705 27319 : y[1] = x[1];
3706 27319 : uel(y,2) = uel(x,2);
3707 27319 : uel(y,3) = uel(x,3);
3708 757670 : for (i=2; i<=d; i++)
3709 : {
3710 730324 : t = Fl_mul(t, i%p, p);
3711 730336 : uel(y,i+2) = Fl_mul(uel(x,i+2), t, p);
3712 : }
3713 27346 : return y;
3714 : }
3715 :
3716 : GEN
3717 5731170 : Flxn_red(GEN a, long n)
3718 : {
3719 5731170 : long i, L, l = lg(a);
3720 : GEN b;
3721 5731170 : if (l == 2 || !n) return zero_Flx(a[1]);
3722 5342108 : L = n+2; if (L > l) L = l;
3723 5342108 : b = cgetg(L, t_VECSMALL); b[1] = a[1];
3724 51678557 : for (i=2; i<L; i++) b[i] = a[i];
3725 5338442 : return Flx_renormalize(b,L);
3726 : }
3727 :
3728 : GEN
3729 4620668 : Flxn_mul_pre(GEN a, GEN b, long n, ulong p, ulong pi)
3730 4620668 : { return Flxn_red(Flx_mul_pre(a, b, p, pi), n); }
3731 : GEN
3732 75203 : Flxn_mul(GEN a, GEN b, long n, ulong p)
3733 75203 : { return Flxn_mul_pre(a, b, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3734 :
3735 : GEN
3736 0 : Flxn_sqr_pre(GEN a, long n, ulong p, ulong pi)
3737 0 : { return Flxn_red(Flx_sqr_pre(a, p, pi), n); }
3738 : GEN
3739 0 : Flxn_sqr(GEN a, long n, ulong p)
3740 0 : { return Flxn_sqr_pre(a, n, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3741 :
3742 : /* (f*g) \/ x^n */
3743 : static GEN
3744 937741 : Flx_mulhigh_i(GEN f, GEN g, long n, ulong p, ulong pi)
3745 937741 : { return Flx_shift(Flx_mul_pre(f, g, p, pi),-n); }
3746 :
3747 : static GEN
3748 516414 : Flxn_mulhigh(GEN f, GEN g, long n2, long n, ulong p, ulong pi)
3749 : {
3750 516414 : GEN F = Flx_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
3751 516099 : return Flx_add(Flx_mulhigh_i(fl, g, n2, p, pi),
3752 : Flxn_mul_pre(fh, g, n - n2, p, pi), p);
3753 : }
3754 :
3755 : /* g==NULL -> assume g==1 */
3756 : GEN
3757 55152 : Flxn_div_pre(GEN g, GEN f, long e, ulong p, ulong pi)
3758 : {
3759 55152 : pari_sp av = avma, av2;
3760 : ulong mask;
3761 : GEN W;
3762 55152 : long n = 1;
3763 55152 : if (lg(f) <= 2) pari_err_INV("Flxn_inv",f);
3764 55152 : W = Fl_to_Flx(Fl_inv(uel(f,2),p), f[1]);
3765 55146 : mask = quadratic_prec_mask(e);
3766 55162 : av2 = avma;
3767 258769 : for (;mask>1;)
3768 : {
3769 : GEN u, fr;
3770 203593 : long n2 = n;
3771 203593 : n<<=1; if (mask & 1) n--;
3772 203593 : mask >>= 1;
3773 203593 : fr = Flxn_red(f, n);
3774 203456 : if (mask>1 || !g)
3775 : {
3776 149381 : u = Flxn_mul_pre(W, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
3777 149707 : W = Flx_sub(W, Flx_shift(u, n2), p);
3778 : } else
3779 : {
3780 54075 : GEN y = Flxn_mul_pre(g, W, n, p, pi), yt = Flxn_red(y, n-n2);
3781 54076 : u = Flxn_mul_pre(yt, Flxn_mulhigh(fr, W, n2, n, p, pi), n-n2, p, pi);
3782 54084 : W = Flx_sub(y, Flx_shift(u, n2), p);
3783 : }
3784 203600 : if (gc_needed(av2,2))
3785 : {
3786 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Flxn_div, e = %ld", n);
3787 0 : W = gerepileupto(av2, W);
3788 : }
3789 : }
3790 55176 : return gerepileupto(av, W);
3791 : }
3792 : GEN
3793 55122 : Flxn_div(GEN g, GEN f, long e, ulong p)
3794 55122 : { return Flxn_div_pre(g, f, e, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
3795 :
3796 : GEN
3797 1032 : Flxn_inv(GEN f, long e, ulong p)
3798 1032 : { return Flxn_div(NULL, f, e, p); }
3799 :
3800 : GEN
3801 109372 : Flxn_expint(GEN h, long e, ulong p)
3802 : {
3803 109372 : pari_sp av = avma, av2;
3804 109372 : long v = h[1], n=1;
3805 109372 : GEN f = pol1_Flx(v), g = pol1_Flx(v);
3806 109336 : ulong mask = quadratic_prec_mask(e), pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
3807 109352 : av2 = avma;
3808 422594 : for (;mask>1;)
3809 : {
3810 : GEN u, w;
3811 422523 : long n2 = n;
3812 422523 : n<<=1; if (mask & 1) n--;
3813 422523 : mask >>= 1;
3814 422523 : u = Flxn_mul_pre(g, Flx_mulhigh_i(f, Flxn_red(h, n2-1), n2-1, p,pi), n-n2, p,pi);
3815 422503 : u = Flx_add(u, Flx_shift(Flxn_red(h, n-1), 1-n2), p);
3816 422546 : w = Flxn_mul_pre(f, Flx_integXn(u, n2-1, p), n-n2, p, pi);
3817 422513 : f = Flx_add(f, Flx_shift(w, n2), p);
3818 422635 : if (mask<=1) break;
3819 313278 : u = Flxn_mul_pre(g, Flxn_mulhigh(f, g, n2, n, p, pi), n-n2, p, pi);
3820 313244 : g = Flx_sub(g, Flx_shift(u, n2), p);
3821 313242 : if (gc_needed(av2,2))
3822 : {
3823 0 : if (DEBUGMEM>1) pari_warn(warnmem,"Flxn_exp, e = %ld", n);
3824 0 : gerepileall(av2, 2, &f, &g);
3825 : }
3826 : }
3827 109428 : return gerepileupto(av, f);
3828 : }
3829 :
3830 : GEN
3831 0 : Flxn_exp(GEN h, long e, ulong p)
3832 : {
3833 0 : if (degpol(h)<1 || uel(h,2)!=0)
3834 0 : pari_err_DOMAIN("Flxn_exp","valuation", "<", gen_1, h);
3835 0 : return Flxn_expint(Flx_deriv(h, p), e, p);
3836 : }
3837 :
3838 : INLINE GEN
3839 217181 : Flxn_recip(GEN x, long n)
3840 : {
3841 217181 : GEN z=Flx_recipspec(x+2,lgpol(x),n);
3842 217035 : z[1]=x[1];
3843 217035 : return z;
3844 : }
3845 :
3846 : GEN
3847 54067 : Flx_Newton(GEN P, long n, ulong p)
3848 : {
3849 54067 : pari_sp av = avma;
3850 54067 : long d = degpol(P);
3851 54064 : GEN dP = Flxn_recip(Flx_deriv(P, p), d);
3852 53974 : GEN Q = Flxn_div(dP, Flxn_recip(P, d+1), n, p);
3853 54041 : return gerepileuptoleaf(av, Q);
3854 : }
3855 :
3856 : GEN
3857 109368 : Flx_fromNewton(GEN P, ulong p)
3858 : {
3859 109368 : pari_sp av = avma;
3860 109368 : ulong n = Flx_constant(P)+1;
3861 109367 : GEN z = Flx_neg(Flx_shift(P, -1), p);
3862 109372 : GEN Q = Flxn_recip(Flxn_expint(z, n, p), n);
3863 109347 : return gerepileuptoleaf(av, Q);
3864 : }
3865 :
3866 : static long
3867 483 : newtonlogint(ulong n, ulong pp)
3868 : {
3869 483 : long s = 0;
3870 2107 : while (n > pp)
3871 : {
3872 1624 : s += ulogint(n, pp);
3873 1624 : n = (n+1)>>1;
3874 : }
3875 483 : return s;
3876 : }
3877 :
3878 : static void
3879 12500 : init_invlaplace(long d, ulong p, GEN *pt_P, GEN *pt_V)
3880 : {
3881 : long i;
3882 : ulong e;
3883 12500 : GEN P = cgetg(d+1, t_VECSMALL);
3884 12500 : GEN V = cgetg(d+1, t_VECSMALL);
3885 1389427 : for (i=1, e=1; i<=d; i++, e++)
3886 : {
3887 1376927 : if (e==p)
3888 : {
3889 455590 : e = 0;
3890 455590 : V[i] = u_lvalrem(i, p, &uel(P,i));
3891 : } else
3892 : {
3893 921337 : V[i] = 0; uel(P,i) = i;
3894 : }
3895 : }
3896 12500 : *pt_P = P; *pt_V = V;
3897 12500 : }
3898 :
3899 : /* return p^val * FpX_invLaplace(1+x+...x^(n-1), q), with q a power of p and
3900 : * val large enough to compensate for the power of p in the factorials */
3901 :
3902 : static GEN
3903 483 : ZpX_invLaplace_init(long n, GEN q, ulong p, long v, long sv)
3904 : {
3905 483 : pari_sp av = avma;
3906 483 : long i, d = n-1, w;
3907 : GEN y, W, E, t;
3908 483 : init_invlaplace(d, p, &E, &W);
3909 483 : t = Fp_inv(FpV_prod(Flv_to_ZV(E), q), q);
3910 483 : w = zv_sum(W);
3911 483 : if (v > w) t = Fp_mul(t, powuu(p, v-w), q);
3912 483 : y = cgetg(d+3,t_POL);
3913 483 : y[1] = evalsigne(1) | sv;
3914 21728 : for (i=d; i>=1; i--)
3915 : {
3916 21245 : gel(y,i+2) = t;
3917 21245 : t = Fp_mulu(t, uel(E,i), q);
3918 21245 : if (uel(W,i)) t = Fp_mul(t, powuu(p, uel(W,i)), q);
3919 : }
3920 483 : gel(y,2) = t;
3921 483 : return gerepilecopy(av, ZX_renormalize(y, d+3));
3922 : }
3923 :
3924 : static GEN
3925 27513 : Flx_composedsum(GEN P, GEN Q, ulong p)
3926 : {
3927 27513 : long n = 1 + degpol(P)*degpol(Q);
3928 27513 : ulong lead = Fl_mul(Fl_powu(Flx_lead(P), degpol(Q), p),
3929 27510 : Fl_powu(Flx_lead(Q), degpol(P), p), p);
3930 : GEN R;
3931 27515 : if (p >= (ulong)n)
3932 : {
3933 27032 : GEN Pl = Flx_invLaplace(Flx_Newton(P,n,p), p);
3934 27040 : GEN Ql = Flx_invLaplace(Flx_Newton(Q,n,p), p);
3935 27037 : GEN L = Flx_Laplace(Flxn_mul(Pl, Ql, n, p), p);
3936 27041 : R = Flx_fromNewton(L, p);
3937 : } else
3938 : {
3939 483 : long v = factorial_lval(n-1, p);
3940 483 : long w = 1 + newtonlogint(n-1, p);
3941 483 : GEN pv = powuu(p, v);
3942 483 : GEN qf = powuu(p, w), q = mulii(pv, qf), q2 = mulii(q, pv);
3943 483 : GEN iL = ZpX_invLaplace_init(n, q, p, v, P[1]);
3944 483 : GEN Pl = FpX_convol(iL, FpX_Newton(Flx_to_ZX(P), n, qf), q);
3945 483 : GEN Ql = FpX_convol(iL, FpX_Newton(Flx_to_ZX(Q), n, qf), q);
3946 483 : GEN Ln = ZX_Z_divexact(FpXn_mul(Pl, Ql, n, q2), pv);
3947 483 : GEN L = ZX_Z_divexact(FpX_Laplace(Ln, q), pv);
3948 483 : R = ZX_to_Flx(FpX_fromNewton(L, qf), p);
3949 : }
3950 27515 : return Flx_Fl_mul(R, lead, p);
3951 : }
3952 :
3953 : GEN
3954 27513 : Flx_direct_compositum(GEN a, GEN b, ulong p)
3955 27513 : { return Flx_composedsum(a, b, p); }
3956 :
3957 : static GEN
3958 3812 : _Flx_direct_compositum(void *E, GEN a, GEN b)
3959 3812 : { return Flx_direct_compositum(a, b, (ulong)E); }
3960 :
3961 : GEN
3962 28769 : FlxV_direct_compositum(GEN V, ulong p)
3963 28769 : { return gen_product(V, (void *)p, &_Flx_direct_compositum); }
3964 :
3965 : /* (x+1)^n mod p; assume 2 <= n < 2p prime */
3966 : static GEN
3967 0 : Fl_Xp1_powu(ulong n, ulong p, long v)
3968 : {
3969 0 : ulong k, d = (n + 1) >> 1;
3970 0 : GEN C, V = identity_zv(d);
3971 :
3972 0 : Flv_inv_inplace(V, p); /* could restrict to odd integers in [3,d] */
3973 0 : C = cgetg(n+3, t_VECSMALL);
3974 0 : C[1] = v;
3975 0 : uel(C,2) = 1UL;
3976 0 : uel(C,3) = n%p;
3977 0 : uel(C,4) = Fl_mul(odd(n)? n: n-1, n >> 1, p);
3978 : /* binom(n,k) = binom(n,k-1) * (n-k+1) / k */
3979 0 : if (SMALL_ULONG(p))
3980 0 : for (k = 3; k <= d; k++)
3981 0 : uel(C,k+2) = Fl_mul(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p);
3982 : else
3983 : {
3984 0 : ulong pi = get_Fl_red(p);
3985 0 : for (k = 3; k <= d; k++)
3986 0 : uel(C,k+2) = Fl_mul_pre(Fl_mul(n-k+1, uel(C,k+1), p), uel(V,k), p, pi);
3987 : }
3988 0 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
3989 0 : return C; /* normalized */
3990 : }
3991 :
3992 : /* p arbitrary */
3993 : GEN
3994 28194 : Flx_translate1_basecase(GEN P, ulong p)
3995 : {
3996 28194 : GEN R = Flx_copy(P);
3997 28194 : long i, k, n = degpol(P);
3998 654732 : for (i = 1; i <= n; i++)
3999 14846488 : for (k = n-i; k < n; k++) uel(R,k+2) = Fl_add(uel(R,k+2), uel(R,k+3), p);
4000 28194 : return R;
4001 : }
4002 :
4003 : static int
4004 41359 : translate_basecase(long n, ulong p)
4005 : {
4006 : #ifdef LONG_IS_64BIT
4007 36066 : if (p <= 19) return n < 40;
4008 29910 : if (p < 1UL<<30) return n < 58;
4009 0 : if (p < 1UL<<59) return n < 100;
4010 0 : if (p < 1UL<<62) return n < 120;
4011 0 : if (p < 1UL<<63) return n < 240;
4012 0 : return n < 250;
4013 : #else
4014 5293 : if (p <= 13) return n < 18;
4015 4136 : if (p <= 17) return n < 22;
4016 4078 : if (p <= 29) return n < 39;
4017 3886 : if (p <= 67) return n < 69;
4018 3667 : if (p < 1UL<< 15) return n < 80;
4019 2047 : if (p < 1UL<< 16) return n < 100;
4020 0 : if (p < 1UL<< 28) return n < 300;
4021 0 : return n < 650;
4022 : #endif
4023 : }
4024 : /* assume p prime */
4025 : GEN
4026 16100 : Flx_translate1(GEN P, ulong p)
4027 : {
4028 16100 : long d, n = degpol(P);
4029 : GEN R, Q, S;
4030 16100 : if (translate_basecase(n, p)) return Flx_translate1_basecase(P, p);
4031 : /* n > 0 */
4032 1148 : d = n >> 1;
4033 1148 : if ((ulong)n < p)
4034 : {
4035 0 : R = Flx_translate1(Flxn_red(P, d), p);
4036 0 : Q = Flx_translate1(Flx_shift(P, -d), p);
4037 0 : S = Fl_Xp1_powu(d, p, P[1]);
4038 0 : return Flx_add(Flx_mul(Q, S, p), R, p);
4039 : }
4040 : else
4041 : {
4042 : ulong q;
4043 1148 : if ((ulong)d > p) (void)ulogintall(d, p, &q); else q = p;
4044 1148 : R = Flx_translate1(Flxn_red(P, q), p);
4045 1148 : Q = Flx_translate1(Flx_shift(P, -q), p);
4046 1148 : S = Flx_add(Flx_shift(Q, q), Q, p);
4047 1148 : return Flx_add(S, R, p); /* P(x+1) = Q(x+1) (x^q+1) + R(x+1) */
4048 : }
4049 : }
4050 :
4051 : static GEN
4052 12017 : zl_Xp1_powu(ulong n, ulong p, ulong q, long e, long vs)
4053 : {
4054 12017 : ulong k, d = n >> 1, c, v = 0;
4055 12017 : GEN C, V, W, U = upowers(p, e-1);
4056 12017 : init_invlaplace(d, p, &V, &W);
4057 12017 : Flv_inv_inplace(V, q);
4058 12017 : C = cgetg(n+3, t_VECSMALL);
4059 12017 : C[1] = vs;
4060 12017 : uel(C,2) = 1UL;
4061 12017 : uel(C,3) = n%q;
4062 12017 : v = u_lvalrem(n, p, &c);
4063 1355682 : for (k = 2; k <= d; k++)
4064 : {
4065 : ulong w;
4066 1343665 : v += u_lvalrem(n-k+1, p, &w) - W[k];
4067 1343665 : c = Fl_mul(Fl_mul(w%q, c, q), uel(V,k), q);
4068 1343665 : uel(C,2+k) = v >= (ulong)e ? 0: v==0 ? c : Fl_mul(c, uel(U, v+1), q);
4069 : }
4070 1374521 : for ( ; k <= n; k++) uel(C,2+k) = uel(C,2+n-k);
4071 12017 : return C; /* normalized */
4072 : }
4073 :
4074 : GEN
4075 25259 : zlx_translate1(GEN P, ulong p, long e)
4076 : {
4077 25259 : ulong d, q = upowuu(p,e), n = degpol(P);
4078 : GEN R, Q, S;
4079 25259 : if (translate_basecase(n, q)) return Flx_translate1_basecase(P, q);
4080 : /* n > 0 */
4081 12017 : d = n >> 1;
4082 12017 : R = zlx_translate1(Flxn_red(P, d), p, e);
4083 12017 : Q = zlx_translate1(Flx_shift(P, -d), p, e);
4084 12017 : S = zl_Xp1_powu(d, p, q, e, P[1]);
4085 12017 : return Flx_add(Flx_mul(Q, S, q), R, q);
4086 : }
4087 :
4088 : /***********************************************************************/
4089 : /** Fl2 **/
4090 : /***********************************************************************/
4091 : /* Fl2 objects are Flv of length 2 [a,b] representing a+bsqrt(D) for
4092 : * a nonsquare D. */
4093 :
4094 : INLINE GEN
4095 6674290 : mkF2(ulong a, ulong b) { return mkvecsmall2(a,b); }
4096 :
4097 : /* allow pi = 0 */
4098 : GEN
4099 1792172 : Fl2_mul_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4100 : {
4101 : ulong xaya, xbyb, Db2, mid, z1, z2;
4102 1792172 : ulong x1 = x[1], x2 = x[2], y1 = y[1], y2 = y[2];
4103 1792172 : if (pi)
4104 : {
4105 1792533 : xaya = Fl_mul_pre(x1,y1,p,pi);
4106 1793041 : if (x2==0 && y2==0) return mkF2(xaya,0);
4107 1734203 : if (x2==0) return mkF2(xaya,Fl_mul_pre(x1,y2,p,pi));
4108 1712136 : if (y2==0) return mkF2(xaya,Fl_mul_pre(x2,y1,p,pi));
4109 1711923 : xbyb = Fl_mul_pre(x2,y2,p,pi);
4110 1711887 : mid = Fl_mul_pre(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p,pi);
4111 1712036 : Db2 = Fl_mul_pre(D, xbyb, p,pi);
4112 : }
4113 0 : else if (p & HIGHMASK)
4114 : {
4115 0 : xaya = Fl_mul(x1,y1,p);
4116 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4117 0 : if (x2==0) return mkF2(xaya,Fl_mul(x1,y2,p));
4118 0 : if (y2==0) return mkF2(xaya,Fl_mul(x2,y1,p));
4119 0 : xbyb = Fl_mul(x2,y2,p);
4120 0 : mid = Fl_mul(Fl_add(x1,x2,p), Fl_add(y1,y2,p),p);
4121 0 : Db2 = Fl_mul(D, xbyb, p);
4122 : }
4123 : else
4124 : {
4125 0 : xaya = (x1 * y1) % p;
4126 0 : if (x2==0 && y2==0) return mkF2(xaya,0);
4127 0 : if (x2==0) return mkF2(xaya, (x1 * y2) % p);
4128 0 : if (y2==0) return mkF2(xaya, (x2 * y1) % p);
4129 0 : xbyb = (x2 * y2) % p;
4130 0 : mid = (Fl_add(x1,x2,p) * Fl_add(y1,y2,p)) % p;
4131 0 : Db2 = (D * xbyb) % p;
4132 : }
4133 1711965 : z1 = Fl_add(xaya,Db2,p);
4134 1711860 : z2 = Fl_sub(mid,Fl_add(xaya,xbyb,p),p);
4135 1711809 : return mkF2(z1,z2);
4136 : }
4137 :
4138 : /* allow pi = 0 */
4139 : GEN
4140 4523983 : Fl2_sqr_pre(GEN x, ulong D, ulong p, ulong pi)
4141 : {
4142 4523983 : ulong a = x[1], b = x[2];
4143 : ulong a2, Db2, ab;
4144 4523983 : if (pi)
4145 : {
4146 4524699 : a2 = Fl_sqr_pre(a,p,pi);
4147 4528162 : if (b==0) return mkF2(a2,0);
4148 4342298 : Db2= Fl_mul_pre(D, Fl_sqr_pre(b,p,pi), p,pi);
4149 4342397 : ab = Fl_mul_pre(a,b,p,pi);
4150 : }
4151 0 : else if (p & HIGHMASK)
4152 : {
4153 0 : a2 = Fl_sqr(a,p);
4154 0 : if (b==0) return mkF2(a2,0);
4155 0 : Db2= Fl_mul(D, Fl_sqr(b,p), p);
4156 0 : ab = Fl_mul(a,b,p);
4157 : }
4158 : else
4159 : {
4160 0 : a2 = (a * a) % p;
4161 0 : if (b==0) return mkF2(a2,0);
4162 0 : Db2= (D * ((b * b) % p)) % p;
4163 0 : ab = (a * b) % p;
4164 : }
4165 4341570 : return mkF2(Fl_add(a2,Db2,p), Fl_double(ab,p));
4166 : }
4167 :
4168 : /* allow pi = 0 */
4169 : ulong
4170 67340 : Fl2_norm_pre(GEN x, ulong D, ulong p, ulong pi)
4171 : {
4172 67340 : ulong a = x[1], b = x[2], a2;
4173 67340 : if (pi)
4174 : {
4175 67340 : a2 = Fl_sqr_pre(a,p,pi);
4176 67340 : return b? Fl_sub(a2, Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p): a2;
4177 : }
4178 0 : else if (p & HIGHMASK)
4179 : {
4180 0 : a2 = Fl_sqr(a,p);
4181 0 : return b? Fl_sub(a2, Fl_mul(D, Fl_sqr(b, p), p), p): a2;
4182 : }
4183 : else
4184 : {
4185 0 : a2 = (a * a) % p;
4186 0 : return b? Fl_sub(a2, (D * ((b * b) % p)) % p, p): a2;
4187 : }
4188 : }
4189 :
4190 : /* allow pi = 0 */
4191 : GEN
4192 177208 : Fl2_inv_pre(GEN x, ulong D, ulong p, ulong pi)
4193 : {
4194 177208 : ulong a = x[1], b = x[2], n, ni;
4195 177208 : if (b == 0) return mkF2(Fl_inv(a,p), 0);
4196 151107 : b = Fl_neg(b, p);
4197 151107 : if (pi)
4198 : {
4199 151107 : n = Fl_sub(Fl_sqr_pre(a, p,pi),
4200 : Fl_mul_pre(D, Fl_sqr_pre(b, p,pi), p,pi), p);
4201 151107 : ni = Fl_inv(n,p);
4202 151112 : return mkF2(Fl_mul_pre(a, ni, p,pi), Fl_mul_pre(b, ni, p,pi));
4203 : }
4204 0 : else if (p & HIGHMASK)
4205 : {
4206 0 : n = Fl_sub(Fl_sqr(a, p), Fl_mul(D, Fl_sqr(b, p), p), p);
4207 0 : ni = Fl_inv(n,p);
4208 0 : return mkF2(Fl_mul(a, ni, p), Fl_mul(b, ni, p));
4209 : }
4210 : else
4211 : {
4212 0 : n = Fl_sub((a * a) % p, (D * ((b * b) % p)) % p, p);
4213 0 : ni = Fl_inv(n,p);
4214 0 : return mkF2((a * ni) % p, (b * ni) % p);
4215 : }
4216 : }
4217 :
4218 : int
4219 404404 : Fl2_equal1(GEN x) { return x[1]==1 && x[2]==0; }
4220 :
4221 : struct _Fl2 {
4222 : ulong p, pi, D;
4223 : };
4224 :
4225 : static GEN
4226 4523019 : _Fl2_sqr(void *data, GEN x)
4227 : {
4228 4523019 : struct _Fl2 *D = (struct _Fl2*)data;
4229 4523019 : return Fl2_sqr_pre(x, D->D, D->p, D->pi);
4230 : }
4231 : static GEN
4232 1763828 : _Fl2_mul(void *data, GEN x, GEN y)
4233 : {
4234 1763828 : struct _Fl2 *D = (struct _Fl2*)data;
4235 1763828 : return Fl2_mul_pre(x,y, D->D, D->p, D->pi);
4236 : }
4237 :
4238 : /* n-Power of x in Z/pZ[X]/(T), as t_VECSMALL; allow pi = 0 */
4239 : GEN
4240 604436 : Fl2_pow_pre(GEN x, GEN n, ulong D, ulong p, ulong pi)
4241 : {
4242 604436 : pari_sp av = avma;
4243 : struct _Fl2 d;
4244 : GEN y;
4245 604436 : long s = signe(n);
4246 604436 : if (!s) return mkF2(1,0);
4247 535611 : if (s < 0)
4248 177208 : x = Fl2_inv_pre(x,D,p,pi);
4249 535602 : if (is_pm1(n)) return s < 0 ? x : zv_copy(x);
4250 394590 : d.p = p; d.pi = pi; d.D=D;
4251 394590 : y = gen_pow_i(x, n, (void*)&d, &_Fl2_sqr, &_Fl2_mul);
4252 394585 : return gerepileuptoleaf(av, y);
4253 : }
4254 :
4255 : static GEN
4256 604406 : _Fl2_pow(void *data, GEN x, GEN n)
4257 : {
4258 604406 : struct _Fl2 *D = (struct _Fl2*)data;
4259 604406 : return Fl2_pow_pre(x, n, D->D, D->p, D->pi);
4260 : }
4261 :
4262 : static GEN
4263 102658 : _Fl2_rand(void *data)
4264 : {
4265 102658 : struct _Fl2 *D = (struct _Fl2*)data;
4266 102658 : ulong a = random_Fl(D->p), b=random_Fl(D->p-1)+1;
4267 102659 : return mkF2(a,b);
4268 : }
4269 :
4270 : static const struct bb_group Fl2_star={_Fl2_mul, _Fl2_pow, _Fl2_rand,
4271 : hash_GEN, zv_equal, Fl2_equal1, NULL};
4272 :
4273 : /* allow pi = 0 */
4274 : GEN
4275 68824 : Fl2_sqrtn_pre(GEN a, GEN n, ulong D, ulong p, ulong pi, GEN *zeta)
4276 : {
4277 : struct _Fl2 E;
4278 : GEN o;
4279 68824 : if (a[1]==0 && a[2]==0)
4280 : {
4281 0 : if (signe(n) < 0) pari_err_INV("Flxq_sqrtn",a);
4282 0 : if (zeta) *zeta=mkF2(1,0);
4283 0 : return zv_copy(a);
4284 : }
4285 68824 : E.p=p; E.pi = pi; E.D = D;
4286 68824 : o = subiu(powuu(p,2), 1);
4287 68825 : return gen_Shanks_sqrtn(a,n,o,zeta,(void*)&E,&Fl2_star);
4288 : }
4289 :
4290 : /* allow pi = 0 */
4291 : GEN
4292 10402 : Flx_Fl2_eval_pre(GEN x, GEN y, ulong D, ulong p, ulong pi)
4293 : {
4294 : GEN p1;
4295 10402 : long i = lg(x)-1;
4296 10402 : if (i <= 2)
4297 2065 : return mkF2(i == 2? x[2]: 0, 0);
4298 8337 : p1 = mkF2(x[i], 0);
4299 36456 : for (i--; i>=2; i--)
4300 : {
4301 28119 : p1 = Fl2_mul_pre(p1, y, D, p, pi);
4302 28119 : uel(p1,1) = Fl_add(uel(p1,1), uel(x,i), p);
4303 : }
4304 8337 : return p1;
4305 : }
4306 :
4307 : /***********************************************************************/
4308 : /** FlxV **/
4309 : /***********************************************************************/
4310 : /* FlxV are t_VEC with Flx coefficients. */
4311 :
4312 : GEN
4313 34482 : FlxV_Flc_mul(GEN V, GEN W, ulong p)
4314 : {
4315 34482 : pari_sp ltop=avma;
4316 : long i;
4317 34482 : GEN z = Flx_Fl_mul(gel(V,1),W[1],p);
4318 257068 : for(i=2;i<lg(V);i++)
4319 222586 : z=Flx_add(z,Flx_Fl_mul(gel(V,i),W[i],p),p);
4320 34482 : return gerepileuptoleaf(ltop,z);
4321 : }
4322 :
4323 : GEN
4324 0 : ZXV_to_FlxV(GEN x, ulong p)
4325 0 : { pari_APPLY_type(t_VEC, ZX_to_Flx(gel(x,i), p)) }
4326 :
4327 : GEN
4328 3768422 : ZXT_to_FlxT(GEN x, ulong p)
4329 : {
4330 3768422 : if (typ(x) == t_POL)
4331 3708428 : return ZX_to_Flx(x, p);
4332 : else
4333 196795 : pari_APPLY_type(t_VEC, ZXT_to_FlxT(gel(x,i), p))
4334 : }
4335 :
4336 : GEN
4337 169704 : FlxV_to_Flm(GEN x, long n)
4338 919882 : { pari_APPLY_type(t_MAT, Flx_to_Flv(gel(x,i), n)) }
4339 :
4340 : GEN
4341 0 : FlxV_red(GEN x, ulong p)
4342 0 : { pari_APPLY_type(t_VEC, Flx_red(gel(x,i), p)) }
4343 :
4344 : GEN
4345 309209 : FlxT_red(GEN x, ulong p)
4346 : {
4347 309209 : if (typ(x) == t_VECSMALL)
4348 208048 : return Flx_red(x, p);
4349 : else
4350 339427 : pari_APPLY_type(t_VEC, FlxT_red(gel(x,i), p))
4351 : }
4352 :
4353 : GEN
4354 113589 : FlxqV_dotproduct_pre(GEN x, GEN y, GEN T, ulong p, ulong pi)
4355 : {
4356 113589 : long i, lx = lg(x);
4357 : pari_sp av;
4358 : GEN c;
4359 113589 : if (lx == 1) return pol0_Flx(get_Flx_var(T));
4360 113589 : av = avma; c = Flx_mul_pre(gel(x,1),gel(y,1), p, pi);
4361 464499 : for (i=2; i<lx; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4362 113589 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4363 : }
4364 : GEN
4365 0 : FlxqV_dotproduct(GEN x, GEN y, GEN T, ulong p)
4366 0 : { return FlxqV_dotproduct_pre(x, y, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
4367 :
4368 : GEN
4369 2108 : FlxqX_dotproduct(GEN x, GEN y, GEN T, ulong p)
4370 : {
4371 2108 : long i, l = minss(lg(x), lg(y));
4372 : ulong pi;
4373 : pari_sp av;
4374 : GEN c;
4375 2108 : if (l == 2) return pol0_Flx(get_Flx_var(T));
4376 2108 : av = avma; pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4377 2108 : c = Flx_mul_pre(gel(x,2),gel(y,2), p, pi);
4378 6780 : for (i=3; i<l; i++) c = Flx_add(c, Flx_mul_pre(gel(x,i),gel(y,i), p, pi), p);
4379 2108 : return gerepileuptoleaf(av, Flx_rem_pre(c,T,p,pi));
4380 : }
4381 :
4382 : /* allow pi = 0 */
4383 : GEN
4384 250755 : FlxC_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4385 : {
4386 250755 : long i, l = lg(z);
4387 250755 : GEN y = cgetg(l, t_VECSMALL);
4388 12724217 : for (i=1; i<l; i++) uel(y,i) = Flx_eval_powers_pre(gel(z,i), x, p, pi);
4389 251009 : return y;
4390 : }
4391 :
4392 : /***********************************************************************/
4393 : /** FlxM **/
4394 : /***********************************************************************/
4395 : /* allow pi = 0 */
4396 : GEN
4397 19404 : FlxM_eval_powers_pre(GEN z, GEN x, ulong p, ulong pi)
4398 : {
4399 19404 : long i, l = lg(z);
4400 19404 : GEN y = cgetg(l, t_MAT);
4401 270162 : for (i=1; i<l; i++) gel(y,i) = FlxC_eval_powers_pre(gel(z,i), x, p, pi);
4402 19404 : return y;
4403 : }
4404 :
4405 : GEN
4406 0 : zero_FlxC(long n, long sv)
4407 : {
4408 0 : GEN x = cgetg(n + 1, t_COL), z = zero_Flx(sv);
4409 : long i;
4410 0 : for (i = 1; i <= n; i++) gel(x, i) = z;
4411 0 : return x;
4412 : }
4413 :
4414 : GEN
4415 0 : FlxC_neg(GEN x, ulong p)
4416 0 : { pari_APPLY_type(t_COL, Flx_neg(gel(x, i), p)) }
4417 :
4418 : GEN
4419 0 : FlxC_sub(GEN x, GEN y, ulong p)
4420 0 : { pari_APPLY_type(t_COL, Flx_sub(gel(x, i), gel(y, i), p)) }
4421 :
4422 : GEN
4423 0 : zero_FlxM(long r, long c, long sv)
4424 : {
4425 0 : GEN x = cgetg(c + 1, t_MAT), z = zero_FlxC(r, sv);
4426 : long j;
4427 0 : for (j = 1; j <= c; j++) gel(x, j) = z;
4428 0 : return x;
4429 : }
4430 :
4431 : GEN
4432 0 : FlxM_neg(GEN x, ulong p)
4433 0 : { pari_APPLY_same(FlxC_neg(gel(x, i), p)) }
4434 :
4435 : GEN
4436 0 : FlxM_sub(GEN x, GEN y, ulong p)
4437 0 : { pari_APPLY_same(FlxC_sub(gel(x, i), gel(y,i), p)) }
4438 :
4439 : GEN
4440 0 : FlxqC_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4441 0 : { pari_APPLY_type(t_COL, Flxq_mul(gel(x, i), y, T, p)) }
4442 :
4443 : GEN
4444 0 : FlxqM_Flxq_mul(GEN x, GEN y, GEN T, ulong p)
4445 0 : { pari_APPLY_same(FlxqC_Flxq_mul(gel(x, i), y, T, p)) }
4446 :
4447 : static GEN
4448 54700 : FlxM_pack_ZM(GEN M, GEN (*pack)(GEN, long)) {
4449 : long i, j, l, lc;
4450 54700 : GEN N = cgetg_copy(M, &l), x;
4451 54700 : if (l == 1)
4452 0 : return N;
4453 54700 : lc = lgcols(M);
4454 269372 : for (j = 1; j < l; j++) {
4455 214672 : gel(N, j) = cgetg(lc, t_COL);
4456 1107480 : for (i = 1; i < lc; i++) {
4457 892808 : x = gcoeff(M, i, j);
4458 892808 : gcoeff(N, i, j) = pack(x + 2, lgpol(x));
4459 : }
4460 : }
4461 54700 : return N;
4462 : }
4463 :
4464 : static GEN
4465 839129 : kron_pack_Flx_spec_half(GEN x, long l) {
4466 839129 : if (l == 0) return gen_0;
4467 480287 : return Flx_to_int_halfspec(x, l);
4468 : }
4469 :
4470 : static GEN
4471 50290 : kron_pack_Flx_spec(GEN x, long l) {
4472 : long i;
4473 : GEN w, y;
4474 50290 : if (l == 0)
4475 9482 : return gen_0;
4476 40808 : y = cgetipos(l + 2);
4477 150922 : for (i = 0, w = int_LSW(y); i < l; i++, w = int_nextW(w))
4478 110114 : *w = x[i];
4479 40808 : return y;
4480 : }
4481 :
4482 : static GEN
4483 3389 : kron_pack_Flx_spec_2(GEN x, long l) { return Flx_eval2BILspec(x, 2, l); }
4484 :
4485 : static GEN
4486 0 : kron_pack_Flx_spec_3(GEN x, long l) { return Flx_eval2BILspec(x, 3, l); }
4487 :
4488 : static GEN
4489 42080 : kron_unpack_Flx(GEN z, ulong p)
4490 : {
4491 42080 : long i, l = lgefint(z);
4492 42080 : GEN x = cgetg(l, t_VECSMALL), w;
4493 198461 : for (w = int_LSW(z), i = 2; i < l; w = int_nextW(w), i++)
4494 156381 : x[i] = ((ulong) *w) % p;
4495 42080 : return Flx_renormalize(x, l);
4496 : }
4497 :
4498 : static GEN
4499 2930 : kron_unpack_Flx_2(GEN x, ulong p) {
4500 2930 : long d = (lgefint(x)-1)/2 - 1;
4501 2930 : return Z_mod2BIL_Flx_2(x, d, p);
4502 : }
4503 :
4504 : static GEN
4505 0 : kron_unpack_Flx_3(GEN x, ulong p) {
4506 0 : long d = lgefint(x)/3 - 1;
4507 0 : return Z_mod2BIL_Flx_3(x, d, p);
4508 : }
4509 :
4510 : static GEN
4511 118539 : FlxM_pack_ZM_bits(GEN M, long b)
4512 : {
4513 : long i, j, l, lc;
4514 118539 : GEN N = cgetg_copy(M, &l), x;
4515 118539 : if (l == 1)
4516 0 : return N;
4517 118539 : lc = lgcols(M);
4518 492013 : for (j = 1; j < l; j++) {
4519 373474 : gel(N, j) = cgetg(lc, t_COL);
4520 5974153 : for (i = 1; i < lc; i++) {
4521 5600679 : x = gcoeff(M, i, j);
4522 5600679 : gcoeff(N, i, j) = kron_pack_Flx_spec_bits(x + 2, b, lgpol(x));
4523 : }
4524 : }
4525 118539 : return N;
4526 : }
4527 :
4528 : static GEN
4529 27350 : ZM_unpack_FlxqM(GEN M, GEN T, ulong p, ulong pi, GEN (*unpack)(GEN, ulong))
4530 : {
4531 27350 : long i, j, l, lc, sv = get_Flx_var(T);
4532 27350 : GEN N = cgetg_copy(M, &l), x;
4533 27350 : if (l == 1)
4534 0 : return N;
4535 27350 : lc = lgcols(M);
4536 161343 : for (j = 1; j < l; j++) {
4537 133993 : gel(N, j) = cgetg(lc, t_COL);
4538 762713 : for (i = 1; i < lc; i++) {
4539 628720 : x = unpack(gcoeff(M, i, j), p);
4540 628720 : x[1] = sv;
4541 628720 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4542 : }
4543 : }
4544 27350 : return N;
4545 : }
4546 :
4547 : static GEN
4548 59310 : ZM_unpack_FlxqM_bits(GEN M, long b, GEN T, ulong p, ulong pi)
4549 : {
4550 59310 : long i, j, l, lc, sv = get_Flx_var(T);
4551 59310 : GEN N = cgetg_copy(M, &l), x;
4552 59310 : if (l == 1)
4553 0 : return N;
4554 59310 : lc = lgcols(M);
4555 59310 : if (b < BITS_IN_LONG) {
4556 202047 : for (j = 1; j < l; j++) {
4557 144420 : gel(N, j) = cgetg(lc, t_COL);
4558 3254377 : for (i = 1; i < lc; i++) {
4559 3109957 : x = kron_unpack_Flx_bits_narrow(gcoeff(M, i, j), b, p);
4560 3109957 : x[1] = sv;
4561 3109957 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4562 : }
4563 : }
4564 : } else {
4565 1683 : ulong pi = get_Fl_red(p);
4566 9844 : for (j = 1; j < l; j++) {
4567 8161 : gel(N, j) = cgetg(lc, t_COL);
4568 175361 : for (i = 1; i < lc; i++) {
4569 167200 : x = kron_unpack_Flx_bits_wide(gcoeff(M, i, j), b, p, pi);
4570 167200 : x[1] = sv;
4571 167200 : gcoeff(N, i, j) = Flx_rem_pre(x, T, p, pi);
4572 : }
4573 : }
4574 : }
4575 59310 : return N;
4576 : }
4577 :
4578 : GEN
4579 86660 : FlxqM_mul_Kronecker(GEN A, GEN B, GEN T, ulong p)
4580 : {
4581 86660 : pari_sp av = avma;
4582 86660 : long b, d = get_Flx_degree(T), n = lg(A) - 1;
4583 : GEN C, D, z;
4584 : GEN (*pack)(GEN, long), (*unpack)(GEN, ulong);
4585 86660 : ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
4586 86660 : int is_sqr = A==B;
4587 :
4588 86660 : z = muliu(muliu(sqru(p - 1), d), n);
4589 86660 : b = expi(z) + 1;
4590 : /* only do expensive bit-packing if it saves at least 1 limb */
4591 86660 : if (b <= BITS_IN_HALFULONG)
4592 82485 : { if (nbits2nlong(d*b) == (d + 1)/2) b = BITS_IN_HALFULONG; }
4593 : else
4594 : {
4595 4175 : long l = lgefint(z) - 2;
4596 4175 : if (nbits2nlong(d*b) == d*l) b = l*BITS_IN_LONG;
4597 : }
4598 86660 : set_avma(av);
4599 :
4600 86660 : switch (b) {
4601 26581 : case BITS_IN_HALFULONG:
4602 26581 : pack = kron_pack_Flx_spec_half;
4603 26581 : unpack = int_to_Flx_half;
4604 26581 : break;
4605 720 : case BITS_IN_LONG:
4606 720 : pack = kron_pack_Flx_spec;
4607 720 : unpack = kron_unpack_Flx;
4608 720 : break;
4609 49 : case 2*BITS_IN_LONG:
4610 49 : pack = kron_pack_Flx_spec_2;
4611 49 : unpack = kron_unpack_Flx_2;
4612 49 : break;
4613 0 : case 3*BITS_IN_LONG:
4614 0 : pack = kron_pack_Flx_spec_3;
4615 0 : unpack = kron_unpack_Flx_3;
4616 0 : break;
4617 59310 : default:
4618 59310 : A = FlxM_pack_ZM_bits(A, b);
4619 59310 : B = is_sqr? A: FlxM_pack_ZM_bits(B, b);
4620 59310 : C = ZM_mul(A, B);
4621 59310 : D = ZM_unpack_FlxqM_bits(C, b, T, p, pi);
4622 59310 : return gerepilecopy(av, D);
4623 : }
4624 27350 : A = FlxM_pack_ZM(A, pack);
4625 27350 : B = is_sqr? A: FlxM_pack_ZM(B, pack);
4626 27350 : C = ZM_mul(A, B);
4627 27350 : D = ZM_unpack_FlxqM(C, T, p, pi, unpack);
4628 27350 : return gerepilecopy(av, D);
4629 : }
|
