Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - Zp.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 764 860 88.8 %
Date: 2022-07-03 07:33:15 Functions: 73 79 92.4 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13