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 4880154 : quadratic_prec_mask(long n)
38 : {
39 4880154 : long a = n, i;
40 4880154 : ulong mask = 0;
41 14582908 : for(i = 1;; i++, mask <<= 1)
42 : {
43 14582908 : mask |= (a&1); a = (a+1)>>1;
44 14582908 : if (a==1) return mask | (1UL << i);
45 : }
46 : }
47 :
48 : /***********************************************************************/
49 : /** **/
50 : /** Zp **/
51 : /** **/
52 : /***********************************************************************/
53 :
54 : static GEN
55 1450706 : Zp_divlift(GEN b, GEN a, GEN x, GEN p, long n)
56 : {
57 1450706 : pari_sp ltop = avma, av;
58 : ulong mask;
59 1450706 : GEN q = p;
60 1450706 : if (n == 1) return gcopy(x);
61 1450167 : mask = quadratic_prec_mask(n);
62 1450167 : av = avma;
63 8858253 : while (mask > 1)
64 : {
65 7408087 : GEN v, q2 = q;
66 7408087 : q = sqri(q);
67 7408086 : if (mask & 1UL) q = diviiexact(q,p);
68 7408085 : mask >>= 1;
69 7408085 : if (mask > 1 || !b)
70 : {
71 7300101 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
72 7300103 : x = Fp_sub(x, Fp_mul(v, x, q), q);
73 : }
74 : else
75 : {
76 107984 : GEN y = Fp_mul(x, b, q), yt = modii(y, q2);
77 107985 : v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
78 107984 : x = Fp_sub(y, Fp_mul(v, yt, q), q);
79 : }
80 7408086 : if (gc_needed(av, 1))
81 : {
82 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_Zp_Newton");
83 0 : (void)gc_all(av, 2, &x, &q);
84 : }
85 : }
86 1450166 : return gc_upto(ltop, x);
87 : }
88 :
89 : GEN
90 1342721 : Zp_invlift(GEN a, GEN x, GEN p, long e)
91 1342721 : { return Zp_divlift(NULL, a, x, p, e); }
92 :
93 : GEN
94 1342721 : Zp_inv(GEN a, GEN p, long e)
95 : {
96 1342721 : pari_sp av=avma;
97 : GEN ai;
98 1342721 : if (lgefint(p)==3)
99 : {
100 1342458 : ulong pp = p[2];
101 1342458 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
102 : } else
103 263 : ai = Fp_inv(modii(a, p), p);
104 1342721 : return gc_upto(av, Zp_invlift(a, ai, p, e));
105 : }
106 :
107 : GEN
108 107987 : Zp_div(GEN b, GEN a, GEN p, long e)
109 : {
110 107987 : pari_sp av=avma;
111 : GEN ai;
112 107987 : if (lgefint(p)==3)
113 : {
114 107924 : ulong pp = p[2];
115 107924 : ai = utoi(Fl_inv(umodiu(a,pp), pp));
116 : } else
117 63 : ai = Fp_inv(modii(a, p), p);
118 107985 : return gc_upto(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 265863 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
134 : {
135 265863 : pari_sp av = avma;
136 : int nis2, pis2;
137 : GEN q, w, n_1;
138 : ulong mask;
139 :
140 265863 : if (e == 1) return icopy(a);
141 194734 : nis2 = equaliu(n, 2); n_1 = nis2? NULL: subiu(n,1);
142 194734 : pis2 = equaliu(p, 2);
143 194734 : mask = quadratic_prec_mask(e);
144 194734 : w = nis2 ? shifti(a,1): Fp_mul(n, Fp_pow(a,n_1,p), p);
145 194734 : w = Fp_inv(w, p);
146 194734 : q = p; /* q = p^e; use e instead of q iff p = 2 */
147 194734 : e = 1;
148 : for(;;)
149 : {
150 296430 : 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 296087 : q = sqri(q); if (mask & 1) q = diviiexact(q, p);
163 296087 : mask >>= 1;
164 296087 : if (lgefint(q) == 3 && lgefint(n) == 3)
165 101414 : {
166 295626 : ulong Q = uel(q,2), N = uel(n,2);
167 295626 : ulong A = umodiu(a, Q);
168 295626 : ulong B = umodiu(b, Q);
169 295626 : ulong W = umodiu(w, Q);
170 :
171 295626 : A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
172 295626 : a = utoi(A);
173 295626 : if (mask == 1) break;
174 101414 : if (nis2)
175 88255 : W = Fl_double(Fl_sub(W, Fl_mul(Fl_sqr(W,Q), A, Q), Q), Q);
176 : else
177 13159 : 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 101414 : w = utoi(W);
180 : }
181 : else
182 : {
183 : /* a -= w (a^n - b) */
184 461 : a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q), b))), q);
185 461 : if (mask == 1) break;
186 : /* w += w - w^2 n a^(n-1)*/
187 149 : if (nis2)
188 15 : w = shifti(subii(w, Fp_mul(Fp_sqr(w,q), a, q)), 1);
189 : else
190 134 : w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
191 : mulii(n, Fp_pow(a,n_1,q)), q));
192 : }
193 : }
194 194734 : if (pis2 && signe(a) < 0) a = addii(a, int2n(e));
195 194734 : return gc_INT(av, a);
196 : }
197 :
198 : /* ZpX_liftroot for the polynomial X^2-b */
199 : GEN
200 156632 : Zp_sqrtlift(GEN b, GEN a, GEN p, long e)
201 156632 : { return Zp_sqrtnlift(b, gen_2, a, p, e); }
202 :
203 : GEN
204 1329879 : Zp_sqrt(GEN x, GEN p, long e)
205 : {
206 : pari_sp av;
207 : GEN z;
208 1329879 : if (absequaliu(p,2)) return Z2_sqrt(x,e);
209 1328241 : av = avma; z = Fp_sqrt(Fp_red(x, p), p);
210 1328243 : if (!z) return NULL;
211 1328201 : if (e > 1) z = Zp_sqrtlift(x, z, p, e);
212 1328201 : return gc_INT(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 196232 : Zp_log(GEN a, GEN p, ulong e)
225 : {
226 196232 : pari_sp av = avma;
227 196232 : ulong v, N, Np, trunc, pp = itou_or_0(p);
228 196232 : GEN pe, pv, trunc_mod, num, den, ans = gen_0;
229 :
230 196232 : if (equali1(a)) return ans; /* very frequent! */
231 : /* First make the argument closer to 1 by raising it to the p^(v-1) */
232 194923 : v = pp? ulogint(e, pp): 0; /* v here is v-1 */
233 194923 : pe = powiu(p,e); pv = powiu(p,v);
234 194923 : a = Fp_pow(a, pv, mulii(pe, pv));
235 194923 : e += v;
236 :
237 : /* Where do we truncate the Taylor expansion */
238 194923 : N = e + v; N /= ++v; /* note the ++v */
239 194923 : Np = N;
240 : while(1)
241 0 : {
242 194923 : ulong e = Np;
243 194923 : if (pp) e += ulogint(N, pp) / v;
244 194923 : if (e == N) break;
245 0 : N = e;
246 : }
247 :
248 194923 : num = cgetg(N+1, t_INT);
249 194923 : den = cgetg(N+1, t_INT);
250 194923 : trunc = v << 1;
251 194923 : trunc_mod = powiu(p, trunc);
252 : while(1)
253 785485 : { /* compute f = 1 - a_i*p^((v+1)*2^i); trunc_mod = p^((v+1)*2^(i+1)) */
254 980408 : GEN f = modii(a, trunc_mod);
255 980408 : if (!equali1(f))
256 : {
257 977317 : ulong i, step = 1;
258 : GEN h, hpow;
259 :
260 977317 : f = subui(2, f);
261 977317 : a = mulii(a, f);
262 :
263 : /* compute the Taylor expansion of log(f), over Q for now */
264 10553481 : for (i = 1; i <= N; i++) { gel(num,i) = gen_1; gel(den,i) = utoipos(i); }
265 977317 : hpow = h = subui(1, f); /* h = a_i*p^(2^i) */
266 : while(1)
267 : {
268 11381609 : for (i = 1; i <= N - step; i += step << 1)
269 : {
270 8598847 : GEN t = mulii(mulii(hpow, gel(num,i+step)), gel(den,i));
271 8598847 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
272 8598847 : gel(num,i) = addii(gel(num,i), t);
273 8598847 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
274 : }
275 2782762 : step <<= 1; if (step >= N) break;
276 1805445 : hpow = sqri(hpow);
277 : }
278 :
279 977317 : if (pp)
280 : { /* simplify the fraction */
281 977061 : GEN d = powuu(pp, factorial_lval(N, pp));
282 977061 : gel(num,1) = diviiexact(gel(num,1), d);
283 977061 : gel(den,1) = diviiexact(gel(den,1), d);
284 : }
285 :
286 977317 : h = diviiexact(h, pv);
287 977317 : ans = addii(ans, mulii(mulii(gel(num,1), h), Zp_inv(gel(den,1), p, e)));
288 : }
289 980408 : if (trunc > e) break;
290 785485 : trunc_mod = sqri(trunc_mod); trunc <<= 1; N >>= 1;
291 : }
292 194923 : return gc_INT(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 107987 : Zp_exp(GEN a, GEN p, ulong e)
301 : {
302 107987 : pari_sp av = avma;
303 107987 : ulong trunc, N = e, pp = itou_or_0(p);
304 107987 : GEN num, den, trunc_mod = NULL, denominator = gen_1, ans = gen_1;
305 107987 : GEN pe = powiu(p, e);
306 107987 : int pis2 = pp == 2;
307 :
308 : /* Where do we truncate the Taylor expansion */
309 107987 : if (!pis2) N += sdivsi(e, subis(p,2));
310 107987 : num = cgetg(N+2, t_VEC);
311 107987 : den = cgetg(N+2, t_VEC);
312 107987 : if (pis2) trunc = 4;
313 : else
314 : {
315 101332 : trunc = 2;
316 101332 : trunc_mod = sqri(p);
317 : }
318 : while(1)
319 200995 : {
320 308982 : GEN h, hpow, f = pis2? remi2n(a, trunc): modii(a, trunc_mod);
321 308981 : a = subii(a, f);
322 308980 : if (signe(f))
323 : {
324 233043 : ulong step = 1, i;
325 : /* Taylor expansion of exp(f), over Q for now */
326 233043 : gel(num,1) = gen_1;
327 233043 : gel(den,1) = gen_1;
328 1068550 : for (i = 2; i <= N+1; i++)
329 : {
330 835502 : gel(num,i) = gen_1;
331 835502 : gel(den,i) = utoipos(i-1);
332 : }
333 233048 : hpow = h = f;
334 : while(1)
335 : {
336 1382373 : for (i = 1; i <= N+1 - step; i += step << 1) {
337 835341 : gel(num,i) = mulii(gel(num,i), gel(den,i+step));
338 835336 : gel(num,i) = addii(gel(num,i), mulii(hpow, gel(num,i+step)));
339 835314 : gel(den,i) = mulii(gel(den,i), gel(den,i+step));
340 : }
341 547032 : step <<= 1; if (step > N) break;
342 313989 : hpow = sqri(hpow);
343 : }
344 :
345 : /* We simplify the fraction */
346 233043 : if (pp)
347 : {
348 232917 : GEN d = powuu(pp, factorial_lval(N, pp));
349 232917 : gel(num,1) = diviiexact(gel(num,1), d);
350 232916 : gel(den,1) = diviiexact(gel(den,1), d);
351 : }
352 : /* We add this contribution to exp(f) */
353 233042 : ans = Fp_mul(ans, gel(num,1), pe);
354 233043 : denominator = Fp_mul(denominator, gel(den,1), pe);
355 : }
356 :
357 308982 : if (trunc > e) break;
358 200995 : if (!pis2) trunc_mod = sqri(trunc_mod);
359 200995 : trunc <<= 1; N >>= 1;
360 : }
361 107987 : return gc_INT(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 201480 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
384 : {
385 201480 : long k = lg(a)-1;
386 : long i, j, s, minp, mind;
387 :
388 852467 : for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
389 :
390 449510 : for (j=1; j <= 2*k-5; j+=2,i++)
391 : {
392 248031 : minp = j;
393 248031 : mind = degpol(gel(V,j));
394 2000410 : for (s=j+1; s<i; s++)
395 1752380 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
396 :
397 248030 : swap(gel(V,j), gel(V,minp));
398 248030 : lswap(link[j], link[minp]);
399 :
400 248030 : minp = j+1;
401 248030 : mind = degpol(gel(V,j+1));
402 1752412 : for (s=j+2; s<i; s++)
403 1504381 : if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
404 :
405 248031 : swap(gel(V,j+1), gel(V,minp));
406 248031 : lswap(link[j+1], link[minp]);
407 :
408 248031 : gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
409 248030 : link[i] = j;
410 : }
411 :
412 650984 : for (j=1; j <= 2*k-3; j+=2)
413 : {
414 : GEN d, u, v;
415 449516 : d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
416 449519 : if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
417 449511 : d = gel(d,2);
418 449511 : if (!gequal1(d))
419 : {
420 404300 : if (typ(d)==t_POL)
421 : {
422 133469 : d = FpXQ_inv(d, T, p);
423 133461 : u = FqX_Fq_mul(u, d, T, p);
424 133464 : v = FqX_Fq_mul(v, d, T, p);
425 : }
426 : else
427 : {
428 270831 : d = Fp_inv(d, p);
429 270828 : u = FqX_Fp_mul(u, d, T,p);
430 270829 : v = FqX_Fp_mul(v, d, T,p);
431 : }
432 : }
433 449505 : gel(W,j) = u;
434 449505 : gel(W,j+1) = v;
435 : }
436 201468 : }
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 1615364 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
442 : {
443 1615364 : pari_sp av = avma;
444 1615364 : long space = lg(f) * lgefint(p1);
445 : GEN a2, b2, g, z, s, t;
446 1615364 : GEN a = gel(V,j), b = gel(V,j+1);
447 1615364 : GEN u = gel(W,j), v = gel(W,j+1);
448 :
449 1615364 : (void)new_chunk(space); /* HACK */
450 1615404 : g = ZX_sub(f, ZX_mul(a,b));
451 1615353 : g = ZX_Z_divexact(g, p0);
452 1615236 : g = FpX_red(g, pd);
453 1615274 : z = FpX_mul(v,g, pd);
454 1615398 : t = FpX_divrem(z,a, pd, &s);
455 1615405 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
456 1615319 : t = FpX_red(t, pd);
457 1615393 : t = ZX_Z_mul(t,p0);
458 1615308 : s = ZX_Z_mul(s,p0);
459 1615331 : set_avma(av);
460 1615328 : a2 = ZX_add(a,s);
461 1615355 : b2 = ZX_add(b,t);
462 :
463 : /* already reduced mod p1 = pd p0 */
464 1615379 : gel(V,j) = a2;
465 1615379 : gel(V,j+1) = b2;
466 1615379 : if (noinv) return;
467 :
468 1302939 : av = avma;
469 1302939 : (void)new_chunk(space); /* HACK */
470 1302964 : g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
471 1302912 : g = Z_ZX_sub(gen_1, g);
472 1302984 : g = ZX_Z_divexact(g, p0);
473 1302862 : g = FpX_red(g, pd);
474 1302851 : z = FpX_mul(v,g, pd);
475 1302973 : t = FpX_divrem(z,a, pd, &s);
476 1302991 : t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
477 1302927 : t = FpX_red(t, pd);
478 1302988 : t = ZX_Z_mul(t,p0);
479 1302920 : s = ZX_Z_mul(s,p0);
480 1302919 : set_avma(av);
481 1302915 : gel(W,j) = ZX_add(u,t);
482 1302934 : gel(W,j+1) = ZX_add(v,s);
483 : }
484 :
485 : static void
486 507236 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
487 : {
488 507236 : pari_sp av = avma;
489 507236 : const long n = degpol(T1), vT = varn(T1);
490 507231 : long space = lg(f) * lgefint(p1) * lg(T1);
491 : GEN a2, b2, g, z, s, t;
492 507231 : GEN a = gel(V,j), b = gel(V,j+1);
493 507231 : GEN u = gel(W,j), v = gel(W,j+1);
494 :
495 507231 : (void)new_chunk(space); /* HACK */
496 507237 : g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
497 507176 : g = FpXQX_red(g, T1, p1);
498 507082 : g = RgX_Rg_divexact(g, p0);
499 507026 : z = FpXQX_mul(v,g, Td,pd);
500 507220 : t = FpXQX_divrem(z,a, Td,pd, &s);
501 507241 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
502 507129 : t = Kronecker_to_ZXX(t, n, vT);
503 507170 : t = FpXQX_red(t, Td, pd);
504 507120 : t = RgX_Rg_mul(t,p0);
505 507103 : s = RgX_Rg_mul(s,p0);
506 507122 : set_avma(av);
507 :
508 507122 : a2 = RgX_add(a,s);
509 507162 : b2 = RgX_add(b,t);
510 : /* already reduced mod p1 = pd p0 */
511 507213 : gel(V,j) = a2;
512 507213 : gel(V,j+1) = b2;
513 507213 : if (noinv) return;
514 :
515 373012 : av = avma;
516 373012 : (void)new_chunk(space); /* HACK */
517 373011 : g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
518 372936 : g = Kronecker_to_ZXX(g, n, vT);
519 373012 : g = Rg_RgX_sub(gen_1, g);
520 373028 : g = FpXQX_red(g, T1, p1);
521 372962 : g = RgX_Rg_divexact(g, p0);
522 372926 : z = FpXQX_mul(v,g, Td,pd);
523 373040 : t = FpXQX_divrem(z,a, Td,pd, &s);
524 373052 : t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
525 372982 : t = Kronecker_to_ZXX(t, n, vT);
526 372999 : t = FpXQX_red(t, Td, pd);
527 372974 : t = RgX_Rg_mul(t,p0);
528 372954 : s = RgX_Rg_mul(s,p0);
529 372960 : set_avma(av);
530 372959 : gel(W,j) = RgX_add(u,t);
531 372977 : 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 3988052 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
538 : GEN f, long j, int noinv)
539 : {
540 3988052 : if (j < 0) return;
541 1615159 : ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
542 1615329 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j) , link[j ], noinv);
543 1615387 : ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
544 : }
545 : static void
546 1146933 : 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 1146933 : if (j < 0) return;
550 507055 : ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
551 507160 : ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j) , link[j ], noinv);
552 507204 : 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 201525 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
562 : {
563 201525 : 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 201525 : if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
569 201525 : if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
570 201525 : if (e0 == 1) return a;
571 :
572 201483 : if (DEBUGLEVEL > 3) timer_start(&Ti);
573 201483 : 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 201483 : e = 1;
583 201483 : v = cgetg(2*k-2 + 1, t_VEC);
584 201483 : w = cgetg(2*k-2 + 1, t_VEC);
585 201483 : link=cgetg(2*k-2 + 1, t_VECSMALL);
586 201483 : BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
587 201477 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
588 : }
589 201477 : mask = quadratic_prec_mask(e0);
590 201477 : eold = 1;
591 201477 : penew = NULL;
592 201477 : Tnew = NULL;
593 201477 : if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
594 1091679 : while (mask > 1)
595 : {
596 890209 : long enew = eold << 1;
597 890209 : if (mask & 1) enew--;
598 890209 : mask >>= 1;
599 890209 : if (enew >= e) { /* mask == 1: last iteration */
600 890209 : GEN peold = penew? penew: powiu(p, eold);
601 890217 : GEN Td = NULL, pd;
602 890217 : long d = enew - eold; /* = eold or eold-1 */
603 : /* lift from p^eold to p^enew */
604 890217 : pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
605 890214 : penew = mulii(peold,pd);
606 890141 : if (T) {
607 132670 : if (Tnew)
608 97972 : Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
609 : else
610 34698 : Td = FpX_red(T, peold);
611 132678 : Tnew = FpX_red(T, penew);
612 132672 : ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
613 : (flag == 0 && mask == 1));
614 : }
615 : else
616 757471 : ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
617 : (flag == 0 && mask == 1));
618 890202 : if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
619 : }
620 890202 : eold = enew;
621 : }
622 :
623 201470 : if (flag)
624 1715 : E = mkvec4(utoipos(e0), link, v, w);
625 : else
626 : {
627 199755 : E = cgetg(k+1, t_VEC);
628 1093036 : for (i = 1; i <= 2*k-2; i++)
629 : {
630 893278 : long t = link[i];
631 893278 : if (t < 0) gel(E,-t) = gel(v,i);
632 : }
633 : }
634 201473 : 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 249595 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
641 : {
642 249595 : pari_sp av = avma;
643 249595 : pol = FpX_normalize(pol, pe);
644 249595 : if (lg(Q) == 2) return mkvec(pol);
645 165105 : return gc_GEN(av, MultiLift(pol, Q, NULL, p, e, 0));
646 : }
647 :
648 : GEN
649 34705 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
650 : {
651 34705 : pari_sp av = avma;
652 34705 : pol = FpXQX_normalize(pol, T, pe);
653 34705 : if (lg(Q) == 2) return mkvec(pol);
654 34705 : return gc_GEN(av, MultiLift(pol, Q, T, p, e, 0));
655 : }
656 :
657 : GEN
658 46570 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
659 46570 : { return T ? ZpXQX_liftfact(f, a, T, pe, p, e): ZpX_liftfact(f, a, pe, p, e); }
660 : GEN
661 112 : ZqX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
662 112 : { 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 7441 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
667 : {
668 : GEN Q, R;
669 7441 : if (j < 0) return;
670 :
671 2863 : Q = FpX_mul(gel(v,j), gel(w,j), pe);
672 2863 : if (U)
673 : {
674 1148 : Q = FpXQ_mul(Q, U, f, pe);
675 1148 : R = FpX_sub(U, Q, pe);
676 : }
677 : else
678 1715 : R = Fp_FpX_sub(gen_1, Q, pe);
679 2863 : gel(w,j+1) = Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
680 2863 : gel(w,j) = R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
681 2863 : BezoutPropagate(link, v, w, pe, R, f, link[j ]);
682 2863 : 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 1736 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
690 : {
691 1736 : pari_sp av = avma;
692 : GEN E, link, v, w, pe;
693 1736 : long i, k = lg(Q)-1;
694 1736 : if (k == 1) retmkvec(pol_1(varn(pol)));
695 1715 : pe = powiu(p, e);
696 1715 : pol = FpX_normalize(pol, pe);
697 1715 : E = MultiLift(pol, Q, NULL, p, e, 1);
698 1715 : link = gel(E,2);
699 1715 : v = gel(E,3);
700 1715 : w = gel(E,4);
701 1715 : BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
702 1715 : E = cgetg(k+1, t_VEC);
703 7441 : for (i = 1; i <= 2*k-2; i++)
704 : {
705 5726 : long t = link[i];
706 5726 : if (t < 0) E[-t] = w[i];
707 : }
708 1715 : return gc_GEN(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 gc_GEN(av, ZqX_liftfact(pol, L, T, powiu(p,N), p, N));
735 : }
736 :
737 : static GEN
738 98138 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
739 : {
740 98138 : long i,l = lg(x);
741 98138 : GEN r = cgetg(l,t_COL);
742 337788 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
743 98134 : return r;
744 : }
745 :
746 : static GEN
747 98076 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
748 : {
749 98076 : pari_sp av = avma;
750 98076 : GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
751 98075 : return gc_upto(av, FqV_roots_from_deg1(y, NULL, q));
752 : }
753 :
754 : GEN
755 97937 : ZpX_roots(GEN F, GEN p, long e)
756 : {
757 97937 : pari_sp av = avma;
758 97937 : GEN q = powiu(p, e);
759 97937 : GEN f = FpX_normalize(F, p);
760 97937 : GEN g = FpX_normalize(FpX_split_part(f, p), p);
761 : GEN S;
762 97936 : if (degpol(g) < degpol(f))
763 : {
764 54830 : GEN h = FpX_div(f, g, p);
765 54830 : F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
766 : }
767 97936 : S = FpX_roots(g, p);
768 97937 : return gc_upto(av, ZpX_liftroots_full(F, S, q, p, e));
769 : }
770 :
771 : static GEN
772 63 : ZpXQX_liftroots_full(GEN f, GEN S, GEN T, GEN q, GEN p, long e)
773 : {
774 63 : pari_sp av = avma;
775 63 : GEN y = ZpXQX_liftfact(f, deg1_from_roots(S, varn(f)), T, q, p, e);
776 63 : return gc_upto(av, FqV_roots_from_deg1(y, T, q));
777 : }
778 :
779 : GEN
780 63 : ZpXQX_roots(GEN F, GEN T, GEN p, long e)
781 : {
782 63 : pari_sp av = avma;
783 63 : GEN q = powiu(p, e);
784 63 : GEN f = FpXQX_normalize(F, T, p);
785 63 : GEN g = FpXQX_normalize(FpXQX_split_part(f, T, p), T, p);
786 : GEN S;
787 63 : 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 63 : S = FpXQX_roots(g, T, p);
793 63 : return gc_upto(av, ZpXQX_liftroots_full(F, S, T, q, p, e));
794 : }
795 :
796 : GEN
797 6941 : ZqX_roots(GEN F, GEN T, GEN p, long e)
798 : {
799 6941 : 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 gc_INT(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 189 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
846 : {
847 189 : pari_sp av = avma, av2;
848 189 : GEN pv = p, q, W, df, Tq, fr, dfr;
849 : ulong mask;
850 189 : a = Fq_red(a, T, p);
851 189 : if (e <= v+1) return a;
852 189 : df = RgX_deriv(f);
853 189 : if (v) { pv = powiu(p,v); df = ZXX_Z_divexact(df, pv); }
854 189 : mask = quadratic_prec_mask(e-v);
855 189 : Tq = FpXT_red(T, p); dfr = FpXQX_red(df, Tq, p);
856 189 : W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
857 189 : q = p;
858 189 : av2 = avma;
859 : for (;;)
860 161 : {
861 : GEN u, fa, qv, q2v, q2, Tq2;
862 350 : q2 = q; q = sqri(q);
863 350 : if (mask & 1) q = diviiexact(q,p);
864 350 : mask >>= 1;
865 350 : if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
866 350 : else { qv = q; q2v = q2; }
867 350 : Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
868 350 : fr = FpXQX_red(f, Tq, qv);
869 350 : fa = FqX_eval(fr, a, Tq, qv);
870 350 : fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
871 350 : a = Fq_sub(a, Fq_mul(Fq_mul(W,fa,Tq2,q2v),q2, Tq,qv), Tq, qv);
872 350 : if (mask == 1) return gc_upto(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 : (void)gc_all(av2, 3, &a, &W, &q);
881 : }
882 : }
883 : }
884 :
885 : GEN
886 189 : 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 21103 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
903 : {
904 21103 : pari_sp av = avma;
905 21103 : long vT = get_FpX_var(T);
906 : GEN bn, bdi;
907 21103 : GEN bd = ZX_Z_add(x, gen_1);
908 21103 : if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
909 : {
910 7138 : bn = ZX_shifti(x,-(k+1));
911 7138 : bdi= ZpXQ_invlift(ZX_shifti(bd,-1), pol_1(vT), T, p, e);
912 : }
913 : else
914 : {
915 13965 : bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
916 13965 : bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
917 : }
918 21103 : return gc_upto(av, FpXQ_mul(bn, bdi, T, pe));
919 : }
920 :
921 : /* Assume p odd, a = 1 [p], return log(a) */
922 : GEN
923 21103 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
924 : {
925 21103 : pari_sp av = avma;
926 : pari_timer ti;
927 21103 : long is2 = absequaliu(p,2);
928 21103 : ulong pp = is2 ? 0: itou_or_0(p);
929 21103 : double lp = is2 ? 1: pp ? log2(pp): expi(p);
930 21103 : long k = maxss(1, (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
931 : GEN ak, s, b, pol;
932 21103 : long e = is2 ? N-1: N;
933 21103 : long i, l = (e-2)/(2*(k+is2));
934 21103 : GEN pe = powiu(p,e);
935 21103 : GEN TNk, pNk = powiu(p,N+k);
936 21103 : if( DEBUGLEVEL>=3) timer_start(&ti);
937 21103 : TNk = FpX_get_red(get_FpX_mod(T), pNk);
938 21103 : ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
939 21103 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
940 21103 : b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
941 21103 : if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
942 21103 : pol= cgetg(l+3,t_POL);
943 21103 : pol[1] = evalsigne(1)|evalvarn(0);
944 58222 : for(i=0; i<=l; i++)
945 : {
946 : GEN g;
947 37119 : ulong z = 2*i+1;
948 37119 : if (pp)
949 : {
950 26341 : long w = u_lvalrem(z, pp, &z);
951 26341 : g = powuu(pp,2*i*k-w);
952 : }
953 10778 : else g = powiu(p,2*i*k);
954 37119 : gel(pol,i+2) = Fp_divu(g, z,pe);
955 : }
956 21103 : if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
957 21103 : s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T, pe);
958 21103 : if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
959 21103 : s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
960 21103 : return gc_upto(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 gc_upto(av, FpC_red(ZM_add(VN2, ZM_Z_mul(VM, q2)), q));
990 : }
991 :
992 : GEN
993 200858 : 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 200858 : pari_sp av = avma;
998 : long N2, M;
999 : GEN VN2, V2, VM, bil;
1000 : GEN q2, qM;
1001 200858 : V = FpX_red(V, q);
1002 200858 : if (N == 1) return invl(E, V);
1003 45829 : N2 = (N + 1)>>1; M = N - N2;
1004 45829 : F = FpXT_red(F, q);
1005 45829 : qM = powiu(p, M);
1006 45829 : q2 = M == N2? qM: mulii(qM, p);
1007 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
1008 45829 : VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
1009 45829 : bil = lin(E, F, VN2, q);
1010 45829 : V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
1011 45829 : VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
1012 45829 : return gc_upto(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 : (void)gc_all(av, 2, &x, &q);
1050 : }
1051 : }
1052 0 : return gc_upto(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 :
1076 : GEN
1077 0 : ZpM_ZpC_gauss(GEN F, GEN V, GEN q, GEN p, long N)
1078 : {
1079 0 : pari_sp av = avma;
1080 : long N2, M;
1081 : GEN VN2, V2, VM, bil;
1082 : GEN q2, qM;
1083 0 : F = FpM_red(F, q);
1084 0 : V = FpC_red(V, q);
1085 0 : if (N == 1) return gc_upto(av, FpM_FpC_gauss(F, V, p));
1086 0 : N2 = (N + 1)>>1; M = N - N2;
1087 0 : qM = powiu(p, M);
1088 0 : q2 = M == N2? qM: mulii(qM, p);
1089 : /* q2 = p^N2, qM = p^M, q = q2 * qM */
1090 0 : VN2 = ZpM_ZpC_gauss(F, V, q2, p, N2);
1091 0 : bil = FpM_FpC_mul(F, VN2, q);
1092 0 : V2 = ZC_Z_divexact(ZC_sub(V, bil), q2);
1093 0 : VM = ZpM_ZpC_gauss(F, V2, qM, p, M);
1094 0 : return gc_upto(av, FpC_red(ZC_add(VN2, ZC_Z_mul(VM, q2)), q));
1095 : }
1096 :
1097 : GEN
1098 266911 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
1099 : GEN eval(void *E, GEN f, GEN q),
1100 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1101 : {
1102 266911 : pari_sp ltop = avma, av;
1103 266911 : long N = 1, N2, M;
1104 : long mask;
1105 266911 : GEN q = p;
1106 266911 : if (n == 1) return gcopy(x);
1107 263082 : mask = quadratic_prec_mask(n);
1108 263082 : av = avma;
1109 930402 : while (mask > 1)
1110 : {
1111 : GEN qM, q2, v, V;
1112 667320 : N2 = N; N <<= 1;
1113 667320 : q2 = q;
1114 667320 : if (mask&1UL) { /* can never happen when q2 = p */
1115 276810 : N--; M = N2-1;
1116 276810 : qM = diviiexact(q2,p); /* > 1 */
1117 276810 : q = mulii(qM,q2);
1118 : } else {
1119 390510 : M = N2;
1120 390510 : qM = q2;
1121 390510 : q = sqri(q2);
1122 : }
1123 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1124 667318 : mask >>= 1;
1125 667318 : v = eval(E, x, q);
1126 667320 : V = ZX_Z_divexact(gel(v,1), q2);
1127 667318 : x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
1128 667319 : if (gc_needed(av, 1))
1129 : {
1130 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
1131 0 : (void)gc_all(av, 2, &x, &q);
1132 : }
1133 : }
1134 263082 : return gc_upto(ltop, x);
1135 : }
1136 :
1137 : struct _ZpXQ_inv
1138 : {
1139 : GEN T, a, p, n;
1140 : };
1141 :
1142 : static GEN
1143 523095 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1144 : {
1145 523095 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1146 523095 : GEN Tq = FpXT_red(d->T, q);
1147 : (void)M;
1148 523096 : return FpXQ_mul(V, gel(v,2), Tq, q);
1149 : }
1150 :
1151 : static GEN
1152 523096 : _inv_eval(void *E, GEN x, GEN q)
1153 : {
1154 523096 : struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
1155 523096 : GEN Tq = FpXT_red(d->T, q);
1156 523095 : GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
1157 523096 : return mkvec2(f, x);
1158 : }
1159 :
1160 : GEN
1161 201822 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
1162 : {
1163 : struct _ZpXQ_inv d;
1164 201822 : d.a = a; d.T = T; d.p = p;
1165 201822 : return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
1166 : }
1167 :
1168 : GEN
1169 180719 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
1170 : {
1171 180719 : pari_sp av=avma;
1172 : GEN ai;
1173 180719 : if (lgefint(p)==3)
1174 : {
1175 180691 : ulong pp = p[2];
1176 180691 : ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
1177 : } else
1178 28 : ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
1179 180719 : return gc_upto(av, ZpXQ_invlift(a, ai, T, p, e));
1180 : }
1181 :
1182 : GEN
1183 31605 : ZpXQ_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
1184 : {
1185 31605 : return FpXQ_mul(a, ZpXQ_inv(b, T, p, e), T, q);
1186 : }
1187 :
1188 : GEN
1189 146440 : ZpXQX_divrem(GEN x, GEN Sp, GEN T, GEN q, GEN p, long e, GEN *pr)
1190 : {
1191 146440 : pari_sp av = avma;
1192 146440 : GEN S = get_FpXQX_mod(Sp);
1193 146440 : GEN b = leading_coeff(S), bi;
1194 : GEN S2, Q;
1195 146440 : if (typ(b)==t_INT) return FpXQX_divrem(x, Sp, T, q, pr);
1196 146440 : bi = ZpXQ_inv(b, T, p, e);
1197 146440 : S2 = FqX_Fq_mul_to_monic(S, bi, T, q);
1198 146440 : Q = FpXQX_divrem(x, S2, T, q, pr);
1199 146440 : if (pr==ONLY_DIVIDES && !Q) { set_avma(av); return NULL; }
1200 146440 : if (pr==ONLY_REM || pr==ONLY_DIVIDES) return gc_upto(av, Q);
1201 146440 : Q = FpXQX_FpXQ_mul(Q, bi, T, q);
1202 146440 : return gc_all(av, 2, &Q, pr);
1203 : }
1204 :
1205 : GEN
1206 623 : ZpXQX_digits(GEN x, GEN B, GEN T, GEN q, GEN p, long e)
1207 : {
1208 623 : pari_sp av = avma;
1209 623 : GEN b = leading_coeff(B), bi;
1210 : GEN B2, P, V, W;
1211 : long i, lV;
1212 623 : if (typ(b)==t_INT) return FpXQX_digits(x, B, T, q);
1213 623 : bi = ZpXQ_inv(b, T, p, e);
1214 623 : B2 = FqX_Fq_mul_to_monic(B, bi, T, q);
1215 623 : V = FpXQX_digits(x, B2, T, q);
1216 623 : lV = lg(V)-1;
1217 623 : P = FpXQ_powers(bi, lV-1, T, q);
1218 623 : W = cgetg(lV+1, t_VEC);
1219 11200 : for(i=1; i<=lV; i++)
1220 10577 : gel(W, i) = FpXQX_FpXQ_mul(gel(V,i), gel(P, i), T, q);
1221 623 : return gc_upto(av, W);
1222 : }
1223 :
1224 : struct _ZpXQ_sqrtn
1225 : {
1226 : GEN T, a, n, ai;
1227 : };
1228 :
1229 : static GEN
1230 1806 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
1231 : {
1232 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1233 1806 : GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
1234 : (void)M;
1235 1806 : return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
1236 : }
1237 :
1238 : static GEN
1239 1806 : _sqrtn_eval(void *E, GEN x, GEN q)
1240 : {
1241 1806 : struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
1242 1806 : GEN Tq = FpX_red(d->T, q);
1243 1806 : GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
1244 1806 : return mkvec2(f, x);
1245 : }
1246 :
1247 : GEN
1248 1155 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1249 : {
1250 : struct _ZpXQ_sqrtn d;
1251 1155 : d.a = a; d.T = T; d.n = n;
1252 1155 : d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
1253 1155 : return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
1254 : }
1255 :
1256 : static GEN
1257 0 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
1258 :
1259 : GEN
1260 0 : Zq_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
1261 : {
1262 0 : return T? ZpXQ_sqrtnlift(to_ZX(a,varn(T)), n, to_ZX(x,varn(T)), T, p, e)
1263 0 : : Zp_sqrtnlift(a, n, x, p, e);
1264 : }
1265 :
1266 : GEN
1267 0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
1268 : {
1269 0 : pari_sp av = avma;
1270 0 : GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
1271 0 : if (!z) return NULL;
1272 0 : if (e <= 1) return gc_upto(av, z);
1273 0 : return gc_upto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
1274 : }
1275 :
1276 : GEN
1277 33756 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
1278 : GEN early(void *E, GEN x, GEN q))
1279 : {
1280 33756 : pari_sp ltop = avma, av;
1281 : long N, r;
1282 : long mask;
1283 : GEN q2, q, W, Q, Tq2, Tq, Pq;
1284 : pari_timer ti;
1285 33756 : T = FpX_get_red(T, powiu(p, n));
1286 33755 : if (n == 1) return gcopy(S);
1287 33755 : mask = quadratic_prec_mask(n);
1288 33756 : av = avma;
1289 33756 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1290 33755 : if (DEBUGLEVEL > 3) timer_start(&ti);
1291 33755 : Tq = FpXT_red(T,q);
1292 33754 : Tq2 = FpXT_red(Tq,q2);
1293 33755 : Pq = FpX_red(P,q);
1294 33754 : W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
1295 33756 : Q = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
1296 33756 : r = brent_kung_optpow(degpol(P), 4, 3);
1297 33756 : if (DEBUGLEVEL > 3)
1298 0 : err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
1299 : for (;;)
1300 38673 : {
1301 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
1302 72429 : H = FpXQ_mul(W, Q, Tq2, q2);
1303 72429 : Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
1304 72429 : if (DEBUGLEVEL > 3)
1305 0 : timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
1306 72429 : if (mask==1)
1307 6440 : return gc_upto(ltop, Sq);
1308 65989 : if (early)
1309 : {
1310 64330 : GEN Se = early(E, Sq, q);
1311 64330 : if (Se) return gc_upto(ltop, Se);
1312 : }
1313 38673 : qq = sqri(q); N <<= 1;
1314 38673 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1315 38673 : mask >>= 1;
1316 38673 : Pqq = FpX_red(P, qq);
1317 38671 : Tqq = FpXT_red(T, qq);
1318 38673 : Spow = FpXQ_powers(Sq, r, Tqq, qq);
1319 38672 : Q = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
1320 38673 : dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
1321 38673 : Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
1322 38672 : Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
1323 38672 : Wq = FpX_sub(W, Wq, q);
1324 38673 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
1325 38673 : if (gc_needed(av, 1))
1326 : {
1327 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_liftroot");
1328 0 : (void)gc_all(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
1329 : }
1330 : }
1331 : }
1332 :
1333 : GEN
1334 1729 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
1335 : {
1336 1729 : return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
1337 : }
1338 :
1339 : GEN
1340 301 : ZpX_Frobenius(GEN T, GEN p, long e)
1341 : {
1342 301 : return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
1343 : }
1344 :
1345 : GEN
1346 147 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
1347 : {
1348 147 : pari_sp av = avma;
1349 147 : GEN xp = ZpX_Frobenius(T, p, e);
1350 147 : GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
1351 147 : return gc_GEN(av, gel(z,2));
1352 : }
1353 :
1354 : GEN
1355 0 : ZpXQX_ZpXQXQ_liftroot(GEN P, GEN S, GEN U, GEN T, GEN p, long n)
1356 : {
1357 0 : pari_sp ltop = avma, av;
1358 : long N, r;
1359 : long mask;
1360 : GEN qn, q2, q, W, Q, Tq2, Tq, Pq, Uq, Uq2;
1361 : pari_timer ti;
1362 0 : qn = powiu(p, n);
1363 0 : T = FpX_get_red(T, qn);
1364 0 : U = FpXQX_get_red(U, T, qn);
1365 0 : if (n == 1) return gcopy(S);
1366 0 : mask = quadratic_prec_mask(n);
1367 0 : av = avma;
1368 0 : q2 = p; q = sqri(p); mask >>= 1; N = 2;
1369 0 : if (DEBUGLEVEL > 3) timer_start(&ti);
1370 0 : Tq = FpXT_red(T,q);
1371 0 : Uq = FpXQXT_red(U, Tq, q);
1372 0 : Tq2 = FpXT_red(Tq, q2);
1373 0 : Uq2 = FpXQXT_red(U, Tq2, q2);
1374 0 : Pq = FpXQX_red(P, Tq, q);
1375 0 : W = FpXQXQ_inv(FpXQX_FpXQXQ_eval(FpXX_deriv(P,q2), S, Uq2, Tq2, q2), Uq2, Tq2, q2);
1376 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQ_eval(Pq, S, Uq, Tq, q), q2);
1377 0 : r = brent_kung_optpow(degpol(P), 4, 3);
1378 0 : if (DEBUGLEVEL > 3)
1379 0 : err_printf("ZpXQX_ZpXQXQ_liftroot: lifting to prec %ld\n",n);
1380 : for (;;)
1381 0 : {
1382 : GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq, Uqq;
1383 0 : H = FpXQXQ_mul(W, Q, Uq2, Tq2, q2);
1384 0 : Sq = FpXX_sub(S, ZXX_Z_mul(H, q2), q);
1385 0 : if (DEBUGLEVEL > 3)
1386 0 : timer_printf(&ti,"ZpXQX_ZpXQXQ_liftroot: reaching prec %ld",N);
1387 0 : if (mask==1)
1388 0 : return gc_upto(ltop, Sq);
1389 0 : qq = sqri(q); N <<= 1;
1390 0 : if (mask&1UL) { qq = diviiexact(qq, p); N--; }
1391 0 : mask >>= 1;
1392 0 : Tqq = FpXT_red(T, qq);
1393 0 : Uqq = FpXQXT_red(U, Tqq, qq);
1394 0 : Pqq = FpXQX_red(P, Tqq, qq);
1395 0 : Spow = FpXQXQ_powers(Sq, r, Uqq, Tqq, qq);
1396 0 : Q = ZXX_Z_divexact(FpXQX_FpXQXQV_eval(Pqq, Spow, Uqq, Tqq, qq), q);
1397 0 : dP = FpXQX_FpXQXQV_eval(FpXX_deriv(Pq, q), FpXQXV_red(Spow, Tq, q), Uq, Tq, q);
1398 0 : Wq = ZXX_Z_divexact(gsub(FpXQXQ_mul(W, dP, Uq, Tq, q), gen_1), q2);
1399 0 : Wq = ZXX_Z_mul(FpXQXQ_mul(W, Wq, Uq2, Tq2, q2), q2);
1400 0 : Wq = FpXX_sub(W, Wq, q);
1401 0 : S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Uq2 = Uq; Uq = Uqq; Pq = Pqq;
1402 0 : if (gc_needed(av, 1))
1403 : {
1404 0 : if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_ZpXQXQ_liftroot");
1405 0 : (void)gc_all(av, 10, &S, &W, &Q, &Uq2, &Uq, &Tq2, &Tq, &Pq, &q, &q2);
1406 : }
1407 : }
1408 : }
1409 :
1410 : GEN
1411 35 : ZqX_ZqXQ_liftroot(GEN f, GEN a, GEN P, GEN T, GEN p, long e)
1412 0 : { return T ? ZpXQX_ZpXQXQ_liftroot(f, a, P, T, p, e)
1413 35 : : ZpX_ZpXQ_liftroot(f, a, P, p, e); }
1414 :
1415 : /* Canonical lift of polynomial */
1416 :
1417 : static GEN
1418 11130 : _can_invl(void *E, GEN V) {(void) E; return V; }
1419 :
1420 : static GEN
1421 3654 : _can_lin(void *E, GEN F, GEN V, GEN q)
1422 : {
1423 3654 : GEN v = RgX_splitting(V, 3);
1424 : (void) E;
1425 3654 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1426 : }
1427 :
1428 : static GEN
1429 7476 : _can_iter(void *E, GEN f, GEN q)
1430 : {
1431 7476 : GEN h = RgX_splitting(f,3);
1432 7476 : GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
1433 7476 : GEN h12 = ZX_mul(gel(h,1), gel(h,2));
1434 7476 : GEN h13 = ZX_mul(gel(h,1), gel(h,3));
1435 7476 : GEN h23 = ZX_mul(gel(h,2), gel(h,3));
1436 7476 : GEN h1c = ZX_mul(gel(h,1), h1s);
1437 7476 : GEN h3c = ZX_mul(gel(h,3), h3s);
1438 7476 : GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
1439 7476 : GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
1440 : (void) E;
1441 7476 : return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
1442 : }
1443 :
1444 : static GEN
1445 7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
1446 : {
1447 7476 : GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
1448 7476 : GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
1449 7476 : GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
1450 : ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
1451 : (void)E;
1452 7476 : return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
1453 : _can_lin, _can_invl);
1454 : }
1455 :
1456 : static GEN
1457 3717 : F3x_frobeniuslift(GEN P, long n)
1458 3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
1459 :
1460 : static GEN
1461 0 : _can5_invl(void *E, GEN V) {(void) E; return V; }
1462 :
1463 : static GEN
1464 0 : _can5_lin(void *E, GEN F, GEN V, GEN q)
1465 : {
1466 0 : ulong p = *(ulong*)E;
1467 0 : GEN v = RgX_splitting(V, p);
1468 0 : return FpX_sub(V,ZXV_dotproduct(v, F), q);
1469 : }
1470 :
1471 : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
1472 : static GEN
1473 0 : _shift(GEN P, long n, ulong p, long v)
1474 : {
1475 0 : long i, l=lg(P);
1476 0 : GEN r = cgetg(l,t_POL); r[1] = P[1];
1477 0 : for(i=2;i<l;i++)
1478 : {
1479 0 : long s = n*(i-2)%p;
1480 0 : GEN ci = gel(P,i);
1481 0 : if (typ(ci)==t_INT)
1482 0 : gel(r,i) = monomial(ci, s, v);
1483 : else
1484 0 : gel(r,i) = RgX_rotate_shallow(ci, s, p);
1485 : }
1486 0 : return FpXX_renormalize(r, l);
1487 : }
1488 :
1489 : static GEN
1490 0 : _can5_iter(void *E, GEN f, GEN q)
1491 : {
1492 0 : pari_sp av = avma;
1493 0 : ulong p = *(ulong*)E;
1494 0 : long i, vT = fetch_var();
1495 : GEN N, P, d, V;
1496 0 : GEN T = ZX_Z_sub(pol_xn(p,vT),gen_1);
1497 0 : N = FpXQX_mul(_shift(f, 1, p, vT), _shift(f, 2, p, vT), T, q);
1498 0 : N = FpXQX_mul(N,_shift(N, 2, p, vT), T, q);
1499 0 : N = ZXX_evalx0(FpXQX_red(N,polcyclo(p,vT),q));
1500 0 : P = FpX_mul(N,f,q);
1501 0 : P = RgX_deflate(P, p);
1502 0 : d = RgX_splitting(N, p);
1503 0 : V = cgetg(p+1,t_VEC);
1504 0 : gel(V,1) = ZX_mulu(gel(d,1), p);
1505 0 : for(i=2; i<= (long)p; i++)
1506 0 : gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
1507 0 : (void)delete_var(); return gc_GEN(av, mkvec2(ZX_sub(f,P),V));
1508 : }
1509 :
1510 : static GEN
1511 0 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
1512 : {
1513 0 : ulong p = *(long*)E;
1514 0 : return gen_ZpX_Dixon(gel(v,2), H, qM, utoipos(p), M, E, _can5_lin, _can5_invl);
1515 : }
1516 :
1517 : struct _canlarge
1518 : {
1519 : ulong p, pi;
1520 : GEN sqx, Tp;
1521 : };
1522 :
1523 : static GEN
1524 29918 : FpXQ_inflate(GEN P, ulong p, GEN P2, GEN q)
1525 : {
1526 29918 : long n = degpol(P);
1527 29918 : return p >= 4*usqrt(n)
1528 12117 : ? FpX_FpXQ_eval(P, FpXQ_powu(pol_x(varn(P)), p, P2, q), P2, q)
1529 42035 : : FpX_rem(RgX_inflate(P, p), P2, q);
1530 : }
1531 :
1532 : static GEN
1533 29918 : _canlarge_invl(void *E, GEN x)
1534 : {
1535 29918 : struct _canlarge *D = (struct _canlarge *) E;
1536 29918 : GEN T = D->Tp;
1537 29918 : ulong p = D->p, pi = D->pi;
1538 29918 : GEN xp = ZX_to_Flx(x, p);
1539 29918 : return Flx_to_ZX_inplace(
1540 29918 : D->sqx ? Flxq_lroot_fast_pre(xp, D->sqx, T, p, pi)
1541 29197 : : Flxq_lroot_pre(xp, T, p, pi));
1542 : }
1543 :
1544 : static GEN
1545 9212 : _canlarge_lin(void *E, GEN F, GEN H, GEN q)
1546 : {
1547 9212 : pari_sp av = avma;
1548 9212 : struct _canlarge * D = (struct _canlarge *) E;
1549 9212 : ulong p = D->p;
1550 9212 : GEN Q = gel(F,2), P = gel(F,3);
1551 9212 : GEN Hp = FpXQ_inflate(H, p, P, q);
1552 9212 : GEN lin = FpX_sub(Hp, FpXQ_mul(H,Q,P,q), q);
1553 9212 : return gc_upto(av, lin);
1554 : }
1555 :
1556 : static GEN
1557 20706 : _canlarge_iter(void *E, GEN P, GEN q)
1558 : {
1559 20706 : pari_sp av = avma;
1560 20706 : struct _canlarge * D = (struct _canlarge *) E;
1561 20706 : ulong p = D->p;
1562 20706 : GEN R2 = FpXQ_inflate(P, p, FpX_sqr(P, q), q);
1563 20706 : GEN Pq = FpX_get_red(P, q);
1564 20706 : GEN R, Q = FpX_divrem(R2, Pq, q, &R);
1565 20706 : return gc_GEN(av, mkvec3(R, FpX_rem(Q, Pq, q), Pq));
1566 : }
1567 :
1568 : static GEN
1569 20706 : _canlarge_invd(void *E, GEN H, GEN v, GEN qM, long M)
1570 : {
1571 20706 : struct _canlarge * D = (struct _canlarge *) E;
1572 20706 : ulong p = D->p;
1573 20706 : return gen_ZpX_Dixon(v, H, qM, utoipos(p), M, E, _canlarge_lin, _canlarge_invl);
1574 : }
1575 :
1576 : static GEN
1577 10269 : Flx_Teichmuller_large(GEN T, ulong p, long n)
1578 : {
1579 : struct _canlarge D;
1580 10269 : ulong pi = get_Fl_red(p);
1581 10269 : D.p = p; D.pi = pi; D.Tp = get_Flx_mod(T);
1582 10269 : if ((ulong)degpol(D.Tp) > p)
1583 : {
1584 42 : GEN lr = Flxq_lroot_pre(polx_Flx(D.Tp[1]), T, p, pi);
1585 42 : D.sqx = Flxq_powers_pre(lr, p-1, T, p, pi);
1586 10227 : } else D.sqx = NULL;
1587 10269 : return gen_ZpX_Newton(Flx_to_ZX(D.Tp),utoipos(p), n, &D, _canlarge_iter, _canlarge_invd);
1588 : }
1589 :
1590 : GEN
1591 13986 : Flx_Teichmuller(GEN P, ulong p, long n)
1592 : {
1593 13986 : long d = degpol(P);
1594 24255 : return p==3 ? F3x_frobeniuslift(P,n):
1595 10192 : (p==5 && d>=311) ?
1596 0 : gen_ZpX_Newton(Flx_to_ZX(P),utoipos(p), n, &p, _can5_iter, _can5_invd)
1597 20461 : : Flx_Teichmuller_large(P,p,n);
1598 : }
1599 :
1600 : GEN
1601 35 : polteichmuller(GEN P, ulong p, long n)
1602 : {
1603 35 : pari_sp av = avma;
1604 35 : GEN q = NULL;
1605 35 : if (typ(P)!=t_POL || !RgX_is_FpX(P,&q)) pari_err_TYPE("polteichmuller",P);
1606 35 : if (q && !equaliu(q,p)) pari_err_MODULUS("polteichmuller",q,utoi(p));
1607 35 : if (n <= 0)
1608 0 : pari_err_DOMAIN("polteichmuller", "precision", "<=",gen_0,stoi(n));
1609 35 : if (signe(P)==0) pari_err_ROOTS0("polteichmuller");
1610 49 : return gc_upto(av, p==2 ? F2x_Teichmuller(RgX_to_F2x(P), n)
1611 21 : : Flx_Teichmuller(RgX_to_Flx(P, p), p, n));
1612 : }
1613 :
1614 : GEN
1615 0 : gen_ZpXX_Newton(GEN x, GEN p, long n, void *E,
1616 : GEN eval(void *E, GEN f, GEN q),
1617 : GEN invd(void *E, GEN V, GEN v, GEN q, long M))
1618 : {
1619 0 : pari_sp ltop = avma, av;
1620 0 : long N = 1, N2, M;
1621 : long mask;
1622 0 : GEN q = p;
1623 0 : if (n == 1) return gcopy(x);
1624 0 : mask = quadratic_prec_mask(n);
1625 0 : av = avma;
1626 0 : while (mask > 1)
1627 : {
1628 : GEN qM, q2, v, V;
1629 0 : N2 = N; N <<= 1;
1630 0 : q2 = q;
1631 0 : if (mask&1UL) { /* can never happen when q2 = p */
1632 0 : N--; M = N2-1;
1633 0 : qM = diviiexact(q2,p); /* > 1 */
1634 0 : q = mulii(qM,q2);
1635 : } else {
1636 0 : M = N2;
1637 0 : qM = q2;
1638 0 : q = sqri(q2);
1639 : }
1640 : /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
1641 0 : mask >>= 1;
1642 0 : v = eval(E, x, q);
1643 0 : V = ZXX_Z_divexact(gel(v,1), q2);
1644 0 : x = FpXX_sub(x, ZXX_Z_mul(invd(E, V, v, qM, M), q2), q);
1645 0 : if (gc_needed(av, 1))
1646 : {
1647 0 : if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
1648 0 : (void)gc_all(av, 2, &x, &q);
1649 : }
1650 : }
1651 0 : return gc_upto(ltop, x);
1652 : }
1653 :
1654 : struct _ZpXQXQ_inv
1655 : {
1656 : GEN S, T, a, p, n;
1657 : };
1658 :
1659 : static GEN
1660 0 : _ZpXQXQ_inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
1661 : {
1662 0 : struct _ZpXQXQ_inv *d = (struct _ZpXQXQ_inv *) E;
1663 0 : GEN Tq = FpXT_red(d->T, q);
1664 0 : GEN Sq = FpXQXT_red(d->S, Tq, q);
1665 : (void)M;
1666 0 : return FpXQXQ_mul(V, gel(v,2), Sq, Tq, q);
1667 : }
1668 :
1669 : GEN FpXX_Fp_sub(GEN x, GEN y, GEN p);
1670 :
1671 : static GEN
1672 0 : _ZpXQXQ_inv_eval(void *E, GEN x, GEN q)
1673 : {
1674 0 : struct _ZpXQXQ_inv *d = (struct _ZpXQXQ_inv *) E;
1675 0 : GEN Tq = FpXT_red(d->T, q);
1676 0 : GEN Sq = FpXQXT_red(d->S, Tq, q);
1677 0 : GEN f = FpXX_Fp_sub(FpXQXQ_mul(x, FpXX_red(d->a, q), Sq, Tq, q), gen_1, q);
1678 0 : return mkvec2(f, x);
1679 : }
1680 :
1681 : GEN
1682 0 : ZpXQXQ_invlift(GEN a, GEN x, GEN S, GEN T, GEN p, long e)
1683 : {
1684 : struct _ZpXQXQ_inv d;
1685 0 : d.a = a; d.S = S; d.T = T; d.p = p;
1686 0 : return gen_ZpXX_Newton(x, p, e, &d, _ZpXQXQ_inv_eval, _ZpXQXQ_inv_invd);
1687 : }
1688 :
1689 : GEN
1690 0 : ZpXQXQ_inv(GEN a, GEN S, GEN T, GEN p, long e)
1691 : {
1692 0 : pari_sp av=avma;
1693 : GEN ai;
1694 0 : if (lgefint(p)==3)
1695 : {
1696 0 : ulong pp = p[2];
1697 0 : GEN Tp = ZXT_to_FlxT(T, pp);
1698 0 : long v = get_FpX_var(Tp);
1699 0 : GEN Sp = ZXXT_to_FlxXT(S, pp, v);
1700 0 : ai = FlxX_to_ZXX(FlxqXQ_inv(ZXX_to_FlxX(a, pp, v), Sp, Tp, pp));
1701 : } else
1702 : {
1703 0 : GEN Tp = FpXT_red(T, p);
1704 0 : GEN Sp = FpXQXT_red(S, Tp, p);
1705 0 : ai = FpXQXQ_inv(FpXX_red(a,p), Sp, Tp, p);
1706 : }
1707 0 : return gc_upto(av, ZpXQXQ_invlift(a, ai, S, T, p, e));
1708 : }
1709 :
1710 : GEN
1711 0 : ZpXQXQ_div(GEN a, GEN b, GEN S, GEN T, GEN q, GEN p, long e)
1712 : {
1713 0 : return FpXQXQ_mul(a, ZpXQXQ_inv(b, S, T, p, e), S, T, q);
1714 : }
|
