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.16.2 lcov report (development 29454-e7b6e21fb9) Lines: 734 860 85.3 %
Date: 2024-07-18 09:03:47 Functions: 69 79 87.3 %
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     4669936 : quadratic_prec_mask(long n)
      38             : {
      39     4669936 :   long a = n, i;
      40     4669936 :   ulong mask = 0;
      41    14109279 :   for(i = 1;; i++, mask <<= 1)
      42             :   {
      43    14109279 :     mask |= (a&1); a = (a+1)>>1;
      44    14109279 :     if (a==1) return mask | (1UL << i);
      45             :   }
      46             : }
      47             : 
      48             : /***********************************************************************/
      49             : /**                                                                   **/
      50             : /**                            Zp                                     **/
      51             : /**                                                                   **/
      52             : /***********************************************************************/
      53             : 
      54             : static GEN
      55     1448447 : Zp_divlift(GEN b, GEN a, GEN x, GEN p, long n)
      56             : {
      57     1448447 :   pari_sp ltop = avma, av;
      58             :   ulong mask;
      59     1448447 :   GEN q = p;
      60     1448447 :   if (n == 1) return gcopy(x);
      61     1447908 :   mask = quadratic_prec_mask(n);
      62     1447908 :   av = avma;
      63     8842102 :   while (mask > 1)
      64             :   {
      65     7394195 :     GEN v, q2 = q;
      66     7394195 :     q = sqri(q);
      67     7394189 :     if (mask & 1UL) q = diviiexact(q,p);
      68     7394188 :     mask >>= 1;
      69     7394188 :     if (mask > 1 || !b)
      70             :     {
      71     7286481 :       v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
      72     7286490 :       x = Fp_sub(x, Fp_mul(v, x, q), q);
      73             :     }
      74             :     else
      75             :     {
      76      107707 :       GEN y = Fp_mul(x, b, q), yt = modii(y, q2);
      77      107709 :       v = Fp_sub(Fp_mul(x, modii(a, q), q), gen_1, q);
      78      107710 :       x = Fp_sub(y, Fp_mul(v, yt, q), q);
      79             :     }
      80     7394196 :     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     1447907 :   return gerepileupto(ltop, x);
      87             : }
      88             : 
      89             : GEN
      90     1340740 : Zp_invlift(GEN a, GEN x, GEN p, long e)
      91     1340740 : { return Zp_divlift(NULL, a, x, p, e); }
      92             : 
      93             : GEN
      94     1340740 : Zp_inv(GEN a, GEN p, long e)
      95             : {
      96     1340740 :   pari_sp av=avma;
      97             :   GEN ai;
      98     1340740 :   if (lgefint(p)==3)
      99             :   {
     100     1340477 :     ulong pp = p[2];
     101     1340477 :     ai = utoi(Fl_inv(umodiu(a,pp), pp));
     102             :   } else
     103         263 :     ai = Fp_inv(modii(a, p), p);
     104     1340740 :   return gerepileupto(av, Zp_invlift(a, ai, p, e));
     105             : }
     106             : 
     107             : GEN
     108      107707 : Zp_div(GEN b, GEN a, GEN p, long e)
     109             : {
     110      107707 :   pari_sp av=avma;
     111             :   GEN ai;
     112      107707 :   if (lgefint(p)==3)
     113             :   {
     114      107644 :     ulong pp = p[2];
     115      107644 :     ai = utoi(Fl_inv(umodiu(a,pp), pp));
     116             :   } else
     117          63 :     ai = Fp_inv(modii(a, p), p);
     118      107707 :   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      248895 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
     134             : {
     135      248895 :   pari_sp av = avma;
     136             :   int nis2, pis2;
     137             :   GEN q, w, n_1;
     138             :   ulong mask;
     139             : 
     140      248895 :   if (e == 1) return icopy(a);
     141      177955 :   nis2 = equaliu(n, 2); n_1 = nis2? NULL: subiu(n,1);
     142      177955 :   pis2 = equaliu(p, 2);
     143      177955 :   mask = quadratic_prec_mask(e);
     144      177955 :   w = nis2 ? shifti(a,1): Fp_mul(n, Fp_pow(a,n_1,p), p);
     145      177955 :   w = Fp_inv(w, p);
     146      177955 :   q = p; /* q = p^e; use e instead of q iff p = 2 */
     147      177955 :   e = 1;
     148             :   for(;;)
     149             :   {
     150      275605 :     if (pis2)
     151             :     {
     152         343 :       e <<= 1; if (mask & 1) e--;
     153         343 :       mask >>= 1;
     154             :       /* a -= w (a^n - b) */
     155         343 :       a = remi2n(subii(a, mulii(w, subii(Fp_pow2n(a, n, e), b))), e);
     156         343 :       if (mask == 1) break;
     157             :       /* w += w - w^2 n a^(n-1)*/
     158         133 :       w = subii(shifti(w,1), remi2n(mulii(remi2n(sqri(w), e),
     159             :                                           mulii(n, Fp_pow2n(a, n_1, e))), e));
     160         133 :       continue;
     161             :     }
     162      275262 :     q = sqri(q); if (mask & 1) q = diviiexact(q, p);
     163      275262 :     mask >>= 1;
     164      275262 :     if (lgefint(q) == 3 && lgefint(n) == 3)
     165       97370 :     {
     166      274818 :       ulong Q = uel(q,2), N = uel(n,2);
     167      274818 :       ulong A = umodiu(a, Q);
     168      274818 :       ulong B = umodiu(b, Q);
     169      274818 :       ulong W = umodiu(w, Q);
     170             : 
     171      274818 :       A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
     172      274818 :       a = utoi(A);
     173      274818 :       if (mask == 1) break;
     174       97370 :       if (nis2)
     175       84293 :         W = Fl_double(Fl_sub(W, Fl_mul(Fl_sqr(W,Q), A, Q), Q), Q);
     176             :       else
     177       13077 :         W = Fl_sub(Fl_double(W,Q),
     178             :                    Fl_mul(Fl_sqr(W,Q), Fl_mul(N,Fl_powu(A, N-1, Q), Q), Q), Q);
     179       97370 :       w = utoi(W);
     180             :     }
     181             :     else
     182             :     {
     183             :       /* a -= w (a^n - b) */
     184         444 :       a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q), b))), q);
     185         444 :       if (mask == 1) break;
     186             :       /* w += w - w^2 n a^(n-1)*/
     187         147 :       if (nis2)
     188          15 :         w = shifti(subii(w, Fp_mul(Fp_sqr(w,q), a, q)), 1);
     189             :       else
     190         132 :         w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
     191             :                                       mulii(n, Fp_pow(a,n_1,q)), q));
     192             :     }
     193             :   }
     194      177955 :   if (pis2 && signe(a) < 0) a = addii(a, int2n(e));
     195      177955 :   return gerepileuptoint(av, a);
     196             : }
     197             : 
     198             : /* ZpX_liftroot for the polynomial X^2-b */
     199             : GEN
     200      139972 : Zp_sqrtlift(GEN b, GEN a, GEN p, long e)
     201      139972 : { return Zp_sqrtnlift(b, gen_2, a, p, e); }
     202             : 
     203             : GEN
     204     1255056 : Zp_sqrt(GEN x, GEN p, long e)
     205             : {
     206             :   pari_sp av;
     207             :   GEN z;
     208     1255056 :   if (absequaliu(p,2)) return Z2_sqrt(x,e);
     209     1253509 :   av = avma; z = Fp_sqrt(Fp_red(x, p), p);
     210     1253510 :   if (!z) return NULL;
     211     1253461 :   if (e > 1) z = Zp_sqrtlift(x, z, p, e);
     212     1253461 :   return gerepileuptoint(av, z);
     213             : }
     214             : 
     215             : /* p-adic logarithm of a = 1 mod p (adapted from a C program of Xavier Caruso)
     216             :  * Algorithm:
     217             :  * 1. raise a at the power p^(v-1)  in order to make it closer to 1
     218             :  * 2. write the new a as a product
     219             :  *      1 / a = (1 - a_0*p^v) (1 - a_1*p^(2*v) (1 - a_2*p^(4*v) ...
     220             :  *    with 0 <= a_i < p^(v*2^i).
     221             :  * 3. compute each log(1 - a_i*p^(v*2^i)) using Taylor expansion
     222             :  *    and binary splitting */
     223             : GEN
     224      195791 : Zp_log(GEN a, GEN p, ulong e)
     225             : {
     226      195791 :   pari_sp av = avma;
     227      195791 :   ulong v, N, Np, trunc, pp = itou_or_0(p);
     228      195791 :   GEN pe, pv, trunc_mod, num, den, ans = gen_0;
     229             : 
     230      195791 :   if (equali1(a)) return ans; /* very frequent! */
     231             :   /* First make the argument closer to 1 by raising it to the p^(v-1) */
     232      194517 :   v = pp? ulogint(e, pp): 0;  /* v here is v-1 */
     233      194517 :   pe = powiu(p,e); pv = powiu(p,v);
     234      194517 :   a = Fp_pow(a, pv, mulii(pe, pv));
     235      194517 :   e += v;
     236             : 
     237             :   /* Where do we truncate the Taylor expansion */
     238      194517 :   N = e + v; N /= ++v; /* note the ++v */
     239      194517 :   Np = N;
     240             :   while(1)
     241           0 :   {
     242      194517 :     ulong e = Np;
     243      194517 :     if (pp) e += ulogint(N, pp) / v;
     244      194517 :     if (e == N) break;
     245           0 :     N = e;
     246             :   }
     247             : 
     248      194517 :   num = cgetg(N+1, t_INT);
     249      194517 :   den = cgetg(N+1, t_INT);
     250      194517 :   trunc = v << 1;
     251      194517 :   trunc_mod = powiu(p, trunc);
     252             :   while(1)
     253      784162 :   { /* compute f = 1 - a_i*p^((v+1)*2^i); trunc_mod = p^((v+1)*2^(i+1)) */
     254      978679 :     GEN f = modii(a, trunc_mod);
     255      978679 :     if (!equali1(f))
     256             :     {
     257      975609 :       ulong i, step = 1;
     258             :       GEN h, hpow;
     259             : 
     260      975609 :       f = subui(2, f);
     261      975608 :       a = mulii(a, f);
     262             : 
     263             :       /* compute the Taylor expansion of log(f), over Q for now */
     264    10526144 :       for (i = 1; i <= N; i++) { gel(num,i) = gen_1; gel(den,i) = utoipos(i); }
     265      975609 :       hpow = h = subui(1, f); /* h = a_i*p^(2^i) */
     266             :       while(1)
     267             :       {
     268    11352643 :         for (i = 1; i <= N - step; i += step << 1)
     269             :         {
     270     8574928 :           GEN t = mulii(mulii(hpow, gel(num,i+step)), gel(den,i));
     271     8574928 :           gel(num,i) = mulii(gel(num,i), gel(den,i+step));
     272     8574928 :           gel(num,i) = addii(gel(num,i), t);
     273     8574928 :           gel(den,i) = mulii(gel(den,i), gel(den,i+step));
     274             :         }
     275     2777715 :         step <<= 1; if (step >= N) break;
     276     1802106 :         hpow = sqri(hpow);
     277             :       }
     278             : 
     279      975609 :       if (pp)
     280             :       { /* simplify the fraction */
     281      975353 :         GEN d = powuu(pp, factorial_lval(N, pp));
     282      975353 :         gel(num,1) = diviiexact(gel(num,1), d);
     283      975353 :         gel(den,1) = diviiexact(gel(den,1), d);
     284             :       }
     285             : 
     286      975609 :       h = diviiexact(h, pv);
     287      975609 :       ans = addii(ans, mulii(mulii(gel(num,1), h), Zp_inv(gel(den,1), p, e)));
     288             :     }
     289      978679 :     if (trunc > e) break;
     290      784162 :     trunc_mod = sqri(trunc_mod); trunc <<= 1; N >>= 1;
     291             :   }
     292      194517 :   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      107709 : Zp_exp(GEN a, GEN p, ulong e)
     301             : {
     302      107709 :   pari_sp av = avma;
     303      107709 :   ulong trunc, N = e, pp = itou_or_0(p);
     304      107709 :   GEN num, den, trunc_mod = NULL, denominator = gen_1, ans = gen_1;
     305      107709 :   GEN pe = powiu(p, e);
     306      107708 :   int pis2 = pp == 2;
     307             : 
     308             :   /* Where do we truncate the Taylor expansion */
     309      107708 :   if (!pis2) N += sdivsi(e, subis(p,2));
     310      107708 :   num = cgetg(N+2, t_VEC);
     311      107708 :   den = cgetg(N+2, t_VEC);
     312      107709 :   if (pis2) trunc = 4;
     313             :   else
     314             :   {
     315      101072 :     trunc = 2;
     316      101072 :     trunc_mod = sqri(p);
     317             :   }
     318             :   while(1)
     319      199490 :   {
     320      307197 :     GEN h, hpow, f = pis2? remi2n(a, trunc): modii(a, trunc_mod);
     321      307195 :     a = subii(a, f);
     322      307189 :     if (signe(f))
     323             :     {
     324      231258 :       ulong step = 1, i;
     325             :       /* Taylor expansion of exp(f), over Q for now */
     326      231258 :       gel(num,1) = gen_1;
     327      231258 :       gel(den,1) = gen_1;
     328     1009899 :       for (i = 2; i <= N+1; i++)
     329             :       {
     330      778634 :         gel(num,i) = gen_1;
     331      778634 :         gel(den,i) = utoipos(i-1);
     332             :       }
     333      231265 :       hpow = h = f;
     334             :       while(1)
     335             :       {
     336     1317121 :         for (i = 1; i <= N+1 - step; i += step << 1) {
     337      778428 :           gel(num,i) = mulii(gel(num,i), gel(den,i+step));
     338      778418 :           gel(num,i) = addii(gel(num,i), mulii(hpow, gel(num,i+step)));
     339      778423 :           gel(den,i) = mulii(gel(den,i), gel(den,i+step));
     340             :         }
     341      538693 :         step <<= 1; if (step > N) break;
     342      307433 :         hpow = sqri(hpow);
     343             :       }
     344             : 
     345             :       /* We simplify the fraction */
     346      231260 :       if (pp)
     347             :       {
     348      231134 :         GEN d = powuu(pp, factorial_lval(N, pp));
     349      231132 :         gel(num,1) = diviiexact(gel(num,1), d);
     350      231133 :         gel(den,1) = diviiexact(gel(den,1), d);
     351             :       }
     352             :       /* We add this contribution to exp(f) */
     353      231256 :       ans = Fp_mul(ans, gel(num,1), pe);
     354      231262 :       denominator = Fp_mul(denominator, gel(den,1), pe);
     355             :     }
     356             : 
     357      307197 :     if (trunc > e) break;
     358      199490 :     if (!pis2) trunc_mod = sqri(trunc_mod);
     359      199490 :     trunc <<= 1; N >>= 1;
     360             :   }
     361      107707 :   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      195908 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
     384             : {
     385      195908 :   long k = lg(a)-1;
     386             :   long i, j, s, minp, mind;
     387             : 
     388      832435 :   for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
     389             : 
     390      440619 :   for (j=1; j <= 2*k-5; j+=2,i++)
     391             :   {
     392      244712 :     minp = j;
     393      244712 :     mind = degpol(gel(V,j));
     394     1961482 :     for (s=j+1; s<i; s++)
     395     1716768 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
     396             : 
     397      244714 :     swap(gel(V,j), gel(V,minp));
     398      244714 :     lswap(link[j], link[minp]);
     399             : 
     400      244714 :     minp = j+1;
     401      244714 :     mind = degpol(gel(V,j+1));
     402     1716773 :     for (s=j+2; s<i; s++)
     403     1472060 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
     404             : 
     405      244713 :     swap(gel(V,j+1), gel(V,minp));
     406      244713 :     lswap(link[j+1], link[minp]);
     407             : 
     408      244713 :     gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
     409      244711 :     link[i] = j;
     410             :   }
     411             : 
     412      636514 :   for (j=1; j <= 2*k-3; j+=2)
     413             :   {
     414             :     GEN d, u, v;
     415      440615 :     d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
     416      440625 :     if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
     417      440618 :     d = gel(d,2);
     418      440618 :     if (!gequal1(d))
     419             :     {
     420      396730 :       if (typ(d)==t_POL)
     421             :       {
     422      132950 :         d = FpXQ_inv(d, T, p);
     423      132945 :         u = FqX_Fq_mul(u, d, T, p);
     424      132944 :         v = FqX_Fq_mul(v, d, T, p);
     425             :       }
     426             :       else
     427             :       {
     428      263780 :         d = Fp_inv(d, p);
     429      263778 :         u = FqX_Fp_mul(u, d, T,p);
     430      263777 :         v = FqX_Fp_mul(v, d, T,p);
     431             :       }
     432             :     }
     433      440607 :     gel(W,j) = u;
     434      440607 :     gel(W,j+1) = v;
     435             :   }
     436      195899 : }
     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     1577843 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
     442             : {
     443     1577843 :   pari_sp av = avma;
     444     1577843 :   long space = lg(f) * lgefint(p1);
     445             :   GEN a2, b2, g, z, s, t;
     446     1577843 :   GEN a = gel(V,j), b = gel(V,j+1);
     447     1577843 :   GEN u = gel(W,j), v = gel(W,j+1);
     448             : 
     449     1577843 :   (void)new_chunk(space); /* HACK */
     450     1577880 :   g = ZX_sub(f, ZX_mul(a,b));
     451     1577837 :   g = ZX_Z_divexact(g, p0);
     452     1577721 :   g = FpX_red(g, pd);
     453     1577743 :   z = FpX_mul(v,g, pd);
     454     1577882 :   t = FpX_divrem(z,a, pd, &s);
     455     1577892 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     456     1577797 :   t = FpX_red(t, pd);
     457     1577804 :   t = ZX_Z_mul(t,p0);
     458     1577772 :   s = ZX_Z_mul(s,p0);
     459     1577774 :   set_avma(av);
     460     1577764 :   a2 = ZX_add(a,s);
     461     1577825 :   b2 = ZX_add(b,t);
     462             : 
     463             :   /* already reduced mod p1 = pd p0 */
     464     1577857 :   gel(V,j)   = a2;
     465     1577857 :   gel(V,j+1) = b2;
     466     1577857 :   if (noinv) return;
     467             : 
     468     1273735 :   av = avma;
     469     1273735 :   (void)new_chunk(space); /* HACK */
     470     1273748 :   g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
     471     1273716 :   g = Z_ZX_sub(gen_1, g);
     472     1273777 :   g = ZX_Z_divexact(g, p0);
     473     1273650 :   g = FpX_red(g, pd);
     474     1273628 :   z = FpX_mul(v,g, pd);
     475     1273777 :   t = FpX_divrem(z,a, pd, &s);
     476     1273777 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     477     1273712 :   t = FpX_red(t, pd);
     478     1273724 :   t = ZX_Z_mul(t,p0);
     479     1273701 :   s = ZX_Z_mul(s,p0);
     480     1273693 :   set_avma(av);
     481     1273689 :   gel(W,j)   = ZX_add(u,t);
     482     1273725 :   gel(W,j+1) = ZX_add(v,s);
     483             : }
     484             : 
     485             : static void
     486      502645 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
     487             : {
     488      502645 :   pari_sp av = avma;
     489      502645 :   const long n = degpol(T1), vT = varn(T1);
     490      502646 :   long space = lg(f) * lgefint(p1) * lg(T1);
     491             :   GEN a2, b2, g, z, s, t;
     492      502646 :   GEN a = gel(V,j), b = gel(V,j+1);
     493      502646 :   GEN u = gel(W,j), v = gel(W,j+1);
     494             : 
     495      502646 :   (void)new_chunk(space); /* HACK */
     496      502656 :   g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
     497      502586 :   g = FpXQX_red(g, T1, p1);
     498      502521 :   g = RgX_Rg_divexact(g, p0);
     499      502471 :   z = FpXQX_mul(v,g, Td,pd);
     500      502664 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     501      502670 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     502      502624 :   t = Kronecker_to_ZXX(t, n, vT);
     503      502620 :   t = FpXQX_red(t, Td, pd);
     504      502601 :   t = RgX_Rg_mul(t,p0);
     505      502565 :   s = RgX_Rg_mul(s,p0);
     506      502541 :   set_avma(av);
     507             : 
     508      502539 :   a2 = RgX_add(a,s);
     509      502600 :   b2 = RgX_add(b,t);
     510             :   /* already reduced mod p1 = pd p0 */
     511      502644 :   gel(V,j)   = a2;
     512      502644 :   gel(V,j+1) = b2;
     513      502644 :   if (noinv) return;
     514             : 
     515      369047 :   av = avma;
     516      369047 :   (void)new_chunk(space); /* HACK */
     517      369046 :   g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
     518      368949 :   g = Kronecker_to_ZXX(g, n, vT);
     519      369042 :   g = Rg_RgX_sub(gen_1, g);
     520      369062 :   g = FpXQX_red(g, T1, p1);
     521      369005 :   g = RgX_Rg_divexact(g, p0);
     522      368959 :   z = FpXQX_mul(v,g, Td,pd);
     523      369077 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     524      369076 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     525      369049 :   t = Kronecker_to_ZXX(t, n, vT);
     526      369052 :   t = FpXQX_red(t, Td, pd);
     527      369031 :   t = RgX_Rg_mul(t,p0);
     528      369025 :   s = RgX_Rg_mul(s,p0);
     529      369028 :   set_avma(av);
     530      369032 :   gel(W,j)   = RgX_add(u,t);
     531      369016 :   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     3889571 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
     538             :                 GEN f, long j, int noinv)
     539             : {
     540     3889571 :   if (j < 0) return;
     541     1577668 :   ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
     542     1577827 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     543     1577893 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
     544             : }
     545             : static void
     546     1137369 : 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     1137369 :   if (j < 0) return;
     550      502437 :   ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
     551      502572 :   ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     552      502615 :   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      195936 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
     562             : {
     563      195936 :   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      195936 :   if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
     569      195936 :   if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
     570      195936 :   if (e0 == 1) return a;
     571             : 
     572      195908 :   if (DEBUGLEVEL > 3) timer_start(&Ti);
     573      195908 :   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      195908 :     e = 1;
     583      195908 :     v = cgetg(2*k-2 + 1, t_VEC);
     584      195908 :     w = cgetg(2*k-2 + 1, t_VEC);
     585      195908 :     link=cgetg(2*k-2 + 1, t_VECSMALL);
     586      195910 :     BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
     587      195901 :     if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
     588             :   }
     589      195901 :   mask = quadratic_prec_mask(e0);
     590      195902 :   eold = 1;
     591      195902 :   penew = NULL;
     592      195902 :   Tnew = NULL;
     593      195902 :   if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
     594     1062247 :   while (mask > 1)
     595             :   {
     596      866349 :     long enew = eold << 1;
     597      866349 :     if (mask & 1) enew--;
     598      866349 :     mask >>= 1;
     599      866349 :     if (enew >= e) { /* mask == 1: last iteration */
     600      866350 :       GEN peold = penew? penew: powiu(p, eold);
     601      866356 :       GEN Td = NULL, pd;
     602      866356 :       long d = enew - eold; /* = eold or eold-1 */
     603             :       /* lift from p^eold to p^enew */
     604      866356 :       pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
     605      866345 :       penew = mulii(peold,pd);
     606      866285 :       if (T) {
     607      132314 :         if (Tnew)
     608       97687 :           Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
     609             :         else
     610       34627 :           Td = FpX_red(T, peold);
     611      132308 :         Tnew = FpX_red(T, penew);
     612      132303 :         ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
     613             :                          (flag == 0 && mask == 1));
     614             :       }
     615             :       else
     616      733971 :         ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
     617             :                         (flag == 0 && mask == 1));
     618      866345 :       if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
     619             :     }
     620      866345 :     eold = enew;
     621             :   }
     622             : 
     623      195898 :   if (flag)
     624        1743 :     E = mkvec4(utoipos(e0), link, v, w);
     625             :   else
     626             :   {
     627      194155 :     E = cgetg(k+1, t_VEC);
     628     1069607 :     for (i = 1; i <= 2*k-2; i++)
     629             :     {
     630      875448 :       long t = link[i];
     631      875448 :       if (t < 0) gel(E,-t) = gel(v,i);
     632             :     }
     633             :   }
     634      195902 :   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      243000 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
     641             : {
     642      243000 :   pari_sp av = avma;
     643      243000 :   pol = FpX_normalize(pol, pe);
     644      243000 :   if (lg(Q) == 2) return mkvec(pol);
     645      159558 :   return gerepilecopy(av, MultiLift(pol, Q, NULL, p, e, 0));
     646             : }
     647             : 
     648             : GEN
     649       34636 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
     650             : {
     651       34636 :   pari_sp av = avma;
     652       34636 :   pol = FpXQX_normalize(pol, T, pe);
     653       34636 :   if (lg(Q) == 2) return mkvec(pol);
     654       34636 :   return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
     655             : }
     656             : 
     657             : GEN
     658       46508 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
     659       46508 : { return T ? ZpXQX_liftfact(f, a, T, pe, p, e): ZpX_liftfact(f, a, pe, p, e); }
     660             : GEN
     661         154 : ZqX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
     662         154 : { return T ? ZpXQX_liftroot(f, a,T , p, e): ZpX_liftroot(f, a, p, e); }
     663             : 
     664             : /* U = NULL treated as 1 */
     665             : static void
     666        7525 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
     667             : {
     668             :   GEN Q, R;
     669        7525 :   if (j < 0) return;
     670             : 
     671        2891 :   Q = FpX_mul(gel(v,j), gel(w,j), pe);
     672        2891 :   if (U)
     673             :   {
     674        1148 :     Q = FpXQ_mul(Q, U, f, pe);
     675        1148 :     R = FpX_sub(U, Q, pe);
     676             :   }
     677             :   else
     678        1743 :     R = Fp_FpX_sub(gen_1, Q, pe);
     679        2891 :   gel(w,j+1) = Q; /*  0 mod U v[j],  1 mod (1-U) v[j+1] */
     680        2891 :   gel(w,j) = R; /*  1 mod U v[j],  0 mod (1-U) v[j+1] */
     681        2891 :   BezoutPropagate(link, v, w, pe, R, f, link[j  ]);
     682        2891 :   BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
     683             : }
     684             : 
     685             : /* as above, but return the Bezout coefficients for the lifted modular factors
     686             :  *   U[i] = 1 mod Qlift[i]
     687             :  *          0 mod Qlift[j], j != i */
     688             : GEN
     689        1764 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
     690             : {
     691        1764 :   pari_sp av = avma;
     692             :   GEN E, link, v, w, pe;
     693        1764 :   long i, k = lg(Q)-1;
     694        1764 :   if (k == 1) retmkvec(pol_1(varn(pol)));
     695        1743 :   pe = powiu(p, e);
     696        1743 :   pol = FpX_normalize(pol, pe);
     697        1743 :   E = MultiLift(pol, Q, NULL, p, e, 1);
     698        1743 :   link = gel(E,2);
     699        1743 :   v    = gel(E,3);
     700        1743 :   w    = gel(E,4);
     701        1743 :   BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
     702        1743 :   E = cgetg(k+1, t_VEC);
     703        7525 :   for (i = 1; i <= 2*k-2; i++)
     704             :   {
     705        5782 :     long t = link[i];
     706        5782 :     if (t < 0) E[-t] = w[i];
     707             :   }
     708        1743 :   return gerepilecopy(av, E);
     709             : }
     710             : 
     711             : /* Front-end for ZpX_liftfact:
     712             :    lift the factorization of pol mod p given by L to p^N (if possible) */
     713             : GEN
     714          28 : polhensellift(GEN pol, GEN L, GEN Tp, long N)
     715             : {
     716             :   GEN T, p;
     717             :   long i, l;
     718          28 :   pari_sp av = avma;
     719             :   void (*chk)(GEN, const char*);
     720             : 
     721          28 :   if (typ(pol) != t_POL) pari_err_TYPE("polhensellift",pol);
     722          28 :   RgX_check_ZXX(pol, "polhensellift");
     723          28 :   if (!is_vec_t(typ(L)) || lg(L) < 3) pari_err_TYPE("polhensellift",L);
     724          28 :   if (N < 1) pari_err_DOMAIN("polhensellift", "precision", "<", gen_1,stoi(N));
     725          28 :   if (!ff_parse_Tp(Tp, &T, &p, 0)) pari_err_TYPE("polhensellift",Tp);
     726          28 :   chk = T? RgX_check_ZXX: RgX_check_ZX;
     727          28 :   l = lg(L); L = leafcopy(L);
     728          70 :   for (i = 1; i < l; i++)
     729             :   {
     730          49 :     GEN q = gel(L,i);
     731          49 :     if (typ(q) != t_POL) gel(L,i) = scalar_ZX_shallow(q, varn(pol));
     732          49 :     else chk(q, "polhensellift");
     733             :   }
     734          21 :   return gerepilecopy(av, ZqX_liftfact(pol, L, T, powiu(p,N), p, N));
     735             : }
     736             : 
     737             : static GEN
     738       95236 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
     739             : {
     740       95236 :   long i,l = lg(x);
     741       95236 :   GEN r = cgetg(l,t_COL);
     742      327770 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
     743       95232 :   return r;
     744             : }
     745             : 
     746             : static GEN
     747       95180 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
     748             : {
     749       95180 :   pari_sp av = avma;
     750       95180 :   GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
     751       95180 :   return gerepileupto(av, FqV_roots_from_deg1(y, NULL, q));
     752             : }
     753             : 
     754             : GEN
     755       95040 : ZpX_roots(GEN F, GEN p, long e)
     756             : {
     757       95040 :   pari_sp av = avma;
     758       95040 :   GEN q = powiu(p, e);
     759       95040 :   GEN f = FpX_normalize(F, p);
     760       95040 :   GEN g = FpX_normalize(FpX_split_part(f, p), p);
     761             :   GEN S;
     762       95041 :   if (degpol(g) < degpol(f))
     763             :   {
     764       53719 :     GEN h = FpX_div(f, g, p);
     765       53719 :     F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
     766             :   }
     767       95041 :   S = FpX_roots(g, p);
     768       95041 :   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        6139 : ZqX_roots(GEN F, GEN T, GEN p, long e)
     798             : {
     799        6139 :   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        6240 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
     810             : {
     811        6240 :   pari_sp av = avma;
     812        6240 :   GEN q = p, fr, W;
     813             :   ulong mask;
     814             : 
     815        6240 :   a = modii(a,q);
     816        6240 :   if (e == 1) return a;
     817        6240 :   mask = quadratic_prec_mask(e);
     818        6240 :   fr = FpX_red(f,q);
     819        6240 :   W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
     820             :   for(;;)
     821             :   {
     822       13400 :     q = sqri(q);
     823       13400 :     if (mask & 1) q = diviiexact(q, p);
     824       13400 :     mask >>= 1;
     825       13400 :     fr = FpX_red(f,q);
     826       13400 :     a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
     827       13400 :     if (mask == 1) return gerepileuptoint(av, a);
     828        7160 :     W = Fp_sub(shifti(W,1), Fp_mul(Fp_sqr(W,q), FpX_eval(ZX_deriv(fr),a,q), q), q);
     829             :   }
     830             : }
     831             : 
     832             : GEN
     833         139 : ZpX_liftroots(GEN f, GEN S, GEN p, long e)
     834             : {
     835         139 :   long i, n = lg(S)-1, d = degpol(f);
     836             :   GEN r;
     837         139 :   if (n == d) return ZpX_liftroots_full(f, S, powiu(p, e), p, e);
     838           0 :   r = cgetg(n+1, typ(S));
     839           0 :   for (i=1; i <= n; i++)
     840           0 :     gel(r,i) = ZpX_liftroot(f, gel(S,i), p, e);
     841           0 :   return r;
     842             : }
     843             : 
     844             : GEN
     845         231 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
     846             : {
     847         231 :   pari_sp av = avma, av2;
     848         231 :   GEN pv = p, q, W, df, Tq, fr, dfr;
     849             :   ulong mask;
     850         231 :   a = Fq_red(a, T, p);
     851         231 :   if (e <= v+1) return a;
     852         231 :   df = RgX_deriv(f);
     853         231 :   if (v) { pv = powiu(p,v); df = ZXX_Z_divexact(df, pv); }
     854         231 :   mask = quadratic_prec_mask(e-v);
     855         231 :   Tq = FpXT_red(T, p); dfr = FpXQX_red(df, Tq, p);
     856         231 :   W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
     857         231 :   q = p;
     858         231 :   av2 = avma;
     859             :   for (;;)
     860         161 :   {
     861             :     GEN u, fa, qv, q2v, q2, Tq2;
     862         392 :     q2 = q; q = sqri(q);
     863         392 :     if (mask & 1) q = diviiexact(q,p);
     864         392 :     mask >>= 1;
     865         392 :     if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
     866         392 :     else { qv = q; q2v = q2; }
     867         392 :     Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
     868         392 :     fr = FpXQX_red(f, Tq, qv);
     869         392 :     fa = FqX_eval(fr, a, Tq, qv);
     870         392 :     fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
     871         392 :     a = Fq_sub(a, Fq_mul(Fq_mul(W,fa,Tq2,q2v),q2, Tq,qv), Tq, qv);
     872         392 :     if (mask == 1) return gerepileupto(av, a);
     873         161 :     dfr = FpXQX_red(df, Tq, q);
     874         161 :     u = Fq_sub(Fq_mul(W,FqX_eval(dfr,a,Tq,q),Tq,q),gen_1,Tq,q);
     875         161 :     u = typ(u)==t_INT? diviiexact(u,q2): ZX_Z_divexact(u,q2);
     876         161 :     W = Fq_sub(W, Fq_mul(Fq_mul(u,W,Tq2,q2),q2, Tq,q),Tq,q);
     877         161 :     if (gc_needed(av2,2))
     878             :     {
     879           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_liftroot, e = %ld", e);
     880           0 :       gerepileall(av2, 3, &a, &W, &q);
     881             :     }
     882             :   }
     883             : }
     884             : 
     885             : GEN
     886         231 : ZpXQX_liftroot(GEN f, GEN a, GEN T, GEN p, long e) { return ZpXQX_liftroot_vald(f,a,0,T,p,e); }
     887             : 
     888             : GEN
     889           0 : ZpXQX_liftroots(GEN f, GEN S, GEN T, GEN p, long e)
     890             : {
     891           0 :   long i, n = lg(S)-1, d = degpol(f);
     892             :   GEN r;
     893           0 :   if (n == d) return ZpXQX_liftroots_full(f, S, T, powiu(p, e), p, e);
     894           0 :   r = cgetg(n+1, typ(S));
     895           0 :   for (i=1; i <= n; i++)
     896           0 :     gel(r,i) = ZpXQX_liftroot(f, gel(S,i), T, p, e);
     897           0 :   return r;
     898             : }
     899             : 
     900             : /* Compute (x-1)/(x+1)/p^k */
     901             : static GEN
     902       21071 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
     903             : {
     904       21071 :   pari_sp av = avma;
     905       21071 :   long vT = get_FpX_var(T);
     906             :   GEN bn, bdi;
     907       21071 :   GEN bd = ZX_Z_add(x, gen_1);
     908       21071 :   if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
     909             :   {
     910        7120 :     bn = ZX_shifti(x,-(k+1));
     911        7120 :     bdi= ZpXQ_invlift(ZX_shifti(bd ,-1), pol_1(vT), T, p, e);
     912             :   }
     913             :   else
     914             :   {
     915       13951 :     bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
     916       13951 :     bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
     917             :   }
     918       21071 :   return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
     919             : }
     920             : 
     921             : /* Assume p odd, a = 1 [p], return log(a) */
     922             : GEN
     923       21071 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
     924             : {
     925       21071 :   pari_sp av = avma;
     926             :   pari_timer ti;
     927       21071 :   long is2 = absequaliu(p,2);
     928       21071 :   ulong pp = is2 ? 0: itou_or_0(p);
     929       21071 :   double lp = is2 ? 1: pp ? log2(pp): expi(p);
     930       21071 :   long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
     931             :   GEN ak, s, b, pol;
     932       21071 :   long e = is2 ? N-1: N;
     933       21071 :   long i, l = (e-2)/(2*(k+is2));
     934       21071 :   GEN pe = powiu(p,e);
     935       21071 :   GEN TNk, pNk = powiu(p,N+k);
     936       21071 :   if( DEBUGLEVEL>=3) timer_start(&ti);
     937       21071 :   TNk = FpX_get_red(get_FpX_mod(T), pNk);
     938       21071 :   ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
     939       21071 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
     940       21071 :   b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
     941       21071 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
     942       21071 :   pol= cgetg(l+3,t_POL);
     943       21071 :   pol[1] = evalsigne(1)|evalvarn(0);
     944       57946 :   for(i=0; i<=l; i++)
     945             :   {
     946             :     GEN g;
     947       36876 :     ulong z = 2*i+1;
     948       36876 :     if (pp)
     949             :     {
     950       26243 :       long w = u_lvalrem(z, pp, &z);
     951       26243 :       g = powuu(pp,2*i*k-w);
     952             :     }
     953       10633 :     else g = powiu(p,2*i*k);
     954       36875 :     gel(pol,i+2) = Fp_divu(g, z,pe);
     955             :   }
     956       21070 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
     957       21070 :   s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T,  pe);
     958       21071 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
     959       21071 :   s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
     960       21071 :   return gerepileupto(av, is2? s: FpX_red(s, pe));
     961             : }
     962             : 
     963             : /***********************************************************************/
     964             : /**                                                                   **/
     965             : /**                 Generic quadratic hensel lift over Zp[X]          **/
     966             : /**                                                                   **/
     967             : /***********************************************************************/
     968             : /* q = p^N */
     969             : GEN
     970           0 : gen_ZpM_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
     971             :                             GEN lin(void *E, GEN F, GEN d, GEN q),
     972             :                             GEN invl(void *E, GEN d))
     973             : {
     974           0 :   pari_sp av = avma;
     975             :   long N2, M;
     976             :   GEN VN2, V2, VM, bil;
     977             :   GEN q2, qM;
     978           0 :   V = FpM_red(V, q);
     979           0 :   if (N == 1) return invl(E, V);
     980           0 :   N2 = (N + 1)>>1; M = N - N2;
     981           0 :   F = FpM_red(F, q);
     982           0 :   qM = powiu(p, M);
     983           0 :   q2 = M == N2? qM: mulii(qM, p);
     984             :   /* q2 = p^N2, qM = p^M, q = q2 * qM */
     985           0 :   VN2 = gen_ZpM_Dixon(F, V, q2, p, N2, E, lin, invl);
     986           0 :   bil = lin(E, F, VN2, q);
     987           0 :   V2 = ZM_Z_divexact(ZM_sub(V, bil), q2);
     988           0 :   VM = gen_ZpM_Dixon(F, V2, qM, p, M, E, lin, invl);
     989           0 :   return gerepileupto(av, FpM_red(ZM_add(VN2, ZM_Z_mul(VM, q2)), q));
     990             : }
     991             : 
     992             : GEN
     993      199955 : gen_ZpX_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
     994             :                             GEN lin(void *E, GEN F, GEN d, GEN q),
     995             :                             GEN invl(void *E, GEN d))
     996             : {
     997      199955 :   pari_sp av = avma;
     998             :   long N2, M;
     999             :   GEN VN2, V2, VM, bil;
    1000             :   GEN q2, qM;
    1001      199955 :   V = FpX_red(V, q);
    1002      199955 :   if (N == 1) return invl(E, V);
    1003       45458 :   N2 = (N + 1)>>1; M = N - N2;
    1004       45458 :   F = FpXT_red(F, q);
    1005       45458 :   qM = powiu(p, M);
    1006       45458 :   q2 = M == N2? qM: mulii(qM, p);
    1007             :   /* q2 = p^N2, qM = p^M, q = q2 * qM */
    1008       45458 :   VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
    1009       45458 :   bil = lin(E, F, VN2, q);
    1010       45458 :   V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
    1011       45458 :   VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
    1012       45458 :   return gerepileupto(av, FpX_red(ZX_add(VN2, ZX_Z_mul(VM, q2)), q));
    1013             : }
    1014             : 
    1015             : GEN
    1016           0 : gen_ZpM_Newton(GEN x, GEN p, long n, void *E,
    1017             :                       GEN eval(void *E, GEN f, GEN q),
    1018             :                       GEN invd(void *E, GEN V, GEN v, GEN q, long M))
    1019             : {
    1020           0 :   pari_sp ltop = avma, av;
    1021           0 :   long N = 1, N2, M;
    1022             :   long mask;
    1023           0 :   GEN q = p;
    1024           0 :   if (n == 1) return gcopy(x);
    1025           0 :   mask = quadratic_prec_mask(n);
    1026           0 :   av = avma;
    1027           0 :   while (mask > 1)
    1028             :   {
    1029             :     GEN qM, q2, v, V;
    1030           0 :     N2 = N; N <<= 1;
    1031           0 :     q2 = q;
    1032           0 :     if (mask&1UL) { /* can never happen when q2 = p */
    1033           0 :       N--; M = N2-1;
    1034           0 :       qM = diviiexact(q2,p); /* > 1 */
    1035           0 :       q = mulii(qM,q2);
    1036             :     } else {
    1037           0 :       M = N2;
    1038           0 :       qM = q2;
    1039           0 :       q = sqri(q2);
    1040             :     }
    1041             :     /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
    1042           0 :     mask >>= 1;
    1043           0 :     v = eval(E, x, q);
    1044           0 :     V = ZM_Z_divexact(gel(v,1), q2);
    1045           0 :     x = FpM_sub(x, ZM_Z_mul(invd(E, V, v, qM, M), q2), q);
    1046           0 :     if (gc_needed(av, 1))
    1047             :     {
    1048           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpM_Newton");
    1049           0 :       gerepileall(av, 2, &x, &q);
    1050             :     }
    1051             :   }
    1052           0 :   return gerepileupto(ltop, x);
    1053             : }
    1054             : 
    1055             : static GEN
    1056           0 : _ZpM_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
    1057             : {
    1058             :   (void)E; (void)M;
    1059           0 :   return FpM_mul(V, gel(v,2), q);
    1060             : }
    1061             : 
    1062             : static GEN
    1063           0 : _ZpM_eval(void *E, GEN x, GEN q)
    1064             : {
    1065           0 :   GEN f = RgM_Rg_sub_shallow(FpM_mul(x, FpM_red((GEN) E, q), q), gen_1);
    1066           0 :   return mkvec2(f, x);
    1067             : }
    1068             : 
    1069             : GEN
    1070           0 : ZpM_invlift(GEN M, GEN C, GEN p, long n)
    1071             : {
    1072           0 :   return gen_ZpM_Newton(C, p, n, (void *)M, _ZpM_eval, _ZpM_invd);
    1073             : }
    1074             : 
    1075             : GEN
    1076      267067 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
    1077             :                       GEN eval(void *E, GEN f, GEN q),
    1078             :                       GEN invd(void *E, GEN V, GEN v, GEN q, long M))
    1079             : {
    1080      267067 :   pari_sp ltop = avma, av;
    1081      267067 :   long N = 1, N2, M;
    1082             :   long mask;
    1083      267067 :   GEN q = p;
    1084      267067 :   if (n == 1) return gcopy(x);
    1085      263238 :   mask = quadratic_prec_mask(n);
    1086      263238 :   av = avma;
    1087      930202 :   while (mask > 1)
    1088             :   {
    1089             :     GEN qM, q2, v, V;
    1090      666964 :     N2 = N; N <<= 1;
    1091      666964 :     q2 = q;
    1092      666964 :     if (mask&1UL) { /* can never happen when q2 = p */
    1093      276561 :       N--; M = N2-1;
    1094      276561 :       qM = diviiexact(q2,p); /* > 1 */
    1095      276558 :       q = mulii(qM,q2);
    1096             :     } else {
    1097      390403 :       M = N2;
    1098      390403 :       qM = q2;
    1099      390403 :       q = sqri(q2);
    1100             :     }
    1101             :     /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
    1102      666962 :     mask >>= 1;
    1103      666962 :     v = eval(E, x, q);
    1104      666965 :     V = ZX_Z_divexact(gel(v,1), q2);
    1105      666958 :     x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
    1106      666964 :     if (gc_needed(av, 1))
    1107             :     {
    1108           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
    1109           0 :       gerepileall(av, 2, &x, &q);
    1110             :     }
    1111             :   }
    1112      263238 :   return gerepileupto(ltop, x);
    1113             : }
    1114             : 
    1115             : struct _ZpXQ_inv
    1116             : {
    1117             :   GEN T, a, p ,n;
    1118             : };
    1119             : 
    1120             : static GEN
    1121      523116 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
    1122             : {
    1123      523116 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
    1124      523116 :   GEN Tq = FpXT_red(d->T, q);
    1125             :   (void)M;
    1126      523117 :   return FpXQ_mul(V, gel(v,2), Tq, q);
    1127             : }
    1128             : 
    1129             : static GEN
    1130      523116 : _inv_eval(void *E, GEN x, GEN q)
    1131             : {
    1132      523116 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
    1133      523116 :   GEN Tq = FpXT_red(d->T, q);
    1134      523116 :   GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
    1135      523117 :   return mkvec2(f, x);
    1136             : }
    1137             : 
    1138             : GEN
    1139      202049 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
    1140             : {
    1141             :   struct _ZpXQ_inv d;
    1142      202049 :   d.a = a; d.T = T; d.p = p;
    1143      202049 :   return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
    1144             : }
    1145             : 
    1146             : GEN
    1147      180978 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
    1148             : {
    1149      180978 :   pari_sp av=avma;
    1150             :   GEN ai;
    1151      180978 :   if (lgefint(p)==3)
    1152             :   {
    1153      180950 :     ulong pp = p[2];
    1154      180950 :     ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
    1155             :   } else
    1156          28 :     ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
    1157      180978 :   return gerepileupto(av, ZpXQ_invlift(a, ai, T, p, e));
    1158             : }
    1159             : 
    1160             : GEN
    1161       31605 : ZpXQ_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
    1162             : {
    1163       31605 :   return FpXQ_mul(a, ZpXQ_inv(b, T, p, e), T, q);
    1164             : }
    1165             : 
    1166             : GEN
    1167      146440 : ZpXQX_divrem(GEN x, GEN Sp, GEN T, GEN q, GEN p, long e, GEN *pr)
    1168             : {
    1169      146440 :   pari_sp av = avma;
    1170      146440 :   GEN S = get_FpXQX_mod(Sp);
    1171      146440 :   GEN b = leading_coeff(S), bi;
    1172             :   GEN S2, Q;
    1173      146440 :   if (typ(b)==t_INT) return FpXQX_divrem(x, Sp, T, q, pr);
    1174      146440 :   bi = ZpXQ_inv(b, T, p, e);
    1175      146440 :   S2 = FqX_Fq_mul_to_monic(S, bi, T, q);
    1176      146440 :   Q = FpXQX_divrem(x, S2, T, q, pr);
    1177      146440 :   if (pr==ONLY_DIVIDES && !Q) { set_avma(av); return NULL; }
    1178      146440 :   if (pr==ONLY_REM || pr==ONLY_DIVIDES) return gerepileupto(av, Q);
    1179      146440 :   Q = FpXQX_FpXQ_mul(Q, bi, T, q);
    1180      146440 :   return gc_all(av, 2, &Q, pr);
    1181             : }
    1182             : 
    1183             : GEN
    1184         623 : ZpXQX_digits(GEN x, GEN B, GEN T, GEN q, GEN p, long e)
    1185             : {
    1186         623 :   pari_sp av = avma;
    1187         623 :   GEN b = leading_coeff(B), bi;
    1188             :   GEN B2, P, V, W;
    1189             :   long i, lV;
    1190         623 :   if (typ(b)==t_INT) return FpXQX_digits(x, B, T, q);
    1191         623 :   bi = ZpXQ_inv(b, T, p, e);
    1192         623 :   B2 = FqX_Fq_mul_to_monic(B, bi, T, q);
    1193         623 :   V = FpXQX_digits(x, B2, T, q);
    1194         623 :   lV = lg(V)-1;
    1195         623 :   P = FpXQ_powers(bi, lV-1, T, q);
    1196         623 :   W = cgetg(lV+1, t_VEC);
    1197       11200 :   for(i=1; i<=lV; i++)
    1198       10577 :     gel(W, i) = FpXQX_FpXQ_mul(gel(V,i), gel(P, i), T, q);
    1199         623 :   return gerepileupto(av, W);
    1200             : }
    1201             : 
    1202             : struct _ZpXQ_sqrtn
    1203             : {
    1204             :   GEN T, a, n, ai;
    1205             : };
    1206             : 
    1207             : static GEN
    1208        1806 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
    1209             : {
    1210        1806 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
    1211        1806 :   GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
    1212             :   (void)M;
    1213        1806 :   return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
    1214             : }
    1215             : 
    1216             : static GEN
    1217        1806 : _sqrtn_eval(void *E, GEN x, GEN q)
    1218             : {
    1219        1806 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
    1220        1806 :   GEN Tq = FpX_red(d->T, q);
    1221        1806 :   GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
    1222        1806 :   return mkvec2(f, x);
    1223             : }
    1224             : 
    1225             : GEN
    1226        1155 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
    1227             : {
    1228             :   struct _ZpXQ_sqrtn d;
    1229        1155 :   d.a = a; d.T = T; d.n = n;
    1230        1155 :   d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
    1231        1155 :   return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
    1232             : }
    1233             : 
    1234             : static GEN
    1235           0 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol_shallow(a,v): a; }
    1236             : 
    1237             : GEN
    1238           0 : Zq_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
    1239             : {
    1240           0 :   return T? ZpXQ_sqrtnlift(to_ZX(a,varn(T)), n, to_ZX(x,varn(T)), T, p, e)
    1241           0 :           : Zp_sqrtnlift(a, n, x, p, e);
    1242             : }
    1243             : 
    1244             : GEN
    1245           0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
    1246             : {
    1247           0 :   pari_sp av = avma;
    1248           0 :   GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
    1249           0 :   if (!z) return NULL;
    1250           0 :   if (e <= 1) return gerepileupto(av, z);
    1251           0 :   return gerepileupto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
    1252             : }
    1253             : 
    1254             : GEN
    1255       33271 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
    1256             :                      GEN early(void *E, GEN x, GEN q))
    1257             : {
    1258       33271 :   pari_sp ltop = avma, av;
    1259             :   long N, r;
    1260             :   long mask;
    1261             :   GEN q2, q, W, Q, Tq2, Tq, Pq;
    1262             :   pari_timer ti;
    1263       33271 :   T = FpX_get_red(T, powiu(p, n));
    1264       33268 :   if (n == 1) return gcopy(S);
    1265       33268 :   mask = quadratic_prec_mask(n);
    1266       33269 :   av = avma;
    1267       33269 :   q2 = p; q = sqri(p); mask >>= 1; N = 2;
    1268       33268 :   if (DEBUGLEVEL > 3) timer_start(&ti);
    1269       33268 :   Tq = FpXT_red(T,q);
    1270       33270 :   Tq2 = FpXT_red(Tq,q2);
    1271       33268 :   Pq = FpX_red(P,q);
    1272       33270 :   W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
    1273       33270 :   Q  = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
    1274       33270 :   r = brent_kung_optpow(degpol(P), 4, 3);
    1275       33270 :   if (DEBUGLEVEL > 3)
    1276           0 :     err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
    1277             :   for (;;)
    1278       38430 :   {
    1279             :     GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
    1280       71700 :     H  = FpXQ_mul(W, Q, Tq2, q2);
    1281       71701 :     Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
    1282       71700 :     if (DEBUGLEVEL > 3)
    1283           0 :       timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
    1284       71700 :     if (mask==1)
    1285        6399 :       return gerepileupto(ltop, Sq);
    1286       65301 :     if (early)
    1287             :     {
    1288       63642 :       GEN Se = early(E, Sq, q);
    1289       63643 :       if (Se) return gerepileupto(ltop, Se);
    1290             :     }
    1291       38430 :     qq = sqri(q); N <<= 1;
    1292       38431 :     if (mask&1UL) { qq = diviiexact(qq, p); N--; }
    1293       38431 :     mask >>= 1;
    1294       38431 :     Pqq  = FpX_red(P, qq);
    1295       38431 :     Tqq  = FpXT_red(T, qq);
    1296       38431 :     Spow = FpXQ_powers(Sq, r, Tqq, qq);
    1297       38431 :     Q  = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
    1298       38430 :     dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
    1299       38430 :     Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
    1300       38429 :     Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
    1301       38430 :     Wq = FpX_sub(W, Wq, q);
    1302       38430 :     S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
    1303       38430 :     if (gc_needed(av, 1))
    1304             :     {
    1305           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_liftroot");
    1306           0 :       gerepileall(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
    1307             :     }
    1308             :   }
    1309             : }
    1310             : 
    1311             : GEN
    1312        1729 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
    1313             : {
    1314        1729 :   return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
    1315             : }
    1316             : 
    1317             : GEN
    1318         301 : ZpX_Frobenius(GEN T, GEN p, long e)
    1319             : {
    1320         301 :   return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
    1321             : }
    1322             : 
    1323             : GEN
    1324         147 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
    1325             : {
    1326         147 :   pari_sp av = avma;
    1327         147 :   GEN xp = ZpX_Frobenius(T, p, e);
    1328         147 :   GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
    1329         147 :   return gerepilecopy(av, gel(z,2));
    1330             : }
    1331             : 
    1332             : GEN
    1333           0 : ZpXQX_ZpXQXQ_liftroot(GEN P, GEN S, GEN U, GEN T, GEN p, long n)
    1334             : {
    1335           0 :   pari_sp ltop = avma, av;
    1336             :   long N, r;
    1337             :   long mask;
    1338             :   GEN  qn, q2, q, W, Q, Tq2, Tq, Pq, Uq, Uq2;
    1339             :   pari_timer ti;
    1340           0 :   qn = powiu(p, n);
    1341           0 :   T = FpX_get_red(T, qn);
    1342           0 :   U = FpXQX_get_red(U, T, qn);
    1343           0 :   if (n == 1) return gcopy(S);
    1344           0 :   mask = quadratic_prec_mask(n);
    1345           0 :   av = avma;
    1346           0 :   q2 = p; q = sqri(p); mask >>= 1; N = 2;
    1347           0 :   if (DEBUGLEVEL > 3) timer_start(&ti);
    1348           0 :   Tq = FpXT_red(T,q);
    1349           0 :   Uq = FpXQXT_red(U, Tq, q);
    1350           0 :   Tq2 = FpXT_red(Tq, q2);
    1351           0 :   Uq2 = FpXQXT_red(U, Tq2, q2);
    1352           0 :   Pq = FpXQX_red(P, Tq, q);
    1353           0 :   W = FpXQXQ_inv(FpXQX_FpXQXQ_eval(FpXX_deriv(P,q2), S, Uq2, Tq2, q2), Uq2, Tq2, q2);
    1354           0 :   Q  = ZXX_Z_divexact(FpXQX_FpXQXQ_eval(Pq, S, Uq, Tq, q), q2);
    1355           0 :   r = brent_kung_optpow(degpol(P), 4, 3);
    1356           0 :   if (DEBUGLEVEL > 3)
    1357           0 :     err_printf("ZpXQX_ZpXQXQ_liftroot: lifting to prec %ld\n",n);
    1358             :   for (;;)
    1359           0 :   {
    1360             :     GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq, Uqq;
    1361           0 :     H  = FpXQXQ_mul(W, Q, Uq2, Tq2, q2);
    1362           0 :     Sq = FpXX_sub(S, ZXX_Z_mul(H, q2), q);
    1363           0 :     if (DEBUGLEVEL > 3)
    1364           0 :       timer_printf(&ti,"ZpXQX_ZpXQXQ_liftroot: reaching prec %ld",N);
    1365           0 :     if (mask==1)
    1366           0 :       return gerepileupto(ltop, Sq);
    1367           0 :     qq = sqri(q); N <<= 1;
    1368           0 :     if (mask&1UL) { qq = diviiexact(qq, p); N--; }
    1369           0 :     mask >>= 1;
    1370           0 :     Tqq  = FpXT_red(T, qq);
    1371           0 :     Uqq  = FpXQXT_red(U, Tqq, qq);
    1372           0 :     Pqq  = FpXQX_red(P, Tqq, qq);
    1373           0 :     Spow = FpXQXQ_powers(Sq, r, Uqq, Tqq, qq);
    1374           0 :     Q  = ZXX_Z_divexact(FpXQX_FpXQXQV_eval(Pqq, Spow, Uqq, Tqq, qq), q);
    1375           0 :     dP = FpXQX_FpXQXQV_eval(FpXX_deriv(Pq, q), FpXQXV_red(Spow, Tq, q), Uq, Tq, q);
    1376           0 :     Wq = ZXX_Z_divexact(gsub(FpXQXQ_mul(W, dP, Uq, Tq, q), gen_1), q2);
    1377           0 :     Wq = ZXX_Z_mul(FpXQXQ_mul(W, Wq, Uq2, Tq2, q2), q2);
    1378           0 :     Wq = FpXX_sub(W, Wq, q);
    1379           0 :     S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Uq2 = Uq; Uq = Uqq;  Pq = Pqq;
    1380           0 :     if (gc_needed(av, 1))
    1381             :     {
    1382           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_ZpXQXQ_liftroot");
    1383           0 :       gerepileall(av, 10, &S, &W, &Q, &Uq2, &Uq, &Tq2, &Tq, &Pq, &q, &q2);
    1384             :     }
    1385             :   }
    1386             : }
    1387             : 
    1388             : GEN
    1389          35 : ZqX_ZqXQ_liftroot(GEN f, GEN a, GEN P, GEN T, GEN p, long e)
    1390          35 : { return T ? ZpXQX_ZpXQXQ_liftroot(f, a, P, T , p, e): ZpX_ZpXQ_liftroot(f, a, P, p, e); }
    1391             : 
    1392             : /* Canonical lift of polynomial */
    1393             : 
    1394       11130 : static GEN _can_invl(void *E, GEN V) {(void) E; return V; }
    1395             : 
    1396        3654 : static GEN _can_lin(void *E, GEN F, GEN V, GEN q)
    1397             : {
    1398        3654 :   GEN v = RgX_splitting(V, 3);
    1399             :   (void) E;
    1400        3654 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
    1401             : }
    1402             : 
    1403             : static GEN
    1404        7476 : _can_iter(void *E, GEN f, GEN q)
    1405             : {
    1406        7476 :   GEN h = RgX_splitting(f,3);
    1407        7476 :   GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
    1408        7476 :   GEN h12 = ZX_mul(gel(h,1), gel(h,2));
    1409        7476 :   GEN h13 = ZX_mul(gel(h,1), gel(h,3));
    1410        7476 :   GEN h23 = ZX_mul(gel(h,2), gel(h,3));
    1411        7476 :   GEN h1c = ZX_mul(gel(h,1), h1s);
    1412        7476 :   GEN h3c = ZX_mul(gel(h,3), h3s);
    1413        7476 :   GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
    1414        7476 :   GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
    1415             :   (void) E;
    1416        7476 :   return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
    1417             : }
    1418             : 
    1419             : static GEN
    1420        7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
    1421             : {
    1422        7476 :   GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
    1423        7476 :   GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
    1424        7476 :   GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
    1425             :                  ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
    1426             :   (void)E;
    1427        7476 :   return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
    1428             :                                                  _can_lin, _can_invl);
    1429             : }
    1430             : 
    1431             : static GEN
    1432        3717 : F3x_frobeniuslift(GEN P, long n)
    1433        3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
    1434             : 
    1435       29778 : static GEN _can5_invl(void *E, GEN V) {(void) E; return V; }
    1436             : 
    1437        9107 : static GEN _can5_lin(void *E, GEN F, GEN V, GEN q)
    1438             : {
    1439        9107 :   ulong p = *(ulong*)E;
    1440        9107 :   GEN v = RgX_splitting(V, p);
    1441        9107 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
    1442             : }
    1443             : 
    1444             : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
    1445             : static GEN
    1446       62314 : _shift(GEN P, long n, ulong p, long v)
    1447             : {
    1448       62314 :   long i, l=lg(P);
    1449       62314 :   GEN r = cgetg(l,t_POL); r[1] = P[1];
    1450      484498 :   for(i=2;i<l;i++)
    1451             :   {
    1452      422184 :     long s = n*(i-2)%p;
    1453      422184 :     GEN ci = gel(P,i);
    1454      422184 :     if (typ(ci)==t_INT)
    1455      105147 :       gel(r,i) = monomial(ci, s, v);
    1456             :     else
    1457      317037 :       gel(r,i) = RgX_rotate_shallow(ci, s, p);
    1458             :   }
    1459       62314 :   return FpXX_renormalize(r, l);
    1460             : }
    1461             : 
    1462             : struct _can_mul
    1463             : {
    1464             :   GEN T, q;
    1465             :   ulong p;
    1466             : };
    1467             : 
    1468             : static GEN
    1469       41643 : _can5_mul(void *E, GEN A, GEN B)
    1470             : {
    1471       41643 :   struct _can_mul *d = (struct _can_mul *)E;
    1472       41643 :   GEN a = gel(A,1), b = gel(B,1);
    1473       41643 :   long n = itos(gel(A,2));
    1474       41643 :   GEN bn = _shift(b, n, d->p, get_FpX_var(d->T));
    1475       41643 :   GEN c = FpXQX_mul(a, bn, d->T, d->q);
    1476       41643 :   return mkvec2(c, addii(gel(A,2), gel(B,2)));
    1477             : }
    1478             : 
    1479             : static GEN
    1480       41447 : _can5_sqr(void *E, GEN A)
    1481             : {
    1482       41447 :   return _can5_mul(E,A,A);
    1483             : }
    1484             : 
    1485             : static GEN
    1486       20671 : _can5_iter(void *E, GEN f, GEN q)
    1487             : {
    1488       20671 :   pari_sp av = avma;
    1489             :   struct _can_mul D;
    1490       20671 :   ulong p = *(ulong*)E;
    1491       20671 :   long i, vT = fetch_var();
    1492             :   GEN N, P, d, V, fs;
    1493       20671 :   D.q = q; D.T = ZX_Z_sub(pol_xn(p,vT),gen_1);
    1494       20671 :   D.p = p;
    1495       20671 :   fs = mkvec2(_shift(f, 1, p, vT), gen_1);
    1496       20671 :   N = gel(gen_powu_i(fs,p-1,(void*)&D,_can5_sqr,_can5_mul),1);
    1497       20671 :   N = ZXX_evalx0(FpXQX_red(N,polcyclo(p,vT),q));
    1498       20671 :   P = FpX_mul(N,f,q);
    1499       20671 :   P = RgX_deflate(P, p);
    1500       20671 :   d = RgX_splitting(N, p);
    1501       20671 :   V = cgetg(p+1,t_VEC);
    1502       20671 :   gel(V,1) = ZX_mulu(gel(d,1), p);
    1503      104265 :   for(i=2; i<= (long)p; i++)
    1504       83594 :     gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
    1505       20671 :   (void)delete_var(); return gerepilecopy(av, mkvec2(ZX_sub(f,P),V));
    1506             : }
    1507             : 
    1508             : static GEN
    1509       20671 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
    1510             : {
    1511       20671 :   ulong p = *(long*)E;
    1512       20671 :   return gen_ZpX_Dixon(gel(v,2), H, qM, utoipos(p), M, E, _can5_lin, _can5_invl);
    1513             : }
    1514             : 
    1515             : GEN
    1516       13979 : Flx_Teichmuller(GEN P, ulong p, long n)
    1517             : {
    1518       24241 :   return p==3 ? F3x_frobeniuslift(P,n):
    1519       10262 :          gen_ZpX_Newton(Flx_to_ZX(P),utoipos(p), n, &p, _can5_iter, _can5_invd);
    1520             : }
    1521             : 
    1522             : GEN
    1523          35 : polteichmuller(GEN P, ulong p, long n)
    1524             : {
    1525          35 :   pari_sp av = avma;
    1526          35 :   GEN q = NULL;
    1527          35 :   if (typ(P)!=t_POL || !RgX_is_FpX(P,&q)) pari_err_TYPE("polteichmuller",P);
    1528          35 :   if (q && !equaliu(q,p)) pari_err_MODULUS("polteichmuller",q,utoi(p));
    1529          35 :   if (n <= 0)
    1530           0 :     pari_err_DOMAIN("polteichmuller", "precision", "<=",gen_0,stoi(n));
    1531          63 :   return gerepileupto(av, p==2 ? F2x_Teichmuller(RgX_to_F2x(P), n)
    1532          28 :                                : Flx_Teichmuller(RgX_to_Flx(P, p), p, n));
    1533             : }

Generated by: LCOV version 1.16