Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21730-50db3b0) Lines: 1520 1637 92.9 %
Date: 2018-01-19 06:18:47 Functions: 134 141 95.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /***********************************************************************/
      18             : /**                                                                   **/
      19             : /**               Factorisation over finite field                     **/
      20             : /**                                                                   **/
      21             : /***********************************************************************/
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : /* Check types and replace F by a monic normalized FpX having the same roots
      29             :  * Don't bother to make constant polynomials monic */
      30             : static void
      31         350 : factmod_init(GEN *F, GEN p)
      32             : {
      33         350 :   if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
      34         350 :   if (signe(p) < 0) pari_err_PRIME("factmod",p);
      35         350 :   if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
      36         350 :   if (lgefint(p) == 3)
      37             :   {
      38         224 :     ulong pp = p[2];
      39         224 :     if (pp < 2) pari_err_PRIME("factmod", p);
      40         224 :     *F = RgX_to_Flx(*F, pp);
      41         224 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42             :   }
      43             :   else
      44             :   {
      45         126 :     *F = RgX_to_FpX(*F, p);
      46         126 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      47             :   }
      48         350 : }
      49             : /* as above, assume p prime and *F a ZX */
      50             : static void
      51      933397 : ZX_factmod_init(GEN *F, GEN p)
      52             : {
      53      933397 :   if (lgefint(p) == 3)
      54             :   {
      55      925968 :     ulong pp = p[2];
      56      925968 :     *F = ZX_to_Flx(*F, pp);
      57      925968 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      58             :   }
      59             :   else
      60             :   {
      61        7429 :     *F = FpX_red(*F, p);
      62        7429 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      63             :   }
      64      933397 : }
      65             : 
      66             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      67             : static GEN
      68          42 : all_roots_mod_p(ulong p, int not_0)
      69             : {
      70             :   GEN r;
      71             :   ulong i;
      72          42 :   if (not_0) {
      73          28 :     r = cgetg(p, t_VECSMALL);
      74          28 :     for (i = 1; i < p; i++) r[i] = i;
      75             :   } else {
      76          14 :     r = cgetg(p+1, t_VECSMALL);
      77          14 :     for (i = 0; i < p; i++) r[i+1] = i;
      78             :   }
      79          42 :   return r;
      80             : }
      81             : 
      82             : /* X^n - 1 */
      83             : static GEN
      84          13 : Flx_Xnm1(long sv, long n, ulong p)
      85             : {
      86          13 :   GEN t = cgetg(n+3, t_VECSMALL);
      87             :   long i;
      88          13 :   t[1] = sv;
      89          13 :   t[2] = p - 1;
      90          13 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      91          13 :   t[i] = 1; return t;
      92             : }
      93             : /* X^n + 1 */
      94             : static GEN
      95           6 : Flx_Xn1(long sv, long n, ulong p)
      96             : {
      97           6 :   GEN t = cgetg(n+3, t_VECSMALL);
      98             :   long i;
      99             :   (void) p;
     100           6 :   t[1] = sv;
     101           6 :   t[2] = 1;
     102           6 :   for (i = 3; i <= n+1; i++) t[i] = 0;
     103           6 :   t[i] = 1; return t;
     104             : }
     105             : 
     106             : static ulong
     107           6 : Fl_nonsquare(ulong p)
     108             : {
     109           6 :   long k = 2;
     110           6 :   for (;; k++)
     111             :   {
     112          12 :     long i = krouu(k, p);
     113          12 :     if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
     114          18 :     if (i < 0) return k;
     115           6 :   }
     116             : }
     117             : 
     118             : /* f monic Flx, f(0) != 0. Return a monic squarefree g with the same
     119             :  * roots as f */
     120             : static GEN
     121          14 : Flx_cut_out_roots(GEN f, ulong p)
     122             : {
     123          14 :   GEN g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
     124          14 :   if (g != f && degpol(g) >= 0) {
     125           7 :     (void)Flx_valrem(g, &g); /* reduction may introduce 0 root */
     126           7 :     g = Flx_gcd(g, Flx_Xnm1(g[1], p-1, p), p);
     127           7 :     g = Flx_normalize(g, p);
     128             :   }
     129          14 :   return g;
     130             : }
     131             : 
     132             : /* by checking f(0..p-1) */
     133             : GEN
     134          14 : Flx_roots_naive(GEN f, ulong p)
     135             : {
     136          14 :   long d, n = 0;
     137          14 :   ulong s = 1UL, r;
     138          14 :   GEN q, y = cgetg(degpol(f) + 1, t_VECSMALL);
     139          14 :   pari_sp av2, av = avma;
     140             : 
     141          14 :   if (Flx_valrem(f, &f)) y[++n] = 0;
     142          14 :   f = Flx_cut_out_roots(f, p);
     143          14 :   d = degpol(f);
     144          14 :   if (d < 0) return all_roots_mod_p(p, n == 0);
     145          14 :   av2 = avma;
     146         308 :   while (d > 1) /* d = current degree of f */
     147             :   {
     148         287 :     q = Flx_div_by_X_x(f, s, p, &r); /* TODO: FFT-type multi-evaluation */
     149         287 :     if (r) avma = av2; else { y[++n] = s; d--; f = q; av2 = avma; }
     150         287 :     if (++s == p) break;
     151             :   }
     152          14 :   if (d == 1)
     153             :   { /* -f[2]/f[3], root of deg 1 polynomial */
     154           7 :     r = Fl_mul(p - Fl_inv(f[3], p), f[2], p);
     155           7 :     if (r >= s) y[++n] = r; /* otherwise double root */
     156             :   }
     157          14 :   avma = av; fixlg(y, n+1); return y;
     158             : }
     159             : static GEN
     160        8659 : Flx_root_mod_2(GEN f)
     161             : {
     162        8659 :   int z1, z0 = !(f[2] & 1);
     163             :   long i,n;
     164             :   GEN y;
     165             : 
     166        8659 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     167        8659 :   z1 = n & 1;
     168        8659 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     169        8659 :   if (z0) y[i++] = 0;
     170        8659 :   if (z1) y[i  ] = 1;
     171        8659 :   return y;
     172             : }
     173             : static ulong
     174          49 : Flx_oneroot_mod_2(GEN f)
     175             : {
     176             :   long i,n;
     177          49 :   if (!(f[2] & 1)) return 0;
     178          49 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     179          49 :   if (n & 1) return 1;
     180          28 :   return 2;
     181             : }
     182             : 
     183             : static GEN FpX_roots_i(GEN f, GEN p);
     184             : static GEN Flx_roots_i(GEN f, ulong p);
     185             : static GEN FpX_Berlekamp_i(GEN f, GEN pp, long flag);
     186             : 
     187             : static int
     188     2869780 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     189             : 
     190             : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
     191             :  * pp is a small prime.
     192             :  * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
     193             :  * assume that f is an FpX, pp a prime and return t_INTs */
     194             : static GEN
     195       70241 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
     196             : {
     197       70241 :   pari_sp av = avma;
     198             :   GEN y;
     199       70241 :   if (gpwrap)
     200          84 :     factmod_init(&f, pp);
     201             :   else
     202       70157 :     ZX_factmod_init(&f, pp);
     203       70241 :   switch(lg(f))
     204             :   {
     205          14 :     case 2: pari_err_ROOTS0("rootmod");
     206          49 :     case 3: avma = av; return cgetg(1,t_COL);
     207             :   }
     208       70178 :   if (typ(f) == t_VECSMALL)
     209             :   {
     210       67324 :     ulong p = pp[2];
     211       67324 :     if (p == 2)
     212        8659 :       y = Flx_root_mod_2(f);
     213             :     else
     214             :     {
     215       58665 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     216       58665 :       y = Roots(f, p);
     217             :     }
     218       67317 :     y = Flc_to_ZC(y);
     219             :   }
     220             :   else
     221        2854 :     y = FpX_roots_i(f, pp);
     222       70164 :   if (gpwrap) y = FpC_to_mod(y, pp);
     223       70164 :   return gerepileupto(av, y);
     224             : }
     225             : /* assume that f is a ZX an pp a prime */
     226             : GEN
     227       70157 : FpX_roots(GEN f, GEN pp)
     228       70157 : { return rootmod_aux(f, pp, Flx_roots_i, 0); }
     229             : /* no assumptions on f and pp */
     230             : GEN
     231          14 : rootmod2(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_naive, 1); }
     232             : GEN
     233          70 : rootmod(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_i, 1); }
     234             : GEN
     235          84 : rootmod0(GEN f, GEN p, long flag)
     236             : {
     237          84 :   switch(flag)
     238             :   {
     239          70 :     case 0: return rootmod(f,p);
     240          14 :     case 1: return rootmod2(f,p);
     241           0 :     default: pari_err_FLAG("polrootsmod");
     242             :   }
     243             :   return NULL; /* LCOV_EXCL_LINE */
     244             : }
     245             : 
     246             : /* assume x reduced mod p > 2, monic. */
     247             : static int
     248          21 : FpX_quad_factortype(GEN x, GEN p)
     249             : {
     250          21 :   GEN b = gel(x,3), c = gel(x,2);
     251          21 :   GEN D = subii(sqri(b), shifti(c,2));
     252          21 :   return kronecker(D,p);
     253             : }
     254             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     255             : GEN
     256        6916 : FpX_quad_root(GEN x, GEN p, int unknown)
     257             : {
     258        6916 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     259             : 
     260        6916 :   if (absequaliu(p, 2)) {
     261           0 :     if (!signe(b)) return c;
     262           0 :     return signe(c)? NULL: gen_1;
     263             :   }
     264        6916 :   D = subii(sqri(b), shifti(c,2));
     265        6916 :   D = remii(D,p);
     266        6916 :   if (unknown && kronecker(D,p) == -1) return NULL;
     267             : 
     268        6823 :   s = Fp_sqrt(D,p);
     269             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     270        6823 :   if (!s) return NULL;
     271        6813 :   return Fp_halve(Fp_sub(s,b, p), p);
     272             : }
     273             : static GEN
     274        2922 : FpX_otherroot(GEN x, GEN r, GEN p)
     275        2922 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     276             : 
     277             : /* disc(x^2+bx+c) = b^2 - 4c */
     278             : static ulong
     279    21210062 : Fl_disc_bc(ulong b, ulong c, ulong p)
     280    21210062 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     281             : /* p > 2 */
     282             : static ulong
     283    20495691 : Flx_quad_root(GEN x, ulong p, int unknown)
     284             : {
     285    20495691 :   ulong s, b = x[3], c = x[2];
     286    20495691 :   ulong D = Fl_disc_bc(b, c, p);
     287    20490841 :   if (unknown && krouu(D,p) == -1) return p;
     288    13688262 :   s = Fl_sqrt(D,p);
     289    13677589 :   if (s==~0UL) return p;
     290    13677564 :   return Fl_halve(Fl_sub(s,b, p), p);
     291             : }
     292             : static ulong
     293    12256653 : Flx_otherroot(GEN x, ulong r, ulong p)
     294    12256653 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     295             : 
     296             : 
     297             : /* 'todo' contains the list of factors to be split.
     298             :  * 'done' the list of finished factors, no longer touched */
     299             : struct split_t { GEN todo, done; };
     300             : static void
     301     4993883 : split_init(struct split_t *S, long max)
     302             : {
     303     4993883 :   S->todo = vectrunc_init(max);
     304     4993999 :   S->done = vectrunc_init(max);
     305     4993893 : }
     306             : #if 0
     307             : /* move todo[i] to done */
     308             : static void
     309             : split_convert(struct split_t *S, long i)
     310             : {
     311             :   long n = lg(S->todo)-1;
     312             :   vectrunc_append(S->done, gel(S->todo,i));
     313             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     314             :   setlg(S->todo, n);
     315             : }
     316             : #endif
     317             : /* append t to todo */
     318             : static void
     319     5232526 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     320             : /* delete todo[i], add t to done */
     321             : static void
     322     5232602 : split_moveto_done(struct split_t *S, long i, GEN t)
     323             : {
     324     5232602 :   long n = lg(S->todo)-1;
     325     5232602 :   vectrunc_append(S->done, t);
     326     5233037 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     327     5233037 :   setlg(S->todo, n);
     328             : 
     329     5233187 : }
     330             : /* append t to done */
     331             : static void
     332      383037 : split_add_done(struct split_t *S, GEN t)
     333      383037 : { vectrunc_append(S->done, t); }
     334             : /* split todo[i] into a and b */
     335             : static void
     336      331536 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     337             : {
     338      331536 :   gel(S->todo, i) = a;
     339      331536 :   split_add(S, b);
     340      331544 : }
     341             : /* split todo[i] into a and b, moved to done */
     342             : static void
     343      364010 : split_done(struct split_t *S, long i, GEN a, GEN b)
     344             : {
     345      364010 :   split_moveto_done(S, i, a);
     346      364017 :   split_add_done(S, b);
     347      364020 : }
     348             : 
     349             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     350             : static GEN
     351        2854 : FpX_roots_i(GEN f, GEN p)
     352             : {
     353             :   GEN pol, pol0, a, q;
     354             :   struct split_t S;
     355             : 
     356        2854 :   split_init(&S, lg(f)-1);
     357        2854 :   settyp(S.done, t_COL);
     358        2854 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     359        2854 :   switch(degpol(f))
     360             :   {
     361           7 :     case 0: return ZC_copy(S.done);
     362          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     363             :     case 2: {
     364        1638 :       GEN s, r = FpX_quad_root(f, p, 1);
     365        1638 :       if (r) {
     366        1638 :         split_add_done(&S, r);
     367        1638 :         s = FpX_otherroot(f,r, p);
     368             :         /* f not known to be square free yet */
     369        1638 :         if (!equalii(r, s)) split_add_done(&S, s);
     370             :       }
     371        1638 :       return sort(S.done);
     372             :     }
     373             :   }
     374             : 
     375        1195 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     376        1195 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     377        1195 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     378        1195 :   a = FpX_gcd(f,a, p);
     379        1195 :   if (!degpol(a)) return ZC_copy(S.done);
     380        1195 :   split_add(&S, FpX_normalize(a,p));
     381             : 
     382        1195 :   q = shifti(p,-1);
     383        1195 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     384        1195 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     385        2502 :   for (pol0[2] = 1;; pol0[2]++)
     386             :   {
     387        2502 :     long j, l = lg(S.todo);
     388        2502 :     if (l == 1) return sort(S.done);
     389        1314 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     390        2733 :     for (j = 1; j < l; j++)
     391             :     {
     392        1426 :       GEN c = gel(S.todo,j);
     393        1426 :       switch(degpol(c))
     394             :       { /* convert linear and quadratics to roots, try to split the rest */
     395             :         case 1:
     396          77 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     397          77 :           j--; l--; break;
     398             :         case 2: {
     399        1230 :           GEN r = FpX_quad_root(c, p, 0), s;
     400        1230 :           if (!r) pari_err_PRIME("polrootsmod",p);
     401        1223 :           s = FpX_otherroot(c,r, p);
     402        1223 :           split_done(&S, j, r, s);
     403        1223 :           j--; l--; break;
     404             :         }
     405             :         default: {
     406             :           /* b = pol^(p-1)/2 - 1 */
     407         119 :           GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
     408             :           long db;
     409         119 :           b = FpX_gcd(c,b, p); db = degpol(b);
     410         119 :           if (db && db < degpol(c))
     411             :           {
     412         112 :             b = FpX_normalize(b, p);
     413         112 :             c = FpX_div(c,b, p);
     414         112 :             split_todo(&S, j, b, c);
     415             :           }
     416             :         }
     417             :       }
     418             :     }
     419        1307 :   }
     420             : }
     421             : 
     422             : /* Assume f is normalized */
     423             : static ulong
     424      171482 : Flx_cubic_root(GEN ff, ulong p)
     425             : {
     426      171482 :   GEN f = Flx_normalize(ff,p);
     427      171482 :   ulong pi = get_Fl_red(p);
     428      171480 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     429      171480 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     430      171484 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     431      171480 :   ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
     432      171481 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     433      171484 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     434      171484 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     435      171483 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     436      171484 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     437      171484 :   if (s!=~0UL)
     438             :   {
     439      108013 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     440      108013 :     if (p%3==2) /* 1 solutions */
     441       18777 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     442             :     else
     443             :     {
     444       89236 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     445       89236 :       if (vS1==~0UL) return p; /*0 solutions*/
     446             :       /*3 solutions*/
     447             :     }
     448       78161 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     449       78161 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     450             :   }
     451             :   else
     452             :   {
     453       63471 :     pari_sp av = avma;
     454       63471 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     455       63471 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     456             :     ulong Sa;
     457       63472 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     458       63472 :     Sa = vS1[1];
     459       63472 :     if (p%3==1) /*1 solutions*/
     460             :     {
     461       23327 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     462       23327 :       if (Fa!=P)
     463       15458 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     464             :     }
     465       63472 :     avma = av;
     466       63472 :     return Fl_sub(Fl_double(Sa,p),t,p);
     467             :   }
     468             : }
     469             : 
     470             : /* assume p > 2 prime */
     471             : static ulong
     472     3107102 : Flx_oneroot_i(GEN f, ulong p, long fl)
     473             : {
     474             :   GEN pol, a;
     475             :   ulong q;
     476             :   long da;
     477             : 
     478     3107102 :   if (Flx_val(f)) return 0;
     479     3106506 :   switch(degpol(f))
     480             :   {
     481       10986 :     case 1: return Fl_neg(f[2], p);
     482     2661436 :     case 2: return Flx_quad_root(f, p, 1);
     483      155930 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     484             :   }
     485             : 
     486      278084 :   if (!fl)
     487             :   {
     488      247501 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     489      247501 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     490      247501 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     491      247501 :     a = Flx_gcd(f,a, p);
     492       30583 :   } else a = f;
     493      278084 :   da = degpol(a);
     494      278085 :   if (!da) return p;
     495      198621 :   a = Flx_normalize(a,p);
     496             : 
     497      198685 :   q = p >> 1;
     498      198685 :   pol = polx_Flx(f[1]);
     499      297216 :   for(pol[2] = 1;; pol[2]++)
     500             :   {
     501      297216 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     502      297152 :     switch(da)
     503             :     {
     504      123819 :       case 1: return Fl_neg(a[2], p);
     505       59385 :       case 2: return Flx_quad_root(a, p, 0);
     506       15552 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     507             :       default: {
     508       98396 :         GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
     509             :         long db;
     510       98570 :         b = Flx_gcd(a,b, p); db = degpol(b);
     511       98563 :         if (db && db < da)
     512             :         {
     513       93651 :           b = Flx_normalize(b, p);
     514       93656 :           if (db <= (da >> 1)) {
     515       57012 :             a = b;
     516       57012 :             da = db;
     517             :           } else {
     518       36644 :             a = Flx_div(a,b, p);
     519       36641 :             da -= db;
     520             :           }
     521             :         }
     522             :       }
     523             :     }
     524       98565 :   }
     525             : }
     526             : 
     527             : /* assume p > 2 prime */
     528             : static GEN
     529        4284 : FpX_oneroot_i(GEN f, GEN p)
     530             : {
     531             :   GEN pol, pol0, a, q;
     532             :   long da;
     533             : 
     534        4284 :   if (ZX_val(f)) return gen_0;
     535        4130 :   switch(degpol(f))
     536             :   {
     537         330 :     case 1: return subii(p, gel(f,2));
     538        3730 :     case 2: return FpX_quad_root(f, p, 1);
     539             :   }
     540             : 
     541          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     542          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     543          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     544          70 :   a = FpX_gcd(f,a, p);
     545          70 :   da = degpol(a);
     546          70 :   if (!da) return NULL;
     547          70 :   a = FpX_normalize(a,p);
     548             : 
     549          70 :   q = shifti(p,-1);
     550          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     551          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     552         224 :   for (pol0[2]=1; ; pol0[2]++)
     553             :   {
     554         224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     555         224 :     switch(da)
     556             :     {
     557          42 :       case 1: return subii(p, gel(a,2));
     558          28 :       case 2: return FpX_quad_root(a, p, 0);
     559             :       default: {
     560         154 :         GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, a,p), gen_1, p);
     561             :         long db;
     562         154 :         b = FpX_gcd(a,b, p); db = degpol(b);
     563         154 :         if (db && db < da)
     564             :         {
     565         147 :           b = FpX_normalize(b, p);
     566         147 :           if (db <= (da >> 1)) {
     567         105 :             a = b;
     568         105 :             da = db;
     569             :           } else {
     570          42 :             a = FpX_div(a,b, p);
     571          42 :             da -= db;
     572             :           }
     573             :         }
     574             :       }
     575             :     }
     576         154 :   }
     577             : }
     578             : 
     579             : ulong
     580     3068468 : Flx_oneroot(GEN f, ulong p)
     581             : {
     582     3068468 :   pari_sp av = avma;
     583             :   ulong r;
     584     3068468 :   switch(lg(f))
     585             :   {
     586           0 :     case 2: return 0;
     587           0 :     case 3: avma = av; return p;
     588             :   }
     589     3068468 :   if (p == 2) return Flx_oneroot_mod_2(f);
     590     3068468 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     591     3068468 :   avma = av; return r;
     592             : }
     593             : 
     594             : ulong
     595       31435 : Flx_oneroot_split(GEN f, ulong p)
     596             : {
     597       31435 :   pari_sp av = avma;
     598             :   ulong r;
     599       31435 :   switch(lg(f))
     600             :   {
     601           0 :     case 2: return 0;
     602           0 :     case 3: avma = av; return p;
     603             :   }
     604       31435 :   if (p == 2) return Flx_oneroot_mod_2(f);
     605       31435 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     606       31541 :   avma = av; return r;
     607             : }
     608             : 
     609             : /* assume that p is prime */
     610             : GEN
     611       11543 : FpX_oneroot(GEN f, GEN pp) {
     612       11543 :   pari_sp av = avma;
     613       11543 :   ZX_factmod_init(&f, pp);
     614       11543 :   switch(lg(f))
     615             :   {
     616           0 :     case 2: avma = av; return gen_0;
     617           0 :     case 3: avma = av; return NULL;
     618             :   }
     619       11543 :   if (typ(f) == t_VECSMALL)
     620             :   {
     621        7259 :     ulong r, p = pp[2];
     622        7259 :     if (p == 2)
     623          49 :       r = Flx_oneroot_mod_2(f);
     624             :     else
     625        7210 :       r = Flx_oneroot_i(f, p, 0);
     626        7259 :     avma = av;
     627        7259 :     return (r == p)? NULL: utoi(r);
     628             :   }
     629        4284 :   f = FpX_oneroot_i(f, pp);
     630        4284 :   if (!f) { avma = av; return NULL; }
     631        4284 :   return gerepileuptoint(av, f);
     632             : }
     633             : 
     634             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     635             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     636             : /* returned in n                                                          */
     637             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     638             : /* twice as big and decrement until it divides p-1.                       */
     639             : static GEN
     640         238 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     641             : {
     642         238 :    pari_sp ltop = avma;
     643             :    GEN pm, factn, power, base, zeta;
     644             :    long n;
     645             : 
     646         238 :    pm = subis (p, 1ul);
     647         238 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     648         238 :    factn = Z_factor(stoi(n));
     649         238 :    power = diviuexact (pm, n);
     650         238 :    base = gen_1;
     651             :    do {
     652         392 :       base = addis (base, 1l);
     653         392 :       zeta = Fp_pow (base, power, p);
     654             :    }
     655         392 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     656         238 :    *pt_n = n;
     657         238 :    return gerepileuptoint (ltop, zeta);
     658             : }
     659             : 
     660             : GEN
     661         665 : FpX_oneroot_split(GEN fact, GEN p)
     662             : {
     663         665 :   pari_sp av = avma;
     664             :   long n, deg_f, i, dmin;
     665             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     666         665 :   fact = FpX_normalize(fact, p);
     667         665 :   deg_f = degpol(fact);
     668         665 :   if (deg_f<=2) return FpX_oneroot(fact, p);
     669         238 :   minfactor = fact; /* factor of minimal degree found so far */
     670         238 :   dmin = degpol(minfactor);
     671         238 :   prim = good_root_of_unity(p, deg_f, 1, &n);
     672         238 :   expo = diviuexact(subiu(p, 1), n);
     673         238 :   xplusa = pol_x(varn(fact));
     674         238 :   zeta = gen_1;
     675         945 :   while (dmin != 1)
     676             :   {
     677             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     678             :     /* zeta varies over the roots of unity in F_p                         */
     679         469 :     fact = minfactor; deg_f = dmin;
     680             :     /* update X+a, avoid a=0 */
     681         469 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     682         469 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     683         833 :     for (i = 0; i < n; i++)
     684             :     {
     685         735 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     686         735 :       long dtmp = degpol(tmp);
     687         735 :       if (dtmp > 0 && dtmp < deg_f)
     688             :       {
     689         427 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     690         427 :         if (dtmp < dmin)
     691             :         {
     692         427 :           minfactor = FpX_normalize (tmp, p);
     693         427 :           dmin = dtmp;
     694         427 :           if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
     695             :             /* stop early to avoid too many gcds */
     696             :             break;
     697             :         }
     698             :       }
     699         364 :       zeta = Fp_mul (zeta, prim, p);
     700             :     }
     701             :   }
     702         238 :   return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                                                                 */
     707             : /*                     FACTORISATION MODULO p                      */
     708             : /*                                                                 */
     709             : /*******************************************************************/
     710             : 
     711             : /* Functions giving information on the factorisation. */
     712             : 
     713             : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
     714             : static GEN
     715         171 : FpX_Berlekamp_ker(GEN u, GEN p)
     716             : {
     717         171 :   pari_sp ltop=avma;
     718         171 :   long j,N = degpol(u);
     719         171 :   GEN Q  = FpX_matFrobenius(u, p);
     720        3025 :   for (j=1; j<=N; j++)
     721        2854 :     gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
     722         171 :   return gerepileupto(ltop, FpM_ker(Q,p));
     723             : }
     724             : 
     725             : static GEN
     726       19455 : F2x_Berlekamp_ker(GEN u)
     727             : {
     728       19455 :   pari_sp ltop=avma;
     729       19455 :   long j,N = F2x_degree(u);
     730             :   GEN Q;
     731             :   pari_timer T;
     732       19455 :   timer_start(&T);
     733       19455 :   Q = F2x_matFrobenius(u);
     734      116502 :   for (j=1; j<=N; j++)
     735       97047 :     F2m_flip(Q,j,j);
     736       19455 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     737       19455 :   Q = F2m_ker_sp(Q,0);
     738       19455 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     739       19455 :   return gerepileupto(ltop,Q);
     740             : }
     741             : 
     742             : static GEN
     743      281150 : Flx_Berlekamp_ker(GEN u, ulong l)
     744             : {
     745      281150 :   pari_sp ltop=avma;
     746      281150 :   long j,N = degpol(u);
     747             :   GEN Q;
     748             :   pari_timer T;
     749      281150 :   timer_start(&T);
     750      281150 :   Q  = Flx_matFrobenius(u, l);
     751     1643958 :   for (j=1; j<=N; j++)
     752     1362808 :     coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
     753      281150 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     754      281150 :   Q = Flm_ker_sp(Q,l,0);
     755      281150 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     756      281150 :   return gerepileupto(ltop,Q);
     757             : }
     758             : 
     759             : /* product of terms of degree 1 in factorization of f */
     760             : GEN
     761      142528 : FpX_split_part(GEN f, GEN p)
     762             : {
     763      142528 :   long n = degpol(f);
     764      142528 :   GEN z, X = pol_x(varn(f));
     765      142528 :   if (n <= 1) return f;
     766      142227 :   f = FpX_red(f, p);
     767      142227 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     768      142227 :   return FpX_gcd(z,f,p);
     769             : }
     770             : 
     771             : /* Compute the number of roots in Fp without counting multiplicity
     772             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     773             : long
     774       99085 : FpX_nbroots(GEN f, GEN p)
     775             : {
     776       99085 :   pari_sp av = avma;
     777       99085 :   GEN z = FpX_split_part(f, p);
     778       99085 :   avma = av; return degpol(z);
     779             : }
     780             : 
     781             : int
     782           0 : FpX_is_totally_split(GEN f, GEN p)
     783             : {
     784           0 :   long n=degpol(f);
     785           0 :   pari_sp av = avma;
     786           0 :   if (n <= 1) return 1;
     787           0 :   if (abscmpui(n, p) > 0) return 0;
     788           0 :   f = FpX_red(f, p);
     789           0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     790             : }
     791             : 
     792             : long
     793     2402593 : Flx_nbroots(GEN f, ulong p)
     794             : {
     795     2402593 :   long n = degpol(f);
     796     2402593 :   pari_sp av = avma;
     797             :   GEN z;
     798     2402593 :   if (n <= 1) return n;
     799     2400486 :   if (n == 2)
     800             :   {
     801             :     ulong D;
     802       11508 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     803       10570 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     804       10570 :     return 1 + krouu(D,p);
     805             :   }
     806     2388978 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     807     2388978 :   z = Flx_gcd(z, f, p);
     808     2388978 :   avma = av; return degpol(z);
     809             : }
     810             : 
     811             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     812             : static GEN
     813        4242 : FpX_ddf(GEN T, GEN XP, GEN p)
     814             : {
     815        4242 :   pari_sp av = avma;
     816             :   GEN b, g, h, F, f, Tr, xq;
     817             :   long i, j, n, v;
     818             :   long B, l, m;
     819             :   pari_timer ti;
     820        4242 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     821        4242 :   if (n == 0) return cgetg(1, t_VEC);
     822        4242 :   if (n == 1) return mkvec(get_FpX_mod(T));
     823        4235 :   B = n/2;
     824        4235 :   l = usqrt(B);
     825        4235 :   m = (B+l-1)/l;
     826        4235 :   T = FpX_get_red(T, p);
     827        4235 :   b = cgetg(l+2, t_VEC);
     828        4235 :   gel(b, 1) = pol_x(v);
     829        4235 :   gel(b, 2) = XP;
     830        4235 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     831        4235 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     832        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
     833       10976 :   for (i = 3; i <= l+1; i++)
     834        6741 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     835        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
     836        4235 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     837        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
     838        4235 :   g = cgetg(m+1, t_VEC);
     839        4235 :   gel(g, 1) = gel(xq, 2);
     840       15379 :   for(i = 2; i <= m; i++)
     841       11144 :     gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     842        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
     843        4235 :   h = cgetg(m+1, t_VEC);
     844       19614 :   for (j = 1; j <= m; j++)
     845             :   {
     846       15379 :     pari_sp av = avma;
     847       15379 :     GEN gj = gel(g, j);
     848       15379 :     GEN e = FpX_sub(gj, gel(b, 1), p);
     849       44079 :     for (i = 2; i <= l; i++)
     850       28700 :       e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
     851       15379 :     gel(h, j) = gerepileupto(av, e);
     852             :   }
     853        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
     854        4235 :   Tr = get_FpX_mod(T);
     855        4235 :   F = cgetg(m+1, t_VEC);
     856       19614 :   for (j = 1; j <= m; j++)
     857             :   {
     858       15379 :     gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
     859       15379 :     Tr = FpX_div(Tr, gel(F,j), p);
     860             :   }
     861        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
     862        4235 :   f = const_vec(n, pol_1(v));
     863       19614 :   for (j = 1; j <= m; j++)
     864             :   {
     865       15379 :     GEN e = gel(F, j);
     866       17416 :     for (i=l-1; i >= 0; i--)
     867             :     {
     868       17416 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     869       17416 :       if (degpol(u))
     870             :       {
     871        2779 :         gel(f, l*j-i) = u;
     872        2779 :         e = FpX_div(e, u, p);
     873             :       }
     874       17416 :       if (!degpol(e)) break;
     875             :     }
     876             :   }
     877        4235 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
     878        4235 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     879        4235 :   return gerepilecopy(av, f);
     880             : }
     881             : 
     882             : static void
     883           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     884             : {
     885           0 :   long n = degpol(Tp), r = n/d;
     886             :   GEN T, f, ff;
     887             :   GEN p2;
     888           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     889           0 :   p2 = shifti(p,-1);
     890           0 :   T = FpX_get_red(Tp, p);
     891           0 :   XP = FpX_rem(XP, T, p);
     892             :   while (1)
     893             :   {
     894           0 :     pari_sp btop = avma;
     895             :     long i;
     896           0 :     GEN g = random_FpX(n, varn(Tp), p);
     897           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     898           0 :     if (signe(t) == 0) continue;
     899           0 :     for(i=1; i<=10; i++)
     900             :     {
     901           0 :       pari_sp btop2 = avma;
     902           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     903           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     904           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     905           0 :       avma = btop2;
     906             :     }
     907           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     908           0 :     avma = btop;
     909           0 :   }
     910           0 :   f = FpX_normalize(f, p);
     911           0 :   ff = FpX_div(Tp, f ,p);
     912           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     913           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     914             : }
     915             : 
     916             : static void
     917         210 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     918             : {
     919             :   pari_sp av;
     920         210 :   GEN Tp = get_FpX_mod(T);
     921         210 :   long n = degpol(hp), vT = varn(Tp);
     922             :   GEN u1, u2, f1, f2;
     923             :   GEN R, h;
     924         210 :   h = FpX_get_red(hp, p);
     925         210 :   t = FpX_rem(t, T, p);
     926         210 :   av = avma;
     927             :   do
     928             :   {
     929         343 :     avma = av;
     930         343 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     931         343 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     932         343 :   } while (degpol(u1)==0 || degpol(u1)==n);
     933         210 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     934         210 :   f1 = FpX_normalize(f1, p);
     935         210 :   u2 = FpX_div(hp, u1, p);
     936         210 :   f2 = FpX_div(Tp, f1, p);
     937         210 :   if (degpol(u1)==1)
     938         105 :     gel(V, idx) = f1;
     939             :   else
     940         105 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     941         210 :   idx += degpol(f1)/d;
     942         210 :   if (degpol(u2)==1)
     943         119 :     gel(V, idx) = f2;
     944             :   else
     945          91 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     946         210 : }
     947             : 
     948             : static void
     949          14 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     950             : {
     951          14 :   long n = degpol(Tp), r = n/d, vT = varn(Tp);
     952             :   GEN T, h, t;
     953             :   pari_timer ti;
     954          28 :   if (r==1) { gel(V, idx) = Tp; return; }
     955          14 :   T = FpX_get_red(Tp, p);
     956          14 :   XP = FpX_rem(XP, T, p);
     957          14 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     958             :   do
     959             :   {
     960          14 :     GEN g = random_FpX(n, vT, p);
     961          14 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     962          14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     963          14 :     h = FpXQ_minpoly(t, T, p);
     964          14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     965          14 :   } while (degpol(h) != r);
     966          14 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     967             : }
     968             : 
     969             : static GEN
     970          14 : FpX_factor_Shoup(GEN T, GEN p)
     971             : {
     972          14 :   long i, n, s = 0;
     973             :   GEN XP, D, V;
     974          14 :   long e = expi(p);
     975             :   pari_timer ti;
     976          14 :   n = get_FpX_degree(T);
     977          14 :   T = FpX_get_red(T, p);
     978          14 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     979          14 :   XP = FpX_Frobenius(T, p);
     980          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     981          14 :   D = FpX_ddf(T, XP, p);
     982          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     983         301 :   for (i = 1; i <= n; i++)
     984         287 :     s += degpol(gel(D,i))/i;
     985          14 :   V = cgetg(s+1, t_COL);
     986         301 :   for (i = 1, s = 1; i <= n; i++)
     987             :   {
     988         287 :     GEN Di = gel(D,i);
     989         287 :     long ni = degpol(Di), ri = ni/i;
     990         287 :     if (ni == 0) continue;
     991          21 :     Di = FpX_normalize(Di, p);
     992          21 :     if (ni == i) { gel(V, s++) = Di; continue; }
     993          14 :     if (ri <= e*expu(e))
     994          14 :       FpX_edf(Di, XP, i, p, V, s);
     995             :     else
     996           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     997          14 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     998          14 :     s += ri;
     999             :   }
    1000          14 :   return V;
    1001             : }
    1002             : 
    1003             : GEN
    1004      352004 : ddf_to_simplefact(GEN D, long n)
    1005             : {
    1006             :   GEN V;
    1007      352004 :   long i, s = 0, j = 1, k;
    1008     2340615 :   for (i = 1; i <= n; i++)
    1009     1988611 :     s += degpol(gel(D,i))/i;
    1010      352004 :   V = cgetg(s+1, t_VECSMALL);
    1011     2340615 :   for (i = 1; i <= n; i++)
    1012             :   {
    1013     1988611 :     long ni = degpol(gel(D,i)), ri = ni/i;
    1014     1988611 :     if (ni == 0) continue;
    1015     1352424 :     for (k = 1; k <= ri; k++)
    1016      905509 :       uel(V, j++) = i;
    1017             :   }
    1018      352004 :   return V;
    1019             : }
    1020             : 
    1021             : static GEN
    1022          14 : FpX_simplefact_Shoup(GEN T, GEN p)
    1023             : {
    1024             :   GEN XP, D;
    1025          14 :   long n = get_FpX_degree(T);
    1026             :   pari_timer ti;
    1027          14 :   T = FpX_get_red(T, p);
    1028          14 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1029          14 :   XP = FpX_Frobenius(T, p);
    1030          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1031          14 :   D = FpX_ddf(T, XP, p);
    1032          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1033          14 :   return ddf_to_simplefact(D, n);
    1034             : }
    1035             : 
    1036             : /* Yun algorithm: Assume p > degpol(T) */
    1037             : static GEN
    1038         453 : FpX_factor_Yun(GEN T, GEN p)
    1039             : {
    1040         453 :   long n = degpol(T);
    1041         453 :   long i = 1;
    1042         453 :   GEN d = FpX_deriv(T, p);
    1043             :   GEN a, b, c;
    1044         453 :   GEN V = cgetg(n+1,t_VEC);
    1045         453 :   a = FpX_gcd(T, d, p);
    1046         453 :   if (degpol(a) == 0) return mkvec(T);
    1047         208 :   b = FpX_div(T, a, p);
    1048             :   do
    1049             :   {
    1050        1830 :     c = FpX_div(d, a, p);
    1051        1830 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1052        1830 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1053        1830 :     gel(V, i++) = a;
    1054        1830 :     b = FpX_div(b, a, p);
    1055        1830 :   } while (degpol(b));
    1056         208 :   setlg(V, i); return V;
    1057             : }
    1058             : GEN
    1059         742 : FpX_factor_squarefree(GEN T, GEN p)
    1060             : {
    1061         742 :   if (abscmpiu(p, degpol(T)) <= 0)
    1062             :   {
    1063         525 :     ulong pp = (ulong)p[2];
    1064         525 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1065         525 :     return FlxV_to_ZXV(u);
    1066             :   }
    1067         217 :   return FpX_factor_Yun(T, p);
    1068             : }
    1069             : 
    1070             : long
    1071         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1072             : {
    1073         168 :   pari_sp av = avma;
    1074             :   GEN lc, F;
    1075         168 :   long i, l, n = degpol(f), v = varn(f);
    1076         168 :   if (n % k) return 0;
    1077         168 :   if (lgefint(p)==3)
    1078             :   {
    1079         126 :     ulong pp = p[2];
    1080         126 :     GEN fp = ZX_to_Flx(f, pp);
    1081         126 :     if (!Flx_ispower(fp, k, pp, pt_r))
    1082          21 :     { avma = av; return 0; }
    1083         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r));
    1084           0 :     else avma = av;
    1085         105 :     return 1;
    1086             :   }
    1087          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1088          42 :   if (!lc) { av = avma; return 0; }
    1089          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1090        1491 :   for(i=1; i <= l; i++)
    1091        1456 :     if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
    1092          35 :   if (pt_r)
    1093             :   {
    1094          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1095        1484 :     for (i=l; i>=1; i--)
    1096             :     {
    1097        1449 :       if (i%k) continue;
    1098         294 :       s = FpX_mul(s, gel(F,i), p);
    1099         294 :       r = FpX_mul(r, s, p);
    1100             :     }
    1101          35 :     *pt_r = gerepileupto(av, r);
    1102           0 :   } else av = avma;
    1103          35 :   return 1;
    1104             : }
    1105             : 
    1106             : /* F / E  a vector of factors / exponents of virtual length l
    1107             :  * (their real lg may be larger). Set their lg to j and return [F,E] */
    1108             : static GEN
    1109      418751 : FE_setlg(GEN F, GEN E, long l)
    1110             : {
    1111      418751 :   setlg(E,l);
    1112      418751 :   setlg(F,l); return mkvec2(F,E);
    1113             : }
    1114             : /* F / E  a vector of vectors of factors / exponents of virtual length l
    1115             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
    1116             : static GEN
    1117      400450 : FE_concat(GEN F, GEN E, long l)
    1118             : {
    1119      400450 :   setlg(E,l); E = shallowconcat1(E);
    1120      400450 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
    1121             : }
    1122             : 
    1123             : static GEN
    1124           7 : FpX_factor_Cantor(GEN T, GEN p)
    1125             : {
    1126           7 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1127           7 :   long i, j, l = lg(V);
    1128           7 :   F = cgetg(l, t_VEC);
    1129           7 :   E = cgetg(l, t_VEC);
    1130          21 :   for (i=1, j=1; i < l; i++)
    1131          14 :     if (degpol(gel(V,i)))
    1132             :     {
    1133          14 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1134          14 :       gel(F, j) = Fj;
    1135          14 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1136          14 :       j++;
    1137             :     }
    1138           7 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1139             : }
    1140             : 
    1141             : static GEN
    1142          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1143             : {
    1144             :   GEN E, F, V;
    1145             :   long i, j, l;
    1146          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p);
    1147          14 :   l = lg(V);
    1148          14 :   E = cgetg(l, t_VEC);
    1149          14 :   F = cgetg(l, t_VEC);
    1150          28 :   for (i=1, j=1; i < l; i++)
    1151          14 :     if (degpol(gel(V,i)))
    1152             :     {
    1153          14 :       GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
    1154          14 :       gel(F, j) = Fj;
    1155          14 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1156          14 :       j++;
    1157             :     }
    1158          14 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1159             : }
    1160             : 
    1161             : static int
    1162           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1163             : {
    1164           0 :   pari_sp av = avma;
    1165             :   pari_timer ti;
    1166             :   long n, d;
    1167           0 :   GEN T = get_FpX_mod(Tp);
    1168           0 :   GEN dT = FpX_deriv(T, p);
    1169             :   GEN XP, D;
    1170           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1171           0 :   n = get_FpX_degree(T);
    1172           0 :   T = FpX_get_red(Tp, p);
    1173           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1174           0 :   XP = FpX_Frobenius(T, p);
    1175           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1176           0 :   D = FpX_ddf(T, XP, p);
    1177           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1178           0 :   d = degpol(gel(D, n));
    1179           0 :   avma = av; return d==n;
    1180             : }
    1181             : 
    1182             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1183             : 
    1184             : /*Assume that p is large and odd*/
    1185             : static GEN
    1186          49 : FpX_factcantor_i(GEN f, GEN pp, long flag)
    1187             : {
    1188          49 :   long d = degpol(f);
    1189          49 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1190          21 :   switch(flag)
    1191             :   {
    1192           7 :     default: return FpX_factor_Cantor(f, pp);
    1193          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1194           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1195             :   }
    1196             : }
    1197             : 
    1198             : long
    1199       85715 : ddf_to_nbfact(GEN D)
    1200             : {
    1201       85715 :   long l = lg(D), i, s = 0;
    1202      793100 :   for(i = 1; i < l; i++)
    1203      707385 :     s += degpol(gel(D,i))/i;
    1204       85715 :   return s;
    1205             : }
    1206             : 
    1207             : long
    1208        4214 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1209             : {
    1210        4214 :   pari_sp av = avma;
    1211        4214 :   long s = ddf_to_nbfact(FpX_ddf(T, XP, p));
    1212        4214 :   avma = av; return s;
    1213             : }
    1214             : 
    1215             : long
    1216          28 : FpX_nbfact(GEN T, GEN p)
    1217             : {
    1218          28 :   pari_sp av = avma;
    1219          28 :   GEN XP = FpX_Frobenius(T, p);
    1220          28 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1221          28 :   avma = av; return n;
    1222             : }
    1223             : 
    1224             : /* p > 2 */
    1225             : static GEN
    1226           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1227             : {
    1228           7 :   switch(d)
    1229             :   {
    1230             :     case -1:
    1231           0 :     case 0: return NULL;
    1232           0 :     case 1: return gen_1;
    1233             :   }
    1234           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1235             : }
    1236             : /* p > 2 */
    1237             : static GEN
    1238          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1239             : {
    1240          14 :   switch(d)
    1241             :   {
    1242           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1243           0 :     case 0: return trivial_fact();
    1244           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1245             :   }
    1246          14 :   switch(FpX_quad_factortype(f, p)) {
    1247           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1248           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1249           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1250             :   }
    1251             : }
    1252             : 
    1253             : GEN
    1254          77 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1255             : GEN
    1256         952 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1257             : 
    1258             : /* Mod(0,p) * x, where x is f's main variable */
    1259             : static GEN
    1260          14 : Mod0pX(GEN f, GEN p)
    1261          14 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
    1262             : static GEN
    1263          14 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
    1264             : 
    1265             : /* not gerepile safe */
    1266             : static GEN
    1267         174 : FpX_factor_2(GEN f, GEN p, long d)
    1268             : {
    1269             :   GEN r, s, R, S;
    1270             :   long v;
    1271             :   int sgn;
    1272         174 :   switch(d)
    1273             :   {
    1274           0 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1275           2 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1276          24 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1277             :   }
    1278         148 :   r = FpX_quad_root(f, p, 1);
    1279         148 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1280          54 :   v = varn(f);
    1281          54 :   s = FpX_otherroot(f, r, p);
    1282          54 :   if (signe(r)) r = subii(p, r);
    1283          54 :   if (signe(s)) s = subii(p, s);
    1284          54 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1285          54 :   R = deg1pol_shallow(gen_1, r, v);
    1286          54 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1287          47 :   S = deg1pol_shallow(gen_1, s, v);
    1288          47 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1289             : }
    1290             : static GEN
    1291         195 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1292             : {
    1293         195 :   switch(flag) {
    1294           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1295          14 :     case 1: return FpX_degfact_2(f, p, d);
    1296         174 :     default: return FpX_factor_2(f, p, d);
    1297             :   }
    1298             : }
    1299             : 
    1300             : static int
    1301       58646 : F2x_quad_factortype(GEN x)
    1302       58646 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1303             : 
    1304             : static GEN
    1305           7 : F2x_is_irred_2(GEN f, long d)
    1306           7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1307             : 
    1308             : static GEN
    1309        7035 : F2x_degfact_2(GEN f, long d)
    1310             : {
    1311        7035 :   if (!d) return trivial_fact();
    1312        7035 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1313        6832 :   switch(F2x_quad_factortype(f)) {
    1314        2254 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1315        2219 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1316        2359 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1317             :   }
    1318             : }
    1319             : 
    1320             : static GEN
    1321       99946 : F2x_factor_2(GEN f, long d)
    1322             : {
    1323       99946 :   long v = f[1];
    1324       99946 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1325       99946 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1326       92246 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1327       45199 :   switch(F2x_quad_factortype(f))
    1328             :   {
    1329        9457 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1330       21364 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1331       14378 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1332             :   }
    1333             : }
    1334             : static GEN
    1335      106988 : F2x_factor_deg2(GEN f, long d, long flag)
    1336             : {
    1337      106988 :   switch(flag) {
    1338           7 :     case 2: return F2x_is_irred_2(f, d);
    1339        7035 :     case 1: return F2x_degfact_2(f, d);
    1340       99946 :     default: return F2x_factor_2(f, d);
    1341             :   }
    1342             : }
    1343             : 
    1344             : static void
    1345          12 : split_squares(struct split_t *S, GEN g, ulong p)
    1346             : {
    1347          12 :   ulong q = p >> 1;
    1348          12 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1349          12 :   if (degpol(a) < 0)
    1350             :   {
    1351             :     ulong i;
    1352           6 :     split_add_done(S, (GEN)1);
    1353           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1354             :   } else {
    1355           6 :     (void)Flx_valrem(a, &a);
    1356           6 :     if (degpol(a))
    1357             :     {
    1358           6 :       a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
    1359           6 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1360             :     }
    1361             :   }
    1362          12 : }
    1363             : static void
    1364          12 : split_nonsquares(struct split_t *S, GEN g, ulong p)
    1365             : {
    1366          12 :   ulong q = p >> 1;
    1367          12 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1368          12 :   if (degpol(a) < 0)
    1369             :   {
    1370           6 :     ulong i, z = Fl_nonsquare(p);
    1371           6 :     split_add_done(S, (GEN)z);
    1372           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1373             :   } else {
    1374           6 :     (void)Flx_valrem(a, &a);
    1375           6 :     if (degpol(a))
    1376             :     {
    1377           6 :       a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
    1378           6 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1379             :     }
    1380             :   }
    1381          12 : }
    1382             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1383             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1384             : static int
    1385     4990974 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1386             : {
    1387     4990974 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1388     4991113 :   long d = degpol(g);
    1389     4991047 :   if (d < 0) return 0;
    1390     4990922 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1391     4991034 :   if (!d) return 1;
    1392     4975445 :   if (p <= 1.4 * (ulong)d) {
    1393             :     /* small p; split further using x^((p-1)/2) +/- 1.
    1394             :      * 30% degree drop makes the extra gcd worth it. */
    1395          12 :     split_squares(S, g, p);
    1396          12 :     split_nonsquares(S, g, p);
    1397             :   } else { /* large p; use x^(p-1) - 1 directly */
    1398     4975433 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1399     4975004 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1400     4975004 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1401     4974703 :     g = Flx_gcd(g,a, p);
    1402     4974976 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1403             :   }
    1404     4975487 :   return 1;
    1405             : }
    1406             : 
    1407             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1408             : static GEN
    1409    21997329 : Flx_roots_i(GEN f, ulong p)
    1410             : {
    1411             :   GEN pol, g;
    1412    21997329 :   long v = Flx_valrem(f, &g);
    1413             :   ulong q;
    1414             :   struct split_t S;
    1415             : 
    1416             :   /* optimization: test for small degree first */
    1417    21998974 :   switch(degpol(g))
    1418             :   {
    1419             :     case 1: {
    1420       25369 :       ulong r = p - g[2];
    1421       25369 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1422             :     }
    1423             :     case 2: {
    1424    16980361 :       ulong r = Flx_quad_root(g, p, 1), s;
    1425    16993915 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1426    11623707 :       s = Flx_otherroot(g,r, p);
    1427    11632888 :       if (r < s)
    1428     2902879 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1429     8730009 :       else if (r > s)
    1430     8729848 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1431             :       else
    1432         161 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1433             :     }
    1434             :   }
    1435     4991083 :   q = p >> 1;
    1436     4991083 :   split_init(&S, lg(f)-1);
    1437     4991058 :   settyp(S.done, t_VECSMALL);
    1438     4991058 :   if (v) split_add_done(&S, (GEN)0);
    1439     4991058 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1440          42 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1441     4990703 :   pol = polx_Flx(f[1]);
    1442    10305261 :   for (pol[2]=1; ; pol[2]++)
    1443             :   {
    1444    10305261 :     long j, l = lg(S.todo);
    1445    10305261 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1446     5314588 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1447    10975360 :     for (j = 1; j < l; j++)
    1448             :     {
    1449     5661331 :       GEN c = gel(S.todo,j);
    1450     5661331 :       switch(degpol(c))
    1451             :       {
    1452             :         case 1:
    1453     4868542 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1454     4868532 :           j--; l--; break;
    1455             :         case 2: {
    1456      362800 :           ulong r = Flx_quad_root(c, p, 0), s;
    1457      362798 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1458      362791 :           s = Flx_otherroot(c,r, p);
    1459      362789 :           split_done(&S, j, (GEN)r, (GEN)s);
    1460      362795 :           j--; l--; break;
    1461             :         }
    1462             :         default: {
    1463      429826 :           GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
    1464             :           long db;
    1465      429822 :           b = Flx_gcd(c,b, p); db = degpol(b);
    1466      429829 :           if (db && db < degpol(c))
    1467             :           {
    1468      331427 :             b = Flx_normalize(b, p);
    1469      331425 :             c = Flx_div(c,b, p);
    1470      331424 :             split_todo(&S, j, b, c);
    1471             :           }
    1472             :         }
    1473             :       }
    1474             :     }
    1475     5314029 :   }
    1476             : }
    1477             : 
    1478             : GEN
    1479    21936451 : Flx_roots(GEN f, ulong p)
    1480             : {
    1481    21936451 :   pari_sp av = avma;
    1482    21936451 :   switch(lg(f))
    1483             :   {
    1484           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1485           0 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1486             :   }
    1487    21947347 :   if (p == 2) return Flx_root_mod_2(f);
    1488    21947347 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1489             : }
    1490             : 
    1491             : /* assume x reduced mod p, monic. */
    1492             : static int
    1493      720734 : Flx_quad_factortype(GEN x, ulong p)
    1494             : {
    1495      720734 :   ulong b = x[3], c = x[2];
    1496      720734 :   return krouu(Fl_disc_bc(b, c, p), p);
    1497             : }
    1498             : static GEN
    1499           0 : Flx_is_irred_2(GEN f, ulong p, long d)
    1500             : {
    1501           0 :   if (!d) return NULL;
    1502           0 :   if (d == 1) return gen_1;
    1503           0 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1504             : }
    1505             : static GEN
    1506      737877 : Flx_degfact_2(GEN f, ulong p, long d)
    1507             : {
    1508      737877 :   if (!d) return trivial_fact();
    1509      737877 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1510      720734 :   switch(Flx_quad_factortype(f, p)) {
    1511      342412 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1512      369894 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1513        8428 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1514             :   }
    1515             : }
    1516             : /* p > 2 */
    1517             : static GEN
    1518      409431 : Flx_factor_2(GEN f, ulong p, long d)
    1519             : {
    1520             :   ulong r, s;
    1521             :   GEN R,S;
    1522      409431 :   long v = f[1];
    1523      409431 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1524      409431 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1525      395483 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1526      297444 :   r = Flx_quad_root(f, p, 1);
    1527      297444 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1528      195135 :   s = Flx_otherroot(f, r, p);
    1529      195135 :   r = Fl_neg(r, p);
    1530      195135 :   s = Fl_neg(s, p);
    1531      195135 :   if (s < r) lswap(s,r);
    1532      195135 :   R = mkvecsmall3(v,r,1);
    1533      195135 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1534      168626 :   S = mkvecsmall3(v,s,1);
    1535      168626 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1536             : }
    1537             : static GEN
    1538     1147308 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1539             : {
    1540     1147308 :   switch(flag) {
    1541           0 :     case 2: return Flx_is_irred_2(f, p, d);
    1542      737877 :     case 1: return Flx_degfact_2(f, p, d);
    1543      409431 :     default: return Flx_factor_2(f, p, d);
    1544             :   }
    1545             : }
    1546             : 
    1547             : void
    1548       18452 : F2xV_to_FlxV_inplace(GEN v)
    1549             : {
    1550             :   long i;
    1551       18452 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1552       18452 : }
    1553             : void
    1554      712526 : FlxV_to_ZXV_inplace(GEN v)
    1555             : {
    1556             :   long i;
    1557      712526 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1558      712526 : }
    1559             : void
    1560      141379 : F2xV_to_ZXV_inplace(GEN v)
    1561             : {
    1562             :   long i;
    1563      141379 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1564      141379 : }
    1565             : 
    1566             : /* Adapted from Shoup NTL */
    1567             : GEN
    1568       76787 : F2x_factor_squarefree(GEN f)
    1569             : {
    1570             :   GEN r, t, v, tv;
    1571       76787 :   long i, q, n = F2x_degree(f);
    1572       76787 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1573      121975 :   for(q = 1;;q *= 2)
    1574             :   {
    1575      121975 :     r = F2x_gcd(f, F2x_deriv(f));
    1576      121973 :     if (F2x_degree(r) == 0)
    1577             :     {
    1578       64185 :       gel(u, q) = f;
    1579       64185 :       break;
    1580             :     }
    1581       57790 :     t = F2x_div(f, r);
    1582       57790 :     if (F2x_degree(t) > 0)
    1583             :     {
    1584             :       long j;
    1585       53699 :       for(j = 1;;j++)
    1586             :       {
    1587       53699 :         v = F2x_gcd(r, t);
    1588       53699 :         tv = F2x_div(t, v);
    1589       53698 :         if (F2x_degree(tv) > 0)
    1590       20374 :           gel(u, j*q) = tv;
    1591       53699 :         if (F2x_degree(v) <= 0) break;
    1592       34722 :         r = F2x_div(r, v);
    1593       34722 :         t = v;
    1594       34722 :       }
    1595       18976 :       if (F2x_degree(r) == 0) break;
    1596             :     }
    1597       45189 :     f = F2x_sqrt(r);
    1598       45188 :   }
    1599      523445 :   for (i = n; i; i--)
    1600      518017 :     if (F2x_degree(gel(u,i))) break;
    1601       76804 :   setlg(u,i+1); return u;
    1602             : }
    1603             : 
    1604             : static GEN
    1605       26396 : F2x_ddf_simple(GEN T, GEN XP)
    1606             : {
    1607       26396 :   pari_sp av = avma, av2;
    1608             :   GEN f, z, Tr, X;
    1609       26396 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1610       26396 :   if (n == 0) return cgetg(1, t_VEC);
    1611       26396 :   if (n == 1) return mkvec(T);
    1612       19828 :   z = XP; Tr = T; X = polx_F2x(v);
    1613       19828 :   f = const_vec(n, pol1_F2x(v));
    1614       19828 :   av2 = avma;
    1615      159553 :   for (j = 1; j <= B; j++)
    1616             :   {
    1617      149635 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1618      149639 :     if (F2x_degree(u))
    1619             :     {
    1620       42255 :       gel(f, j) = u;
    1621       42255 :       Tr = F2x_div(Tr, u);
    1622       42250 :       av2 = avma;
    1623      107419 :     } else z = gerepileuptoleaf(av2, z);
    1624      149685 :     if (!F2x_degree(Tr)) break;
    1625      139742 :     z = F2xq_sqr(z, Tr);
    1626             :   }
    1627       19827 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1628       19827 :   return gerepilecopy(av, f);
    1629             : }
    1630             : 
    1631             : static GEN
    1632        8329 : F2xq_frobtrace(GEN a, long d, GEN T)
    1633             : {
    1634        8329 :   pari_sp av = avma;
    1635             :   long i;
    1636        8329 :   GEN x = a;
    1637       54422 :   for(i=1; i<d; i++)
    1638             :   {
    1639       46094 :     x = F2x_add(a, F2xq_sqr(x,T));
    1640       46092 :     if (gc_needed(av, 2))
    1641           0 :       x = gerepileuptoleaf(av, x);
    1642             :   }
    1643        8328 :   return x;
    1644             : }
    1645             : 
    1646             : static void
    1647       12795 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1648             : {
    1649       12795 :   long n = F2x_degree(Tp), r = n/d;
    1650             :   GEN T, f, ff;
    1651       25590 :   if (r==1) { gel(V, idx) = Tp; return; }
    1652        4386 :   T = Tp;
    1653        4386 :   XP = F2x_rem(XP, T);
    1654             :   while (1)
    1655             :   {
    1656        8328 :     pari_sp btop = avma;
    1657             :     long df;
    1658        8328 :     GEN g = random_F2x(n, Tp[1]);
    1659        8328 :     GEN t = F2xq_frobtrace(g, d, T);
    1660        8328 :     if (lgpol(t) == 0) continue;
    1661        6403 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1662        6403 :     if (df > 0 && df < n) break;
    1663        2017 :     avma = btop;
    1664        3942 :   }
    1665        4386 :   ff = F2x_div(Tp, f);
    1666        4386 :   F2x_edf_simple(f, XP, d, V, idx);
    1667        4386 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1668             : }
    1669             : 
    1670             : static GEN
    1671       24451 : F2x_factor_Shoup(GEN T)
    1672             : {
    1673       24451 :   long i, n, s = 0;
    1674             :   GEN XP, D, V;
    1675             :   pari_timer ti;
    1676       24451 :   n = F2x_degree(T);
    1677       24451 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1678       24451 :   XP = F2x_Frobenius(T);
    1679       24450 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1680       24450 :   D = F2x_ddf_simple(T, XP);
    1681       24450 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1682      363130 :   for (i = 1; i <= n; i++)
    1683      338680 :     s += F2x_degree(gel(D,i))/i;
    1684       24450 :   V = cgetg(s+1, t_COL);
    1685      363097 :   for (i = 1, s = 1; i <= n; i++)
    1686             :   {
    1687      338648 :     GEN Di = gel(D,i);
    1688      338648 :     long ni = F2x_degree(Di), ri = ni/i;
    1689      338632 :     if (ni == 0) continue;
    1690       56085 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1691        4023 :     F2x_edf_simple(Di, XP, i, V, s);
    1692        4023 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1693        4023 :     s += ri;
    1694             :   }
    1695       24449 :   return V;
    1696             : }
    1697             : 
    1698             : static GEN
    1699        1778 : F2x_simplefact_Shoup(GEN T)
    1700             : {
    1701        1778 :   long i, n, s = 0, j = 1, k;
    1702             :   GEN XP, D, V;
    1703             :   pari_timer ti;
    1704        1778 :   n = F2x_degree(T);
    1705        1778 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1706        1778 :   XP = F2x_Frobenius(T);
    1707        1778 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1708        1778 :   D = F2x_ddf_simple(T, XP);
    1709        1778 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1710        9065 :   for (i = 1; i <= n; i++)
    1711        7287 :     s += F2x_degree(gel(D,i))/i;
    1712        1778 :   V = cgetg(s+1, t_VECSMALL);
    1713        9065 :   for (i = 1; i <= n; i++)
    1714             :   {
    1715        7287 :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1716        7287 :     if (ni == 0) continue;
    1717        4333 :     for (k = 1; k <= ri; k++)
    1718        2219 :       V[j++] = i;
    1719             :   }
    1720        1778 :   return V;
    1721             : }
    1722             : 
    1723             : static GEN
    1724       17791 : F2x_factor_Cantor(GEN T)
    1725             : {
    1726       17791 :   GEN E, F, V = F2x_factor_squarefree(T);
    1727       17791 :   long i, j, l = lg(V);
    1728       17791 :   E = cgetg(l, t_VEC);
    1729       17791 :   F = cgetg(l, t_VEC);
    1730       46733 :   for (i=1, j=1; i < l; i++)
    1731       28942 :     if (F2x_degree(gel(V,i)))
    1732             :     {
    1733       24451 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1734       24449 :       gel(F, j) = Fj;
    1735       24449 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1736       24451 :       j++;
    1737             :     }
    1738       17791 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1739             : }
    1740             : 
    1741             : static GEN
    1742        1736 : F2x_simplefact_Cantor(GEN T)
    1743             : {
    1744        1736 :   GEN E, F, V = F2x_factor_squarefree(T);
    1745        1736 :   long i, j, l = lg(V);
    1746        1736 :   F = cgetg(l, t_VEC);
    1747        1736 :   E = cgetg(l, t_VEC);
    1748        4690 :   for (i=1, j=1; i < l; i++)
    1749        2954 :     if (F2x_degree(gel(V,i)))
    1750             :     {
    1751        1778 :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1752        1778 :       gel(F, j) = Fj;
    1753        1778 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1754        1778 :       j++;
    1755             :     }
    1756        1736 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1757             : }
    1758             : 
    1759             : static int
    1760         238 : F2x_isirred_Cantor(GEN T)
    1761             : {
    1762         238 :   pari_sp av = avma;
    1763             :   pari_timer ti;
    1764             :   long n, d;
    1765         238 :   GEN dT = F2x_deriv(T);
    1766             :   GEN XP, D;
    1767         238 :   if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
    1768         168 :   n = F2x_degree(T);
    1769         168 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1770         168 :   XP = F2x_Frobenius(T);
    1771         168 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1772         168 :   D = F2x_ddf_simple(T, XP);
    1773         168 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1774         168 :   d = F2x_degree(gel(D, n));
    1775         168 :   avma = av; return d==n;
    1776             : }
    1777             : 
    1778             : static GEN
    1779       30671 : F2x_factcantor_i(GEN f, long flag)
    1780             : {
    1781       30671 :   long d = F2x_degree(f);
    1782       30671 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1783       19765 :   switch(flag)
    1784             :   {
    1785       17791 :     default: return F2x_factor_Cantor(f);
    1786        1736 :     case 1: return F2x_simplefact_Cantor(f);
    1787         238 :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1788             :   }
    1789             : }
    1790             : 
    1791             : GEN
    1792       15152 : F2x_factcantor(GEN f, long flag)
    1793             : {
    1794       15152 :   pari_sp av = avma;
    1795       15152 :   GEN z = F2x_factcantor_i(f, flag);
    1796       15152 :   if (flag == 2) { avma = av; return z; }
    1797       15152 :   return gerepilecopy(av, z);
    1798             : }
    1799             : 
    1800             : GEN
    1801           0 : F2x_degfact(GEN f)
    1802             : {
    1803           0 :   pari_sp av = avma;
    1804           0 :   GEN z = F2x_factcantor_i(f, 1);
    1805           0 :   return gerepilecopy(av, z);
    1806             : }
    1807             : 
    1808             : int
    1809         238 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f, 2); }
    1810             : 
    1811             : /* Adapted from Shoup NTL */
    1812             : GEN
    1813      742871 : Flx_factor_squarefree(GEN f, ulong p)
    1814             : {
    1815      742871 :   long i, q, n = degpol(f);
    1816      742871 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1817      796175 :   for(q = 1;;q *= p)
    1818             :   {
    1819      796175 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1820      796175 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1821       89182 :     t = Flx_div(f, r, p);
    1822       89182 :     if (degpol(t) > 0)
    1823             :     {
    1824             :       long j;
    1825      115504 :       for(j = 1;;j++)
    1826             :       {
    1827      115504 :         v = Flx_gcd(r, t, p);
    1828      115504 :         tv = Flx_div(t, v, p);
    1829      115504 :         if (degpol(tv) > 0)
    1830       52515 :           gel(u, j*q) = Flx_normalize(tv, p);
    1831      115504 :         if (degpol(v) <= 0) break;
    1832       79158 :         r = Flx_div(r, v, p);
    1833       79158 :         t = v;
    1834       79158 :       }
    1835       36346 :       if (degpol(r) == 0) break;
    1836             :     }
    1837       53304 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1838       53304 :   }
    1839     3822487 :   for (i = n; i; i--)
    1840     3819806 :     if (degpol(gel(u,i))) break;
    1841      742871 :   setlg(u,i+1); return u;
    1842             : }
    1843             : 
    1844             : long
    1845         126 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1846             : {
    1847         126 :   pari_sp av = avma;
    1848             :   ulong lc;
    1849             :   GEN F;
    1850         126 :   long i, n = degpol(f), v = f[1], l;
    1851         126 :   if (n % k) return 0;
    1852         126 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    1853         126 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    1854         126 :   F = Flx_factor_squarefree(f, p); l = lg(F)-1;
    1855        6405 :   for (i = 1; i <= l; i++)
    1856        6300 :     if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
    1857         105 :   if (pt_r)
    1858             :   {
    1859         105 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    1860        6384 :     for(i = l; i >= 1; i--)
    1861             :     {
    1862        6279 :       if (i%k) continue;
    1863        1274 :       s = Flx_mul(s, gel(F,i), p);
    1864        1274 :       r = Flx_mul(r, s, p);
    1865             :     }
    1866         105 :     *pt_r = gerepileuptoleaf(av, r);
    1867           0 :   } else av = avma;
    1868         105 :   return 1;
    1869             : }
    1870             : 
    1871             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1872             : static GEN
    1873      653450 : Flx_ddf(GEN T, GEN XP, ulong p)
    1874             : {
    1875      653450 :   pari_sp av = avma;
    1876             :   GEN b, g, h, F, f, Tr, xq;
    1877             :   long i, j, n, v, bo, ro;
    1878             :   long B, l, m;
    1879             :   pari_timer ti;
    1880      653450 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1881      653450 :   if (n == 0) return cgetg(1, t_VEC);
    1882      653450 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1883      642116 :   B = n/2;
    1884      642116 :   l = usqrt(B);
    1885      642116 :   m = (B+l-1)/l;
    1886      642116 :   T = Flx_get_red(T, p);
    1887      642116 :   b = cgetg(l+2, t_VEC);
    1888      642116 :   gel(b, 1) = polx_Flx(v);
    1889      642116 :   gel(b, 2) = XP;
    1890      642116 :   bo = brent_kung_optpow(n, l-1, 1);
    1891      642116 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    1892      642116 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1893      642116 :   if (expu(p) <= ro)
    1894      151347 :     for (i = 3; i <= l+1; i++)
    1895       85470 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1896             :   else
    1897             :   {
    1898      576239 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1899      576239 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
    1900      675981 :     for (i = 3; i <= l+1; i++)
    1901       99742 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1902             :   }
    1903      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
    1904      642116 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1905      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
    1906      642116 :   g = cgetg(m+1, t_VEC);
    1907      642116 :   gel(g, 1) = gel(xq, 2);
    1908     1251721 :   for(i = 2; i <= m; i++)
    1909      609605 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1910      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
    1911      642116 :   h = cgetg(m+1, t_VEC);
    1912     1893836 :   for (j = 1; j <= m; j++)
    1913             :   {
    1914     1251720 :     pari_sp av = avma;
    1915     1251720 :     GEN gj = gel(g, j);
    1916     1251720 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1917     1923868 :     for (i = 2; i <= l; i++)
    1918      672150 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1919     1251718 :     gel(h, j) = gerepileupto(av, e);
    1920             :   }
    1921      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
    1922      642116 :   Tr = get_Flx_mod(T);
    1923      642116 :   F = cgetg(m+1, t_VEC);
    1924     1893837 :   for (j = 1; j <= m; j++)
    1925             :   {
    1926     1251721 :     gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
    1927     1251721 :     Tr = Flx_div(Tr, gel(F,j), p);
    1928             :   }
    1929      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
    1930      642116 :   f = const_vec(n, pol1_Flx(v));
    1931     1893836 :   for (j = 1; j <= m; j++)
    1932             :   {
    1933     1251720 :     GEN e = gel(F, j);
    1934     1390015 :     for (i=l-1; i >= 0; i--)
    1935             :     {
    1936     1390016 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    1937     1390016 :       if (degpol(u))
    1938             :       {
    1939      525134 :         gel(f, l*j-i) = u;
    1940      525134 :         e = Flx_div(e, u, p);
    1941             :       }
    1942     1390014 :       if (!degpol(e)) break;
    1943             :     }
    1944             :   }
    1945      642116 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
    1946      642116 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    1947      642116 :   return gerepilecopy(av, f);
    1948             : }
    1949             : 
    1950             : static void
    1951       45743 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1952             : {
    1953       45743 :   long n = degpol(Tp), r = n/d;
    1954             :   GEN T, f, ff;
    1955             :   ulong p2;
    1956       91486 :   if (r==1) { gel(V, idx) = Tp; return; }
    1957       20646 :   p2 = p>>1;
    1958       20646 :   T = Flx_get_red(Tp, p);
    1959       20646 :   XP = Flx_rem(XP, T, p);
    1960             :   while (1)
    1961             :   {
    1962       22195 :     pari_sp btop = avma;
    1963             :     long i;
    1964       22195 :     GEN g = random_Flx(n, Tp[1], p);
    1965       22195 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1966       22195 :     if (lgpol(t) == 0) continue;
    1967       43697 :     for(i=1; i<=10; i++)
    1968             :     {
    1969       42511 :       pari_sp btop2 = avma;
    1970       42511 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    1971       42510 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    1972       42511 :       if (degpol(f) > 0 && degpol(f) < n) break;
    1973       21865 :       avma = btop2;
    1974             :     }
    1975       21832 :     if (degpol(f) > 0 && degpol(f) < n) break;
    1976        1186 :     avma = btop;
    1977        1549 :   }
    1978       20646 :   f = Flx_normalize(f, p);
    1979       20646 :   ff = Flx_div(Tp, f ,p);
    1980       20646 :   Flx_edf_simple(f, XP, d, p, V, idx);
    1981       20646 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    1982             : }
    1983             : static void
    1984             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    1985             : 
    1986             : static void
    1987       32701 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    1988             : {
    1989             :   pari_sp av;
    1990       32701 :   GEN Tp = get_Flx_mod(T);
    1991       32701 :   long n = degpol(hp), vT = Tp[1];
    1992             :   GEN u1, u2, f1, f2;
    1993       32701 :   ulong p2 = p>>1;
    1994             :   GEN R, h;
    1995       32701 :   h = Flx_get_red(hp, p);
    1996       32701 :   t = Flx_rem(t, T, p);
    1997       32701 :   av = avma;
    1998             :   do
    1999             :   {
    2000       55865 :     avma = av;
    2001       55865 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    2002       55865 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    2003       55865 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2004       32701 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    2005       32701 :   f1 = Flx_normalize(f1, p);
    2006       32701 :   u2 = Flx_div(hp, u1, p);
    2007       32701 :   f2 = Flx_div(Tp, f1, p);
    2008       32701 :   if (degpol(u1)==1)
    2009             :   {
    2010       25850 :     if (degpol(f1)==d)
    2011       25489 :       gel(V, idx) = f1;
    2012             :     else
    2013         361 :       Flx_edf(f1, XP, d, p, V, idx);
    2014             :   }
    2015             :   else
    2016        6851 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    2017       32701 :   idx += degpol(f1)/d;
    2018       32701 :   if (degpol(u2)==1)
    2019             :   {
    2020       25728 :     if (degpol(f2)==d)
    2021       25278 :       gel(V, idx) = f2;
    2022             :     else
    2023         450 :       Flx_edf(f2, XP, d, p, V, idx);
    2024             :   }
    2025             :   else
    2026        6973 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    2027       32701 : }
    2028             : 
    2029             : static void
    2030       18877 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2031             : {
    2032       18877 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2033             :   GEN T, h, t;
    2034             :   pari_timer ti;
    2035       37754 :   if (r==1) { gel(V, idx) = Tp; return; }
    2036       18877 :   T = Flx_get_red(Tp, p);
    2037       18877 :   XP = Flx_rem(XP, T, p);
    2038       18877 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2039             :   do
    2040             :   {
    2041       20817 :     GEN g = random_Flx(n, vT, p);
    2042       20817 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2043       20817 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2044       20817 :     h = Flxq_minpoly(t, T, p);
    2045       20817 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2046       20817 :   } while (degpol(h) <= 1);
    2047       18877 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    2048             : }
    2049             : 
    2050             : static GEN
    2051       33012 : Flx_factor_Shoup(GEN T, ulong p)
    2052             : {
    2053       33012 :   long i, n, s = 0;
    2054             :   GEN XP, D, V;
    2055       33012 :   long e = expu(p);
    2056             :   pari_timer ti;
    2057       33012 :   n = get_Flx_degree(T);
    2058       33012 :   T = Flx_get_red(T, p);
    2059       33012 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2060       33012 :   XP = Flx_Frobenius(T, p);
    2061       33012 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2062       33012 :   D = Flx_ddf(T, XP, p);
    2063       33012 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2064      339031 :   for (i = 1; i <= n; i++)
    2065      306019 :     s += degpol(gel(D,i))/i;
    2066       33012 :   V = cgetg(s+1, t_COL);
    2067      339031 :   for (i = 1, s = 1; i <= n; i++)
    2068             :   {
    2069      306019 :     GEN Di = gel(D,i);
    2070      306019 :     long ni = degpol(Di), ri = ni/i;
    2071      306019 :     if (ni == 0) continue;
    2072       49685 :     Di = Flx_normalize(Di, p);
    2073       49685 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2074       22517 :     if (ri <= e*expu(e))
    2075       18066 :       Flx_edf(Di, XP, i, p, V, s);
    2076             :     else
    2077        4451 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2078       22517 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2079       22517 :     s += ri;
    2080             :   }
    2081       33012 :   return V;
    2082             : }
    2083             : 
    2084             : static GEN
    2085      351990 : Flx_simplefact_Shoup(GEN T, ulong p)
    2086             : {
    2087             :   GEN XP, D;
    2088             :   pari_timer ti;
    2089      351990 :   long n = get_Flx_degree(T);
    2090      351990 :   T = Flx_get_red(T, p);
    2091      351990 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2092      351990 :   XP = Flx_Frobenius(T, p);
    2093      351990 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2094      351990 :   D = Flx_ddf(T, XP, p);
    2095      351990 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2096      351990 :   return ddf_to_simplefact(D, n);
    2097             : }
    2098             : 
    2099             : static GEN
    2100       29894 : Flx_factor_Cantor(GEN T, ulong p)
    2101             : {
    2102       29894 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2103       29894 :   long i, j, l = lg(V);
    2104       29894 :   F = cgetg(l, t_VEC);
    2105       29894 :   E = cgetg(l, t_VEC);
    2106       65946 :   for (i=1, j=1; i < l; i++)
    2107       36052 :     if (degpol(gel(V,i)))
    2108             :     {
    2109       33012 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2110       33012 :       gel(F, j) = Fj;
    2111       33012 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2112       33012 :       j++;
    2113             :     }
    2114       29894 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2115             : }
    2116             : 
    2117             : static GEN
    2118      351008 : Flx_simplefact_Cantor(GEN T, ulong p)
    2119             : {
    2120      351008 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2121      351008 :   long i, j, l = lg(V);
    2122      351008 :   F = cgetg(l, t_VEC);
    2123      351008 :   E = cgetg(l, t_VEC);
    2124      705001 :   for (i=1, j=1; i < l; i++)
    2125      353993 :     if (degpol(gel(V,i)))
    2126             :     {
    2127      351990 :       GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
    2128      351990 :       gel(F, j) = Fj;
    2129      351990 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2130      351990 :       j++;
    2131             :     }
    2132      351008 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    2133             : }
    2134             : 
    2135             : static int
    2136         581 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2137             : {
    2138         581 :   pari_sp av = avma;
    2139             :   pari_timer ti;
    2140             :   long n, d;
    2141         581 :   GEN T = get_Flx_mod(Tp);
    2142         581 :   GEN dT = Flx_deriv(T, p);
    2143             :   GEN XP, D;
    2144         581 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    2145         441 :   n = get_Flx_degree(T);
    2146         441 :   T = Flx_get_red(Tp, p);
    2147         441 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2148         441 :   XP = Flx_Frobenius(T, p);
    2149         441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2150         441 :   D = Flx_ddf(T, XP, p);
    2151         441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2152         441 :   d = degpol(gel(D, n));
    2153         441 :   avma = av; return d==n;
    2154             : }
    2155             : 
    2156             : static GEN
    2157     1138992 : Flx_factcantor_i(GEN f, ulong pp, long flag)
    2158             : {
    2159             :   long d;
    2160     1138992 :   if (pp==2) { /*We need to handle 2 specially */
    2161       15246 :     GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2162       15246 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2163       15246 :     return F;
    2164             :   }
    2165     1123746 :   d = degpol(f);
    2166     1123746 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2167      381483 :   switch(flag)
    2168             :   {
    2169       29894 :     default: return Flx_factor_Cantor(f, pp);
    2170      351008 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2171         581 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2172             :   }
    2173             : }
    2174             : 
    2175             : GEN
    2176        7622 : Flx_factcantor(GEN f, ulong p, long flag)
    2177             : {
    2178        7622 :   pari_sp av = avma;
    2179        7622 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
    2180        7622 :   if (flag == 2) { avma = av; return z; }
    2181        7622 :   return gerepilecopy(av, z);
    2182             : }
    2183             : 
    2184             : GEN
    2185     1097362 : Flx_degfact(GEN f, ulong p)
    2186             : {
    2187     1097362 :   pari_sp av = avma;
    2188     1097362 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
    2189     1097362 :   return gerepilecopy(av, z);
    2190             : }
    2191             : 
    2192             : /* T must be squarefree mod p*/
    2193             : GEN
    2194      186688 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2195             : {
    2196             :   GEN XP, D;
    2197             :   pari_timer ti;
    2198      186688 :   long i, s, n = get_Flx_degree(T);
    2199      186688 :   GEN V = const_vecsmall(n, 0);
    2200      186688 :   pari_sp av = avma;
    2201      186688 :   T = Flx_get_red(T, p);
    2202      186688 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2203      186688 :   XP = Flx_Frobenius(T, p);
    2204      186688 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2205      186688 :   D = Flx_ddf(T, XP, p);
    2206      186688 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2207     1203945 :   for (i = 1, s = 0; i <= n; i++)
    2208             :   {
    2209     1017257 :     V[i] = degpol(gel(D,i))/i;
    2210     1017257 :     s += V[i];
    2211             :   }
    2212      186688 :   *nb = s;
    2213      186688 :   avma = av; return V;
    2214             : }
    2215             : 
    2216             : long
    2217       81319 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2218             : {
    2219       81319 :   pari_sp av = avma;
    2220       81319 :   long s = ddf_to_nbfact(Flx_ddf(T, XP, p));
    2221       81319 :   avma = av; return s;
    2222             : }
    2223             : 
    2224             : /* T must be squarefree mod p*/
    2225             : long
    2226       81319 : Flx_nbfact(GEN T, ulong p)
    2227             : {
    2228       81319 :   pari_sp av = avma;
    2229       81319 :   GEN XP = Flx_Frobenius(T, p);
    2230       81319 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2231       81319 :   avma = av; return n;
    2232             : }
    2233             : 
    2234             : int
    2235         581 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
    2236             : 
    2237             : /* factor f (FpX or Flx) mod pp.
    2238             :  * flag = 1: return the degrees, not the factors
    2239             :  * flag = 2: return NULL if f is not irreducible.
    2240             :  * Not gerepile-safe */
    2241             : static GEN
    2242         112 : factcantor_i(GEN f, GEN pp, long flag)
    2243             : {
    2244         112 :   if (typ(f) == t_VECSMALL)
    2245             :   { /* lgefint(pp) = 3 */
    2246             :     GEN F;
    2247          63 :     ulong p = pp[2];
    2248          63 :     if (p==2) {
    2249          35 :       F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2250          35 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2251             :     } else {
    2252          28 :       F = Flx_factcantor_i(f,p,flag);
    2253          28 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2254             :     }
    2255          63 :     return F;
    2256             :   }
    2257          49 :   return FpX_factcantor_i(f, pp, flag);
    2258             : }
    2259             : GEN
    2260           0 : FpX_factcantor(GEN f, GEN pp, long flag)
    2261             : {
    2262           0 :   pari_sp av = avma;
    2263             :   GEN z;
    2264           0 :   ZX_factmod_init(&f,pp);
    2265           0 :   z = factcantor_i(f,pp,flag);
    2266           0 :   if (flag == 2) { avma = av; return z; }
    2267           0 :   return gerepilecopy(av, z);
    2268             : }
    2269             : 
    2270             : static GEN
    2271         266 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
    2272             : {
    2273         266 :   pari_sp av = avma;
    2274             :   long j, lfact;
    2275             :   GEN z, t, E, y, u, v;
    2276             : 
    2277         266 :   factmod_init(&f, p);
    2278         266 :   switch(lg(f))
    2279             :   {
    2280          49 :     case 3: avma = av; return trivial_fact();
    2281          14 :     case 2: return gerepileupto(av, zero_fact_intmod(f, p));
    2282             :   }
    2283         203 :   z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
    2284         203 :   lfact = lg(t); y = cgetg(3, t_MAT);
    2285         203 :   gel(y,1) = u = cgetg(lfact,t_COL);
    2286         203 :   gel(y,2) = v = cgetg(lfact,t_COL);
    2287         203 :   if (flag)
    2288         119 :     for (j=1; j<lfact; j++)
    2289             :     {
    2290          77 :       gel(u,j) = utoi(uel(t,j));
    2291          77 :       gel(v,j) = utoi(uel(E,j));
    2292             :     }
    2293             :   else
    2294        8078 :     for (j=1; j<lfact; j++)
    2295             :     {
    2296        7917 :       gel(u,j) = FpX_to_mod(gel(t,j), p);
    2297        7917 :       gel(v,j) = utoi(uel(E,j));
    2298             :     }
    2299         203 :   return gerepileupto(av, y);
    2300             : }
    2301             : GEN
    2302         112 : factcantor0(GEN f, GEN p, long flag)
    2303         112 : { return factmod_aux(f, p, &factcantor_i, flag); }
    2304             : GEN
    2305         154 : factmod(GEN f, GEN p)
    2306         154 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
    2307             : 
    2308             : /* Use this function when you think f is reducible, and that there are lots of
    2309             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2310             : int
    2311          14 : FpX_is_irred(GEN f, GEN p) {
    2312          14 :   ZX_factmod_init(&f,p);
    2313          14 :   return !!factcantor_i(f,p,2);
    2314             : }
    2315             : GEN
    2316           0 : FpX_degfact(GEN f, GEN p) {
    2317           0 :   pari_sp av = avma;
    2318             :   GEN z;
    2319           0 :   ZX_factmod_init(&f,p);
    2320           0 :   z = factcantor_i(f,p,1);
    2321           0 :   return gerepilecopy(av, z);
    2322             : }
    2323             : GEN
    2324          70 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
    2325             : GEN
    2326          42 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
    2327             : 
    2328             : /* set x <-- x + c*y mod p */
    2329             : /* x is not required to be normalized.*/
    2330             : static void
    2331     1246893 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2332             : {
    2333             :   long i, lx, ly;
    2334     1246893 :   ulong *x=(ulong *)gx;
    2335     1246893 :   ulong *y=(ulong *)gy;
    2336     2493786 :   if (!c) return;
    2337     1213188 :   lx = lg(gx);
    2338     1213188 :   ly = lg(gy);
    2339     1213188 :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2340     1213188 :   if (SMALL_ULONG(p))
    2341     1210710 :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2342             :   else
    2343        2478 :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2344             : }
    2345             : 
    2346             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    2347             : /* assume x1 != 0 */
    2348             : static GEN
    2349      151412 : deg1_Flx(ulong x1, ulong x0, ulong sv)
    2350             : {
    2351      151412 :   return mkvecsmall3(sv, x0, x1);
    2352             : }
    2353             : 
    2354             : static long
    2355       52920 : F2x_split_Berlekamp(GEN *t)
    2356             : {
    2357       52920 :   GEN u = *t, a, b, vker;
    2358       52920 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    2359             : 
    2360       52920 :   if (du == 1) return 1;
    2361       26049 :   if (du == 2)
    2362             :   {
    2363        6594 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    2364             :     {
    2365           0 :       t[0] = mkvecsmall2(sv, 2);
    2366           0 :       t[1] = mkvecsmall2(sv, 3);
    2367           0 :       return 2;
    2368             :     }
    2369        6594 :     return 1;
    2370             :   }
    2371             : 
    2372       19455 :   vker = F2x_Berlekamp_ker(u);
    2373       19455 :   lb = lgcols(vker);
    2374       19455 :   d = lg(vker)-1;
    2375       19455 :   ir = 0;
    2376             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2377       45980 :   for (L=1; L<d; )
    2378             :   {
    2379             :     GEN pol;
    2380        7070 :     if (d == 2)
    2381        6139 :       pol = F2v_to_F2x(gel(vker,2), sv);
    2382             :     else
    2383             :     {
    2384         931 :       GEN v = zero_zv(lb);
    2385         931 :       v[1] = du;
    2386         931 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    2387        3045 :       for (i=2; i<=d; i++)
    2388        2114 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    2389         931 :       pol = F2v_to_F2x(v, sv);
    2390             :     }
    2391       14959 :     for (i=ir; i<L && L<d; i++)
    2392             :     {
    2393        7889 :       a = t[i]; la = F2x_degree(a);
    2394        7889 :       if (la == 1) { set_irred(i); }
    2395        7868 :       else if (la == 2)
    2396             :       {
    2397          14 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    2398             :         {
    2399           0 :           t[i] = mkvecsmall2(sv, 2);
    2400           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    2401             :         }
    2402          14 :         set_irred(i);
    2403             :       }
    2404             :       else
    2405             :       {
    2406        7854 :         pari_sp av = avma;
    2407             :         long lb;
    2408        7854 :         b = F2x_rem(pol, a);
    2409        7854 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    2410        6748 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    2411        6748 :         if (lb && lb < la)
    2412             :         {
    2413        6748 :           t[L] = F2x_div(a,b);
    2414        6748 :           t[i]= b; L++;
    2415             :         }
    2416           0 :         else avma = av;
    2417             :       }
    2418             :     }
    2419             :   }
    2420       19455 :   return d;
    2421             : }
    2422             : 
    2423             : /* p != 2 */
    2424             : static long
    2425      370782 : Flx_split_Berlekamp(GEN *t, ulong p)
    2426             : {
    2427      370782 :   GEN u = *t, a,b,vker;
    2428      370782 :   long d, i, ir, L, la, lb, sv = u[1];
    2429      370782 :   long l = lg(u);
    2430             :   ulong po2;
    2431             : 
    2432      370782 :   if (p == 2)
    2433             :   {
    2434           0 :     *t = Flx_to_F2x(*t);
    2435           0 :     d = F2x_split_Berlekamp(t);
    2436           0 :     for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
    2437           0 :     return d;
    2438             :   }
    2439      370782 :   la = degpol(u);
    2440      370782 :   if (la == 1) return 1;
    2441      296106 :   if (la == 2)
    2442             :   {
    2443       14956 :     ulong r = Flx_quad_root(u,p,1);
    2444       14956 :     if (r != p)
    2445             :     {
    2446        7102 :       t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
    2447        7102 :       t[1] = deg1_Flx(1, p - r, sv);
    2448        7102 :       return 2;
    2449             :     }
    2450        7854 :     return 1;
    2451             :   }
    2452             : 
    2453      281150 :   vker = Flx_Berlekamp_ker(u,p);
    2454      281150 :   vker = Flm_to_FlxV(vker, sv);
    2455      281150 :   d = lg(vker)-1;
    2456      281150 :   po2 = p >> 1; /* (p-1) / 2 */
    2457      281150 :   ir = 0;
    2458             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2459     1058345 :   for (L=1; L<d; )
    2460             :   {
    2461      496045 :     GEN pol = zero_zv(l-2);
    2462      496045 :     pol[1] = sv;
    2463      496045 :     pol[2] = random_Fl(p); /*Assume vker[1]=1*/
    2464     1742938 :     for (i=2; i<=d; i++)
    2465     1246893 :       Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
    2466      496045 :     (void)Flx_renormalize(pol,l-1);
    2467             : 
    2468     1292416 :     for (i=ir; i<L && L<d; i++)
    2469             :     {
    2470      796371 :       a = t[i]; la = degpol(a);
    2471      796371 :       if (la == 1) { set_irred(i); }
    2472      738828 :       else if (la == 2)
    2473             :       {
    2474      129324 :         ulong r = Flx_quad_root(a,p,1);
    2475      129324 :         if (r != p)
    2476             :         {
    2477       68604 :           t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
    2478       68604 :           t[L] = deg1_Flx(1, p - r, sv); L++;
    2479             :         }
    2480      129324 :         set_irred(i);
    2481             :       }
    2482             :       else
    2483             :       {
    2484      609504 :         pari_sp av = avma;
    2485      609504 :         b = Flx_rem(pol, a, p);
    2486      609504 :         if (degpol(b) <= 0) { avma=av; continue; }
    2487      559871 :         b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
    2488      559871 :         b = Flx_gcd(a,b, p); lb = degpol(b);
    2489      559871 :         if (lb && lb < la)
    2490             :         {
    2491      311835 :           b = Flx_normalize(b, p);
    2492      311835 :           t[L] = Flx_div(a,b,p);
    2493      311835 :           t[i]= b; L++;
    2494             :         }
    2495      248036 :         else avma = av;
    2496             :       }
    2497             :     }
    2498             :   }
    2499      281150 :   return d;
    2500             : }
    2501             : 
    2502             : static long
    2503         177 : FpX_split_Berlekamp(GEN *t, GEN p)
    2504             : {
    2505         177 :   GEN u = *t, a,b,po2,vker;
    2506         177 :   long d, i, ir, L, la, lb, vu = varn(u);
    2507         177 :   if (lgefint(p) == 3)
    2508             :   {
    2509           0 :     ulong up = p[2];
    2510           0 :     if (up == 2)
    2511             :     {
    2512           0 :       *t = ZX_to_F2x(*t);
    2513           0 :       d = F2x_split_Berlekamp(t);
    2514           0 :       for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
    2515             :     }
    2516             :     else
    2517             :     {
    2518           0 :       *t = ZX_to_Flx(*t, up);
    2519           0 :       d = Flx_split_Berlekamp(t, up);
    2520           0 :       for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
    2521             :     }
    2522           0 :     return d;
    2523             :   }
    2524         177 :   la = degpol(u);
    2525         177 :   if (la == 1) return 1;
    2526         173 :   if (la == 2)
    2527             :   {
    2528           2 :     GEN r = FpX_quad_root(u,p,1);
    2529           2 :     if (r)
    2530             :     {
    2531           0 :       t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
    2532           0 :       t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
    2533           0 :       return 2;
    2534             :     }
    2535           2 :     return 1;
    2536             :   }
    2537         171 :   vker = FpX_Berlekamp_ker(u,p);
    2538         171 :   vker = RgM_to_RgXV(vker,vu);
    2539         171 :   d = lg(vker)-1;
    2540         171 :   po2 = shifti(p, -1); /* (p-1) / 2 */
    2541         171 :   ir = 0;
    2542             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2543         384 :   for (L=1; L<d; )
    2544             :   {
    2545          42 :     GEN pol = scalar_ZX_shallow(randomi(p), vu);
    2546         164 :     for (i=2; i<=d; i++)
    2547         122 :       pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
    2548          42 :     pol = FpX_red(pol,p);
    2549         116 :     for (i=ir; i<L && L<d; i++)
    2550             :     {
    2551          74 :       a = t[i]; la = degpol(a);
    2552          74 :       if (la == 1) { set_irred(i); }
    2553          65 :       else if (la == 2)
    2554             :       {
    2555           7 :         GEN r = FpX_quad_root(a,p,1);
    2556           7 :         if (r)
    2557             :         {
    2558           7 :           t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
    2559           7 :           t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
    2560             :         }
    2561           7 :         set_irred(i);
    2562             :       }
    2563             :       else
    2564             :       {
    2565          58 :         pari_sp av = avma;
    2566          58 :         b = FpX_rem(pol, a, p);
    2567          58 :         if (degpol(b) <= 0) { avma=av; continue; }
    2568          58 :         b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
    2569          58 :         b = FpX_gcd(a,b, p); lb = degpol(b);
    2570          58 :         if (lb && lb < la)
    2571             :         {
    2572          46 :           b = FpX_normalize(b, p);
    2573          46 :           t[L] = FpX_div(a,b,p);
    2574          46 :           t[i]= b; L++;
    2575             :         }
    2576          12 :         else avma = av;
    2577             :       }
    2578             :     }
    2579             :   }
    2580         171 :   return d;
    2581             : }
    2582             : 
    2583             : static GEN
    2584      153342 : F2x_Berlekamp_i(GEN f, long flag)
    2585             : {
    2586      153342 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    2587             :   GEN y, E, t, V;
    2588             : 
    2589      153342 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    2590             : 
    2591       57260 :   val = F2x_valrem(f, &f);
    2592       57260 :   if (flag == 2 && val > 1) return NULL;
    2593       57260 :   V = F2x_factor_squarefree(f); lV = lg(V);
    2594       57260 :   if (flag == 2 && lV > 2) return NULL;
    2595             : 
    2596             :   /* to hold factors and exponents */
    2597       57260 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2598       57260 :   E = cgetg(d+1,t_VECSMALL);
    2599       57260 :   lfact = 1;
    2600       57260 :   if (val) {
    2601       18646 :     if (flag == 1)
    2602           0 :       t[1] = 1;
    2603             :     else
    2604       18646 :       gel(t,1) = polx_F2x(f[1]);
    2605       18646 :     E[1] = val; lfact++;
    2606             :   }
    2607             : 
    2608      222018 :   for (k=1; k<lV; k++)
    2609             :   {
    2610      164758 :     if (F2x_degree(gel(V, k))==0) continue;
    2611       52920 :     gel(t,lfact) = gel(V, k);
    2612       52920 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    2613       52920 :     if (flag == 2 && d != 1) return NULL;
    2614       52920 :     if (flag == 1)
    2615           0 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    2616       52920 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2617       52920 :     lfact += d;
    2618             :   }
    2619       57260 :   if (flag == 2) return gen_1; /* irreducible */
    2620       57260 :   y = FE_setlg(t,E, lfact);
    2621       57260 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2622       57260 :               : sort_factor_pol(y, cmpGuGu);
    2623             : }
    2624             : 
    2625             : static GEN
    2626      778333 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
    2627             : {
    2628      778333 :   long lfact, val, d = degpol(f), j, k, lV;
    2629             :   GEN y, E, t, V;
    2630             : 
    2631      778333 :   if (p == 2)
    2632             :   {
    2633       11970 :     GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
    2634       11970 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2635       11970 :     return F;
    2636             :   }
    2637      766363 :   if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
    2638      361318 :   val = Flx_valrem(f, &f);
    2639      361318 :   if (flag == 2 && val > 1) return NULL;
    2640      361318 :   V = Flx_factor_squarefree(f, p); lV = lg(V);
    2641      361318 :   if (flag == 2 && lV > 2) return NULL;
    2642             : 
    2643             :   /* to hold factors and exponents */
    2644      361318 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2645      361318 :   E = cgetg(d+1,t_VECSMALL);
    2646      361318 :   lfact = 1;
    2647      361318 :   if (val) {
    2648       14552 :     if (flag == 1)
    2649           0 :       t[1] = 1;
    2650             :     else
    2651       14552 :       gel(t,1) = polx_Flx(f[1]);
    2652       14552 :     E[1] = val; lfact++;
    2653             :   }
    2654             : 
    2655      906598 :   for (k=1; k<lV; k++)
    2656             :   {
    2657      545280 :     if (degpol(gel(V, k))==0) continue;
    2658      370782 :     gel(t,lfact) = Flx_normalize(gel(V, k), p);
    2659      370782 :     d = Flx_split_Berlekamp(&gel(t,lfact), p);
    2660      370782 :     if (flag == 2 && d != 1) return NULL;
    2661      370782 :     if (flag == 1)
    2662           0 :       for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
    2663      370782 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2664      370782 :     lfact += d;
    2665             :   }
    2666      361318 :   if (flag == 2) return gen_1; /* irreducible */
    2667      361318 :   y = FE_setlg(t,E,lfact);
    2668      361318 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2669      361318 :               : sort_factor_pol(y, cmpGuGu);
    2670             : }
    2671             : 
    2672             : /* f an FpX or an Flx */
    2673             : static GEN
    2674      851788 : FpX_Berlekamp_i(GEN f, GEN p, long flag)
    2675             : {
    2676      851788 :   long lfact, val, d = degpol(f), j, k, lV;
    2677             :   GEN y, E, t ,V;
    2678             : 
    2679      851788 :   if (typ(f) == t_VECSMALL)
    2680             :   {/* lgefint(p) == 3 */
    2681      851448 :     ulong pp = p[2];
    2682             :     GEN F;
    2683      851448 :     if (pp == 2) {
    2684      141358 :       F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
    2685      141358 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2686             :     } else {
    2687      710090 :       F = Flx_Berlekamp_i(f, pp, flag);
    2688      710090 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2689             :     }
    2690      851448 :     return F;
    2691             :   }
    2692             :   /* p is large (and odd) */
    2693         340 :   if (d <= 2) return FpX_factor_deg2(f, p, d, flag);
    2694         173 :   val = ZX_valrem(f, &f);
    2695         173 :   if (flag == 2 && val > 1) return NULL;
    2696         173 :   V = FpX_factor_Yun(f, p); lV = lg(V);
    2697         173 :   if (flag == 2 && lg(V) > 2) return NULL;
    2698             : 
    2699             :   /* to hold factors and exponents */
    2700         173 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2701         173 :   E = cgetg(d+1,t_VECSMALL);
    2702         173 :   lfact = 1;
    2703         173 :   if (val) {
    2704           2 :     if (flag == 1)
    2705           0 :       t[1] = 1;
    2706             :     else
    2707           2 :       gel(t,1) = pol_x(varn(f));
    2708           2 :     E[1] = val; lfact++;
    2709             :   }
    2710             : 
    2711         351 :   for (k=1; k<lV; k++)
    2712             :   {
    2713         178 :     if (degpol(gel(V,k))==0) continue;
    2714         177 :     gel(t,lfact) = FpX_normalize(gel(V, k), p);
    2715         177 :     d = FpX_split_Berlekamp(&gel(t,lfact), p);
    2716         177 :     if (flag == 2 && d != 1) return NULL;
    2717         177 :     if (flag == 1)
    2718           0 :       for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
    2719         177 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2720         177 :     lfact += d;
    2721             :   }
    2722         173 :   if (flag == 2) return gen_1; /* irreducible */
    2723         173 :   y = FE_setlg(t,E, lfact);
    2724         173 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2725         173 :               : sort_factor_pol(y, cmpii);
    2726             : }
    2727             : GEN
    2728      851683 : FpX_factor(GEN f, GEN p)
    2729             : {
    2730      851683 :   pari_sp av = avma;
    2731      851683 :   ZX_factmod_init(&f, p);
    2732      851683 :   return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
    2733             : }
    2734             : GEN
    2735      101642 : Flx_factor(GEN f, ulong p)
    2736             : {
    2737      101642 :   pari_sp av = avma;
    2738      101642 :   GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
    2739      101642 :   return gerepilecopy(av, F);
    2740             : }
    2741             : GEN
    2742          14 : F2x_factor(GEN f)
    2743             : {
    2744          14 :   pari_sp av = avma;
    2745          14 :   return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
    2746             : }
    2747             : 
    2748             : GEN
    2749         126 : factormod0(GEN f, GEN p, long flag)
    2750             : {
    2751         126 :   switch(flag)
    2752             :   {
    2753          84 :     case 0: return factmod(f,p);
    2754          42 :     case 1: return simplefactmod(f,p);
    2755           0 :     default: pari_err_FLAG("factormod");
    2756             :   }
    2757             :   return NULL; /* LCOV_EXCL_LINE */
    2758             : }

Generated by: LCOV version 1.11