Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

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

Generated by: LCOV version 1.14