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 - QX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16912-212c0f0) Lines: 709 765 92.7 %
Date: 2014-10-20 Functions: 40 43 93.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 455 588 77.4 %

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

Generated by: LCOV version 1.9