Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - nffactor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16937-4bd9b4e) Lines: 1086 1189 91.3 %
Date: 2014-10-24 Functions: 65 69 94.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 616 822 74.9 %

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

Generated by: LCOV version 1.9