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 - Hensel.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24961-b01324313d) Lines: 636 683 93.1 %
Date: 2020-01-21 05:55:55 Functions: 65 68 95.6 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : 
      16             : /***********************************************************************/
      17             : /**                                                                   **/
      18             : /**       QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL)         **/
      19             : /**                                                                   **/
      20             : /***********************************************************************/
      21             : /* Setup for divide/conquer quadratic Hensel lift
      22             :  *  a = set of k t_POL in Z[X] = factors over Fp (T=NULL) or Fp[Y]/(T)
      23             :  *  V = set of products of factors built as follows
      24             :  *  1) V[1..k] = initial a
      25             :  *  2) iterate:
      26             :  *      append to V the two smallest factors (minimal degree) in a, remove them
      27             :  *      from a and replace them by their product [net loss for a = 1 factor]
      28             :  *
      29             :  * W = bezout coeffs W[i]V[i] + W[i+1]V[i+1] = 1
      30             :  *
      31             :  * link[i] = -j if V[i] = a[j]
      32             :  *            j if V[i] = V[j] * V[j+1]
      33             :  * Arrays (link, V, W) pre-allocated for 2k - 2 elements */
      34             : static void
      35       59306 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
      36             : {
      37       59306 :   long k = lg(a)-1;
      38             :   long i, j, s, minp, mind;
      39             : 
      40       59306 :   for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
      41             : 
      42       98433 :   for (j=1; j <= 2*k-5; j+=2,i++)
      43             :   {
      44       39127 :     minp = j;
      45       39127 :     mind = degpol(gel(V,j));
      46      484318 :     for (s=j+1; s<i; s++)
      47      445191 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
      48             : 
      49       39127 :     swap(gel(V,j), gel(V,minp));
      50       39127 :     lswap(link[j], link[minp]);
      51             : 
      52       39127 :     minp = j+1;
      53       39127 :     mind = degpol(gel(V,j+1));
      54      445191 :     for (s=j+2; s<i; s++)
      55      406064 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
      56             : 
      57       39127 :     swap(gel(V,j+1), gel(V,minp));
      58       39127 :     lswap(link[j+1], link[minp]);
      59             : 
      60       39127 :     gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
      61       39127 :     link[i] = j;
      62             :   }
      63             : 
      64      157732 :   for (j=1; j <= 2*k-3; j+=2)
      65             :   {
      66             :     GEN d, u, v;
      67       98433 :     d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
      68       98433 :     if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
      69       98426 :     d = gel(d,2);
      70       98426 :     if (!gequal1(d))
      71             :     {
      72       89100 :       if (typ(d)==t_POL)
      73             :       {
      74        8309 :         d = FpXQ_inv(d, T, p);
      75        8309 :         u = FqX_Fq_mul(u, d, T, p);
      76        8309 :         v = FqX_Fq_mul(v, d, T, p);
      77             :       }
      78             :       else
      79             :       {
      80       80791 :         d = Fp_inv(d, p);
      81       80791 :         u = FqX_Fp_mul(u, d, T,p);
      82       80791 :         v = FqX_Fp_mul(v, d, T,p);
      83             :       }
      84             :     }
      85       98426 :     gel(W,j) = u;
      86       98426 :     gel(W,j+1) = v;
      87             :   }
      88       59299 : }
      89             : 
      90             : /* au + bv = 1 (p0), ab = f (p0). Lift mod p1 = p0 pd (<= p0^2).
      91             :  * If noinv is set, don't lift the inverses u and v */
      92             : static void
      93      412103 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
      94             : {
      95      412103 :   pari_sp av = avma;
      96      412103 :   long space = lg(f) * lgefint(p1);
      97             :   GEN a2, b2, g, z, s, t;
      98      412103 :   GEN a = gel(V,j), b = gel(V,j+1);
      99      412103 :   GEN u = gel(W,j), v = gel(W,j+1);
     100             : 
     101      412103 :   (void)new_chunk(space); /* HACK */
     102      412103 :   g = ZX_sub(f, ZX_mul(a,b));
     103      412103 :   g = ZX_Z_divexact(g, p0);
     104      412103 :   g = FpX_red(g, pd);
     105      412103 :   z = FpX_mul(v,g, pd);
     106      412103 :   t = FpX_divrem(z,a, pd, &s);
     107      412103 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     108      412103 :   t = FpX_red(t, pd);
     109      412103 :   t = ZX_Z_mul(t,p0);
     110      412103 :   s = ZX_Z_mul(s,p0);
     111      412103 :   set_avma(av);
     112      412103 :   a2 = ZX_add(a,s);
     113      412103 :   b2 = ZX_add(b,t);
     114             : 
     115             :   /* already reduced mod p1 = pd p0 */
     116      412103 :   gel(V,j)   = a2;
     117      412103 :   gel(V,j+1) = b2;
     118      412103 :   if (noinv) return;
     119             : 
     120      323421 :   av = avma;
     121      323421 :   (void)new_chunk(space); /* HACK */
     122      323421 :   g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
     123      323421 :   g = Z_ZX_sub(gen_1, g);
     124      323421 :   g = ZX_Z_divexact(g, p0);
     125      323421 :   g = FpX_red(g, pd);
     126      323421 :   z = FpX_mul(v,g, pd);
     127      323421 :   t = FpX_divrem(z,a, pd, &s);
     128      323421 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     129      323421 :   t = FpX_red(t, pd);
     130      323421 :   t = ZX_Z_mul(t,p0);
     131      323421 :   s = ZX_Z_mul(s,p0);
     132      323421 :   set_avma(av);
     133      323421 :   gel(W,j)   = ZX_add(u,t);
     134      323421 :   gel(W,j+1) = ZX_add(v,s);
     135             : }
     136             : 
     137             : static void
     138       36757 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
     139             : {
     140       36757 :   pari_sp av = avma;
     141       36757 :   const long n = degpol(T1), vT = varn(T1);
     142       36757 :   long space = lg(f) * lgefint(p1) * lg(T1);
     143             :   GEN a2, b2, g, z, s, t;
     144       36757 :   GEN a = gel(V,j), b = gel(V,j+1);
     145       36757 :   GEN u = gel(W,j), v = gel(W,j+1);
     146             : 
     147       36757 :   (void)new_chunk(space); /* HACK */
     148       36757 :   g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
     149       36757 :   g = FpXQX_red(g, T1, p1);
     150       36757 :   g = RgX_Rg_divexact(g, p0);
     151       36757 :   z = FpXQX_mul(v,g, Td,pd);
     152       36757 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     153       36757 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     154       36757 :   t = Kronecker_to_ZXX(t, n, vT);
     155       36757 :   t = FpXQX_red(t, Td, pd);
     156       36757 :   t = RgX_Rg_mul(t,p0);
     157       36757 :   s = RgX_Rg_mul(s,p0);
     158       36757 :   set_avma(av);
     159             : 
     160       36757 :   a2 = RgX_add(a,s);
     161       36757 :   b2 = RgX_add(b,t);
     162             :   /* already reduced mod p1 = pd p0 */
     163       36757 :   gel(V,j)   = a2;
     164       36757 :   gel(V,j+1) = b2;
     165       36757 :   if (noinv) return;
     166             : 
     167       27860 :   av = avma;
     168       27860 :   (void)new_chunk(space); /* HACK */
     169       27860 :   g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
     170       27860 :   g = Kronecker_to_ZXX(g, n, vT);
     171       27860 :   g = Rg_RgX_sub(gen_1, g);
     172       27860 :   g = FpXQX_red(g, T1, p1);
     173       27860 :   g = RgX_Rg_divexact(g, p0);
     174       27860 :   z = FpXQX_mul(v,g, Td,pd);
     175       27860 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     176       27860 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     177       27860 :   t = Kronecker_to_ZXX(t, n, vT);
     178       27860 :   t = FpXQX_red(t, Td, pd);
     179       27860 :   t = RgX_Rg_mul(t,p0);
     180       27860 :   s = RgX_Rg_mul(s,p0);
     181       27860 :   set_avma(av);
     182       27860 :   gel(W,j)   = RgX_add(u,t);
     183       27860 :   gel(W,j+1) = RgX_add(v,s);
     184             : }
     185             : 
     186             : /* v list of factors, w list of inverses.  f = v[j] v[j+1]
     187             :  * Lift v[j] and v[j+1] mod p0 pd (possibly mod T), then all their divisors */
     188             : static void
     189     1061084 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
     190             :                 GEN f, long j, int noinv)
     191             : {
     192     1061084 :   if (j < 0) return;
     193      412103 :   ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
     194      412103 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     195      412103 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
     196             : }
     197             : static void
     198       77014 : ZpXQ_RecTreeLift(GEN link, GEN v, GEN w, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1,
     199             :                  GEN f, long j, int noinv)
     200             : {
     201       77014 :   if (j < 0) return;
     202       36757 :   ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
     203       36757 :   ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     204       36757 :   ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
     205             : }
     206             : 
     207             : /* Assume n > 0. We want to go to accuracy n, starting from accuracy 1, using
     208             :  * a quadratically convergent algorithm. Goal: 9 -> 1,2,3,5,9 instead of
     209             :  * 1,2,4,8,9 (sequence of accuracies).
     210             :  *
     211             :  * Let a0 = 1, a1 = 2, a2, ... ak = n, the sequence of accuracies. To obtain
     212             :  * it, work backwards:
     213             :  *   a(k) = n, a(i-1) = (a(i) + 1) \ 2,
     214             :  * but we do not want to store a(i) explicitly, even as a t_VECSMALL, since
     215             :  * this would leave an object on the stack. We store a(i) implicitly in a
     216             :  * MASK: let a(0) = 1, if the i-bit of MASK is set, set a(i+1) = 2 a(i) - 1,
     217             :  * and 2a(i) otherwise.
     218             :  *
     219             :  * In fact, we do something a little more complicated to simplify the
     220             :  * function interface and avoid returning k and MASK separately: we return
     221             :  * MASK + 2^(k+1), so the highest bit of the mask indicates the length of the
     222             :  * sequence, and the following ones are as above. */
     223             : ulong
     224     2716746 : quadratic_prec_mask(long n)
     225             : {
     226     2716746 :   long a = n, i;
     227     2716746 :   ulong mask = 0;
     228     6565900 :   for(i = 1;; i++, mask <<= 1)
     229             :   {
     230    10415054 :     mask |= (a&1); a = (a+1)>>1;
     231     9282646 :     if (a==1) return mask | (1UL << i);
     232             :   }
     233             : }
     234             : 
     235             : /* Lift to precision p^e0.
     236             :  * a = modular factors of f mod (p,T) [possibly T=NULL]
     237             :  *  OR a TreeLift structure [e, link, v, w]: go on lifting
     238             :  * flag = 0: standard.
     239             :  * flag = 1: return TreeLift structure */
     240             : static GEN
     241       59334 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
     242             : {
     243       59334 :   long i, eold, e, k = lg(a) - 1;
     244             :   GEN E, v, w, link, penew, Tnew;
     245             :   ulong mask;
     246             :   pari_timer Ti;
     247             : 
     248       59334 :   if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
     249       59334 :   if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
     250       59334 :   if (e0 == 1) return a;
     251             : 
     252       59306 :   if (DEBUGLEVEL > 3) timer_start(&Ti);
     253       59306 :   if (typ(gel(a,1)) == t_INT)
     254             :   { /* a = TreeLift structure */
     255           0 :     e = itos(gel(a,1));
     256           0 :     link = gel(a,2);
     257           0 :     v    = gel(a,3);
     258           0 :     w    = gel(a,4);
     259             :   }
     260             :   else
     261             :   {
     262       59306 :     e = 1;
     263       59306 :     v = cgetg(2*k-2 + 1, t_VEC);
     264       59306 :     w = cgetg(2*k-2 + 1, t_VEC);
     265       59306 :     link=cgetg(2*k-2 + 1, t_VECSMALL);
     266       59306 :     BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
     267       59299 :     if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
     268             :   }
     269       59299 :   mask = quadratic_prec_mask(e0);
     270       59299 :   eold = 1;
     271       59299 :   penew = NULL;
     272       59299 :   Tnew = NULL;
     273       59299 :   if (DEBUGLEVEL > 3) err_printf("lifting to prec %ld\n", e0);
     274      358976 :   while (mask > 1)
     275             :   {
     276      240378 :     long enew = eold << 1;
     277      240378 :     if (mask & 1) enew--;
     278      240378 :     mask >>= 1;
     279      240378 :     if (enew >= e) { /* mask == 1: last iteration */
     280      240378 :       GEN peold = penew? penew: powiu(p, eold);
     281      240378 :       GEN Td = NULL, pd;
     282      240378 :       long d = enew - eold; /* = eold or eold-1 */
     283             :       /* lift from p^eold to p^enew */
     284      240378 :       pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
     285      240378 :       penew = mulii(peold,pd);
     286      240378 :       if (T) {
     287        3500 :         if (Tnew)
     288        2534 :           Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
     289             :         else
     290         966 :           Td = FpX_red(T, peold);
     291        3500 :         Tnew = FpX_red(T, penew);
     292        3500 :         ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
     293             :                          (flag == 0 && mask == 1));
     294             :       }
     295             :       else
     296      236878 :         ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
     297             :                         (flag == 0 && mask == 1));
     298      240378 :       if (DEBUGLEVEL > 3) timer_printf(&Ti, "reaching prec %ld", enew);
     299             :     }
     300      240378 :     eold = enew;
     301             :   }
     302             : 
     303       59299 :   if (flag)
     304         483 :     E = mkvec4(utoipos(e0), link, v, w);
     305             :   else
     306             :   {
     307       58816 :     E = cgetg(k+1, t_VEC);
     308      253974 :     for (i = 1; i <= 2*k-2; i++)
     309             :     {
     310      195158 :       long t = link[i];
     311      195158 :       if (t < 0) gel(E,-t) = gel(v,i);
     312             :     }
     313             :   }
     314       59299 :   return E;
     315             : }
     316             : 
     317             : /* Q list of (coprime, monic) factors of pol mod (T,p). Lift mod p^e = pe.
     318             :  * T may be NULL */
     319             : GEN
     320      101128 : ZpX_liftfact(GEN pol, GEN Q, GEN pe, GEN p, long e)
     321             : {
     322      101128 :   pari_sp av = avma;
     323      101128 :   pol = FpX_normalize(pol, pe);
     324      101128 :   if (lg(Q) == 2) return mkvec(pol);
     325       57878 :   return gerepilecopy(av, MultiLift(pol, Q, NULL, p, e, 0));
     326             : }
     327             : 
     328             : GEN
     329         973 : ZpXQX_liftfact(GEN pol, GEN Q, GEN T, GEN pe, GEN p, long e)
     330             : {
     331         973 :   pari_sp av = avma;
     332         973 :   pol = FpXQX_normalize(pol, T, pe);
     333         973 :   if (lg(Q) == 2) return mkvec(pol);
     334         973 :   return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
     335             : }
     336             : 
     337             : GEN
     338        2324 : ZqX_liftfact(GEN f, GEN a, GEN T, GEN pe, GEN p, long e)
     339        2324 : { return T ? ZpXQX_liftfact(f, a, T, pe, p, e): ZpX_liftfact(f, a, pe, p, e); }
     340             : GEN
     341          42 : ZqX_liftroot(GEN f, GEN a, GEN T, GEN p, long e)
     342          42 : { return T ? ZpXQX_liftroot(f, a,T , p, e): ZpX_liftroot(f, a, p, e); }
     343             : 
     344             : 
     345             : /* U = NULL treated as 1 */
     346             : static void
     347        2177 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
     348             : {
     349             :   GEN Q, R;
     350        2177 :   if (j < 0) return;
     351             : 
     352         847 :   Q = FpX_mul(gel(v,j), gel(w,j), pe);
     353         847 :   if (U)
     354             :   {
     355         364 :     Q = FpXQ_mul(Q, U, f, pe);
     356         364 :     R = FpX_sub(U, Q, pe);
     357             :   }
     358             :   else
     359         483 :     R = Fp_FpX_sub(gen_1, Q, pe);
     360         847 :   gel(w,j+1) = Q; /*  0 mod U v[j],  1 mod (1-U) v[j+1] */
     361         847 :   gel(w,j) = R; /*  1 mod U v[j],  0 mod (1-U) v[j+1] */
     362         847 :   BezoutPropagate(link, v, w, pe, R, f, link[j  ]);
     363         847 :   BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
     364             : }
     365             : 
     366             : /* as above, but return the Bezout coefficients for the lifted modular factors
     367             :  *   U[i] = 1 mod Qlift[i]
     368             :  *          0 mod Qlift[j], j != i */
     369             : GEN
     370         504 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
     371             : {
     372         504 :   pari_sp av = avma;
     373             :   GEN E, link, v, w, pe;
     374         504 :   long i, k = lg(Q)-1;
     375         504 :   if (k == 1) retmkvec(pol_1(varn(pol)));
     376         483 :   pe = powiu(p, e);
     377         483 :   pol = FpX_normalize(pol, pe);
     378         483 :   E = MultiLift(pol, Q, NULL, p, e, 1);
     379         483 :   link = gel(E,2);
     380         483 :   v    = gel(E,3);
     381         483 :   w    = gel(E,4);
     382         483 :   BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
     383         483 :   E = cgetg(k+1, t_VEC);
     384        2177 :   for (i = 1; i <= 2*k-2; i++)
     385             :   {
     386        1694 :     long t = link[i];
     387        1694 :     if (t < 0) E[-t] = w[i];
     388             :   }
     389         483 :   return gerepilecopy(av, E);
     390             : }
     391             : 
     392             : /* Front-end for ZpX_liftfact:
     393             :    lift the factorization of pol mod p given by L to p^N (if possible) */
     394             : GEN
     395          28 : polhensellift(GEN pol, GEN L, GEN Tp, long N)
     396             : {
     397             :   GEN T, p;
     398             :   long i, l;
     399          28 :   pari_sp av = avma;
     400             :   void (*chk)(GEN, const char*);
     401             : 
     402          28 :   if (typ(pol) != t_POL) pari_err_TYPE("polhensellift",pol);
     403          28 :   RgX_check_ZXX(pol, "polhensellift");
     404          28 :   if (!is_vec_t(typ(L)) || lg(L) < 3) pari_err_TYPE("polhensellift",L);
     405          28 :   if (N < 1) pari_err_DOMAIN("polhensellift", "precision", "<", gen_1,stoi(N));
     406          28 :   if (!ff_parse_Tp(Tp, &T, &p, 0)) pari_err_TYPE("polhensellift",Tp);
     407          28 :   chk = T? RgX_check_ZXX: RgX_check_ZX;
     408          28 :   l = lg(L); L = leafcopy(L);
     409          70 :   for (i = 1; i < l; i++)
     410             :   {
     411          49 :     GEN q = gel(L,i);
     412          49 :     if (typ(q) != t_POL) gel(L,i) = scalar_ZX_shallow(q, varn(pol));
     413          49 :     else chk(q, "polhensellift");
     414             :   }
     415          21 :   return gerepilecopy(av, ZqX_liftfact(pol, L, T, powiu(p,N), p, N));
     416             : }
     417             : 
     418             : static GEN
     419       47951 : FqV_roots_from_deg1(GEN x, GEN T, GEN p)
     420             : {
     421       47951 :   long i,l = lg(x);
     422       47951 :   GEN r = cgetg(l,t_COL);
     423       47951 :   for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = Fq_neg(gel(P,2), T, p); }
     424       47951 :   return r;
     425             : }
     426             : 
     427             : static GEN
     428       47909 : ZpX_liftroots_full(GEN f, GEN S, GEN q, GEN p, long e)
     429             : {
     430       47909 :   pari_sp av = avma;
     431       47909 :   GEN y = ZpX_liftfact(f, deg1_from_roots(S, varn(f)), q, p, e);
     432       47909 :   return gerepileupto(av, FqV_roots_from_deg1(y, NULL, q));
     433             : }
     434             : 
     435             : GEN
     436       47751 : ZpX_roots(GEN F, GEN p, long e)
     437             : {
     438       47751 :   pari_sp av = avma;
     439       47751 :   GEN q = powiu(p, e);
     440       47751 :   GEN f = FpX_normalize(F, p);
     441       47751 :   GEN g = FpX_normalize(FpX_split_part(f, p), p);
     442             :   GEN S;
     443       47751 :   if (degpol(g) < degpol(f))
     444             :   {
     445       41991 :     GEN h = FpX_div(f, g, p);
     446       41991 :     F = gel(ZpX_liftfact(F, mkvec2(g, h), q, p, e), 1);
     447             :   }
     448       47751 :   S = FpX_roots(g, p);
     449       47751 :   return gerepileupto(av, ZpX_liftroots_full(F, S, q, p, e));
     450             : }
     451             : 
     452             : static GEN
     453          42 : ZpXQX_liftroots_full(GEN f, GEN S, GEN T, GEN q, GEN p, long e)
     454             : {
     455          42 :   pari_sp av = avma;
     456          42 :   GEN y = ZpXQX_liftfact(f, deg1_from_roots(S, varn(f)), T, q, p, e);
     457          42 :   return gerepileupto(av, FqV_roots_from_deg1(y, T, q));
     458             : }
     459             : 
     460             : GEN
     461          42 : ZpXQX_roots(GEN F, GEN T, GEN p, long e)
     462             : {
     463          42 :   pari_sp av = avma;
     464          42 :   GEN q = powiu(p, e);
     465          42 :   GEN f = FpXQX_normalize(F, T, p);
     466          42 :   GEN g = FpXQX_normalize(FpXQX_split_part(f, T, p), T, p);
     467             :   GEN S;
     468          42 :   if (degpol(g) < degpol(f))
     469             :   {
     470          14 :     GEN h = FpXQX_div(f, g, T, p);
     471          14 :     F = gel(ZpXQX_liftfact(F, mkvec2(g, h), T, q, p, e), 1);
     472             :   }
     473          42 :   S = FpXQX_roots(g, T, p);
     474          42 :   return gerepileupto(av, ZpXQX_liftroots_full(F, S, T, q, p, e));
     475             : }
     476             : 
     477             : GEN
     478        2429 : ZqX_roots(GEN F, GEN T, GEN p, long e)
     479             : {
     480        2429 :   return T ? ZpXQX_roots(F, T, p, e): ZpX_roots(F, p, e);
     481             : }
     482             : 
     483             : /* SPEC:
     484             :  * p is a t_INT > 1, e >= 1
     485             :  * f is a ZX with leading term prime to p.
     486             :  * a is a simple root mod l for all l|p.
     487             :  * Return roots of f mod p^e, as integers (implicitly mod p^e)
     488             :  * STANDARD USE: p is a prime power */
     489             : GEN
     490         406 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
     491             : {
     492         406 :   pari_sp av = avma;
     493         406 :   GEN q = p, fr, W;
     494             :   ulong mask;
     495             : 
     496         406 :   a = modii(a,q);
     497         406 :   if (e == 1) return a;
     498         406 :   mask = quadratic_prec_mask(e);
     499         406 :   fr = FpX_red(f,q);
     500         406 :   W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
     501             :   for(;;)
     502             :   {
     503        2632 :     q = sqri(q);
     504        1519 :     if (mask & 1) q = diviiexact(q, p);
     505        1519 :     mask >>= 1;
     506        1519 :     fr = FpX_red(f,q);
     507        1519 :     a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
     508        1519 :     if (mask == 1) return gerepileuptoint(av, a);
     509        1113 :     W = Fp_sub(shifti(W,1), Fp_mul(Fp_sqr(W,q), FpX_eval(ZX_deriv(fr),a,q), q), q);
     510             :   }
     511             : }
     512             : 
     513             : GEN
     514         158 : ZpX_liftroots(GEN f, GEN S, GEN p, long e)
     515             : {
     516         158 :   long i, n = lg(S)-1, d = degpol(f);
     517             :   GEN r;
     518         158 :   if (n == d) return ZpX_liftroots_full(f, S, powiu(p, e), p, e);
     519           0 :   r = cgetg(n+1, typ(S));
     520           0 :   for (i=1; i <= n; i++)
     521           0 :     gel(r,i) = ZpX_liftroot(f, gel(S,i), p, e);
     522           0 :   return r;
     523             : }
     524             : 
     525             : GEN
     526         119 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
     527             : {
     528         119 :   pari_sp av = avma, av2;
     529         119 :   GEN pv = p, q, W, df, Tq, fr, dfr;
     530             :   ulong mask;
     531         119 :   a = Fq_red(a, T, p);
     532         119 :   if (e <= v+1) return a;
     533         119 :   df = RgX_deriv(f);
     534         119 :   if (v) { pv = powiu(p,v); df = ZXX_Z_divexact(df, pv); }
     535         119 :   mask = quadratic_prec_mask(e-v);
     536         119 :   Tq = FpXT_red(T, p); dfr = FpXQX_red(df, Tq, p);
     537         119 :   W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
     538         119 :   q = p;
     539         119 :   av2 = avma;
     540             :   for (;;)
     541         175 :   {
     542             :     GEN u, fa, qv, q2v, q2, Tq2;
     543         294 :     q2 = q; q = sqri(q);
     544         294 :     if (mask & 1) q = diviiexact(q,p);
     545         294 :     mask >>= 1;
     546         294 :     if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
     547         294 :     else { qv = q; q2v = q2; }
     548         294 :     Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
     549         294 :     fr = FpXQX_red(f, Tq, qv);
     550         294 :     fa = FqX_eval(fr, a, Tq, qv);
     551         294 :     fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
     552         294 :     a = Fq_sub(a, Fq_mul(Fq_mul(W,fa,Tq2,q2v),q2, Tq,qv), Tq, qv);
     553         294 :     if (mask == 1) return gerepileupto(av, a);
     554         175 :     dfr = FpXQX_red(df, Tq, q);
     555         175 :     u = Fq_sub(Fq_mul(W,FqX_eval(dfr,a,Tq,q),Tq,q),gen_1,Tq,q);
     556         175 :     u = typ(u)==t_INT? diviiexact(u,q2): ZX_Z_divexact(u,q2);
     557         175 :     W = Fq_sub(W, Fq_mul(Fq_mul(u,W,Tq2,q2),q2, Tq,q),Tq,q);
     558         175 :     if (gc_needed(av2,2))
     559             :     {
     560           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_liftroot, e = %ld", e);
     561           0 :       gerepileall(av2, 3, &a, &W, &q);
     562             :     }
     563             :   }
     564             : }
     565             : 
     566             : GEN
     567         119 : ZpXQX_liftroot(GEN f, GEN a, GEN T, GEN p, long e) { return ZpXQX_liftroot_vald(f,a,0,T,p,e); }
     568             : 
     569             : GEN
     570           0 : ZpXQX_liftroots(GEN f, GEN S, GEN T, GEN p, long e)
     571             : {
     572           0 :   long i, n = lg(S)-1, d = degpol(f);
     573             :   GEN r;
     574           0 :   if (n == d) return ZpXQX_liftroots_full(f, S, T, powiu(p, e), p, e);
     575           0 :   r = cgetg(n+1, typ(S));
     576           0 :   for (i=1; i <= n; i++)
     577           0 :     gel(r,i) = ZpXQX_liftroot(f, gel(S,i), T, p, e);
     578           0 :   return r;
     579             : }
     580             : 
     581             : /* Same as ZpX_liftroot for the polynomial X^n-b*/
     582             : GEN
     583      177303 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
     584             : {
     585      177303 :   pari_sp ltop=avma;
     586             :   GEN q, w, n_1;
     587             :   ulong mask;
     588      177303 :   long pis2 = equalii(n, gen_2)? 1: 0;
     589      177303 :   if (e == 1) return icopy(a);
     590      172431 :   n_1 = subiu(n,1);
     591      172431 :   mask = quadratic_prec_mask(e);
     592      172431 :   w = Fp_inv(pis2 ? shifti(a,1): Fp_mul(n,Fp_pow(a,n_1,p), p), p);
     593      172431 :   q = p;
     594             :   for(;;)
     595             :   {
     596      366023 :     q = sqri(q);
     597      269227 :     if (mask & 1) q = diviiexact(q, p);
     598      269227 :     mask >>= 1;
     599      269227 :     if (lgefint(q) == 3 && lgefint(n) == 3)
     600       96638 :     {
     601      268628 :       ulong Q = uel(q,2), N = uel(n,2);
     602      268628 :       ulong A = umodiu(a, Q);
     603      268628 :       ulong B = umodiu(b, Q);
     604      268628 :       ulong W = umodiu(w, Q);
     605             : 
     606      268628 :       A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
     607      268628 :       a = utoi(A);
     608      268628 :       if (mask == 1) break;
     609       96638 :       W = Fl_sub(Fl_add(W,W,Q),
     610             :                  Fl_mul(Fl_sqr(W,Q), Fl_mul(N,Fl_powu(A, N-1, Q), Q), Q), Q);
     611       96638 :       w = utoi(W);
     612             :     }
     613             :     else
     614             :     {
     615             :       /* a -= w (a^n - b) */
     616         599 :       a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q),b))), q);
     617         599 :       if (mask == 1) break;
     618             :       /* w += w - w^2 n a^(n-1)*/
     619         301 :       w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
     620         143 :                            pis2? shifti(a,1): mulii(n,Fp_pow(a,n_1,q)), q));
     621             :     }
     622             :   }
     623      172431 :   return gerepileuptoint(ltop,a);
     624             : }
     625             : 
     626             : 
     627             : /* Same as ZpX_liftroot for the polynomial X^2-b */
     628             : GEN
     629      139769 : Zp_sqrtlift(GEN b, GEN a, GEN p, long e)
     630             : {
     631      139769 :   return Zp_sqrtnlift(b, gen_2, a, p, e);
     632             : }
     633             : 
     634             : GEN
     635     1352379 : Zp_sqrt(GEN x, GEN p, long e)
     636             : {
     637             :   pari_sp av;
     638             :   GEN z;
     639     1352379 :   if (absequaliu(p,2)) return Z2_sqrt(x,e);
     640     1350874 :   av = avma;
     641     1350874 :   z = Fp_sqrt(Fp_red(x, p), p);
     642     1350874 :   if (!z) return NULL;
     643     1350825 :   if (e > 1) z = Zp_sqrtlift(x, z, p, e);
     644     1350825 :   return gerepileuptoint(av, z);
     645             : }
     646             : 
     647             : /* Compute (x-1)/(x+1)/p^k */
     648             : static GEN
     649       20880 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
     650             : {
     651       20880 :   pari_sp av = avma;
     652       20880 :   long vT = get_FpX_var(T);
     653             :   GEN bn, bdi;
     654       20880 :   GEN bd = ZX_Z_add(x, gen_1);
     655       20879 :   if (absequaliu(p,2)) /*For p=2, we need to simplify by 2*/
     656             :   {
     657        6949 :     bn = ZX_shifti(x,-(k+1));
     658        6950 :     bdi= ZpXQ_invlift(ZX_shifti(bd ,-1), pol_1(vT), T, p, e);
     659             :   }
     660             :   else
     661             :   {
     662       13930 :     bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
     663       13930 :     bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
     664             :   }
     665       20880 :   return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
     666             : }
     667             : 
     668             : /* Assume p odd, a = 1 [p], return log(a) */
     669             : GEN
     670       20880 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
     671             : {
     672       20880 :   pari_sp av = avma;
     673             :   pari_timer ti;
     674       20880 :   long is2 = absequaliu(p,2);
     675       20880 :   ulong pp = is2 ? 0: itou_or_0(p);
     676       20880 :   double lp = is2 ? 1: pp ? log2(pp): expi(p);
     677       20880 :   long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
     678             :   GEN ak, s, b, pol;
     679       20880 :   long e = is2 ? N-1: N;
     680       20880 :   long i, l = (e-2)/(2*(k+is2));
     681       20880 :   GEN pe = powiu(p,e);
     682       20880 :   GEN TNk, pNk = powiu(p,N+k);
     683       20880 :   if( DEBUGLEVEL>=3) timer_start(&ti);
     684       20880 :   TNk = FpX_get_red(get_FpX_mod(T), pNk);
     685       20880 :   ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
     686       20880 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
     687       20880 :   b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
     688       20880 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
     689       20880 :   pol= cgetg(l+3,t_POL);
     690       20880 :   pol[1] = evalsigne(1)|evalvarn(0);
     691       56376 :   for(i=0; i<=l; i++)
     692             :   {
     693             :     GEN g;
     694       35496 :     ulong z = 2*i+1;
     695       35496 :     if (pp)
     696             :     {
     697       26222 :       long w = u_lvalrem(z, pp, &z);
     698       26222 :       g = powuu(pp,2*i*k-w);
     699             :     }
     700        9274 :     else g = powiu(p,2*i*k);
     701       35496 :     gel(pol,i+2) = Fp_divu(g, z,pe);
     702             :   }
     703       20880 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
     704       20880 :   s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T,  pe);
     705       20880 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
     706       20880 :   s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
     707       20880 :   return gerepileupto(av, is2? s: FpX_red(s, pe));
     708             : }
     709             : 
     710             : /***********************************************************************/
     711             : /**                                                                   **/
     712             : /**                 Generic quadratic hensel lift over Zp[X]          **/
     713             : /**                                                                   **/
     714             : /***********************************************************************/
     715             : /* q = p^N */
     716             : 
     717             : GEN
     718           0 : gen_ZpM_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
     719             :                             GEN lin(void *E, GEN F, GEN d, GEN q),
     720             :                             GEN invl(void *E, GEN d))
     721             : {
     722           0 :   pari_sp av = avma;
     723             :   long N2, M;
     724             :   GEN VN2, V2, VM, bil;
     725             :   GEN q2, qM;
     726           0 :   V = FpM_red(V, q);
     727           0 :   if (N == 1) return invl(E, V);
     728           0 :   N2 = (N + 1)>>1; M = N - N2;
     729           0 :   F = FpM_red(F, q);
     730           0 :   qM = powiu(p, M);
     731           0 :   q2 = M == N2? qM: mulii(qM, p);
     732             :   /* q2 = p^N2, qM = p^M, q = q2 * qM */
     733           0 :   VN2 = gen_ZpM_Dixon(F, V, q2, p, N2, E, lin, invl);
     734           0 :   bil = lin(E, F, VN2, q);
     735           0 :   V2 = ZM_Z_divexact(ZM_sub(V, bil), q2);
     736           0 :   VM = gen_ZpM_Dixon(F, V2, qM, p, M, E, lin, invl);
     737           0 :   return gerepileupto(av, FpM_red(ZM_add(VN2, ZM_Z_mul(VM, q2)), q));
     738             : }
     739             : 
     740             : GEN
     741      199829 : gen_ZpX_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
     742             :                             GEN lin(void *E, GEN F, GEN d, GEN q),
     743             :                             GEN invl(void *E, GEN d))
     744             : {
     745      199829 :   pari_sp av = avma;
     746             :   long N2, M;
     747             :   GEN VN2, V2, VM, bil;
     748             :   GEN q2, qM;
     749      199829 :   V = FpX_red(V, q);
     750      199829 :   if (N == 1) return invl(E, V);
     751       45458 :   N2 = (N + 1)>>1; M = N - N2;
     752       45458 :   F = FpXT_red(F, q);
     753       45458 :   qM = powiu(p, M);
     754       45458 :   q2 = M == N2? qM: mulii(qM, p);
     755             :   /* q2 = p^N2, qM = p^M, q = q2 * qM */
     756       45458 :   VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
     757       45458 :   bil = lin(E, F, VN2, q);
     758       45458 :   V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
     759       45458 :   VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
     760       45458 :   return gerepileupto(av, FpX_red(ZX_add(VN2, ZX_Z_mul(VM, q2)), q));
     761             : }
     762             : 
     763             : GEN
     764      525353 : gen_ZpM_Newton(GEN x, GEN p, long n, void *E,
     765             :                       GEN eval(void *E, GEN f, GEN q),
     766             :                       GEN invd(void *E, GEN V, GEN v, GEN q, long M))
     767             : {
     768      525353 :   pari_sp ltop = avma, av;
     769      525353 :   long N = 1, N2, M;
     770             :   long mask;
     771      525353 :   GEN q = p;
     772      525353 :   if (n == 1) return gcopy(x);
     773      525353 :   mask = quadratic_prec_mask(n);
     774      525353 :   av = avma;
     775     1606074 :   while (mask > 1)
     776             :   {
     777             :     GEN qM, q2, v, V;
     778      555368 :     N2 = N; N <<= 1;
     779      555368 :     q2 = q;
     780      555368 :     if (mask&1UL) { /* can never happen when q2 = p */
     781       21910 :       N--; M = N2-1;
     782       21910 :       qM = diviiexact(q2,p); /* > 1 */
     783       21910 :       q = mulii(qM,q2);
     784             :     } else {
     785      533458 :       M = N2;
     786      533458 :       qM = q2;
     787      533458 :       q = sqri(q2);
     788             :     }
     789             :     /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
     790      555368 :     mask >>= 1;
     791      555368 :     v = eval(E, x, q);
     792      555368 :     V = ZM_Z_divexact(gel(v,1), q2);
     793      555368 :     x = FpM_sub(x, ZM_Z_mul(invd(E, V, v, qM, M), q2), q);
     794      555368 :     if (gc_needed(av, 1))
     795             :     {
     796           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
     797           0 :       gerepileall(av, 2, &x, &q);
     798             :     }
     799             :   }
     800      525353 :   return gerepileupto(ltop, x);
     801             : }
     802             : 
     803             : static GEN
     804      555368 : _ZpM_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
     805             : {
     806             :   (void)E; (void)M;
     807      555368 :   return FpM_mul(V, gel(v,2), q);
     808             : }
     809             : 
     810             : static GEN
     811      555368 : _ZpM_eval(void *E, GEN x, GEN q)
     812             : {
     813      555368 :   GEN f = RgM_Rg_sub_shallow(FpM_mul(x, FpM_red((GEN) E, q), q), gen_1);
     814      555368 :   return mkvec2(f, x);
     815             : }
     816             : 
     817             : GEN
     818      525353 : ZpM_invlift(GEN M, GEN C, GEN p, long n)
     819             : {
     820      525353 :   return gen_ZpM_Newton(C, p, n, (void *)M, _ZpM_eval, _ZpM_invd);
     821             : }
     822             : 
     823             : GEN
     824      265856 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
     825             :                       GEN eval(void *E, GEN f, GEN q),
     826             :                       GEN invd(void *E, GEN V, GEN v, GEN q, long M))
     827             : {
     828      265856 :   pari_sp ltop = avma, av;
     829      265856 :   long N = 1, N2, M;
     830             :   long mask;
     831      265856 :   GEN q = p;
     832      265856 :   if (n == 1) return gcopy(x);
     833      262027 :   mask = quadratic_prec_mask(n);
     834      262027 :   av = avma;
     835     1187441 :   while (mask > 1)
     836             :   {
     837             :     GEN qM, q2, v, V;
     838      663387 :     N2 = N; N <<= 1;
     839      663387 :     q2 = q;
     840      663387 :     if (mask&1UL) { /* can never happen when q2 = p */
     841      275742 :       N--; M = N2-1;
     842      275742 :       qM = diviiexact(q2,p); /* > 1 */
     843      275743 :       q = mulii(qM,q2);
     844             :     } else {
     845      387645 :       M = N2;
     846      387645 :       qM = q2;
     847      387645 :       q = sqri(q2);
     848             :     }
     849             :     /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
     850      663379 :     mask >>= 1;
     851      663379 :     v = eval(E, x, q);
     852      663391 :     V = ZX_Z_divexact(gel(v,1), q2);
     853      663383 :     x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
     854      663387 :     if (gc_needed(av, 1))
     855             :     {
     856           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
     857           0 :       gerepileall(av, 2, &x, &q);
     858             :     }
     859             :   }
     860      262027 :   return gerepileupto(ltop, x);
     861             : }
     862             : 
     863             : struct _ZpXQ_inv
     864             : {
     865             :   GEN T, a, p ,n;
     866             : };
     867             : 
     868             : static GEN
     869      521479 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
     870             : {
     871      521479 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
     872      521479 :   GEN Tq = FpXT_red(d->T, q);
     873             :   (void)M;
     874      521479 :   return FpXQ_mul(V, gel(v,2), Tq, q);
     875             : }
     876             : 
     877             : static GEN
     878      521477 : _inv_eval(void *E, GEN x, GEN q)
     879             : {
     880      521477 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
     881      521477 :   GEN Tq = FpXT_red(d->T, q);
     882      521477 :   GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
     883      521479 :   return mkvec2(f, x);
     884             : }
     885             : 
     886             : GEN
     887      201193 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
     888             : {
     889             :   struct _ZpXQ_inv d;
     890      201193 :   d.a = a; d.T = T; d.p = p;
     891      201193 :   return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
     892             : }
     893             : 
     894             : GEN
     895      180313 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
     896             : {
     897      180313 :   pari_sp av=avma;
     898             :   GEN ai;
     899      180313 :   if (lgefint(p)==3)
     900             :   {
     901      180285 :     ulong pp = p[2];
     902      180285 :     ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
     903             :   } else
     904          28 :     ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
     905      180313 :   return gerepileupto(av, ZpXQ_invlift(a, ai, T, p, e));
     906             : }
     907             : 
     908             : GEN
     909       31584 : ZpXQ_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
     910             : {
     911       31584 :   return FpXQ_mul(a, ZpXQ_inv(b, T, p, e), T, q);
     912             : }
     913             : 
     914             : GEN
     915      146440 : ZpXQX_divrem(GEN x, GEN Sp, GEN T, GEN q, GEN p, long e, GEN *pr)
     916             : {
     917      146440 :   pari_sp av = avma;
     918      146440 :   GEN S = get_FpXQX_mod(Sp);
     919      146440 :   GEN b = leading_coeff(S), bi;
     920             :   GEN S2, Q;
     921      146440 :   if (typ(b)==t_INT) return FpXQX_divrem(x, Sp, T, q, pr);
     922      146440 :   bi = ZpXQ_inv(b, T, p, e);
     923      146440 :   S2 = FqX_Fq_mul_to_monic(S, bi, T, q);
     924      146440 :   Q = FpXQX_divrem(x, S2, T, q, pr);
     925      146440 :   if (pr==ONLY_DIVIDES && !Q) { set_avma(av); return NULL; }
     926      146440 :   if (pr==ONLY_REM || pr==ONLY_DIVIDES) return gerepileupto(av, Q);
     927      146440 :   Q = FpXQX_FpXQ_mul(Q, bi, T, q);
     928      146440 :   gerepileall(av, 2, &Q, pr);
     929      146440 :   return Q;
     930             : }
     931             : 
     932             : GEN
     933         623 : ZpXQX_digits(GEN x, GEN B, GEN T, GEN q, GEN p, long e)
     934             : {
     935         623 :   pari_sp av = avma;
     936         623 :   GEN b = leading_coeff(B), bi;
     937             :   GEN B2, P, V, W;
     938             :   long i, lV;
     939         623 :   if (typ(b)==t_INT) return FpXQX_digits(x, B, T, q);
     940         623 :   bi = ZpXQ_inv(b, T, p, e);
     941         623 :   B2 = FqX_Fq_mul_to_monic(B, bi, T, q);
     942         623 :   V = FpXQX_digits(x, B2, T, q);
     943         623 :   lV = lg(V)-1;
     944         623 :   P = FpXQ_powers(bi, lV-1, T, q);
     945         623 :   W = cgetg(lV+1, t_VEC);
     946       11200 :   for(i=1; i<=lV; i++)
     947       10577 :     gel(W, i) = FpXQX_FpXQ_mul(gel(V,i), gel(P, i), T, q);
     948         623 :   return gerepileupto(av, W);
     949             : }
     950             : 
     951             : struct _ZpXQ_sqrtn
     952             : {
     953             :   GEN T, a, n, ai;
     954             : };
     955             : 
     956             : static GEN
     957        2037 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
     958             : {
     959        2037 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
     960        2037 :   GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
     961             :   (void)M;
     962        2037 :   return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
     963             : }
     964             : 
     965             : static GEN
     966        2037 : _sqrtn_eval(void *E, GEN x, GEN q)
     967             : {
     968        2037 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
     969        2037 :   GEN Tq = FpX_red(d->T, q);
     970        2037 :   GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
     971        2037 :   return mkvec2(f, x);
     972             : }
     973             : 
     974             : GEN
     975        1204 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
     976             : {
     977             :   struct _ZpXQ_sqrtn d;
     978        1204 :   d.a = a; d.T = T; d.n = n;
     979        1204 :   d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
     980        1204 :   return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
     981             : }
     982             : 
     983             : static GEN
     984          98 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
     985             : 
     986             : GEN
     987         665 : Zq_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
     988             : {
     989          49 :   return T? ZpXQ_sqrtnlift(to_ZX(a,varn(T)), n, to_ZX(x,varn(T)), T, p, e)
     990         714 :           : Zp_sqrtnlift(a, n, x, p, e);
     991             : }
     992             : 
     993             : GEN
     994           0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
     995             : {
     996           0 :   pari_sp av = avma;
     997           0 :   GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
     998           0 :   if (!z) return NULL;
     999           0 :   if (e <= 1) return gerepileupto(av, z);
    1000           0 :   return gerepileupto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
    1001             : }
    1002             : 
    1003             : GEN
    1004        5614 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
    1005             :                      GEN early(void *E, GEN x, GEN q))
    1006             : {
    1007        5614 :   pari_sp ltop = avma, av;
    1008             :   long N, r;
    1009             :   long mask;
    1010             :   GEN q2, q, W, Q, Tq2, Tq, Pq;
    1011             :   pari_timer ti;
    1012        5614 :   T = FpX_get_red(T, powiu(p, n));
    1013        5614 :   if (n == 1) return gcopy(S);
    1014        5614 :   mask = quadratic_prec_mask(n);
    1015        5614 :   av = avma;
    1016        5614 :   q2 = p; q = sqri(p); mask >>= 1; N = 2;
    1017        5614 :   if (DEBUGLEVEL > 3) timer_start(&ti);
    1018        5614 :   Tq = FpXT_red(T,q);
    1019        5614 :   Tq2 = FpXT_red(Tq,q2);
    1020        5614 :   Pq = FpX_red(P,q);
    1021        5614 :   W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
    1022        5614 :   Q  = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
    1023        5614 :   r = brent_kung_optpow(degpol(P), 4, 3);
    1024        5614 :   if (DEBUGLEVEL > 3)
    1025           0 :     err_printf("ZpX_ZpXQ_liftroot: lifting to prec %ld\n",n);
    1026             :   for (;;)
    1027       12559 :   {
    1028             :     GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
    1029       18173 :     H  = FpXQ_mul(W, Q, Tq2, q2);
    1030       18173 :     Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
    1031       18173 :     if (DEBUGLEVEL > 3)
    1032           0 :       timer_printf(&ti,"ZpX_ZpXQ_liftroot: reaching prec %ld",N);
    1033       18173 :     if (mask==1)
    1034        3529 :       return gerepileupto(ltop, Sq);
    1035       14644 :     if (early)
    1036             :     {
    1037       13118 :       GEN Se = early(E, Sq, q);
    1038       13118 :       if (Se) return gerepileupto(ltop, Se);
    1039             :     }
    1040       12559 :     qq = sqri(q); N <<= 1;
    1041       12559 :     if (mask&1UL) { qq = diviiexact(qq, p); N--; }
    1042       12559 :     mask >>= 1;
    1043       12559 :     Pqq  = FpX_red(P, qq);
    1044       12559 :     Tqq  = FpXT_red(T, qq);
    1045       12559 :     Spow = FpXQ_powers(Sq, r, Tqq, qq);
    1046       12559 :     Q  = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
    1047       12559 :     dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
    1048       12559 :     Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
    1049       12559 :     Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
    1050       12559 :     Wq = FpX_sub(W, Wq, q);
    1051       12559 :     S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
    1052       12559 :     if (gc_needed(av, 1))
    1053             :     {
    1054           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_Newton");
    1055           0 :       gerepileall(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
    1056             :     }
    1057             :   }
    1058             : }
    1059             : 
    1060             : GEN
    1061        1694 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
    1062             : {
    1063        1694 :   return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
    1064             : }
    1065             : 
    1066             : GEN
    1067         301 : ZpX_Frobenius(GEN T, GEN p, long e)
    1068             : {
    1069         301 :   return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
    1070             : }
    1071             : 
    1072             : GEN
    1073         147 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
    1074             : {
    1075         147 :   pari_sp av = avma;
    1076         147 :   GEN xp = ZpX_Frobenius(T, p, e);
    1077         147 :   GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
    1078         147 :   return gerepilecopy(av, gel(z,2));
    1079             : }
    1080             : 
    1081             : /* Canonical lift of polynomial */
    1082             : 
    1083       11130 : static GEN _can_invl(void *E, GEN V) {(void) E; return V; }
    1084             : 
    1085        3654 : static GEN _can_lin(void *E, GEN F, GEN V, GEN q)
    1086             : {
    1087        3654 :   GEN v = RgX_splitting(V, 3);
    1088             :   (void) E;
    1089        3654 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
    1090             : }
    1091             : 
    1092             : static GEN
    1093        7476 : _can_iter(void *E, GEN f, GEN q)
    1094             : {
    1095        7476 :   GEN h = RgX_splitting(f,3);
    1096        7476 :   GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
    1097        7476 :   GEN h12 = ZX_mul(gel(h,1), gel(h,2));
    1098        7476 :   GEN h13 = ZX_mul(gel(h,1), gel(h,3));
    1099        7476 :   GEN h23 = ZX_mul(gel(h,2), gel(h,3));
    1100        7476 :   GEN h1c = ZX_mul(gel(h,1), h1s);
    1101        7476 :   GEN h3c = ZX_mul(gel(h,3), h3s);
    1102        7476 :   GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
    1103        7476 :   GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
    1104             :   (void) E;
    1105        7476 :   return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
    1106             : }
    1107             : 
    1108             : static GEN
    1109        7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
    1110             : {
    1111        7476 :   GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
    1112        7476 :   GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
    1113        7476 :   GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
    1114             :                  ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
    1115             :   (void)E;
    1116        7476 :   return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
    1117             :                                                  _can_lin, _can_invl);
    1118             : }
    1119             : 
    1120             : static GEN
    1121        3717 : F3x_frobeniuslift(GEN P, long n)
    1122        3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
    1123             : 
    1124       29736 : static GEN _can5_invl(void *E, GEN V) {(void) E; return V; }
    1125             : 
    1126        9107 : static GEN _can5_lin(void *E, GEN F, GEN V, GEN q)
    1127             : {
    1128        9107 :   ulong p = *(ulong*)E;
    1129        9107 :   GEN v = RgX_splitting(V, p);
    1130        9107 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
    1131             : }
    1132             : 
    1133             : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
    1134             : static GEN
    1135       62146 : _shift(GEN P, long n, ulong p, long v)
    1136             : {
    1137       62146 :   long i, l=lg(P);
    1138       62146 :   GEN r = cgetg(l,t_POL); r[1] = P[1];
    1139      483238 :   for(i=2;i<l;i++)
    1140             :   {
    1141      421092 :     long s = n*(i-2)%p;
    1142      421092 :     GEN ci = gel(P,i);
    1143      421092 :     if (typ(ci)==t_INT)
    1144      104979 :       gel(r,i) = monomial(ci, s, v);
    1145             :     else
    1146      316113 :       gel(r,i) = RgX_rotate_shallow(ci, s, p);
    1147             :   }
    1148       62146 :   return FpXX_renormalize(r, l);
    1149             : }
    1150             : 
    1151             : struct _can_mul
    1152             : {
    1153             :   GEN T, q;
    1154             :   ulong p;
    1155             : };
    1156             : 
    1157             : static GEN
    1158       41517 : _can5_mul(void *E, GEN A, GEN B)
    1159             : {
    1160       41517 :   struct _can_mul *d = (struct _can_mul *)E;
    1161       41517 :   GEN a = gel(A,1), b = gel(B,1);
    1162       41517 :   long n = itos(gel(A,2));
    1163       41517 :   GEN bn = _shift(b, n, d->p, get_FpX_var(d->T));
    1164       41517 :   GEN c = FpXQX_mul(a, bn, d->T, d->q);
    1165       41517 :   return mkvec2(c, addii(gel(A,2), gel(B,2)));
    1166             : }
    1167             : 
    1168             : static GEN
    1169       41349 : _can5_sqr(void *E, GEN A)
    1170             : {
    1171       41349 :   return _can5_mul(E,A,A);
    1172             : }
    1173             : 
    1174             : static GEN
    1175       20629 : _can5_iter(void *E, GEN f, GEN q)
    1176             : {
    1177       20629 :   pari_sp av = avma;
    1178             :   struct _can_mul D;
    1179       20629 :   ulong p = *(ulong*)E;
    1180       20629 :   long i, vT = fetch_var();
    1181             :   GEN N, P, d, V, fs;
    1182       20629 :   D.q = q; D.T = ZX_Z_sub(pol_xn(p,vT),gen_1);
    1183       20629 :   D.p = p;
    1184       20629 :   fs = mkvec2(_shift(f, 1, p, vT), gen_1);
    1185       20629 :   N = gel(gen_powu_i(fs,p-1,(void*)&D,_can5_sqr,_can5_mul),1);
    1186       20629 :   N = ZXX_evalx0(FpXQX_red(N,polcyclo(p,vT),q));
    1187       20629 :   P = FpX_mul(N,f,q);
    1188       20629 :   P = RgX_deflate(P, p);
    1189       20629 :   d = RgX_splitting(N, p);
    1190       20629 :   V = cgetg(p+1,t_VEC);
    1191       20629 :   gel(V,1) = ZX_mulu(gel(d,1), p);
    1192      103915 :   for(i=2; i<= (long)p; i++)
    1193       83286 :     gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
    1194       20629 :   (void)delete_var(); return gerepilecopy(av, mkvec2(ZX_sub(f,P),V));
    1195             : }
    1196             : 
    1197             : static GEN
    1198       20629 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
    1199             : {
    1200       20629 :   ulong p = *(long*)E;
    1201       20629 :   return gen_ZpX_Dixon(gel(v,2), H, qM, utoi(p), M, E, _can5_lin, _can5_invl);
    1202             : }
    1203             : 
    1204             : GEN
    1205       13958 : Flx_Teichmuller(GEN P, ulong p, long n)
    1206             : {
    1207       24199 :   return p==3 ? F3x_frobeniuslift(P,n):
    1208       10241 :          gen_ZpX_Newton(Flx_to_ZX(P),utoi(p), n, &p, _can5_iter, _can5_invd);
    1209             : }
    1210             : 
    1211             : GEN
    1212          35 : polteichmuller(GEN P, ulong p, long n)
    1213             : {
    1214          35 :   pari_sp av = avma;
    1215          35 :   GEN q = NULL;
    1216          35 :   if (typ(P)!=t_POL || !RgX_is_FpX(P,&q))
    1217           0 :     pari_err_TYPE("polteichmuller",P);
    1218          35 :   if (q && !equaliu(q,p))
    1219           0 :     pari_err_MODULUS("polteichmuller",q,utoi(p));
    1220          35 :   if (n <= 0)
    1221           0 :     pari_err_DOMAIN("polteichmuller", "precision", "<=",gen_0,stoi(n));
    1222          63 :   return gerepileupto(av, p==2 ? F2x_Teichmuller(RgX_to_F2x(P), n)
    1223          28 :                                : Flx_Teichmuller(RgX_to_Flx(P, p), p, n));
    1224             : }

Generated by: LCOV version 1.13