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 - nffactor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23332-367b47754) Lines: 1155 1241 93.1 %
Date: 2018-12-10 05:41:52 Functions: 69 73 94.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2004  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             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*            POLYNOMIAL FACTORIZATION IN A NUMBER FIELD           */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : static GEN nfsqff(GEN nf,GEN pol,long fl,GEN den);
      23             : static int nfsqff_use_Trager(long n, long dpol);
      24             : 
      25             : enum { FACTORS = 0, ROOTS, ROOTS_SPLIT };
      26             : 
      27             : /* for nf_bestlift: reconstruction of algebraic integers known mod P^k,
      28             :  * P maximal ideal above p */
      29             : typedef struct {
      30             :   long k;    /* input known mod P^k */
      31             :   GEN p, pk; /* p^k = denom(prk^-1) [ assume pr unramified ]*/
      32             :   GEN prk;   /* |.|^2 LLL-reduced basis (b_i) of P^k  (NOT T2-reduced) */
      33             :   GEN iprk;  /* den * prk^-1 */
      34             :   GEN GSmin; /* min |b_i^*|^2 */
      35             : 
      36             :   GEN Tp; /* Tpk mod p */
      37             :   GEN Tpk;
      38             :   GEN ZqProj;/* projector to Zp / P^k = Z/p^k[X] / Tpk */
      39             : 
      40             :   GEN tozk;
      41             :   GEN topow;
      42             :   GEN topowden; /* topow x / topowden = basistoalg(x) */
      43             :   GEN dn; /* NULL (we trust nf.zk) or a t_INT > 1 (an alg. integer has
      44             :              denominator dividing dn, when expressed on nf.zk */
      45             : } nflift_t;
      46             : 
      47             : typedef struct
      48             : {
      49             :   nflift_t *L;
      50             :   GEN nf;
      51             :   GEN pol, polbase; /* leading coeff is a t_INT */
      52             :   GEN fact;
      53             :   GEN Br, bound, ZC, BS_2;
      54             : } nfcmbf_t;
      55             : 
      56             : /*******************************************************************/
      57             : /*              RATIONAL RECONSTRUCTION (use ratlift)              */
      58             : /*******************************************************************/
      59             : /* NOT stack clean. a, b stay on the stack */
      60             : static GEN
      61    14175337 : lift_to_frac_tdenom(GEN t, GEN mod, GEN amax, GEN bmax, GEN denom, GEN tdenom)
      62             : {
      63             :   GEN a, b;
      64    14175337 :   if (signe(t) < 0) t = addii(t, mod); /* in case t is a centerlift */
      65    14175337 :   if (tdenom)
      66             :   {
      67    11988154 :     pari_sp av = avma;
      68    11988154 :     a = Fp_center_i(Fp_mul(t, tdenom, mod), mod, shifti(mod,-1));
      69    11988154 :     if (abscmpii(a, amax) < 0)
      70             :     {
      71    11432164 :       GEN d = gcdii(a, tdenom);
      72    11432164 :       a = diviiexact(a, d);
      73    11432164 :       b = diviiexact(tdenom, d);
      74    11432164 :       if (is_pm1(b)) { return gerepileuptoint(av, a); }
      75     1530935 :       return gerepilecopy(av, mkfrac(a, b));
      76             :     }
      77      555990 :     set_avma(av);
      78             :   }
      79     2743173 :   if (!Fp_ratlift(t, mod, amax,bmax, &a,&b)
      80     2616287 :      || (denom && !dvdii(denom,b))
      81     2614512 :      || !is_pm1(gcdii(a,b))) return NULL;
      82     2614416 :   if (is_pm1(b)) { cgiv(b); return a; }
      83     1233459 :   return mkfrac(a, b);
      84             : }
      85             : 
      86             : static GEN
      87      152101 : lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN denom)
      88             : {
      89      152101 :   return lift_to_frac_tdenom(t, mod, amax, bmax, denom, NULL);
      90             : }
      91             : 
      92             : /* Compute rational lifting for all the components of M modulo mod.
      93             :  * Assume all Fp_ratlift preconditions are met; we allow centerlifts but in
      94             :  * that case are no longer stack clean. If one component fails, return NULL.
      95             :  * If denom != NULL, check that the denominators divide denom.
      96             :  *
      97             :  * We suppose gcd(mod, denom) = 1, then a and b are coprime; so we can use
      98             :  * mkfrac rather than gdiv */
      99             : GEN
     100     2035082 : FpC_ratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)
     101             : {
     102     2035082 :   pari_sp ltop = avma;
     103             :   long j, l;
     104     2035082 :   GEN a, d, tdenom = NULL, Q = cgetg_copy(P, &l);
     105     2035082 :   if (l==1) return Q;
     106    15938395 :   for (j = 1; j < l; ++j)
     107             :   {
     108    14023236 :     a = lift_to_frac_tdenom(gel(P,j), mod, amax, bmax, denom, tdenom);
     109    14023236 :     if (!a) return gc_NULL(ltop);
     110    13903313 :     d = Q_denom(a);
     111    13903313 :     tdenom = tdenom ? cmpii(tdenom, d)<0? d: tdenom : d;
     112    13903313 :     gel(Q,j) = a;
     113             :   }
     114     1915159 :   return Q;
     115             : }
     116             : 
     117             : GEN
     118     1003102 : FpM_ratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom)
     119             : {
     120     1003102 :   pari_sp av = avma;
     121     1003102 :   long j, l = lg(M);
     122     1003102 :   GEN N = cgetg_copy(M, &l);
     123     1003102 :   if (l == 1) return N;
     124     2874084 :   for (j = 1; j < l; ++j)
     125             :   {
     126     2035082 :     GEN a = FpC_ratlift(gel(M, j), mod, amax, bmax, denom);
     127     2035082 :     if (!a) return gc_NULL(av);
     128     1915159 :     gel(N,j) = a;
     129             :   }
     130      839002 :   return N;
     131             : }
     132             : 
     133             : GEN
     134       62947 : FpX_ratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)
     135             : {
     136       62947 :   pari_sp ltop = avma;
     137             :   long j, l;
     138       62947 :   GEN a, Q = cgetg_copy(P, &l);
     139       62947 :   Q[1] = P[1];
     140      206214 :   for (j = 2; j < l; ++j)
     141             :   {
     142      152101 :     a = lift_to_frac(gel(P,j), mod, amax,bmax,denom);
     143      152101 :     if (!a) return gc_NULL(ltop);
     144      143267 :     gel(Q,j) = a;
     145             :   }
     146       54113 :   return Q;
     147             : }
     148             : 
     149             : /*******************************************************************/
     150             : /*              GCD in K[X], K NUMBER FIELD                        */
     151             : /*******************************************************************/
     152             : /* P,Q in Z[X,Y], T in Z[Y] irreducible. compute GCD in Q[Y]/(T)[X].
     153             :  *
     154             :  * M. Encarnacion "On a modular Algorithm for computing GCDs of polynomials
     155             :  * over number fields" (ISSAC'94).
     156             :  *
     157             :  * We procede as follows
     158             :  *  1:compute the gcd modulo primes discarding bad primes as they are detected.
     159             :  *  2:reconstruct the result via FpM_ratlift, stoping as soon as we get weird
     160             :  *    denominators.
     161             :  *  3:if FpM_ratlift succeeds, try the full division.
     162             :  * Suppose accuracy is insufficient to get the result right: FpM_ratlift will
     163             :  * rarely succeed, and even if it does the polynomial we get has sensible
     164             :  * coefficients, so the full division will not be too costly.
     165             :  *
     166             :  * If not NULL, den must be a multiple of the denominator of the gcd,
     167             :  * for example the discriminant of T.
     168             :  *
     169             :  * NOTE: if T is not irreducible, nfgcd may loop forever, esp. if gcd | T */
     170             : GEN
     171       11558 : nfgcd_all(GEN P, GEN Q, GEN T, GEN den, GEN *Pnew)
     172             : {
     173       11558 :   pari_sp btop, ltop = avma;
     174       11558 :   GEN lP, lQ, M, dsol, R, bo, sol, mod = NULL, lden = NULL;
     175       11558 :   long vP = varn(P), vT = varn(T), dT = degpol(T), dM = 0, dR;
     176             :   forprime_t S;
     177             : 
     178       11558 :   if (!signe(P)) { if (Pnew) *Pnew = pol_0(vT); return gcopy(Q); }
     179       11558 :   if (!signe(Q)) { if (Pnew) *Pnew = pol_1(vT);   return gcopy(P); }
     180             :   /*Compute denominators*/
     181       11474 :   lP = leading_coeff(P);
     182       11474 :   lQ = leading_coeff(Q);
     183       11474 :   if ( !((typ(lP)==t_INT && is_pm1(lP)) || (typ(lQ)==t_INT && is_pm1(lQ))) )
     184             :   {
     185        3521 :     lden = gcdii(ZX_resultant(lP, T), ZX_resultant(lQ, T));
     186        3521 :     if (den) den = mulii(den, lden);
     187             :   }
     188             : 
     189       11474 :   init_modular_small(&S);
     190       11474 :   btop = avma;
     191             :   for(;;)
     192        1468 :   {
     193       12942 :     ulong p = u_forprime_next(&S);
     194             :     GEN Tp;
     195       12942 :     if (!p) pari_err_OVERFLOW("nfgcd [ran out of primes]");
     196             :     /*Discard primes dividing disc(T) or lc(PQ) */
     197       12942 :     if (lden && !umodiu(lden, p)) continue;
     198       12942 :     Tp = ZX_to_Flx(T,p);
     199       12942 :     if (!Flx_is_squarefree(Tp, p)) continue;
     200             :     /*Discard primes when modular gcd does not exist*/
     201       12942 :     if ((R = FlxqX_safegcd(ZXX_to_FlxX(P,p,vT),
     202             :                            ZXX_to_FlxX(Q,p,vT),
     203           0 :                            Tp, p)) == NULL) continue;
     204       12942 :     dR = degpol(R);
     205       12942 :     if (dR == 0) { set_avma(ltop); if (Pnew) *Pnew = P; return pol_1(vP); }
     206        2469 :     if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
     207             : 
     208        2469 :     R = FlxX_to_Flm(R, dT);
     209             :     /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
     210        2469 :     if (!mod || dR < dM) { M = ZM_init_CRT(R, p); mod = utoipos(p); dM = dR; continue; }
     211        1468 :     (void)ZM_incremental_CRT(&M,R, &mod,p);
     212        1468 :     if (gc_needed(btop, 1))
     213             :     {
     214           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"nfgcd");
     215           0 :       gerepileall(btop, 2, &M, &mod);
     216             :     }
     217             :     /* I suspect it must be better to take amax > bmax*/
     218        1468 :     bo = sqrti(shifti(mod, -1));
     219        1468 :     if ((sol = FpM_ratlift(M, mod, bo, bo, den)) == NULL) continue;
     220        1001 :     sol = RgM_to_RgXX(sol,vP,vT);
     221        1001 :     dsol = Q_primpart(sol);
     222             : 
     223        1001 :     if (!ZXQX_dvd(Q, dsol, T)) continue;
     224        1001 :     if (Pnew)
     225             :     {
     226         175 :       *Pnew = RgXQX_pseudodivrem(P, dsol, T, &R);
     227         175 :       if (signe(R)) continue;
     228             :     }
     229             :     else
     230             :     {
     231         826 :       if (!ZXQX_dvd(P, dsol, T)) continue;
     232             :     }
     233        1001 :     gerepileall(ltop, Pnew? 2: 1, &dsol, Pnew);
     234        1001 :     return dsol; /* both remainders are 0 */
     235             :   }
     236             : }
     237             : GEN
     238        5236 : nfgcd(GEN P, GEN Q, GEN T, GEN den)
     239        5236 : { return nfgcd_all(P,Q,T,den,NULL); }
     240             : 
     241             : int
     242        3094 : nfissquarefree(GEN nf, GEN x)
     243             : {
     244        3094 :   pari_sp av = avma;
     245        3094 :   GEN g, y = RgX_deriv(x);
     246        3094 :   if (RgX_is_rational(x)) g = QX_gcd(x, y);
     247             :   else
     248             :   {
     249        2121 :     GEN T = get_nfpol(nf,&nf);
     250        2121 :     x = Q_primpart( liftpol_shallow(x) );
     251        2121 :     y = Q_primpart( liftpol_shallow(y) );
     252        2121 :     g = nfgcd(x, y, T, nf? nf_get_index(nf): NULL);
     253             :   }
     254        3094 :   return gc_bool(av, degpol(g) == 0);
     255             : }
     256             : 
     257             : /*******************************************************************/
     258             : /*             FACTOR OVER (Z_K/pr)[X] --> FqX_factor              */
     259             : /*******************************************************************/
     260             : GEN
     261           7 : nffactormod(GEN nf, GEN x, GEN pr)
     262             : {
     263           7 :   long j, l, vx = varn(x), vn;
     264           7 :   pari_sp av = avma;
     265             :   GEN F, E, rep, xrd, modpr, T, p;
     266             : 
     267           7 :   nf = checknf(nf);
     268           7 :   vn = nf_get_varn(nf);
     269           7 :   if (typ(x)!=t_POL) pari_err_TYPE("nffactormod",x);
     270           7 :   if (varncmp(vx,vn) >= 0) pari_err_PRIORITY("nffactormod", x, ">=", vn);
     271             : 
     272           7 :   modpr = nf_to_Fq_init(nf, &pr, &T, &p);
     273           7 :   xrd = nfX_to_FqX(x, nf, modpr);
     274           7 :   rep = FqX_factor(xrd,T,p);
     275           7 :   settyp(rep, t_MAT);
     276           7 :   F = gel(rep,1); l = lg(F);
     277           7 :   E = gel(rep,2); settyp(E, t_COL);
     278          14 :   for (j = 1; j < l; j++) {
     279           7 :     gel(F,j) = FqX_to_nfX(gel(F,j), modpr);
     280           7 :     gel(E,j) = stoi(E[j]);
     281             :   }
     282           7 :   return gerepilecopy(av, rep);
     283             : }
     284             : 
     285             : /*******************************************************************/
     286             : /*               MAIN ROUTINES nfroots / nffactor                  */
     287             : /*******************************************************************/
     288             : static GEN
     289        7183 : QXQX_normalize(GEN P, GEN T)
     290             : {
     291        7183 :   GEN P0 = leading_coeff(P);
     292        7183 :   long t = typ(P0);
     293        7183 :   if (t == t_POL)
     294             :   {
     295        1197 :     if (degpol(P0)) return RgXQX_RgXQ_mul(P, QXQ_inv(P0,T), T);
     296         441 :     P0 = gel(P0,2); t = typ(P0);
     297             :   }
     298             :   /* t = t_INT/t_FRAC */
     299        6427 :   if (t == t_INT && is_pm1(P0) && signe(P0) > 0) return P; /* monic */
     300        2660 :   return RgX_Rg_div(P, P0);
     301             : }
     302             : /* assume leading term of P is an integer */
     303             : static GEN
     304        2940 : RgX_int_normalize(GEN P)
     305             : {
     306        2940 :   GEN P0 = leading_coeff(P);
     307             :   /* cater for t_POL */
     308        2940 :   if (typ(P0) == t_POL)
     309             :   {
     310          87 :     P0 = gel(P0,2); /* non-0 constant */
     311          87 :     P = shallowcopy(P);
     312          87 :     gel(P,lg(P)-1) = P0; /* now leading term is a t_INT */
     313             :   }
     314        2940 :   if (typ(P0) != t_INT) pari_err_BUG("RgX_int_normalize");
     315        2940 :   if (is_pm1(P0)) return signe(P0) > 0? P: RgX_neg(P);
     316        1883 :   return RgX_Rg_div(P, P0);
     317             : }
     318             : 
     319             : /* discard change of variable if nf is of the form [nf,c] as return by nfinit
     320             :  * for non-monic polynomials */
     321             : static GEN
     322         952 : proper_nf(GEN nf)
     323         952 : { return (lg(nf) == 3)? gel(nf,1): nf; }
     324             : 
     325             : /* if *pnf = NULL replace if by a "quick" K = nfinit(T), ensuring maximality
     326             :  * by small primes only. Return a multiplicative bound for the denominator of
     327             :  * algebraic integers in Z_K in terms of K.zk */
     328             : static GEN
     329        6119 : fix_nf(GEN *pnf, GEN *pT, GEN *pA)
     330             : {
     331        6119 :   GEN nf, NF, fa, P, Q, q, D, T = *pT;
     332             :   nfmaxord_t S;
     333             :   long i, l;
     334             : 
     335        6119 :   if (*pnf) return gen_1;
     336         952 :   nfmaxord(&S, T, nf_PARTIALFACT);
     337         952 :   NF = nfinit_complete(&S, 0, DEFAULTPREC);
     338         952 :   *pnf = nf = proper_nf(NF);
     339         952 :   if (nf != NF) { /* t_POL defining base field changed (not monic) */
     340          35 :     GEN A = *pA, a = cgetg_copy(A, &l);
     341          35 :     GEN rev = gel(NF,2), pow, dpow;
     342             : 
     343          35 :     *pT = T = nf_get_pol(nf); /* need to update T */
     344          35 :     pow = QXQ_powers(lift_shallow(rev), degpol(T)-1, T);
     345          35 :     pow = Q_remove_denom(pow, &dpow);
     346          35 :     a[1] = A[1];
     347         154 :     for (i=2; i<l; i++) {
     348         119 :       GEN c = gel(A,i);
     349         119 :       if (typ(c) == t_POL) c = QX_ZXQV_eval(c, pow, dpow);
     350         119 :       gel(a,i) = c;
     351             :     }
     352          35 :     *pA = Q_primpart(a); /* need to update A */
     353             :   }
     354             : 
     355         952 :   D = nf_get_disc(nf);
     356         952 :   if (is_pm1(D)) return gen_1;
     357         945 :   fa = absZ_factor_limit(D, 0);
     358         945 :   P = gel(fa,1); q = gel(P, lg(P)-1);
     359         945 :   if (BPSW_psp(q)) return gen_1;
     360             :   /* nf_get_disc(nf) may be incorrect */
     361          14 :   P = nf_get_ramified_primes(nf);
     362          14 :   l = lg(P);
     363          14 :   Q = q; q = gen_1;
     364          70 :   for (i = 1; i < l; i++)
     365             :   {
     366          56 :     GEN p = gel(P,i);
     367          56 :     if (Z_pvalrem(Q, p, &Q) && !BPSW_psp(p)) q = mulii(q, p);
     368             :   }
     369          14 :   return q;
     370             : }
     371             : 
     372             : /* lt(A) is an integer; ensure it is not a constant t_POL. In place */
     373             : static void
     374        6259 : ensure_lt_INT(GEN A)
     375             : {
     376        6259 :   long n = lg(A)-1;
     377        6259 :   GEN lt = gel(A,n);
     378        6259 :   while (typ(lt) != t_INT) gel(A,n) = lt = gel(lt,2);
     379        6259 : }
     380             : 
     381             : /* set B = A/gcd(A,A'), squarefree */
     382             : static GEN
     383        6245 : get_nfsqff_data(GEN *pnf, GEN *pT, GEN *pA, GEN *pB, GEN *ptbad)
     384             : {
     385        6245 :   GEN den, bad, D, B, A = *pA, T = *pT;
     386        6245 :   long n = degpol(T);
     387             : 
     388        6245 :   A = Q_primpart( QXQX_normalize(A, T) );
     389        6245 :   if (nfsqff_use_Trager(n, degpol(A)))
     390             :   {
     391         210 :     *pnf = T;
     392         210 :     bad = den = ZX_disc(T);
     393         210 :     if (is_pm1(leading_coeff(T))) den = indexpartial(T, den);
     394             :   }
     395             :   else
     396             :   {
     397        6035 :     den = fix_nf(pnf, &T, &A);
     398        6035 :     bad = nf_get_index(*pnf);
     399        6035 :     if (den != gen_1) bad = mulii(bad, den);
     400             :   }
     401        6245 :   D = nfgcd_all(A, RgX_deriv(A), T, bad, &B);
     402        6245 :   if (degpol(D)) B = Q_primpart( QXQX_normalize(B, T) );
     403        6245 :   if (ptbad) *ptbad = bad;
     404        6245 :   *pA = A;
     405        6245 :   *pB = B; ensure_lt_INT(B);
     406        6245 :   *pT = T; return den;
     407             : }
     408             : 
     409             : /* return the roots of pol in nf */
     410             : GEN
     411        7064 : nfroots(GEN nf,GEN pol)
     412             : {
     413        7064 :   pari_sp av = avma;
     414             :   GEN z, A, B, T, den;
     415             :   long d, dT;
     416             : 
     417        7064 :   if (!nf) return nfrootsQ(pol);
     418        5048 :   T = get_nfpol(nf, &nf);
     419        5048 :   RgX_check_ZX(T,"nfroots");
     420        5048 :   A = RgX_nffix("nfroots", T,pol,1);
     421        5048 :   d = degpol(A);
     422        5048 :   if (d < 0) pari_err_ROOTS0("nfroots");
     423        5048 :   if (d == 0) return cgetg(1,t_VEC);
     424        5048 :   if (d == 1)
     425             :   {
     426          14 :     A = QXQX_normalize(A,T);
     427          14 :     A = mkpolmod(gneg_i(gel(A,2)), T);
     428          14 :     return gerepilecopy(av, mkvec(A));
     429             :   }
     430        5034 :   dT = degpol(T);
     431        5034 :   if (dT == 1) return gerepileupto(av, nfrootsQ(simplify_shallow(A)));
     432             : 
     433        4936 :   den = get_nfsqff_data(&nf, &T, &A, &B, NULL);
     434        4936 :   if (RgX_is_ZX(B))
     435             :   {
     436        1485 :     GEN v = gel(ZX_factor(B), 1);
     437        1485 :     long i, l = lg(v), p = mael(factoru(dT),1,1); /* smallest prime divisor */
     438        1485 :     z = cgetg(1, t_VEC);
     439        4041 :     for (i = 1; i < l; i++)
     440             :     {
     441        2556 :       GEN b = gel(v,i); /* irreducible / Q */
     442        2556 :       long db = degpol(b);
     443        2556 :       if (db != 1 && degpol(b) < p) continue;
     444        2556 :       z = shallowconcat(z, nfsqff(nf, b, ROOTS, den));
     445             :     }
     446             :   }
     447             :   else
     448        3451 :     z = nfsqff(nf,B, ROOTS, den);
     449        4936 :   z = gerepileupto(av, QXQV_to_mod(z, T));
     450        4936 :   gen_sort_inplace(z, (void*)&cmp_RgX, &cmp_nodata, NULL);
     451        4936 :   return z;
     452             : }
     453             : 
     454             : static GEN
     455      177450 : _norml2(GEN x) { return RgC_fpnorml2(x, DEFAULTPREC); }
     456             : 
     457             : /* return a minimal lift of elt modulo id, as a ZC */
     458             : static GEN
     459       47172 : nf_bestlift(GEN elt, GEN bound, nflift_t *L)
     460             : {
     461             :   GEN u;
     462       47172 :   long i,l = lg(L->prk), t = typ(elt);
     463       47172 :   if (t != t_INT)
     464             :   {
     465       10266 :     if (t == t_POL) elt = ZM_ZX_mul(L->tozk, elt);
     466       10266 :     u = ZM_ZC_mul(L->iprk,elt);
     467       10266 :     for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
     468             :   }
     469             :   else
     470             :   {
     471       36906 :     u = ZC_Z_mul(gel(L->iprk,1), elt);
     472       36906 :     for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
     473       36906 :     elt = scalarcol(elt, l-1);
     474             :   }
     475       47172 :   u = ZC_sub(elt, ZM_ZC_mul(L->prk, u));
     476       47172 :   if (bound && gcmp(_norml2(u), bound) > 0) u = NULL;
     477       47172 :   return u;
     478             : }
     479             : 
     480             : /* Warning: return L->topowden * (best lift). */
     481             : static GEN
     482       33452 : nf_bestlift_to_pol(GEN elt, GEN bound, nflift_t *L)
     483             : {
     484       33452 :   pari_sp av = avma;
     485       33452 :   GEN u,v = nf_bestlift(elt,bound,L);
     486       33452 :   if (!v) return NULL;
     487       27922 :   if (ZV_isscalar(v))
     488             :   {
     489       17850 :     if (L->topowden)
     490       17850 :       u = mulii(L->topowden, gel(v,1));
     491             :     else
     492           0 :       u = icopy(gel(v,1));
     493       17850 :     u = gerepileuptoint(av, u);
     494             :   }
     495             :   else
     496             :   {
     497       10072 :     v = gclone(v); set_avma(av);
     498       10072 :     u = RgV_dotproduct(L->topow, v);
     499       10072 :     gunclone(v);
     500             :   }
     501       27922 :   return u;
     502             : }
     503             : 
     504             : /* return the T->powden * (lift of pol with coefficients of T2-norm <= C)
     505             :  * if it exists. */
     506             : static GEN
     507        7427 : nf_pol_lift(GEN pol, GEN bound, nflift_t *L)
     508             : {
     509        7427 :   long i, l = lg(pol);
     510        7427 :   GEN x = cgetg(l,t_POL);
     511             : 
     512        7427 :   x[1] = pol[1];
     513        7427 :   gel(x,l-1) = mul_content(gel(pol,l-1), L->topowden);
     514       29988 :   for (i=l-2; i>1; i--)
     515             :   {
     516       28091 :     GEN t = nf_bestlift_to_pol(gel(pol,i), bound, L);
     517       28091 :     if (!t) return NULL;
     518       22561 :     gel(x,i) = t;
     519             :   }
     520        1897 :   return x;
     521             : }
     522             : 
     523             : static GEN
     524           0 : zerofact(long v)
     525             : {
     526           0 :   GEN z = cgetg(3, t_MAT);
     527           0 :   gel(z,1) = mkcol(pol_0(v));
     528           0 :   gel(z,2) = mkcol(gen_1); return z;
     529             : }
     530             : 
     531             : /* Return the factorization of A in Q[X]/(T) in rep [pre-allocated with
     532             :  * cgetg(3,t_MAT)], reclaiming all memory between avma and rep.
     533             :  * y is the vector of irreducible factors of B = Q_primpart( A/gcd(A,A') ).
     534             :  * Bad primes divide 'bad' */
     535             : static void
     536        1323 : fact_from_sqff(GEN rep, GEN A, GEN B, GEN y, GEN T, GEN bad)
     537             : {
     538        1323 :   pari_sp av = (pari_sp)rep;
     539        1323 :   long n = lg(y)-1;
     540             :   GEN ex;
     541             : 
     542        1323 :   if (A != B)
     543             :   { /* not squarefree */
     544          77 :     if (n == 1)
     545             :     { /* perfect power, simple ! */
     546          14 :       long e = degpol(A) / degpol(gel(y,1));
     547          14 :       y = gerepileupto(av, QXQXV_to_mod(y, T));
     548          14 :       ex = mkcol(utoipos(e));
     549             :     }
     550             :     else
     551             :     { /* compute valuations mod a prime of degree 1 (avoid coeff explosion) */
     552          63 :       GEN quo, p, r, Bp, lb = leading_coeff(B), E = cgetalloc(t_VECSMALL,n+1);
     553          63 :       pari_sp av1 = avma;
     554             :       ulong pp;
     555             :       long j;
     556             :       forprime_t S;
     557          63 :       u_forprime_init(&S, degpol(T), ULONG_MAX);
     558         168 :       for (; ; set_avma(av1))
     559             :       {
     560         399 :         pp = u_forprime_next(&S);
     561         231 :         if (! umodiu(bad,pp) || !umodiu(lb, pp)) continue;
     562         217 :         p = utoipos(pp);
     563         217 :         r = FpX_oneroot(T, p);
     564         217 :         if (!r) continue;
     565         112 :         Bp = FpXY_evalx(B, r, p);
     566         112 :         if (FpX_is_squarefree(Bp, p)) break;
     567             :       }
     568             : 
     569          63 :       quo = FpXY_evalx(Q_primpart(A), r, p);
     570         140 :       for (j=n; j>=2; j--)
     571             :       {
     572          77 :         GEN junk, fact = Q_remove_denom(gel(y,j), &junk);
     573          77 :         long e = 0;
     574          77 :         fact = FpXY_evalx(fact, r, p);
     575         182 :         for(;; e++)
     576         182 :         {
     577         259 :           GEN q = FpX_divrem(quo,fact,p,ONLY_DIVIDES);
     578         259 :           if (!q) break;
     579         182 :           quo = q;
     580             :         }
     581          77 :         E[j] = e;
     582             :       }
     583          63 :       E[1] = degpol(quo) / degpol(gel(y,1));
     584          63 :       y = gerepileupto(av, QXQXV_to_mod(y, T));
     585          63 :       ex = zc_to_ZC(E); pari_free((void*)E);
     586             :     }
     587             :   }
     588             :   else
     589             :   {
     590        1246 :     y = gerepileupto(av, QXQXV_to_mod(y, T));
     591        1246 :     ex = const_col(n, gen_1);
     592             :   }
     593        1323 :   gel(rep,1) = y; settyp(y, t_COL);
     594        1323 :   gel(rep,2) = ex;
     595        1323 : }
     596             : 
     597             : /* return the factorization of x in nf */
     598             : GEN
     599        1498 : nffactor(GEN nf,GEN pol)
     600             : {
     601        1498 :   GEN bad, A, B, y, T, den, rep = cgetg(3, t_MAT);
     602        1498 :   pari_sp av = avma;
     603             :   long dA;
     604             :   pari_timer ti;
     605             : 
     606        1498 :   if (DEBUGLEVEL>2) { timer_start(&ti); err_printf("\nEntering nffactor:\n"); }
     607        1498 :   T = get_nfpol(nf, &nf);
     608        1498 :   RgX_check_ZX(T,"nffactor");
     609        1498 :   A = RgX_nffix("nffactor",T,pol,1);
     610        1498 :   dA = degpol(A);
     611        1498 :   if (dA <= 0) {
     612           7 :     set_avma((pari_sp)(rep + 3));
     613           7 :     return (dA == 0)? trivial_fact(): zerofact(varn(pol));
     614             :   }
     615        1491 :   if (dA == 1) {
     616             :     GEN c;
     617         112 :     A = Q_primpart( QXQX_normalize(A, T) );
     618         112 :     A = gerepilecopy(av, A); c = gel(A,2);
     619         112 :     if (typ(c) == t_POL && degpol(c) > 0) gel(A,2) = mkpolmod(c, ZX_copy(T));
     620         112 :     gel(rep,1) = mkcol(A);
     621         112 :     gel(rep,2) = mkcol(gen_1); return rep;
     622             :   }
     623        1379 :   if (degpol(T) == 1) return gerepileupto(av, QX_factor(simplify_shallow(A)));
     624             : 
     625        1309 :   den = get_nfsqff_data(&nf, &T, &A, &B, &bad);
     626        1309 :   if (DEBUGLEVEL>2) timer_printf(&ti, "squarefree test");
     627        1309 :   if (RgX_is_ZX(B))
     628             :   {
     629         924 :     GEN v = gel(ZX_factor(B), 1);
     630         924 :     long i, l = lg(v);
     631         924 :     y = cgetg(1, t_VEC);
     632        1897 :     for (i = 1; i < l; i++)
     633             :     {
     634         973 :       GEN b = gel(v,i); /* irreducible / Q */
     635         973 :       y = shallowconcat(y, nfsqff(nf, b, 0, den));
     636             :     }
     637             :   }
     638             :   else
     639         385 :     y = nfsqff(nf,B, 0, den);
     640        1309 :   if (DEBUGLEVEL>3) err_printf("number of factor(s) found: %ld\n", lg(y)-1);
     641             : 
     642        1309 :   fact_from_sqff(rep, A, B, y, T, bad);
     643        1309 :   return sort_factor_pol(rep, cmp_RgX);
     644             : }
     645             : 
     646             : /* assume x scalar or t_COL, G t_MAT */
     647             : static GEN
     648       23926 : arch_for_T2(GEN G, GEN x)
     649             : {
     650       23926 :   return (typ(x) == t_COL)? RgM_RgC_mul(G,x)
     651       23926 :                           : RgC_Rg_mul(gel(G,1),x);
     652             : }
     653             : 
     654             : /* polbase a zkX with t_INT leading coeff; return a bound for T_2(P),
     655             :  * P | polbase in C[X]. NB: Mignotte bound: A | S ==>
     656             :  *  |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S)
     657             :  *
     658             :  * Apply to sigma(S) for all embeddings sigma, then take the L_2 norm over
     659             :  * sigma, then take the sup over i */
     660             : static GEN
     661        1064 : nf_Mignotte_bound(GEN nf, GEN polbase)
     662        1064 : { GEN lS = leading_coeff(polbase); /* t_INT */
     663             :   GEN p1, C, N2, binlS, bin;
     664        1064 :   long prec = nf_get_prec(nf), n = nf_get_degree(nf), r1 = nf_get_r1(nf);
     665        1064 :   long i, j, d = degpol(polbase);
     666             : 
     667        1064 :   binlS = bin = vecbinomial(d-1);
     668        1064 :   if (!isint1(lS)) binlS = ZC_Z_mul(bin,lS);
     669             : 
     670        1064 :   N2 = cgetg(n+1, t_VEC);
     671             :   for (;;)
     672           0 :   {
     673        1064 :     GEN G = nf_get_G(nf), matGS = cgetg(d+2, t_MAT);
     674             : 
     675        1064 :     for (j=0; j<=d; j++) gel(matGS,j+1) = arch_for_T2(G, gel(polbase,j+2));
     676        1064 :     matGS = shallowtrans(matGS);
     677        2688 :     for (j=1; j <= r1; j++) /* N2[j] = || sigma_j(S) ||_2 */
     678             :     {
     679        1624 :       GEN c = sqrtr( _norml2(gel(matGS,j)) );
     680        1624 :       gel(N2,j) = c; if (!signe(c)) goto PRECPB;
     681             :     }
     682        2485 :     for (   ; j <= n; j+=2)
     683             :     {
     684        1421 :       GEN q1 = _norml2(gel(matGS, j));
     685        1421 :       GEN q2 = _norml2(gel(matGS, j+1));
     686        1421 :       GEN c = sqrtr( gmul2n(addrr(q1, q2), -1) );
     687        1421 :       gel(N2,j) = gel(N2,j+1) = c; if (!signe(c)) goto PRECPB;
     688             :     }
     689        1064 :     break; /* done */
     690             : PRECPB:
     691           0 :     prec = precdbl(prec);
     692           0 :     nf = nfnewprec_shallow(nf, prec);
     693           0 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
     694             :   }
     695             : 
     696             :   /* Take sup over 0 <= i <= d of
     697             :    * sum_j | binom(d-1, i-1) ||sigma_j(S)||_2 + binom(d-1,i) lc(S) |^2 */
     698             : 
     699             :   /* i = 0: n lc(S)^2 */
     700        1064 :   C = mului(n, sqri(lS));
     701             :   /* i = d: sum_sigma ||sigma(S)||_2^2 */
     702        1064 :   p1 = gnorml2(N2); if (gcmp(C, p1) < 0) C = p1;
     703       13440 :   for (i = 1; i < d; i++)
     704             :   {
     705       12376 :     GEN B = gel(bin,i), L = gel(binlS,i+1);
     706       12376 :     GEN s = sqrr(addri(mulir(B, gel(N2,1)),  L)); /* j=1 */
     707       12376 :     for (j = 2; j <= n; j++) s = addrr(s, sqrr(addri(mulir(B, gel(N2,j)), L)));
     708       12376 :     if (mpcmp(C, s) < 0) C = s;
     709             :   }
     710        1064 :   return C;
     711             : }
     712             : 
     713             : /* return a bound for T_2(P), P | polbase
     714             :  * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2,
     715             :  * where [P]_2 is Bombieri's 2-norm
     716             :  * Sum over conjugates */
     717             : static GEN
     718        1064 : nf_Beauzamy_bound(GEN nf, GEN polbase)
     719             : {
     720             :   GEN lt, C, s, POL, bin;
     721        1064 :   long d = degpol(polbase), n = nf_get_degree(nf), prec = nf_get_prec(nf);
     722        1064 :   bin = vecbinomial(d);
     723        1064 :   POL = polbase + 2;
     724             :   /* compute [POL]_2 */
     725             :   for (;;)
     726           0 :   {
     727        1064 :     GEN G = nf_get_G(nf);
     728             :     long i;
     729             : 
     730        1064 :     s = real_0(prec);
     731       15568 :     for (i=0; i<=d; i++)
     732             :     {
     733       14504 :       GEN c = gel(POL,i);
     734       14504 :       if (gequal0(c)) continue;
     735        9422 :       c = _norml2(arch_for_T2(G,c));
     736        9422 :       if (!signe(c)) goto PRECPB;
     737             :       /* s += T2(POL[i]) / binomial(d,i) */
     738        9422 :       s = addrr(s, divri(c, gel(bin,i+1)));
     739             :     }
     740        1064 :     break;
     741             : PRECPB:
     742           0 :     prec = precdbl(prec);
     743           0 :     nf = nfnewprec_shallow(nf, prec);
     744           0 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
     745             :   }
     746        1064 :   lt = leading_coeff(polbase);
     747        1064 :   s = mulri(s, muliu(sqri(lt), n));
     748        1064 :   C = powruhalf(utor(3,DEFAULTPREC), 3 + 2*d); /* 3^{3/2 + d} */
     749        1064 :   return divrr(mulrr(C, s), mulur(d, mppi(DEFAULTPREC)));
     750             : }
     751             : 
     752             : static GEN
     753        1064 : nf_factor_bound(GEN nf, GEN polbase)
     754             : {
     755        1064 :   pari_sp av = avma;
     756        1064 :   GEN a = nf_Mignotte_bound(nf, polbase);
     757        1064 :   GEN b = nf_Beauzamy_bound(nf, polbase);
     758        1064 :   if (DEBUGLEVEL>2)
     759             :   {
     760           0 :     err_printf("Mignotte bound: %Ps\n",a);
     761           0 :     err_printf("Beauzamy bound: %Ps\n",b);
     762             :   }
     763        1064 :   return gerepileupto(av, gmin(a, b));
     764             : }
     765             : 
     766             : /* True nf; return Bs: if r a root of sigma_i(P), |r| < Bs[i] */
     767             : static GEN
     768        3417 : nf_root_bounds(GEN nf, GEN P)
     769             : {
     770             :   long lR, i, j, l, prec, r1;
     771             :   GEN Ps, R, V;
     772             : 
     773        3417 :   if (RgX_is_rational(P)) return polrootsbound(P, NULL);
     774        2275 :   r1 = nf_get_r1(nf);
     775        2275 :   P = Q_primpart(P);
     776        2275 :   prec = ZXX_max_lg(P) + 1;
     777        2275 :   l = lg(P);
     778        2275 :   if (nf_get_prec(nf) >= prec)
     779        1952 :     R = nf_get_roots(nf);
     780             :   else
     781         323 :     R = QX_complex_roots(nf_get_pol(nf), prec);
     782        2275 :   lR = lg(R);
     783        2275 :   V = cgetg(lR, t_VEC);
     784        2275 :   Ps = cgetg(l, t_POL); /* sigma (P) */
     785        2275 :   Ps[1] = P[1];
     786        6664 :   for (j=1; j<lg(R); j++)
     787             :   {
     788        4389 :     GEN r = gel(R,j);
     789        4389 :     for (i=2; i<l; i++) gel(Ps,i) = poleval(gel(P,i), r);
     790        4389 :     gel(V,j) = polrootsbound(Ps, NULL);
     791             :   }
     792        2275 :   return mkvec2(vecslice(V,1,r1), vecslice(V,r1+1,lg(V)-1));
     793             : }
     794             : 
     795             : /* return B such that, if x = sum x_i K.zk[i] in O_K, then ||x||_2^2 <= B T_2(x)
     796             :  * den = multiplicative bound for denom(x) [usually NULL, for 1, but when we
     797             :  * use nf_PARTIALFACT K.zk may not generate O_K] */
     798             : static GEN
     799        3846 : L2_bound(GEN nf, GEN den)
     800             : {
     801        3846 :   GEN M, L, prep, T = nf_get_pol(nf), tozk = nf_get_invzk(nf);
     802        3846 :   long prec = ZM_max_lg(tozk) + ZX_max_lg(T) + nbits2prec(degpol(T));
     803        3846 :   (void)initgaloisborne(nf, den? den: gen_1, prec, &L, &prep, NULL);
     804        3846 :   M = vandermondeinverse(L, RgX_gtofp(T,prec), den, prep);
     805        3846 :   return RgM_fpnorml2(RgM_mul(tozk,M), DEFAULTPREC);
     806             : }
     807             : 
     808             : /* sum_i L[i]^p */
     809             : static GEN
     810        8596 : normlp(GEN L, long p)
     811             : {
     812        8596 :   long i, l = lg(L);
     813             :   GEN z;
     814        8596 :   if (l == 1) return gen_0;
     815        4543 :   z = gpowgs(gel(L,1), p);
     816        4543 :   for (i=2; i<l; i++) z = gadd(z, gpowgs(gel(L,i), p));
     817        4543 :   return z;
     818             : }
     819             : /* \sum_i deg(sigma_i) L[i]^p in dimension n (L may be a scalar
     820             :  * or [L1,L2], where Ld corresponds to the archimedean places of degree d) */
     821             : static GEN
     822        6015 : normTp(GEN L, long p, long n)
     823             : {
     824        6015 :   if (typ(L) != t_VEC) return gmulsg(n, gpowgs(L, p));
     825        4298 :   return gadd(normlp(gel(L,1),p), gmul2n(normlp(gel(L,2),p), 1));
     826             : }
     827             : 
     828             : /* S = S0 + tS1, P = P0 + tP1 (Euclidean div. by t integer). For a true
     829             :  * factor (vS, vP), we have:
     830             :  *    | S vS + P vP |^2 < Btra
     831             :  * This implies | S1 vS + P1 vP |^2 < Bhigh, assuming t > sqrt(Btra).
     832             :  * d = dimension of low part (= [nf:Q])
     833             :  * n0 = bound for |vS|^2
     834             :  * */
     835             : static double
     836        1113 : get_Bhigh(long n0, long d)
     837             : {
     838        1113 :   double sqrtd = sqrt((double)d);
     839        1113 :   double z = n0*sqrtd + sqrtd/2 * (d * (n0+1));
     840        1113 :   z = 1. + 0.5 * z; return z * z;
     841             : }
     842             : 
     843             : typedef struct {
     844             :   GEN d;
     845             :   GEN dPinvS;   /* d P^(-1) S   [ integral ] */
     846             :   double **PinvSdbl; /* P^(-1) S as double */
     847             :   GEN S1, P1;   /* S = S0 + S1 q, idem P */
     848             : } trace_data;
     849             : 
     850             : /* S1 * u - P1 * round(P^-1 S u). K non-zero coords in u given by ind */
     851             : static GEN
     852      135471 : get_trace(GEN ind, trace_data *T)
     853             : {
     854      135471 :   long i, j, l, K = lg(ind)-1;
     855             :   GEN z, s, v;
     856             : 
     857      135471 :   s = gel(T->S1, ind[1]);
     858      135471 :   if (K == 1) return s;
     859             : 
     860             :   /* compute s = S1 u */
     861      132797 :   for (j=2; j<=K; j++) s = ZC_add(s, gel(T->S1, ind[j]));
     862             : 
     863             :   /* compute v := - round(P^1 S u) */
     864      132797 :   l = lg(s);
     865      132797 :   v = cgetg(l, t_VECSMALL);
     866     1862749 :   for (i=1; i<l; i++)
     867             :   {
     868     1729952 :     double r, t = 0.;
     869             :     /* quick approximate computation */
     870     1729952 :     for (j=1; j<=K; j++) t += T->PinvSdbl[ ind[j] ][i];
     871     1729952 :     r = floor(t + 0.5);
     872     1729952 :     if (fabs(t + 0.5 - r) < 0.0001)
     873             :     { /* dubious, compute exactly */
     874         196 :       z = gen_0;
     875         196 :       for (j=1; j<=K; j++) z = addii(z, ((GEN**)T->dPinvS)[ ind[j] ][i]);
     876         196 :       v[i] = - itos( diviiround(z, T->d) );
     877             :     }
     878             :     else
     879     1729756 :       v[i] = - (long)r;
     880             :   }
     881      132797 :   return ZC_add(s, ZM_zc_mul(T->P1, v));
     882             : }
     883             : 
     884             : static trace_data *
     885        2128 : init_trace(trace_data *T, GEN S, nflift_t *L, GEN q)
     886             : {
     887        2128 :   long e = gexpo(S), i,j, l,h;
     888             :   GEN qgood, S1, invd;
     889             : 
     890        2128 :   if (e < 0) return NULL; /* S = 0 */
     891             : 
     892        1974 :   qgood = int2n(e - 32); /* single precision check */
     893        1974 :   if (cmpii(qgood, q) > 0) q = qgood;
     894             : 
     895        1974 :   S1 = gdivround(S, q);
     896        1974 :   if (gequal0(S1)) return NULL;
     897             : 
     898         371 :   invd = invr(itor(L->pk, DEFAULTPREC));
     899             : 
     900         371 :   T->dPinvS = ZM_mul(L->iprk, S);
     901         371 :   l = lg(S);
     902         371 :   h = lgcols(T->dPinvS);
     903         371 :   T->PinvSdbl = (double**)cgetg(l, t_MAT);
     904        4109 :   for (j = 1; j < l; j++)
     905             :   {
     906        3738 :     double *t = (double *) stack_malloc_align(h * sizeof(double), sizeof(double));
     907        3738 :     GEN c = gel(T->dPinvS,j);
     908        3738 :     pari_sp av = avma;
     909        3738 :     T->PinvSdbl[j] = t;
     910        3738 :     for (i=1; i < h; i++) t[i] = rtodbl(mulri(invd, gel(c,i)));
     911        3738 :     set_avma(av);
     912             :   }
     913             : 
     914         371 :   T->d  = L->pk;
     915         371 :   T->P1 = gdivround(L->prk, q);
     916         371 :   T->S1 = S1; return T;
     917             : }
     918             : 
     919             : static void
     920       17374 : update_trace(trace_data *T, long k, long i)
     921             : {
     922       17374 :   if (!T) return;
     923        9513 :   gel(T->S1,k)     = gel(T->S1,i);
     924        9513 :   gel(T->dPinvS,k) = gel(T->dPinvS,i);
     925        9513 :   T->PinvSdbl[k]   = T->PinvSdbl[i];
     926             : }
     927             : 
     928             : /* reduce coeffs mod (T,pk), then center mod pk */
     929             : static GEN
     930       16926 : FqX_centermod(GEN z, GEN T, GEN pk, GEN pks2)
     931             : {
     932             :   long i, l;
     933             :   GEN y;
     934       16926 :   if (!T) return centermod_i(z, pk, pks2);
     935       14266 :   y = FpXQX_red(z, T, pk); l = lg(y);
     936      134547 :   for (i = 2; i < l; i++)
     937             :   {
     938      120281 :     GEN c = gel(y,i);
     939      120281 :     if (typ(c) == t_INT)
     940       80248 :       c = Fp_center_i(c, pk, pks2);
     941             :     else
     942       40033 :       c = FpX_center_i(c, pk, pks2);
     943      120281 :     gel(y,i) = c;
     944             :   }
     945       14266 :   return y;
     946             : }
     947             : 
     948             : typedef struct {
     949             :   GEN lt, C, Clt, C2lt, C2ltpol;
     950             : } div_data;
     951             : 
     952             : static void
     953        3888 : init_div_data(div_data *D, GEN pol, nflift_t *L)
     954             : {
     955        3888 :   GEN C2lt, Clt, C = mul_content(L->topowden, L->dn);
     956        3888 :   GEN lc = leading_coeff(pol), lt = is_pm1(lc)? NULL: absi_shallow(lc);
     957        3888 :   if (C)
     958             :   {
     959        3888 :     GEN C2 = sqri(C);
     960        3888 :     if (lt) {
     961        1036 :       C2lt = mulii(C2, lt);
     962        1036 :       Clt = mulii(C,lt);
     963             :     } else {
     964        2852 :       C2lt = C2;
     965        2852 :       Clt = C;
     966             :     }
     967             :   }
     968             :   else
     969           0 :     C2lt = Clt = lt;
     970        3888 :   D->lt = lt;
     971        3888 :   D->C = C;
     972        3888 :   D->Clt = Clt;
     973        3888 :   D->C2lt = C2lt;
     974        3888 :   D->C2ltpol = C2lt? RgX_Rg_mul(pol, C2lt): pol;
     975        3888 : }
     976             : static void
     977        1799 : update_target(div_data *D, GEN pol)
     978        1799 : { D->C2ltpol = D->Clt? RgX_Rg_mul(pol, D->Clt): pol; }
     979             : 
     980             : /* nb = number of modular factors; return a "good" K such that naive
     981             :  * recombination of up to maxK modular factors is not too costly */
     982             : long
     983       14155 : cmbf_maxK(long nb)
     984             : {
     985       14155 :   if (nb >  10) return 3;
     986       13588 :   return nb-1;
     987             : }
     988             : /* Naive recombination of modular factors: combine up to maxK modular
     989             :  * factors, degree <= klim
     990             :  *
     991             :  * target = polynomial we want to factor
     992             :  * famod = array of modular factors.  Product should be congruent to
     993             :  * target/lc(target) modulo p^a
     994             :  * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
     995             : /* set *done = 1 if factorisation is known to be complete */
     996             : static GEN
     997        1064 : nfcmbf(nfcmbf_t *T, long klim, long *pmaxK, int *done)
     998             : {
     999        1064 :   GEN nf = T->nf, famod = T->fact, bound = T->bound;
    1000        1064 :   GEN ltdn, nfpol = nf_get_pol(nf);
    1001        1064 :   long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol);
    1002        1064 :   pari_sp av0 = avma;
    1003        1064 :   GEN Tpk = T->L->Tpk, pk = T->L->pk, pks2 = shifti(pk,-1);
    1004        1064 :   GEN ind      = cgetg(lfamod+1, t_VECSMALL);
    1005        1064 :   GEN deg      = cgetg(lfamod+1, t_VECSMALL);
    1006        1064 :   GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
    1007        1064 :   GEN fa       = cgetg(lfamod+1, t_VEC);
    1008        1064 :   const double Bhigh = get_Bhigh(lfamod, dnf);
    1009             :   trace_data _T1, _T2, *T1, *T2;
    1010             :   div_data D;
    1011             :   pari_timer ti;
    1012             : 
    1013        1064 :   timer_start(&ti);
    1014             : 
    1015        1064 :   *pmaxK = cmbf_maxK(lfamod);
    1016        1064 :   init_div_data(&D, T->pol, T->L);
    1017        1064 :   ltdn = mul_content(D.lt, T->L->dn);
    1018             :   {
    1019        1064 :     GEN q = ceil_safe(sqrtr(T->BS_2));
    1020        1064 :     GEN t1,t2, lt2dn = mul_content(ltdn, D.lt);
    1021        1064 :     GEN trace1   = cgetg(lfamod+1, t_MAT);
    1022        1064 :     GEN trace2   = cgetg(lfamod+1, t_MAT);
    1023        5635 :     for (i=1; i <= lfamod; i++)
    1024             :     {
    1025        4571 :       pari_sp av = avma;
    1026        4571 :       GEN P = gel(famod,i);
    1027        4571 :       long d = degpol(P);
    1028             : 
    1029        4571 :       deg[i] = d; P += 2;
    1030        4571 :       t1 = gel(P,d-1);/* = - S_1 */
    1031        4571 :       t2 = Fq_sqr(t1, Tpk, pk);
    1032        4571 :       if (d > 1) t2 = Fq_sub(t2, gmul2n(gel(P,d-2), 1), Tpk, pk);
    1033             :       /* t2 = S_2 Newton sum */
    1034        4571 :       if (ltdn)
    1035             :       {
    1036         294 :         t1 = Fq_Fp_mul(t1, ltdn, Tpk, pk);
    1037         294 :         t2 = Fq_Fp_mul(t2, lt2dn, Tpk, pk);
    1038             :       }
    1039        4571 :       gel(trace1,i) = gclone( nf_bestlift(t1, NULL, T->L) );
    1040        4571 :       gel(trace2,i) = gclone( nf_bestlift(t2, NULL, T->L) ); set_avma(av);
    1041             :     }
    1042        1064 :     T1 = init_trace(&_T1, trace1, T->L, q);
    1043        1064 :     T2 = init_trace(&_T2, trace2, T->L, q);
    1044        5635 :     for (i=1; i <= lfamod; i++) {
    1045        4571 :       gunclone(gel(trace1,i));
    1046        4571 :       gunclone(gel(trace2,i));
    1047             :     }
    1048             :   }
    1049        1064 :   degsofar[0] = 0; /* sentinel */
    1050             : 
    1051             :   /* ind runs through strictly increasing sequences of length K,
    1052             :    * 1 <= ind[i] <= lfamod */
    1053             : nextK:
    1054        1547 :   if (K > *pmaxK || 2*K > lfamod) goto END;
    1055        1316 :   if (DEBUGLEVEL > 3)
    1056           0 :     err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
    1057        1316 :   setlg(ind, K+1); ind[1] = 1;
    1058        1316 :   i = 1; curdeg = deg[ind[1]];
    1059             :   for(;;)
    1060             :   { /* try all combinations of K factors */
    1061      297332 :     for (j = i; j < K; j++)
    1062             :     {
    1063       16870 :       degsofar[j] = curdeg;
    1064       16870 :       ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
    1065             :     }
    1066      140889 :     if (curdeg <= klim) /* trial divide */
    1067             :     {
    1068             :       GEN t, y, q;
    1069             :       pari_sp av;
    1070             : 
    1071      140889 :       av = avma;
    1072      140889 :       if (T1)
    1073             :       { /* d-1 test */
    1074       52353 :         t = get_trace(ind, T1);
    1075       52353 :         if (rtodbl(_norml2(t)) > Bhigh)
    1076             :         {
    1077       51268 :           if (DEBUGLEVEL>6) err_printf(".");
    1078       51268 :           set_avma(av); goto NEXT;
    1079             :         }
    1080             :       }
    1081       89621 :       if (T2)
    1082             :       { /* d-2 test */
    1083       83118 :         t = get_trace(ind, T2);
    1084       83118 :         if (rtodbl(_norml2(t)) > Bhigh)
    1085             :         {
    1086       82285 :           if (DEBUGLEVEL>3) err_printf("|");
    1087       82285 :           set_avma(av); goto NEXT;
    1088             :         }
    1089             :       }
    1090        7336 :       set_avma(av);
    1091        7336 :       y = ltdn; /* full computation */
    1092       24262 :       for (i=1; i<=K; i++)
    1093             :       {
    1094       16926 :         GEN q = gel(famod, ind[i]);
    1095       16926 :         if (y) q = gmul(y, q);
    1096       16926 :         y = FqX_centermod(q, Tpk, pk, pks2);
    1097             :       }
    1098        7336 :       y = nf_pol_lift(y, bound, T->L);
    1099        7336 :       if (!y)
    1100             :       {
    1101        5530 :         if (DEBUGLEVEL>3) err_printf("@");
    1102        5530 :         set_avma(av); goto NEXT;
    1103             :       }
    1104             :       /* y = topowden*dn*lt*\prod_{i in ind} famod[i] is apparently in O_K[X],
    1105             :        * in fact in (Z[Y]/nf.pol)[X] due to multiplication by C = topowden*dn.
    1106             :        * Try out this candidate factor */
    1107        1806 :       q = RgXQX_divrem(D.C2ltpol, y, nfpol, ONLY_DIVIDES);
    1108        1806 :       if (!q)
    1109             :       {
    1110          70 :         if (DEBUGLEVEL>3) err_printf("*");
    1111          70 :         set_avma(av); goto NEXT;
    1112             :       }
    1113             :       /* Original T->pol in O_K[X] with leading coeff lt in Z,
    1114             :        * y = C*lt \prod famod[i] is in O_K[X] with leading coeff in Z
    1115             :        * q = C^2*lt*pol / y = C * (lt*pol) / (lt*\prod famod[i]) is a
    1116             :        * K-rational factor, in fact in Z[Y]/nf.pol)[X] as above, with
    1117             :        * leading term C*lt. */
    1118        1736 :       update_target(&D, q);
    1119        1736 :       gel(fa,cnt++) = D.C2lt? RgX_int_normalize(y): y; /* make monic */
    1120       12474 :       for (i=j=k=1; i <= lfamod; i++)
    1121             :       { /* remove used factors */
    1122       10738 :         if (j <= K && i == ind[j]) j++;
    1123             :         else
    1124             :         {
    1125        8687 :           gel(famod,k) = gel(famod,i);
    1126        8687 :           update_trace(T1, k, i);
    1127        8687 :           update_trace(T2, k, i);
    1128        8687 :           deg[k] = deg[i]; k++;
    1129             :         }
    1130             :       }
    1131        1736 :       lfamod -= K;
    1132        1736 :       *pmaxK = cmbf_maxK(lfamod);
    1133        1736 :       if (lfamod < 2*K) goto END;
    1134         903 :       i = 1; curdeg = deg[ind[1]];
    1135         903 :       if (DEBUGLEVEL > 2)
    1136             :       {
    1137           0 :         err_printf("\n"); timer_printf(&ti, "to find factor %Ps",gel(fa,cnt-1));
    1138           0 :         err_printf("remaining modular factor(s): %ld\n", lfamod);
    1139             :       }
    1140         903 :       continue;
    1141             :     }
    1142             : 
    1143             : NEXT:
    1144      139153 :     for (i = K+1;;)
    1145             :     {
    1146      173229 :       if (--i == 0) { K++; goto nextK; }
    1147      155708 :       if (++ind[i] <= lfamod - K + i)
    1148             :       {
    1149      138670 :         curdeg = degsofar[i-1] + deg[ind[i]];
    1150      138670 :         if (curdeg <= klim) break;
    1151             :       }
    1152             :     }
    1153             :   }
    1154             : END:
    1155        1064 :   *done = 1;
    1156        1064 :   if (degpol(D.C2ltpol) > 0)
    1157             :   { /* leftover factor */
    1158        1064 :     GEN q = D.C2ltpol;
    1159        1064 :     if (D.C2lt) q = RgX_int_normalize(q);
    1160        1064 :     if (lfamod >= 2*K)
    1161             :     { /* restore leading coefficient [#930] */
    1162          49 :       if (D.lt) q = RgX_Rg_mul(q, D.lt);
    1163          49 :       *done = 0; /* ... may still be reducible */
    1164             :     }
    1165        1064 :     setlg(famod, lfamod+1);
    1166        1064 :     gel(fa,cnt++) = q;
    1167             :   }
    1168        1064 :   if (DEBUGLEVEL>6) err_printf("\n");
    1169        1064 :   setlg(fa, cnt);
    1170        1064 :   return gerepilecopy(av0, fa);
    1171             : }
    1172             : 
    1173             : static GEN
    1174          42 : nf_chk_factors(nfcmbf_t *T, GEN P, GEN M_L, GEN famod, GEN pk)
    1175             : {
    1176          42 :   GEN nf = T->nf, bound = T->bound;
    1177          42 :   GEN nfT = nf_get_pol(nf);
    1178             :   long i, r;
    1179          42 :   GEN pol = P, list, piv, y;
    1180          42 :   GEN Tpk = T->L->Tpk;
    1181             :   div_data D;
    1182             : 
    1183          42 :   piv = ZM_hnf_knapsack(M_L);
    1184          42 :   if (!piv) return NULL;
    1185          28 :   if (DEBUGLEVEL>3) err_printf("ZM_hnf_knapsack output:\n%Ps\n",piv);
    1186             : 
    1187          28 :   r  = lg(piv)-1;
    1188          28 :   list = cgetg(r+1, t_VEC);
    1189          28 :   init_div_data(&D, pol, T->L);
    1190          28 :   for (i = 1;;)
    1191          63 :   {
    1192          91 :     pari_sp av = avma;
    1193          91 :     if (DEBUGLEVEL) err_printf("nf_LLL_cmbf: checking factor %ld\n", i);
    1194          91 :     y = chk_factors_get(D.lt, famod, gel(piv,i), Tpk, pk);
    1195             : 
    1196          91 :     if (! (y = nf_pol_lift(y, bound, T->L)) ) return NULL;
    1197          91 :     y = gerepilecopy(av, y);
    1198             :     /* y is the candidate factor */
    1199          91 :     pol = RgXQX_divrem(D.C2ltpol, y, nfT, ONLY_DIVIDES);
    1200          91 :     if (!pol) return NULL;
    1201             : 
    1202          91 :     if (D.C2lt) y = RgX_int_normalize(y);
    1203          91 :     gel(list,i) = y;
    1204          91 :     if (++i >= r) break;
    1205             : 
    1206          63 :     update_target(&D, pol);
    1207             :   }
    1208          28 :   gel(list,i) = RgX_int_normalize(pol); return list;
    1209             : }
    1210             : 
    1211             : static GEN
    1212       26221 : nf_to_Zq(GEN x, GEN T, GEN pk, GEN pks2, GEN proj)
    1213             : {
    1214             :   GEN y;
    1215       26221 :   if (typ(x) != t_COL) return centermodii(x, pk, pks2);
    1216        6258 :   if (!T)
    1217             :   {
    1218        6062 :     y = ZV_dotproduct(proj, x);
    1219        6062 :     return centermodii(y, pk, pks2);
    1220             :   }
    1221         196 :   y = ZM_ZC_mul(proj, x);
    1222         196 :   y = RgV_to_RgX(y, varn(T));
    1223         196 :   return FpX_center_i(FpX_rem(y, T, pk), pk, pks2);
    1224             : }
    1225             : 
    1226             : /* Assume P in nfX form, lc(P) != 0 mod p. Reduce P to Zp[X]/(T) mod p^a */
    1227             : static GEN
    1228        3424 : ZqX(GEN P, GEN pk, GEN T, GEN proj)
    1229             : {
    1230        3424 :   long i, l = lg(P);
    1231        3424 :   GEN z, pks2 = shifti(pk,-1);
    1232             : 
    1233        3424 :   z = cgetg(l,t_POL); z[1] = P[1];
    1234        3424 :   for (i=2; i<l; i++) gel(z,i) = nf_to_Zq(gel(P,i),T,pk,pks2,proj);
    1235        3424 :   return normalizepol_lg(z, l);
    1236             : }
    1237             : 
    1238             : static GEN
    1239        3424 : ZqX_normalize(GEN P, GEN lt, nflift_t *L)
    1240             : {
    1241        3424 :   GEN R = lt? RgX_Rg_mul(P, Fp_inv(lt, L->pk)): P;
    1242        3424 :   return ZqX(R, L->pk, L->Tpk, L->ZqProj);
    1243             : }
    1244             : 
    1245             : /* k allowing to reconstruct x, |x|^2 < C, from x mod pr^k */
    1246             : /* return log [  2sqrt(C/d) * ( (3/2)sqrt(gamma) )^(d-1) ] ^d / log N(pr)
    1247             :  * cf. Belabas relative van Hoeij algorithm, lemma 3.12 */
    1248             : static double
    1249        3431 : bestlift_bound(GEN C, long d, double alpha, GEN p, long f)
    1250             : {
    1251        3431 :   const double g = 1 / (alpha - 0.25); /* = 2 if alpha = 3/4 */
    1252        3431 :   GEN C4 = shiftr(gtofp(C,DEFAULTPREC), 2);
    1253        3431 :   double t, logp = log(gtodouble(p));
    1254        3431 :   if (f == d)
    1255             :   { /* p inert, no LLL fudge factor: p^(2k) / 4 > C */
    1256          35 :     t = 0.5 * rtodbl(mplog(C4));
    1257          35 :     return ceil(t / logp);
    1258             :   }
    1259             :   /* (1/2)log (4C/d) + (d-1)(log 3/2 sqrt(gamma)) */
    1260        3396 :   t = 0.5 * rtodbl(mplog(divru(C4,d))) + (d-1) * log(1.5 * sqrt(g));
    1261        3396 :   return ceil((t * d) / (logp * f));
    1262             : }
    1263             : 
    1264             : static GEN
    1265        3895 : get_R(GEN M)
    1266             : {
    1267             :   GEN R;
    1268        3895 :   long i, l, prec = nbits2prec( gexpo(M) + 64 );
    1269             : 
    1270             :   for(;;)
    1271             :   {
    1272        3895 :     R = gaussred_from_QR(M, prec);
    1273        3895 :     if (R) break;
    1274           0 :     prec = precdbl(prec);
    1275             :   }
    1276        3895 :   l = lg(R);
    1277        3895 :   for (i=1; i<l; i++) gcoeff(R,i,i) = gen_1;
    1278        3895 :   return R;
    1279             : }
    1280             : 
    1281             : static void
    1282        3853 : init_proj(nflift_t *L, GEN prkHNF, GEN nfT)
    1283             : {
    1284        3853 :   if (degpol(L->Tp)>1)
    1285             :   {
    1286         140 :     GEN coTp = FpX_div(FpX_red(nfT, L->p), L->Tp,  L->p); /* Tp's cofactor */
    1287             :     GEN z, proj;
    1288         140 :     z = ZpX_liftfact(nfT, mkvec2(L->Tp, coTp), L->pk, L->p, L->k);
    1289         140 :     L->Tpk = gel(z,1);
    1290         140 :     proj = QXQV_to_FpM(L->topow, L->Tpk, L->pk);
    1291         140 :     if (L->topowden)
    1292         140 :       proj = FpM_red(ZM_Z_mul(proj, Fp_inv(L->topowden, L->pk)), L->pk);
    1293         140 :     L->ZqProj = proj;
    1294             :   }
    1295             :   else
    1296             :   {
    1297        3713 :     L->Tpk = NULL;
    1298        3713 :     L->ZqProj = dim1proj(prkHNF);
    1299             :   }
    1300        3853 : }
    1301             : 
    1302             : /* Square of the radius of largest ball inscript in PRK's fundamental domain,
    1303             :  *   whose orthogonalized vector's norms are the Bi
    1304             :  * Rmax ^2 =  min 1/4T_i where T_i = sum_j ( s_ij^2 / B_j)
    1305             :  * For p inert, S = Id, T_i = 1 / p^{2k} and Rmax = p^k / 2 */
    1306             : static GEN
    1307        3895 : max_radius(GEN PRK, GEN B)
    1308             : {
    1309        3895 :   GEN S, smax = gen_0;
    1310        3895 :   pari_sp av = avma;
    1311        3895 :   long i, j, d = lg(PRK)-1;
    1312             : 
    1313        3895 :   S = RgM_inv( get_R(PRK) ); if (!S) pari_err_PREC("max_radius");
    1314       18074 :   for (i=1; i<=d; i++)
    1315             :   {
    1316       14179 :     GEN s = gen_0;
    1317      157786 :     for (j=1; j<=d; j++)
    1318      143607 :       s = mpadd(s, mpdiv( mpsqr(gcoeff(S,i,j)), gel(B,j)));
    1319       14179 :     if (mpcmp(s, smax) > 0) smax = s;
    1320             :   }
    1321        3895 :   return gerepileupto(av, ginv(gmul2n(smax, 2)));
    1322             : }
    1323             : 
    1324             : static void
    1325        3853 : bestlift_init(long a, GEN nf, GEN C, nflift_t *L)
    1326             : {
    1327        3853 :   const double alpha = 0.99; /* LLL parameter */
    1328        3853 :   const long d = nf_get_degree(nf);
    1329        3853 :   pari_sp av = avma, av2;
    1330        3853 :   GEN prk, PRK, iPRK, GSmin, T = L->Tp, p = L->p;
    1331        3853 :   long f = degpol(T);
    1332             :   pari_timer ti;
    1333             : 
    1334        3853 :   if (f == d)
    1335             :   { /* inert p, much simpler */
    1336          35 :     long a0 = bestlift_bound(C, d, alpha, p, f);
    1337             :     GEN q;
    1338          35 :     if (a < a0) a = a0; /* guarantees GSmin >= C */
    1339          35 :     if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
    1340          35 :     q = powiu(p,a);
    1341          35 :     PRK = prk = scalarmat_shallow(q, d);
    1342          35 :     GSmin = shiftr(itor(q, DEFAULTPREC), -1);
    1343          35 :     iPRK = matid(d); goto END;
    1344             :   }
    1345        3818 :   timer_start(&ti);
    1346        3818 :   if (!a) a = (long)bestlift_bound(C, d, alpha, p, f);
    1347          77 :   for (;; set_avma(av), a += (a==1)? 1: (a>>1)) /* roughly a *= 1.5 */
    1348          77 :   {
    1349        3895 :     GEN B, q = powiu(p,a), Tq = FpXQ_powu(T, a, FpX_red(nf_get_pol(nf), q), q);
    1350        3895 :     if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
    1351        3895 :     prk = idealhnf_two(nf, mkvec2(q, Tq));
    1352        3895 :     av2 = avma;
    1353        3895 :     PRK = ZM_lll_norms(prk, alpha, LLL_INPLACE, &B);
    1354        3895 :     GSmin = max_radius(PRK, B);
    1355        3895 :     if (gcmp(GSmin, C) >= 0) break;
    1356             :   }
    1357        3818 :   gerepileall(av2, 2, &PRK, &GSmin);
    1358        3818 :   iPRK = ZM_inv(PRK, NULL);
    1359        3818 :   if (DEBUGLEVEL>2)
    1360           0 :     err_printf("for this exponent, GSmin = %Ps\nTime reduction: %ld\n",
    1361             :                GSmin, timer_delay(&ti));
    1362             : END:
    1363        3853 :   L->k = a;
    1364        3853 :   L->pk = gcoeff(prk,1,1);
    1365        3853 :   L->prk = PRK;
    1366        3853 :   L->iprk = iPRK;
    1367        3853 :   L->GSmin= GSmin;
    1368        3853 :   init_proj(L, prk, nf_get_pol(nf));
    1369        3853 : }
    1370             : 
    1371             : /* Let X = Tra * M_L, Y = bestlift(X) return V s.t Y = X - PRK V
    1372             :  * and set *eT2 = gexpo(Y)  [cf nf_bestlift, but memory efficient] */
    1373             : static GEN
    1374         266 : get_V(GEN Tra, GEN M_L, GEN PRK, GEN PRKinv, GEN pk, long *eT2)
    1375             : {
    1376         266 :   long i, e = 0, l = lg(M_L);
    1377         266 :   GEN V = cgetg(l, t_MAT);
    1378         266 :   *eT2 = 0;
    1379        3731 :   for (i = 1; i < l; i++)
    1380             :   { /* cf nf_bestlift(Tra * c) */
    1381        3465 :     pari_sp av = avma, av2;
    1382        3465 :     GEN v, T2 = ZM_ZC_mul(Tra, gel(M_L,i));
    1383             : 
    1384        3465 :     v = gdivround(ZM_ZC_mul(PRKinv, T2), pk); /* small */
    1385        3465 :     av2 = avma;
    1386        3465 :     T2 = ZC_sub(T2, ZM_ZC_mul(PRK, v));
    1387        3465 :     e = gexpo(T2); if (e > *eT2) *eT2 = e;
    1388        3465 :     set_avma(av2);
    1389        3465 :     gel(V,i) = gerepileupto(av, v); /* small */
    1390             :   }
    1391         266 :   return V;
    1392             : }
    1393             : 
    1394             : static GEN
    1395          49 : nf_LLL_cmbf(nfcmbf_t *T, long rec)
    1396             : {
    1397          49 :   const double BitPerFactor = 0.4; /* nb bits / modular factor */
    1398          49 :   nflift_t *L = T->L;
    1399          49 :   GEN famod = T->fact, ZC = T->ZC, Br = T->Br, P = T->pol, dn = T->L->dn;
    1400          49 :   long dnf = nf_get_degree(T->nf), dP = degpol(P);
    1401             :   long i, C, tmax, n0;
    1402             :   GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO, Btra;
    1403             :   double Bhigh;
    1404             :   pari_sp av, av2;
    1405          49 :   long ti_LLL = 0, ti_CF = 0;
    1406             :   pari_timer ti2, TI;
    1407             : 
    1408          49 :   lP = absi_shallow(leading_coeff(P));
    1409          49 :   if (is_pm1(lP)) lP = NULL;
    1410             : 
    1411          49 :   n0 = lg(famod) - 1;
    1412             :  /* Lattice: (S PRK), small vector (vS vP). To find k bound for the image,
    1413             :   * write S = S1 q + S0, P = P1 q + P0
    1414             :   * |S1 vS + P1 vP|^2 <= Bhigh for all (vS,vP) assoc. to true factors */
    1415          49 :   Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2, dnf)));
    1416          49 :   Bhigh = get_Bhigh(n0, dnf);
    1417          49 :   C = (long)ceil(sqrt(Bhigh/n0)) + 1; /* C^2 n0 ~ Bhigh */
    1418          49 :   Bnorm = dbltor( n0 * C * C + Bhigh );
    1419          49 :   ZERO = zeromat(n0, dnf);
    1420             : 
    1421          49 :   av = avma;
    1422          49 :   TT = const_vec(n0, NULL);
    1423          49 :   Tra  = cgetg(n0+1, t_MAT);
    1424          49 :   CM_L = scalarmat_s(C, n0);
    1425             :   /* tmax = current number of traces used (and computed so far) */
    1426         196 :   for(tmax = 0;; tmax++)
    1427         147 :   {
    1428         196 :     long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1;
    1429             :     GEN M_L, q, CM_Lp, oldCM_L, S1, P1, VV;
    1430         196 :     int first = 1;
    1431             : 
    1432             :     /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */
    1433         196 :     Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2*tnew, dnf)));
    1434         196 :     bmin = logint(ceil_safe(sqrtr(Btra)), gen_2) + 1;
    1435         196 :     if (DEBUGLEVEL>2)
    1436           0 :       err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
    1437             :                  r, tmax, bmin);
    1438             : 
    1439             :     /* compute Newton sums (possibly relifting first) */
    1440         196 :     if (gcmp(L->GSmin, Btra) < 0)
    1441             :     {
    1442             :       GEN polred;
    1443             : 
    1444           7 :       bestlift_init((L->k)<<1, T->nf, Btra, L);
    1445           7 :       polred = ZqX_normalize(T->polbase, lP, L);
    1446           7 :       famod = ZqX_liftfact(polred, famod, L->Tpk, L->pk, L->p, L->k);
    1447           7 :       for (i=1; i<=n0; i++) gel(TT,i) = NULL;
    1448             :     }
    1449        4774 :     for (i=1; i<=n0; i++)
    1450             :     {
    1451        4578 :       GEN h, lPpow = lP? powiu(lP, tnew): NULL;
    1452        4578 :       GEN z = polsym_gen(gel(famod,i), gel(TT,i), tnew, L->Tpk, L->pk);
    1453        4578 :       gel(TT,i) = z;
    1454        4578 :       h = gel(z,tnew+1);
    1455             :       /* make Newton sums integral */
    1456        4578 :       lPpow = mul_content(lPpow, dn);
    1457        4578 :       if (lPpow)
    1458         126 :         h = (typ(h) == t_INT)? Fp_mul(h, lPpow, L->pk): FpX_Fp_mul(h, lPpow, L->pk);
    1459        4578 :       gel(Tra,i) = nf_bestlift(h, NULL, L); /* S_tnew(famod) */
    1460             :     }
    1461             : 
    1462             :     /* compute truncation parameter */
    1463         196 :     if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
    1464         196 :     oldCM_L = CM_L;
    1465         196 :     av2 = avma;
    1466         196 :     b = delta = 0; /* -Wall */
    1467             : AGAIN:
    1468         266 :     M_L = Q_div_to_int(CM_L, utoipos(C));
    1469         266 :     VV = get_V(Tra, M_L, L->prk, L->iprk, L->pk, &a);
    1470         266 :     if (first)
    1471             :     { /* initialize lattice, using few p-adic digits for traces */
    1472         196 :       bgood = (long)(a - maxss(32, (long)(BitPerFactor * r)));
    1473         196 :       b = maxss(bmin, bgood);
    1474         196 :       delta = a - b;
    1475             :     }
    1476             :     else
    1477             :     { /* add more p-adic digits and continue reduction */
    1478          70 :       if (a < b) b = a;
    1479          70 :       b = maxss(b-delta, bmin);
    1480          70 :       if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
    1481             :     }
    1482             : 
    1483             :     /* restart with truncated entries */
    1484         266 :     q = int2n(b);
    1485         266 :     P1 = gdivround(L->prk, q);
    1486         266 :     S1 = gdivround(Tra, q);
    1487         266 :     T2 = ZM_sub(ZM_mul(S1, M_L), ZM_mul(P1, VV));
    1488         266 :     m = vconcat( CM_L, T2 );
    1489         266 :     if (first)
    1490             :     {
    1491         196 :       first = 0;
    1492         196 :       m = shallowconcat( m, vconcat(ZERO, P1) );
    1493             :       /*     [ C M_L   0  ]
    1494             :        * m = [            ]   square matrix
    1495             :        *     [  T2'   PRK ]   T2' = Tra * M_L  truncated
    1496             :        */
    1497             :     }
    1498         266 :     CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
    1499         266 :     if (DEBUGLEVEL>2)
    1500           0 :       err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
    1501           0 :                  a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
    1502         315 :     if (!CM_L) { list = mkcol(RgX_int_normalize(P)); break; }
    1503         245 :     if (b > bmin)
    1504             :     {
    1505          70 :       CM_L = gerepilecopy(av2, CM_L);
    1506          70 :       goto AGAIN;
    1507             :     }
    1508         175 :     if (DEBUGLEVEL>2) timer_printf(&ti2, "for this trace");
    1509             : 
    1510         175 :     i = lg(CM_L) - 1;
    1511         175 :     if (i == r && ZM_equal(CM_L, oldCM_L))
    1512             :     {
    1513          91 :       CM_L = oldCM_L;
    1514          91 :       set_avma(av2); continue;
    1515             :     }
    1516             : 
    1517          84 :     CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
    1518          84 :     if (lg(CM_Lp) != lg(CM_L))
    1519             :     {
    1520           0 :       if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
    1521           0 :       CM_L = ZM_hnf(CM_L);
    1522             :     }
    1523             : 
    1524          84 :     if (i <= r && i*rec < n0)
    1525             :     {
    1526             :       pari_timer ti;
    1527          42 :       if (DEBUGLEVEL>2) timer_start(&ti);
    1528          42 :       list = nf_chk_factors(T, P, Q_div_to_int(CM_L,utoipos(C)), famod, L->pk);
    1529          42 :       if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
    1530          42 :       if (list) break;
    1531             :     }
    1532          56 :     if (gc_needed(av,1))
    1533             :     {
    1534           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"nf_LLL_cmbf");
    1535           0 :       gerepileall(av, L->Tpk? 9: 8,
    1536             :                       &CM_L,&TT,&Tra,&famod,&L->GSmin,&L->pk,&L->prk,&L->iprk,
    1537             :                       &L->Tpk);
    1538             :     }
    1539          56 :     else CM_L = gerepilecopy(av2, CM_L);
    1540             :   }
    1541          49 :   if (DEBUGLEVEL>2)
    1542           0 :     err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
    1543          49 :   return list;
    1544             : }
    1545             : 
    1546             : static GEN
    1547        1064 : nf_combine_factors(nfcmbf_t *T, GEN polred, long klim)
    1548             : {
    1549        1064 :   nflift_t *L = T->L;
    1550             :   GEN res;
    1551             :   long maxK;
    1552             :   int done;
    1553             :   pari_timer ti;
    1554             : 
    1555        1064 :   if (DEBUGLEVEL>2) timer_start(&ti);
    1556        1064 :   T->fact = ZqX_liftfact(polred, T->fact, L->Tpk, L->pk, L->p, L->k);
    1557        1064 :   if (DEBUGLEVEL>2) timer_printf(&ti, "Hensel lift");
    1558        1064 :   res = nfcmbf(T, klim, &maxK, &done);
    1559        1064 :   if (DEBUGLEVEL>2) timer_printf(&ti, "Naive recombination");
    1560        1064 :   if (!done)
    1561             :   {
    1562          49 :     long l = lg(res)-1;
    1563             :     GEN v;
    1564          49 :     if (l > 1)
    1565             :     {
    1566             :       GEN den;
    1567          14 :       T->pol = gel(res,l);
    1568          14 :       T->polbase = Q_remove_denom(RgX_to_nfX(T->nf, T->pol), &den);
    1569          14 :       if (den) { T->Br = gmul(T->Br, den); T->pol = RgX_Rg_mul(T->pol, den); }
    1570             :     }
    1571          49 :     v = nf_LLL_cmbf(T, maxK);
    1572             :     /* remove last elt, possibly unfactored. Add all new ones. */
    1573          49 :     setlg(res, l); res = shallowconcat(res, v);
    1574             :   }
    1575        1064 :   return res;
    1576             : }
    1577             : 
    1578             : static GEN
    1579        2353 : nf_DDF_roots(GEN pol, GEN polred, GEN nfpol, long fl, nflift_t *L)
    1580             : {
    1581             :   GEN z, Cltx_r, ltdn;
    1582             :   long i, m, lz;
    1583             :   div_data D;
    1584             : 
    1585        2353 :   init_div_data(&D, pol, L);
    1586        2353 :   ltdn = mul_content(D.lt, L->dn);
    1587        2353 :   z = ZqX_roots(polred, L->Tpk, L->p, L->k);
    1588        2353 :   Cltx_r = deg1pol_shallow(D.Clt? D.Clt: gen_1, NULL, varn(pol));
    1589        2353 :   lz = lg(z);
    1590        2353 :   if (DEBUGLEVEL > 3) err_printf("Checking %ld roots:",lz-1);
    1591        7271 :   for (m=1,i=1; i<lz; i++)
    1592             :   {
    1593        4918 :     GEN r = gel(z,i);
    1594             :     int dvd;
    1595             :     pari_sp av;
    1596        4918 :     if (DEBUGLEVEL > 3) err_printf(" %ld",i);
    1597             :     /* lt*dn*topowden * r = Clt * r */
    1598        4918 :     r = nf_bestlift_to_pol(ltdn? gmul(ltdn,r): r, NULL, L);
    1599        4918 :     av = avma;
    1600        4918 :     gel(Cltx_r,2) = gneg(r); /* check P(r) == 0 */
    1601        4918 :     dvd = ZXQX_dvd(D.C2ltpol, Cltx_r, nfpol); /* integral */
    1602        4918 :     set_avma(av);
    1603             :     /* don't go on with q, usually much larger that C2ltpol */
    1604        4918 :     if (dvd) {
    1605        4603 :       if (D.Clt) r = gdiv(r, D.Clt);
    1606        4603 :       gel(z,m++) = r;
    1607             :     }
    1608         315 :     else if (fl == ROOTS_SPLIT) return cgetg(1, t_VEC);
    1609             :   }
    1610        2353 :   if (DEBUGLEVEL > 3) err_printf(" done\n");
    1611        2353 :   z[0] = evaltyp(t_VEC) | evallg(m);
    1612        2353 :   return z;
    1613             : }
    1614             : 
    1615             : /* returns a few factors of T in Fp of degree <= maxf, NULL if none exist */
    1616             : static GEN
    1617       86994 : get_good_factor(GEN T, ulong p, long maxf)
    1618             : {
    1619       86994 :   pari_sp av = avma;
    1620       86994 :   GEN r, R = gel(Flx_factor(ZX_to_Flx(T,p),p), 1);
    1621       86994 :   if (maxf == 1)
    1622             :   { /* degree 1 factors are best */
    1623       83599 :     r = gel(R,1);
    1624       83599 :     if (degpol(r) == 1) return mkvec(r);
    1625             :   }
    1626             :   else
    1627             :   { /* otherwise, pick factor of largish degree */
    1628        3395 :     long i, j, dr, d = 0, l = lg(R);
    1629             :     GEN v;
    1630        3395 :     if (l == 2) return mkvec(gel(R,1)); /* inert is fine */
    1631        3031 :     v = cgetg(l, t_VEC);
    1632       21399 :     for (i = j = 1; i < l; i++)
    1633             :     {
    1634       18949 :       r = gel(R,i); dr = degpol(r);
    1635       18949 :       if (dr > maxf) break;
    1636       18368 :       if (dr != d) { gel(v,j++) = r; d = dr; }
    1637             :     }
    1638        3031 :     setlg(v,j); if (j > 1) return v;
    1639             :   }
    1640       55927 :   return gc_NULL(av); /* failure */
    1641             : }
    1642             : 
    1643             : /* n = number of modular factors, f = residue degree; nold/fold current best
    1644             :  * return 1 if new values are "better" than old ones */
    1645             : static int
    1646       24107 : record(long nold, long n, long fold, long f)
    1647             : {
    1648       24107 :   if (!nold) return 1; /* uninitialized */
    1649       19479 :   if (fold == f) return n < nold;
    1650             :   /* if f increases, allow increasing n a little */
    1651        1939 :   if (fold < f) return n <= 20 || n < 1.1*nold;
    1652             :   /* f decreases, only allow if decreasing n a lot */
    1653         840 :   return n < 0.7*nold;
    1654             : }
    1655             : /* Optimization problem: factorization of polynomials over large Fq is slow,
    1656             :  * BUT bestlift correspondingly faster.
    1657             :  * Return maximal residue degree to be considered when picking a prime ideal */
    1658             : static long
    1659        6037 : get_maxf(long nfdeg)
    1660             : {
    1661        6037 :   long maxf = 1;
    1662        6037 :   if      (nfdeg >= 45) maxf =32;
    1663        6023 :   else if (nfdeg >= 30) maxf =16;
    1664        6009 :   else if (nfdeg >= 15) maxf = 8;
    1665        6037 :   return maxf;
    1666             : }
    1667             : /* number of maximal ideals to test before settling on best prime and number
    1668             :  * of factors; B = [K:Q]*deg(P) */
    1669             : static long
    1670        5608 : get_nbprimes(long B)
    1671             : {
    1672        5608 :   if (B <= 128) return 5;
    1673         197 :   if (B <= 1024) return 20;
    1674          35 :   if (B <= 2048) return 65;
    1675           0 :   return 100;
    1676             : }
    1677             : /* Select a prime ideal pr over which to factor pol.
    1678             :  * Return the number of factors (or roots, according to flag fl) mod pr.
    1679             :  * Set:
    1680             :  *   lt: leading term of polbase (t_INT or NULL [ for 1 ])
    1681             :  *   pr: a suitable maximal ideal
    1682             :  *   Fa: factors found mod pr
    1683             :  *   Tp: polynomial defining Fq/Fp */
    1684             : static long
    1685        5608 : nf_pick_prime(GEN nf, GEN pol, long fl, GEN *lt, GEN *Tp, ulong *pp)
    1686             : {
    1687        5608 :   GEN nfpol = nf_get_pol(nf), bad = mulii(nf_get_disc(nf), nf_get_index(nf));
    1688        5608 :   long nfdeg = degpol(nfpol), dpol = degpol(pol), nold = 0, fold = 1;
    1689        5608 :   long maxf = get_maxf(nfdeg), ct = get_nbprimes(nfdeg * dpol);
    1690             :   ulong p;
    1691             :   forprime_t S;
    1692             :   pari_timer ti_pr;
    1693             : 
    1694        5608 :   if (DEBUGLEVEL>3) timer_start(&ti_pr);
    1695        5608 :   *lt  = leading_coeff(pol); /* t_INT */
    1696        5608 :   if (gequal1(*lt)) *lt = NULL;
    1697        5608 :   *pp = 0;
    1698        5608 :   *Tp = NULL;
    1699        5608 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
    1700             :   /* select pr such that pol has the smallest number of factors, ct attempts */
    1701        5608 :   while ((p = u_forprime_next(&S)))
    1702             :   {
    1703             :     GEN vT;
    1704       92594 :     long n, i, l, ok = 0;
    1705       92594 :     ulong ltp = 0;
    1706             : 
    1707       92594 :     if (! umodiu(bad,p)) continue;
    1708       84708 :     if (*lt) { ltp = umodiu(*lt, p); if (!ltp) continue; }
    1709       83567 :     vT = get_good_factor(nfpol, p, maxf);
    1710       83567 :     if (!vT) continue;
    1711       30638 :     l = lg(vT);
    1712       59792 :     for (i = 1; i < l; i++)
    1713             :     {
    1714       31352 :       pari_sp av2 = avma;
    1715       31352 :       GEN T = gel(vT,i), red = RgX_to_FlxqX(pol, T, p);
    1716       31352 :       long f = degpol(T);
    1717       31352 :       if (f == 1)
    1718             :       { /* degree 1 */
    1719       28160 :         red = FlxX_to_Flx(red);
    1720       28160 :         if (ltp) red = Flx_normalize(red, p);
    1721       28160 :         if (!Flx_is_squarefree(red, p)) { set_avma(av2); continue; }
    1722       23407 :         ok = 1;
    1723       23407 :         n = (fl == FACTORS)? Flx_nbfact(red,p): Flx_nbroots(red,p);
    1724             :       }
    1725             :       else
    1726             :       {
    1727        3192 :         if (ltp) red = FlxqX_normalize(red, T, p);
    1728        3192 :         if (!FlxqX_is_squarefree(red, T, p)) { set_avma(av2); continue; }
    1729        2898 :         ok = 1;
    1730        2898 :         n = (fl == FACTORS)? FlxqX_nbfact(red,T,p): FlxqX_nbroots(red,T,p);
    1731             :       }
    1732       26305 :       if (fl == ROOTS_SPLIT && n < dpol) return n; /* not split */
    1733       26263 :       if (n <= 1)
    1734             :       {
    1735        6258 :         if (fl == FACTORS) return n; /* irreducible */
    1736        6027 :         if (!n) return 0; /* no root */
    1737             :       }
    1738       24114 :       if (DEBUGLEVEL>3)
    1739           0 :         err_printf("%3ld %s at prime (%ld,x^%ld+...)\n Time: %ld\n",
    1740             :             n, (fl == FACTORS)? "factors": "roots", p,f, timer_delay(&ti_pr));
    1741             : 
    1742       24114 :       if (fl == ROOTS && f==nfdeg) { *Tp = T; *pp = p; return n; }
    1743       24107 :       if (record(nold, n, fold, f)) { nold = n; fold = f; *Tp = T; *pp = p; }
    1744       18387 :       else set_avma(av2);
    1745             :     }
    1746       28440 :     if (ok && --ct <= 0) break;
    1747             :   }
    1748        3410 :   if (!nold) pari_err_OVERFLOW("nf_pick_prime [ran out of primes]");
    1749        3410 :   return nold;
    1750             : }
    1751             : 
    1752             : /* Assume lt(T) is a t_INT and T square free. Return t_VEC of irred. factors */
    1753             : static GEN
    1754         210 : nfsqff_trager(GEN u, GEN T, GEN dent)
    1755             : {
    1756         210 :   long k = 0, i, lx;
    1757         210 :   GEN U, P, x0, mx0, fa, n = ZX_ZXY_rnfequation(T, u, &k);
    1758             :   int tmonic;
    1759         210 :   if (DEBUGLEVEL>4) err_printf("nfsqff_trager: choosing k = %ld\n",k);
    1760             : 
    1761             :   /* n guaranteed to be squarefree */
    1762         210 :   fa = ZX_DDF(Q_primpart(n)); lx = lg(fa);
    1763         210 :   if (lx == 2) return mkvec(u);
    1764             : 
    1765         126 :   tmonic = is_pm1(leading_coeff(T));
    1766         126 :   P = cgetg(lx,t_VEC);
    1767         126 :   x0 = deg1pol_shallow(stoi(-k), gen_0, varn(T));
    1768         126 :   mx0 = deg1pol_shallow(stoi(k), gen_0, varn(T));
    1769         126 :   U = RgXQX_translate(u, mx0, T);
    1770         126 :   if (!tmonic) U = Q_primpart(U);
    1771         490 :   for (i=lx-1; i>0; i--)
    1772             :   {
    1773         364 :     GEN f = gel(fa,i), F = nfgcd(U, f, T, dent);
    1774         364 :     F = RgXQX_translate(F, x0, T);
    1775             :     /* F = gcd(f, u(t - x0)) [t + x0] = gcd(f(t + x0), u), more efficient */
    1776         364 :     if (typ(F) != t_POL || degpol(F) == 0)
    1777           0 :       pari_err_IRREDPOL("factornf [modulus]",T);
    1778         364 :     gel(P,i) = QXQX_normalize(F, T);
    1779             :   }
    1780         126 :   gen_sort_inplace(P, (void*)&cmp_RgX, &gen_cmp_RgX, NULL);
    1781         126 :   return P;
    1782             : }
    1783             : 
    1784             : /* Factor polynomial a on the number field defined by polynomial T, using
    1785             :  * Trager's trick */
    1786             : GEN
    1787          14 : polfnf(GEN a, GEN T)
    1788             : {
    1789          14 :   GEN rep = cgetg(3, t_MAT), A, B, y, dent, bad;
    1790             :   long dA;
    1791             :   int tmonic;
    1792             : 
    1793          14 :   if (typ(a)!=t_POL) pari_err_TYPE("polfnf",a);
    1794          14 :   if (typ(T)!=t_POL) pari_err_TYPE("polfnf",T);
    1795          14 :   T = Q_primpart(T); tmonic = is_pm1(leading_coeff(T));
    1796          14 :   RgX_check_ZX(T,"polfnf");
    1797          14 :   A = Q_primpart( QXQX_normalize(RgX_nffix("polfnf",T,a,1), T) );
    1798          14 :   dA = degpol(A);
    1799          14 :   if (dA <= 0)
    1800             :   {
    1801           0 :     set_avma((pari_sp)(rep + 3));
    1802           0 :     return (dA == 0)? trivial_fact(): zerofact(varn(A));
    1803             :   }
    1804          14 :   bad = dent = ZX_disc(T);
    1805          14 :   if (tmonic) dent = indexpartial(T, dent);
    1806          14 :   (void)nfgcd_all(A,RgX_deriv(A), T, dent, &B);
    1807          14 :   if (degpol(B) != dA) B = Q_primpart( QXQX_normalize(B, T) );
    1808          14 :   ensure_lt_INT(B);
    1809          14 :   y = nfsqff_trager(B, T, dent);
    1810          14 :   fact_from_sqff(rep, A, B, y, T, bad);
    1811          14 :   return sort_factor_pol(rep, cmp_RgX);
    1812             : }
    1813             : 
    1814             : static int
    1815       11874 : nfsqff_use_Trager(long n, long dpol)
    1816             : {
    1817       11874 :   return dpol*3<n;
    1818             : }
    1819             : 
    1820             : /* return the factorization of the square-free polynomial pol. Not memory-clean
    1821             :    The coeffs of pol are in Z_nf and its leading term is a rational integer.
    1822             :    deg(pol) > 0, deg(nfpol) > 1
    1823             :    fl is either FACTORS (return factors), or ROOTS / ROOTS_SPLIT (return roots):
    1824             :      - ROOTS, return only the roots of x in nf
    1825             :      - ROOTS_SPLIT, as ROOTS if pol splits, [] otherwise
    1826             :    den is usually 1, otherwise nf.zk is doubtful, and den bounds the
    1827             :    denominator of an arbitrary element of Z_nf on nf.zk */
    1828             : static GEN
    1829        7449 : nfsqff(GEN nf, GEN pol, long fl, GEN den)
    1830             : {
    1831        7449 :   long n, nbf, dpol = degpol(pol);
    1832             :   GEN C0, polbase;
    1833        7449 :   GEN N2, res, polred, lt, nfpol = typ(nf)==t_POL?nf:nf_get_pol(nf);
    1834             :   ulong pp;
    1835             :   nfcmbf_t T;
    1836             :   nflift_t L;
    1837             :   pari_timer ti, ti_tot;
    1838             : 
    1839        7449 :   if (DEBUGLEVEL>2) { timer_start(&ti); timer_start(&ti_tot); }
    1840        7449 :   n = degpol(nfpol);
    1841             :   /* deg = 1 => irreducible */
    1842        7449 :   if (dpol == 1) {
    1843        1645 :     if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol));
    1844        1610 :     return mkvec(gneg(gdiv(gel(pol,2),gel(pol,3))));
    1845             :   }
    1846        5804 :   if (typ(nf)==t_POL || nfsqff_use_Trager(n,dpol))
    1847             :   {
    1848             :     GEN z;
    1849         196 :     if (DEBUGLEVEL>2) err_printf("Using Trager's method\n");
    1850         196 :     if (typ(nf) != t_POL) den =  mulii(den, nf_get_index(nf));
    1851         196 :     z = nfsqff_trager(Q_primpart(pol), nfpol, den);
    1852         196 :     if (fl != FACTORS) {
    1853         168 :       long i, l = lg(z);
    1854         378 :       for (i = 1; i < l; i++)
    1855             :       {
    1856         301 :         GEN LT, t = gel(z,i); if (degpol(t) > 1) break;
    1857         210 :         LT = gel(t,3);
    1858         210 :         if (typ(LT) == t_POL) LT = gel(LT,2); /* constant */
    1859         210 :         gel(z,i) = gdiv(gel(t,2), negi(LT));
    1860             :       }
    1861         168 :       setlg(z, i);
    1862         168 :       if (fl == ROOTS_SPLIT && i != l) return cgetg(1,t_VEC);
    1863             :     }
    1864         196 :     return z;
    1865             :   }
    1866             : 
    1867        5608 :   polbase = RgX_to_nfX(nf, pol);
    1868        5608 :   nbf = nf_pick_prime(nf, pol, fl, &lt, &L.Tp, &pp);
    1869        5608 :   if (L.Tp)
    1870             :   {
    1871        4635 :     L.Tp = Flx_to_ZX(L.Tp);
    1872        4635 :     L.p = utoi(pp);
    1873             :   }
    1874             : 
    1875        5608 :   if (fl == ROOTS_SPLIT && nbf < dpol) return cgetg(1,t_VEC);
    1876        5566 :   if (nbf <= 1)
    1877             :   {
    1878        2954 :     if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol)); /* irred. */
    1879        2723 :     if (!nbf) return cgetg(1,t_VEC); /* no root */
    1880             :   }
    1881             : 
    1882        3417 :   if (DEBUGLEVEL>2) {
    1883           0 :     timer_printf(&ti, "choice of a prime ideal");
    1884           0 :     err_printf("Prime ideal chosen: (%lu,x^%ld+...)\n", pp, degpol(L.Tp));
    1885             :   }
    1886        3417 :   L.tozk = nf_get_invzk(nf);
    1887        3417 :   L.topow= nf_get_zkprimpart(nf);
    1888        3417 :   L.topowden = nf_get_zkden(nf);
    1889        3417 :   if (is_pm1(den)) den = NULL;
    1890        3417 :   L.dn = den;
    1891        3417 :   T.ZC = L2_bound(nf, den);
    1892        3417 :   T.Br = nf_root_bounds(nf, pol); if (lt) T.Br = gmul(T.Br, lt);
    1893             : 
    1894             :   /* C0 = bound for T_2(Q_i), Q | P */
    1895        3417 :   if (fl != FACTORS) C0 = normTp(T.Br, 2, n);
    1896        1064 :   else               C0 = nf_factor_bound(nf, polbase);
    1897        3417 :   T.bound = mulrr(T.ZC, C0); /* bound for |Q_i|^2 in Z^n on chosen Z-basis */
    1898             : 
    1899        3417 :   N2 = mulur(dpol*dpol, normTp(T.Br, 4, n)); /* bound for T_2(lt * S_2) */
    1900        3417 :   T.BS_2 = mulrr(T.ZC, N2); /* bound for |S_2|^2 on chosen Z-basis */
    1901             : 
    1902        3417 :   if (DEBUGLEVEL>2) {
    1903           0 :     timer_printf(&ti, "bound computation");
    1904           0 :     err_printf("  1) T_2 bound for %s: %Ps\n",
    1905             :                fl == FACTORS?"factor": "root", C0);
    1906           0 :     err_printf("  2) Conversion from T_2 --> | |^2 bound : %Ps\n", T.ZC);
    1907           0 :     err_printf("  3) Final bound: %Ps\n", T.bound);
    1908             :   }
    1909             : 
    1910        3417 :   bestlift_init(0, nf, T.bound, &L);
    1911        3417 :   if (DEBUGLEVEL>2) timer_start(&ti);
    1912        3417 :   polred = ZqX_normalize(polbase, lt, &L); /* monic */
    1913             : 
    1914        3417 :   if (fl != FACTORS) {
    1915        2353 :     GEN z = nf_DDF_roots(pol, polred, nfpol, fl, &L);
    1916        2353 :     if (lg(z) == 1) return cgetg(1, t_VEC);
    1917        2297 :     return z;
    1918             :   }
    1919             : 
    1920        1064 :   T.fact = gel(FqX_factor(polred, L.Tp, L.p), 1);
    1921        1064 :   if (DEBUGLEVEL>2)
    1922           0 :     timer_printf(&ti, "splitting mod %Ps^%ld", L.p, degpol(L.Tp));
    1923        1064 :   T.L  = &L;
    1924        1064 :   T.polbase = polbase;
    1925        1064 :   T.pol   = pol;
    1926        1064 :   T.nf    = nf;
    1927        1064 :   res = nf_combine_factors(&T, polred, dpol-1);
    1928        1064 :   if (DEBUGLEVEL>2)
    1929           0 :     err_printf("Total Time: %ld\n===========\n", timer_delay(&ti_tot));
    1930        1064 :   return res;
    1931             : }
    1932             : 
    1933             : /* assume pol monic in nf.zk[X] */
    1934             : GEN
    1935          84 : nfroots_if_split(GEN *pnf, GEN pol)
    1936             : {
    1937          84 :   GEN T = get_nfpol(*pnf,pnf), den = fix_nf(pnf, &T, &pol);
    1938          84 :   pari_sp av = avma;
    1939          84 :   GEN z = nfsqff(*pnf, pol, ROOTS_SPLIT, den);
    1940          84 :   if (lg(z) == 1) return gc_NULL(av);
    1941          42 :   return gerepilecopy(av, z);
    1942             : }
    1943             : 
    1944             : /*******************************************************************/
    1945             : /*                                                                 */
    1946             : /*              Roots of unity in a number field                   */
    1947             : /*     (alternative to nfrootsof1 using factorization in K[X])     */
    1948             : /*                                                                 */
    1949             : /*******************************************************************/
    1950             : /* Code adapted from nffactor. Structure of the algorithm; only step 1 is
    1951             :  * specific to roots of unity.
    1952             :  *
    1953             :  * [Step 1]: guess roots via ramification. If trivial output this.
    1954             :  * [Step 2]: select prime [p] unramified and ideal [pr] above
    1955             :  * [Step 3]: evaluate the maximal exponent [k] such that the fondamental domain
    1956             :  *           of a LLL-reduction of [prk] = pr^k contains a ball of radius larger
    1957             :  *           than the norm of any root of unity.
    1958             :  * [Step 3]: select a heuristic exponent,
    1959             :  *           LLL reduce prk=pr^k and verify the exponent is sufficient,
    1960             :  *           otherwise try a larger one.
    1961             :  * [Step 4]: factor the cyclotomic polynomial mod [pr],
    1962             :  *           Hensel lift to pr^k and find the representative in the ball
    1963             :  *           If there is it is a primitive root */
    1964             : 
    1965             : /* Choose prime ideal unramified with "large" inertia degree */
    1966             : static void
    1967         429 : nf_pick_prime_for_units(GEN nf, nflift_t *L)
    1968             : {
    1969         429 :   GEN nfpol = nf_get_pol(nf), bad = mulii(nf_get_disc(nf), nf_get_index(nf));
    1970         429 :   GEN ap = NULL, r = NULL;
    1971         429 :   long nfdeg = degpol(nfpol), maxf = get_maxf(nfdeg);
    1972             :   ulong pp;
    1973             :   forprime_t S;
    1974             : 
    1975         429 :   (void)u_forprime_init(&S, 2, ULONG_MAX);
    1976         429 :   while ( (pp = u_forprime_next(&S)) )
    1977             :   {
    1978        4577 :     if (! umodiu(bad,pp)) continue;
    1979        3427 :     r = get_good_factor(nfpol, pp, maxf);
    1980        3427 :     if (r) break;
    1981             :   }
    1982         429 :   if (!r) pari_err_OVERFLOW("nf_pick_prime [ran out of primes]");
    1983         429 :   r = gel(r,lg(r)-1); /* largest inertia degree */
    1984         429 :   ap = utoipos(pp);
    1985         429 :   L->p = ap;
    1986         429 :   L->Tp = Flx_to_ZX(r);
    1987         429 :   L->tozk = nf_get_invzk(nf);
    1988         429 :   L->topow = nf_get_zkprimpart(nf);
    1989         429 :   L->topowden = nf_get_zkden(nf);
    1990         429 : }
    1991             : 
    1992             : /* *Heuristic* exponent k such that the fundamental domain of pr^k
    1993             :  * should contain the ball of radius C */
    1994             : static double
    1995         429 : mybestlift_bound(GEN C)
    1996             : {
    1997         429 :   C = gtofp(C,DEFAULTPREC);
    1998         429 :   return ceil(log(gtodouble(C)) / 0.2) + 3;
    1999             : }
    2000             : 
    2001             : /* simplified nf_DDF_roots: polcyclo(n) monic in ZX either splits or has no
    2002             :  * root in nf.
    2003             :  * Return a root or NULL (no root) */
    2004             : static GEN
    2005         443 : nfcyclo_root(long n, GEN nfpol, nflift_t *L)
    2006             : {
    2007         443 :   GEN q, r, Cltx_r, pol = polcyclo(n,0), gn = utoipos(n);
    2008             :   div_data D;
    2009             : 
    2010         443 :   init_div_data(&D, pol, L);
    2011         443 :   (void)Fq_sqrtn(gen_1, gn, L->Tp, L->p, &r);
    2012             :   /* r primitive n-th root of 1 in Fq */
    2013         443 :   r = Zq_sqrtnlift(gen_1, gn, r, L->Tpk, L->p, L->k);
    2014             :   /* lt*dn*topowden * r = Clt * r */
    2015         443 :   r = nf_bestlift_to_pol(r, NULL, L);
    2016         443 :   Cltx_r = deg1pol_shallow(D.Clt? D.Clt: gen_1, gneg(r), varn(pol));
    2017             :   /* check P(r) == 0 */
    2018         443 :   q = RgXQX_divrem(D.C2ltpol, Cltx_r, nfpol, ONLY_DIVIDES); /* integral */
    2019         443 :   if (!q) return NULL;
    2020         415 :   if (D.Clt) r = gdiv(r, D.Clt);
    2021         415 :   return r;
    2022             : }
    2023             : 
    2024             : /* Guesses the number of roots of unity in number field [nf].
    2025             :  * Computes gcd of N(P)-1 for some primes. The value returned is a proven
    2026             :  * multiple of the correct value. */
    2027             : static long
    2028        8150 : guess_roots(GEN nf)
    2029             : {
    2030        8150 :   long c = 0, nfdegree = nf_get_degree(nf), B = nfdegree + 20, l;
    2031        8150 :   ulong p = 2;
    2032        8150 :   GEN T = nf_get_pol(nf), D = nf_get_disc(nf), index = nf_get_index(nf);
    2033        8150 :   GEN nbroots = NULL;
    2034             :   forprime_t S;
    2035             :   pari_sp av;
    2036             : 
    2037        8150 :   (void)u_forprime_init(&S, 3, ULONG_MAX);
    2038        8150 :   av = avma;
    2039             :   /* result must be stationary (counter c) for at least B loops */
    2040      224989 :   for (l=1; (p = u_forprime_next(&S)); l++)
    2041             :   {
    2042             :     GEN old, F, pf_1, Tp;
    2043      224989 :     long i, nb, gcdf = 0;
    2044             : 
    2045      224989 :     if (!umodiu(D,p) || !umodiu(index,p)) continue;
    2046      215535 :     Tp = ZX_to_Flx(T,p); /* squarefree */
    2047      215535 :     F = Flx_nbfact_by_degree(Tp, &nb, p);
    2048             :     /* the gcd of the p^f - 1 is p^(gcd of the f's) - 1 */
    2049      647032 :     for (i = 1; i <= nfdegree; i++)
    2050      530443 :       if (F[i]) {
    2051      215913 :         gcdf = gcdf? ugcd(gcdf, i): i;
    2052      215913 :         if (gcdf == 1) break;
    2053             :       }
    2054      215535 :     pf_1 = subiu(powuu(p, gcdf), 1);
    2055      215535 :     old = nbroots;
    2056      215535 :     nbroots = nbroots? gcdii(pf_1, nbroots): pf_1;
    2057      215535 :     if (DEBUGLEVEL>5)
    2058           0 :       err_printf("p=%lu; gcf(f(P/p))=%ld; nbroots | %Ps",p, gcdf, nbroots);
    2059             :     /* if same result go on else reset the stop counter [c] */
    2060      215535 :     if (old && equalii(nbroots,old))
    2061      200136 :     { if (!is_bigint(nbroots) && ++c > B) break; }
    2062             :     else
    2063       15399 :       c = 0;
    2064             :   }
    2065        8150 :   if (!nbroots) pari_err_OVERFLOW("guess_roots [ran out of primes]");
    2066        8150 :   if (DEBUGLEVEL>5) err_printf("%ld loops\n",l);
    2067        8150 :   return gc_long(av, itos(nbroots));
    2068             : }
    2069             : 
    2070             : /* T(x) an irreducible ZX. Is it of the form Phi_n(c \pm x) ?
    2071             :  * Return NULL if not, and a root of 1 of maximal order in Z[x]/(T) otherwise
    2072             :  *
    2073             :  * N.B. Set n_squarefree = 1 if n is squarefree, and 0 otherwise.
    2074             :  * This last parameter is inconvenient, but it allows a cheap
    2075             :  * stringent test. (n guessed from guess_roots())*/
    2076             : static GEN
    2077        1246 : ZXirred_is_cyclo_translate(GEN T, long n_squarefree)
    2078             : {
    2079        1246 :   long r, m, d = degpol(T);
    2080        1246 :   GEN T1, c = divis_rem(gel(T, d+1), d, &r); /* d-1 th coeff divisible by d ? */
    2081             :   /* The trace coefficient of polcyclo(n) is \pm1 if n is square free, and 0
    2082             :    * otherwise. */
    2083        1246 :   if (!n_squarefree)
    2084         623 :   { if (r) return NULL; }
    2085             :   else
    2086             :   {
    2087         623 :     if (r < -1)
    2088             :     {
    2089           0 :       r += d;
    2090           0 :       c = subiu(c, 1);
    2091             :     }
    2092         623 :     else if (r == d-1)
    2093             :     {
    2094          35 :       r = -1;
    2095          35 :       c = addiu(c, 1);
    2096             :     }
    2097         623 :     if (r != 1 && r != -1) return NULL;
    2098             :   }
    2099        1211 :   if (signe(c)) /* presumably Phi_guess(c \pm x) */
    2100          35 :     T = RgX_translate(T, negi(c));
    2101        1211 :   if (!n_squarefree) T = RgX_deflate_max(T, &m);
    2102             :   /* presumably Phi_core(guess)(\pm x), cyclotomic iff original T was */
    2103        1211 :   T1 = ZX_graeffe(T);
    2104        1211 :   if (ZX_equal(T, T1)) /* T = Phi_n, n odd */
    2105          35 :     return deg1pol_shallow(gen_m1, negi(c), varn(T));
    2106        1176 :   else if (ZX_equal(T1, ZX_z_unscale(T, -1))) /* T = Phi_{2n}, nodd */
    2107        1155 :     return deg1pol_shallow(gen_1, c, varn(T));
    2108          21 :   return NULL;
    2109             : }
    2110             : 
    2111             : static GEN
    2112        9911 : trivroots(void) { return mkvec2(gen_2, gen_m1); }
    2113             : /* Number of roots of unity in number field [nf]. */
    2114             : GEN
    2115       11530 : rootsof1(GEN nf)
    2116             : {
    2117             :   nflift_t L;
    2118             :   GEN T, q, fa, LP, LE, C0, z, disc;
    2119             :   pari_timer ti;
    2120             :   long i, l, nbguessed, nbroots, nfdegree;
    2121             :   pari_sp av;
    2122             : 
    2123       11530 :   nf = checknf(nf);
    2124       11530 :   if (nf_get_r1(nf)) return trivroots();
    2125             : 
    2126             :   /* Step 1 : guess number of roots and discard trivial case 2 */
    2127        8150 :   if (DEBUGLEVEL>2) timer_start(&ti);
    2128        8150 :   nbguessed = guess_roots(nf);
    2129        8150 :   if (DEBUGLEVEL>2)
    2130           0 :     timer_printf(&ti, "guessing roots of 1 [guess = %ld]", nbguessed);
    2131        8150 :   if (nbguessed == 2) return trivroots();
    2132             : 
    2133        1619 :   nfdegree = nf_get_degree(nf);
    2134        1619 :   fa = factoru(nbguessed);
    2135        1619 :   LP = gel(fa,1); l = lg(LP);
    2136        1619 :   LE = gel(fa,2);
    2137        1619 :   disc = nf_get_disc(nf);
    2138        4205 :   for (i = 1; i < l; i++)
    2139             :   {
    2140        2586 :     long p = LP[i];
    2141             :     /* Degree and ramification test: find largest k such that Q(zeta_{p^k})
    2142             :      * may be a subfield of K. Q(zeta_p^k) has degree (p-1)p^(k-1)
    2143             :      * and v_p(discriminant) = ((p-1)k-1)p^(k-1); so we must have
    2144             :      * v_p(disc_K) >= ((p-1)k-1) * n / (p-1) = kn - q, where q = n/(p-1) */
    2145        2586 :     if (p == 2)
    2146             :     { /* the test simplifies a little in that case */
    2147             :       long v, vnf, k;
    2148        1619 :       if (LE[i] == 1) continue;
    2149         708 :       vnf = vals(nfdegree);
    2150         708 :       v = vali(disc);
    2151         736 :       for (k = minss(LE[i], vnf+1); k >= 1; k--)
    2152         736 :         if (v >= nfdegree*(k-1)) { nbguessed >>= LE[i]-k; LE[i] = k; break; }
    2153             :       /* N.B the test above always works for k = 1: LE[i] >= 1 */
    2154             :     }
    2155             :     else
    2156             :     {
    2157             :       long v, vnf, k;
    2158         967 :       ulong r, q = udivuu_rem(nfdegree, p-1, &r);
    2159         967 :       if (r) { nbguessed /= upowuu(p, LE[i]); LE[i] = 0; continue; }
    2160             :       /* q = n/(p-1) */
    2161         967 :       vnf = u_lval(q, p);
    2162         967 :       v = Z_lval(disc, p);
    2163         967 :       for (k = minss(LE[i], vnf+1); k >= 0; k--)
    2164         967 :         if (v >= nfdegree*k-(long)q)
    2165         967 :         { nbguessed /= upowuu(p, LE[i]-k); LE[i] = k; break; }
    2166             :       /* N.B the test above always works for k = 0: LE[i] >= 0 */
    2167             :     }
    2168             :   }
    2169        1619 :   if (DEBUGLEVEL>2)
    2170           0 :     timer_printf(&ti, "after ramification conditions [guess = %ld]", nbguessed);
    2171        1619 :   if (nbguessed == 2) return trivroots();
    2172        1619 :   av = avma;
    2173             : 
    2174             :   /* Step 1.5 : test if nf.pol == subst(polcyclo(nbguessed), x, \pm x+c) */
    2175        1619 :   T = nf_get_pol(nf);
    2176        1619 :   if (eulerphiu_fact(fa) == (ulong)nfdegree)
    2177             :   {
    2178        1246 :     z = ZXirred_is_cyclo_translate(T, uissquarefree_fact(fa));
    2179        1246 :     if (DEBUGLEVEL>2) timer_printf(&ti, "checking for cyclotomic polynomial");
    2180        1246 :     if (z)
    2181             :     {
    2182        1190 :       z = nf_to_scalar_or_basis(nf,z);
    2183        1190 :       return gerepilecopy(av, mkvec2(utoipos(nbguessed), z));
    2184             :     }
    2185          56 :     set_avma(av);
    2186             :   }
    2187             : 
    2188             :   /* Step 2 : choose a prime ideal for local lifting */
    2189         429 :   nf_pick_prime_for_units(nf, &L);
    2190         429 :   if (DEBUGLEVEL>2)
    2191           0 :     timer_printf(&ti, "choosing prime %Ps, degree %ld",
    2192           0 :              L.p, L.Tp? degpol(L.Tp): 1);
    2193             : 
    2194             :   /* Step 3 : compute a reduced pr^k allowing lifting of local solutions */
    2195             :   /* evaluate maximum L2 norm of a root of unity in nf */
    2196         429 :   C0 = gmulsg(nfdegree, L2_bound(nf, gen_1));
    2197             :   /* lift and reduce pr^k */
    2198         429 :   if (DEBUGLEVEL>2) err_printf("Lift pr^k; GSmin wanted: %Ps\n",C0);
    2199         429 :   bestlift_init((long)mybestlift_bound(C0), nf, C0, &L);
    2200         429 :   L.dn = NULL;
    2201         429 :   if (DEBUGLEVEL>2) timer_start(&ti);
    2202             : 
    2203             :   /* Step 4 : actual computation of roots */
    2204         429 :   nbroots = 2; z = gen_m1;
    2205         429 :   q = powiu(L.p,degpol(L.Tp));
    2206        1216 :   for (i = 1; i < l; i++)
    2207             :   { /* for all prime power factors of nbguessed, find a p^k-th root of unity */
    2208         787 :     long k, p = LP[i];
    2209        1159 :     for (k = minss(LE[i], Z_lval(subiu(q,1UL),p)); k > 0; k--)
    2210             :     { /* find p^k-th roots */
    2211         787 :       pari_sp av = avma;
    2212         787 :       long pk = upowuu(p,k);
    2213             :       GEN r;
    2214         787 :       if (pk==2) continue; /* no need to test second roots ! */
    2215         443 :       r = nfcyclo_root(pk, T, &L);
    2216         443 :       if (DEBUGLEVEL>2) timer_printf(&ti, "for factoring Phi_%ld^%ld", p,k);
    2217         443 :       if (r) {
    2218         415 :         if (DEBUGLEVEL>2) err_printf("  %s root of unity found\n",uordinal(pk));
    2219         415 :         if (p==2) { nbroots = pk; z = r; }
    2220         337 :         else     { nbroots *= pk; z = nfmul(nf, z,r); }
    2221         415 :         break;
    2222             :       }
    2223          28 :       set_avma(av);
    2224          28 :       if (DEBUGLEVEL) pari_warn(warner,"rootsof1: wrong guess");
    2225             :     }
    2226             :   }
    2227         429 :   return gerepilecopy(av, mkvec2(utoi(nbroots), z));
    2228             : }
    2229             : 
    2230             : static long
    2231           0 : zk_equal1(GEN y)
    2232             : {
    2233           0 :   if (typ(y) == t_INT) return equali1(y);
    2234           0 :   return equali1(gel(y,1)) && ZV_isscalar(y);
    2235             : }
    2236             : /* x^w = 1 */
    2237             : static GEN
    2238           0 : is_primitive_root(GEN nf, GEN fa, GEN x, long w)
    2239             : {
    2240           0 :   GEN P = gel(fa,1);
    2241           0 :   long i, l = lg(P);
    2242             : 
    2243           0 :   for (i = 1; i < l; i++)
    2244             :   {
    2245           0 :     long p = itos(gel(P,i));
    2246           0 :     GEN y = nfpow_u(nf,x, w/p);
    2247           0 :     if (zk_equal1(y) > 0) /* y = 1 */
    2248             :     {
    2249           0 :       if (p != 2 || !equali1(gcoeff(fa,i,2))) return NULL;
    2250           0 :       x = gneg_i(x);
    2251             :     }
    2252             :   }
    2253           0 :   return x;
    2254             : }
    2255             : GEN
    2256           0 : rootsof1_kannan(GEN nf)
    2257             : {
    2258           0 :   pari_sp av = avma;
    2259             :   long N, k, i, ws, prec;
    2260             :   GEN z, y, d, list, w;
    2261             : 
    2262           0 :   nf = checknf(nf);
    2263           0 :   if ( nf_get_r1(nf) ) return trivroots();
    2264             : 
    2265           0 :   N = nf_get_degree(nf); prec = nf_get_prec(nf);
    2266             :   for (;;)
    2267           0 :   {
    2268           0 :     GEN R = R_from_QR(nf_get_G(nf), prec);
    2269           0 :     if (R)
    2270             :     {
    2271           0 :       y = fincke_pohst(mkvec(R), utoipos(N), N * N, 0, NULL);
    2272           0 :       if (y) break;
    2273             :     }
    2274           0 :     prec = precdbl(prec);
    2275           0 :     if (DEBUGLEVEL) pari_warn(warnprec,"rootsof1",prec);
    2276           0 :     nf = nfnewprec_shallow(nf,prec);
    2277             :   }
    2278           0 :   if (itos(ground(gel(y,2))) != N) pari_err_BUG("rootsof1 (bug1)");
    2279           0 :   w = gel(y,1); ws = itos(w);
    2280           0 :   if (ws == 2) { set_avma(av); return trivroots(); }
    2281             : 
    2282           0 :   d = Z_factor(w); list = gel(y,3); k = lg(list);
    2283           0 :   for (i=1; i<k; i++)
    2284             :   {
    2285           0 :     z = is_primitive_root(nf, d, gel(list,i), ws);
    2286           0 :     if (z) return gerepilecopy(av, mkvec2(w, z));
    2287             :   }
    2288           0 :   pari_err_BUG("rootsof1");
    2289             :   return NULL; /* LCOV_EXCL_LINE */
    2290             : }

Generated by: LCOV version 1.13