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 4687826 : quadratic_prec_mask(long n)
38 : {
39 4687826 : long a = n, i;
40 4687826 : ulong mask = 0;
41 14413473 : for(i = 1;; i++, mask <<= 1)
42 : {
43 14413473 : mask |= (a&1); a = (a+1)>>1;
44 14413473 : if (a==1) return mask | (1UL << i);
45 : }
46 : }
47 :
48 : /***********************************************************************/
49 : /** **/
50 : /** Zp **/
51 : /** **/
52 : /***********************************************************************/
53 :
54 : static GEN
55 1477598 : Zp_divlift(GEN b, GEN a, GEN x, GEN p, long n)
56 : {
57 1477598 : pari_sp ltop = avma, av;
58 : ulong mask;
59 1477598 : GEN q = p;
60 1477598 : if (n == 1) return gcopy(x);
61 1477059 : mask = quadratic_prec_mask(n);
62 1477057 : av = avma;
63 9192087 : while (mask > 1)
64 : {
65 7715029 : GEN v, q2 = q;
66 7715029 : q = sqri(q);
67 7715023 : if (mask & 1UL) q = diviiexact(q,p);
68 7715022 : mask >>= 1;
69 7715022 : if (mask > 1 || !b)
70 : {
71 7607811 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
72 7607816 : x = Fp_sub(x, Fp_mul(v, x, q), q);
73 : }
74 : else
75 : {
76 107211 : GEN y = Fp_mul(x, b, q), yt = modii(y, q2);
77 107215 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
78 107215 : x = Fp_sub(y, Fp_mul(v, yt, q), q);
79 : }
80 7715033 : 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 1477058 : return gerepileupto(ltop, x);
87 : }
88 :
89 : GEN
90 1370382 : Zp_invlift(GEN a, GEN x, GEN p, long e)
91 1370382 : { return Zp_divlift(NULL, a, x, p, e); }
92 :
93 : GEN
94 1370382 : Zp_inv(GEN a, GEN p, long e)
95 : {
96 1370382 : pari_sp av=avma;
97 : GEN ai;
98 1370382 : if (lgefint(p)==3)
99 : {
100 1370119 : ulong pp = p[2];
101 1370119 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
102 : } else
103 263 : ai = Fp_inv(modii(a, p), p);
104 1370382 : return gerepileupto(av, Zp_invlift(a, ai, p, e));
105 : }
106 :
107 : GEN
108 107214 : Zp_div(GEN b, GEN a, GEN p, long e)
109 : {
110 107214 : pari_sp av=avma;
111 : GEN ai;
112 107214 : if (lgefint(p)==3)
113 : {
114 107151 : ulong pp = p[2];
115 107151 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
116 : } else
117 63 : ai = Fp_inv(modii(a, p), p);
118 107216 : 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 248888 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
134 : {
135 248888 : pari_sp av = avma;
136 : int nis2, pis2;
137 : GEN q, w, n_1;
138 : ulong mask;
139 :
140 248888 : if (e == 1) return icopy(a);
141 177955 : nis2 = equaliu(n, 2); n_1 = nis2? NULL: subiu(n,1);
142 177955 : pis2 = equaliu(p, 2);
143 177955 : mask = quadratic_prec_mask(e);
144 177955 : w = nis2 ? shifti(a,1): Fp_mul(n, Fp_pow(a,n_1,p), p);
145 177955 : w = Fp_inv(w, p);
146 177955 : q = p; /* q = p^e; use e instead of q iff p = 2 */
147 177955 : e = 1;
148 : for(;;)
149 : {
150 275605 : 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 275262 : q = sqri(q); if (mask & 1) q = diviiexact(q, p);
163 275262 : mask >>= 1;
164 275262 : if (lgefint(q) == 3 && lgefint(n) == 3)
165 97370 : {
166 274818 : ulong Q = uel(q,2), N = uel(n,2);
167 274818 : ulong A = umodiu(a, Q);
168 274818 : ulong B = umodiu(b, Q);
169 274818 : ulong W = umodiu(w, Q);
170 :
171 274818 : A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
172 274818 : a = utoi(A);
173 274818 : if (mask == 1) break;
174 97370 : if (nis2)
175 84293 : W = Fl_double(Fl_sub(W, Fl_mul(Fl_sqr(W,Q), A, Q), Q), Q);
176 : else
177 13077 : 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 97370 : w = utoi(W);
180 : }
181 : else
182 : {
183 : /* a -= w (a^n - b) */
184 444 : a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q), b))), q);
185 444 : if (mask == 1) break;
186 : /* w += w - w^2 n a^(n-1)*/
187 147 : if (nis2)
188 15 : w = shifti(subii(w, Fp_mul(Fp_sqr(w,q), a, q)), 1);
189 : else
190 132 : w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
191 : mulii(n, Fp_pow(a,n_1,q)), q));
192 : }
193 : }
194 177955 : if (pis2 && signe(a) < 0) a = addii(a, int2n(e));
195 177955 : 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 1255058 : Zp_sqrt(GEN x, GEN p, long e)
205 : {
206 : pari_sp av;
207 : GEN z;
208 1255058 : if (absequaliu(p,2)) return Z2_sqrt(x,e);
209 1253511 : av = avma; z = Fp_sqrt(Fp_red(x, p), p);
210 1253510 : if (!z) return NULL;
211 1253461 : if (e > 1) z = Zp_sqrtlift(x, z, p, e);
212 1253461 : 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 196464 : Zp_log(GEN a, GEN p, ulong e)
225 : {
226 196464 : pari_sp av = avma;
227 196464 : ulong v, N, Np, trunc, pp = itou_or_0(p);
228 196464 : GEN pe, pv, trunc_mod, num, den, ans = gen_0;
229 :
230 196464 : if (equali1(a)) return ans; /* very frequent! */
231 : /* First make the argument closer to 1 by raising it to the p^(v-1) */
232 195435 : v = pp? ulogint(e, pp): 0; /* v here is v-1 */
233 195435 : pe = powiu(p,e); pv = powiu(p,v);
234 195435 : a = Fp_pow(a, pv, mulii(pe, pv));
235 195435 : e += v;
236 :
237 : /* Where do we truncate the Taylor expansion */
238 195435 : N = e + v; N /= ++v; /* note the ++v */
239 195435 : Np = N;
240 : while(1)
241 0 : {
242 195435 : ulong e = Np;
243 195435 : if (pp) e += ulogint(N, pp) / v;
244 195435 : if (e == N) break;
245 0 : N = e;
246 : }
247 :
248 195435 : num = cgetg(N+1, t_INT);
249 195435 : den = cgetg(N+1, t_INT);
250 195435 : trunc = v << 1;
251 195435 : trunc_mod = powiu(p, trunc);
252 : while(1)
253 824434 : { /* compute f = 1 - a_i*p^((v+1)*2^i); trunc_mod = p^((v+1)*2^(i+1)) */
254 1019869 : GEN f = modii(a, trunc_mod);
255 1019869 : if (!equali1(f))
256 : {
257 1017082 : ulong i, step = 1;
258 : GEN h, hpow;
259 :
260 1017082 : f = subui(2, f);
261 1017082 : a = mulii(a, f);
262 :
263 : /* compute the Taylor expansion of log(f), over Q for now */
264 11324418 : for (i = 1; i <= N; i++) { gel(num,i) = gen_1; gel(den,i) = utoipos(i); }
265 1017082 : hpow = h = subui(1, f); /* h = a_i*p^(2^i) */
266 : while(1)
267 : {
268 12220094 : for (i = 1; i <= N - step; i += step << 1)
269 : {
270 9290254 : GEN t = mulii(mulii(hpow, gel(num,i+step)), gel(den,i));
271 9290254 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
272 9290254 : gel(num,i) = addii(gel(num,i), t);
273 9290254 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
274 : }
275 2929840 : step <<= 1; if (step >= N) break;
276 1912759 : hpow = sqri(hpow);
277 : }
278 :
279 1017081 : if (pp)
280 : { /* simplify the fraction */
281 1016825 : GEN d = powuu(pp, factorial_lval(N, pp));
282 1016825 : gel(num,1) = diviiexact(gel(num,1), d);
283 1016825 : gel(den,1) = diviiexact(gel(den,1), d);
284 : }
285 :
286 1017081 : h = diviiexact(h, pv);
287 1017081 : ans = addii(ans, mulii(mulii(gel(num,1), h), Zp_inv(gel(den,1), p, e)));
288 : }
289 1019869 : if (trunc > e) break;
290 824434 : trunc_mod = sqri(trunc_mod); trunc <<= 1; N >>= 1;
291 : }
292 195435 : 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 107213 : Zp_exp(GEN a, GEN p, ulong e)
301 : {
302 107213 : pari_sp av = avma;
303 107213 : ulong trunc, N = e, pp = itou_or_0(p);
304 107213 : GEN num, den, trunc_mod = NULL, denominator = gen_1, ans = gen_1;
305 107213 : GEN pe = powiu(p, e);
306 107216 : int pis2 = pp == 2;
307 :
308 : /* Where do we truncate the Taylor expansion */
309 107216 : if (!pis2) N += sdivsi(e, subis(p,2));
310 107216 : num = cgetg(N+2, t_VEC);
311 107216 : den = cgetg(N+2, t_VEC);
312 107216 : if (pis2) trunc = 4;
313 : else
314 : {
315 100721 : trunc = 2;
316 100721 : trunc_mod = sqri(p);
317 : }
318 : while(1)
319 198217 : {
320 305430 : GEN h, hpow, f = pis2? remi2n(a, trunc): modii(a, trunc_mod);
321 305428 : a = subii(a, f);
322 305430 : if (signe(f))
323 : {
324 229675 : ulong step = 1, i;
325 : /* Taylor expansion of exp(f), over Q for now */
326 229675 : gel(num,1) = gen_1;
327 229675 : gel(den,1) = gen_1;
328 999498 : for (i = 2; i <= N+1; i++)
329 : {
330 769810 : gel(num,i) = gen_1;
331 769810 : gel(den,i) = utoipos(i-1);
332 : }
333 229688 : hpow = h = f;
334 : while(1)
335 : {
336 1304250 : for (i = 1; i <= N+1 - step; i += step << 1) {
337 769532 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
338 769517 : gel(num,i) = addii(gel(num,i), mulii(hpow, gel(num,i+step)));
339 769483 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
340 : }
341 534718 : step <<= 1; if (step > N) break;
342 305041 : hpow = sqri(hpow);
343 : }
344 :
345 : /* We simplify the fraction */
346 229677 : if (pp)
347 : {
348 229551 : GEN d = powuu(pp, factorial_lval(N, pp));
349 229551 : gel(num,1) = diviiexact(gel(num,1), d);
350 229551 : gel(den,1) = diviiexact(gel(den,1), d);
351 : }
352 : /* We add this contribution to exp(f) */
353 229674 : ans = Fp_mul(ans, gel(num,1), pe);
354 229675 : denominator = Fp_mul(denominator, gel(den,1), pe);
355 : }
356 :
357 305430 : if (trunc > e) break;
358 198216 : if (!pis2) trunc_mod = sqri(trunc_mod);
359 198217 : trunc <<= 1; N >>= 1;
360 : }
361 107214 : 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 190790 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
384 : {
385 190790 : long k = lg(a)-1;
386 : long i, j, s, minp, mind;
387 :
388 816731 : for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
389 :
390 435148 : for (j=1; j <= 2*k-5; j+=2,i++)
391 : {
392 244358 : minp = j;
393 244358 : mind = degpol(gel(V,j));
394 1960124 : for (s=j+1; s<i; s++)
395 1715766 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
396 :
397 244358 : swap(gel(V,j), gel(V,minp));
398 244358 : lswap(link[j], link[minp]);
399 :
400 244358 : minp = j+1;
401 244358 : mind = degpol(gel(V,j+1));
402 1715745 : for (s=j+2; s<i; s++)
403 1471384 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
404 :
405 244361 : swap(gel(V,j+1), gel(V,minp));
406 244361 : lswap(link[j+1], link[minp]);
407 :
408 244361 : gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
409 244358 : link[i] = j;
410 : }
411 :
412 625923 : for (j=1; j <= 2*k-3; j+=2)
413 : {
414 : GEN d, u, v;
415 435143 : d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
416 435155 : if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
417 435148 : d = gel(d,2);
418 435148 : if (!gequal1(d))
419 : {
420 393298 : if (typ(d)==t_POL)
421 : {
422 132829 : d = FpXQ_inv(d, T, p);
423 132828 : u = FqX_Fq_mul(u, d, T, p);
424 132821 : v = FqX_Fq_mul(v, d, T, p);
425 : }
426 : else
427 : {
428 260469 : d = Fp_inv(d, p);
429 260467 : u = FqX_Fp_mul(u, d, T,p);
430 260461 : v = FqX_Fp_mul(v, d, T,p);
431 : }
432 : }
433 435133 : gel(W,j) = u;
434 435133 : gel(W,j+1) = v;
435 : }
436 190780 : }
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 1561988 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
442 : {
443 1561988 : pari_sp av = avma;
444 1561988 : long space = lg(f) * lgefint(p1);
445 : GEN a2, b2, g, z, s, t;
446 1561988 : GEN a = gel(V,j), b = gel(V,j+1);
447 1561988 : GEN u = gel(W,j), v = gel(W,j+1);
448 :
449 1561988 : (void)new_chunk(space); /* HACK */
450 1562024 : g = ZX_sub(f, ZX_mul(a,b));
451 1561998 : g = ZX_Z_divexact(g, p0);
452 1561864 : g = FpX_red(g, pd);
453 1561933 : z = FpX_mul(v,g, pd);
454 1562011 : t = FpX_divrem(z,a, pd, &s);
455 1562036 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
456 1561948 : t = FpX_red(t, pd);
457 1561967 : t = ZX_Z_mul(t,p0);
458 1561922 : s = ZX_Z_mul(s,p0);
459 1561927 : set_avma(av);
460 1561916 : a2 = ZX_add(a,s);
461 1561981 : b2 = ZX_add(b,t);
462 :
463 : /* already reduced mod p1 = pd p0 */
464 1562012 : gel(V,j) = a2;
465 1562012 : gel(V,j+1) = b2;
466 1562012 : if (noinv) return;
467 :
468 1263238 : av = avma;
469 1263238 : (void)new_chunk(space); /* HACK */
470 1263244 : g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
471 1263201 : g = Z_ZX_sub(gen_1, g);
472 1263265 : g = ZX_Z_divexact(g, p0);
473 1263158 : g = FpX_red(g, pd);
474 1263126 : z = FpX_mul(v,g, pd);
475 1263257 : t = FpX_divrem(z,a, pd, &s);
476 1263273 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
477 1263202 : t = FpX_red(t, pd);
478 1263230 : t = ZX_Z_mul(t,p0);
479 1263183 : s = ZX_Z_mul(s,p0);
480 1263170 : set_avma(av);
481 1263155 : gel(W,j) = ZX_add(u,t);
482 1263196 : gel(W,j+1) = ZX_add(v,s);
483 : }
484 :
485 : static void
486 502472 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
487 : {
488 502472 : pari_sp av = avma;
489 502472 : const long n = degpol(T1), vT = varn(T1);
490 502466 : long space = lg(f) * lgefint(p1) * lg(T1);
491 : GEN a2, b2, g, z, s, t;
492 502466 : GEN a = gel(V,j), b = gel(V,j+1);
493 502466 : GEN u = gel(W,j), v = gel(W,j+1);
494 :
495 502466 : (void)new_chunk(space); /* HACK */
496 502474 : g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
497 502405 : g = FpXQX_red(g, T1, p1);
498 502340 : g = RgX_Rg_divexact(g, p0);
499 502302 : z = FpXQX_mul(v,g, Td,pd);
500 502470 : t = FpXQX_divrem(z,a, Td,pd, &s);
501 502468 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
502 502416 : t = Kronecker_to_ZXX(t, n, vT);
503 502430 : t = FpXQX_red(t, Td, pd);
504 502365 : t = RgX_Rg_mul(t,p0);
505 502363 : s = RgX_Rg_mul(s,p0);
506 502374 : set_avma(av);
507 :
508 502367 : a2 = RgX_add(a,s);
509 502410 : b2 = RgX_add(b,t);
510 : /* already reduced mod p1 = pd p0 */
511 502458 : gel(V,j) = a2;
512 502458 : gel(V,j+1) = b2;
513 502458 : if (noinv) return;
514 :
515 368979 : av = avma;
516 368979 : (void)new_chunk(space); /* HACK */
517 368979 : g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
518 368874 : g = Kronecker_to_ZXX(g, n, vT);
519 368964 : g = Rg_RgX_sub(gen_1, g);
520 368998 : g = FpXQX_red(g, T1, p1);
521 368946 : g = RgX_Rg_divexact(g, p0);
522 368884 : z = FpXQX_mul(v,g, Td,pd);
523 369009 : t = FpXQX_divrem(z,a, Td,pd, &s);
524 368996 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
525 368966 : t = Kronecker_to_ZXX(t, n, vT);
526 368980 : t = FpXQX_red(t, Td, pd);
527 368953 : t = RgX_Rg_mul(t,p0);
528 368943 : s = RgX_Rg_mul(s,p0);
529 368977 : set_avma(av);
530 368976 : gel(W,j) = RgX_add(u,t);
531 368972 : 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 3842568 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
538 : GEN f, long j, int noinv)
539 : {
540 3842568 : if (j < 0) return;
541 1561813 : ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
542 1561965 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j) , link[j ], noinv);
543 1562014 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
544 : }
545 : static void
546 1137097 : 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 1137097 : if (j < 0) return;
550 502317 : ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
551 502419 : ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j) , link[j ], noinv);
552 502455 : 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 190818 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
562 : {
563 190818 : 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 190818 : if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
569 190818 : if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
570 190818 : if (e0 == 1) return a;
571 :
572 190790 : if (DEBUGLEVEL > 3) timer_start(&Ti);
573 190790 : 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 190790 : e = 1;
583 190790 : v = cgetg(2*k-2 + 1, t_VEC);
584 190790 : w = cgetg(2*k-2 + 1, t_VEC);
585 190790 : link=cgetg(2*k-2 + 1, t_VECSMALL);
586 190789 : BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
587 190783 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
588 : }
589 190783 : mask = quadratic_prec_mask(e0);
590 190784 : eold = 1;
591 190784 : penew = NULL;
592 190784 : Tnew = NULL;
593 190784 : if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
594 1041870 : while (mask > 1)
595 : {
596 851096 : long enew = eold << 1;
597 851096 : if (mask & 1) enew--;
598 851096 : mask >>= 1;
599 851096 : if (enew >= e) { /* mask == 1: last iteration */
600 851096 : GEN peold = penew? penew: powiu(p, eold);
601 851106 : GEN Td = NULL, pd;
602 851106 : long d = enew - eold; /* = eold or eold-1 */
603 : /* lift from p^eold to p^enew */
604 851106 : pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
605 851100 : penew = mulii(peold,pd);
606 851034 : if (T) {
607 132322 : if (Tnew)
608 97718 : Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
609 : else
610 34604 : Td = FpX_red(T, peold);
611 132330 : Tnew = FpX_red(T, penew);
612 132321 : ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
613 : (flag == 0 && mask == 1));
614 : }
615 : else
616 718712 : ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
617 : (flag == 0 && mask == 1));
618 851087 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
619 : }
620 851086 : eold = enew;
621 : }
622 :
623 190774 : if (flag)
624 1743 : E = mkvec4(utoipos(e0), link, v, w);
625 : else
626 : {
627 189031 : E = cgetg(k+1, t_VEC);
628 1053548 : for (i = 1; i <= 2*k-2; i++)
629 : {
630 864508 : long t = link[i];
631 864508 : if (t < 0) gel(E,-t) = gel(v,i);
632 : }
633 : }
634 190783 : 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 234577 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
641 : {
642 234577 : pari_sp av = avma;
643 234577 : pol = FpX_normalize(pol, pe);
644 234577 : if (lg(Q) == 2) return mkvec(pol);
645 154462 : return gerepilecopy(av, MultiLift(pol, Q, NULL, p, e, 0));
646 : }
647 :
648 : GEN
649 34614 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
650 : {
651 34614 : pari_sp av = avma;
652 34614 : pol = FpXQX_normalize(pol, T, pe);
653 34613 : if (lg(Q) == 2) return mkvec(pol);
654 34613 : return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
655 : }
656 :
657 : GEN
658 46478 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
659 46478 : { return T ? ZpXQX_liftfact(f, a, T, pe, p, e): ZpX_liftfact(f, a, pe, p, e); }
660 : GEN
661 154 : ZqX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
662 154 : { 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 7525 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
667 : {
668 : GEN Q, R;
669 7525 : if (j < 0) return;
670 :
671 2891 : Q = FpX_mul(gel(v,j), gel(w,j), pe);
672 2891 : if (U)
673 : {
674 1148 : Q = FpXQ_mul(Q, U, f, pe);
675 1148 : R = FpX_sub(U, Q, pe);
676 : }
677 : else
678 1743 : R = Fp_FpX_sub(gen_1, Q, pe);
679 2891 : gel(w,j+1) = Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
680 2891 : gel(w,j) = R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
681 2891 : BezoutPropagate(link, v, w, pe, R, f, link[j ]);
682 2891 : 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 1764 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
690 : {
691 1764 : pari_sp av = avma;
692 : GEN E, link, v, w, pe;
693 1764 : long i, k = lg(Q)-1;
694 1764 : if (k == 1) retmkvec(pol_1(varn(pol)));
695 1743 : pe = powiu(p, e);
696 1743 : pol = FpX_normalize(pol, pe);
697 1743 : E = MultiLift(pol, Q, NULL, p, e, 1);
698 1743 : link = gel(E,2);
699 1743 : v = gel(E,3);
700 1743 : w = gel(E,4);
701 1743 : BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
702 1743 : E = cgetg(k+1, t_VEC);
703 7525 : for (i = 1; i <= 2*k-2; i++)
704 : {
705 5782 : long t = link[i];
706 5782 : if (t < 0) E[-t] = w[i];
707 : }
708 1743 : 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 90329 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
739 : {
740 90329 : long i,l = lg(x);
741 90329 : GEN r = cgetg(l,t_COL);
742 316186 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
743 90325 : return r;
744 : }
745 :
746 : static GEN
747 90273 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
748 : {
749 90273 : pari_sp av = avma;
750 90273 : GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
751 90273 : return gerepileupto(av, FqV_roots_from_deg1(y, NULL, q));
752 : }
753 :
754 : GEN
755 90134 : ZpX_roots(GEN F, GEN p, long e)
756 : {
757 90134 : pari_sp av = avma;
758 90134 : GEN q = powiu(p, e);
759 90134 : GEN f = FpX_normalize(F, p);
760 90134 : GEN g = FpX_normalize(FpX_split_part(f, p), p);
761 : GEN S;
762 90134 : if (degpol(g) < degpol(f))
763 : {
764 50338 : GEN h = FpX_div(f, g, p);
765 50338 : F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
766 : }
767 90134 : S = FpX_roots(g, p);
768 90134 : 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 6153 : ZqX_roots(GEN F, GEN T, GEN p, long e)
798 : {
799 6153 : 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 6233 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
810 : {
811 6233 : pari_sp av = avma;
812 6233 : GEN q = p, fr, W;
813 : ulong mask;
814 :
815 6233 : a = modii(a,q);
816 6233 : if (e == 1) return a;
817 6233 : mask = quadratic_prec_mask(e);
818 6233 : fr = FpX_red(f,q);
819 6233 : W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
820 : for(;;)
821 : {
822 13393 : q = sqri(q);
823 13393 : if (mask & 1) q = diviiexact(q, p);
824 13393 : mask >>= 1;
825 13393 : fr = FpX_red(f,q);
826 13393 : a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
827 13393 : if (mask == 1) return gerepileuptoint(av, a);
828 7160 : 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 231 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
846 : {
847 231 : pari_sp av = avma, av2;
848 231 : GEN pv = p, q, W, df, Tq, fr, dfr;
849 : ulong mask;
850 231 : a = Fq_red(a, T, p);
851 231 : if (e <= v+1) return a;
852 231 : df = RgX_deriv(f);
853 231 : if (v) { pv = powiu(p,v); df = ZXX_Z_divexact(df, pv); }
854 231 : mask = quadratic_prec_mask(e-v);
855 231 : Tq = FpXT_red(T, p); dfr = FpXQX_red(df, Tq, p);
856 231 : W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
857 231 : q = p;
858 231 : av2 = avma;
859 : for (;;)
860 161 : {
861 : GEN u, fa, qv, q2v, q2, Tq2;
862 392 : q2 = q; q = sqri(q);
863 392 : if (mask & 1) q = diviiexact(q,p);
864 392 : mask >>= 1;
865 392 : if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
866 392 : else { qv = q; q2v = q2; }
867 392 : Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
868 392 : fr = FpXQX_red(f, Tq, qv);
869 392 : fa = FqX_eval(fr, a, Tq, qv);
870 392 : fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
871 392 : a = Fq_sub(a, Fq_mul(Fq_mul(W,fa,Tq2,q2v),q2, Tq,qv), Tq, qv);
872 392 : 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 231 : 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 21069 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
903 : {
904 21069 : pari_sp av = avma;
905 21069 : long vT = get_FpX_var(T);
906 : GEN bn, bdi;
907 21069 : GEN bd = ZX_Z_add(x, gen_1);
908 21069 : if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
909 : {
910 7118 : bn = ZX_shifti(x,-(k+1));
911 7118 : bdi= ZpXQ_invlift(ZX_shifti(bd ,-1), pol_1(vT), T, p, e);
912 : }
913 : else
914 : {
915 13951 : bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
916 13951 : bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
917 : }
918 21069 : return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
919 : }
920 :
921 : /* Assume p odd, a = 1 [p], return log(a) */
922 : GEN
923 21069 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
924 : {
925 21069 : pari_sp av = avma;
926 : pari_timer ti;
927 21069 : long is2 = absequaliu(p,2);
928 21069 : ulong pp = is2 ? 0: itou_or_0(p);
929 21069 : double lp = is2 ? 1: pp ? log2(pp): expi(p);
930 21069 : long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
931 : GEN ak, s, b, pol;
932 21069 : long e = is2 ? N-1: N;
933 21069 : long i, l = (e-2)/(2*(k+is2));
934 21069 : GEN pe = powiu(p,e);
935 21069 : GEN TNk, pNk = powiu(p,N+k);
936 21069 : if( DEBUGLEVEL>=3) timer_start(&ti);
937 21069 : TNk = FpX_get_red(get_FpX_mod(T), pNk);
938 21069 : ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
939 21069 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
940 21069 : b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
941 21069 : if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
942 21069 : pol= cgetg(l+3,t_POL);
943 21069 : pol[1] = evalsigne(1)|evalvarn(0);
944 57930 : for(i=0; i<=l; i++)
945 : {
946 : GEN g;
947 36861 : ulong z = 2*i+1;
948 36861 : if (pp)
949 : {
950 26243 : long w = u_lvalrem(z, pp, &z);
951 26243 : g = powuu(pp,2*i*k-w);
952 : }
953 10618 : else g = powiu(p,2*i*k);
954 36859 : gel(pol,i+2) = Fp_divu(g, z,pe);
955 : }
956 21069 : if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
957 21069 : s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T, pe);
958 21069 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
959 21069 : s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
960 21069 : 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 : GEN
970 0 : gen_ZpM_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
971 : GEN lin(void *E, GEN F, GEN d, GEN q),
972 : GEN invl(void *E, GEN d))
973 : {
974 0 : pari_sp av = avma;
975 : long N2, M;
976 : GEN VN2, V2, VM, bil;
977 : GEN q2, qM;
978 0 : V = FpM_red(V, q);
979 0 : if (N == 1) return invl(E, V);
980 0 : N2 = (N + 1)>>1; M = N - N2;
981 0 : F = FpM_red(F, q);
982 0 : qM = powiu(p, M);
983 0 : q2 = M == N2? qM: mulii(qM, p);
984 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
985 0 : VN2 = gen_ZpM_Dixon(F, V, q2, p, N2, E, lin, invl);
986 0 : bil = lin(E, F, VN2, q);
987 0 : V2 = ZM_Z_divexact(ZM_sub(V, bil), q2);
988 0 : VM = gen_ZpM_Dixon(F, V2, qM, p, M, E, lin, invl);
989 0 : return gerepileupto(av, FpM_red(ZM_add(VN2, ZM_Z_mul(VM, q2)), q));
990 : }
991 :
992 : GEN
993 199955 : gen_ZpX_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
994 : GEN lin(void *E, GEN F, GEN d, GEN q),
995 : GEN invl(void *E, GEN d))
996 : {
997 199955 : pari_sp av = avma;
998 : long N2, M;
999 : GEN VN2, V2, VM, bil;
1000 : GEN q2, qM;
1001 199955 : V = FpX_red(V, q);
1002 199955 : if (N == 1) return invl(E, V);
1003 45458 : N2 = (N + 1)>>1; M = N - N2;
1004 45458 : F = FpXT_red(F, q);
1005 45458 : qM = powiu(p, M);
1006 45458 : q2 = M == N2? qM: mulii(qM, p);
1007 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
1008 45458 : VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
1009 45458 : bil = lin(E, F, VN2, q);
1010 45458 : V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
1011 45458 : VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
1012 45458 : return gerepileupto(av, FpX_red(ZX_add(VN2, ZX_Z_mul(VM, q2)), q));
1013 : }
1014 :
1015 : GEN
1016 0 : gen_ZpM_Newton(GEN x, GEN p, long n, void *E,
1017 : GEN eval(void *E, GEN f, GEN q),
1018 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1019 : {
1020 0 : pari_sp ltop = avma, av;
1021 0 : long N = 1, N2, M;
1022 : long mask;
1023 0 : GEN q = p;
1024 0 : if (n == 1) return gcopy(x);
1025 0 : mask = quadratic_prec_mask(n);
1026 0 : av = avma;
1027 0 : while (mask > 1)
1028 : {
1029 : GEN qM, q2, v, V;
1030 0 : N2 = N; N <<= 1;
1031 0 : q2 = q;
1032 0 : if (mask&1UL) { /* can never happen when q2 = p */
1033 0 : N--; M = N2-1;
1034 0 : qM = diviiexact(q2,p); /* > 1 */
1035 0 : q = mulii(qM,q2);
1036 : } else {
1037 0 : M = N2;
1038 0 : qM = q2;
1039 0 : q = sqri(q2);
1040 : }
1041 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1042 0 : mask >>= 1;
1043 0 : v = eval(E, x, q);
1044 0 : V = ZM_Z_divexact(gel(v,1), q2);
1045 0 : x = FpM_sub(x, ZM_Z_mul(invd(E, V, v, qM, M), q2), q);
1046 0 : if (gc_needed(av, 1))
1047 : {
1048 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Newton");
1049 0 : gerepileall(av, 2, &x, &q);
1050 : }
1051 : }
1052 0 : return gerepileupto(ltop, x);
1053 : }
1054 :
1055 : static GEN
1056 0 : _ZpM_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1057 : {
1058 : (void)E; (void)M;
1059 0 : return FpM_mul(V, gel(v,2), q);
1060 : }
1061 :
1062 : static GEN
1063 0 : _ZpM_eval(void *E, GEN x, GEN q)
1064 : {
1065 0 : GEN f = RgM_Rg_sub_shallow(FpM_mul(x, FpM_red((GEN) E, q), q), gen_1);
1066 0 : return mkvec2(f, x);
1067 : }
1068 :
1069 : GEN
1070 0 : ZpM_invlift(GEN M, GEN C, GEN p, long n)
1071 : {
1072 0 : return gen_ZpM_Newton(C, p, n, (void *)M, _ZpM_eval, _ZpM_invd);
1073 : }
1074 :
1075 : GEN
1076 267061 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
1077 : GEN eval(void *E, GEN f, GEN q),
1078 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1079 : {
1080 267061 : pari_sp ltop = avma, av;
1081 267061 : long N = 1, N2, M;
1082 : long mask;
1083 267061 : GEN q = p;
1084 267061 : if (n == 1) return gcopy(x);
1085 263232 : mask = quadratic_prec_mask(n);
1086 263232 : av = avma;
1087 930157 : while (mask > 1)
1088 : {
1089 : GEN qM, q2, v, V;
1090 666925 : N2 = N; N <<= 1;
1091 666925 : q2 = q;
1092 666925 : if (mask&1UL) { /* can never happen when q2 = p */
1093 276550 : N--; M = N2-1;
1094 276550 : qM = diviiexact(q2,p); /* > 1 */
1095 276550 : q = mulii(qM,q2);
1096 : } else {
1097 390375 : M = N2;
1098 390375 : qM = q2;
1099 390375 : q = sqri(q2);
1100 : }
1101 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1102 666924 : mask >>= 1;
1103 666924 : v = eval(E, x, q);
1104 666927 : V = ZX_Z_divexact(gel(v,1), q2);
1105 666918 : x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
1106 666923 : if (gc_needed(av, 1))
1107 : {
1108 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
1109 0 : gerepileall(av, 2, &x, &q);
1110 : }
1111 : }
1112 263232 : return gerepileupto(ltop, x);
1113 : }
1114 :
1115 : struct _ZpXQ_inv
1116 : {
1117 : GEN T, a, p ,n;
1118 : };
1119 :
1120 : static GEN
1121 523101 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1122 : {
1123 523101 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1124 523101 : GEN Tq = FpXT_red(d->T, q);
1125 : (void)M;
1126 523103 : return FpXQ_mul(V, gel(v,2), Tq, q);
1127 : }
1128 :
1129 : static GEN
1130 523103 : _inv_eval(void *E, GEN x, GEN q)
1131 : {
1132 523103 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1133 523103 : GEN Tq = FpXT_red(d->T, q);
1134 523100 : GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
1135 523103 : return mkvec2(f, x);
1136 : }
1137 :
1138 : GEN
1139 202047 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
1140 : {
1141 : struct _ZpXQ_inv d;
1142 202047 : d.a = a; d.T = T; d.p = p;
1143 202047 : return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
1144 : }
1145 :
1146 : GEN
1147 180978 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
1148 : {
1149 180978 : pari_sp av=avma;
1150 : GEN ai;
1151 180978 : if (lgefint(p)==3)
1152 : {
1153 180950 : ulong pp = p[2];
1154 180950 : ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
1155 : } else
1156 28 : ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
1157 180978 : return gerepileupto(av, ZpXQ_invlift(a, ai, T, p, e));
1158 : }
1159 :
1160 : GEN
1161 31605 : ZpXQ_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
1162 : {
1163 31605 : return FpXQ_mul(a, ZpXQ_inv(b, T, p, e), T, q);
1164 : }
1165 :
1166 : GEN
1167 146440 : ZpXQX_divrem(GEN x, GEN Sp, GEN T, GEN q, GEN p, long e, GEN *pr)
1168 : {
1169 146440 : pari_sp av = avma;
1170 146440 : GEN S = get_FpXQX_mod(Sp);
1171 146440 : GEN b = leading_coeff(S), bi;
1172 : GEN S2, Q;
1173 146440 : if (typ(b)==t_INT) return FpXQX_divrem(x, Sp, T, q, pr);
1174 146440 : bi = ZpXQ_inv(b, T, p, e);
1175 146440 : S2 = FqX_Fq_mul_to_monic(S, bi, T, q);
1176 146440 : Q = FpXQX_divrem(x, S2, T, q, pr);
1177 146440 : if (pr==ONLY_DIVIDES && !Q) { set_avma(av); return NULL; }
1178 146440 : if (pr==ONLY_REM || pr==ONLY_DIVIDES) return gerepileupto(av, Q);
1179 146440 : Q = FpXQX_FpXQ_mul(Q, bi, T, q);
1180 146440 : return gc_all(av, 2, &Q, pr);
1181 : }
1182 :
1183 : GEN
1184 623 : ZpXQX_digits(GEN x, GEN B, GEN T, GEN q, GEN p, long e)
1185 : {
1186 623 : pari_sp av = avma;
1187 623 : GEN b = leading_coeff(B), bi;
1188 : GEN B2, P, V, W;
1189 : long i, lV;
1190 623 : if (typ(b)==t_INT) return FpXQX_digits(x, B, T, q);
1191 623 : bi = ZpXQ_inv(b, T, p, e);
1192 623 : B2 = FqX_Fq_mul_to_monic(B, bi, T, q);
1193 623 : V = FpXQX_digits(x, B2, T, q);
1194 623 : lV = lg(V)-1;
1195 623 : P = FpXQ_powers(bi, lV-1, T, q);
1196 623 : W = cgetg(lV+1, t_VEC);
1197 11200 : for(i=1; i<=lV; i++)
1198 10577 : gel(W, i) = FpXQX_FpXQ_mul(gel(V,i), gel(P, i), T, q);
1199 623 : return gerepileupto(av, W);
1200 : }
1201 :
1202 : struct _ZpXQ_sqrtn
1203 : {
1204 : GEN T, a, n, ai;
1205 : };
1206 :
1207 : static GEN
1208 1806 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
1209 : {
1210 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1211 1806 : GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
1212 : (void)M;
1213 1806 : return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
1214 : }
1215 :
1216 : static GEN
1217 1806 : _sqrtn_eval(void *E, GEN x, GEN q)
1218 : {
1219 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1220 1806 : GEN Tq = FpX_red(d->T, q);
1221 1806 : GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
1222 1806 : return mkvec2(f, x);
1223 : }
1224 :
1225 : GEN
1226 1155 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1227 : {
1228 : struct _ZpXQ_sqrtn d;
1229 1155 : d.a = a; d.T = T; d.n = n;
1230 1155 : d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
1231 1155 : return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
1232 : }
1233 :
1234 : static GEN
1235 0 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
1236 :
1237 : GEN
1238 0 : Zq_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1239 : {
1240 0 : return T? ZpXQ_sqrtnlift(to_ZX(a,varn(T)), n, to_ZX(x,varn(T)), T, p, e)
1241 0 : : Zp_sqrtnlift(a, n, x, p, e);
1242 : }
1243 :
1244 : GEN
1245 0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
1246 : {
1247 0 : pari_sp av = avma;
1248 0 : GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
1249 0 : if (!z) return NULL;
1250 0 : if (e <= 1) return gerepileupto(av, z);
1251 0 : return gerepileupto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
1252 : }
1253 :
1254 : GEN
1255 33390 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
1256 : GEN early(void *E, GEN x, GEN q))
1257 : {
1258 33390 : pari_sp ltop = avma, av;
1259 : long N, r;
1260 : long mask;
1261 : GEN q2, q, W, Q, Tq2, Tq, Pq;
1262 : pari_timer ti;
1263 33390 : T = FpX_get_red(T, powiu(p, n));
1264 33390 : if (n == 1) return gcopy(S);
1265 33390 : mask = quadratic_prec_mask(n);
1266 33390 : av = avma;
1267 33390 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1268 33388 : if (DEBUGLEVEL > 3) timer_start(&ti);
1269 33388 : Tq = FpXT_red(T,q);
1270 33390 : Tq2 = FpXT_red(Tq,q2);
1271 33390 : Pq = FpX_red(P,q);
1272 33389 : W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
1273 33389 : Q = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
1274 33390 : r = brent_kung_optpow(degpol(P), 4, 3);
1275 33390 : if (DEBUGLEVEL > 3)
1276 0 : err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
1277 : for (;;)
1278 38693 : {
1279 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
1280 72083 : H = FpXQ_mul(W, Q, Tq2, q2);
1281 72082 : Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
1282 72083 : if (DEBUGLEVEL > 3)
1283 0 : timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
1284 72083 : if (mask==1)
1285 6462 : return gerepileupto(ltop, Sq);
1286 65621 : if (early)
1287 : {
1288 63962 : GEN Se = early(E, Sq, q);
1289 63961 : if (Se) return gerepileupto(ltop, Se);
1290 : }
1291 38692 : qq = sqri(q); N <<= 1;
1292 38692 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1293 38692 : mask >>= 1;
1294 38692 : Pqq = FpX_red(P, qq);
1295 38693 : Tqq = FpXT_red(T, qq);
1296 38693 : Spow = FpXQ_powers(Sq, r, Tqq, qq);
1297 38693 : Q = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
1298 38693 : dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
1299 38693 : Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
1300 38692 : Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
1301 38692 : Wq = FpX_sub(W, Wq, q);
1302 38693 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
1303 38693 : if (gc_needed(av, 1))
1304 : {
1305 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_liftroot");
1306 0 : gerepileall(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
1307 : }
1308 : }
1309 : }
1310 :
1311 : GEN
1312 1729 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
1313 : {
1314 1729 : return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
1315 : }
1316 :
1317 : GEN
1318 301 : ZpX_Frobenius(GEN T, GEN p, long e)
1319 : {
1320 301 : return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
1321 : }
1322 :
1323 : GEN
1324 147 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
1325 : {
1326 147 : pari_sp av = avma;
1327 147 : GEN xp = ZpX_Frobenius(T, p, e);
1328 147 : GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
1329 147 : return gerepilecopy(av, gel(z,2));
1330 : }
1331 :
1332 : GEN
1333 0 : ZpXQX_ZpXQXQ_liftroot(GEN P, GEN S, GEN U, GEN T, GEN p, long n)
1334 : {
1335 0 : pari_sp ltop = avma, av;
1336 : long N, r;
1337 : long mask;
1338 : GEN qn, q2, q, W, Q, Tq2, Tq, Pq, Uq, Uq2;
1339 : pari_timer ti;
1340 0 : qn = powiu(p, n);
1341 0 : T = FpX_get_red(T, qn);
1342 0 : U = FpXQX_get_red(U, T, qn);
1343 0 : if (n == 1) return gcopy(S);
1344 0 : mask = quadratic_prec_mask(n);
1345 0 : av = avma;
1346 0 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1347 0 : if (DEBUGLEVEL > 3) timer_start(&ti);
1348 0 : Tq = FpXT_red(T,q);
1349 0 : Uq = FpXQXT_red(U, Tq, q);
1350 0 : Tq2 = FpXT_red(Tq, q2);
1351 0 : Uq2 = FpXQXT_red(U, Tq2, q2);
1352 0 : Pq = FpXQX_red(P, Tq, q);
1353 0 : W = FpXQXQ_inv(FpXQX_FpXQXQ_eval(FpXX_deriv(P,q2), S, Uq2, Tq2, q2), Uq2, Tq2, q2);
1354 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQ_eval(Pq, S, Uq, Tq, q), q2);
1355 0 : r = brent_kung_optpow(degpol(P), 4, 3);
1356 0 : if (DEBUGLEVEL > 3)
1357 0 : err_printf("ZpXQX_ZpXQXQ_liftroot: lifting to prec %ld\n",n);
1358 : for (;;)
1359 0 : {
1360 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq, Uqq;
1361 0 : H = FpXQXQ_mul(W, Q, Uq2, Tq2, q2);
1362 0 : Sq = FpXX_sub(S, ZXX_Z_mul(H, q2), q);
1363 0 : if (DEBUGLEVEL > 3)
1364 0 : timer_printf(&ti,"ZpXQX_ZpXQXQ_liftroot: reaching prec %ld",N);
1365 0 : if (mask==1)
1366 0 : return gerepileupto(ltop, Sq);
1367 0 : qq = sqri(q); N <<= 1;
1368 0 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1369 0 : mask >>= 1;
1370 0 : Tqq = FpXT_red(T, qq);
1371 0 : Uqq = FpXQXT_red(U, Tqq, qq);
1372 0 : Pqq = FpXQX_red(P, Tqq, qq);
1373 0 : Spow = FpXQXQ_powers(Sq, r, Uqq, Tqq, qq);
1374 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQV_eval(Pqq, Spow, Uqq, Tqq, qq), q);
1375 0 : dP = FpXQX_FpXQXQV_eval(FpXX_deriv(Pq, q), FpXQXV_red(Spow, Tq, q), Uq, Tq, q);
1376 0 : Wq = ZXX_Z_divexact(gsub(FpXQXQ_mul(W, dP, Uq, Tq, q), gen_1), q2);
1377 0 : Wq = ZXX_Z_mul(FpXQXQ_mul(W, Wq, Uq2, Tq2, q2), q2);
1378 0 : Wq = FpXX_sub(W, Wq, q);
1379 0 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Uq2 = Uq; Uq = Uqq; Pq = Pqq;
1380 0 : if (gc_needed(av, 1))
1381 : {
1382 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_ZpXQXQ_liftroot");
1383 0 : gerepileall(av, 10, &S, &W, &Q, &Uq2, &Uq, &Tq2, &Tq, &Pq, &q, &q2);
1384 : }
1385 : }
1386 : }
1387 :
1388 : GEN
1389 35 : ZqX_ZqXQ_liftroot(GEN f, GEN a, GEN P, GEN T, GEN p, long e)
1390 35 : { return T ? ZpXQX_ZpXQXQ_liftroot(f, a, P, T , p, e): ZpX_ZpXQ_liftroot(f, a, P, p, e); }
1391 :
1392 : /* Canonical lift of polynomial */
1393 :
1394 11130 : static GEN _can_invl(void *E, GEN V) {(void) E; return V; }
1395 :
1396 3654 : static GEN _can_lin(void *E, GEN F, GEN V, GEN q)
1397 : {
1398 3654 : GEN v = RgX_splitting(V, 3);
1399 : (void) E;
1400 3654 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1401 : }
1402 :
1403 : static GEN
1404 7476 : _can_iter(void *E, GEN f, GEN q)
1405 : {
1406 7476 : GEN h = RgX_splitting(f,3);
1407 7476 : GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
1408 7476 : GEN h12 = ZX_mul(gel(h,1), gel(h,2));
1409 7476 : GEN h13 = ZX_mul(gel(h,1), gel(h,3));
1410 7476 : GEN h23 = ZX_mul(gel(h,2), gel(h,3));
1411 7476 : GEN h1c = ZX_mul(gel(h,1), h1s);
1412 7476 : GEN h3c = ZX_mul(gel(h,3), h3s);
1413 7476 : GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
1414 7476 : GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
1415 : (void) E;
1416 7476 : return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
1417 : }
1418 :
1419 : static GEN
1420 7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
1421 : {
1422 7476 : GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
1423 7476 : GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
1424 7476 : GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
1425 : ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
1426 : (void)E;
1427 7476 : return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
1428 : _can_lin, _can_invl);
1429 : }
1430 :
1431 : static GEN
1432 3717 : F3x_frobeniuslift(GEN P, long n)
1433 3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
1434 :
1435 29778 : static GEN _can5_invl(void *E, GEN V) {(void) E; return V; }
1436 :
1437 9107 : static GEN _can5_lin(void *E, GEN F, GEN V, GEN q)
1438 : {
1439 9107 : ulong p = *(ulong*)E;
1440 9107 : GEN v = RgX_splitting(V, p);
1441 9107 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1442 : }
1443 :
1444 : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
1445 : static GEN
1446 62314 : _shift(GEN P, long n, ulong p, long v)
1447 : {
1448 62314 : long i, l=lg(P);
1449 62314 : GEN r = cgetg(l,t_POL); r[1] = P[1];
1450 484498 : for(i=2;i<l;i++)
1451 : {
1452 422184 : long s = n*(i-2)%p;
1453 422184 : GEN ci = gel(P,i);
1454 422184 : if (typ(ci)==t_INT)
1455 105147 : gel(r,i) = monomial(ci, s, v);
1456 : else
1457 317037 : gel(r,i) = RgX_rotate_shallow(ci, s, p);
1458 : }
1459 62314 : return FpXX_renormalize(r, l);
1460 : }
1461 :
1462 : struct _can_mul
1463 : {
1464 : GEN T, q;
1465 : ulong p;
1466 : };
1467 :
1468 : static GEN
1469 41643 : _can5_mul(void *E, GEN A, GEN B)
1470 : {
1471 41643 : struct _can_mul *d = (struct _can_mul *)E;
1472 41643 : GEN a = gel(A,1), b = gel(B,1);
1473 41643 : long n = itos(gel(A,2));
1474 41643 : GEN bn = _shift(b, n, d->p, get_FpX_var(d->T));
1475 41643 : GEN c = FpXQX_mul(a, bn, d->T, d->q);
1476 41643 : return mkvec2(c, addii(gel(A,2), gel(B,2)));
1477 : }
1478 :
1479 : static GEN
1480 41447 : _can5_sqr(void *E, GEN A)
1481 : {
1482 41447 : return _can5_mul(E,A,A);
1483 : }
1484 :
1485 : static GEN
1486 20671 : _can5_iter(void *E, GEN f, GEN q)
1487 : {
1488 20671 : pari_sp av = avma;
1489 : struct _can_mul D;
1490 20671 : ulong p = *(ulong*)E;
1491 20671 : long i, vT = fetch_var();
1492 : GEN N, P, d, V, fs;
1493 20671 : D.q = q; D.T = ZX_Z_sub(pol_xn(p,vT),gen_1);
1494 20671 : D.p = p;
1495 20671 : fs = mkvec2(_shift(f, 1, p, vT), gen_1);
1496 20671 : N = gel(gen_powu_i(fs,p-1,(void*)&D,_can5_sqr,_can5_mul),1);
1497 20671 : N = ZXX_evalx0(FpXQX_red(N,polcyclo(p,vT),q));
1498 20671 : P = FpX_mul(N,f,q);
1499 20671 : P = RgX_deflate(P, p);
1500 20671 : d = RgX_splitting(N, p);
1501 20671 : V = cgetg(p+1,t_VEC);
1502 20671 : gel(V,1) = ZX_mulu(gel(d,1), p);
1503 104265 : for(i=2; i<= (long)p; i++)
1504 83594 : gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
1505 20671 : (void)delete_var(); return gerepilecopy(av, mkvec2(ZX_sub(f,P),V));
1506 : }
1507 :
1508 : static GEN
1509 20671 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
1510 : {
1511 20671 : ulong p = *(long*)E;
1512 20671 : return gen_ZpX_Dixon(gel(v,2), H, qM, utoipos(p), M, E, _can5_lin, _can5_invl);
1513 : }
1514 :
1515 : GEN
1516 13979 : Flx_Teichmuller(GEN P, ulong p, long n)
1517 : {
1518 24241 : return p==3 ? F3x_frobeniuslift(P,n):
1519 10262 : gen_ZpX_Newton(Flx_to_ZX(P),utoipos(p), n, &p, _can5_iter, _can5_invd);
1520 : }
1521 :
1522 : GEN
1523 35 : polteichmuller(GEN P, ulong p, long n)
1524 : {
1525 35 : pari_sp av = avma;
1526 35 : GEN q = NULL;
1527 35 : if (typ(P)!=t_POL || !RgX_is_FpX(P,&q)) pari_err_TYPE("polteichmuller",P);
1528 35 : if (q && !equaliu(q,p)) pari_err_MODULUS("polteichmuller",q,utoi(p));
1529 35 : if (n <= 0)
1530 0 : pari_err_DOMAIN("polteichmuller", "precision", "<=",gen_0,stoi(n));
1531 63 : return gerepileupto(av, p==2 ? F2x_Teichmuller(RgX_to_F2x(P), n)
1532 28 : : Flx_Teichmuller(RgX_to_Flx(P, p), p, n));
1533 : }
|
