Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - QX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29796-c9bd12ad19) Lines: 773 801 96.5 %
Date: 2024-12-09 09:09:46 Functions: 44 45 97.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_factor
      18             : 
      19             : /* x,y two ZX, y non constant. Return q = x/y if y divides x in Z[X] and NULL
      20             :  * otherwise. If not NULL, B is a t_INT upper bound for ||q||_oo. */
      21             : static GEN
      22     6569998 : ZX_divides_i(GEN x, GEN y, GEN B)
      23             : {
      24             :   long dx, dy, dz, i, j;
      25             :   pari_sp av;
      26             :   GEN z,p1,y_lead;
      27             : 
      28     6569998 :   dy=degpol(y);
      29     6569999 :   dx=degpol(x);
      30     6570000 :   dz=dx-dy; if (dz<0) return NULL;
      31     6568936 :   z=cgetg(dz+3,t_POL); z[1] = x[1];
      32     6568936 :   x += 2; y += 2; z += 2;
      33     6568936 :   y_lead = gel(y,dy);
      34     6568936 :   if (equali1(y_lead)) y_lead = NULL;
      35             : 
      36     6568936 :   p1 = gel(x,dx);
      37     6568936 :   if (y_lead) {
      38             :     GEN r;
      39       30135 :     p1 = dvmdii(p1,y_lead, &r);
      40       30135 :     if (r != gen_0) return NULL;
      41             :   }
      42     6538801 :   else p1 = icopy(p1);
      43     6564775 :   gel(z,dz) = p1;
      44     8228016 :   for (i=dx-1; i>=dy; i--)
      45             :   {
      46     1668417 :     av = avma; p1 = gel(x,i);
      47     5790512 :     for (j=i-dy+1; j<=i && j<=dz; j++)
      48     4122124 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
      49     1668388 :     if (y_lead) {
      50             :       GEN r;
      51       52822 :       p1 = dvmdii(p1,y_lead, &r);
      52       52822 :       if (r != gen_0) return NULL;
      53             :     }
      54     1666889 :     if (B && abscmpii(p1, B) > 0) return NULL;
      55     1663214 :     p1 = gerepileuptoint(av, p1);
      56     1663241 :     gel(z,i-dy) = p1;
      57             :   }
      58     6559599 :   av = avma;
      59    17095563 :   for (; i >= 0; i--)
      60             :   {
      61    10579116 :     p1 = gel(x,i);
      62             :     /* we always enter this loop at least once */
      63    23806364 :     for (j=0; j<=i && j<=dz; j++)
      64    13227259 :       p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
      65    10579105 :     if (signe(p1)) return NULL;
      66    10535963 :     set_avma(av);
      67             :   }
      68     6516447 :   return z - 2;
      69             : }
      70             : static GEN
      71     6541347 : ZX_divides(GEN x, GEN y) { return ZX_divides_i(x,y,NULL); }
      72             : 
      73             : #if 0
      74             : /* cf Beauzamy et al: upper bound for
      75             :  *      lc(x) * [2^(5/8) / pi^(3/8)] e^(1/4n) 2^(n/2) sqrt([x]_2)/ n^(3/8)
      76             :  * where [x]_2 = sqrt(\sum_i=0^n x[i]^2 / binomial(n,i)). One factor has
      77             :  * all coeffs less than then bound */
      78             : static GEN
      79             : two_factor_bound(GEN x)
      80             : {
      81             :   long i, j, n = lg(x) - 3;
      82             :   pari_sp av = avma;
      83             :   GEN *invbin, c, r = cgetr(LOWDEFAULTPREC), z;
      84             : 
      85             :   x += 2; invbin = (GEN*)new_chunk(n+1);
      86             :   z = real_1(LOWDEFAULTPREC); /* invbin[i] = 1 / binomial(n, i) */
      87             :   for (i=0,j=n; j >= i; i++,j--)
      88             :   {
      89             :     invbin[i] = invbin[j] = z;
      90             :     z = divru(mulru(z, i+1), n-i);
      91             :   }
      92             :   z = invbin[0]; /* = 1 */
      93             :   for (i=0; i<=n; i++)
      94             :   {
      95             :     c = gel(x,i); if (!signe(c)) continue;
      96             :     affir(c, r);
      97             :     z = addrr(z, mulrr(sqrr(r), invbin[i]));
      98             :   }
      99             :   z = shiftr(sqrtr(z), n);
     100             :   z = divrr(z, dbltor(pow((double)n, 0.75)));
     101             :   z = roundr_safe(sqrtr(z));
     102             :   z = mulii(z, absi_shallow(gel(x,n)));
     103             :   return gerepileuptoint(av, shifti(z, 1));
     104             : }
     105             : #endif
     106             : 
     107             : /* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */
     108             : static GEN
     109       60134 : Mignotte_bound(GEN S)
     110             : {
     111       60134 :   long i, d = degpol(S);
     112       60134 :   GEN C, N2, t, binlS, lS = leading_coeff(S), bin = vecbinomial(d-1);
     113             : 
     114       60134 :   N2 = sqrtr(RgX_fpnorml2(S,DEFAULTPREC));
     115       60133 :   binlS = is_pm1(lS)? bin: ZC_Z_mul(bin, lS);
     116             : 
     117             :   /* i = 0 */
     118       60133 :   C = gel(binlS,1);
     119             :   /* i = d */
     120       60133 :   t = N2; if (gcmp(C, t) < 0) C = t;
     121      520512 :   for (i = 1; i < d; i++)
     122             :   {
     123      460379 :     t = addri(mulir(gel(bin,i), N2), gel(binlS,i+1));
     124      460379 :     if (mpcmp(C, t) < 0) C = t;
     125             :   }
     126       60133 :   return C;
     127             : }
     128             : /* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,
     129             :  * where [P]_2 is Bombieri's 2-norm */
     130             : static GEN
     131       60133 : Beauzamy_bound(GEN S)
     132             : {
     133       60133 :   const long prec = DEFAULTPREC;
     134       60133 :   long i, d = degpol(S);
     135             :   GEN bin, lS, s, C;
     136       60133 :   bin = vecbinomial(d);
     137             : 
     138       60133 :   s = real_0(prec);
     139      640782 :   for (i=0; i<=d; i++)
     140             :   {
     141      580648 :     GEN c = gel(S,i+2);
     142      580648 :     if (gequal0(c)) continue;
     143             :     /* s += P_i^2 / binomial(d,i) */
     144      488619 :     s = addrr(s, divri(itor(sqri(c), prec), gel(bin,i+1)));
     145             :   }
     146             :   /* s = [S]_2^2 */
     147       60134 :   C = powruhalf(utor(3,prec), 3 + 2*d); /* 3^{3/2 + d} */
     148       60134 :   C = divrr(mulrr(C, s), mulur(4*d, mppi(prec)));
     149       60134 :   lS = absi_shallow(leading_coeff(S));
     150       60134 :   return mulir(lS, sqrtr(C));
     151             : }
     152             : 
     153             : static GEN
     154       60134 : factor_bound(GEN S)
     155             : {
     156       60134 :   pari_sp av = avma;
     157       60134 :   GEN a = Mignotte_bound(S);
     158       60133 :   GEN b = Beauzamy_bound(S);
     159       60134 :   if (DEBUGLEVEL>2)
     160             :   {
     161           0 :     err_printf("Mignotte bound: %Ps\n",a);
     162           0 :     err_printf("Beauzamy bound: %Ps\n",b);
     163             :   }
     164       60134 :   return gerepileupto(av, ceil_safe(gmin_shallow(a, b)));
     165             : }
     166             : 
     167             : /* Naive recombination of modular factors: combine up to maxK modular
     168             :  * factors, degree <= klim
     169             :  *
     170             :  * target = polynomial we want to factor
     171             :  * famod = array of modular factors.  Product should be congruent to
     172             :  * target/lc(target) modulo p^a
     173             :  * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
     174             : static GEN
     175       49623 : cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,
     176             :      long klim, long *pmaxK, int *done)
     177             : {
     178       49623 :   long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
     179             :   ulong spa_b, spa_bs2, Sbound;
     180       49623 :   GEN lc, lcpol, pa = powiu(p,a), pas2 = shifti(pa,-1);
     181       49623 :   GEN trace1   = cgetg(lfamod+1, t_VECSMALL);
     182       49623 :   GEN trace2   = cgetg(lfamod+1, t_VECSMALL);
     183       49623 :   GEN ind      = cgetg(lfamod+1, t_VECSMALL);
     184       49623 :   GEN deg      = cgetg(lfamod+1, t_VECSMALL);
     185       49623 :   GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
     186       49623 :   GEN listmod  = cgetg(lfamod+1, t_VEC);
     187       49623 :   GEN fa       = cgetg(lfamod+1, t_VEC);
     188             : 
     189       49623 :   *pmaxK = cmbf_maxK(lfamod);
     190       49623 :   lc = absi_shallow(leading_coeff(pol));
     191       49623 :   if (equali1(lc)) lc = NULL;
     192       49623 :   lcpol = lc? ZX_Z_mul(pol, lc): pol;
     193             : 
     194             :   {
     195       49623 :     GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;
     196             : 
     197       49623 :     pa_b = powiu(p, a-b); /* < 2^31 */
     198       49623 :     pa_bs2 = shifti(pa_b,-1);
     199       49623 :     pb= powiu(p, b);
     200      175700 :     for (i=1; i <= lfamod; i++)
     201             :     {
     202      126077 :       GEN T1,T2, P = gel(famod,i);
     203      126077 :       long d = degpol(P);
     204             : 
     205      126077 :       deg[i] = d; P += 2;
     206      126077 :       T1 = gel(P,d-1);/* = - S_1 */
     207      126077 :       T2 = sqri(T1);
     208      126077 :       if (d > 1) T2 = subii(T2, shifti(gel(P,d-2),1));
     209      126077 :       T2 = modii(T2, pa); /* = S_2 Newton sum */
     210      126077 :       if (lc)
     211             :       {
     212        4193 :         T1 = Fp_mul(lc, T1, pa);
     213        4193 :         T2 = Fp_mul(lc2,T2, pa);
     214             :       }
     215      126077 :       uel(trace1,i) = itou(diviiround(T1, pb));
     216      126077 :       uel(trace2,i) = itou(diviiround(T2, pb));
     217             :     }
     218       49623 :     spa_b   = uel(pa_b,2); /* < 2^31 */
     219       49623 :     spa_bs2 = uel(pa_bs2,2); /* < 2^31 */
     220             :   }
     221       49623 :   degsofar[0] = 0; /* sentinel */
     222             : 
     223             :   /* ind runs through strictly increasing sequences of length K,
     224             :    * 1 <= ind[i] <= lfamod */
     225       90540 : nextK:
     226       90540 :   if (K > *pmaxK || 2*K > lfamod) goto END;
     227       55200 :   if (DEBUGLEVEL > 3)
     228           0 :     err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
     229       55200 :   setlg(ind, K+1); ind[1] = 1;
     230       55200 :   Sbound = (ulong) ((K+1)>>1);
     231       55200 :   i = 1; curdeg = deg[ind[1]];
     232             :   for(;;)
     233             :   { /* try all combinations of K factors */
     234      622734 :     for (j = i; j < K; j++)
     235             :     {
     236       85886 :       degsofar[j] = curdeg;
     237       85886 :       ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
     238             :     }
     239      536848 :     if (curdeg <= klim) /* trial divide */
     240       10504 :     {
     241             :       GEN y, q, list;
     242             :       pari_sp av;
     243             :       ulong t;
     244             : 
     245             :       /* d - 1 test */
     246     1353090 :       for (t=uel(trace1,ind[1]),i=2; i<=K; i++)
     247      816242 :         t = Fl_add(t, uel(trace1,ind[i]), spa_b);
     248      536848 :       if (t > spa_bs2) t = spa_b - t;
     249      536848 :       if (t > Sbound)
     250             :       {
     251      428220 :         if (DEBUGLEVEL>6) err_printf(".");
     252      428220 :         goto NEXT;
     253             :       }
     254             :       /* d - 2 test */
     255      233889 :       for (t=uel(trace2,ind[1]),i=2; i<=K; i++)
     256      125261 :         t = Fl_add(t, uel(trace2,ind[i]), spa_b);
     257      108628 :       if (t > spa_bs2) t = spa_b - t;
     258      108628 :       if (t > Sbound)
     259             :       {
     260       57582 :         if (DEBUGLEVEL>6) err_printf("|");
     261       57582 :         goto NEXT;
     262             :       }
     263             : 
     264       51046 :       av = avma;
     265             :       /* check trailing coeff */
     266       51046 :       y = lc;
     267      150217 :       for (i=1; i<=K; i++)
     268             :       {
     269       99171 :         GEN q = constant_coeff(gel(famod,ind[i]));
     270       99171 :         if (y) q = mulii(y, q);
     271       99171 :         y = centermodii(q, pa, pas2);
     272             :       }
     273       51046 :       if (!signe(y) || !dvdii(constant_coeff(lcpol), y))
     274             :       {
     275       22556 :         if (DEBUGLEVEL>3) err_printf("T");
     276       22556 :         set_avma(av); goto NEXT;
     277             :       }
     278       28490 :       y = lc; /* full computation */
     279       63602 :       for (i=1; i<=K; i++)
     280             :       {
     281       35112 :         GEN q = gel(famod,ind[i]);
     282       35112 :         if (y) q = gmul(y, q);
     283       35112 :         y = centermod_i(q, pa, pas2);
     284             :       }
     285             : 
     286             :       /* y is the candidate factor */
     287       28490 :       if (! (q = ZX_divides_i(lcpol,y,bound)) )
     288             :       {
     289        3703 :         if (DEBUGLEVEL>3) err_printf("*");
     290        3703 :         set_avma(av); goto NEXT;
     291             :       }
     292             :       /* found a factor */
     293       24787 :       list = cgetg(K+1, t_VEC);
     294       24787 :       gel(listmod,cnt) = list;
     295       50288 :       for (i=1; i<=K; i++) list[i] = famod[ind[i]];
     296             : 
     297       24787 :       y = Q_primpart(y);
     298       24787 :       gel(fa,cnt++) = y;
     299             :       /* fix up pol */
     300       24787 :       pol = q;
     301       24787 :       if (lc) pol = Q_div_to_int(pol, leading_coeff(y));
     302       91942 :       for (i=j=k=1; i <= lfamod; i++)
     303             :       { /* remove used factors */
     304       67155 :         if (j <= K && i == ind[j]) j++;
     305             :         else
     306             :         {
     307       41654 :           gel(famod,k) = gel(famod,i);
     308       41654 :           uel(trace1,k) = uel(trace1,i);
     309       41654 :           uel(trace2,k) = uel(trace2,i);
     310       41654 :           deg[k] = deg[i]; k++;
     311             :         }
     312             :       }
     313       24787 :       lfamod -= K;
     314       24787 :       *pmaxK = cmbf_maxK(lfamod);
     315       24787 :       if (lfamod < 2*K) goto END;
     316       10504 :       i = 1; curdeg = deg[ind[1]];
     317       10504 :       bound = factor_bound(pol);
     318       10504 :       if (lc) lc = absi_shallow(leading_coeff(pol));
     319       10504 :       lcpol = lc? ZX_Z_mul(pol, lc): pol;
     320       10504 :       if (DEBUGLEVEL>3)
     321           0 :         err_printf("\nfound factor %Ps\nremaining modular factor(s): %ld\n",
     322             :                    y, lfamod);
     323       10504 :       continue;
     324             :     }
     325             : 
     326           0 : NEXT:
     327      512061 :     for (i = K+1;;)
     328             :     {
     329      638150 :       if (--i == 0) { K++; goto nextK; }
     330      597233 :       if (++ind[i] <= lfamod - K + i)
     331             :       {
     332      471144 :         curdeg = degsofar[i-1] + deg[ind[i]];
     333      471144 :         if (curdeg <= klim) break;
     334             :       }
     335             :     }
     336             :   }
     337       49623 : END:
     338       49623 :   *done = 1;
     339       49623 :   if (degpol(pol) > 0)
     340             :   { /* leftover factor */
     341       49623 :     if (signe(leading_coeff(pol)) < 0) pol = ZX_neg(pol);
     342       49623 :     if (lfamod >= 2*K) *done = 0;
     343             : 
     344       49623 :     setlg(famod, lfamod+1);
     345       49623 :     gel(listmod,cnt) = leafcopy(famod);
     346       49623 :     gel(fa,cnt++) = pol;
     347             :   }
     348       49623 :   if (DEBUGLEVEL>6) err_printf("\n");
     349       49623 :   setlg(listmod, cnt);
     350       49623 :   setlg(fa, cnt); return mkvec2(fa, listmod);
     351             : }
     352             : 
     353             : /* recombination of modular factors: van Hoeij's algorithm */
     354             : 
     355             : /* Q in Z[X], return Q(2^n) */
     356             : static GEN
     357      183288 : shifteval(GEN Q, long n)
     358             : {
     359      183288 :   pari_sp av = avma;
     360      183288 :   long i, l = lg(Q);
     361             :   GEN s;
     362             : 
     363      183288 :   if (!signe(Q)) return gen_0;
     364      183288 :   s = gel(Q,l-1);
     365      970775 :   for (i = l-2; i > 1; i--)
     366             :   {
     367      787488 :     s = addii(gel(Q,i), shifti(s, n));
     368      787487 :     if (gc_needed(av,1)) s = gerepileuptoint(av, s);
     369             :   }
     370      183287 :   return s;
     371             : }
     372             : 
     373             : /* return integer y such that all |a| <= y if P(a) = 0 */
     374             : static GEN
     375      109860 : root_bound(GEN P0)
     376             : {
     377      109860 :   GEN Q = leafcopy(P0), lP = absi_shallow(leading_coeff(Q)), x,y,z;
     378      109860 :   long k, d = degpol(Q);
     379             : 
     380             :   /* P0 = lP x^d + Q, deg Q < d */
     381      109860 :   Q = normalizepol_lg(Q, d+2);
     382      691755 :   for (k=lg(Q)-1; k>1; k--) gel(Q,k) = absi_shallow(gel(Q,k));
     383      109860 :   k = (long)(fujiwara_bound(P0));
     384      184114 :   for (  ; k >= 0; k--)
     385             :   {
     386      183288 :     pari_sp av = avma;
     387             :     /* y = 2^k; Q(y) >= lP y^d ? */
     388      183288 :     if (cmpii(shifteval(Q,k), shifti(lP, d*k)) >= 0) break;
     389       74254 :     set_avma(av);
     390             :   }
     391      109859 :   if (k < 0) k = 0;
     392      109859 :   y = int2n(k+1);
     393      109859 :   if (d > 2000) return y; /* likely to be expensive, don't bother */
     394      109859 :   x = int2n(k);
     395      109859 :   for(k=0; ; k++)
     396             :   {
     397      598704 :     z = shifti(addii(x,y), -1);
     398      598702 :     if (equalii(x,z) || k > 5) break;
     399      488844 :     if (cmpii(ZX_Z_eval(Q,z), mulii(lP, powiu(z, d))) < 0)
     400      261595 :       y = z;
     401             :     else
     402      227250 :       x = z;
     403             :   }
     404      109860 :   return y;
     405             : }
     406             : 
     407             : GEN
     408         350 : chk_factors_get(GEN lt, GEN famod, GEN c, GEN T, GEN N)
     409             : {
     410         350 :   long i = 1, j, l = lg(famod);
     411         350 :   GEN V = cgetg(l, t_VEC);
     412        8414 :   for (j = 1; j < l; j++)
     413        8064 :     if (signe(gel(c,j))) gel(V,i++) = gel(famod,j);
     414         350 :   if (lt && i > 1) gel(V,1) = RgX_Rg_mul(gel(V,1), lt);
     415         350 :   setlg(V, i);
     416         350 :   return T? FpXQXV_prod(V, T, N): FpXV_prod(V,N);
     417             : }
     418             : 
     419             : static GEN
     420         140 : chk_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)
     421             : {
     422             :   long i, r;
     423         140 :   GEN pol = P, list, piv, y, ltpol, lt, paov2;
     424             : 
     425         140 :   piv = ZM_hnf_knapsack(M_L);
     426         140 :   if (!piv) return NULL;
     427          70 :   if (DEBUGLEVEL>7) err_printf("ZM_hnf_knapsack output:\n%Ps\n",piv);
     428             : 
     429          70 :   r  = lg(piv)-1;
     430          70 :   list = cgetg(r+1, t_VEC);
     431          70 :   lt = absi_shallow(leading_coeff(pol));
     432          70 :   if (equali1(lt)) lt = NULL;
     433          70 :   ltpol = lt? ZX_Z_mul(pol, lt): pol;
     434          70 :   paov2 = shifti(pa,-1);
     435          70 :   for (i = 1;;)
     436             :   {
     437         161 :     if (DEBUGLEVEL) err_printf("LLL_cmbf: checking factor %ld\n",i);
     438         161 :     y = chk_factors_get(lt, famod, gel(piv,i), NULL, pa);
     439         161 :     y = FpX_center_i(y, pa, paov2);
     440         161 :     if (! (pol = ZX_divides_i(ltpol,y,bound)) ) return NULL;
     441         133 :     if (lt) y = Q_primpart(y);
     442         133 :     gel(list,i) = y;
     443         133 :     if (++i >= r) break;
     444             : 
     445          91 :     if (lt)
     446             :     {
     447          35 :       pol = ZX_Z_divexact(pol, leading_coeff(y));
     448          35 :       lt = absi_shallow(leading_coeff(pol));
     449          35 :       ltpol = ZX_Z_mul(pol, lt);
     450             :     }
     451             :     else
     452          56 :       ltpol = pol;
     453             :   }
     454          42 :   y = Q_primpart(pol);
     455          42 :   gel(list,i) = y; return list;
     456             : }
     457             : 
     458             : GEN
     459        1547 : LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, long *ti_LLL)
     460             : {
     461             :   GEN norm, u;
     462             :   long i, R;
     463             :   pari_timer T;
     464             : 
     465        1547 :   if (DEBUGLEVEL>2) timer_start(&T);
     466        1547 :   u = ZM_lll_norms(m, final? 0.999: 0.75, LLL_INPLACE | LLL_NOFLATTER, &norm);
     467        1547 :   if (DEBUGLEVEL>2) *ti_LLL += timer_delay(&T);
     468       10542 :   for (R=lg(m)-1; R > 0; R--)
     469       10542 :     if (cmprr(gel(norm,R), Bnorm) < 0) break;
     470       16156 :   for (i=1; i<=R; i++) setlg(u[i], n0+1);
     471        1547 :   if (R <= 1)
     472             :   {
     473         105 :     if (!R) pari_err_BUG("LLL_cmbf [no factor]");
     474         105 :     return NULL; /* irreducible */
     475             :   }
     476        1442 :   setlg(u, R+1); return u;
     477             : }
     478             : 
     479             : static ulong
     480          14 : next2pow(ulong a)
     481             : {
     482          14 :   ulong b = 1;
     483         112 :   while (b < a) b <<= 1;
     484          14 :   return b;
     485             : }
     486             : 
     487             : /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
     488             :  * van Hoeij's knapsack
     489             :  *
     490             :  * P = squarefree in Z[X].
     491             :  * famod = array of (lifted) modular factors mod p^a
     492             :  * bound = Mignotte bound for the size of divisors of P (for the sup norm)
     493             :  * previously recombined all set of factors with less than rec elts */
     494             : static GEN
     495         126 : LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
     496             : {
     497         126 :   const long N0 = 1; /* # of traces added at each step */
     498         126 :   double BitPerFactor = 0.4; /* nb bits in p^(a-b) / modular factor */
     499         126 :   long i,j,tmax,n0,C, dP = degpol(P);
     500         126 :   double logp = log((double)itos(p)), LOGp2 = M_LN2/logp;
     501         126 :   double b0 = log((double)dP*2) / logp, logBr;
     502             :   GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
     503             :   pari_sp av, av2;
     504         126 :   long ti_LLL = 0, ti_CF  = 0;
     505             : 
     506         126 :   lP = absi_shallow(leading_coeff(P));
     507         126 :   if (equali1(lP)) lP = NULL;
     508         126 :   Br = root_bound(P);
     509         126 :   if (lP) Br = mulii(lP, Br);
     510         126 :   logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;
     511             : 
     512         126 :   n0 = lg(famod) - 1;
     513         126 :   C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */
     514         126 :   Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);
     515         126 :   ZERO = zeromat(n0, N0);
     516             : 
     517         126 :   av = avma;
     518         126 :   TT = cgetg(n0+1, t_VEC);
     519         126 :   Tra  = cgetg(n0+1, t_MAT);
     520        2401 :   for (i=1; i<=n0; i++)
     521             :   {
     522        2275 :     TT[i]  = 0;
     523        2275 :     gel(Tra,i) = cgetg(N0+1, t_COL);
     524             :   }
     525         126 :   CM_L = scalarmat_s(C, n0);
     526             :   /* tmax = current number of traces used (and computed so far) */
     527         126 :   for (tmax = 0;; tmax += N0)
     528         469 :   {
     529         595 :     long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;
     530             :     GEN M_L, q, CM_Lp, oldCM_L;
     531         595 :     int first = 1;
     532             :     pari_timer ti2, TI;
     533             : 
     534         595 :     bmin = (long)ceil(b0 + tnew*logBr);
     535         595 :     if (DEBUGLEVEL>2)
     536           0 :       err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
     537             :                  r, tmax, bmin);
     538             : 
     539             :     /* compute Newton sums (possibly relifting first) */
     540         595 :     if (a <= bmin)
     541             :     {
     542          14 :       a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */
     543          14 :       a = (long)next2pow((ulong)a);
     544             : 
     545          14 :       pa = powiu(p,a);
     546          14 :       famod = ZpX_liftfact(P, famod, pa, p, a);
     547         266 :       for (i=1; i<=n0; i++) TT[i] = 0;
     548             :     }
     549       11557 :     for (i=1; i<=n0; i++)
     550             :     {
     551       10962 :       GEN p1 = gel(Tra,i);
     552       10962 :       GEN p2 = polsym_gen(gel(famod,i), gel(TT,i), tnew, NULL, pa);
     553       10962 :       gel(TT,i) = p2;
     554       10962 :       p2 += 1+tmax; /* ignore traces number 0...tmax */
     555       21924 :       for (j=1; j<=N0; j++) gel(p1,j) = gel(p2,j);
     556       10962 :       if (lP)
     557             :       { /* make Newton sums integral */
     558        1848 :         GEN lPpow = powiu(lP, tmax);
     559        3696 :         for (j=1; j<=N0; j++)
     560             :         {
     561        1848 :           lPpow = mulii(lPpow,lP);
     562        1848 :           gel(p1,j) = mulii(gel(p1,j), lPpow);
     563             :         }
     564             :       }
     565             :     }
     566             : 
     567             :     /* compute truncation parameter */
     568         595 :     if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
     569         595 :     oldCM_L = CM_L;
     570         595 :     av2 = avma;
     571         595 :     delta = b = 0; /* -Wall */
     572        1071 : AGAIN:
     573        1071 :     M_L = Q_div_to_int(CM_L, utoipos(C));
     574        1071 :     T2 = centermod( ZM_mul(Tra, M_L), pa );
     575        1071 :     if (first)
     576             :     { /* initialize lattice, using few p-adic digits for traces */
     577         595 :       double t = gexpo(T2) - maxdd(32.0, BitPerFactor*r);
     578         595 :       bgood = (long) (t * LOGp2);
     579         595 :       b = maxss(bmin, bgood);
     580         595 :       delta = a - b;
     581             :     }
     582             :     else
     583             :     { /* add more p-adic digits and continue reduction */
     584         476 :       long b0 = (long)(gexpo(T2) * LOGp2);
     585         476 :       if (b0 < b) b = b0;
     586         476 :       b = maxss(b-delta, bmin);
     587         476 :       if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
     588             :     }
     589             : 
     590        1071 :     q = powiu(p, b);
     591        1071 :     m = vconcat( CM_L, gdivround(T2, q) );
     592        1071 :     if (first)
     593             :     {
     594         595 :       GEN P1 = scalarmat(powiu(p, a-b), N0);
     595         595 :       first = 0;
     596         595 :       m = shallowconcat( m, vconcat(ZERO, P1) );
     597             :       /*     [ C M_L        0     ]
     598             :        * m = [                    ]   square matrix
     599             :        *     [  T2'  p^(a-b) I_N0 ]   T2' = Tra * M_L  truncated
     600             :        */
     601             :     }
     602             : 
     603        1071 :     CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
     604        1071 :     if (DEBUGLEVEL>2)
     605           0 :       err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
     606           0 :                  a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
     607        1071 :     if (!CM_L) { list = mkvec(P); break; }
     608         987 :     if (b > bmin)
     609             :     {
     610         476 :       CM_L = gerepilecopy(av2, CM_L);
     611         476 :       goto AGAIN;
     612             :     }
     613         511 :     if (DEBUGLEVEL>2) timer_printf(&ti2, "for this block of traces");
     614             : 
     615         511 :     i = lg(CM_L) - 1;
     616         511 :     if (i == r && ZM_equal(CM_L, oldCM_L))
     617             :     {
     618          63 :       CM_L = oldCM_L;
     619          63 :       set_avma(av2); continue;
     620             :     }
     621             : 
     622         448 :     CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
     623         448 :     if (lg(CM_Lp) != lg(CM_L))
     624             :     {
     625           7 :       if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
     626           7 :       CM_L = ZM_hnf(CM_L);
     627             :     }
     628             : 
     629         448 :     if (i <= r && i*rec < n0)
     630             :     {
     631             :       pari_timer ti;
     632         140 :       if (DEBUGLEVEL>2) timer_start(&ti);
     633         140 :       list = chk_factors(P, Q_div_to_int(CM_L,utoipos(C)), bound, famod, pa);
     634         140 :       if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
     635         140 :       if (list) break;
     636          98 :       if (DEBUGLEVEL>2) err_printf("LLL_cmbf: chk_factors failed");
     637             :     }
     638         406 :     CM_L = gerepilecopy(av2, CM_L);
     639         406 :     if (gc_needed(av,1))
     640             :     {
     641           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"LLL_cmbf");
     642           0 :       gerepileall(av, 5, &CM_L, &TT, &Tra, &famod, &pa);
     643             :     }
     644             :   }
     645         126 :   if (DEBUGLEVEL>2)
     646           0 :     err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
     647         126 :   return list;
     648             : }
     649             : 
     650             : /* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */
     651             : static int
     652       49623 : cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)
     653             : {
     654       49623 :   long a,b,amin,d = (long)(31 * M_LN2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);
     655       49623 :   int fl = 0;
     656             : 
     657       49623 :   b = logintall(B, q, qb) + 1;
     658       49623 :   *qb = mulii(*qb, q);
     659       49623 :   amin = b + d;
     660       49623 :   if (gcmp(powiu(q, amin), A) <= 0)
     661             :   {
     662       14216 :     a = logintall(A, q, qa) + 1;
     663       14216 :     *qa = mulii(*qa, q);
     664       14216 :     b = a - d; *qb = powiu(q, b);
     665             :   }
     666             :   else
     667             :   { /* not enough room */
     668       35407 :     a = amin;  *qa = powiu(q, a);
     669       35407 :     fl = 1;
     670             :   }
     671       49623 :   if (DEBUGLEVEL > 3) {
     672           0 :     err_printf("S_2   bound: %Ps^%ld\n", q,b);
     673           0 :     err_printf("coeff bound: %Ps^%ld\n", q,a);
     674             :   }
     675       49623 :   *pta = a;
     676       49623 :   *ptb = b; return fl;
     677             : }
     678             : 
     679             : /* use van Hoeij's knapsack algorithm */
     680             : static GEN
     681       49623 : combine_factors(GEN target, GEN famod, GEN p, long klim)
     682             : {
     683             :   GEN la, B, A, res, L, pa, pb, listmod;
     684       49623 :   long a,b, l, maxK, n = degpol(target);
     685             :   int done;
     686             :   pari_timer T;
     687             : 
     688       49623 :   A = factor_bound(target);
     689             : 
     690       49623 :   la = absi_shallow(leading_coeff(target));
     691       49623 :   B = mului(n, sqri(mulii(la, root_bound(target)))); /* = bound for S_2 */
     692             : 
     693       49623 :   (void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);
     694             : 
     695       49623 :   if (DEBUGLEVEL>2) timer_start(&T);
     696       49623 :   famod = ZpX_liftfact(target, famod, pa, p, a);
     697       49623 :   if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %Ps^%ld)", p,a);
     698       49623 :   L = cmbf(target, famod, A, p, a, b, klim, &maxK, &done);
     699       49623 :   if (DEBUGLEVEL>2) timer_printf(&T, "Naive recombination");
     700             : 
     701       49623 :   res     = gel(L,1);
     702       49623 :   listmod = gel(L,2); l = lg(listmod)-1;
     703       49623 :   famod = gel(listmod,l);
     704       49623 :   if (maxK > 0 && lg(famod)-1 > 2*maxK)
     705             :   {
     706         126 :     if (l!=1) A = factor_bound(gel(res,l));
     707         126 :     if (DEBUGLEVEL > 4) err_printf("last factor still to be checked\n");
     708         126 :     L = LLL_cmbf(gel(res,l), famod, p, pa, A, a, maxK);
     709         126 :     if (DEBUGLEVEL>2) timer_printf(&T,"Knapsack");
     710             :     /* remove last elt, possibly unfactored. Add all new ones. */
     711         126 :     setlg(res, l); res = shallowconcat(res, L);
     712             :   }
     713       49623 :   return res;
     714             : }
     715             : 
     716             : /* Assume 'a' a squarefree ZX; return 0 if no root (fl=1) / irreducible (fl=0).
     717             :  * Otherwise return prime p such that a mod p has fewest roots / factors */
     718             : static ulong
     719     1824042 : pick_prime(GEN a, long fl, pari_timer *T)
     720             : {
     721     1824042 :   pari_sp av = avma, av1;
     722     1824042 :   const long MAXNP = 7, da = degpol(a);
     723     1824042 :   long nmax = da+1, np;
     724     1824042 :   ulong chosenp = 0;
     725     1824042 :   GEN lead = gel(a,da+2);
     726             :   forprime_t S;
     727     1824042 :   if (equali1(lead)) lead = NULL;
     728     1824041 :   u_forprime_init(&S, 2, ULONG_MAX);
     729     1824046 :   av1 = avma;
     730     7402312 :   for (np = 0; np < MAXNP; set_avma(av1))
     731             :   {
     732     7292574 :     ulong p = u_forprime_next(&S);
     733             :     long nfacp;
     734             :     GEN z;
     735             : 
     736     7292570 :     if (!p) pari_err_OVERFLOW("DDF [out of small primes]");
     737     7292570 :     if (lead && !umodiu(lead,p)) continue;
     738     7226136 :     z = ZX_to_Flx(a, p);
     739     7226124 :     if (!Flx_is_squarefree(z, p)) continue;
     740             : 
     741     4494249 :     if (fl==1)
     742             :     {
     743     3909263 :       nfacp = Flx_nbroots(z, p);
     744     3909263 :       if (!nfacp) { chosenp = 0; break; } /* no root */
     745             :     }
     746      584986 :     else if(fl==0)
     747             :     {
     748      584237 :       nfacp = Flx_nbfact(z, p);
     749      584245 :       if (nfacp == 1) { chosenp = 0; break; } /* irreducible */
     750             :     } else
     751             :     {
     752         749 :       GEN f = gel(Flx_degfact(z, p),1);
     753         749 :       nfacp = lg(f)-1;
     754         749 :       if (f[1] > fl) { chosenp = 0; break; } /* no small factors */
     755             :     }
     756     2779946 :     if (DEBUGLEVEL>4)
     757           0 :       err_printf("...tried prime %3lu (%-3ld %s). Time = %ld\n",
     758             :                   p, nfacp, fl==1? "roots": "factors", timer_delay(T));
     759     2779948 :     if (nfacp < nmax)
     760             :     {
     761      992458 :       nmax = nfacp; chosenp = p;
     762      992458 :       if (da > 100 && nmax < 5) break; /* large degree, few factors. Enough */
     763             :     }
     764     2779941 :     np++;
     765             :   }
     766     1824054 :   return gc_ulong(av, chosenp);
     767             : }
     768             : 
     769             : /* Assume A a squarefree ZX; return the vector of its rational roots */
     770             : static GEN
     771     1621056 : DDF_roots(GEN A)
     772             : {
     773             :   GEN p, lc, lcpol, z, pe, pes2, bound;
     774             :   long i, m, e, lz;
     775             :   ulong pp;
     776             :   pari_sp av;
     777             :   pari_timer T;
     778             : 
     779     1621056 :   if (DEBUGLEVEL>2) timer_start(&T);
     780     1621056 :   pp = pick_prime(A, 1, &T);
     781     1621058 :   if (!pp) return cgetg(1,t_COL); /* no root */
     782       60111 :   p = utoipos(pp);
     783       60111 :   lc = leading_coeff(A);
     784       60111 :   if (is_pm1(lc))
     785       52599 :   { lc = NULL; lcpol = A; }
     786             :   else
     787        7512 :   { lc = absi_shallow(lc); lcpol = ZX_Z_mul(A, lc); }
     788       60111 :   bound = root_bound(A); if (lc) bound = mulii(lc, bound);
     789       60111 :   e = logintall(addiu(shifti(bound, 1), 1), p, &pe) + 1;
     790       60111 :   pe = mulii(pe, p);
     791       60111 :   pes2 = shifti(pe, -1);
     792       60111 :   if (DEBUGLEVEL>2) timer_printf(&T, "Root bound");
     793       60111 :   av = avma;
     794       60111 :   z = ZpX_roots(A, p, e); lz = lg(z);
     795       60111 :   z = deg1_from_roots(z, varn(A));
     796       60111 :   if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %lu^%ld)", pp,e);
     797      130547 :   for (m=1, i=1; i < lz; i++)
     798             :   {
     799       70436 :     GEN q, r, y = gel(z,i);
     800       70436 :     if (lc) y = ZX_Z_mul(y, lc);
     801       70436 :     y = centermod_i(y, pe, pes2);
     802       70436 :     if (! (q = ZX_divides(lcpol, y)) ) continue;
     803             : 
     804       25946 :     lcpol = q;
     805       25946 :     r = negi( constant_coeff(y) );
     806       25946 :     if (lc) {
     807        8277 :       r = gdiv(r,lc);
     808        8277 :       lcpol = Q_primpart(lcpol);
     809        8277 :       lc = absi_shallow( leading_coeff(lcpol) );
     810        8277 :       if (is_pm1(lc)) lc = NULL; else lcpol = ZX_Z_mul(lcpol, lc);
     811             :     }
     812       25946 :     gel(z,m++) = r;
     813       25946 :     if (gc_needed(av,2))
     814             :     {
     815           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"DDF_roots, m = %ld", m);
     816           0 :       gerepileall(av, lc? 3:2, &z, &lcpol, &lc);
     817             :     }
     818             :   }
     819       60111 :   if (DEBUGLEVEL>2) timer_printf(&T, "Recombination");
     820       60111 :   setlg(z, m); return z;
     821             : }
     822             : 
     823             : /* Assume a squarefree ZX, deg(a) > 0, return rational factors.
     824             :  * In fact, a(0) != 0 but we don't use this
     825             :  * if dmax>0, Only look for factor of degree at most dmax */
     826             : GEN
     827      202986 : ZX_DDF_max(GEN a, long dmax)
     828             : {
     829             :   GEN ap, prime, famod, z;
     830      202986 :   long ti = 0;
     831      202986 :   ulong p = 0;
     832      202986 :   pari_sp av = avma;
     833             :   pari_timer T, T2;
     834             : 
     835      202986 :   if (DEBUGLEVEL>2) { timer_start(&T); timer_start(&T2); }
     836      202986 :   p = pick_prime(a, dmax, &T2);
     837      202987 :   if (!p) return mkvec(a);
     838       49623 :   prime = utoipos(p);
     839       49623 :   ap = Flx_normalize(ZX_to_Flx(a, p), p);
     840       49623 :   famod = gel(Flx_factor(ap, p), 1);
     841       49623 :   if (DEBUGLEVEL>2)
     842             :   {
     843           0 :     if (DEBUGLEVEL>4) timer_printf(&T2, "splitting mod p = %lu", p);
     844           0 :     ti = timer_delay(&T);
     845           0 :     err_printf("Time setup: %ld\n", ti);
     846             :   }
     847       49623 :   z = combine_factors(a, FlxV_to_ZXV(famod), prime, degpol(a)-1);
     848       49623 :   if (DEBUGLEVEL>2)
     849           0 :     err_printf("Total Time: %ld\n===========\n", ti + timer_delay(&T));
     850       49623 :   return gerepilecopy(av, z);
     851             : }
     852             : 
     853             : /* Distinct Degree Factorization (deflating first)
     854             :  * Assume x squarefree, degree(x) > 0, x(0) != 0 */
     855             : GEN
     856      151491 : ZX_DDF(GEN x)
     857             : {
     858             :   GEN L;
     859             :   long m;
     860      151491 :   x = ZX_deflate_max(x, &m);
     861      151491 :   L = ZX_DDF_max(x,0);
     862      151491 :   if (m > 1)
     863             :   {
     864       49195 :     GEN e, v, fa = factoru(m);
     865             :     long i,j,k, l;
     866             : 
     867       49195 :     e = gel(fa,2); k = 0;
     868       49195 :     fa= gel(fa,1); l = lg(fa);
     869       98733 :     for (i=1; i<l; i++) k += e[i];
     870       49195 :     v = cgetg(k+1, t_VECSMALL); k = 1;
     871       98733 :     for (i=1; i<l; i++)
     872      100501 :       for (j=1; j<=e[i]; j++) v[k++] = fa[i];
     873      100159 :     for (k--; k; k--)
     874             :     {
     875       50963 :       GEN L2 = cgetg(1,t_VEC);
     876      102291 :       for (i=1; i < lg(L); i++)
     877       51327 :               L2 = shallowconcat(L2, ZX_DDF_max(RgX_inflate(gel(L,i), v[k]),0));
     878       50964 :       L = L2;
     879             :     }
     880             :   }
     881      151492 :   return L;
     882             : }
     883             : 
     884             : /* SquareFree Factorization in Z[X] (char 0 is enough, if ZX_gcd -> RgX_gcd)
     885             :  * f = prod Q[i]^E[i], E[1] < E[2] < ..., and Q[i] squarefree and coprime.
     886             :  * Return Q, set *pE = E. For efficiency, caller should have used ZX_valrem
     887             :  * so that f(0) != 0 */
     888             : GEN
     889      332424 : ZX_squff(GEN f, GEN *pE)
     890             : {
     891             :   GEN T, V, P, E;
     892      332424 :   long i, k, n = 1 + degpol(f);
     893             : 
     894      332424 :   if (signe(leading_coeff(f)) < 0) f = ZX_neg(f);
     895      332424 :   E = cgetg(n, t_VECSMALL);
     896      332424 :   P = cgetg(n, t_COL);
     897      332424 :   T = ZX_gcd_all(f, ZX_deriv(f), &V);
     898      332424 :   for (k = i = 1;; k++)
     899        3641 :   {
     900      336065 :     GEN W = ZX_gcd_all(T,V, &T); /* V and W are squarefree */
     901      336065 :     long dW = degpol(W), dV = degpol(V);
     902             :     /* f = prod_i T_i^{e_i}
     903             :      * W = prod_{i: e_i > k} T_i,
     904             :      * V = prod_{i: e_i >= k} T_i,
     905             :      * T = prod_{i: e_i > k} T_i^{e_i - k} */
     906      336065 :     if (!dW)
     907             :     {
     908      332424 :       if (dV) { gel(P,i) = Q_primpart(V); E[i] = k; i++; }
     909      332424 :       break;
     910             :     }
     911        3641 :     if (dW == dV)
     912             :     {
     913             :       GEN U;
     914        1638 :       while ( (U = ZX_divides(T, V)) ) { k++; T = U; }
     915             :     }
     916             :     else
     917             :     {
     918        2388 :       gel(P,i) = Q_primpart(RgX_div(V,W));
     919        2388 :       E[i] = k; i++; V = W;
     920             :     }
     921             :   }
     922      332424 :   setlg(P,i);
     923      332424 :   setlg(E,i); *pE = E; return P;
     924             : }
     925             : 
     926             : static GEN
     927       50323 : fact_from_DDF(GEN Q, GEN E, long n)
     928             : {
     929       50323 :   GEN v,w, y = cgetg(3, t_MAT);
     930       50323 :   long i,j,k, l = lg(Q);
     931             : 
     932       50323 :   v = cgetg(n+1, t_COL); gel(y,1) = v;
     933       50323 :   w = cgetg(n+1, t_COL); gel(y,2) = w;
     934      102445 :   for (k = i = 1; i < l; i++)
     935             :   {
     936       52122 :     GEN L = gel(Q,i), e = utoipos(E[i]);
     937       52122 :     long J = lg(L);
     938      128716 :     for (j = 1; j < J; j++,k++)
     939             :     {
     940       76594 :       gel(v,k) = ZX_copy(gel(L,j));
     941       76594 :       gel(w,k) = e;
     942             :     }
     943             :   }
     944       50323 :   return y;
     945             : }
     946             : 
     947             : /* Factor T in Z[x] */
     948             : static GEN
     949       50330 : ZX_factor_i(GEN T)
     950             : {
     951             :   GEN Q, E, y;
     952             :   long n, i, l, v;
     953             : 
     954       50330 :   if (!signe(T)) return prime_fact(T);
     955       50323 :   v = ZX_valrem(T, &T);
     956       50323 :   Q = ZX_squff(T, &E); l = lg(Q);
     957      100884 :   for (i = 1, n = 0; i < l; i++)
     958             :   {
     959       50561 :     gel(Q,i) = ZX_DDF(gel(Q,i));
     960       50561 :     n += lg(gel(Q,i)) - 1;
     961             :   }
     962       50323 :   if (v)
     963             :   {
     964        1561 :     Q = vec_append(Q, mkvec(pol_x(varn(T))));
     965        1561 :     E = vecsmall_append(E, v); n++;
     966             :   }
     967       50323 :   y = fact_from_DDF(Q, E, n);
     968       50323 :   return sort_factor_pol(y, cmpii);
     969             : }
     970             : GEN
     971       49553 : ZX_factor(GEN x)
     972             : {
     973       49553 :   pari_sp av = avma;
     974       49553 :   return gerepileupto(av, ZX_factor_i(x));
     975             : }
     976             : GEN
     977         777 : QX_factor(GEN x)
     978             : {
     979         777 :   pari_sp av = avma;
     980         777 :   return gerepileupto(av, ZX_factor_i(Q_primpart(x)));
     981             : }
     982             : 
     983             : long
     984      101762 : ZX_is_irred(GEN x)
     985             : {
     986      101762 :   pari_sp av = avma;
     987      101762 :   long l = lg(x);
     988             :   GEN y;
     989      101762 :   if (l <= 3) return 0; /* degree < 1 */
     990      101762 :   if (l == 4) return 1; /* degree 1 */
     991       98252 :   if (ZX_val(x)) return 0;
     992       98042 :   if (!ZX_is_squarefree(x)) return 0;
     993       97899 :   y = ZX_DDF(x); set_avma(av);
     994       97900 :   return (lg(y) == 2);
     995             : }
     996             : 
     997             : GEN
     998     1621055 : nfrootsQ(GEN x)
     999             : {
    1000     1621055 :   pari_sp av = avma;
    1001             :   GEN z;
    1002             :   long val;
    1003             : 
    1004     1621055 :   if (typ(x)!=t_POL) pari_err_TYPE("nfrootsQ",x);
    1005     1621055 :   if (!signe(x)) pari_err_ROOTS0("nfrootsQ");
    1006     1621055 :   x = Q_primpart(x);
    1007     1621054 :   RgX_check_ZX(x,"nfrootsQ");
    1008     1621054 :   val = ZX_valrem(x, &x);
    1009     1621054 :   z = DDF_roots( ZX_radical(x) );
    1010     1621059 :   if (val) z = vec_append(z, gen_0);
    1011     1621059 :   return gerepileupto(av, sort(z));
    1012             : }
    1013             : 
    1014             : /************************************************************************
    1015             :  *                   GCD OVER Z[X] / Q[X]                               *
    1016             :  ************************************************************************/
    1017             : int
    1018      197094 : ZX_is_squarefree(GEN x)
    1019             : {
    1020      197094 :   pari_sp av = avma;
    1021             :   GEN d;
    1022             :   long m;
    1023      197094 :   if (lg(x) == 2) return 0;
    1024      197094 :   m = ZX_deflate_order(x);
    1025      197094 :   if (m > 1)
    1026             :   {
    1027       85202 :     if (!signe(gel(x,2))) return 0;
    1028       84964 :     x = RgX_deflate(x, m);
    1029             :   }
    1030      196856 :   d = ZX_gcd(x,ZX_deriv(x));
    1031      196870 :   return gc_bool(av, lg(d) == 3);
    1032             : }
    1033             : 
    1034             : static int
    1035      114559 : ZX_gcd_filter(GEN *pt_A, GEN *pt_P)
    1036             : {
    1037      114559 :   GEN A = *pt_A, P = *pt_P;
    1038      114559 :   long i, j, l = lg(A), n = 1, d = degpol(gel(A,1));
    1039             :   GEN B, Q;
    1040      230124 :   for (i=2; i<l; i++)
    1041             :   {
    1042      115565 :     long di = degpol(gel(A,i));
    1043      115565 :     if (di==d) n++;
    1044          36 :     else if (d > di)
    1045          36 :     { n=1; d = di; }
    1046             :   }
    1047      114559 :   if (n == l-1)
    1048      114523 :     return 0;
    1049          36 :   B = cgetg(n+1, t_VEC);
    1050          36 :   Q = cgetg(n+1, typ(P));
    1051         156 :   for (i=1, j=1; i<l; i++)
    1052             :   {
    1053         120 :     if (degpol(gel(A,i))==d)
    1054             :     {
    1055          84 :       gel(B,j) = gel(A,i);
    1056          84 :       Q[j] = P[i];
    1057          84 :       j++;
    1058             :     }
    1059             :   }
    1060          36 :   *pt_A = B; *pt_P = Q; return 1;
    1061             : }
    1062             : 
    1063             : static GEN
    1064     3348155 : ZX_gcd_Flx(GEN a, GEN b, ulong g, ulong p)
    1065             : {
    1066     3348155 :   GEN H = Flx_gcd(a, b, p);
    1067     3348154 :   if (!g)
    1068     3326282 :     return Flx_normalize(H, p);
    1069             :   else
    1070             :   {
    1071       21872 :     ulong t = Fl_mul(g, Fl_inv(Flx_lead(H), p), p);
    1072       21872 :     return Flx_Fl_mul(H, t, p);
    1073             :   }
    1074             : }
    1075             : 
    1076             : static GEN
    1077     3342241 : ZX_gcd_slice(GEN A, GEN B, GEN g, GEN P, GEN *mod)
    1078             : {
    1079     3342241 :   pari_sp av = avma;
    1080     3342241 :   long i, n = lg(P)-1;
    1081             :   GEN H, T;
    1082     3342241 :   if (n == 1)
    1083             :   {
    1084     3337307 :     ulong p = uel(P,1), gp = g ? umodiu(g, p): 0;
    1085     3337307 :     GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
    1086     3337307 :     GEN Hp = ZX_gcd_Flx(a, b, gp, p);
    1087     3337306 :     H = gerepileupto(av, Flx_to_ZX(Hp));
    1088     3337306 :     *mod = utoi(p);
    1089     3337306 :     return H;
    1090             :   }
    1091        4934 :   T = ZV_producttree(P);
    1092        4934 :   A = ZX_nv_mod_tree(A, P, T);
    1093        4934 :   B = ZX_nv_mod_tree(B, P, T);
    1094        4934 :   g = g ?  Z_ZV_mod_tree(g, P, T): NULL;
    1095        4934 :   H = cgetg(n+1, t_VEC);
    1096       15782 :   for(i=1; i <= n; i++)
    1097             :   {
    1098       10848 :     ulong p = P[i];
    1099       10848 :     GEN a = gel(A,i), b = gel(B,i);
    1100       10848 :     gel(H,i) = ZX_gcd_Flx(a, b, g? g[i]: 0, p);
    1101             :   }
    1102        4934 :   if (ZX_gcd_filter(&H, &P))
    1103          12 :     T = ZV_producttree(P);
    1104        4934 :   H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
    1105        4934 :   *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
    1106             : }
    1107             : 
    1108             : GEN
    1109     3342241 : ZX_gcd_worker(GEN P, GEN A, GEN B, GEN g)
    1110             : {
    1111     3342241 :   GEN V = cgetg(3, t_VEC);
    1112     3342241 :   gel(V,1) = ZX_gcd_slice(A, B, equali1(g)? NULL: g , P, &gel(V,2));
    1113     3342240 :   return V;
    1114             : }
    1115             : 
    1116             : static GEN
    1117      109625 : ZX_gcd_chinese(GEN A, GEN P, GEN *mod)
    1118             : {
    1119      109625 :   ZX_gcd_filter(&A, &P);
    1120      109625 :   return nxV_chinese_center(A, P, mod);
    1121             : }
    1122             : 
    1123             : GEN
    1124    13099008 : ZX_gcd_all(GEN A, GEN B, GEN *Anew)
    1125             : {
    1126    13099008 :   pari_sp av = avma;
    1127    13099008 :   long k, valH, valA, valB, vA = varn(A), dA = degpol(A), dB = degpol(B);
    1128    13098975 :   GEN worker, c, cA, cB, g, Ag, Bg, H = NULL, mod = gen_1, R;
    1129             :   GEN Ap, Bp, Hp;
    1130             :   forprime_t S;
    1131             :   ulong pp;
    1132    13098975 :   if (dA < 0) { if (Anew) *Anew = pol_0(vA); return ZX_copy(B); }
    1133    13098954 :   if (dB < 0) { if (Anew) *Anew = pol_1(vA); return ZX_copy(A); }
    1134    13097806 :   A = Q_primitive_part(A, &cA);
    1135    13097928 :   B = Q_primitive_part(B, &cB);
    1136    13097892 :   valA = ZX_valrem(A, &A); dA -= valA;
    1137    13097888 :   valB = ZX_valrem(B, &B); dB -= valB;
    1138    13097932 :   valH = minss(valA, valB);
    1139    13097950 :   valA -= valH; /* valuation(Anew) */
    1140    13097950 :   c = (cA && cB)? gcdii(cA, cB): NULL; /* content(gcd) */
    1141    13097957 :   if (!dA || !dB)
    1142             :   {
    1143     6482112 :     if (Anew) *Anew = RgX_shift_shallow(A, valA);
    1144     6482112 :     return monomial(c? c: gen_1, valH, vA);
    1145             :   }
    1146     6615845 :   g = gcdii(leading_coeff(A), leading_coeff(B)); /* multiple of lead(gcd) */
    1147     6615765 :   if (is_pm1(g)) {
    1148     6409435 :     g = NULL;
    1149     6409435 :     Ag = A;
    1150     6409435 :     Bg = B;
    1151             :   } else {
    1152      206321 :     Ag = ZX_Z_mul(A,g);
    1153      206321 :     Bg = ZX_Z_mul(B,g);
    1154             :   }
    1155     6615756 :   init_modular_big(&S);
    1156             :   do {
    1157     6615848 :     pp = u_forprime_next(&S);
    1158     6615850 :     Ap = ZX_to_Flx(Ag, pp);
    1159     6615880 :     Bp = ZX_to_Flx(Bg, pp);
    1160     6615880 :   } while (degpol(Ap) != dA || degpol(Bp) != dB);
    1161     6615862 :   if (degpol(Flx_gcd(Ap, Bp, pp)) == 0)
    1162             :   {
    1163     3383244 :     if (Anew) *Anew = RgX_shift_shallow(A, valA);
    1164     3383245 :     return monomial(c? c: gen_1, valH, vA);
    1165             :   }
    1166     3232590 :   worker = snm_closure(is_entry("_ZX_gcd_worker"), mkvec3(A, B, g? g: gen_1));
    1167     3232590 :   av = avma;
    1168     3340561 :   for (k = 1; ;k *= 2)
    1169             :   {
    1170     3340561 :     gen_inccrt_i("ZX_gcd", worker, g, (k+1)>>1, 0, &S, &H, &mod, ZX_gcd_chinese, NULL);
    1171     3340562 :     gerepileall(av, 2, &H, &mod);
    1172     3340562 :     Hp = ZX_to_Flx(H, pp);
    1173     3340562 :     if (lgpol(Flx_rem(Ap, Hp, pp)) || lgpol(Flx_rem(Bp, Hp, pp))) continue;
    1174     3236656 :     if (!ZX_divides(Bg, H)) continue;
    1175     3232618 :     R = ZX_divides(Ag, H);
    1176     3232618 :     if (R) break;
    1177             :   }
    1178             :   /* lead(H) = g */
    1179     3232589 :   if (g) H = Q_primpart(H);
    1180     3232589 :   if (c) H = ZX_Z_mul(H,c);
    1181     3232589 :   if (DEBUGLEVEL>5) err_printf("done\n");
    1182     3232589 :   if (Anew) *Anew = RgX_shift_shallow(R, valA);
    1183     3232589 :   return valH? RgX_shift_shallow(H, valH): H;
    1184             : }
    1185             : 
    1186             : #if 0
    1187             : /* ceil( || p ||_oo / lc(p) ) */
    1188             : static GEN
    1189             : maxnorm(GEN p)
    1190             : {
    1191             :   long i, n = degpol(p), av = avma;
    1192             :   GEN x, m = gen_0;
    1193             : 
    1194             :   p += 2;
    1195             :   for (i=0; i<n; i++)
    1196             :   {
    1197             :     x = gel(p,i);
    1198             :     if (abscmpii(x,m) > 0) m = x;
    1199             :   }
    1200             :   m = divii(m, gel(p,n));
    1201             :   return gerepileuptoint(av, addiu(absi_shallow(m),1));
    1202             : }
    1203             : #endif
    1204             : 
    1205             : GEN
    1206     9767031 : ZX_gcd(GEN A, GEN B)
    1207             : {
    1208     9767031 :   pari_sp av = avma;
    1209     9767031 :   return gerepilecopy(av, ZX_gcd_all(A,B,NULL));
    1210             : }
    1211             : 
    1212             : GEN
    1213     2659225 : ZX_radical(GEN A) { GEN B; (void)ZX_gcd_all(A,ZX_deriv(A),&B); return B; }
    1214             : 
    1215             : static GEN
    1216       18634 : _gcd(GEN a, GEN b)
    1217             : {
    1218       18634 :   if (!a) a = gen_1;
    1219       18634 :   if (!b) b = gen_1;
    1220       18634 :   return Q_gcd(a,b);
    1221             : }
    1222             : /* A0 and B0 in Q[X] */
    1223             : GEN
    1224       18634 : QX_gcd(GEN A0, GEN B0)
    1225             : {
    1226             :   GEN a, b, D;
    1227       18634 :   pari_sp av = avma, av2;
    1228             : 
    1229       18634 :   D = ZX_gcd(Q_primitive_part(A0, &a), Q_primitive_part(B0, &b));
    1230       18634 :   av2 = avma; a = _gcd(a,b);
    1231       18634 :   if (isint1(a)) set_avma(av2); else D = ZX_Q_mul(D, a);
    1232       18634 :   return gerepileupto(av, D);
    1233             : }
    1234             : 
    1235             : /*****************************************************************************
    1236             :  * Variants of the Bradford-Davenport algorithm: look for cyclotomic         *
    1237             :  * factors, and decide whether a ZX is cyclotomic or a product of cyclotomic *
    1238             :  *****************************************************************************/
    1239             : /* f of degree 1, return a cyclotomic factor (Phi_1 or Phi_2) or NULL */
    1240             : static GEN
    1241           0 : BD_deg1(GEN f)
    1242             : {
    1243           0 :   GEN a = gel(f,3), b = gel(f,2); /* f = ax + b */
    1244           0 :   if (!absequalii(a,b)) return NULL;
    1245           0 :   return polcyclo((signe(a) == signe(b))? 2: 1, varn(f));
    1246             : }
    1247             : 
    1248             : /* f a squarefree ZX; not divisible by any Phi_n, n even */
    1249             : static GEN
    1250         406 : BD_odd(GEN f)
    1251             : {
    1252         406 :   while(degpol(f) > 1)
    1253             :   {
    1254         406 :     GEN f1 = ZX_graeffe(f); /* contain all cyclotomic divisors of f */
    1255         406 :     if (ZX_equal(f1, f)) return f; /* product of cyclotomics */
    1256           0 :     f = ZX_gcd(f, f1);
    1257             :   }
    1258           0 :   if (degpol(f) == 1) return BD_deg1(f);
    1259           0 :   return NULL; /* no cyclotomic divisor */
    1260             : }
    1261             : 
    1262             : static GEN
    1263        2310 : myconcat(GEN v, GEN x)
    1264             : {
    1265        2310 :   if (typ(x) != t_VEC) x = mkvec(x);
    1266        2310 :   if (!v) return x;
    1267        1470 :   return shallowconcat(v, x);
    1268             : }
    1269             : 
    1270             : /* Bradford-Davenport algorithm.
    1271             :  * f a squarefree ZX of degree > 0, return NULL or a vector of coprime
    1272             :  * cyclotomic factors of f [ possibly reducible ] */
    1273             : static GEN
    1274        2359 : BD(GEN f)
    1275             : {
    1276        2359 :   GEN G = NULL, Gs = NULL, Gp = NULL, Gi = NULL;
    1277             :   GEN fs2, fp, f2, f1, fe, fo, fe1, fo1;
    1278        2359 :   RgX_even_odd(f, &fe, &fo);
    1279        2359 :   fe1 = ZX_eval1(fe);
    1280        2359 :   fo1 = ZX_eval1(fo);
    1281        2359 :   if (absequalii(fe1, fo1)) /* f(1) = 0 or f(-1) = 0 */
    1282             :   {
    1283        1519 :     long i, v = varn(f);
    1284        1519 :     if (!signe(fe1))
    1285         371 :       G = mkvec2(polcyclo(1, v), polcyclo(2, v)); /* both 0 */
    1286        1148 :     else if (signe(fe1) == signe(fo1))
    1287         693 :       G = mkvec(polcyclo(2, v)); /*f(-1) = 0*/
    1288             :     else
    1289         455 :       G = mkvec(polcyclo(1, v)); /*f(1) = 0*/
    1290        3409 :     for (i = lg(G)-1; i; i--) f = RgX_div(f, gel(G,i));
    1291             :   }
    1292             :   /* f no longer divisible by Phi_1 or Phi_2 */
    1293        2359 :   if (degpol(f) <= 1) return G;
    1294        2058 :   f1 = ZX_graeffe(f); /* has at most square factors */
    1295        2058 :   if (ZX_equal(f1, f)) return myconcat(G,f); /* f = product of Phi_n, n odd */
    1296             : 
    1297        1183 :   fs2 = ZX_gcd_all(f1, ZX_deriv(f1), &f2); /* fs2 squarefree */
    1298        1183 :   if (degpol(fs2))
    1299             :   { /* fs contains all Phi_n | f, 4 | n; and only those */
    1300             :     /* In that case, Graeffe(Phi_n) = Phi_{n/2}^2, and Phi_n = Phi_{n/2}(x^2) */
    1301        1029 :     GEN fs = RgX_inflate(fs2, 2);
    1302        1029 :     (void)ZX_gcd_all(f, fs, &f); /* remove those Phi_n | f, 4 | n */
    1303        1029 :     Gs = BD(fs2);
    1304        1029 :     if (Gs)
    1305             :     {
    1306             :       long i;
    1307        2555 :       for (i = lg(Gs)-1; i; i--) gel(Gs,i) = RgX_inflate(gel(Gs,i), 2);
    1308             :       /* prod Gs[i] is the product of all Phi_n | f, 4 | n */
    1309        1029 :       G = myconcat(G, Gs);
    1310             :     }
    1311             :     /* f2 = f1 / fs2 */
    1312        1029 :     f1 = RgX_div(f2, fs2); /* f1 / fs2^2 */
    1313             :   }
    1314        1183 :   fp = ZX_gcd(f, f1); /* contains all Phi_n | f, n > 1 odd; and only those */
    1315        1183 :   if (degpol(fp))
    1316             :   {
    1317         196 :     Gp = BD_odd(fp);
    1318             :     /* Gp is the product of all Phi_n | f, n odd */
    1319         196 :     if (Gp) G = myconcat(G, Gp);
    1320         196 :     f = RgX_div(f, fp);
    1321             :   }
    1322        1183 :   if (degpol(f))
    1323             :   { /* contains all Phi_n originally dividing f, n = 2 mod 4, n > 2;
    1324             :      * and only those
    1325             :      * In that case, Graeffe(Phi_n) = Phi_{n/2}, and Phi_n = Phi_{n/2}(-x) */
    1326         210 :     Gi = BD_odd(ZX_z_unscale(f, -1));
    1327         210 :     if (Gi)
    1328             :     { /* N.B. Phi_2 does not divide f */
    1329         210 :       Gi = ZX_z_unscale(Gi, -1);
    1330             :       /* Gi is the product of all Phi_n | f, n = 2 mod 4 */
    1331         210 :       G = myconcat(G, Gi);
    1332             :     }
    1333             :   }
    1334        1183 :   return G;
    1335             : }
    1336             : 
    1337             : /* Let f be a nonzero QX, return the (squarefree) product of cyclotomic
    1338             :  * divisors of f */
    1339             : GEN
    1340         315 : polcyclofactors(GEN f)
    1341             : {
    1342         315 :   pari_sp av = avma;
    1343         315 :   if (typ(f) != t_POL || !signe(f)) pari_err_TYPE("polcyclofactors",f);
    1344         315 :   (void)RgX_valrem(f, &f);
    1345         315 :   f = Q_primpart(f);
    1346         315 :   RgX_check_ZX(f,"polcyclofactors");
    1347         315 :   if (degpol(f))
    1348             :   {
    1349         315 :     f = BD(ZX_radical(f));
    1350         315 :     if (f) return gerepilecopy(av, f);
    1351             :   }
    1352           0 :   set_avma(av); return cgetg(1,t_VEC);
    1353             : }
    1354             : 
    1355             : /* list of all squarefree odd x such that phi(x) = n, P^-(x) > m. Unsorted */
    1356             : static GEN
    1357       19547 : invphi(ulong n, ulong m)
    1358             : {
    1359             :   GEN C, D;
    1360             :   long l, i;
    1361       19547 :   if (n == 1) return mkvecsmall(1);
    1362       14280 :   D = divisorsu(n); l = lg(D);
    1363       14280 :   C = cgetg(1, t_VECSMALL);
    1364       39676 :   for (i = 2; i < l; i++) /* skip 1 */
    1365             :   {
    1366       25396 :     ulong d = D[i], p;
    1367       25396 :     if (d < m) continue;
    1368       20239 :     p = d + 1; if (!uisprime(p)) continue;
    1369       10507 :     C = vecsmall_concat(C, zv_z_mul(invphi(D[l-i], p), p));
    1370             :   }
    1371       14280 :   return C;
    1372             : }
    1373             : 
    1374             : long
    1375       98901 : poliscyclo(GEN f)
    1376             : {
    1377       98901 :   const ulong p = 2147483647; /* prime */
    1378             :   pari_sp av;
    1379             :   long i, n, e, l;
    1380             :   ulong f3, fm3;
    1381             :   GEN D, fp, _3;
    1382       98901 :   if (typ(f) != t_POL) pari_err_TYPE("poliscyclo", f);
    1383       98894 :   n = degpol(f);
    1384       98894 :   if (n <= 0 || !RgX_is_ZX(f)) return 0;
    1385       98887 :   if (!equali1(gel(f,n+2)) || !is_pm1(gel(f,2))) return 0;
    1386        9138 :   if (n == 1) return signe(gel(f,2)) > 0? 2: 1;
    1387        9040 :   av = avma;
    1388        9040 :   f = ZX_deflate_max(f, &e); if (e != 1) n = degpol(f);
    1389        9040 :   D = invphi(n, 1); /* squareefree odd d s.t. phi(d) = n */
    1390        9040 :   l = lg(D); _3 = gmodulss(3, p);
    1391        9040 :   fp = ZX_to_Flx(f, p);
    1392        9040 :   f3 = Flx_eval(fp, 3, p);
    1393        9040 :   fm3 = Flx_eval(fp, p-3, p);
    1394             :   /* f(x^e) is cyclotomic (= Phi_{de}) iff f = Phi_d, where all prime dividing
    1395             :    * e also divide d. */
    1396       11885 :   for (i = 1; i < l; i++)
    1397             :   {
    1398        5134 :     long d = D[i]; /* squarefree odd */
    1399        5134 :     if (odd(e))
    1400             :     {
    1401        4075 :       if (e == 1 || u_ppo(e, d) == 1)
    1402             :       { /* early abort: check whether f(3) = Phi_d(3) or Phi_2d(3) = Phi_d(-3)
    1403             :          * mod p before checking in Z. N.B. phi(d) and value at 3 mod p
    1404             :          * determine Phi_d for all d <= 10^7 */
    1405        3833 :         ulong F3 = Rg_to_Fl(polcyclo_eval(d, _3), p);
    1406        3833 :         if (F3 == f3 && ZX_equal(f, polcyclo(d, varn(f))))
    1407         973 :           return gc_long(av, d * e);
    1408        2860 :         if (F3 == fm3 && ZX_equal(f, polcyclo(2*d, varn(f))))
    1409         749 :           return gc_long(av, 2* d * e);
    1410             :       }
    1411             :     }
    1412             :     else
    1413             :     {
    1414        1059 :       if (u_ppo(e, 2*d) == 1)
    1415             :       { /* early abort: check whether f(3) = Phi_2d(3) mod p */
    1416        1052 :         ulong F3 = Rg_to_Fl(polcyclo_eval(2*d, _3), p);
    1417        1052 :         if (F3 == f3 && ZX_equal(f, polcyclo(2*d, varn(f))))
    1418         567 :           return gc_long(av, 2* d * e);
    1419             :       }
    1420             :     }
    1421             :   }
    1422        6751 :   return gc_long(av, 0);
    1423             : }
    1424             : 
    1425             : long
    1426        1029 : poliscycloprod(GEN f)
    1427             : {
    1428        1029 :   pari_sp av = avma;
    1429        1029 :   long i, d = degpol(f);
    1430        1029 :   if (typ(f) != t_POL) pari_err_TYPE("poliscycloprod",f);
    1431        1029 :   if (!RgX_is_ZX(f)) return 0;
    1432        1029 :   if (!ZX_is_monic(f) || !is_pm1(constant_coeff(f))) return 0;
    1433        1029 :   if (d < 2) return (d == 1);
    1434        1022 :   if ( degpol(ZX_gcd_all(f, ZX_deriv(f), &f)) )
    1435             :   {
    1436          14 :     d = degpol(f);
    1437          14 :     if (d == 1) return 1;
    1438             :   }
    1439        1015 :   f = BD(f); if (!f) return 0;
    1440        3619 :   for (i = lg(f)-1; i; i--) d -= degpol(gel(f,i));
    1441        1015 :   return gc_long(av, d == 0);
    1442             : }

Generated by: LCOV version 1.16