Line data Source code
1 : /* Copyright (C) 2000 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 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define DEBUGLEVEL DEBUGLEVEL_hensel
18 :
19 : /* Assume n > 0. We want to go to accuracy n, starting from accuracy 1, using
20 : * a quadratically convergent algorithm. Goal: 9 -> 1,2,3,5,9 instead of
21 : * 1,2,4,8,9 (sequence of accuracies).
22 : *
23 : * Let a0 = 1, a1 = 2, a2, ... ak = n, the sequence of accuracies. To obtain
24 : * it, work backwards:
25 : * a(k) = n, a(i-1) = (a(i) + 1) \ 2,
26 : * but we do not want to store a(i) explicitly, even as a t_VECSMALL, since
27 : * this would leave an object on the stack. We store a(i) implicitly in a
28 : * MASK: let a(0) = 1, if the i-bit of MASK is set, set a(i+1) = 2 a(i) - 1,
29 : * and 2a(i) otherwise.
30 : *
31 : * In fact, we do something a little more complicated to simplify the
32 : * function interface and avoid returning k and MASK separately: we return
33 : * MASK + 2^(k+1), so the highest bit of the mask indicates the length of the
34 : * sequence, and the following ones are as above. */
35 :
36 : ulong
37 6649354 : quadratic_prec_mask(long n)
38 : {
39 6649354 : long a = n, i;
40 6649354 : ulong mask = 0;
41 6649354 : for(i = 1;; i++, mask <<= 1)
42 : {
43 11195620 : mask |= (a&1); a = (a+1)>>1;
44 17844974 : if (a==1) return mask | (1UL << i);
45 : }
46 : }
47 :
48 : /***********************************************************************/
49 : /** **/
50 : /** Zp **/
51 : /** **/
52 : /***********************************************************************/
53 :
54 : static GEN
55 1963003 : Zp_divlift(GEN b, GEN a, GEN x, GEN p, long n)
56 : {
57 1963003 : pari_sp ltop = avma, av;
58 : ulong mask;
59 1963003 : GEN q = p;
60 1963003 : if (n == 1) return gcopy(x);
61 1962702 : mask = quadratic_prec_mask(n);
62 1962703 : av = avma;
63 12037472 : while (mask > 1)
64 : {
65 10074770 : GEN v, q2 = q;
66 10074770 : q = sqri(q);
67 10074765 : if (mask & 1UL) q = diviiexact(q,p);
68 10074766 : mask >>= 1;
69 10074766 : if (mask > 1 || !b)
70 : {
71 9967538 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
72 9967540 : x = Fp_sub(x, Fp_mul(v, x, q), q);
73 : }
74 : else
75 : {
76 107228 : GEN y = Fp_mul(x, b, q), yt = modii(y, q2);
77 107228 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
78 107229 : x = Fp_sub(y, Fp_mul(v, yt, q), q);
79 : }
80 10074771 : if (gc_needed(av, 1))
81 : {
82 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_Zp_Newton");
83 0 : gerepileall(av, 2, &x, &q);
84 : }
85 : }
86 1962702 : return gerepileupto(ltop, x);
87 : }
88 :
89 : GEN
90 1855774 : Zp_invlift(GEN a, GEN x, GEN p, long e)
91 1855774 : { return Zp_divlift(NULL, a, x, p, e); }
92 :
93 : GEN
94 1855774 : Zp_inv(GEN a, GEN p, long e)
95 : {
96 1855774 : pari_sp av=avma;
97 : GEN ai;
98 1855774 : if (lgefint(p)==3)
99 : {
100 1855511 : ulong pp = p[2];
101 1855511 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
102 : } else
103 263 : ai = Fp_inv(modii(a, p), p);
104 1855774 : return gerepileupto(av, Zp_invlift(a, ai, p, e));
105 : }
106 :
107 : GEN
108 107228 : Zp_div(GEN b, GEN a, GEN p, long e)
109 : {
110 107228 : pari_sp av=avma;
111 : GEN ai;
112 107228 : if (lgefint(p)==3)
113 : {
114 107165 : ulong pp = p[2];
115 107165 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
116 : } else
117 63 : ai = Fp_inv(modii(a, p), p);
118 107229 : return gerepileupto(av, Zp_divlift(b, a, ai, p, e));
119 : }
120 :
121 : static GEN
122 616 : mul2n(void *E, GEN x, GEN y) { return remi2n(mulii(x, y), (long)E); }
123 : static GEN
124 749 : sqr2n(void *E, GEN x) { return remi2n(sqri(x), (long)E); }
125 :
126 : /* a^n mod 2^e using remi2n (result has the same sign as a) */
127 : static GEN
128 476 : Fp_pow2n(GEN a, GEN n, long e)
129 476 : { return gen_pow(a, n, (void*)e, &sqr2n, &mul2n); }
130 :
131 : /* Same as ZpX_liftroot for the polynomial X^n-b*/
132 : GEN
133 173915 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
134 : {
135 173915 : pari_sp av = avma;
136 : int nis2, pis2;
137 : GEN q, w, n_1;
138 : ulong mask;
139 :
140 173915 : if (e == 1) return icopy(a);
141 168944 : nis2 = equaliu(n, 2); n_1 = nis2? NULL: subiu(n,1);
142 168944 : pis2 = equaliu(p, 2);
143 168944 : mask = quadratic_prec_mask(e);
144 168944 : w = nis2 ? shifti(a,1): Fp_mul(n, Fp_pow(a,n_1,p), p);
145 168944 : w = Fp_inv(w, p);
146 168944 : q = p; /* q = p^e; use e instead of q iff p = 2 */
147 168944 : e = 1;
148 : for(;;)
149 : {
150 266496 : if (pis2)
151 : {
152 343 : e <<= 1; if (mask & 1) e--;
153 343 : mask >>= 1;
154 : /* a -= w (a^n - b) */
155 343 : a = remi2n(subii(a, mulii(w, subii(Fp_pow2n(a, n, e), b))), e);
156 343 : if (mask == 1) break;
157 : /* w += w - w^2 n a^(n-1)*/
158 133 : w = subii(shifti(w,1), remi2n(mulii(remi2n(sqri(w), e),
159 : mulii(n, Fp_pow2n(a, n_1, e))), e));
160 133 : continue;
161 : }
162 266153 : q = sqri(q); if (mask & 1) q = diviiexact(q, p);
163 266153 : mask >>= 1;
164 266153 : if (lgefint(q) == 3 && lgefint(n) == 3)
165 97281 : {
166 265732 : ulong Q = uel(q,2), N = uel(n,2);
167 265732 : ulong A = umodiu(a, Q);
168 265732 : ulong B = umodiu(b, Q);
169 265732 : ulong W = umodiu(w, Q);
170 :
171 265732 : A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
172 265732 : a = utoi(A);
173 265732 : if (mask == 1) break;
174 97281 : if (nis2)
175 84293 : W = Fl_double(Fl_sub(W, Fl_mul(Fl_sqr(W,Q), A, Q), Q), Q);
176 : else
177 12988 : W = Fl_sub(Fl_double(W,Q),
178 : Fl_mul(Fl_sqr(W,Q), Fl_mul(N,Fl_powu(A, N-1, Q), Q), Q), Q);
179 97281 : w = utoi(W);
180 : }
181 : else
182 : {
183 : /* a -= w (a^n - b) */
184 421 : a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q), b))), q);
185 421 : if (mask == 1) break;
186 : /* w += w - w^2 n a^(n-1)*/
187 138 : if (nis2)
188 15 : w = shifti(subii(w, Fp_mul(Fp_sqr(w,q), a, q)), 1);
189 : else
190 123 : w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
191 : mulii(n, Fp_pow(a,n_1,q)), q));
192 : }
193 : }
194 168944 : if (pis2 && signe(a) < 0) a = addii(a, int2n(e));
195 168944 : return gerepileuptoint(av, a);
196 : }
197 :
198 : /* ZpX_liftroot for the polynomial X^2-b */
199 : GEN
200 139972 : Zp_sqrtlift(GEN b, GEN a, GEN p, long e)
201 139972 : { return Zp_sqrtnlift(b, gen_2, a, p, e); }
202 :
203 : GEN
204 1379248 : Zp_sqrt(GEN x, GEN p, long e)
205 : {
206 : pari_sp av;
207 : GEN z;
208 1379248 : if (absequaliu(p,2)) return Z2_sqrt(x,e);
209 1377702 : av = avma; z = Fp_sqrt(Fp_red(x, p), p);
210 1377714 : if (!z) return NULL;
211 1377651 : if (e > 1) z = Zp_sqrtlift(x, z, p, e);
212 1377651 : return gerepileuptoint(av, z);
213 : }
214 :
215 : /* p-adic logarithm of a = 1 mod p (adapted from a C program of Xavier Caruso)
216 : * Algorithm:
217 : * 1. raise a at the power p^(v-1) in order to make it closer to 1
218 : * 2. write the new a as a product
219 : * 1 / a = (1 - a_0*p^v) (1 - a_1*p^(2*v) (1 - a_2*p^(4*v) ...
220 : * with 0 <= a_i < p^(v*2^i).
221 : * 3. compute each log(1 - a_i*p^(v*2^i)) using Taylor expansion
222 : * and binary splitting */
223 : GEN
224 317707 : Zp_log(GEN a, GEN p, ulong e)
225 : {
226 317707 : pari_sp av = avma;
227 317707 : ulong v, N, Np, trunc, pp = itou_or_0(p);
228 317707 : GEN pe, pv, trunc_mod, num, den, ans = gen_0;
229 :
230 317707 : if (equali1(a)) return ans; /* very frequent! */
231 : /* First make the argument closer to 1 by raising it to the p^(v-1) */
232 316685 : v = pp? ulogint(e, pp): 0; /* v here is v-1 */
233 316685 : pe = powiu(p,e); pv = powiu(p,v);
234 316686 : a = Fp_pow(a, pv, mulii(pe, pv));
235 316686 : e += v;
236 :
237 : /* Where do we truncate the Taylor expansion */
238 316686 : N = e + v; N /= ++v; /* note the ++v */
239 316686 : Np = N;
240 : while(1)
241 0 : {
242 316686 : ulong e = Np;
243 316686 : if (pp) e += ulogint(N, pp) / v;
244 316686 : if (e == N) break;
245 0 : N = e;
246 : }
247 :
248 316686 : num = cgetg(N+1, t_INT);
249 316686 : den = cgetg(N+1, t_INT);
250 316686 : trunc = v << 1;
251 316686 : trunc_mod = powiu(p, trunc);
252 : while(1)
253 1190958 : { /* compute f = 1 - a_i*p^((v+1)*2^i); trunc_mod = p^((v+1)*2^(i+1)) */
254 1507644 : GEN f = modii(a, trunc_mod);
255 1507643 : if (!equali1(f))
256 : {
257 1502855 : ulong i, step = 1;
258 : GEN h, hpow;
259 :
260 1502855 : f = subui(2, f);
261 1502855 : a = mulii(a, f);
262 :
263 : /* compute the Taylor expansion of log(f), over Q for now */
264 12738482 : for (i = 1; i <= N; i++) { gel(num,i) = gen_1; gel(den,i) = utoipos(i); }
265 1502856 : hpow = h = subui(1, f); /* h = a_i*p^(2^i) */
266 : while(1)
267 : {
268 13822488 : for (i = 1; i <= N - step; i += step << 1)
269 : {
270 9732770 : GEN t = mulii(mulii(hpow, gel(num,i+step)), gel(den,i));
271 9732770 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
272 9732769 : gel(num,i) = addii(gel(num,i), t);
273 9732770 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
274 : }
275 4089718 : step <<= 1; if (step >= N) break;
276 2586863 : hpow = sqri(hpow);
277 : }
278 :
279 1502855 : if (pp)
280 : { /* simplify the fraction */
281 1502599 : GEN d = powuu(pp, factorial_lval(N, pp));
282 1502599 : gel(num,1) = diviiexact(gel(num,1), d);
283 1502599 : gel(den,1) = diviiexact(gel(den,1), d);
284 : }
285 :
286 1502855 : h = diviiexact(h, pv);
287 1502855 : ans = addii(ans, mulii(mulii(gel(num,1), h), Zp_inv(gel(den,1), p, e)));
288 : }
289 1507644 : if (trunc > e) break;
290 1190958 : trunc_mod = sqri(trunc_mod); trunc <<= 1; N >>= 1;
291 : }
292 316686 : return gerepileuptoint(av, modii(ans, pe));
293 : }
294 :
295 : /* p-adic exponential of a = 0 (mod 2p)
296 : * 1. write a as a sum a = a_0*p + a_1*p^2 + a_2*p^4 + ...
297 : * with 0 <= a_i < p^(2^i).
298 : * 2. compute exp(a_i*p^(2^i)) using Taylor expansion and binary splitting */
299 : GEN
300 107227 : Zp_exp(GEN a, GEN p, ulong e)
301 : {
302 107227 : pari_sp av = avma;
303 107227 : ulong trunc, N = e, pp = itou_or_0(p);
304 107227 : GEN num, den, trunc_mod = NULL, denominator = gen_1, ans = gen_1;
305 107227 : GEN pe = powiu(p, e);
306 107227 : int pis2 = pp == 2;
307 :
308 : /* Where do we truncate the Taylor expansion */
309 107227 : if (!pis2) N += sdivsi(e, subis(p,2));
310 107227 : num = cgetg(N+2, t_VEC);
311 107226 : den = cgetg(N+2, t_VEC);
312 107227 : if (pis2) trunc = 4;
313 : else
314 : {
315 100700 : trunc = 2;
316 100700 : trunc_mod = sqri(p);
317 : }
318 : while(1)
319 198350 : {
320 305576 : GEN h, hpow, f = pis2? remi2n(a, trunc): modii(a, trunc_mod);
321 305575 : a = subii(a, f);
322 305576 : if (signe(f))
323 : {
324 229834 : ulong step = 1, i;
325 : /* Taylor expansion of exp(f), over Q for now */
326 229834 : gel(num,1) = gen_1;
327 229834 : gel(den,1) = gen_1;
328 1003563 : for (i = 2; i <= N+1; i++)
329 : {
330 773729 : gel(num,i) = gen_1;
331 773729 : gel(den,i) = utoipos(i-1);
332 : }
333 229834 : hpow = h = f;
334 : while(1)
335 : {
336 1308799 : for (i = 1; i <= N+1 - step; i += step << 1) {
337 773275 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
338 773267 : gel(num,i) = addii(gel(num,i), mulii(hpow, gel(num,i+step)));
339 773258 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
340 : }
341 535524 : step <<= 1; if (step > N) break;
342 305688 : hpow = sqri(hpow);
343 : }
344 :
345 : /* We simplify the fraction */
346 229836 : if (pp)
347 : {
348 229710 : GEN d = powuu(pp, factorial_lval(N, pp));
349 229710 : gel(num,1) = diviiexact(gel(num,1), d);
350 229709 : gel(den,1) = diviiexact(gel(den,1), d);
351 : }
352 : /* We add this contribution to exp(f) */
353 229836 : ans = Fp_mul(ans, gel(num,1), pe);
354 229837 : denominator = Fp_mul(denominator, gel(den,1), pe);
355 : }
356 :
357 305578 : if (trunc > e) break;
358 198350 : if (!pis2) trunc_mod = sqri(trunc_mod);
359 198350 : trunc <<= 1; N >>= 1;
360 : }
361 107228 : return gerepileuptoint(av, Zp_div(ans, denominator, p, e));
362 : }
363 :
364 : /***********************************************************************/
365 : /** **/
366 : /** QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL) **/
367 : /** **/
368 : /***********************************************************************/
369 : /* Setup for divide/conquer quadratic Hensel lift
370 : * a = set of k t_POL in Z[X] = factors over Fp (T=NULL) or Fp[Y]/(T)
371 : * V = set of products of factors built as follows
372 : * 1) V[1..k] = initial a
373 : * 2) iterate:
374 : * append to V the two smallest factors (minimal degree) in a, remove them
375 : * from a and replace them by their product [net loss for a = 1 factor]
376 : *
377 : * W = bezout coeffs W[i]V[i] + W[i+1]V[i+1] = 1
378 : *
379 : * link[i] = -j if V[i] = a[j]
380 : * j if V[i] = V[j] * V[j+1]
381 : * Arrays (link, V, W) pre-allocated for 2k - 2 elements */
382 : static void
383 181864 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
384 : {
385 181864 : long k = lg(a)-1;
386 : long i, j, s, minp, mind;
387 :
388 789725 : for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
389 :
390 425995 : for (j=1; j <= 2*k-5; j+=2,i++)
391 : {
392 244135 : minp = j;
393 244135 : mind = degpol(gel(V,j));
394 1958882 : for (s=j+1; s<i; s++)
395 1714749 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
396 :
397 244133 : swap(gel(V,j), gel(V,minp));
398 244133 : lswap(link[j], link[minp]);
399 :
400 244133 : minp = j+1;
401 244133 : mind = degpol(gel(V,j+1));
402 1714747 : for (s=j+2; s<i; s++)
403 1470616 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
404 :
405 244131 : swap(gel(V,j+1), gel(V,minp));
406 244131 : lswap(link[j+1], link[minp]);
407 :
408 244131 : gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
409 244131 : link[i] = j;
410 : }
411 :
412 607840 : for (j=1; j <= 2*k-3; j+=2)
413 : {
414 : GEN d, u, v;
415 425995 : d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
416 426003 : if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
417 425996 : d = gel(d,2);
418 425996 : if (!gequal1(d))
419 : {
420 387061 : if (typ(d)==t_POL)
421 : {
422 133017 : d = FpXQ_inv(d, T, p);
423 133010 : u = FqX_Fq_mul(u, d, T, p);
424 133008 : v = FqX_Fq_mul(v, d, T, p);
425 : }
426 : else
427 : {
428 254044 : d = Fp_inv(d, p);
429 254039 : u = FqX_Fp_mul(u, d, T,p);
430 254040 : v = FqX_Fp_mul(v, d, T,p);
431 : }
432 : }
433 425980 : gel(W,j) = u;
434 425980 : gel(W,j+1) = v;
435 : }
436 181845 : }
437 :
438 : /* au + bv = 1 (p0), ab = f (p0). Lift mod p1 = p0 pd (<= p0^2).
439 : * If noinv is set, don't lift the inverses u and v */
440 : static void
441 1521921 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
442 : {
443 1521921 : pari_sp av = avma;
444 1521921 : long space = lg(f) * lgefint(p1);
445 : GEN a2, b2, g, z, s, t;
446 1521921 : GEN a = gel(V,j), b = gel(V,j+1);
447 1521921 : GEN u = gel(W,j), v = gel(W,j+1);
448 :
449 1521921 : (void)new_chunk(space); /* HACK */
450 1521938 : g = ZX_sub(f, ZX_mul(a,b));
451 1521875 : g = ZX_Z_divexact(g, p0);
452 1521703 : g = FpX_red(g, pd);
453 1521821 : z = FpX_mul(v,g, pd);
454 1521902 : t = FpX_divrem(z,a, pd, &s);
455 1521902 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
456 1521840 : t = FpX_red(t, pd);
457 1521922 : t = ZX_Z_mul(t,p0);
458 1521731 : s = ZX_Z_mul(s,p0);
459 1521701 : set_avma(av);
460 1521739 : a2 = ZX_add(a,s);
461 1521882 : b2 = ZX_add(b,t);
462 :
463 : /* already reduced mod p1 = pd p0 */
464 1521905 : gel(V,j) = a2;
465 1521905 : gel(V,j+1) = b2;
466 1521905 : if (noinv) return;
467 :
468 1232448 : av = avma;
469 1232448 : (void)new_chunk(space); /* HACK */
470 1232463 : g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
471 1232417 : g = Z_ZX_sub(gen_1, g);
472 1232495 : g = ZX_Z_divexact(g, p0);
473 1232337 : g = FpX_red(g, pd);
474 1232381 : z = FpX_mul(v,g, pd);
475 1232448 : t = FpX_divrem(z,a, pd, &s);
476 1232474 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
477 1232413 : t = FpX_red(t, pd);
478 1232470 : t = ZX_Z_mul(t,p0);
479 1232343 : s = ZX_Z_mul(s,p0);
480 1232351 : set_avma(av);
481 1232371 : gel(W,j) = ZX_add(u,t);
482 1232425 : gel(W,j+1) = ZX_add(v,s);
483 : }
484 :
485 : static void
486 502087 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
487 : {
488 502087 : pari_sp av = avma;
489 502087 : const long n = degpol(T1), vT = varn(T1);
490 502091 : long space = lg(f) * lgefint(p1) * lg(T1);
491 : GEN a2, b2, g, z, s, t;
492 502091 : GEN a = gel(V,j), b = gel(V,j+1);
493 502091 : GEN u = gel(W,j), v = gel(W,j+1);
494 :
495 502091 : (void)new_chunk(space); /* HACK */
496 502077 : g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
497 502020 : g = FpXQX_red(g, T1, p1);
498 501988 : g = RgX_Rg_divexact(g, p0);
499 501717 : z = FpXQX_mul(v,g, Td,pd);
500 502066 : t = FpXQX_divrem(z,a, Td,pd, &s);
501 502110 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
502 502036 : t = Kronecker_to_ZXX(t, n, vT);
503 502041 : t = FpXQX_red(t, Td, pd);
504 501903 : t = RgX_Rg_mul(t,p0);
505 501965 : s = RgX_Rg_mul(s,p0);
506 501974 : set_avma(av);
507 :
508 501958 : a2 = RgX_add(a,s);
509 502014 : b2 = RgX_add(b,t);
510 : /* already reduced mod p1 = pd p0 */
511 502059 : gel(V,j) = a2;
512 502059 : gel(V,j+1) = b2;
513 502059 : if (noinv) return;
514 :
515 368423 : av = avma;
516 368423 : (void)new_chunk(space); /* HACK */
517 368419 : g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
518 368387 : g = Kronecker_to_ZXX(g, n, vT);
519 368443 : g = Rg_RgX_sub(gen_1, g);
520 368454 : g = FpXQX_red(g, T1, p1);
521 368409 : g = RgX_Rg_divexact(g, p0);
522 368283 : z = FpXQX_mul(v,g, Td,pd);
523 368445 : t = FpXQX_divrem(z,a, Td,pd, &s);
524 368458 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
525 368425 : t = Kronecker_to_ZXX(t, n, vT);
526 368441 : t = FpXQX_red(t, Td, pd);
527 368380 : t = RgX_Rg_mul(t,p0);
528 368399 : s = RgX_Rg_mul(s,p0);
529 368398 : set_avma(av);
530 368387 : gel(W,j) = RgX_add(u,t);
531 368397 : gel(W,j+1) = RgX_add(v,s);
532 : }
533 :
534 : /* v list of factors, w list of inverses. f = v[j] v[j+1]
535 : * Lift v[j] and v[j+1] mod p0 pd (possibly mod T), then all their divisors */
536 : static void
537 3722721 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
538 : GEN f, long j, int noinv)
539 : {
540 3722721 : if (j < 0) return;
541 1521687 : ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
542 1521855 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j) , link[j ], noinv);
543 1521896 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
544 : }
545 : static void
546 1136532 : ZpXQ_RecTreeLift(GEN link, GEN v, GEN w, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1,
547 : GEN f, long j, int noinv)
548 : {
549 1136532 : if (j < 0) return;
550 501956 : ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
551 502000 : ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j) , link[j ], noinv);
552 502044 : ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
553 : }
554 :
555 : /* Lift to precision p^e0.
556 : * a = modular factors of f mod (p,T) [possibly T=NULL]
557 : * OR a TreeLift structure [e, link, v, w]: go on lifting
558 : * flag = 0: standard.
559 : * flag = 1: return TreeLift structure */
560 : static GEN
561 181893 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
562 : {
563 181893 : long i, eold, e, k = lg(a) - 1;
564 : GEN E, v, w, link, penew, Tnew;
565 : ulong mask;
566 : pari_timer Ti;
567 :
568 181893 : if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
569 181893 : if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
570 181893 : if (e0 == 1) return a;
571 :
572 181865 : if (DEBUGLEVEL > 3) timer_start(&Ti);
573 181865 : if (typ(gel(a,1)) == t_INT)
574 : { /* a = TreeLift structure */
575 0 : e = itos(gel(a,1));
576 0 : link = gel(a,2);
577 0 : v = gel(a,3);
578 0 : w = gel(a,4);
579 : }
580 : else
581 : {
582 181865 : e = 1;
583 181865 : v = cgetg(2*k-2 + 1, t_VEC);
584 181866 : w = cgetg(2*k-2 + 1, t_VEC);
585 181867 : link=cgetg(2*k-2 + 1, t_VECSMALL);
586 181867 : BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
587 181859 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
588 : }
589 181859 : mask = quadratic_prec_mask(e0);
590 181862 : eold = 1;
591 181862 : penew = NULL;
592 181862 : Tnew = NULL;
593 181862 : if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
594 993601 : while (mask > 1)
595 : {
596 811743 : long enew = eold << 1;
597 811743 : if (mask & 1) enew--;
598 811743 : mask >>= 1;
599 811743 : if (enew >= e) { /* mask == 1: last iteration */
600 811739 : GEN peold = penew? penew: powiu(p, eold);
601 811751 : GEN Td = NULL, pd;
602 811751 : long d = enew - eold; /* = eold or eold-1 */
603 : /* lift from p^eold to p^enew */
604 811751 : pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
605 811729 : penew = mulii(peold,pd);
606 811601 : if (T) {
607 132543 : if (Tnew)
608 97853 : Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
609 : else
610 34690 : Td = FpX_red(T, peold);
611 132555 : Tnew = FpX_red(T, penew);
612 132554 : ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
613 : (flag == 0 && mask == 1));
614 : }
615 : else
616 679058 : ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
617 : (flag == 0 && mask == 1));
618 811739 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
619 : }
620 811739 : eold = enew;
621 : }
622 :
623 181858 : if (flag)
624 1736 : E = mkvec4(utoipos(e0), link, v, w);
625 : else
626 : {
627 180122 : E = cgetg(k+1, t_VEC);
628 1026344 : for (i = 1; i <= 2*k-2; i++)
629 : {
630 846219 : long t = link[i];
631 846219 : if (t < 0) gel(E,-t) = gel(v,i);
632 : }
633 : }
634 181861 : return E;
635 : }
636 :
637 : /* Q list of (coprime, monic) factors of pol mod (T,p). Lift mod p^e = pe.
638 : * T may be NULL */
639 : GEN
640 219255 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
641 : {
642 219255 : pari_sp av = avma;
643 219255 : pol = FpX_normalize(pol, pe);
644 219255 : if (lg(Q) == 2) return mkvec(pol);
645 145459 : return gerepilecopy(av, MultiLift(pol, Q, NULL, p, e, 0));
646 : }
647 :
648 : GEN
649 34699 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
650 : {
651 34699 : pari_sp av = avma;
652 34699 : pol = FpXQX_normalize(pol, T, pe);
653 34698 : if (lg(Q) == 2) return mkvec(pol);
654 34698 : return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
655 : }
656 :
657 : GEN
658 46515 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
659 46515 : { return T ? ZpXQX_liftfact(f, a, T, pe, p, e): ZpX_liftfact(f, a, pe, p, e); }
660 : GEN
661 42 : ZqX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
662 42 : { return T ? ZpXQX_liftroot(f, a,T , p, e): ZpX_liftroot(f, a, p, e); }
663 :
664 : /* U = NULL treated as 1 */
665 : static void
666 7504 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
667 : {
668 : GEN Q, R;
669 7504 : if (j < 0) return;
670 :
671 2884 : Q = FpX_mul(gel(v,j), gel(w,j), pe);
672 2884 : if (U)
673 : {
674 1148 : Q = FpXQ_mul(Q, U, f, pe);
675 1148 : R = FpX_sub(U, Q, pe);
676 : }
677 : else
678 1736 : R = Fp_FpX_sub(gen_1, Q, pe);
679 2884 : gel(w,j+1) = Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
680 2884 : gel(w,j) = R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
681 2884 : BezoutPropagate(link, v, w, pe, R, f, link[j ]);
682 2884 : BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
683 : }
684 :
685 : /* as above, but return the Bezout coefficients for the lifted modular factors
686 : * U[i] = 1 mod Qlift[i]
687 : * 0 mod Qlift[j], j != i */
688 : GEN
689 1757 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
690 : {
691 1757 : pari_sp av = avma;
692 : GEN E, link, v, w, pe;
693 1757 : long i, k = lg(Q)-1;
694 1757 : if (k == 1) retmkvec(pol_1(varn(pol)));
695 1736 : pe = powiu(p, e);
696 1736 : pol = FpX_normalize(pol, pe);
697 1736 : E = MultiLift(pol, Q, NULL, p, e, 1);
698 1736 : link = gel(E,2);
699 1736 : v = gel(E,3);
700 1736 : w = gel(E,4);
701 1736 : BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
702 1736 : E = cgetg(k+1, t_VEC);
703 7504 : for (i = 1; i <= 2*k-2; i++)
704 : {
705 5768 : long t = link[i];
706 5768 : if (t < 0) E[-t] = w[i];
707 : }
708 1736 : return gerepilecopy(av, E);
709 : }
710 :
711 : /* Front-end for ZpX_liftfact:
712 : lift the factorization of pol mod p given by L to p^N (if possible) */
713 : GEN
714 28 : polhensellift(GEN pol, GEN L, GEN Tp, long N)
715 : {
716 : GEN T, p;
717 : long i, l;
718 28 : pari_sp av = avma;
719 : void (*chk)(GEN, const char*);
720 :
721 28 : if (typ(pol) != t_POL) pari_err_TYPE("polhensellift",pol);
722 28 : RgX_check_ZXX(pol, "polhensellift");
723 28 : if (!is_vec_t(typ(L)) || lg(L) < 3) pari_err_TYPE("polhensellift",L);
724 28 : if (N < 1) pari_err_DOMAIN("polhensellift", "precision", "<", gen_1,stoi(N));
725 28 : if (!ff_parse_Tp(Tp, &T, &p, 0)) pari_err_TYPE("polhensellift",Tp);
726 28 : chk = T? RgX_check_ZXX: RgX_check_ZX;
727 28 : l = lg(L); L = leafcopy(L);
728 70 : for (i = 1; i < l; i++)
729 : {
730 49 : GEN q = gel(L,i);
731 49 : if (typ(q) != t_POL) gel(L,i) = scalar_ZX_shallow(q, varn(pol));
732 49 : else chk(q, "polhensellift");
733 : }
734 21 : return gerepilecopy(av, ZqX_liftfact(pol, L, T, powiu(p,N), p, N));
735 : }
736 :
737 : static GEN
738 83614 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
739 : {
740 83614 : long i,l = lg(x);
741 83614 : GEN r = cgetg(l,t_COL);
742 302204 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
743 83608 : return r;
744 : }
745 :
746 : static GEN
747 83559 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
748 : {
749 83559 : pari_sp av = avma;
750 83559 : GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
751 83558 : return gerepileupto(av, FqV_roots_from_deg1(y, NULL, q));
752 : }
753 :
754 : GEN
755 83420 : ZpX_roots(GEN F, GEN p, long e)
756 : {
757 83420 : pari_sp av = avma;
758 83420 : GEN q = powiu(p, e);
759 83420 : GEN f = FpX_normalize(F, p);
760 83420 : GEN g = FpX_normalize(FpX_split_part(f, p), p);
761 : GEN S;
762 83421 : if (degpol(g) < degpol(f))
763 : {
764 43870 : GEN h = FpX_div(f, g, p);
765 43870 : F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
766 : }
767 83421 : S = FpX_roots(g, p);
768 83420 : return gerepileupto(av, ZpX_liftroots_full(F, S, q, p, e));
769 : }
770 :
771 : static GEN
772 56 : ZpXQX_liftroots_full(GEN f, GEN S, GEN T, GEN q, GEN p, long e)
773 : {
774 56 : pari_sp av = avma;
775 56 : GEN y = ZpXQX_liftfact(f, deg1_from_roots(S, varn(f)), T, q, p, e);
776 56 : return gerepileupto(av, FqV_roots_from_deg1(y, T, q));
777 : }
778 :
779 : GEN
780 56 : ZpXQX_roots(GEN F, GEN T, GEN p, long e)
781 : {
782 56 : pari_sp av = avma;
783 56 : GEN q = powiu(p, e);
784 56 : GEN f = FpXQX_normalize(F, T, p);
785 56 : GEN g = FpXQX_normalize(FpXQX_split_part(f, T, p), T, p);
786 : GEN S;
787 56 : if (degpol(g) < degpol(f))
788 : {
789 21 : GEN h = FpXQX_div(f, g, T, p);
790 21 : F = gel(ZpXQX_liftfact(F, mkvec2(g, h), T, q, p, e), 1);
791 : }
792 56 : S = FpXQX_roots(g, T, p);
793 56 : return gerepileupto(av, ZpXQX_liftroots_full(F, S, T, q, p, e));
794 : }
795 :
796 : GEN
797 3598 : ZqX_roots(GEN F, GEN T, GEN p, long e)
798 : {
799 3598 : return T ? ZpXQX_roots(F, T, p, e): ZpX_roots(F, p, e);
800 : }
801 :
802 : /* SPEC:
803 : * p is a t_INT > 1, e >= 1
804 : * f is a ZX with leading term prime to p.
805 : * a is a simple root mod l for all l|p.
806 : * Return roots of f mod p^e, as integers (implicitly mod p^e)
807 : * STANDARD USE: p is a prime power */
808 : GEN
809 3696 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
810 : {
811 3696 : pari_sp av = avma;
812 3696 : GEN q = p, fr, W;
813 : ulong mask;
814 :
815 3696 : a = modii(a,q);
816 3696 : if (e == 1) return a;
817 3696 : mask = quadratic_prec_mask(e);
818 3696 : fr = FpX_red(f,q);
819 3696 : W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
820 : for(;;)
821 : {
822 8015 : q = sqri(q);
823 8015 : if (mask & 1) q = diviiexact(q, p);
824 8015 : mask >>= 1;
825 8015 : fr = FpX_red(f,q);
826 8015 : a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
827 8015 : if (mask == 1) return gerepileuptoint(av, a);
828 4319 : W = Fp_sub(shifti(W,1), Fp_mul(Fp_sqr(W,q), FpX_eval(ZX_deriv(fr),a,q), q), q);
829 : }
830 : }
831 :
832 : GEN
833 139 : ZpX_liftroots(GEN f, GEN S, GEN p, long e)
834 : {
835 139 : long i, n = lg(S)-1, d = degpol(f);
836 : GEN r;
837 139 : if (n == d) return ZpX_liftroots_full(f, S, powiu(p, e), p, e);
838 0 : r = cgetg(n+1, typ(S));
839 0 : for (i=1; i <= n; i++)
840 0 : gel(r,i) = ZpX_liftroot(f, gel(S,i), p, e);
841 0 : return r;
842 : }
843 :
844 : GEN
845 119 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
846 : {
847 119 : pari_sp av = avma, av2;
848 119 : GEN pv = p, q, W, df, Tq, fr, dfr;
849 : ulong mask;
850 119 : a = Fq_red(a, T, p);
851 119 : if (e <= v+1) return a;
852 119 : df = RgX_deriv(f);
853 119 : if (v) { pv = powiu(p,v); df = ZXX_Z_divexact(df, pv); }
854 119 : mask = quadratic_prec_mask(e-v);
855 119 : Tq = FpXT_red(T, p); dfr = FpXQX_red(df, Tq, p);
856 119 : W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
857 119 : q = p;
858 119 : av2 = avma;
859 : for (;;)
860 161 : {
861 : GEN u, fa, qv, q2v, q2, Tq2;
862 280 : q2 = q; q = sqri(q);
863 280 : if (mask & 1) q = diviiexact(q,p);
864 280 : mask >>= 1;
865 280 : if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
866 280 : else { qv = q; q2v = q2; }
867 280 : Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
868 280 : fr = FpXQX_red(f, Tq, qv);
869 280 : fa = FqX_eval(fr, a, Tq, qv);
870 280 : fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
871 280 : a = Fq_sub(a, Fq_mul(Fq_mul(W,fa,Tq2,q2v),q2, Tq,qv), Tq, qv);
872 280 : if (mask == 1) return gerepileupto(av, a);
873 161 : dfr = FpXQX_red(df, Tq, q);
874 161 : u = Fq_sub(Fq_mul(W,FqX_eval(dfr,a,Tq,q),Tq,q),gen_1,Tq,q);
875 161 : u = typ(u)==t_INT? diviiexact(u,q2): ZX_Z_divexact(u,q2);
876 161 : W = Fq_sub(W, Fq_mul(Fq_mul(u,W,Tq2,q2),q2, Tq,q),Tq,q);
877 161 : if (gc_needed(av2,2))
878 : {
879 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_liftroot, e = %ld", e);
880 0 : gerepileall(av2, 3, &a, &W, &q);
881 : }
882 : }
883 : }
884 :
885 : GEN
886 119 : ZpXQX_liftroot(GEN f, GEN a, GEN T, GEN p, long e) { return ZpXQX_liftroot_vald(f,a,0,T,p,e); }
887 :
888 : GEN
889 0 : ZpXQX_liftroots(GEN f, GEN S, GEN T, GEN p, long e)
890 : {
891 0 : long i, n = lg(S)-1, d = degpol(f);
892 : GEN r;
893 0 : if (n == d) return ZpXQX_liftroots_full(f, S, T, powiu(p, e), p, e);
894 0 : r = cgetg(n+1, typ(S));
895 0 : for (i=1; i <= n; i++)
896 0 : gel(r,i) = ZpXQX_liftroot(f, gel(S,i), T, p, e);
897 0 : return r;
898 : }
899 :
900 : /* Compute (x-1)/(x+1)/p^k */
901 : static GEN
902 21080 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
903 : {
904 21080 : pari_sp av = avma;
905 21080 : long vT = get_FpX_var(T);
906 : GEN bn, bdi;
907 21080 : GEN bd = ZX_Z_add(x, gen_1);
908 21080 : if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
909 : {
910 7150 : bn = ZX_shifti(x,-(k+1));
911 7150 : bdi= ZpXQ_invlift(ZX_shifti(bd ,-1), pol_1(vT), T, p, e);
912 : }
913 : else
914 : {
915 13930 : bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
916 13930 : bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
917 : }
918 21080 : return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
919 : }
920 :
921 : /* Assume p odd, a = 1 [p], return log(a) */
922 : GEN
923 21080 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
924 : {
925 21080 : pari_sp av = avma;
926 : pari_timer ti;
927 21080 : long is2 = absequaliu(p,2);
928 21080 : ulong pp = is2 ? 0: itou_or_0(p);
929 21080 : double lp = is2 ? 1: pp ? log2(pp): expi(p);
930 21080 : long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
931 : GEN ak, s, b, pol;
932 21080 : long e = is2 ? N-1: N;
933 21080 : long i, l = (e-2)/(2*(k+is2));
934 21080 : GEN pe = powiu(p,e);
935 21080 : GEN TNk, pNk = powiu(p,N+k);
936 21080 : if( DEBUGLEVEL>=3) timer_start(&ti);
937 21080 : TNk = FpX_get_red(get_FpX_mod(T), pNk);
938 21080 : ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
939 21080 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
940 21080 : b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
941 21080 : if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
942 21080 : pol= cgetg(l+3,t_POL);
943 21080 : pol[1] = evalsigne(1)|evalvarn(0);
944 58175 : for(i=0; i<=l; i++)
945 : {
946 : GEN g;
947 37096 : ulong z = 2*i+1;
948 37096 : if (pp)
949 : {
950 26222 : long w = u_lvalrem(z, pp, &z);
951 26222 : g = powuu(pp,2*i*k-w);
952 : }
953 10874 : else g = powiu(p,2*i*k);
954 37092 : gel(pol,i+2) = Fp_divu(g, z,pe);
955 : }
956 21079 : if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
957 21079 : s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T, pe);
958 21080 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
959 21080 : s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
960 21080 : return gerepileupto(av, is2? s: FpX_red(s, pe));
961 : }
962 :
963 : /***********************************************************************/
964 : /** **/
965 : /** Generic quadratic hensel lift over Zp[X] **/
966 : /** **/
967 : /***********************************************************************/
968 : /* q = p^N */
969 :
970 : GEN
971 0 : gen_ZpM_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
972 : GEN lin(void *E, GEN F, GEN d, GEN q),
973 : GEN invl(void *E, GEN d))
974 : {
975 0 : pari_sp av = avma;
976 : long N2, M;
977 : GEN VN2, V2, VM, bil;
978 : GEN q2, qM;
979 0 : V = FpM_red(V, q);
980 0 : if (N == 1) return invl(E, V);
981 0 : N2 = (N + 1)>>1; M = N - N2;
982 0 : F = FpM_red(F, q);
983 0 : qM = powiu(p, M);
984 0 : q2 = M == N2? qM: mulii(qM, p);
985 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
986 0 : VN2 = gen_ZpM_Dixon(F, V, q2, p, N2, E, lin, invl);
987 0 : bil = lin(E, F, VN2, q);
988 0 : V2 = ZM_Z_divexact(ZM_sub(V, bil), q2);
989 0 : VM = gen_ZpM_Dixon(F, V2, qM, p, M, E, lin, invl);
990 0 : return gerepileupto(av, FpM_red(ZM_add(VN2, ZM_Z_mul(VM, q2)), q));
991 : }
992 :
993 : GEN
994 199829 : gen_ZpX_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
995 : GEN lin(void *E, GEN F, GEN d, GEN q),
996 : GEN invl(void *E, GEN d))
997 : {
998 199829 : pari_sp av = avma;
999 : long N2, M;
1000 : GEN VN2, V2, VM, bil;
1001 : GEN q2, qM;
1002 199829 : V = FpX_red(V, q);
1003 199829 : if (N == 1) return invl(E, V);
1004 45458 : N2 = (N + 1)>>1; M = N - N2;
1005 45458 : F = FpXT_red(F, q);
1006 45458 : qM = powiu(p, M);
1007 45458 : q2 = M == N2? qM: mulii(qM, p);
1008 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
1009 45458 : VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
1010 45458 : bil = lin(E, F, VN2, q);
1011 45458 : V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
1012 45458 : VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
1013 45458 : return gerepileupto(av, FpX_red(ZX_add(VN2, ZX_Z_mul(VM, q2)), q));
1014 : }
1015 :
1016 : GEN
1017 948418 : gen_ZpM_Newton(GEN x, GEN p, long n, void *E,
1018 : GEN eval(void *E, GEN f, GEN q),
1019 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1020 : {
1021 948418 : pari_sp ltop = avma, av;
1022 948418 : long N = 1, N2, M;
1023 : long mask;
1024 948418 : GEN q = p;
1025 948418 : if (n == 1) return gcopy(x);
1026 948418 : mask = quadratic_prec_mask(n);
1027 948425 : av = avma;
1028 1923510 : while (mask > 1)
1029 : {
1030 : GEN qM, q2, v, V;
1031 975090 : N2 = N; N <<= 1;
1032 975090 : q2 = q;
1033 975090 : if (mask&1UL) { /* can never happen when q2 = p */
1034 18644 : N--; M = N2-1;
1035 18644 : qM = diviiexact(q2,p); /* > 1 */
1036 18644 : q = mulii(qM,q2);
1037 : } else {
1038 956446 : M = N2;
1039 956446 : qM = q2;
1040 956446 : q = sqri(q2);
1041 : }
1042 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1043 975069 : mask >>= 1;
1044 975069 : v = eval(E, x, q);
1045 975096 : V = ZM_Z_divexact(gel(v,1), q2);
1046 975074 : x = FpM_sub(x, ZM_Z_mul(invd(E, V, v, qM, M), q2), q);
1047 975083 : if (gc_needed(av, 1))
1048 : {
1049 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Newton");
1050 0 : gerepileall(av, 2, &x, &q);
1051 : }
1052 : }
1053 948420 : return gerepileupto(ltop, x);
1054 : }
1055 :
1056 : static GEN
1057 975082 : _ZpM_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1058 : {
1059 : (void)E; (void)M;
1060 975082 : return FpM_mul(V, gel(v,2), q);
1061 : }
1062 :
1063 : static GEN
1064 975071 : _ZpM_eval(void *E, GEN x, GEN q)
1065 : {
1066 975071 : GEN f = RgM_Rg_sub_shallow(FpM_mul(x, FpM_red((GEN) E, q), q), gen_1);
1067 975079 : return mkvec2(f, x);
1068 : }
1069 :
1070 : GEN
1071 948410 : ZpM_invlift(GEN M, GEN C, GEN p, long n)
1072 : {
1073 948410 : return gen_ZpM_Newton(C, p, n, (void *)M, _ZpM_eval, _ZpM_invd);
1074 : }
1075 :
1076 : GEN
1077 266338 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
1078 : GEN eval(void *E, GEN f, GEN q),
1079 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1080 : {
1081 266338 : pari_sp ltop = avma, av;
1082 266338 : long N = 1, N2, M;
1083 : long mask;
1084 266338 : GEN q = p;
1085 266338 : if (n == 1) return gcopy(x);
1086 262509 : mask = quadratic_prec_mask(n);
1087 262509 : av = avma;
1088 929119 : while (mask > 1)
1089 : {
1090 : GEN qM, q2, v, V;
1091 666610 : N2 = N; N <<= 1;
1092 666610 : q2 = q;
1093 666610 : if (mask&1UL) { /* can never happen when q2 = p */
1094 276637 : N--; M = N2-1;
1095 276637 : qM = diviiexact(q2,p); /* > 1 */
1096 276636 : q = mulii(qM,q2);
1097 : } else {
1098 389973 : M = N2;
1099 389973 : qM = q2;
1100 389973 : q = sqri(q2);
1101 : }
1102 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1103 666598 : mask >>= 1;
1104 666598 : v = eval(E, x, q);
1105 666611 : V = ZX_Z_divexact(gel(v,1), q2);
1106 666601 : x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
1107 666610 : if (gc_needed(av, 1))
1108 : {
1109 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
1110 0 : gerepileall(av, 2, &x, &q);
1111 : }
1112 : }
1113 262509 : return gerepileupto(ltop, x);
1114 : }
1115 :
1116 : struct _ZpXQ_inv
1117 : {
1118 : GEN T, a, p ,n;
1119 : };
1120 :
1121 : static GEN
1122 522529 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1123 : {
1124 522529 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1125 522529 : GEN Tq = FpXT_red(d->T, q);
1126 : (void)M;
1127 522529 : return FpXQ_mul(V, gel(v,2), Tq, q);
1128 : }
1129 :
1130 : static GEN
1131 522529 : _inv_eval(void *E, GEN x, GEN q)
1132 : {
1133 522529 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1134 522529 : GEN Tq = FpXT_red(d->T, q);
1135 522529 : GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
1136 522529 : return mkvec2(f, x);
1137 : }
1138 :
1139 : GEN
1140 201323 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
1141 : {
1142 : struct _ZpXQ_inv d;
1143 201323 : d.a = a; d.T = T; d.p = p;
1144 201323 : return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
1145 : }
1146 :
1147 : GEN
1148 180243 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
1149 : {
1150 180243 : pari_sp av=avma;
1151 : GEN ai;
1152 180243 : if (lgefint(p)==3)
1153 : {
1154 180215 : ulong pp = p[2];
1155 180215 : ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
1156 : } else
1157 28 : ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
1158 180243 : return gerepileupto(av, ZpXQ_invlift(a, ai, T, p, e));
1159 : }
1160 :
1161 : GEN
1162 31584 : ZpXQ_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
1163 : {
1164 31584 : return FpXQ_mul(a, ZpXQ_inv(b, T, p, e), T, q);
1165 : }
1166 :
1167 : GEN
1168 146440 : ZpXQX_divrem(GEN x, GEN Sp, GEN T, GEN q, GEN p, long e, GEN *pr)
1169 : {
1170 146440 : pari_sp av = avma;
1171 146440 : GEN S = get_FpXQX_mod(Sp);
1172 146440 : GEN b = leading_coeff(S), bi;
1173 : GEN S2, Q;
1174 146440 : if (typ(b)==t_INT) return FpXQX_divrem(x, Sp, T, q, pr);
1175 146440 : bi = ZpXQ_inv(b, T, p, e);
1176 146440 : S2 = FqX_Fq_mul_to_monic(S, bi, T, q);
1177 146440 : Q = FpXQX_divrem(x, S2, T, q, pr);
1178 146440 : if (pr==ONLY_DIVIDES && !Q) { set_avma(av); return NULL; }
1179 146440 : if (pr==ONLY_REM || pr==ONLY_DIVIDES) return gerepileupto(av, Q);
1180 146440 : Q = FpXQX_FpXQ_mul(Q, bi, T, q);
1181 146440 : return gc_all(av, 2, &Q, pr);
1182 : }
1183 :
1184 : GEN
1185 623 : ZpXQX_digits(GEN x, GEN B, GEN T, GEN q, GEN p, long e)
1186 : {
1187 623 : pari_sp av = avma;
1188 623 : GEN b = leading_coeff(B), bi;
1189 : GEN B2, P, V, W;
1190 : long i, lV;
1191 623 : if (typ(b)==t_INT) return FpXQX_digits(x, B, T, q);
1192 623 : bi = ZpXQ_inv(b, T, p, e);
1193 623 : B2 = FqX_Fq_mul_to_monic(B, bi, T, q);
1194 623 : V = FpXQX_digits(x, B2, T, q);
1195 623 : lV = lg(V)-1;
1196 623 : P = FpXQ_powers(bi, lV-1, T, q);
1197 623 : W = cgetg(lV+1, t_VEC);
1198 11200 : for(i=1; i<=lV; i++)
1199 10577 : gel(W, i) = FpXQX_FpXQ_mul(gel(V,i), gel(P, i), T, q);
1200 623 : return gerepileupto(av, W);
1201 : }
1202 :
1203 : struct _ZpXQ_sqrtn
1204 : {
1205 : GEN T, a, n, ai;
1206 : };
1207 :
1208 : static GEN
1209 1806 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
1210 : {
1211 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1212 1806 : GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
1213 : (void)M;
1214 1806 : return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
1215 : }
1216 :
1217 : static GEN
1218 1806 : _sqrtn_eval(void *E, GEN x, GEN q)
1219 : {
1220 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1221 1806 : GEN Tq = FpX_red(d->T, q);
1222 1806 : GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
1223 1806 : return mkvec2(f, x);
1224 : }
1225 :
1226 : GEN
1227 1155 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1228 : {
1229 : struct _ZpXQ_sqrtn d;
1230 1155 : d.a = a; d.T = T; d.n = n;
1231 1155 : d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
1232 1155 : return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
1233 : }
1234 :
1235 : static GEN
1236 0 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
1237 :
1238 : GEN
1239 0 : Zq_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1240 : {
1241 0 : return T? ZpXQ_sqrtnlift(to_ZX(a,varn(T)), n, to_ZX(x,varn(T)), T, p, e)
1242 0 : : Zp_sqrtnlift(a, n, x, p, e);
1243 : }
1244 :
1245 : GEN
1246 0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
1247 : {
1248 0 : pari_sp av = avma;
1249 0 : GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
1250 0 : if (!z) return NULL;
1251 0 : if (e <= 1) return gerepileupto(av, z);
1252 0 : return gerepileupto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
1253 : }
1254 :
1255 : GEN
1256 33662 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
1257 : GEN early(void *E, GEN x, GEN q))
1258 : {
1259 33662 : pari_sp ltop = avma, av;
1260 : long N, r;
1261 : long mask;
1262 : GEN q2, q, W, Q, Tq2, Tq, Pq;
1263 : pari_timer ti;
1264 33662 : T = FpX_get_red(T, powiu(p, n));
1265 33661 : if (n == 1) return gcopy(S);
1266 33661 : mask = quadratic_prec_mask(n);
1267 33662 : av = avma;
1268 33662 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1269 33658 : if (DEBUGLEVEL > 3) timer_start(&ti);
1270 33658 : Tq = FpXT_red(T,q);
1271 33662 : Tq2 = FpXT_red(Tq,q2);
1272 33662 : Pq = FpX_red(P,q);
1273 33663 : W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
1274 33663 : Q = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
1275 33662 : r = brent_kung_optpow(degpol(P), 4, 3);
1276 33662 : if (DEBUGLEVEL > 3)
1277 0 : err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
1278 : for (;;)
1279 38834 : {
1280 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
1281 72496 : H = FpXQ_mul(W, Q, Tq2, q2);
1282 72495 : Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
1283 72496 : if (DEBUGLEVEL > 3)
1284 0 : timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
1285 72496 : if (mask==1)
1286 6538 : return gerepileupto(ltop, Sq);
1287 65958 : if (early)
1288 : {
1289 64299 : GEN Se = early(E, Sq, q);
1290 64299 : if (Se) return gerepileupto(ltop, Se);
1291 : }
1292 38834 : qq = sqri(q); N <<= 1;
1293 38832 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1294 38833 : mask >>= 1;
1295 38833 : Pqq = FpX_red(P, qq);
1296 38833 : Tqq = FpXT_red(T, qq);
1297 38833 : Spow = FpXQ_powers(Sq, r, Tqq, qq);
1298 38834 : Q = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
1299 38834 : dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
1300 38834 : Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
1301 38832 : Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
1302 38834 : Wq = FpX_sub(W, Wq, q);
1303 38834 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
1304 38834 : if (gc_needed(av, 1))
1305 : {
1306 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_liftroot");
1307 0 : gerepileall(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
1308 : }
1309 : }
1310 : }
1311 :
1312 : GEN
1313 1729 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
1314 : {
1315 1729 : return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
1316 : }
1317 :
1318 : GEN
1319 301 : ZpX_Frobenius(GEN T, GEN p, long e)
1320 : {
1321 301 : return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
1322 : }
1323 :
1324 : GEN
1325 147 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
1326 : {
1327 147 : pari_sp av = avma;
1328 147 : GEN xp = ZpX_Frobenius(T, p, e);
1329 147 : GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
1330 147 : return gerepilecopy(av, gel(z,2));
1331 : }
1332 :
1333 : GEN
1334 0 : ZpXQX_ZpXQXQ_liftroot(GEN P, GEN S, GEN U, GEN T, GEN p, long n)
1335 : {
1336 0 : pari_sp ltop = avma, av;
1337 : long N, r;
1338 : long mask;
1339 : GEN qn, q2, q, W, Q, Tq2, Tq, Pq, Uq, Uq2;
1340 : pari_timer ti;
1341 0 : qn = powiu(p, n);
1342 0 : T = FpX_get_red(T, qn);
1343 0 : U = FpXQX_get_red(U, T, qn);
1344 0 : if (n == 1) return gcopy(S);
1345 0 : mask = quadratic_prec_mask(n);
1346 0 : av = avma;
1347 0 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1348 0 : if (DEBUGLEVEL > 3) timer_start(&ti);
1349 0 : Tq = FpXT_red(T,q);
1350 0 : Uq = FpXQXT_red(U, Tq, q);
1351 0 : Tq2 = FpXT_red(Tq, q2);
1352 0 : Uq2 = FpXQXT_red(U, Tq2, q2);
1353 0 : Pq = FpXQX_red(P, Tq, q);
1354 0 : W = FpXQXQ_inv(FpXQX_FpXQXQ_eval(FpXX_deriv(P,q2), S, Uq2, Tq2, q2), Uq2, Tq2, q2);
1355 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQ_eval(Pq, S, Uq, Tq, q), q2);
1356 0 : r = brent_kung_optpow(degpol(P), 4, 3);
1357 0 : if (DEBUGLEVEL > 3)
1358 0 : err_printf("ZpXQX_ZpXQXQ_liftroot: lifting to prec %ld\n",n);
1359 : for (;;)
1360 0 : {
1361 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq, Uqq;
1362 0 : H = FpXQXQ_mul(W, Q, Uq2, Tq2, q2);
1363 0 : Sq = FpXX_sub(S, ZXX_Z_mul(H, q2), q);
1364 0 : if (DEBUGLEVEL > 3)
1365 0 : timer_printf(&ti,"ZpXQX_ZpXQXQ_liftroot: reaching prec %ld",N);
1366 0 : if (mask==1)
1367 0 : return gerepileupto(ltop, Sq);
1368 0 : qq = sqri(q); N <<= 1;
1369 0 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1370 0 : mask >>= 1;
1371 0 : Tqq = FpXT_red(T, qq);
1372 0 : Uqq = FpXQXT_red(U, Tqq, qq);
1373 0 : Pqq = FpXQX_red(P, Tqq, qq);
1374 0 : Spow = FpXQXQ_powers(Sq, r, Uqq, Tqq, qq);
1375 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQV_eval(Pqq, Spow, Uqq, Tqq, qq), q);
1376 0 : dP = FpXQX_FpXQXQV_eval(FpXX_deriv(Pq, q), FpXQXV_red(Spow, Tq, q), Uq, Tq, q);
1377 0 : Wq = ZXX_Z_divexact(gsub(FpXQXQ_mul(W, dP, Uq, Tq, q), gen_1), q2);
1378 0 : Wq = ZXX_Z_mul(FpXQXQ_mul(W, Wq, Uq2, Tq2, q2), q2);
1379 0 : Wq = FpXX_sub(W, Wq, q);
1380 0 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Uq2 = Uq; Uq = Uqq; Pq = Pqq;
1381 0 : if (gc_needed(av, 1))
1382 : {
1383 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_ZpXQXQ_liftroot");
1384 0 : gerepileall(av, 10, &S, &W, &Q, &Uq2, &Uq, &Tq2, &Tq, &Pq, &q, &q2);
1385 : }
1386 : }
1387 : }
1388 :
1389 : GEN
1390 35 : ZqX_ZqXQ_liftroot(GEN f, GEN a, GEN P, GEN T, GEN p, long e)
1391 35 : { return T ? ZpXQX_ZpXQXQ_liftroot(f, a, P, T , p, e): ZpX_ZpXQ_liftroot(f, a, P, p, e); }
1392 :
1393 : /* Canonical lift of polynomial */
1394 :
1395 11130 : static GEN _can_invl(void *E, GEN V) {(void) E; return V; }
1396 :
1397 3654 : static GEN _can_lin(void *E, GEN F, GEN V, GEN q)
1398 : {
1399 3654 : GEN v = RgX_splitting(V, 3);
1400 : (void) E;
1401 3654 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1402 : }
1403 :
1404 : static GEN
1405 7476 : _can_iter(void *E, GEN f, GEN q)
1406 : {
1407 7476 : GEN h = RgX_splitting(f,3);
1408 7476 : GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
1409 7476 : GEN h12 = ZX_mul(gel(h,1), gel(h,2));
1410 7476 : GEN h13 = ZX_mul(gel(h,1), gel(h,3));
1411 7476 : GEN h23 = ZX_mul(gel(h,2), gel(h,3));
1412 7476 : GEN h1c = ZX_mul(gel(h,1), h1s);
1413 7476 : GEN h3c = ZX_mul(gel(h,3), h3s);
1414 7476 : GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
1415 7476 : GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
1416 : (void) E;
1417 7476 : return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
1418 : }
1419 :
1420 : static GEN
1421 7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
1422 : {
1423 7476 : GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
1424 7476 : GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
1425 7476 : GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
1426 : ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
1427 : (void)E;
1428 7476 : return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
1429 : _can_lin, _can_invl);
1430 : }
1431 :
1432 : static GEN
1433 3717 : F3x_frobeniuslift(GEN P, long n)
1434 3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
1435 :
1436 29736 : static GEN _can5_invl(void *E, GEN V) {(void) E; return V; }
1437 :
1438 9107 : static GEN _can5_lin(void *E, GEN F, GEN V, GEN q)
1439 : {
1440 9107 : ulong p = *(ulong*)E;
1441 9107 : GEN v = RgX_splitting(V, p);
1442 9107 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1443 : }
1444 :
1445 : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
1446 : static GEN
1447 62146 : _shift(GEN P, long n, ulong p, long v)
1448 : {
1449 62146 : long i, l=lg(P);
1450 62146 : GEN r = cgetg(l,t_POL); r[1] = P[1];
1451 483238 : for(i=2;i<l;i++)
1452 : {
1453 421092 : long s = n*(i-2)%p;
1454 421092 : GEN ci = gel(P,i);
1455 421092 : if (typ(ci)==t_INT)
1456 104979 : gel(r,i) = monomial(ci, s, v);
1457 : else
1458 316113 : gel(r,i) = RgX_rotate_shallow(ci, s, p);
1459 : }
1460 62146 : return FpXX_renormalize(r, l);
1461 : }
1462 :
1463 : struct _can_mul
1464 : {
1465 : GEN T, q;
1466 : ulong p;
1467 : };
1468 :
1469 : static GEN
1470 41517 : _can5_mul(void *E, GEN A, GEN B)
1471 : {
1472 41517 : struct _can_mul *d = (struct _can_mul *)E;
1473 41517 : GEN a = gel(A,1), b = gel(B,1);
1474 41517 : long n = itos(gel(A,2));
1475 41517 : GEN bn = _shift(b, n, d->p, get_FpX_var(d->T));
1476 41517 : GEN c = FpXQX_mul(a, bn, d->T, d->q);
1477 41517 : return mkvec2(c, addii(gel(A,2), gel(B,2)));
1478 : }
1479 :
1480 : static GEN
1481 41349 : _can5_sqr(void *E, GEN A)
1482 : {
1483 41349 : return _can5_mul(E,A,A);
1484 : }
1485 :
1486 : static GEN
1487 20629 : _can5_iter(void *E, GEN f, GEN q)
1488 : {
1489 20629 : pari_sp av = avma;
1490 : struct _can_mul D;
1491 20629 : ulong p = *(ulong*)E;
1492 20629 : long i, vT = fetch_var();
1493 : GEN N, P, d, V, fs;
1494 20629 : D.q = q; D.T = ZX_Z_sub(pol_xn(p,vT),gen_1);
1495 20629 : D.p = p;
1496 20629 : fs = mkvec2(_shift(f, 1, p, vT), gen_1);
1497 20629 : N = gel(gen_powu_i(fs,p-1,(void*)&D,_can5_sqr,_can5_mul),1);
1498 20629 : N = ZXX_evalx0(FpXQX_red(N,polcyclo(p,vT),q));
1499 20629 : P = FpX_mul(N,f,q);
1500 20629 : P = RgX_deflate(P, p);
1501 20629 : d = RgX_splitting(N, p);
1502 20629 : V = cgetg(p+1,t_VEC);
1503 20629 : gel(V,1) = ZX_mulu(gel(d,1), p);
1504 103915 : for(i=2; i<= (long)p; i++)
1505 83286 : gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
1506 20629 : (void)delete_var(); return gerepilecopy(av, mkvec2(ZX_sub(f,P),V));
1507 : }
1508 :
1509 : static GEN
1510 20629 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
1511 : {
1512 20629 : ulong p = *(long*)E;
1513 20629 : return gen_ZpX_Dixon(gel(v,2), H, qM, utoipos(p), M, E, _can5_lin, _can5_invl);
1514 : }
1515 :
1516 : GEN
1517 13958 : Flx_Teichmuller(GEN P, ulong p, long n)
1518 : {
1519 24199 : return p==3 ? F3x_frobeniuslift(P,n):
1520 10241 : gen_ZpX_Newton(Flx_to_ZX(P),utoipos(p), n, &p, _can5_iter, _can5_invd);
1521 : }
1522 :
1523 : GEN
1524 35 : polteichmuller(GEN P, ulong p, long n)
1525 : {
1526 35 : pari_sp av = avma;
1527 35 : GEN q = NULL;
1528 35 : if (typ(P)!=t_POL || !RgX_is_FpX(P,&q)) pari_err_TYPE("polteichmuller",P);
1529 35 : if (q && !equaliu(q,p)) pari_err_MODULUS("polteichmuller",q,utoi(p));
1530 35 : if (n <= 0)
1531 0 : pari_err_DOMAIN("polteichmuller", "precision", "<=",gen_0,stoi(n));
1532 63 : return gerepileupto(av, p==2 ? F2x_Teichmuller(RgX_to_F2x(P), n)
1533 28 : : Flx_Teichmuller(RgX_to_Flx(P, p), p, n));
1534 : }
|
