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