Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - QX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19588-1c9967d) Lines: 731 762 95.9 %
Date: 2016-09-24 05:54:29 Functions: 41 42 97.6 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : 
      16             : /* x,y two ZX, y non constant. Return q = x/y if y divides x in Z[X] and NULL
      17             :  * otherwise. If not NULL, B is a t_INT upper bound for ||q||_oo. */
      18             : static GEN
      19     6264059 : ZX_divides_i(GEN x, GEN y, GEN B)
      20             : {
      21             :   long dx, dy, dz, i, j;
      22             :   pari_sp av;
      23             :   GEN z,p1,y_lead;
      24             : 
      25     6264059 :   dy=degpol(y);
      26     6264059 :   dx=degpol(x);
      27     6264059 :   dz=dx-dy; if (dz<0) return NULL;
      28     6263870 :   z=cgetg(dz+3,t_POL); z[1] = x[1];
      29     6263870 :   x += 2; y += 2; z += 2;
      30     6263870 :   y_lead = gel(y,dy);
      31     6263870 :   if (equali1(y_lead)) y_lead = NULL;
      32             : 
      33     6263870 :   p1 = gel(x,dx);
      34     6263870 :   if (y_lead) {
      35             :     GEN r;
      36        3474 :     p1 = dvmdii(p1,y_lead, &r);
      37        3474 :     if (r != gen_0) return NULL;
      38             :   }
      39     6260396 :   else p1 = icopy(p1);
      40     6263870 :   gel(z,dz) = p1;
      41     6962303 :   for (i=dx-1; i>=dy; i--)
      42             :   {
      43      698727 :     av = avma; p1 = gel(x,i);
      44     2371889 :     for (j=i-dy+1; j<=i && j<=dz; j++)
      45     1673162 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
      46      698727 :     if (y_lead) {
      47             :       GEN r;
      48       14126 :       p1 = dvmdii(p1,y_lead, &r);
      49       14126 :       if (r != gen_0) return NULL;
      50             :     }
      51      698692 :     if (B && abscmpii(p1, B) > 0) return NULL;
      52      698433 :     p1 = gerepileuptoint(av, p1);
      53      698433 :     gel(z,i-dy) = p1;
      54             :   }
      55     6263576 :   av = avma;
      56    15730937 :   for (; i >= 0; i--)
      57             :   {
      58     9505189 :     p1 = gel(x,i);
      59             :     /* we always enter this loop at least once */
      60    19991675 :     for (j=0; j<=i && j<=dz; j++)
      61    10486486 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
      62     9505189 :     if (signe(p1)) return NULL;
      63     9467361 :     avma = av;
      64             :   }
      65     6225748 :   return z - 2;
      66             : }
      67             : static GEN
      68     6261511 : ZX_divides(GEN x, GEN y) { return ZX_divides_i(x,y,NULL); }
      69             : 
      70             : #if 0
      71             : /* cf Beauzamy et al: upper bound for
      72             :  *      lc(x) * [2^(5/8) / pi^(3/8)] e^(1/4n) 2^(n/2) sqrt([x]_2)/ n^(3/8)
      73             :  * where [x]_2 = sqrt(\sum_i=0^n x[i]^2 / binomial(n,i)). One factor has
      74             :  * all coeffs less than then bound */
      75             : static GEN
      76             : two_factor_bound(GEN x)
      77             : {
      78             :   long i, j, n = lg(x) - 3;
      79             :   pari_sp av = avma;
      80             :   GEN *invbin, c, r = cgetr(3), z;
      81             : 
      82             :   x += 2; invbin = (GEN*)new_chunk(n+1);
      83             :   z = real_1(LOWDEFAULTPREC); /* invbin[i] = 1 / binomial(n, i) */
      84             :   for (i=0,j=n; j >= i; i++,j--)
      85             :   {
      86             :     invbin[i] = invbin[j] = z;
      87             :     z = divru(mulru(z, i+1), n-i);
      88             :   }
      89             :   z = invbin[0]; /* = 1 */
      90             :   for (i=0; i<=n; i++)
      91             :   {
      92             :     c = gel(x,i); if (!signe(c)) continue;
      93             :     affir(c, r);
      94             :     z = addrr(z, mulrr(sqrr(r), invbin[i]));
      95             :   }
      96             :   z = shiftr(sqrtr(z), n);
      97             :   z = divrr(z, dbltor(pow((double)n, 0.75)));
      98             :   z = roundr_safe(sqrtr(z));
      99             :   z = mulii(z, absi(gel(x,n)));
     100             :   return gerepileuptoint(av, shifti(z, 1));
     101             : }
     102             : #endif
     103             : 
     104             : /* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */
     105             : static GEN
     106        6063 : Mignotte_bound(GEN S)
     107             : {
     108        6063 :   long i, d = degpol(S);
     109        6063 :   GEN C, N2, t, binlS, lS = leading_coeff(S), bin = vecbinome(d-1);
     110             : 
     111        6063 :   N2 = sqrtr(RgX_fpnorml2(S,DEFAULTPREC));
     112        6063 :   binlS = is_pm1(lS)? bin: ZC_Z_mul(bin, lS);
     113             : 
     114             :   /* i = 0 */
     115        6063 :   C = gel(binlS,1);
     116             :   /* i = d */
     117        6063 :   t = N2; if (gcmp(C, t) < 0) C = t;
     118       84845 :   for (i = 1; i < d; i++)
     119             :   {
     120       78782 :     t = addri(mulir(gel(bin,i), N2), gel(binlS,i+1));
     121       78782 :     if (mpcmp(C, t) < 0) C = t;
     122             :   }
     123        6063 :   return C;
     124             : }
     125             : /* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,
     126             :  * where [P]_2 is Bombieri's 2-norm */
     127             : static GEN
     128        6063 : Beauzamy_bound(GEN S)
     129             : {
     130        6063 :   const long prec = DEFAULTPREC;
     131        6063 :   long i, d = degpol(S);
     132             :   GEN bin, lS, s, C;
     133        6063 :   bin = vecbinome(d);
     134             : 
     135        6063 :   s = real_0(prec);
     136       96971 :   for (i=0; i<=d; i++)
     137             :   {
     138       90908 :     GEN c = gel(S,i+2);
     139       90908 :     if (gequal0(c)) continue;
     140             :     /* s += P_i^2 / binomial(d,i) */
     141       74885 :     s = addrr(s, divri(itor(sqri(c), prec), gel(bin,i+1)));
     142             :   }
     143             :   /* s = [S]_2^2 */
     144        6063 :   C = powruhalf(stor(3,prec), 3 + 2*d); /* 3^{3/2 + d} */
     145        6063 :   C = divrr(mulrr(C, s), mulur(4*d, mppi(prec)));
     146        6063 :   lS = absi(leading_coeff(S));
     147        6063 :   return mulir(lS, sqrtr(C));
     148             : }
     149             : 
     150             : static GEN
     151        6063 : factor_bound(GEN S)
     152             : {
     153        6063 :   pari_sp av = avma;
     154        6063 :   GEN a = Mignotte_bound(S);
     155        6063 :   GEN b = Beauzamy_bound(S);
     156        6063 :   if (DEBUGLEVEL>2)
     157             :   {
     158           0 :     err_printf("Mignotte bound: %Ps\n",a);
     159           0 :     err_printf("Beauzamy bound: %Ps\n",b);
     160             :   }
     161        6063 :   return gerepileupto(av, ceil_safe(gmin(a, b)));
     162             : }
     163             : 
     164             : /* Naive recombination of modular factors: combine up to maxK modular
     165             :  * factors, degree <= klim
     166             :  *
     167             :  * target = polynomial we want to factor
     168             :  * famod = array of modular factors.  Product should be congruent to
     169             :  * target/lc(target) modulo p^a
     170             :  * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
     171             : static GEN
     172        4887 : cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,
     173             :      long klim, long *pmaxK, int *done)
     174             : {
     175        4887 :   long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
     176             :   ulong spa_b, spa_bs2, Sbound;
     177        4887 :   GEN lc, lcpol, pa = powiu(p,a), pas2 = shifti(pa,-1);
     178        4887 :   GEN trace1   = cgetg(lfamod+1, t_VECSMALL);
     179        4887 :   GEN trace2   = cgetg(lfamod+1, t_VECSMALL);
     180        4887 :   GEN ind      = cgetg(lfamod+1, t_VECSMALL);
     181        4887 :   GEN deg      = cgetg(lfamod+1, t_VECSMALL);
     182        4887 :   GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
     183        4887 :   GEN listmod  = cgetg(lfamod+1, t_VEC);
     184        4887 :   GEN fa       = cgetg(lfamod+1, t_VEC);
     185             : 
     186        4887 :   *pmaxK = cmbf_maxK(lfamod);
     187        4887 :   lc = absi(leading_coeff(pol));
     188        4887 :   if (is_pm1(lc)) lc = NULL;
     189        4887 :   lcpol = lc? ZX_Z_mul(pol, lc): pol;
     190             : 
     191             :   {
     192        4887 :     GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;
     193             : 
     194        4887 :     pa_b = powiu(p, a-b); /* < 2^31 */
     195        4887 :     pa_bs2 = shifti(pa_b,-1);
     196        4887 :     pb= powiu(p, b);
     197       19526 :     for (i=1; i <= lfamod; i++)
     198             :     {
     199       14639 :       GEN T1,T2, P = gel(famod,i);
     200       14639 :       long d = degpol(P);
     201             : 
     202       14639 :       deg[i] = d; P += 2;
     203       14639 :       T1 = gel(P,d-1);/* = - S_1 */
     204       14639 :       T2 = sqri(T1);
     205       14639 :       if (d > 1) T2 = subii(T2, shifti(gel(P,d-2),1));
     206       14639 :       T2 = modii(T2, pa); /* = S_2 Newton sum */
     207       14639 :       if (lc)
     208             :       {
     209         658 :         T1 = Fp_mul(lc, T1, pa);
     210         658 :         T2 = Fp_mul(lc2,T2, pa);
     211             :       }
     212       14639 :       uel(trace1,i) = itou(diviiround(T1, pb));
     213       14639 :       uel(trace2,i) = itou(diviiround(T2, pb));
     214             :     }
     215        4887 :     spa_b   = uel(pa_b,2); /* < 2^31 */
     216        4887 :     spa_bs2 = uel(pa_bs2,2); /* < 2^31 */
     217             :   }
     218        4887 :   degsofar[0] = 0; /* sentinel */
     219             : 
     220             :   /* ind runs through strictly increasing sequences of length K,
     221             :    * 1 <= ind[i] <= lfamod */
     222             : nextK:
     223        9830 :   if (K > *pmaxK || 2*K > lfamod) goto END;
     224        5930 :   if (DEBUGLEVEL > 3)
     225           0 :     err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
     226        5930 :   setlg(ind, K+1); ind[1] = 1;
     227        5930 :   Sbound = (ulong) ((K+1)>>1);
     228        5930 :   i = 1; curdeg = deg[ind[1]];
     229             :   for(;;)
     230             :   { /* try all combinations of K factors */
     231      365423 :     for (j = i; j < K; j++)
     232             :     {
     233       44639 :       degsofar[j] = curdeg;
     234       44639 :       ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
     235             :     }
     236      320784 :     if (curdeg <= klim) /* trial divide */
     237             :     {
     238             :       GEN y, q, list;
     239             :       pari_sp av;
     240             :       ulong t;
     241             : 
     242             :       /* d - 1 test */
     243      933067 :       for (t=uel(trace1,ind[1]),i=2; i<=K; i++)
     244      612283 :         t = Fl_add(t, uel(trace1,ind[i]), spa_b);
     245      320784 :       if (t > spa_bs2) t = spa_b - t;
     246      320784 :       if (t > Sbound)
     247             :       {
     248      285497 :         if (DEBUGLEVEL>6) err_printf(".");
     249      285497 :         goto NEXT;
     250             :       }
     251             :       /* d - 2 test */
     252      112805 :       for (t=uel(trace2,ind[1]),i=2; i<=K; i++)
     253       77518 :         t = Fl_add(t, uel(trace2,ind[i]), spa_b);
     254       35287 :       if (t > spa_bs2) t = spa_b - t;
     255       35287 :       if (t > Sbound)
     256             :       {
     257       22134 :         if (DEBUGLEVEL>6) err_printf("|");
     258       22134 :         goto NEXT;
     259             :       }
     260             : 
     261       13153 :       av = avma;
     262             :       /* check trailing coeff */
     263       13153 :       y = lc;
     264       56392 :       for (i=1; i<=K; i++)
     265             :       {
     266       43239 :         GEN q = constant_coeff(gel(famod,ind[i]));
     267       43239 :         if (y) q = mulii(y, q);
     268       43239 :         y = centermodii(q, pa, pas2);
     269             :       }
     270       13153 :       if (!signe(y) || remii(constant_coeff(lcpol), y) != gen_0)
     271             :       {
     272       10745 :         if (DEBUGLEVEL>3) err_printf("T");
     273       10745 :         avma = av; goto NEXT;
     274             :       }
     275        2408 :       y = lc; /* full computation */
     276        5306 :       for (i=1; i<=K; i++)
     277             :       {
     278        2898 :         GEN q = gel(famod,ind[i]);
     279        2898 :         if (y) q = gmul(y, q);
     280        2898 :         y = centermod_i(q, pa, pas2);
     281             :       }
     282             : 
     283             :       /* y is the candidate factor */
     284        2408 :       if (! (q = ZX_divides_i(lcpol,y,bound)) )
     285             :       {
     286         252 :         if (DEBUGLEVEL>3) err_printf("*");
     287         252 :         avma = av; goto NEXT;
     288             :       }
     289             :       /* found a factor */
     290        2156 :       list = cgetg(K+1, t_VEC);
     291        2156 :       gel(listmod,cnt) = list;
     292        2156 :       for (i=1; i<=K; i++) list[i] = famod[ind[i]];
     293             : 
     294        2156 :       y = Q_primpart(y);
     295        2156 :       gel(fa,cnt++) = y;
     296             :       /* fix up pol */
     297        2156 :       pol = q;
     298        2156 :       if (lc) pol = Q_div_to_int(pol, leading_coeff(y));
     299       10234 :       for (i=j=k=1; i <= lfamod; i++)
     300             :       { /* remove used factors */
     301        8078 :         if (j <= K && i == ind[j]) j++;
     302             :         else
     303             :         {
     304        5628 :           gel(famod,k) = gel(famod,i);
     305        5628 :           uel(trace1,k) = uel(trace1,i);
     306        5628 :           uel(trace2,k) = uel(trace2,i);
     307        5628 :           deg[k] = deg[i]; k++;
     308             :         }
     309             :       }
     310        2156 :       lfamod -= K;
     311        2156 :       *pmaxK = cmbf_maxK(lfamod);
     312        2156 :       if (lfamod < 2*K) goto END;
     313        1169 :       i = 1; curdeg = deg[ind[1]];
     314        1169 :       bound = factor_bound(pol);
     315        1169 :       if (lc) lc = absi(leading_coeff(pol));
     316        1169 :       lcpol = lc? ZX_Z_mul(pol, lc): pol;
     317        1169 :       if (DEBUGLEVEL>3)
     318           0 :         err_printf("\nfound factor %Ps\nremaining modular factor(s): %ld\n",
     319             :                    y, lfamod);
     320        1169 :       continue;
     321             :     }
     322             : 
     323             : NEXT:
     324      318628 :     for (i = K+1;;)
     325             :     {
     326      367916 :       if (--i == 0) { K++; goto nextK; }
     327      362973 :       if (++ind[i] <= lfamod - K + i)
     328             :       {
     329      313685 :         curdeg = degsofar[i-1] + deg[ind[i]];
     330      313685 :         if (curdeg <= klim) break;
     331             :       }
     332       49288 :     }
     333      314854 :   }
     334             : END:
     335        4887 :   *done = 1;
     336        4887 :   if (degpol(pol) > 0)
     337             :   { /* leftover factor */
     338        4887 :     if (signe(leading_coeff(pol)) < 0) pol = ZX_neg(pol);
     339        4887 :     if (lfamod >= 2*K) *done = 0;
     340             : 
     341        4887 :     setlg(famod, lfamod+1);
     342        4887 :     gel(listmod,cnt) = leafcopy(famod);
     343        4887 :     gel(fa,cnt++) = pol;
     344             :   }
     345        4887 :   if (DEBUGLEVEL>6) err_printf("\n");
     346        4887 :   setlg(listmod, cnt);
     347        4887 :   setlg(fa, cnt); return mkvec2(fa, listmod);
     348             : }
     349             : 
     350             : /* recombination of modular factors: van Hoeij's algorithm */
     351             : 
     352             : /* Q in Z[X], return Q(2^n) */
     353             : static GEN
     354       75980 : shifteval(GEN Q, long n)
     355             : {
     356       75980 :   long i, l = lg(Q);
     357             :   GEN s;
     358             : 
     359       75980 :   if (!signe(Q)) return gen_0;
     360       75980 :   s = gel(Q,l-1);
     361       75980 :   for (i = l-2; i > 1; i--) s = addii(gel(Q,i), shifti(s, n));
     362       75980 :   return s;
     363             : }
     364             : 
     365             : /* return integer y such that all |a| <= y if P(a) = 0 */
     366             : static GEN
     367       44241 : root_bound(GEN P0)
     368             : {
     369       44241 :   GEN Q = leafcopy(P0), lP = absi(leading_coeff(Q)), x,y,z;
     370       44241 :   long k, d = degpol(Q);
     371             : 
     372             :   /* P0 = lP x^d + Q, deg Q < d */
     373       44241 :   Q = normalizepol_lg(Q, d+2);
     374       44241 :   for (k=lg(Q)-1; k>1; k--) gel(Q,k) = absi(gel(Q,k));
     375       44241 :   k = (long)(fujiwara_bound(P0));
     376       76071 :   for (  ; k >= 0; k--)
     377             :   {
     378       75980 :     pari_sp av = avma;
     379             :     /* y = 2^k; Q(y) >= lP y^d ? */
     380       75980 :     if (cmpii(shifteval(Q,k), shifti(lP, d*k)) >= 0) break;
     381       31830 :     avma = av;
     382             :   }
     383       44241 :   if (k < 0) k = 0;
     384       44241 :   x = int2n(k);
     385       44241 :   y = int2n(k+1);
     386      291743 :   for(k=0; ; k++)
     387             :   {
     388      291743 :     z = shifti(addii(x,y), -1);
     389      291743 :     if (equalii(x,z) || k > 5) break;
     390      247502 :     if (cmpii(poleval(Q,z), mulii(lP, powiu(z, d))) < 0)
     391      131958 :       y = z;
     392             :     else
     393      115544 :       x = z;
     394      247502 :   }
     395       44241 :   return y;
     396             : }
     397             : 
     398             : GEN
     399         175 : special_pivot(GEN x)
     400             : {
     401         175 :   GEN t, perm, H = ZM_hnfperm(x,NULL,&perm);
     402         175 :   long i,j, l = lg(H), h = lgcols(H);
     403        1855 :   for (i=1; i<h; i++)
     404             :   {
     405        1778 :     int fl = 0;
     406       10955 :     for (j=1; j<l; j++)
     407             :     {
     408        9275 :       t = gcoeff(H,i,j);
     409        9275 :       if (signe(t))
     410             :       {
     411        1813 :         if (!is_pm1(t) || fl) return NULL;
     412        1715 :         fl = 1;
     413             :       }
     414             :     }
     415             :   }
     416          77 :   return rowpermute(H, perm_inv(perm));
     417             : }
     418             : 
     419             : GEN
     420         266 : chk_factors_get(GEN lt, GEN famod, GEN c, GEN T, GEN N)
     421             : {
     422         266 :   long i = 1, j, l = lg(famod);
     423         266 :   GEN V = cgetg(l, t_VEC);
     424        7476 :   for (j = 1; j < l; j++)
     425        7210 :     if (signe(gel(c,j))) gel(V,i++) = gel(famod,j);
     426         266 :   if (lt && i > 1) gel(V,1) = RgX_Rg_mul(gel(V,1), lt);
     427         266 :   setlg(V, i);
     428         266 :   return T? FpXQXV_prod(V, T, N): FpXV_prod(V,N);
     429             : }
     430             : 
     431             : static GEN
     432         119 : chk_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)
     433             : {
     434             :   long i, r;
     435         119 :   GEN pol = P, list, piv, y, ltpol, lt, paov2;
     436             : 
     437         119 :   piv = special_pivot(M_L);
     438         119 :   if (!piv) return NULL;
     439          49 :   if (DEBUGLEVEL>7) err_printf("special_pivot output:\n%Ps\n",piv);
     440             : 
     441          49 :   r  = lg(piv)-1;
     442          49 :   list = cgetg(r+1, t_VEC);
     443          49 :   lt = absi(leading_coeff(pol));
     444          49 :   if (is_pm1(lt)) lt = NULL;
     445          49 :   ltpol = lt? ZX_Z_mul(pol, lt): pol;
     446          49 :   paov2 = shifti(pa,-1);
     447          49 :   for (i = 1;;)
     448             :   {
     449         140 :     if (DEBUGLEVEL) err_printf("LLL_cmbf: checking factor %ld\n",i);
     450         140 :     y = chk_factors_get(lt, famod, gel(piv,i), NULL, pa);
     451         140 :     y = FpX_center(y, pa, paov2);
     452         140 :     if (! (pol = ZX_divides_i(ltpol,y,bound)) ) return NULL;
     453         133 :     if (lt) y = Q_primpart(y);
     454         133 :     gel(list,i) = y;
     455         133 :     if (++i >= r) break;
     456             : 
     457          91 :     if (lt)
     458             :     {
     459          35 :       pol = ZX_Z_divexact(pol, leading_coeff(y));
     460          35 :       lt = absi(leading_coeff(pol));
     461          35 :       ltpol = ZX_Z_mul(pol, lt);
     462             :     }
     463             :     else
     464          56 :       ltpol = pol;
     465          91 :   }
     466          42 :   y = Q_primpart(pol);
     467          42 :   gel(list,i) = y; return list;
     468             : }
     469             : 
     470             : GEN
     471        1330 : LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, long *ti_LLL)
     472             : {
     473             :   GEN norm, u;
     474             :   long i, R;
     475             :   pari_timer T;
     476             : 
     477        1330 :   if (DEBUGLEVEL>2) timer_start(&T);
     478        1330 :   u = ZM_lll_norms(m, final? 0.999: 0.75, LLL_INPLACE, &norm);
     479        1330 :   if (DEBUGLEVEL>2) *ti_LLL += timer_delay(&T);
     480        8897 :   for (R=lg(m)-1; R > 0; R--)
     481        8897 :     if (cmprr(gel(norm,R), Bnorm) < 0) break;
     482        1330 :   for (i=1; i<=R; i++) setlg(u[i], n0+1);
     483        1330 :   if (R <= 1)
     484             :   {
     485         119 :     if (!R) pari_err_BUG("LLL_cmbf [no factor]");
     486         119 :     return NULL; /* irreducible */
     487             :   }
     488        1211 :   setlg(u, R+1); return u;
     489             : }
     490             : 
     491             : static ulong
     492          14 : next2pow(ulong a)
     493             : {
     494          14 :   ulong b = 1;
     495          14 :   while (b < a) b <<= 1;
     496          14 :   return b;
     497             : }
     498             : 
     499             : /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
     500             :  * van Hoeij's knapsack
     501             :  *
     502             :  * P = squarefree in Z[X].
     503             :  * famod = array of (lifted) modular factors mod p^a
     504             :  * bound = Mignotte bound for the size of divisors of P (for the sup norm)
     505             :  * previously recombined all set of factors with less than rec elts */
     506             : static GEN
     507          91 : LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
     508             : {
     509          91 :   const long N0 = 1; /* # of traces added at each step */
     510          91 :   double BitPerFactor = 0.4; /* nb bits in p^(a-b) / modular factor */
     511          91 :   long i,j,tmax,n0,C, dP = degpol(P);
     512          91 :   double logp = log((double)itos(p)), LOGp2 = LOG2/logp;
     513          91 :   double b0 = log((double)dP*2) / logp, logBr;
     514             :   GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
     515             :   pari_sp av, av2;
     516          91 :   long ti_LLL = 0, ti_CF  = 0;
     517             : 
     518          91 :   lP = absi(leading_coeff(P));
     519          91 :   if (is_pm1(lP)) lP = NULL;
     520          91 :   Br = root_bound(P);
     521          91 :   if (lP) Br = mulii(lP, Br);
     522          91 :   logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;
     523             : 
     524          91 :   n0 = lg(famod) - 1;
     525          91 :   C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */
     526          91 :   Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);
     527          91 :   ZERO = zeromat(n0, N0);
     528             : 
     529          91 :   av = avma;
     530          91 :   TT = cgetg(n0+1, t_VEC);
     531          91 :   Tra  = cgetg(n0+1, t_MAT);
     532        1792 :   for (i=1; i<=n0; i++)
     533             :   {
     534        1701 :     TT[i]  = 0;
     535        1701 :     gel(Tra,i) = cgetg(N0+1, t_COL);
     536             :   }
     537          91 :   CM_L = scalarmat_s(C, n0);
     538             :   /* tmax = current number of traces used (and computed so far) */
     539         441 :   for (tmax = 0;; tmax += N0)
     540             :   {
     541         441 :     long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;
     542             :     GEN M_L, q, CM_Lp, oldCM_L;
     543         441 :     int first = 1;
     544             :     pari_timer ti2, TI;
     545             : 
     546         441 :     bmin = (long)ceil(b0 + tnew*logBr);
     547         441 :     if (DEBUGLEVEL>2)
     548           0 :       err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
     549             :                  r, tmax, bmin);
     550             : 
     551             :     /* compute Newton sums (possibly relifting first) */
     552         441 :     if (a <= bmin)
     553             :     {
     554          14 :       a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */
     555          14 :       a = (long)next2pow((ulong)a);
     556             : 
     557          14 :       pa = powiu(p,a);
     558          14 :       famod = ZpX_liftfact(P, famod, pa, p, a);
     559          14 :       for (i=1; i<=n0; i++) TT[i] = 0;
     560             :     }
     561        8365 :     for (i=1; i<=n0; i++)
     562             :     {
     563        7924 :       GEN p1 = gel(Tra,i);
     564        7924 :       GEN p2 = polsym_gen(gel(famod,i), gel(TT,i), tnew, NULL, pa);
     565        7924 :       gel(TT,i) = p2;
     566        7924 :       p2 += 1+tmax; /* ignore traces number 0...tmax */
     567        7924 :       for (j=1; j<=N0; j++) gel(p1,j) = gel(p2,j);
     568        7924 :       if (lP)
     569             :       { /* make Newton sums integral */
     570        1848 :         GEN lPpow = powiu(lP, tmax);
     571        3696 :         for (j=1; j<=N0; j++)
     572             :         {
     573        1848 :           lPpow = mulii(lPpow,lP);
     574        1848 :           gel(p1,j) = mulii(gel(p1,j), lPpow);
     575             :         }
     576             :       }
     577             :     }
     578             : 
     579             :     /* compute truncation parameter */
     580         441 :     if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
     581         441 :     oldCM_L = CM_L;
     582         441 :     av2 = avma;
     583         441 :     delta = b = 0; /* -Wall */
     584             : AGAIN:
     585         840 :     M_L = Q_div_to_int(CM_L, utoipos(C));
     586         840 :     T2 = centermod( ZM_mul(Tra, M_L), pa );
     587         840 :     if (first)
     588             :     { /* initialize lattice, using few p-adic digits for traces */
     589         441 :       double t = gexpo(T2) - maxdd(32.0, BitPerFactor*r);
     590         441 :       bgood = (long) (t * LOGp2);
     591         441 :       b = maxss(bmin, bgood);
     592         441 :       delta = a - b;
     593             :     }
     594             :     else
     595             :     { /* add more p-adic digits and continue reduction */
     596         399 :       long b0 = (long)(gexpo(T2) * LOGp2);
     597         399 :       if (b0 < b) b = b0;
     598         399 :       b = maxss(b-delta, bmin);
     599         399 :       if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
     600             :     }
     601             : 
     602         840 :     q = powiu(p, b);
     603         840 :     m = vconcat( CM_L, gdivround(T2, q) );
     604         840 :     if (first)
     605             :     {
     606         441 :       GEN P1 = scalarmat(powiu(p, a-b), N0);
     607         441 :       first = 0;
     608         441 :       m = shallowconcat( m, vconcat(ZERO, P1) );
     609             :       /*     [ C M_L        0     ]
     610             :        * m = [                    ]   square matrix
     611             :        *     [  T2'  p^(a-b) I_N0 ]   T2' = Tra * M_L  truncated
     612             :        */
     613             :     }
     614             : 
     615         840 :     CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
     616         840 :     if (DEBUGLEVEL>2)
     617           0 :       err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
     618           0 :                  a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
     619         931 :     if (!CM_L) { list = mkvec(P); break; }
     620         791 :     if (b > bmin)
     621             :     {
     622         399 :       CM_L = gerepilecopy(av2, CM_L);
     623         399 :       goto AGAIN;
     624             :     }
     625         392 :     if (DEBUGLEVEL>2) timer_printf(&ti2, "for this block of traces");
     626             : 
     627         392 :     i = lg(CM_L) - 1;
     628         392 :     if (i == r && ZM_equal(CM_L, oldCM_L))
     629             :     {
     630          21 :       CM_L = oldCM_L;
     631          21 :       avma = av2; continue;
     632             :     }
     633             : 
     634         371 :     CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
     635         371 :     if (lg(CM_Lp) != lg(CM_L))
     636             :     {
     637           7 :       if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
     638           7 :       CM_L = ZM_hnf(CM_L);
     639             :     }
     640             : 
     641         371 :     if (i <= r && i*rec < n0)
     642             :     {
     643             :       pari_timer ti;
     644         119 :       if (DEBUGLEVEL>2) timer_start(&ti);
     645         119 :       list = chk_factors(P, Q_div_to_int(CM_L,utoipos(C)), bound, famod, pa);
     646         119 :       if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
     647         119 :       if (list) break;
     648          77 :       if (DEBUGLEVEL>2) err_printf("LLL_cmbf: chk_factors failed");
     649             :     }
     650         329 :     CM_L = gerepilecopy(av2, CM_L);
     651         329 :     if (gc_needed(av,1))
     652             :     {
     653           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"LLL_cmbf");
     654           0 :       gerepileall(av, 5, &CM_L, &TT, &Tra, &famod, &pa);
     655             :     }
     656         350 :   }
     657          91 :   if (DEBUGLEVEL>2)
     658           0 :     err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
     659          91 :   return list;
     660             : }
     661             : 
     662             : /* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */
     663             : static int
     664        4887 : cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)
     665             : {
     666        4887 :   long a,b,amin,d = (long)(31 * LOG2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);
     667        4887 :   int fl = 0;
     668             : 
     669        4887 :   b = logintall(B, q, qb) + 1;
     670        4887 :   *qb = mulii(*qb, q);
     671        4887 :   amin = b + d;
     672        4887 :   if (gcmp(powiu(q, amin), A) <= 0)
     673             :   {
     674        1337 :     a = logintall(A, q, qa) + 1;
     675        1337 :     *qa = mulii(*qa, q);
     676        1337 :     b = a - d; *qb = powiu(q, b);
     677             :   }
     678             :   else
     679             :   { /* not enough room */
     680        3550 :     a = amin;  *qa = powiu(q, a);
     681        3550 :     fl = 1;
     682             :   }
     683        4887 :   if (DEBUGLEVEL > 3) {
     684           0 :     err_printf("S_2   bound: %Ps^%ld\n", q,b);
     685           0 :     err_printf("coeff bound: %Ps^%ld\n", q,a);
     686             :   }
     687        4887 :   *pta = a;
     688        4887 :   *ptb = b; return fl;
     689             : }
     690             : 
     691             : /* use van Hoeij's knapsack algorithm */
     692             : static GEN
     693        4887 : combine_factors(GEN target, GEN famod, GEN p, long klim)
     694             : {
     695             :   GEN la, B, A, res, L, pa, pb, listmod;
     696        4887 :   long a,b, l, maxK, n = degpol(target);
     697             :   int done;
     698             :   pari_timer T;
     699             : 
     700        4887 :   A = factor_bound(target);
     701             : 
     702        4887 :   la = absi(leading_coeff(target));
     703        4887 :   B = mului(n, sqri(mulii(la, root_bound(target)))); /* = bound for S_2 */
     704             : 
     705        4887 :   (void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);
     706             : 
     707        4887 :   if (DEBUGLEVEL>2) timer_start(&T);
     708        4887 :   famod = ZpX_liftfact(target, famod, pa, p, a);
     709        4887 :   if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %Ps^%ld)", p,a);
     710        4887 :   L = cmbf(target, famod, A, p, a, b, klim, &maxK, &done);
     711        4887 :   if (DEBUGLEVEL>2) timer_printf(&T, "Naive recombination");
     712             : 
     713        4887 :   res     = gel(L,1);
     714        4887 :   listmod = gel(L,2); l = lg(listmod)-1;
     715        4887 :   famod = gel(listmod,l);
     716        4887 :   if (maxK > 0 && lg(famod)-1 > 2*maxK)
     717             :   {
     718          91 :     if (l!=1) A = factor_bound(gel(res,l));
     719          91 :     if (DEBUGLEVEL > 4) err_printf("last factor still to be checked\n");
     720          91 :     L = LLL_cmbf(gel(res,l), famod, p, pa, A, a, maxK);
     721          91 :     if (DEBUGLEVEL>2) timer_printf(&T,"Knapsack");
     722             :     /* remove last elt, possibly unfactored. Add all new ones. */
     723          91 :     setlg(res, l); res = shallowconcat(res, L);
     724             :   }
     725        4887 :   return res;
     726             : }
     727             : 
     728             : /* Assume 'a' a squarefree ZX; return 0 if no root (fl=1) / irreducible (fl=0).
     729             :  * Otherwise return prime p such that a mod p has fewest roots / factors */
     730             : static ulong
     731      597458 : pick_prime(GEN a, long fl, pari_timer *T)
     732             : {
     733      597458 :   pari_sp av = avma, av1;
     734      597458 :   const long MAXNP = 7, da = degpol(a);
     735      597458 :   long nmax = da+1, np;
     736      597458 :   ulong chosenp = 0;
     737      597458 :   GEN lead = gel(a,da+2);
     738             :   forprime_t S;
     739      597458 :   if (equali1(lead)) lead = NULL;
     740      597458 :   u_forprime_init(&S, 2, ULONG_MAX);
     741      597458 :   av1 = avma;
     742     3410608 :   for (np = 0; np < MAXNP; avma = av1)
     743             :   {
     744     3366472 :     ulong p = u_forprime_next(&S);
     745             :     long nfacp;
     746             :     GEN z;
     747             : 
     748     3366472 :     if (!p) pari_err_OVERFLOW("DDF [out of small primes]");
     749     3366472 :     if (lead && !umodiu(lead,p)) continue;
     750     3363966 :     z = ZX_to_Flx(a, p);
     751     3363966 :     if (!Flx_is_squarefree(z, p)) continue;
     752             : 
     753     2119866 :     if (fl)
     754             :     {
     755     2070089 :       nfacp = Flx_nbroots(z, p);
     756     2070089 :       if (!nfacp) { chosenp = 0; break; } /* no root */
     757             :     }
     758             :     else
     759             :     {
     760       49777 :       nfacp = Flx_nbfact(z, p);
     761       49777 :       if (nfacp == 1) { chosenp = 0; break; } /* irreducible */
     762             :     }
     763     1566558 :     if (DEBUGLEVEL>4)
     764           0 :       err_printf("...tried prime %3lu (%-3ld %s). Time = %ld\n",
     765             :                   p, nfacp, fl? "roots": "factors", timer_delay(T));
     766     1566558 :     if (nfacp < nmax)
     767             :     {
     768      531631 :       nmax = nfacp; chosenp = p;
     769      531631 :       if (da > 100 && nmax < 5) break; /* large degree, few factors. Enough */
     770             :     }
     771     1566544 :     np++;
     772             :   }
     773      597458 :   avma = av; return chosenp;
     774             : }
     775             : 
     776             : /* Assume A a squarefree ZX; return the vector of its rational roots */
     777             : static GEN
     778      581994 : DDF_roots(GEN A)
     779             : {
     780             :   GEN p, lc, lcpol, z, pe, pes2, bound;
     781             :   long i, m, e, lz;
     782             :   ulong pp;
     783             :   pari_sp av;
     784             :   pari_timer T;
     785             : 
     786      581994 :   if (DEBUGLEVEL>2) timer_start(&T);
     787      581994 :   pp = pick_prime(A, 1, &T);
     788      581994 :   if (!pp) return cgetg(1,t_VEC); /* no root */
     789       39263 :   p = utoipos(pp);
     790       39263 :   lc = leading_coeff(A);
     791       39263 :   if (is_pm1(lc))
     792       38332 :   { lc = NULL; lcpol = A; }
     793             :   else
     794         931 :   { lc = absi_shallow(lc); lcpol = ZX_Z_mul(A, lc); }
     795       39263 :   bound = root_bound(A); if (lc) bound = mulii(lc, bound);
     796       39263 :   e = logintall(addiu(shifti(bound, 1), 1), p, &pe) + 1;
     797       39263 :   pe = mulii(pe, p);
     798       39263 :   pes2 = shifti(pe, -1);
     799       39263 :   if (DEBUGLEVEL>2) timer_printf(&T, "Root bound");
     800       39263 :   av = avma;
     801       39263 :   z = ZpX_roots(A, p, e); lz = lg(z);
     802       39263 :   z = deg1_from_roots(z, varn(A));
     803       39263 :   if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %lu^%ld)", pp,e);
     804       79366 :   for (m=1, i=1; i < lz; i++)
     805             :   {
     806       40103 :     GEN q, r, y = gel(z,i);
     807       40103 :     if (lc) y = ZX_Z_mul(y, lc);
     808       40103 :     y = centermod_i(y, pe, pes2);
     809       40103 :     if (! (q = ZX_divides(lcpol, y)) ) continue;
     810             : 
     811        2268 :     lcpol = q;
     812        2268 :     r = negi( constant_coeff(y) );
     813        2268 :     if (lc) {
     814        1323 :       r = gdiv(r,lc);
     815        1323 :       lcpol = Q_primpart(lcpol);
     816        1323 :       lc = absi_shallow( leading_coeff(lcpol) );
     817        1323 :       if (is_pm1(lc)) lc = NULL; else lcpol = ZX_Z_mul(lcpol, lc);
     818             :     }
     819        2268 :     gel(z,m++) = r;
     820        2268 :     if (gc_needed(av,2))
     821             :     {
     822           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"DDF_roots, m = %ld", m);
     823           0 :       gerepileall(av, lc? 3:2, &z, &lcpol, &lc);
     824             :     }
     825             :   }
     826       39263 :   if (DEBUGLEVEL>2) timer_printf(&T, "Recombination");
     827       39263 :   z[0] = evaltyp(t_VEC) | evallg(m); return z;
     828             : }
     829             : 
     830             : /* Assume a squarefree ZX, deg(a) > 0, return rational factors.
     831             :  * In fact, a(0) != 0 but we don't use this */
     832             : static GEN
     833       15464 : DDF(GEN a)
     834             : {
     835             :   GEN ap, prime, famod, z;
     836       15464 :   long ti = 0;
     837       15464 :   ulong p = 0;
     838       15464 :   pari_sp av = avma;
     839             :   pari_timer T, T2;
     840             : 
     841       15464 :   if (DEBUGLEVEL>2) { timer_start(&T); timer_start(&T2); }
     842       15464 :   p = pick_prime(a, 0, &T2);
     843       15464 :   if (!p) return mkvec(a);
     844        4887 :   prime = utoipos(p);
     845        4887 :   ap = Flx_normalize(ZX_to_Flx(a, p), p);
     846        4887 :   famod = gel(Flx_factor(ap, p), 1);
     847        4887 :   if (DEBUGLEVEL>2)
     848             :   {
     849           0 :     if (DEBUGLEVEL>4) timer_printf(&T2, "splitting mod p = %lu", p);
     850           0 :     ti = timer_delay(&T);
     851           0 :     err_printf("Time setup: %ld\n", ti);
     852             :   }
     853        4887 :   z = combine_factors(a, FlxV_to_ZXV(famod), prime, degpol(a)-1);
     854        4887 :   if (DEBUGLEVEL>2)
     855           0 :     err_printf("Total Time: %ld\n===========\n", ti + timer_delay(&T));
     856        4887 :   return gerepilecopy(av, z);
     857             : }
     858             : 
     859             : /* Distinct Degree Factorization (deflating first)
     860             :  * Assume x squarefree, degree(x) > 0, x(0) != 0 */
     861             : GEN
     862       10935 : ZX_DDF(GEN x)
     863             : {
     864             :   GEN L;
     865             :   long m;
     866       10935 :   x = ZX_deflate_max(x, &m);
     867       10935 :   L = DDF(x);
     868       10935 :   if (m > 1)
     869             :   {
     870        3710 :     GEN e, v, fa = factoru(m);
     871             :     long i,j,k, l;
     872             : 
     873        3710 :     e = gel(fa,2); k = 0;
     874        3710 :     fa= gel(fa,1); l = lg(fa);
     875        3710 :     for (i=1; i<l; i++) k += e[i];
     876        3710 :     v = cgetg(k+1, t_VECSMALL); k = 1;
     877        7553 :     for (i=1; i<l; i++)
     878        3843 :       for (j=1; j<=e[i]; j++) v[k++] = fa[i];
     879        8162 :     for (k--; k; k--)
     880             :     {
     881        4452 :       GEN L2 = cgetg(1,t_VEC);
     882        8981 :       for (i=1; i < lg(L); i++)
     883        4529 :               L2 = shallowconcat(L2, DDF(RgX_inflate(gel(L,i), v[k])));
     884        4452 :       L = L2;
     885             :     }
     886             :   }
     887       10935 :   return L;
     888             : }
     889             : 
     890             : /* SquareFree Factorization. f = prod P^e, all e distinct, in Z[X] (char 0
     891             :  * would be enough, if ZX_gcd --> ggcd). Return (P), set *ex = (e) */
     892             : GEN
     893        8061 : ZX_squff(GEN f, GEN *ex)
     894             : {
     895             :   GEN T, V, P, e;
     896             :   long i, k, n, val;
     897             : 
     898        8061 :   if (signe(leading_coeff(f)) < 0) f = gneg_i(f);
     899        8061 :   val = ZX_valrem(f, &f);
     900        8061 :   n = 1 + degpol(f); if (val) n++;
     901        8061 :   e = cgetg(n,t_VECSMALL);
     902        8061 :   P = cgetg(n,t_COL);
     903             : 
     904        8061 :   T = ZX_gcd_all(f, ZX_deriv(f), &V);
     905        8404 :   for (k=i=1;; k++)
     906             :   {
     907        8404 :     pari_sp av = avma;
     908        8404 :     GEN W = ZX_gcd_all(T,V, &T);
     909        8404 :     long dW = degpol(W);
     910             :     /* W = prod P^e, e > k; V = prod P^e, e >= k */
     911        8404 :     if (dW == degpol(V)) /* V | T */
     912             :     {
     913             :       GEN U;
     914        1225 :       if (!dW) { avma = av; break; }
     915         217 :       while ( (U = ZX_divides(T, V)) ) { k++; T = U; }
     916         217 :       T = gerepileupto(av, T);
     917             :     }
     918             :     else
     919             :     {
     920        7179 :       gel(P,i) = Q_primpart(RgX_div(V,W));
     921        7179 :       e[i] = k; i++;
     922        7179 :       if (!dW) break;
     923         126 :       V = W;
     924             :     }
     925         343 :   }
     926        8061 :   if (val) { gel(P,i) = pol_x(varn(f)); e[i] = val; i++;}
     927        8061 :   setlg(P,i);
     928        8061 :   setlg(e,i); *ex = e; return P;
     929             : }
     930             : 
     931             : static GEN
     932        2773 : fact_from_DDF(GEN fa, GEN e, long n)
     933             : {
     934        2773 :   GEN v,w, y = cgetg(3, t_MAT);
     935        2773 :   long i,j,k, l = lg(fa);
     936             : 
     937        2773 :   v = cgetg(n+1,t_COL); gel(y,1) = v;
     938        2773 :   w = cgetg(n+1,t_COL); gel(y,2) = w;
     939        5812 :   for (k=i=1; i<l; i++)
     940             :   {
     941        3039 :     GEN L = gel(fa,i), ex = utoipos(e[i]);
     942        3039 :     long J = lg(L);
     943        7954 :     for (j=1; j<J; j++,k++)
     944             :     {
     945        4915 :       gel(v,k) = gcopy(gel(L,j));
     946        4915 :       gel(w,k) = ex;
     947             :     }
     948             :   }
     949        2773 :   return y;
     950             : }
     951             : 
     952             : /* Factor x in Z[t] */
     953             : static GEN
     954        2773 : ZX_factor_i(GEN x)
     955             : {
     956             :   GEN fa,ex,y;
     957             :   long n,i,l;
     958             : 
     959        2773 :   if (!signe(x)) return prime_fact(x);
     960        2773 :   fa = ZX_squff(x, &ex);
     961        2773 :   l = lg(fa); n = 0;
     962        5812 :   for (i=1; i<l; i++)
     963             :   {
     964        3039 :     gel(fa,i) = ZX_DDF(gel(fa,i));
     965        3039 :     n += lg(gel(fa,i))-1;
     966             :   }
     967        2773 :   y = fact_from_DDF(fa,ex,n);
     968        2773 :   return sort_factor_pol(y, cmpii);
     969             : }
     970             : GEN
     971        2458 : ZX_factor(GEN x)
     972             : {
     973        2458 :   pari_sp av = avma;
     974        2458 :   return gerepileupto(av, ZX_factor_i(x));
     975             : }
     976             : GEN
     977         315 : QX_factor(GEN x)
     978             : {
     979         315 :   pari_sp av = avma;
     980         315 :   return gerepileupto(av, ZX_factor_i(Q_primpart(x)));
     981             : }
     982             : 
     983             : long
     984        8897 : ZX_is_irred(GEN x)
     985             : {
     986        8897 :   pari_sp av = avma;
     987        8897 :   long l = lg(x);
     988             :   GEN y;
     989        8897 :   if (l <= 3) return 0; /* degree < 1 */
     990        8897 :   if (l == 4) return 1; /* degree 1 */
     991        7637 :   if (ZX_val(x)) return 0;
     992        7434 :   if (!ZX_is_squarefree(x)) return 0;
     993        7385 :   y = ZX_DDF(x); avma = av;
     994        7385 :   return (lg(y) == 2);
     995             : }
     996             : 
     997             : GEN
     998      581994 : nfrootsQ(GEN x)
     999             : {
    1000      581994 :   pari_sp av = avma;
    1001             :   GEN z;
    1002             :   long val;
    1003             : 
    1004      581994 :   if (typ(x)!=t_POL) pari_err_TYPE("nfrootsQ",x);
    1005      581994 :   if (!signe(x)) pari_err_ROOTS0("nfrootsQ");
    1006      581994 :   x = Q_primpart(x);
    1007      581994 :   RgX_check_ZX(x,"nfrootsQ");
    1008      581994 :   val = ZX_valrem(x, &x);
    1009      581994 :   (void)ZX_gcd_all(x, ZX_deriv(x), &x);
    1010      581994 :   z = DDF_roots(x);
    1011      581994 :   if (val) z = shallowconcat(z, gen_0);
    1012      581994 :   return gerepileupto(av, sort(z));
    1013             : }
    1014             : 
    1015             : /************************************************************************
    1016             :  *                   GCD OVER Z[X] / Q[X]                               *
    1017             :  ************************************************************************/
    1018             : int
    1019       16793 : ZX_is_squarefree(GEN x)
    1020             : {
    1021       16793 :   pari_sp av = avma;
    1022             :   GEN d;
    1023             :   long m;
    1024             :   int r;
    1025       16793 :   if (lg(x) == 2) return 0;
    1026       16793 :   m = ZX_deflate_order(x);
    1027       16793 :   if (m > 1)
    1028             :   {
    1029        5089 :     if (!signe(gel(x,2))) return 0;
    1030        5026 :     x = RgX_deflate(x, m);
    1031             :   }
    1032       16730 :   d = ZX_gcd(x,ZX_deriv(x));
    1033       16730 :   r = (lg(d) == 3); avma = av; return r;
    1034             : }
    1035             : 
    1036             : #if 0
    1037             : /* ceil( || p ||_oo / lc(p) ) */
    1038             : static GEN
    1039             : maxnorm(GEN p)
    1040             : {
    1041             :   long i, n = degpol(p), av = avma;
    1042             :   GEN x, m = gen_0;
    1043             : 
    1044             :   p += 2;
    1045             :   for (i=0; i<n; i++)
    1046             :   {
    1047             :     x = gel(p,i);
    1048             :     if (abscmpii(x,m) > 0) m = x;
    1049             :   }
    1050             :   m = divii(m, gel(p,n));
    1051             :   return gerepileuptoint(av, addis(absi(m),1));
    1052             : }
    1053             : #endif
    1054             : 
    1055             : /* A, B in Z[X] */
    1056             : GEN
    1057     4098478 : ZX_gcd_all(GEN A, GEN B, GEN *Anew)
    1058             : {
    1059             :   GEN R, a, b, q, H, Hp, g, Ag, Bg;
    1060     4098478 :   long m, n, valX, valA, vA = varn(A);
    1061             :   ulong p;
    1062             :   int small;
    1063             :   pari_sp ltop, av;
    1064             :   forprime_t S;
    1065             : 
    1066     4098478 :   if (!signe(A)) { if (Anew) *Anew = pol_0(vA); return ZX_copy(B); }
    1067     4098478 :   if (!signe(B)) { if (Anew) *Anew = pol_1(vA); return ZX_copy(A); }
    1068     4097386 :   valA = ZX_valrem(A, &A);
    1069     4097386 :   valX = minss(valA, ZX_valrem(B, &B));
    1070     4097386 :   ltop = avma;
    1071             : 
    1072     4097386 :   n = 1 + minss(degpol(A), degpol(B)); /* > degree(gcd) */
    1073     4097386 :   g = gcdii(leading_coeff(A), leading_coeff(B)); /* multiple of lead(gcd) */
    1074     4097386 :   if (is_pm1(g)) {
    1075     4090881 :     g = NULL;
    1076     4090881 :     Ag = A;
    1077     4090881 :     Bg = B;
    1078             :   } else {
    1079        6505 :     Ag = ZX_Z_mul(A,g);
    1080        6505 :     Bg = ZX_Z_mul(B,g);
    1081             :   }
    1082     4097386 :   small = (ZX_max_lg(A) == 3 && ZX_max_lg(B) == 3);
    1083     4097386 :   init_modular_big(&S);
    1084     4097386 :   av = avma;
    1085     4097386 :   R = NULL;/*-Wall*/
    1086     4097386 :   H = NULL;
    1087     8195222 :   while ((p = u_forprime_next(&S)))
    1088             :   {
    1089     4097836 :     if (g && !umodiu(g,p)) continue;
    1090     4097836 :     a = ZX_to_Flx(A, p);
    1091     4097836 :     b = ZX_to_Flx(B, p); Hp = Flx_gcd(a,b, p);
    1092     4097836 :     m = degpol(Hp);
    1093     4097836 :     if (m == 0) { /* coprime. DONE */
    1094      986829 :       avma = ltop;
    1095      986829 :       if (Anew) {
    1096      604137 :         if (valA != valX) A = RgX_shift(A, valA - valX);
    1097      604137 :         *Anew = A;
    1098             :       }
    1099      986829 :       return monomial(gen_1, valX, vA);
    1100             :     }
    1101     3111007 :     if (m > n) continue; /* p | Res(A/G, B/G). Discard */
    1102             : 
    1103     3111007 :     if (!g) /* make sure lead(H) = g mod p */
    1104     3110036 :       Hp = Flx_normalize(Hp, p);
    1105             :     else
    1106             :     {
    1107         971 :       ulong t = Fl_mul(umodiu(g, p), Fl_inv(Hp[m+2],p), p);
    1108         971 :       Hp = Flx_Fl_mul(Hp, t, p);
    1109             :     }
    1110     3111007 :     if (DEBUGLEVEL>5) err_printf("gcd mod %lu (bound 2^%ld)\n", p,expi(q));
    1111     3111007 :     if (m < n)
    1112             :     { /* First time or degree drop [all previous p were as above; restart]. */
    1113     3110557 :       H = ZX_init_CRT(Hp,p,vA);
    1114     3110557 :       q = utoipos(p); n = m;
    1115     3110557 :       if (!small) continue; /* if gcd is small, try our luck */
    1116             :     }
    1117             :     else
    1118         450 :       if (!ZX_incremental_CRT(&H, Hp, &q, p)) continue;
    1119             :     /* H stable: check divisibility */
    1120     3110557 :     if (!ZX_divides(Bg, H)) continue;
    1121     3110557 :     R = ZX_divides(Ag, H);
    1122     3110557 :     if (R) break;
    1123             : 
    1124           0 :     if (gc_needed(av,1))
    1125             :     {
    1126           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"QX_gcd");
    1127           0 :       gerepileall(av, 3, &H, &q, &Hp);
    1128             :     }
    1129             :   }
    1130     3110557 :   if (!p) pari_err_OVERFLOW("ZX_gcd_all [ran out of primes]");
    1131     3110557 :   if (Anew) {
    1132       11306 :     A = R;
    1133       11306 :     if (valA != valX) A = RgX_shift(A, valA - valX);
    1134       11306 :     *Anew = A;
    1135             :   }
    1136     3110557 :   return valX ? RgX_shift(H, valX): H;
    1137             : }
    1138             : GEN
    1139     3481950 : ZX_gcd(GEN A, GEN B) { return ZX_gcd_all(A,B,NULL); }
    1140             : 
    1141             : static GEN
    1142     3453444 : _gcd(GEN a, GEN b)
    1143             : {
    1144     3453444 :   if (!a) a = gen_1;
    1145     3453444 :   if (!b) b = gen_1;
    1146     3453444 :   return Q_gcd(a,b);
    1147             : }
    1148             : /* A0 and B0 in Q[X] */
    1149             : GEN
    1150     3453444 : QX_gcd(GEN A0, GEN B0)
    1151             : {
    1152             :   GEN a, b, D;
    1153     3453444 :   pari_sp av = avma, av2;
    1154             : 
    1155     3453444 :   D = ZX_gcd(Q_primitive_part(A0, &a), Q_primitive_part(B0, &b));
    1156     3453444 :   av2 = avma; a = _gcd(a,b);
    1157     3453444 :   if (isint1(a)) avma = av2; else D = RgX_Rg_mul(D, a);
    1158     3453444 :   return gerepileupto(av, D);
    1159             : }
    1160             : 
    1161             : /*****************************************************************************
    1162             :  * Variants of the Bradford-Davenport algorithm: look for cyclotomic         *
    1163             :  * factors, and decide whether a ZX is cyclotomic or a product of cyclotomic *
    1164             :  *****************************************************************************/
    1165             : /* f of degree 1, return a cyclotomic factor (Phi_1 or Phi_2) or NULL */
    1166             : static GEN
    1167           0 : BD_deg1(GEN f)
    1168             : {
    1169           0 :   GEN a = gel(f,3), b = gel(f,2); /* f = ax + b */
    1170           0 :   if (!absequalii(a,b)) return NULL;
    1171           0 :   return polcyclo((signe(a) == signe(b))? 2: 1, varn(f));
    1172             : }
    1173             : 
    1174             : /* f a squarefree ZX; not divisible by any Phi_n, n even */
    1175             : static GEN
    1176         406 : BD_odd(GEN f)
    1177             : {
    1178         812 :   while(degpol(f) > 1)
    1179             :   {
    1180         406 :     GEN f1 = ZX_graeffe(f); /* contain all cyclotomic divisors of f */
    1181         406 :     if (ZX_equal(f1, f)) return f; /* product of cyclotomics */
    1182           0 :     f = ZX_gcd(f, f1);
    1183             :   }
    1184           0 :   if (degpol(f) == 1) return BD_deg1(f);
    1185           0 :   return NULL; /* no cyclotomic divisor */
    1186             : }
    1187             : 
    1188             : static GEN
    1189        2310 : myconcat(GEN v, GEN x)
    1190             : {
    1191        2310 :   if (typ(x) != t_VEC) x = mkvec(x);
    1192        2310 :   if (!v) return x;
    1193        1470 :   return shallowconcat(v, x);
    1194             : }
    1195             : 
    1196             : /* Bradford-Davenport algorithm.
    1197             :  * f a squarefree ZX of degree > 0, return NULL or a vector of coprime
    1198             :  * cyclotomic factors of f [ possibly reducible ] */
    1199             : static GEN
    1200        2359 : BD(GEN f)
    1201             : {
    1202        2359 :   GEN G = NULL, Gs = NULL, Gp = NULL, Gi = NULL;
    1203             :   GEN fs2, fp, f2, f1, fe, fo, fe1, fo1;
    1204        2359 :   RgX_even_odd(f, &fe, &fo);
    1205        2359 :   fe1 = ZX_eval1(fe);
    1206        2359 :   fo1 = ZX_eval1(fo);
    1207        2359 :   if (absequalii(fe1, fo1)) /* f(1) = 0 or f(-1) = 0 */
    1208             :   {
    1209        1519 :     long i, v = varn(f);
    1210        1519 :     if (!signe(fe1))
    1211         371 :       G = mkvec2(polcyclo(1, v), polcyclo(2, v)); /* both 0 */
    1212        1148 :     else if (signe(fe1) == signe(fo1))
    1213         693 :       G = mkvec(polcyclo(2, v)); /*f(-1) = 0*/
    1214             :     else
    1215         455 :       G = mkvec(polcyclo(1, v)); /*f(1) = 0*/
    1216        1519 :     for (i = lg(G)-1; i; i--) f = RgX_div(f, gel(G,i));
    1217             :   }
    1218             :   /* f no longer divisible by Phi_1 or Phi_2 */
    1219        2359 :   if (degpol(f) <= 1) return G;
    1220        2058 :   f1 = ZX_graeffe(f); /* has at most square factors */
    1221        2058 :   if (ZX_equal(f1, f)) return myconcat(G,f); /* f = product of Phi_n, n odd */
    1222             : 
    1223        1183 :   fs2 = ZX_gcd_all(f1, ZX_deriv(f1), &f2); /* fs2 squarefree */
    1224        1183 :   if (degpol(fs2))
    1225             :   { /* fs contains all Phi_n | f, 4 | n; and only those */
    1226             :     /* In that case, Graeffe(Phi_n) = Phi_{n/2}^2, and Phi_n = Phi_{n/2}(x^2) */
    1227        1029 :     GEN fs = RgX_inflate(fs2, 2);
    1228        1029 :     (void)ZX_gcd_all(f, fs, &f); /* remove those Phi_n | f, 4 | n */
    1229        1029 :     Gs = BD(fs2);
    1230        1029 :     if (Gs)
    1231             :     {
    1232             :       long i;
    1233        1029 :       for (i = lg(Gs)-1; i; i--) gel(Gs,i) = RgX_inflate(gel(Gs,i), 2);
    1234             :       /* prod Gs[i] is the product of all Phi_n | f, 4 | n */
    1235        1029 :       G = myconcat(G, Gs);
    1236             :     }
    1237             :     /* f2 = f1 / fs2 */
    1238        1029 :     f1 = RgX_div(f2, fs2); /* f1 / fs2^2 */
    1239             :   }
    1240        1183 :   fp = ZX_gcd(f, f1); /* contains all Phi_n | f, n > 1 odd; and only those */
    1241        1183 :   if (degpol(fp))
    1242             :   {
    1243         196 :     Gp = BD_odd(fp);
    1244             :     /* Gp is the product of all Phi_n | f, n odd */
    1245         196 :     if (Gp) G = myconcat(G, Gp);
    1246         196 :     f = RgX_div(f, fp);
    1247             :   }
    1248        1183 :   if (degpol(f))
    1249             :   { /* contains all Phi_n originally dividing f, n = 2 mod 4, n > 2;
    1250             :      * and only those
    1251             :      * In that case, Graeffe(Phi_n) = Phi_{n/2}, and Phi_n = Phi_{n/2}(-x) */
    1252         210 :     Gi = BD_odd(ZX_unscale(f, gen_m1));
    1253         210 :     if (Gi)
    1254             :     { /* N.B. Phi_2 does not divide f */
    1255         210 :       Gi = ZX_unscale(Gi, gen_m1);
    1256             :       /* Gi is the product of all Phi_n | f, n = 2 mod 4 */
    1257         210 :       G = myconcat(G, Gi);
    1258             :     }
    1259             :   }
    1260        1183 :   return G;
    1261             : }
    1262             : 
    1263             : /* Let f be a non-zero QX, return the (squarefree) product of cyclotomic
    1264             :  * divisors of f */
    1265             : GEN
    1266         315 : polcyclofactors(GEN f)
    1267             : {
    1268         315 :   pari_sp av = avma;
    1269         315 :   if (typ(f) != t_POL || !signe(f)) pari_err_TYPE("polcyclofactors",f);
    1270         315 :   (void)RgX_valrem(f, &f);
    1271         315 :   f = Q_primpart(f);
    1272         315 :   RgX_check_ZX(f,"polcyclofactors");
    1273         315 :   if (degpol(f))
    1274             :   {
    1275         315 :     (void)ZX_gcd_all(f, ZX_deriv(f), &f);
    1276         315 :     f = BD(f);
    1277         315 :     if (f) return gerepilecopy(av, f);
    1278             :   }
    1279           0 :   avma = av; return cgetg(1,t_VEC);
    1280             : }
    1281             : 
    1282             : /* return t*x mod T(x), T a monic ZX. Assume deg(t) < deg(T) */
    1283             : static GEN
    1284       46018 : ZXQ_mul_by_X(GEN t, GEN T)
    1285             : {
    1286             :   GEN lt;
    1287       46018 :   t = RgX_shift_shallow(t, 1);
    1288       46018 :   if (degpol(t) < degpol(T)) return t;
    1289        3822 :   lt = leading_coeff(t);
    1290        3822 :   if (is_pm1(lt)) return signe(lt) > 0 ? ZX_sub(t, T): ZX_add(t, T);
    1291         217 :   return ZX_sub(t, ZX_Z_mul(T, leading_coeff(t)));
    1292             : }
    1293             : /* f a product of Phi_n, all n odd; deg f > 1. Is it irreducible ? */
    1294             : static long
    1295         840 : BD_odd_iscyclo(GEN f)
    1296             : {
    1297             :   pari_sp av;
    1298             :   long d, e, n, bound;
    1299             :   GEN t;
    1300         840 :   f = ZX_deflate_max(f, &e);
    1301         840 :   av = avma;
    1302             :   /* The original f is cyclotomic (= Phi_{ne}) iff the present one is Phi_n,
    1303             :    * where all prime dividing e also divide n. If current f is Phi_n,
    1304             :    * then n is odd and squarefree */
    1305         840 :   d = degpol(f); /* = phi(n) */
    1306             :   /* Let e > 0, g multiplicative such that
    1307             :        g(p) = p / (p-1)^(1+e) < 1 iff p < (p-1)^(1+e)
    1308             :      For all squarefree odd n, we have g(n) < C, hence n < C phi(n)^(1+e), where
    1309             :        C = \prod_{p odd | p > (p-1)^(1+e)} g(p)
    1310             :      For e = 1/10,   we obtain p = 3, 5 and C < 1.523
    1311             :      For e = 1/100,  we obtain p = 3, 5, ..., 29 and C < 2.573
    1312             :      In fact, for n <= 10^7 odd & squarefree, we have n < 2.92 * phi(n)
    1313             :      By the above, n<10^7 covers all d <= (10^7/2.573)^(1/(1+1/100)) < 3344391.
    1314             :   */
    1315         840 :   if (d <= 3344391)
    1316         840 :     bound = (long)(2.92 * d);
    1317             :   else
    1318           0 :     bound = (long)(2.573 * pow(d,1.01));
    1319             :   /* IF f = Phi_n, n squarefree odd, then n <= bound */
    1320         840 :   t = monomial(gen_1, d-1, varn(f));
    1321       46053 :   for (n = d; n <= bound; n++)
    1322             :   {
    1323       46018 :     t = ZXQ_mul_by_X(t, f);
    1324             :     /* t = (X mod f(X))^d */
    1325       46018 :     if (degpol(t) == 0) break;
    1326       45213 :     if (gc_needed(av,1))
    1327             :     {
    1328         454 :       if(DEBUGMEM>1) pari_warn(warnmem,"BD_odd_iscyclo");
    1329         454 :       t = gerepilecopy(av, t);
    1330             :     }
    1331             :   }
    1332         840 :   if (n > bound || eulerphiu(n) != (ulong)d) return 0;
    1333             : 
    1334         777 :   if (e > 1) return (ucoprime_part(e, n) == 1)? e * n : 0;
    1335         651 :   return n;
    1336             : }
    1337             : 
    1338             : /* Checks if f, monic squarefree ZX with |constant coeff| = 1, is a cyclotomic
    1339             :  * polynomial. Returns n if f = Phi_n, and 0 otherwise */
    1340             : static long
    1341        2373 : BD_iscyclo(GEN f)
    1342             : {
    1343        2373 :   pari_sp av = avma;
    1344             :   GEN f2, fn, f1;
    1345             : 
    1346        2373 :   if (degpol(f) == 1) return isint1(gel(f,2))? 2: 1;
    1347        2240 :   f1 = ZX_graeffe(f);
    1348             :   /* f = product of Phi_n, n odd */
    1349        2240 :   if (ZX_equal(f, f1)) { avma = av; return BD_odd_iscyclo(f); }
    1350             : 
    1351        1778 :   fn = ZX_unscale(f, gen_m1); /* f(-x) */
    1352             :   /* f = product of Phi_n, n = 2 mod 4 */
    1353        1778 :   if (ZX_equal(f1, fn)) return 2*BD_odd_iscyclo(fn);
    1354             : 
    1355        1400 :   if (issquareall(f1, &f2))
    1356             :   {
    1357         595 :     GEN lt = leading_coeff(f2);
    1358             :     long c;
    1359         595 :     if (signe(lt) < 0) f2 = ZX_neg(f2);
    1360         595 :     c = BD_iscyclo(f2);
    1361         595 :     return odd(c)? 0: 2*c;
    1362             :   }
    1363         805 :   avma = av; return 0;
    1364             : }
    1365             : long
    1366        3647 : poliscyclo(GEN f)
    1367             : {
    1368             :   long d;
    1369        3647 :   if (typ(f) != t_POL) pari_err_TYPE("poliscyclo", f);
    1370        3640 :   d = degpol(f);
    1371        3640 :   if (d <= 0 || !RgX_is_ZX(f)) return 0;
    1372        3633 :   if (!equali1(gel(f,d+2)) || !is_pm1(gel(f,2))) return 0;
    1373        1813 :   if (d == 1) return signe(gel(f,2)) > 0? 2: 1;
    1374        1778 :   return ZX_is_squarefree(f)? BD_iscyclo(f): 0;
    1375             : }
    1376             : 
    1377             : long
    1378        1029 : poliscycloprod(GEN f)
    1379             : {
    1380        1029 :   pari_sp av = avma;
    1381        1029 :   long i, d = degpol(f);
    1382        1029 :   if (typ(f) != t_POL) pari_err_TYPE("poliscycloprod",f);
    1383        1029 :   if (!RgX_is_ZX(f)) return 0;
    1384        1029 :   if (!equali1(leading_coeff(f)) || !is_pm1(constant_coeff(f))) return 0;
    1385        1029 :   if (d < 2) return (d == 1);
    1386        1022 :   if ( degpol(ZX_gcd_all(f, ZX_deriv(f), &f)) )
    1387             :   {
    1388          14 :     d = degpol(f);
    1389          14 :     if (d == 1) return 1;
    1390             :   }
    1391        1015 :   f = BD(f); if (!f) return 0;
    1392        1015 :   for (i = lg(f)-1; i; i--) d -= degpol(gel(f,i));
    1393        1015 :   avma = av; return d == 0;
    1394             : }

Generated by: LCOV version 1.11