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 4682286 : quadratic_prec_mask(long n)
38 : {
39 4682286 : long a = n, i;
40 4682286 : ulong mask = 0;
41 14143837 : for(i = 1;; i++, mask <<= 1)
42 : {
43 14143837 : mask |= (a&1); a = (a+1)>>1;
44 14143837 : if (a==1) return mask | (1UL << i);
45 : }
46 : }
47 :
48 : /***********************************************************************/
49 : /** **/
50 : /** Zp **/
51 : /** **/
52 : /***********************************************************************/
53 :
54 : static GEN
55 1450569 : Zp_divlift(GEN b, GEN a, GEN x, GEN p, long n)
56 : {
57 1450569 : pari_sp ltop = avma, av;
58 : ulong mask;
59 1450569 : GEN q = p;
60 1450569 : if (n == 1) return gcopy(x);
61 1450030 : mask = quadratic_prec_mask(n);
62 1450030 : av = avma;
63 8857304 : while (mask > 1)
64 : {
65 7407271 : GEN v, q2 = q;
66 7407271 : q = sqri(q);
67 7407271 : if (mask & 1UL) q = diviiexact(q,p);
68 7407270 : mask >>= 1;
69 7407270 : if (mask > 1 || !b)
70 : {
71 7299316 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
72 7299316 : x = Fp_sub(x, Fp_mul(v, x, q), q);
73 : }
74 : else
75 : {
76 107954 : GEN y = Fp_mul(x, b, q), yt = modii(y, q2);
77 107955 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
78 107956 : x = Fp_sub(y, Fp_mul(v, yt, q), q);
79 : }
80 7407273 : 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 1450033 : return gerepileupto(ltop, x);
87 : }
88 :
89 : GEN
90 1342616 : Zp_invlift(GEN a, GEN x, GEN p, long e)
91 1342616 : { return Zp_divlift(NULL, a, x, p, e); }
92 :
93 : GEN
94 1342616 : Zp_inv(GEN a, GEN p, long e)
95 : {
96 1342616 : pari_sp av=avma;
97 : GEN ai;
98 1342616 : if (lgefint(p)==3)
99 : {
100 1342353 : ulong pp = p[2];
101 1342353 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
102 : } else
103 263 : ai = Fp_inv(modii(a, p), p);
104 1342616 : return gerepileupto(av, Zp_invlift(a, ai, p, e));
105 : }
106 :
107 : GEN
108 107953 : Zp_div(GEN b, GEN a, GEN p, long e)
109 : {
110 107953 : pari_sp av=avma;
111 : GEN ai;
112 107953 : if (lgefint(p)==3)
113 : {
114 107890 : ulong pp = p[2];
115 107890 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
116 : } else
117 63 : ai = Fp_inv(modii(a, p), p);
118 107953 : 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 248972 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
134 : {
135 248972 : pari_sp av = avma;
136 : int nis2, pis2;
137 : GEN q, w, n_1;
138 : ulong mask;
139 :
140 248972 : if (e == 1) return icopy(a);
141 177983 : nis2 = equaliu(n, 2); n_1 = nis2? NULL: subiu(n,1);
142 177983 : pis2 = equaliu(p, 2);
143 177983 : mask = quadratic_prec_mask(e);
144 177983 : w = nis2 ? shifti(a,1): Fp_mul(n, Fp_pow(a,n_1,p), p);
145 177983 : w = Fp_inv(w, p);
146 177983 : q = p; /* q = p^e; use e instead of q iff p = 2 */
147 177983 : e = 1;
148 : for(;;)
149 : {
150 275633 : 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 275290 : q = sqri(q); if (mask & 1) q = diviiexact(q, p);
163 275290 : mask >>= 1;
164 275290 : if (lgefint(q) == 3 && lgefint(n) == 3)
165 97370 : {
166 274846 : ulong Q = uel(q,2), N = uel(n,2);
167 274846 : ulong A = umodiu(a, Q);
168 274846 : ulong B = umodiu(b, Q);
169 274846 : ulong W = umodiu(w, Q);
170 :
171 274846 : A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
172 274846 : a = utoi(A);
173 274846 : 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 177983 : if (pis2 && signe(a) < 0) a = addii(a, int2n(e));
195 177983 : 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 1255057 : Zp_sqrt(GEN x, GEN p, long e)
205 : {
206 : pari_sp av;
207 : GEN z;
208 1255057 : if (absequaliu(p,2)) return Z2_sqrt(x,e);
209 1253510 : 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 196190 : Zp_log(GEN a, GEN p, ulong e)
225 : {
226 196190 : pari_sp av = avma;
227 196190 : ulong v, N, Np, trunc, pp = itou_or_0(p);
228 196190 : GEN pe, pv, trunc_mod, num, den, ans = gen_0;
229 :
230 196190 : if (equali1(a)) return ans; /* very frequent! */
231 : /* First make the argument closer to 1 by raising it to the p^(v-1) */
232 194881 : v = pp? ulogint(e, pp): 0; /* v here is v-1 */
233 194881 : pe = powiu(p,e); pv = powiu(p,v);
234 194881 : a = Fp_pow(a, pv, mulii(pe, pv));
235 194881 : e += v;
236 :
237 : /* Where do we truncate the Taylor expansion */
238 194881 : N = e + v; N /= ++v; /* note the ++v */
239 194881 : Np = N;
240 : while(1)
241 0 : {
242 194881 : ulong e = Np;
243 194881 : if (pp) e += ulogint(N, pp) / v;
244 194881 : if (e == N) break;
245 0 : N = e;
246 : }
247 :
248 194881 : num = cgetg(N+1, t_INT);
249 194881 : den = cgetg(N+1, t_INT);
250 194881 : trunc = v << 1;
251 194881 : trunc_mod = powiu(p, trunc);
252 : while(1)
253 785401 : { /* compute f = 1 - a_i*p^((v+1)*2^i); trunc_mod = p^((v+1)*2^(i+1)) */
254 980282 : GEN f = modii(a, trunc_mod);
255 980282 : if (!equali1(f))
256 : {
257 977212 : ulong i, step = 1;
258 : GEN h, hpow;
259 :
260 977212 : f = subui(2, f);
261 977212 : a = mulii(a, f);
262 :
263 : /* compute the Taylor expansion of log(f), over Q for now */
264 10553061 : for (i = 1; i <= N; i++) { gel(num,i) = gen_1; gel(den,i) = utoipos(i); }
265 977212 : hpow = h = subui(1, f); /* h = a_i*p^(2^i) */
266 : while(1)
267 : {
268 11381210 : for (i = 1; i <= N - step; i += step << 1)
269 : {
270 8598637 : GEN t = mulii(mulii(hpow, gel(num,i+step)), gel(den,i));
271 8598637 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
272 8598637 : gel(num,i) = addii(gel(num,i), t);
273 8598637 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
274 : }
275 2782573 : step <<= 1; if (step >= N) break;
276 1805361 : hpow = sqri(hpow);
277 : }
278 :
279 977212 : if (pp)
280 : { /* simplify the fraction */
281 976956 : GEN d = powuu(pp, factorial_lval(N, pp));
282 976956 : gel(num,1) = diviiexact(gel(num,1), d);
283 976956 : gel(den,1) = diviiexact(gel(den,1), d);
284 : }
285 :
286 977212 : h = diviiexact(h, pv);
287 977212 : ans = addii(ans, mulii(mulii(gel(num,1), h), Zp_inv(gel(den,1), p, e)));
288 : }
289 980282 : if (trunc > e) break;
290 785401 : trunc_mod = sqri(trunc_mod); trunc <<= 1; N >>= 1;
291 : }
292 194881 : 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 107955 : Zp_exp(GEN a, GEN p, ulong e)
301 : {
302 107955 : pari_sp av = avma;
303 107955 : ulong trunc, N = e, pp = itou_or_0(p);
304 107955 : GEN num, den, trunc_mod = NULL, denominator = gen_1, ans = gen_1;
305 107955 : GEN pe = powiu(p, e);
306 107954 : int pis2 = pp == 2;
307 :
308 : /* Where do we truncate the Taylor expansion */
309 107954 : if (!pis2) N += sdivsi(e, subis(p,2));
310 107954 : num = cgetg(N+2, t_VEC);
311 107955 : den = cgetg(N+2, t_VEC);
312 107955 : if (pis2) trunc = 4;
313 : else
314 : {
315 101317 : trunc = 2;
316 101317 : trunc_mod = sqri(p);
317 : }
318 : while(1)
319 200860 : {
320 308814 : GEN h, hpow, f = pis2? remi2n(a, trunc): modii(a, trunc_mod);
321 308813 : a = subii(a, f);
322 308813 : if (signe(f))
323 : {
324 232883 : ulong step = 1, i;
325 : /* Taylor expansion of exp(f), over Q for now */
326 232883 : gel(num,1) = gen_1;
327 232883 : gel(den,1) = gen_1;
328 1065928 : for (i = 2; i <= N+1; i++)
329 : {
330 833056 : gel(num,i) = gen_1;
331 833056 : gel(den,i) = utoipos(i-1);
332 : }
333 232872 : hpow = h = f;
334 : while(1)
335 : {
336 1379277 : for (i = 1; i <= N+1 - step; i += step << 1) {
337 832855 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
338 832863 : gel(num,i) = addii(gel(num,i), mulii(hpow, gel(num,i+step)));
339 832852 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
340 : }
341 546422 : step <<= 1; if (step > N) break;
342 313539 : hpow = sqri(hpow);
343 : }
344 :
345 : /* We simplify the fraction */
346 232883 : if (pp)
347 : {
348 232757 : GEN d = powuu(pp, factorial_lval(N, pp));
349 232757 : gel(num,1) = diviiexact(gel(num,1), d);
350 232757 : gel(den,1) = diviiexact(gel(den,1), d);
351 : }
352 : /* We add this contribution to exp(f) */
353 232881 : ans = Fp_mul(ans, gel(num,1), pe);
354 232880 : denominator = Fp_mul(denominator, gel(den,1), pe);
355 : }
356 :
357 308813 : if (trunc > e) break;
358 200860 : if (!pis2) trunc_mod = sqri(trunc_mod);
359 200860 : trunc <<= 1; N >>= 1;
360 : }
361 107953 : 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 197493 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
384 : {
385 197493 : long k = lg(a)-1;
386 : long i, j, s, minp, mind;
387 :
388 837864 : for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
389 :
390 442878 : for (j=1; j <= 2*k-5; j+=2,i++)
391 : {
392 245384 : minp = j;
393 245384 : mind = degpol(gel(V,j));
394 1964171 : for (s=j+1; s<i; s++)
395 1718785 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
396 :
397 245386 : swap(gel(V,j), gel(V,minp));
398 245386 : lswap(link[j], link[minp]);
399 :
400 245386 : minp = j+1;
401 245386 : mind = degpol(gel(V,j+1));
402 1718790 : for (s=j+2; s<i; s++)
403 1473404 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
404 :
405 245386 : swap(gel(V,j+1), gel(V,minp));
406 245386 : lswap(link[j+1], link[minp]);
407 :
408 245386 : gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
409 245385 : link[i] = j;
410 : }
411 :
412 640353 : for (j=1; j <= 2*k-3; j+=2)
413 : {
414 : GEN d, u, v;
415 442869 : d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
416 442880 : if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
417 442872 : d = gel(d,2);
418 442872 : if (!gequal1(d))
419 : {
420 398465 : if (typ(d)==t_POL)
421 : {
422 133136 : d = FpXQ_inv(d, T, p);
423 133132 : u = FqX_Fq_mul(u, d, T, p);
424 133128 : v = FqX_Fq_mul(v, d, T, p);
425 : }
426 : else
427 : {
428 265329 : d = Fp_inv(d, p);
429 265327 : u = FqX_Fp_mul(u, d, T,p);
430 265327 : v = FqX_Fp_mul(v, d, T,p);
431 : }
432 : }
433 442859 : gel(W,j) = u;
434 442859 : gel(W,j+1) = v;
435 : }
436 197484 : }
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 1587678 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
442 : {
443 1587678 : pari_sp av = avma;
444 1587678 : long space = lg(f) * lgefint(p1);
445 : GEN a2, b2, g, z, s, t;
446 1587678 : GEN a = gel(V,j), b = gel(V,j+1);
447 1587678 : GEN u = gel(W,j), v = gel(W,j+1);
448 :
449 1587678 : (void)new_chunk(space); /* HACK */
450 1587714 : g = ZX_sub(f, ZX_mul(a,b));
451 1587677 : g = ZX_Z_divexact(g, p0);
452 1587574 : g = FpX_red(g, pd);
453 1587611 : z = FpX_mul(v,g, pd);
454 1587698 : t = FpX_divrem(z,a, pd, &s);
455 1587721 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
456 1587632 : t = FpX_red(t, pd);
457 1587670 : t = ZX_Z_mul(t,p0);
458 1587601 : s = ZX_Z_mul(s,p0);
459 1587586 : set_avma(av);
460 1587574 : a2 = ZX_add(a,s);
461 1587679 : b2 = ZX_add(b,t);
462 :
463 : /* already reduced mod p1 = pd p0 */
464 1587698 : gel(V,j) = a2;
465 1587698 : gel(V,j+1) = b2;
466 1587698 : if (noinv) return;
467 :
468 1281515 : av = avma;
469 1281515 : (void)new_chunk(space); /* HACK */
470 1281532 : g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
471 1281486 : g = Z_ZX_sub(gen_1, g);
472 1281539 : g = ZX_Z_divexact(g, p0);
473 1281445 : g = FpX_red(g, pd);
474 1281417 : z = FpX_mul(v,g, pd);
475 1281548 : t = FpX_divrem(z,a, pd, &s);
476 1281554 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
477 1281490 : t = FpX_red(t, pd);
478 1281520 : t = ZX_Z_mul(t,p0);
479 1281470 : s = ZX_Z_mul(s,p0);
480 1281466 : set_avma(av);
481 1281460 : gel(W,j) = ZX_add(u,t);
482 1281483 : gel(W,j+1) = ZX_add(v,s);
483 : }
484 :
485 : static void
486 503419 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
487 : {
488 503419 : pari_sp av = avma;
489 503419 : const long n = degpol(T1), vT = varn(T1);
490 503388 : long space = lg(f) * lgefint(p1) * lg(T1);
491 : GEN a2, b2, g, z, s, t;
492 503388 : GEN a = gel(V,j), b = gel(V,j+1);
493 503388 : GEN u = gel(W,j), v = gel(W,j+1);
494 :
495 503388 : (void)new_chunk(space); /* HACK */
496 503424 : g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
497 503374 : g = FpXQX_red(g, T1, p1);
498 503264 : g = RgX_Rg_divexact(g, p0);
499 503216 : z = FpXQX_mul(v,g, Td,pd);
500 503419 : t = FpXQX_divrem(z,a, Td,pd, &s);
501 503428 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
502 503356 : t = Kronecker_to_ZXX(t, n, vT);
503 503380 : t = FpXQX_red(t, Td, pd);
504 503266 : t = RgX_Rg_mul(t,p0);
505 503253 : s = RgX_Rg_mul(s,p0);
506 503298 : set_avma(av);
507 :
508 503294 : a2 = RgX_add(a,s);
509 503366 : b2 = RgX_add(b,t);
510 : /* already reduced mod p1 = pd p0 */
511 503393 : gel(V,j) = a2;
512 503393 : gel(V,j+1) = b2;
513 503393 : if (noinv) return;
514 :
515 369615 : av = avma;
516 369615 : (void)new_chunk(space); /* HACK */
517 369608 : g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
518 369560 : g = Kronecker_to_ZXX(g, n, vT);
519 369633 : g = Rg_RgX_sub(gen_1, g);
520 369641 : g = FpXQX_red(g, T1, p1);
521 369578 : g = RgX_Rg_divexact(g, p0);
522 369556 : z = FpXQX_mul(v,g, Td,pd);
523 369641 : t = FpXQX_divrem(z,a, Td,pd, &s);
524 369665 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
525 369617 : t = Kronecker_to_ZXX(t, n, vT);
526 369626 : t = FpXQX_red(t, Td, pd);
527 369576 : t = RgX_Rg_mul(t,p0);
528 369599 : s = RgX_Rg_mul(s,p0);
529 369599 : set_avma(av);
530 369597 : gel(W,j) = RgX_add(u,t);
531 369622 : 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 3916520 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
538 : GEN f, long j, int noinv)
539 : {
540 3916520 : if (j < 0) return;
541 1587488 : ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
542 1587636 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j) , link[j ], noinv);
543 1587716 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
544 : }
545 : static void
546 1139179 : 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 1139179 : if (j < 0) return;
550 503199 : ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
551 503333 : ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j) , link[j ], noinv);
552 503389 : 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 197521 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
562 : {
563 197521 : 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 197521 : if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
569 197521 : if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
570 197521 : if (e0 == 1) return a;
571 :
572 197493 : if (DEBUGLEVEL > 3) timer_start(&Ti);
573 197493 : 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 197493 : e = 1;
583 197493 : v = cgetg(2*k-2 + 1, t_VEC);
584 197493 : w = cgetg(2*k-2 + 1, t_VEC);
585 197493 : link=cgetg(2*k-2 + 1, t_VECSMALL);
586 197493 : BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
587 197486 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
588 : }
589 197486 : mask = quadratic_prec_mask(e0);
590 197486 : eold = 1;
591 197486 : penew = NULL;
592 197486 : Tnew = NULL;
593 197486 : if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
594 1071429 : while (mask > 1)
595 : {
596 873945 : long enew = eold << 1;
597 873945 : if (mask & 1) enew--;
598 873945 : mask >>= 1;
599 873945 : if (enew >= e) { /* mask == 1: last iteration */
600 873945 : GEN peold = penew? penew: powiu(p, eold);
601 873953 : GEN Td = NULL, pd;
602 873953 : long d = enew - eold; /* = eold or eold-1 */
603 : /* lift from p^eold to p^enew */
604 873953 : pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
605 873951 : penew = mulii(peold,pd);
606 873879 : if (T) {
607 132604 : if (Tnew)
608 97914 : Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
609 : else
610 34690 : Td = FpX_red(T, peold);
611 132608 : Tnew = FpX_red(T, penew);
612 132595 : ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
613 : (flag == 0 && mask == 1));
614 : }
615 : else
616 741275 : ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
617 : (flag == 0 && mask == 1));
618 873941 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
619 : }
620 873943 : eold = enew;
621 : }
622 :
623 197484 : if (flag)
624 1750 : E = mkvec4(utoipos(e0), link, v, w);
625 : else
626 : {
627 195734 : E = cgetg(k+1, t_VEC);
628 1075671 : for (i = 1; i <= 2*k-2; i++)
629 : {
630 879936 : long t = link[i];
631 879936 : if (t < 0) gel(E,-t) = gel(v,i);
632 : }
633 : }
634 197486 : 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 244907 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
641 : {
642 244907 : pari_sp av = avma;
643 244907 : pol = FpX_normalize(pol, pe);
644 244906 : if (lg(Q) == 2) return mkvec(pol);
645 161073 : return gerepilecopy(av, MultiLift(pol, Q, NULL, p, e, 0));
646 : }
647 :
648 : GEN
649 34698 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
650 : {
651 34698 : pari_sp av = avma;
652 34698 : 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 46654 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
659 46654 : { 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 7546 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
667 : {
668 : GEN Q, R;
669 7546 : if (j < 0) return;
670 :
671 2898 : Q = FpX_mul(gel(v,j), gel(w,j), pe);
672 2898 : if (U)
673 : {
674 1148 : Q = FpXQ_mul(Q, U, f, pe);
675 1148 : R = FpX_sub(U, Q, pe);
676 : }
677 : else
678 1750 : R = Fp_FpX_sub(gen_1, Q, pe);
679 2898 : gel(w,j+1) = Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
680 2898 : gel(w,j) = R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
681 2898 : BezoutPropagate(link, v, w, pe, R, f, link[j ]);
682 2898 : 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 1771 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
690 : {
691 1771 : pari_sp av = avma;
692 : GEN E, link, v, w, pe;
693 1771 : long i, k = lg(Q)-1;
694 1771 : if (k == 1) retmkvec(pol_1(varn(pol)));
695 1750 : pe = powiu(p, e);
696 1750 : pol = FpX_normalize(pol, pe);
697 1750 : E = MultiLift(pol, Q, NULL, p, e, 1);
698 1750 : link = gel(E,2);
699 1750 : v = gel(E,3);
700 1750 : w = gel(E,4);
701 1750 : BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
702 1750 : E = cgetg(k+1, t_VEC);
703 7546 : for (i = 1; i <= 2*k-2; i++)
704 : {
705 5796 : long t = link[i];
706 5796 : if (t < 0) E[-t] = w[i];
707 : }
708 1750 : 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 96241 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
739 : {
740 96241 : long i,l = lg(x);
741 96241 : GEN r = cgetg(l,t_COL);
742 330797 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
743 96241 : return r;
744 : }
745 :
746 : static GEN
747 96185 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
748 : {
749 96185 : pari_sp av = avma;
750 96185 : GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
751 96185 : return gerepileupto(av, FqV_roots_from_deg1(y, NULL, q));
752 : }
753 :
754 : GEN
755 96046 : ZpX_roots(GEN F, GEN p, long e)
756 : {
757 96046 : pari_sp av = avma;
758 96046 : GEN q = powiu(p, e);
759 96046 : GEN f = FpX_normalize(F, p);
760 96046 : GEN g = FpX_normalize(FpX_split_part(f, p), p);
761 : GEN S;
762 96046 : if (degpol(g) < degpol(f))
763 : {
764 54118 : GEN h = FpX_div(f, g, p);
765 54118 : F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
766 : }
767 96046 : S = FpX_roots(g, p);
768 96046 : 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 6249 : ZqX_roots(GEN F, GEN T, GEN p, long e)
798 : {
799 6249 : 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 6397 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
810 : {
811 6397 : pari_sp av = avma;
812 6397 : GEN q = p, fr, W;
813 : ulong mask;
814 :
815 6397 : a = modii(a,q);
816 6397 : if (e == 1) return a;
817 6397 : mask = quadratic_prec_mask(e);
818 6397 : fr = FpX_red(f,q);
819 6397 : W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
820 : for(;;)
821 : {
822 13644 : q = sqri(q);
823 13644 : if (mask & 1) q = diviiexact(q, p);
824 13644 : mask >>= 1;
825 13644 : fr = FpX_red(f,q);
826 13644 : a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
827 13644 : if (mask == 1) return gerepileuptoint(av, a);
828 7247 : 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 21072 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
903 : {
904 21072 : pari_sp av = avma;
905 21072 : long vT = get_FpX_var(T);
906 : GEN bn, bdi;
907 21072 : GEN bd = ZX_Z_add(x, gen_1);
908 21072 : if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
909 : {
910 7121 : bn = ZX_shifti(x,-(k+1));
911 7121 : 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 21072 : return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
919 : }
920 :
921 : /* Assume p odd, a = 1 [p], return log(a) */
922 : GEN
923 21072 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
924 : {
925 21072 : pari_sp av = avma;
926 : pari_timer ti;
927 21072 : long is2 = absequaliu(p,2);
928 21072 : ulong pp = is2 ? 0: itou_or_0(p);
929 21072 : double lp = is2 ? 1: pp ? log2(pp): expi(p);
930 21072 : long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
931 : GEN ak, s, b, pol;
932 21072 : long e = is2 ? N-1: N;
933 21072 : long i, l = (e-2)/(2*(k+is2));
934 21072 : GEN pe = powiu(p,e);
935 21072 : GEN TNk, pNk = powiu(p,N+k);
936 21072 : if( DEBUGLEVEL>=3) timer_start(&ti);
937 21072 : TNk = FpX_get_red(get_FpX_mod(T), pNk);
938 21072 : ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
939 21072 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
940 21072 : b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
941 21072 : if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
942 21072 : pol= cgetg(l+3,t_POL);
943 21072 : pol[1] = evalsigne(1)|evalvarn(0);
944 57957 : for(i=0; i<=l; i++)
945 : {
946 : GEN g;
947 36885 : ulong z = 2*i+1;
948 36885 : if (pp)
949 : {
950 26243 : long w = u_lvalrem(z, pp, &z);
951 26243 : g = powuu(pp,2*i*k-w);
952 : }
953 10642 : else g = powiu(p,2*i*k);
954 36884 : gel(pol,i+2) = Fp_divu(g, z,pe);
955 : }
956 21072 : if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
957 21072 : s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T, pe);
958 21072 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
959 21072 : s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
960 21072 : 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 267070 : 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 267070 : pari_sp ltop = avma, av;
1081 267070 : long N = 1, N2, M;
1082 : long mask;
1083 267070 : GEN q = p;
1084 267070 : if (n == 1) return gcopy(x);
1085 263241 : mask = quadratic_prec_mask(n);
1086 263241 : av = avma;
1087 930222 : while (mask > 1)
1088 : {
1089 : GEN qM, q2, v, V;
1090 666982 : N2 = N; N <<= 1;
1091 666982 : q2 = q;
1092 666982 : if (mask&1UL) { /* can never happen when q2 = p */
1093 276567 : N--; M = N2-1;
1094 276567 : qM = diviiexact(q2,p); /* > 1 */
1095 276567 : q = mulii(qM,q2);
1096 : } else {
1097 390415 : M = N2;
1098 390415 : qM = q2;
1099 390415 : q = sqri(q2);
1100 : }
1101 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1102 666981 : mask >>= 1;
1103 666981 : v = eval(E, x, q);
1104 666982 : V = ZX_Z_divexact(gel(v,1), q2);
1105 666975 : x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
1106 666980 : 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 263240 : return gerepileupto(ltop, x);
1113 : }
1114 :
1115 : struct _ZpXQ_inv
1116 : {
1117 : GEN T, a, p ,n;
1118 : };
1119 :
1120 : static GEN
1121 523123 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1122 : {
1123 523123 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1124 523123 : GEN Tq = FpXT_red(d->T, q);
1125 : (void)M;
1126 523123 : return FpXQ_mul(V, gel(v,2), Tq, q);
1127 : }
1128 :
1129 : static GEN
1130 523122 : _inv_eval(void *E, GEN x, GEN q)
1131 : {
1132 523122 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1133 523122 : GEN Tq = FpXT_red(d->T, q);
1134 523124 : GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
1135 523124 : return mkvec2(f, x);
1136 : }
1137 :
1138 : GEN
1139 202050 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
1140 : {
1141 : struct _ZpXQ_inv d;
1142 202050 : d.a = a; d.T = T; d.p = p;
1143 202050 : 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 33457 : 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 33457 : 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 33457 : T = FpX_get_red(T, powiu(p, n));
1264 33458 : if (n == 1) return gcopy(S);
1265 33458 : mask = quadratic_prec_mask(n);
1266 33458 : av = avma;
1267 33458 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1268 33458 : if (DEBUGLEVEL > 3) timer_start(&ti);
1269 33458 : Tq = FpXT_red(T,q);
1270 33457 : Tq2 = FpXT_red(Tq,q2);
1271 33457 : Pq = FpX_red(P,q);
1272 33457 : W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
1273 33458 : Q = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
1274 33458 : r = brent_kung_optpow(degpol(P), 4, 3);
1275 33458 : if (DEBUGLEVEL > 3)
1276 0 : err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
1277 : for (;;)
1278 38631 : {
1279 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
1280 72089 : H = FpXQ_mul(W, Q, Tq2, q2);
1281 72089 : Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
1282 72088 : if (DEBUGLEVEL > 3)
1283 0 : timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
1284 72088 : if (mask==1)
1285 6440 : return gerepileupto(ltop, Sq);
1286 65648 : if (early)
1287 : {
1288 63989 : GEN Se = early(E, Sq, q);
1289 63989 : if (Se) return gerepileupto(ltop, Se);
1290 : }
1291 38630 : qq = sqri(q); N <<= 1;
1292 38630 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1293 38630 : mask >>= 1;
1294 38630 : Pqq = FpX_red(P, qq);
1295 38630 : Tqq = FpXT_red(T, qq);
1296 38631 : Spow = FpXQ_powers(Sq, r, Tqq, qq);
1297 38631 : Q = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
1298 38630 : dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
1299 38630 : Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
1300 38631 : Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
1301 38631 : Wq = FpX_sub(W, Wq, q);
1302 38631 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
1303 38631 : 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 : }
|
