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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16962-5a32637) Lines: 487 508 95.9 %
Date: 2014-10-29 Functions: 35 36 97.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 186 232 80.2 %

           Branch data     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                 :       7014 : BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
      36                 :            : {
      37                 :       7014 :   long k = lg(a)-1;
      38                 :            :   long i, j, s, minp, mind;
      39                 :            : 
      40         [ +  + ]:      43477 :   for (i=1; i<=k; i++) { gel(V,i) = gel(a,i); link[i] = -i; }
      41                 :            : 
      42         [ +  + ]:      29449 :   for (j=1; j <= 2*k-5; j+=2,i++)
      43                 :            :   {
      44                 :      22435 :     minp = j;
      45                 :      22435 :     mind = degpol(gel(V,j));
      46         [ +  + ]:     358050 :     for (s=j+1; s<i; s++)
      47         [ +  + ]:     335615 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
      48                 :            : 
      49                 :      22435 :     swap(gel(V,j), gel(V,minp));
      50                 :      22435 :     lswap(link[j], link[minp]);
      51                 :            : 
      52                 :      22435 :     minp = j+1;
      53                 :      22435 :     mind = degpol(gel(V,j+1));
      54         [ +  + ]:     335615 :     for (s=j+2; s<i; s++)
      55         [ +  + ]:     313180 :       if (degpol(gel(V,s)) < mind) { minp = s; mind = degpol(gel(V,s)); }
      56                 :            : 
      57                 :      22435 :     swap(gel(V,j+1), gel(V,minp));
      58                 :      22435 :     lswap(link[j+1], link[minp]);
      59                 :            : 
      60                 :      22435 :     gel(V,i) = FqX_mul(gel(V,j), gel(V,j+1), T, p);
      61                 :      22435 :     link[i] = j;
      62                 :            :   }
      63                 :            : 
      64         [ +  + ]:      36456 :   for (j=1; j <= 2*k-3; j+=2)
      65                 :            :   {
      66                 :            :     GEN d, u, v;
      67                 :      29449 :     d = FqX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
      68         [ +  + ]:      29449 :     if (degpol(d) > 0) pari_err_COPRIME("BuildTree", gel(V,j), gel(V,j+1));
      69                 :      29442 :     d = gel(d,2);
      70         [ +  + ]:      29442 :     if (!gequal1(d))
      71                 :            :     {
      72         [ +  + ]:      26218 :       if (typ(d)==t_POL)
      73                 :            :       {
      74                 :       4242 :         d = FpXQ_inv(d, T, p);
      75                 :       4242 :         u = FqX_Fq_mul(u, d, T, p);
      76                 :       4242 :         v = FqX_Fq_mul(v, d, T, p);
      77                 :            :       }
      78                 :            :       else
      79                 :            :       {
      80                 :      21976 :         d = Fp_inv(d, p);
      81                 :      21976 :         u = FqX_Fp_mul(u, d, T,p);
      82                 :      21976 :         v = FqX_Fp_mul(v, d, T,p);
      83                 :            :       }
      84                 :            :     }
      85                 :      29442 :     gel(W,j) = u;
      86                 :      29442 :     gel(W,j+1) = v;
      87                 :            :   }
      88                 :       7007 : }
      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                 :     145736 : ZpX_HenselLift(GEN V, GEN W, long j, GEN f, GEN pd, GEN p0, GEN p1, int noinv)
      94                 :            : {
      95                 :     145736 :   pari_sp av = avma;
      96                 :     145736 :   long space = lg(f) * lgefint(p1);
      97                 :            :   GEN a2, b2, g, z, s, t;
      98                 :     145736 :   GEN a = gel(V,j), b = gel(V,j+1);
      99                 :     145736 :   GEN u = gel(W,j), v = gel(W,j+1);
     100                 :            : 
     101                 :     145736 :   (void)new_chunk(space); /* HACK */
     102                 :     145736 :   g = ZX_sub(f, ZX_mul(a,b));
     103                 :     145736 :   g = ZX_Z_divexact(g, p0);
     104                 :     145736 :   g = FpX_red(g, pd);
     105                 :     145736 :   z = FpX_mul(v,g, pd);
     106                 :     145736 :   t = FpX_divrem(z,a, pd, &s);
     107                 :     145736 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     108                 :     145736 :   t = FpX_red(t, pd);
     109                 :     145736 :   t = ZX_Z_mul(t,p0);
     110                 :     145736 :   s = ZX_Z_mul(s,p0);
     111                 :     145736 :   avma = av;
     112                 :     145736 :   a2 = ZX_add(a,s);
     113                 :     145736 :   b2 = ZX_add(b,t);
     114                 :            : 
     115                 :            :   /* already reduced mod p1 = pd p0 */
     116                 :     145736 :   gel(V,j)   = a2;
     117                 :     145736 :   gel(V,j+1) = b2;
     118         [ +  + ]:     291472 :   if (noinv) return;
     119                 :            : 
     120                 :     123546 :   av = avma;
     121                 :     123546 :   (void)new_chunk(space); /* HACK */
     122                 :     123546 :   g = ZX_add(ZX_mul(u,a2), ZX_mul(v,b2));
     123                 :     123546 :   g = Z_ZX_sub(gen_1, g);
     124                 :     123546 :   g = ZX_Z_divexact(g, p0);
     125                 :     123546 :   g = FpX_red(g, pd);
     126                 :     123546 :   z = FpX_mul(v,g, pd);
     127                 :     123546 :   t = FpX_divrem(z,a, pd, &s);
     128                 :     123546 :   t = ZX_add(ZX_mul(u,g), ZX_mul(t,b));
     129                 :     123546 :   t = FpX_red(t, pd);
     130                 :     123546 :   t = ZX_Z_mul(t,p0);
     131                 :     123546 :   s = ZX_Z_mul(s,p0);
     132                 :     123546 :   avma = av;
     133                 :     123546 :   gel(W,j)   = ZX_add(u,t);
     134                 :     145736 :   gel(W,j+1) = ZX_add(v,s);
     135                 :            : }
     136                 :            : 
     137                 :            : static void
     138                 :      21455 : ZpXQ_HenselLift(GEN V, GEN W, long j, GEN f, GEN Td, GEN T1, GEN pd, GEN p0, GEN p1, int noinv)
     139                 :            : {
     140                 :      21455 :   pari_sp av = avma;
     141                 :      21455 :   const long n = degpol(T1), vT = varn(T1);
     142                 :      21455 :   long space = lg(f) * lgefint(p1) * lg(T1);
     143                 :            :   GEN a2, b2, g, z, s, t;
     144                 :      21455 :   GEN a = gel(V,j), b = gel(V,j+1);
     145                 :      21455 :   GEN u = gel(W,j), v = gel(W,j+1);
     146                 :            : 
     147                 :      21455 :   (void)new_chunk(space); /* HACK */
     148                 :      21455 :   g = RgX_sub(f, Kronecker_to_ZXX(ZXX_mul_Kronecker(a,b,n), n, vT));
     149                 :      21455 :   g = FpXQX_red(g, T1, p1);
     150                 :      21455 :   g = RgX_Rg_divexact(g, p0);
     151                 :      21455 :   z = FpXQX_mul(v,g, Td,pd);
     152                 :      21455 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     153                 :      21455 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     154                 :      21455 :   t = Kronecker_to_ZXX(t, n, vT);
     155                 :      21455 :   t = FpXQX_red(t, Td, pd);
     156                 :      21455 :   t = RgX_Rg_mul(t,p0);
     157                 :      21455 :   s = RgX_Rg_mul(s,p0);
     158                 :      21455 :   avma = av;
     159                 :            : 
     160                 :      21455 :   a2 = RgX_add(a,s);
     161                 :      21455 :   b2 = RgX_add(b,t);
     162                 :            :   /* already reduced mod p1 = pd p0 */
     163                 :      21455 :   gel(V,j)   = a2;
     164                 :      21455 :   gel(V,j+1) = b2;
     165         [ +  + ]:      42910 :   if (noinv) return;
     166                 :            : 
     167                 :      16625 :   av = avma;
     168                 :      16625 :   (void)new_chunk(space); /* HACK */
     169                 :      16625 :   g = ZX_add(ZXX_mul_Kronecker(u,a2,n), ZXX_mul_Kronecker(v,b2,n));
     170                 :      16625 :   g = Kronecker_to_ZXX(g, n, vT);
     171                 :      16625 :   g = Rg_RgX_sub(gen_1, g);
     172                 :      16625 :   g = FpXQX_red(g, T1, p1);
     173                 :      16625 :   g = RgX_Rg_divexact(g, p0);
     174                 :      16625 :   z = FpXQX_mul(v,g, Td,pd);
     175                 :      16625 :   t = FpXQX_divrem(z,a, Td,pd, &s);
     176                 :      16625 :   t = ZX_add(ZXX_mul_Kronecker(u,g,n), ZXX_mul_Kronecker(t,b,n));
     177                 :      16625 :   t = Kronecker_to_ZXX(t, n, vT);
     178                 :      16625 :   t = FpXQX_red(t, Td, pd);
     179                 :      16625 :   t = RgX_Rg_mul(t,p0);
     180                 :      16625 :   s = RgX_Rg_mul(s,p0);
     181                 :      16625 :   avma = av;
     182                 :      16625 :   gel(W,j)   = RgX_add(u,t);
     183                 :      21455 :   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                 :     322948 : ZpX_RecTreeLift(GEN link, GEN v, GEN w, GEN pd, GEN p0, GEN p1,
     190                 :            :                 GEN f, long j, int noinv)
     191                 :            : {
     192         [ +  + ]:     468684 :   if (j < 0) return;
     193                 :     145736 :   ZpX_HenselLift(v, w, j, f, pd, p0,p1, noinv);
     194                 :     145736 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     195                 :     145736 :   ZpX_RecTreeLift(link, v, w, pd, p0,p1, gel(v,j+1), link[j+1], noinv);
     196                 :            : }
     197                 :            : static void
     198                 :      45528 : 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         [ +  + ]:      66983 :   if (j < 0) return;
     202                 :      21455 :   ZpXQ_HenselLift(v, w, j, f, Td,T1, pd, p0,p1, noinv);
     203                 :      21455 :   ZpXQ_RecTreeLift(link, v, w, Td,T1, pd, p0,p1, gel(v,j)  , link[j  ], noinv);
     204                 :      21455 :   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                 :     225147 : quadratic_prec_mask(long n)
     225                 :            : {
     226                 :     225147 :   long a = n, i;
     227                 :     225147 :   ulong mask = 0;
     228                 :     225147 :   for(i = 1;; i++, mask <<= 1)
     229                 :            :   {
     230                 :     518520 :     mask |= (a&1); a = (a+1)>>1;
     231         [ +  + ]:     518520 :     if (a==1) return mask | (1UL << i);
     232                 :     293373 :   }
     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                 :       7021 : MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
     242                 :            : {
     243                 :       7021 :   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         [ -  + ]:       7021 :   if (k < 2) pari_err_DOMAIN("MultiLift", "#(modular factors)", "<", gen_2,a);
     249         [ -  + ]:       7021 :   if (e0 < 1) pari_err_DOMAIN("MultiLift", "precision", "<", gen_1,stoi(e0));
     250         [ +  + ]:       7021 :   if (e0 == 1) return a;
     251                 :            : 
     252         [ -  + ]:       7014 :   if (DEBUGLEVEL > 3) timer_start(&Ti);
     253         [ -  + ]:       7014 :   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                 :       7014 :     e = 1;
     263                 :       7014 :     v = cgetg(2*k-2 + 1, t_VEC);
     264                 :       7014 :     w = cgetg(2*k-2 + 1, t_VEC);
     265                 :       7014 :     link=cgetg(2*k-2 + 1, t_VECSMALL);
     266         [ +  + ]:       7014 :     BuildTree(link, v, w, a, T? FpX_red(T,p): NULL, p);
     267         [ -  + ]:       7007 :     if (DEBUGLEVEL > 3) timer_printf(&Ti, "building tree");
     268                 :            :   }
     269                 :       7007 :   mask = quadratic_prec_mask(e0);
     270                 :       7007 :   eold = 1;
     271                 :       7007 :   penew = NULL;
     272                 :       7007 :   Tnew = NULL;
     273         [ +  + ]:      41101 :   while (mask > 1)
     274                 :            :   {
     275                 :      34094 :     long enew = eold << 1;
     276         [ +  + ]:      34094 :     if (mask & 1) enew--;
     277                 :      34094 :     mask >>= 1;
     278         [ +  - ]:      34094 :     if (enew >= e) { /* mask == 1: last iteration */
     279         [ +  + ]:      34094 :       GEN peold = penew? penew: powiu(p, eold);
     280                 :      34094 :       GEN Td = NULL, pd;
     281                 :      34094 :       long d = enew - eold; /* = eold or eold-1 */
     282                 :            :       /* lift from p^eold to p^enew */
     283         [ +  + ]:      34094 :       pd = (d == eold)? peold: diviiexact(peold, p); /* p^d */
     284                 :      34094 :       penew = mulii(peold,pd);
     285         [ +  + ]:      34094 :       if (T) {
     286         [ +  + ]:       2618 :         if (Tnew)
     287         [ +  + ]:       1904 :           Td = (d == eold)? Tnew: FpX_red(Tnew,pd);
     288                 :            :         else
     289                 :        714 :           Td = FpX_red(T, peold);
     290                 :       2618 :         Tnew = FpX_red(T, penew);
     291 [ +  - ][ +  + ]:       2618 :         ZpXQ_RecTreeLift(link, v, w, Td, Tnew, pd, peold, penew, f, lg(v)-2,
     292                 :            :                          (flag == 0 && mask == 1));
     293                 :            :       }
     294                 :            :       else
     295 [ +  + ][ +  + ]:      31476 :         ZpX_RecTreeLift(link, v, w, pd, peold, penew, f, lg(v)-2,
     296                 :            :                         (flag == 0 && mask == 1));
     297         [ -  + ]:      34094 :       if (DEBUGLEVEL > 3) timer_printf(&Ti, "lifting to prec %ld", enew);
     298                 :            :     }
     299                 :      34094 :     eold = enew;
     300                 :            :   }
     301                 :            : 
     302         [ +  + ]:       7007 :   if (flag)
     303                 :        525 :     E = mkvec4(utoipos(e0), link, v, w);
     304                 :            :   else
     305                 :            :   {
     306                 :       6482 :     E = cgetg(k+1, t_VEC);
     307         [ +  + ]:      60522 :     for (i = 1; i <= 2*k-2; i++)
     308                 :            :     {
     309                 :      54040 :       long t = link[i];
     310         [ +  + ]:      54040 :       if (t < 0) gel(E,-t) = gel(v,i);
     311                 :            :     }
     312                 :            :   }
     313                 :       7014 :   return E;
     314                 :            : }
     315                 :            : 
     316                 :            : /* Q list of (coprime, monic) factors of pol mod (T,p). Lift mod p^e = pe.
     317                 :            :  * T may be NULL */
     318                 :            : GEN
     319                 :       6685 : ZpX_liftfact(GEN pol, GEN Q, GEN T, GEN p, long e, GEN pe)
     320                 :            : {
     321                 :       6685 :   pari_sp av = avma;
     322         [ +  + ]:       6685 :   if (lg(Q) == 2) return mkvec(pol);
     323                 :       6496 :   pol = FqX_normalize(pol, T, pe);
     324                 :       6685 :   return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
     325                 :            : }
     326                 :            : 
     327                 :            : /* U = NULL treated as 1 */
     328                 :            : static void
     329                 :       5369 : BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
     330                 :            : {
     331                 :            :   GEN Q, R;
     332         [ +  + ]:       7791 :   if (j < 0) return;
     333                 :            : 
     334                 :       2422 :   Q = FpX_mul(gel(v,j), gel(w,j), pe);
     335         [ +  + ]:       2422 :   if (U)
     336                 :            :   {
     337                 :       1897 :     Q = FpXQ_mul(Q, U, f, pe);
     338                 :       1897 :     R = FpX_sub(U, Q, pe);
     339                 :            :   }
     340                 :            :   else
     341                 :        525 :     R = Fp_FpX_sub(gen_1, Q, pe);
     342                 :       2422 :   gel(w,j+1) = Q; /*  0 mod U v[j],  1 mod (1-U) v[j+1] */
     343                 :       2422 :   gel(w,j) = R; /*  1 mod U v[j],  0 mod (1-U) v[j+1] */
     344                 :       2422 :   BezoutPropagate(link, v, w, pe, R, f, link[j  ]);
     345                 :       2422 :   BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
     346                 :            : }
     347                 :            : 
     348                 :            : /* as above, but return the Bezout coefficients for the lifted modular factors
     349                 :            :  *   U[i] = 1 mod Qlift[i]
     350                 :            :  *          0 mod Qlift[j], j != i */
     351                 :            : GEN
     352                 :        532 : bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
     353                 :            : {
     354                 :        532 :   pari_sp av = avma;
     355                 :            :   GEN E, link, v, w, pe;
     356                 :        532 :   long i, k = lg(Q)-1;
     357         [ +  + ]:        532 :   if (k == 1) return mkvec(pol);
     358                 :        525 :   pe = powiu(p, e);
     359                 :        525 :   pol = FpX_normalize(pol, pe);
     360                 :        525 :   E = MultiLift(pol, Q, NULL, p, e, 1);
     361                 :        525 :   link = gel(E,2);
     362                 :        525 :   v    = gel(E,3);
     363                 :        525 :   w    = gel(E,4);
     364                 :        525 :   BezoutPropagate(link, v, w, pe, NULL, pol, lg(v)-2);
     365                 :        525 :   E = cgetg(k+1, t_VEC);
     366         [ +  + ]:       5369 :   for (i = 1; i <= 2*k-2; i++)
     367                 :            :   {
     368                 :       4844 :     long t = link[i];
     369         [ +  + ]:       4844 :     if (t < 0) E[-t] = w[i];
     370                 :            :   }
     371                 :        532 :   return gerepilecopy(av, E);
     372                 :            : }
     373                 :            : 
     374                 :            : /* Front-end for ZpX_liftfact:
     375                 :            :    lift the factorization of pol mod p given by L to p^N (if possible) */
     376                 :            : GEN
     377                 :         21 : polhensellift(GEN pol, GEN L, GEN p, long N)
     378                 :            : {
     379                 :         21 :   GEN T = NULL;
     380                 :            :   long i, l, t;
     381                 :         21 :   pari_sp av = avma;
     382                 :            : 
     383         [ -  + ]:         21 :   if (typ(pol) != t_POL) pari_err_TYPE("polhensellift",pol);
     384                 :         21 :   RgX_check_ZXX(pol, "polhensellift");
     385 [ +  - ][ -  + ]:         21 :   if (!is_vec_t(typ(L)) || lg(L) < 3) pari_err_TYPE("polhensellift",L);
     386                 :         21 :   t = typ(p);
     387         [ +  + ]:         21 :   if (t == t_VEC) /* [p, T] */
     388                 :            :   {
     389                 :          7 :     T = gel(p,2);
     390         [ -  + ]:          7 :     if (typ(T) != t_POL) pari_err_TYPE("polhensellift",pol);
     391                 :          7 :     RgX_check_ZX(T, "polhensellift");
     392                 :          7 :     p = gel(p,1); t = typ(p);
     393                 :            :   }
     394         [ -  + ]:         21 :   if (t != t_INT) pari_err_TYPE("polhensellift",p);
     395         [ -  + ]:         21 :   if (N < 1) pari_err_DOMAIN("polhensellift", "precision", "<", gen_1,stoi(N));
     396                 :            : 
     397                 :         21 :   l = lg(L); L = leafcopy(L);
     398         [ +  + ]:         63 :   for (i = 1; i < l; i++)
     399                 :            :   {
     400         [ -  + ]:         42 :     if (typ(gel(L,i)) != t_POL)
     401                 :          0 :       gel(L,i) = scalar_ZX_shallow(gel(L,i), varn(pol));
     402                 :         42 :     RgX_check_ZXX(gel(L,i), "polhensellift");
     403                 :            :   }
     404                 :         21 :   return gerepilecopy(av, ZpX_liftfact(pol, L, T, p, N, powiu(p,N)));
     405                 :            : }
     406                 :            : 
     407                 :            : static GEN
     408                 :         35 : ZpX_liftroots_full(GEN f, GEN S, GEN p, long e)
     409                 :            : {
     410                 :         35 :   long i, n = lg(S)-1;
     411                 :         35 :   GEN r = cgetg(n+1, typ(S));
     412                 :         35 :   pari_sp av = avma;
     413                 :         35 :   long v = varn(f);
     414                 :            :   GEN y, q;
     415         [ +  + ]:        833 :   for (i=1; i<=n; i++)
     416                 :        798 :     gel(r,i) = deg1pol(gen_1, Fp_neg(gel(S, i), p), v);
     417                 :         35 :   q = powiu(p, e);
     418                 :         35 :   y = ZpX_liftfact(f, r, NULL, p, e, q);
     419                 :         35 :   r = cgetg(n+1 ,t_COL);
     420         [ +  + ]:        833 :   for (i=1; i<=n; i++)
     421                 :        798 :     gel(r,i) = Fp_neg(gmael(y, i, 2), q);
     422                 :         35 :   return gerepileupto(av, r);
     423                 :            : }
     424                 :            : 
     425                 :            : GEN
     426                 :      39949 : ZpX_roots(GEN F, GEN p, long e)
     427                 :            : {
     428                 :      39949 :   pari_sp av = avma;
     429                 :      39949 :   long i, v = varn(F);
     430                 :            :   GEN y, r, q;
     431                 :      39949 :   GEN f = FpX_normalize(F, p);
     432                 :      39949 :   GEN g = FpX_normalize(FpX_split_part(f, p), p);
     433                 :      39949 :   GEN S = FpX_roots(g, p);
     434         [ +  + ]:      39949 :   long l = lg(S)-1, n = l < degpol(f)? l+1: l;
     435         [ -  + ]:      39949 :   if (l==0) return cgetg(1, t_COL);
     436         [ +  + ]:      39949 :   if (l==1) return mkcol(ZpX_liftroot(F,gel(S,1),p,e));
     437                 :       1967 :   r = cgetg(n+1 ,t_COL);
     438         [ +  + ]:      17990 :   for (i=1; i<=l; i++)
     439                 :      16023 :     gel(r,i) = deg1pol_shallow(gen_1, Fp_neg(gel(S,i), p), v);
     440         [ +  + ]:       1967 :   if (l < degpol(f)) gel(r, n) = FpX_div(f, g, p);
     441                 :       1967 :   q = powiu(p, e);
     442                 :       1967 :   y = ZpX_liftfact(F, r, NULL, p, e, q);
     443                 :       1967 :   r = cgetg(l+1 ,t_COL);
     444         [ +  + ]:      17990 :   for (i=1; i<=l; i++) gel(r,i) = Fp_neg(gmael(y,i,2), q);
     445                 :      39949 :   return gerepileupto(av, r);
     446                 :            : }
     447                 :            : 
     448                 :            : /* SPEC:
     449                 :            :  * p is a t_INT > 1, e >= 1
     450                 :            :  * f is a ZX with leading term prime to p.
     451                 :            :  * a is a simple root mod l for all l|p.
     452                 :            :  * Return roots of f mod p^e, as integers (implicitly mod p^e)
     453                 :            :  * STANDARD USE: p is a prime power */
     454                 :            : GEN
     455                 :      38241 : ZpX_liftroot(GEN f, GEN a, GEN p, long e)
     456                 :            : {
     457                 :      38241 :   pari_sp av = avma;
     458                 :      38241 :   GEN q = p, fr, W;
     459                 :            :   ulong mask;
     460                 :            : 
     461                 :      38241 :   a = modii(a,q);
     462         [ +  + ]:      38241 :   if (e == 1) return a;
     463                 :      38227 :   mask = quadratic_prec_mask(e);
     464                 :      38227 :   fr = FpX_red(f,q);
     465                 :      38227 :   W = Fp_inv(FpX_eval(ZX_deriv(fr), a, q), q); /* 1/f'(a) mod p */
     466                 :            :   for(;;)
     467                 :            :   {
     468                 :     145536 :     q = sqri(q);
     469         [ +  + ]:     145536 :     if (mask & 1) q = diviiexact(q, p);
     470                 :     145536 :     mask >>= 1;
     471                 :     145536 :     fr = FpX_red(f,q);
     472                 :     145536 :     a = Fp_sub(a, Fp_mul(W, FpX_eval(fr, a,q), q), q);
     473         [ +  + ]:     145536 :     if (mask == 1) return gerepileuptoint(av, a);
     474                 :     107309 :     W = Fp_sub(shifti(W,1), Fp_mul(Fp_sqr(W,q), FpX_eval(ZX_deriv(fr),a,q), q), q);
     475                 :     145550 :   }
     476                 :            : }
     477                 :            : 
     478                 :            : GEN
     479                 :         35 : ZpX_liftroots(GEN f, GEN S, GEN p, long e)
     480                 :            : {
     481                 :         35 :   long i, n = lg(S)-1, d = degpol(f);
     482                 :            :   GEN r;
     483         [ +  - ]:         35 :   if (n == d) return ZpX_liftroots_full(f, S, p, e);
     484                 :          0 :   r = cgetg(n+1, typ(S));
     485         [ #  # ]:          0 :   for (i=1; i <= n; i++)
     486                 :          0 :     gel(r,i) = ZpX_liftroot(f, gel(S,i), p, e);
     487                 :         35 :   return r;
     488                 :            : }
     489                 :            : 
     490                 :            : GEN
     491                 :         91 : ZpXQX_liftroot_vald(GEN f, GEN a, long v, GEN T, GEN p, long e)
     492                 :            : {
     493                 :         91 :   pari_sp av = avma, av2;
     494                 :         91 :   GEN pv = p, q, qv, W, df, Tq, fr, dfr;
     495                 :            :   ulong mask;
     496                 :         91 :   a = Fq_red(a, T, p);
     497         [ -  + ]:         91 :   if (e <= v+1) return a;
     498                 :         91 :   df = RgX_deriv(f);
     499         [ -  + ]:         91 :   if (v) { pv = powiu(p,v); qv = mulii(pv,p); df = ZXX_Z_divexact(df, pv); }
     500                 :         91 :   else qv = p;
     501                 :         91 :   mask = quadratic_prec_mask(e-v);
     502                 :         91 :   Tq = FpXT_red(T, qv); dfr = FpXQX_red(df, Tq, p);
     503                 :         91 :   W = Fq_inv(FqX_eval(dfr, a, Tq, p), Tq, p); /* 1/f'(a) mod (T,p) */
     504                 :         91 :   q = p;
     505                 :         91 :   av2 = avma;
     506                 :            :   for (;;)
     507                 :            :   {
     508                 :            :     GEN u, fa, qv, q2v, q2, Tq2;
     509                 :        161 :     q2 = q; q = sqri(q);
     510         [ +  + ]:        161 :     if (mask & 1) q = diviiexact(q,p);
     511                 :        161 :     mask >>= 1;
     512         [ -  + ]:        161 :     if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
     513                 :        161 :     else { qv = q; q2v = q2; }
     514                 :        161 :     Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
     515                 :        161 :     fr = FpXQX_red(f, Tq, qv);
     516                 :        161 :     fa = FqX_eval(fr, a, Tq, qv);
     517         [ -  + ]:        161 :     fa = typ(fa)==t_INT? diviiexact(fa,q2v): ZX_Z_divexact(fa, q2v);
     518                 :        161 :     a = Fq_sub(a, ZX_Z_mul(Fq_mul(W, fa, Tq2, q2v), q2), Tq, qv);
     519         [ +  + ]:        161 :     if (mask == 1) return gerepileupto(av, a);
     520                 :         70 :     dfr = FpXQX_red(df, Tq, q);
     521                 :         70 :     u = ZX_Z_divexact(FpX_Fp_sub(Fq_mul(W,FqX_eval(dfr,a,Tq,q),Tq,q),gen_1,q),q2);
     522                 :         70 :     W = Fq_sub(W,ZX_Z_mul(Fq_mul(u,W,Tq2,q2),q2),Tq,q);
     523         [ -  + ]:         70 :     if (gc_needed(av2,2))
     524                 :            :     {
     525         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_liftroot, e = %ld", e);
     526                 :          0 :       gerepileall(av2, 3, &a, &W, &q);
     527                 :            :     }
     528                 :        161 :   }
     529                 :            : }
     530                 :            : 
     531                 :            : GEN
     532                 :         91 : ZpXQX_liftroot(GEN f, GEN a, GEN T, GEN p, long e) { return ZpXQX_liftroot_vald(f,a,0,T,p,e); }
     533                 :            : 
     534                 :            : /* Same as ZpX_liftroot for the polynomial X^n-b*/
     535                 :            : GEN
     536                 :      28105 : Zp_sqrtnlift(GEN b, GEN n, GEN a, GEN p, long e)
     537                 :            : {
     538                 :      28105 :   pari_sp ltop=avma;
     539                 :            :   GEN q, w, n_1;
     540                 :            :   ulong mask;
     541                 :      28105 :   long pis2 = equalii(n, gen_2)? 1: 0;
     542         [ +  + ]:      28105 :   if (e == 1) return icopy(a);
     543                 :      25893 :   n_1 = subis(n,1);
     544                 :      25893 :   mask = quadratic_prec_mask(e);
     545         [ +  + ]:      25893 :   w = Fp_inv(pis2 ? shifti(a,1): Fp_mul(n,Fp_pow(a,n_1,p), p), p);
     546                 :      25893 :   q = p;
     547                 :            :   for(;;)
     548                 :            :   {
     549                 :      57442 :     q = sqri(q);
     550         [ +  + ]:      57442 :     if (mask & 1) q = diviiexact(q, p);
     551                 :      57442 :     mask >>= 1;
     552 [ +  + ][ +  - ]:      57442 :     if (lgefint(q) == 3 && lgefint(n) == 3)
     553                 :      31413 :     {
     554                 :      57059 :       ulong Q = uel(q,2), N = uel(n,2);
     555                 :      57059 :       ulong A = umodiu(a, Q);
     556                 :      57059 :       ulong B = umodiu(b, Q);
     557                 :      57059 :       ulong W = umodiu(w, Q);
     558                 :            : 
     559                 :      57059 :       A = Fl_sub(A, Fl_mul(W, Fl_sub(Fl_powu(A,N,Q), B, Q), Q), Q);
     560                 :      57059 :       a = utoi(A);
     561         [ +  + ]:      57059 :       if (mask == 1) break;
     562                 :      31413 :       W = Fl_sub(Fl_add(W,W,Q),
     563                 :            :                  Fl_mul(Fl_sqr(W,Q), Fl_mul(N,Fl_powu(A, N-1, Q), Q), Q), Q);
     564                 :      31413 :       w = utoi(W);
     565                 :            :     }
     566                 :            :     else
     567                 :            :     {
     568                 :            :       /* a -= w (a^n - b) */
     569                 :        383 :       a = modii(subii(a, mulii(w, subii(Fp_pow(a,n,q),b))), q);
     570         [ +  + ]:        383 :       if (mask == 1) break;
     571                 :            :       /* w += w - w^2 n a^(n-1)*/
     572         [ +  + ]:        257 :       w = subii(shifti(w,1), Fp_mul(Fp_sqr(w,q),
     573                 :        121 :                            pis2? shifti(a,1): mulii(n,Fp_pow(a,n_1,q)), q));
     574                 :            :     }
     575                 :      31549 :   }
     576                 :      28105 :   return gerepileuptoint(ltop,a);
     577                 :            : }
     578                 :            : 
     579                 :            : 
     580                 :            : /* Same as ZpX_liftroot for the polynomial X^2-b */
     581                 :            : GEN
     582                 :       7028 : Zp_sqrtlift(GEN b, GEN a, GEN p, long e)
     583                 :            : {
     584                 :       7028 :   return Zp_sqrtnlift(b, gen_2, a, p, e);
     585                 :            : }
     586                 :            : 
     587                 :            : GEN
     588                 :       1568 : Zp_sqrt(GEN x, GEN p, long e)
     589                 :            : {
     590                 :       1568 :   pari_sp av = avma;
     591                 :       1568 :   GEN z = Fp_sqrt(Fp_red(x, p), p);
     592         [ +  + ]:       1568 :   if (!z) return NULL;
     593         [ +  + ]:       1547 :   if (e > 1) z = Zp_sqrtlift(x, z, p, e);
     594                 :       1568 :   return gerepileuptoint(av, z);
     595                 :            : }
     596                 :            : 
     597                 :            : /* Compute (x-1)/(x+1)/p^k */
     598                 :            : static GEN
     599                 :      10262 : ZpXQ_log_to_ath(GEN x, long k, GEN T, GEN p, long e, GEN pe)
     600                 :            : {
     601                 :      10262 :   pari_sp av = avma;
     602                 :      10262 :   long vT = get_FpX_var(T);
     603                 :            :   GEN bn, bdi;
     604                 :      10262 :   GEN bd = ZX_Z_add(x, gen_1);
     605         [ +  + ]:      10262 :   if (equaliu(p,2)) /*For p=2, we need to simplify by 2*/
     606                 :            :   {
     607                 :       3332 :     bn = ZX_shifti(x,-(k+1));
     608                 :       3332 :     bdi= ZpXQ_invlift(ZX_shifti(bd ,-1), pol_1(vT), T, p, e);
     609                 :            :   }
     610                 :            :   else
     611                 :            :   {
     612                 :       6930 :     bn = ZX_Z_divexact(ZX_Z_sub(x, gen_1),powiu(p,k));
     613                 :       6930 :     bdi= ZpXQ_invlift(bd, scalarpol(Fp_inv(gen_2,p),vT), T, p, e);
     614                 :            :   }
     615                 :      10262 :   return gerepileupto(av, FpXQ_mul(bn, bdi, T, pe));
     616                 :            : }
     617                 :            : 
     618                 :            : /* Assume p odd, a = 1 [p], return log(a) */
     619                 :            : GEN
     620                 :      10262 : ZpXQ_log(GEN a, GEN T, GEN p, long N)
     621                 :            : {
     622                 :      10262 :   pari_sp av = avma;
     623                 :            :   pari_timer ti;
     624                 :      10262 :   long is2 = equaliu(p,2);
     625         [ +  + ]:      10262 :   ulong pp = is2 ? 0: itou_or_0(p);
     626 [ +  + ][ +  - ]:      10262 :   double lp = is2 ? 1: pp ? log2(pp): expi(p);
     627                 :      10262 :   long k = maxss(1 , (long) .5+pow((double)(N>>1)/(lp*lp), 1./3));
     628                 :            :   GEN ak, s, b, pol;
     629         [ +  + ]:      10262 :   long e = is2 ? N-1: N;
     630                 :      10262 :   long i, l = (e-2)/(2*(k+is2));
     631                 :      10262 :   GEN pe = powiu(p,e);
     632                 :      10262 :   GEN TNk, pNk = powiu(p,N+k);
     633         [ -  + ]:      10262 :   if( DEBUGLEVEL>=3) timer_start(&ti);
     634                 :      10262 :   TNk = FpX_get_red(get_FpX_mod(T), pNk);
     635                 :      10262 :   ak = FpXQ_pow(a, powiu(p,k), TNk, pNk);
     636         [ -  + ]:      10262 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpXQ_pow(%ld)",k);
     637                 :      10262 :   b = ZpXQ_log_to_ath(ak, k, T, p, e, pe);
     638         [ -  + ]:      10262 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"ZpXQ_log_to_ath");
     639                 :      10262 :   pol= cgetg(l+3,t_POL);
     640                 :      10262 :   pol[1] = evalsigne(1)|evalvarn(0);
     641         [ +  + ]:      26852 :   for(i=0; i<=l; i++)
     642                 :            :   {
     643                 :            :     GEN g;
     644                 :      16590 :     ulong z = 2*i+1;
     645         [ +  + ]:      16590 :     if (pp)
     646                 :            :     {
     647                 :      13132 :       long w = u_lvalrem(z, pp, &z);
     648                 :      13132 :       g = powuu(pp,2*i*k-w);
     649                 :            :     }
     650                 :       3458 :     else g = powiu(p,2*i*k);
     651                 :      16590 :     gel(pol,i+2) = Fp_div(g, utoi(z),pe);
     652                 :            :   }
     653         [ -  + ]:      10262 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"pol(%ld)",l);
     654                 :      10262 :   s = FpX_FpXQ_eval(pol, FpXQ_sqr(b, T, pe), T,  pe);
     655         [ -  + ]:      10262 :   if( DEBUGLEVEL>=3) timer_printf(&ti,"FpX_FpXQ_eval");
     656                 :      10262 :   s = ZX_shifti(FpXQ_mul(b, s, T, pe), 1);
     657         [ +  + ]:      10262 :   return gerepileupto(av, is2? s: FpX_red(s, pe));
     658                 :            : }
     659                 :            : 
     660                 :            : /***********************************************************************/
     661                 :            : /**                                                                   **/
     662                 :            : /**                 Generic quadratic hensel lift over Zp[X]          **/
     663                 :            : /**                                                                   **/
     664                 :            : /***********************************************************************/
     665                 :            : /* q = p^N */
     666                 :            : 
     667                 :            : GEN
     668                 :      99918 : gen_ZpX_Dixon(GEN F, GEN V, GEN q, GEN p, long N, void *E,
     669                 :            :                             GEN lin(void *E, GEN F, GEN d, GEN q),
     670                 :            :                             GEN invl(void *E, GEN d))
     671                 :            : {
     672                 :      99918 :   pari_sp av = avma;
     673                 :            :   long N2, M;
     674                 :            :   GEN VN2, V2, VM, bil;
     675                 :            :   GEN q2, qM;
     676                 :      99918 :   V = FpX_red(V, q);
     677         [ +  + ]:      99918 :   if (N == 1) return invl(E, V);
     678                 :      22869 :   N2 = (N + 1)>>1; M = N - N2;
     679                 :      22869 :   F = FpXT_red(F, q);
     680                 :      22869 :   qM = powiu(p, M);
     681         [ +  + ]:      22869 :   q2 = M == N2? qM: mulii(qM, p);
     682                 :            :   /* q2 = p^N2, qM = p^M, q = q2 * qM */
     683                 :      22869 :   VN2 = gen_ZpX_Dixon(F, V, q2, p, N2, E, lin, invl);
     684                 :      22869 :   bil = lin(E, F, VN2, q);
     685                 :      22869 :   V2 = ZX_Z_divexact(ZX_sub(V, bil), q2);
     686                 :      22869 :   VM = gen_ZpX_Dixon(F, V2, qM, p, M, E, lin, invl);
     687                 :      99918 :   return gerepileupto(av, FpX_red(ZX_add(VN2, ZX_Z_mul(VM, q2)), q));
     688                 :            : }
     689                 :            : 
     690                 :            : GEN
     691                 :      56819 : gen_ZpX_Newton(GEN x, GEN p, long n, void *E,
     692                 :            :                       GEN eval(void *E, GEN f, GEN q),
     693                 :            :                       GEN invd(void *E, GEN V, GEN v, GEN q, long M))
     694                 :            : {
     695                 :      56819 :   pari_sp ltop = avma, av;
     696                 :      56819 :   long N = 1, N2, M;
     697                 :            :   long mask;
     698                 :      56819 :   GEN q = p;
     699         [ +  + ]:      56819 :   if (n == 1) return gcopy(x);
     700                 :      56301 :   mask = quadratic_prec_mask(n);
     701                 :      56301 :   av = avma;
     702         [ +  + ]:     180684 :   while (mask > 1)
     703                 :            :   {
     704                 :            :     GEN qM, q2, v, V;
     705                 :     124383 :     N2 = N; N <<= 1;
     706                 :     124383 :     q2 = q;
     707         [ +  + ]:     124383 :     if (mask&1UL) { /* can never happen when q2 = p */
     708                 :      32424 :       N--; M = N2-1;
     709                 :      32424 :       qM = diviiexact(q2,p); /* > 1 */
     710                 :      32424 :       q = mulii(qM,q2);
     711                 :            :     } else {
     712                 :      91959 :       M = N2;
     713                 :      91959 :       qM = q2;
     714                 :      91959 :       q = sqri(q2);
     715                 :            :     }
     716                 :            :     /* q2 = p^N2, qM = p^M, q = p^N = q2 * qM */
     717                 :     124383 :     mask >>= 1;
     718                 :     124383 :     v = eval(E, x, q);
     719                 :     124383 :     V = ZX_Z_divexact(gel(v,1), q2);
     720                 :     124383 :     x = FpX_sub(x, ZX_Z_mul(invd(E, V, v, qM, M), q2), q);
     721         [ -  + ]:     124383 :     if (gc_needed(av, 1))
     722                 :            :     {
     723         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gen_ZpX_Newton");
     724                 :          0 :       gerepileall(av, 2, &x, &q);
     725                 :            :     }
     726                 :            :   }
     727                 :      56819 :   return gerepileupto(ltop, x);
     728                 :            : }
     729                 :            : 
     730                 :            : struct _ZpXQ_inv
     731                 :            : {
     732                 :            :   GEN T, a, p ,n;
     733                 :            : };
     734                 :            : 
     735                 :            : static GEN
     736                 :      54691 : _inv_invd(void *E, GEN V, GEN v, GEN q, long M/*unused*/)
     737                 :            : {
     738                 :      54691 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
     739                 :      54691 :   GEN Tq = FpXT_red(d->T, q);
     740                 :            :   (void)M;
     741                 :      54691 :   return FpXQ_mul(V, gel(v,2), Tq, q);
     742                 :            : }
     743                 :            : 
     744                 :            : static GEN
     745                 :      54691 : _inv_eval(void *E, GEN x, GEN q)
     746                 :            : {
     747                 :      54691 :   struct _ZpXQ_inv *d = (struct _ZpXQ_inv *) E;
     748                 :      54691 :   GEN Tq = FpXT_red(d->T, q);
     749                 :      54691 :   GEN f = FpX_Fp_sub(FpXQ_mul(x, FpX_red(d->a, q), Tq, q), gen_1, q);
     750                 :      54691 :   return mkvec2(f, x);
     751                 :            : }
     752                 :            : 
     753                 :            : GEN
     754                 :      24360 : ZpXQ_invlift(GEN a, GEN x, GEN T, GEN p, long e)
     755                 :            : {
     756                 :            :   struct _ZpXQ_inv d;
     757                 :      24360 :   d.a = a; d.T = T; d.p = p;
     758                 :      24360 :   return gen_ZpX_Newton(x, p, e, &d, _inv_eval, _inv_invd);
     759                 :            : }
     760                 :            : 
     761                 :            : GEN
     762                 :      14098 : ZpXQ_inv(GEN a, GEN T, GEN p, long e)
     763                 :            : {
     764                 :      14098 :   pari_sp av=avma;
     765                 :            :   GEN ai;
     766         [ +  + ]:      14098 :   if (lgefint(p)==3)
     767                 :            :   {
     768                 :      14070 :     ulong pp = p[2];
     769                 :      14070 :     ai = Flx_to_ZX(Flxq_inv(ZX_to_Flx(a,pp), ZXT_to_FlxT(T, pp), pp));
     770                 :            :   } else
     771                 :         28 :     ai = FpXQ_inv(FpX_red(a,p), FpXT_red(T,p),p);
     772                 :      14098 :   return gerepileupto(av, ZpXQ_invlift(a, ai, T, p, e));
     773                 :            : }
     774                 :            : 
     775                 :            : struct _ZpXQ_sqrtn
     776                 :            : {
     777                 :            :   GEN T, a, n, ai;
     778                 :            : };
     779                 :            : 
     780                 :            : static GEN
     781                 :       1778 : _sqrtn_invd(void *E, GEN V, GEN v, GEN q, long M)
     782                 :            : {
     783                 :       1778 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
     784                 :       1778 :   GEN Tq = FpX_red(d->T, q), aiq = FpX_red(d->ai, q);
     785                 :            :   (void)M;
     786                 :       1778 :   return FpXQ_mul(FpXQ_mul(V, gel(v,2), Tq, q), aiq, Tq, q);
     787                 :            : }
     788                 :            : 
     789                 :            : static GEN
     790                 :       1778 : _sqrtn_eval(void *E, GEN x, GEN q)
     791                 :            : {
     792                 :       1778 :   struct _ZpXQ_sqrtn *d = (struct _ZpXQ_sqrtn *) E;
     793                 :       1778 :   GEN Tq = FpX_red(d->T, q);
     794                 :       1778 :   GEN f = FpX_sub(FpXQ_pow(x, d->n, Tq, q), d->a, q);
     795                 :       1778 :   return mkvec2(f, x);
     796                 :            : }
     797                 :            : 
     798                 :            : GEN
     799                 :       1141 : ZpXQ_sqrtnlift(GEN a, GEN n, GEN x, GEN T, GEN p, long e)
     800                 :            : {
     801                 :            :   struct _ZpXQ_sqrtn d;
     802                 :       1141 :   d.a = a; d.T = T; d.n = n;
     803                 :       1141 :   d.ai = ZpXQ_inv(ZX_Z_mul(a, n),T,p,(e+1)>>1);
     804                 :       1141 :   return gen_ZpX_Newton(x, p, e, &d, _sqrtn_eval, _sqrtn_invd);
     805                 :            : }
     806                 :            : 
     807                 :            : GEN
     808                 :          0 : ZpXQ_sqrt(GEN a, GEN T, GEN p, long e)
     809                 :            : {
     810                 :          0 :   pari_sp av = avma;
     811                 :          0 :   GEN z = FpXQ_sqrt(FpX_red(a, p), T, p);
     812         [ #  # ]:          0 :   if (!z) return NULL;
     813         [ #  # ]:          0 :   if (e <= 1) return gerepileupto(av, z);
     814                 :          0 :   return gerepileupto(av, ZpXQ_sqrtnlift(a, gen_2, z, T, p, e));
     815                 :            : }
     816                 :            : 
     817                 :            : GEN
     818                 :       1799 : ZpX_ZpXQ_liftroot_ea(GEN P, GEN S, GEN T, GEN p, long n, void *E,
     819                 :            :                      int early(void *E, GEN x, GEN q))
     820                 :            : {
     821                 :       1799 :   pari_sp ltop = avma, av;
     822                 :            :   long N, r;
     823                 :            :   long mask;
     824                 :            :   GEN q2, q, W, Q, Tq2, Tq, Pq;
     825                 :            :   pari_timer ti;
     826                 :       1799 :   T = FpX_get_red(T, powiu(p, n));
     827         [ -  + ]:       1799 :   if (n == 1) return gcopy(S);
     828                 :       1799 :   mask = quadratic_prec_mask(n);
     829                 :       1799 :   av = avma;
     830                 :       1799 :   q2 = p; q = sqri(p); mask >>= 1; N = 2;
     831         [ -  + ]:       1799 :   if (DEBUGLEVEL > 3) timer_start(&ti);
     832                 :       1799 :   Tq = FpXT_red(T,q);
     833                 :       1799 :   Tq2 = FpXT_red(Tq,q2);
     834                 :       1799 :   Pq = FpX_red(P,q);
     835                 :       1799 :   W = FpXQ_inv(FpX_FpXQ_eval(FpX_deriv(P,q2), S, Tq2, q2), Tq2, q2);
     836                 :       1799 :   Q  = ZX_Z_divexact(FpX_FpXQ_eval(Pq, S, Tq, q), q2);
     837                 :       1799 :   r = brent_kung_optpow(degpol(P), 4, 3);
     838                 :            :   for (;;)
     839                 :            :   {
     840                 :            :     GEN H, Sq, Wq, Spow, dP, qq, Pqq, Tqq;
     841                 :       6400 :     H  = FpXQ_mul(W, Q, Tq2, q2);
     842                 :       6400 :     Sq = FpX_sub(S, ZX_Z_mul(H, q2), q);
     843         [ -  + ]:       6400 :     if (DEBUGLEVEL > 3)
     844                 :          0 :       timer_printf(&ti,"ZpX_ZpXQ_liftroot: lift to prec %ld",N);
     845 [ +  + ][ +  + ]:       6400 :     if (mask==1 || (early && early(E, Sq, q)))
                 [ +  + ]
     846                 :       1799 :       return gerepileupto(ltop, Sq);
     847                 :       4601 :     qq = sqri(q); N <<= 1;
     848         [ +  + ]:       4601 :     if (mask&1UL) { qq = diviiexact(qq, p); N--; }
     849                 :       4601 :     mask >>= 1;
     850                 :       4601 :     Pqq  = FpX_red(P, qq);
     851                 :       4601 :     Tqq  = FpXT_red(T, qq);
     852                 :       4601 :     Spow = FpXQ_powers(Sq, r, Tqq, qq);
     853                 :       4601 :     Q  = ZX_Z_divexact(FpX_FpXQV_eval(Pqq, Spow, Tqq, qq), q);
     854                 :       4601 :     dP = FpX_FpXQV_eval(FpX_deriv(Pq, q), FpXV_red(Spow, q), Tq, q);
     855                 :       4601 :     Wq = ZX_Z_divexact(FpX_Fp_sub(FpXQ_mul(W, dP, Tq, q), gen_1, q), q2);
     856                 :       4601 :     Wq = ZX_Z_mul(FpXQ_mul(W, Wq, Tq2, q2), q2);
     857                 :       4601 :     Wq = FpX_sub(W, Wq, q);
     858                 :       4601 :     S = Sq; W = Wq; q2 = q; q = qq; Tq2 = Tq; Tq = Tqq; Pq = Pqq;
     859         [ -  + ]:       4601 :     if (gc_needed(av, 1))
     860                 :            :     {
     861         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpX_ZpXQ_Newton");
     862                 :          0 :       gerepileall(av, 8, &S, &W, &Q, &Tq2, &Tq, &Pq, &q, &q2);
     863                 :            :     }
     864                 :       6400 :   }
     865                 :            : }
     866                 :            : 
     867                 :            : GEN
     868                 :        210 : ZpX_ZpXQ_liftroot(GEN P, GEN S, GEN T, GEN p, long n)
     869                 :            : {
     870                 :        210 :   return ZpX_ZpXQ_liftroot_ea(P, S, T, p, n, NULL, NULL);
     871                 :            : }
     872                 :            : 
     873                 :            : GEN
     874                 :         56 : ZpX_Frobenius(GEN T, GEN p, long e)
     875                 :            : {
     876                 :         56 :   return ZpX_ZpXQ_liftroot(get_FpX_mod(T), FpX_Frobenius(T, p), T, p, e);
     877                 :            : }
     878                 :            : 
     879                 :            : GEN
     880                 :         28 : ZpXQM_prodFrobenius(GEN M, GEN T, GEN p, long e)
     881                 :            : {
     882                 :         28 :   pari_sp av = avma;
     883                 :         28 :   GEN xp = ZpX_Frobenius(T, p, e);
     884                 :         28 :   GEN z = FpXQM_autsum(mkvec2(xp, M), get_FpX_degree(T), T, powiu(p,e));
     885                 :         28 :   return gerepilecopy(av, gel(z,2));
     886                 :            : }

Generated by: LCOV version 1.9