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

Generated by: LCOV version 1.11