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