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.8.0 lcov report (development 19207-2ed2f69) Lines: 1471 1614 91.1 %
Date: 2016-07-25 07:10:32 Functions: 123 132 93.2 %
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         266 : factmod_init(GEN *F, GEN p)
      32             : {
      33         266 :   if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
      34         266 :   if (signe(p) < 0) pari_err_PRIME("factmod",p);
      35         266 :   if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
      36         266 :   if (lgefint(p) == 3)
      37             :   {
      38         161 :     ulong pp = p[2];
      39         161 :     if (pp < 2) pari_err_PRIME("factmod", p);
      40         161 :     *F = RgX_to_Flx(*F, pp);
      41         161 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42             :   }
      43             :   else
      44             :   {
      45         105 :     *F = RgX_to_FpX(*F, p);
      46         105 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      47             :   }
      48         266 : }
      49             : /* as above, assume p prime and *F a ZX */
      50             : static void
      51      791075 : ZX_factmod_init(GEN *F, GEN p)
      52             : {
      53      791075 :   if (lgefint(p) == 3)
      54             :   {
      55      784546 :     ulong pp = p[2];
      56      784546 :     *F = ZX_to_Flx(*F, pp);
      57      784546 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      58             :   }
      59             :   else
      60             :   {
      61        6529 :     *F = FpX_red(*F, p);
      62        6529 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      63             :   }
      64      791075 : }
      65             : 
      66             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      67             : static GEN
      68          35 : all_roots_mod_p(ulong p, int not_0)
      69             : {
      70             :   GEN r;
      71             :   ulong i;
      72          35 :   if (not_0) {
      73          21 :     r = cgetg(p, t_VECSMALL);
      74          21 :     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          35 :   return r;
      80             : }
      81             : 
      82             : /* X^n - 1 */
      83             : static GEN
      84          20 : Flx_Xnm1(long sv, long n, ulong p)
      85             : {
      86          20 :   GEN t = cgetg(n+3, t_VECSMALL);
      87             :   long i;
      88          20 :   t[1] = sv;
      89          20 :   t[2] = p - 1;
      90          20 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      91          20 :   t[i] = 1; return t;
      92             : }
      93             : /* X^n + 1 */
      94             : static GEN
      95          13 : Flx_Xn1(long sv, long n, ulong p)
      96             : {
      97          13 :   GEN t = cgetg(n+3, t_VECSMALL);
      98             :   long i;
      99             :   (void) p;
     100          13 :   t[1] = sv;
     101          13 :   t[2] = 1;
     102          13 :   for (i = 3; i <= n+1; i++) t[i] = 0;
     103          13 :   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        6377 : Flx_root_mod_2(GEN f)
     161             : {
     162        6377 :   int z1, z0 = !(f[2] & 1);
     163             :   long i,n;
     164             :   GEN y;
     165             : 
     166        6377 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     167        6377 :   z1 = n & 1;
     168        6377 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     169        6377 :   if (z0) y[i++] = 0;
     170        6377 :   if (z1) y[i  ] = 1;
     171        6377 :   return y;
     172             : }
     173             : static ulong
     174          14 : Flx_oneroot_mod_2(GEN f)
     175             : {
     176             :   long i,n;
     177          14 :   if (!(f[2] & 1)) return 0;
     178          14 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     179          14 :   if (n & 1) return 1;
     180           7 :   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     1939761 : 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       60886 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
     196             : {
     197       60886 :   pari_sp av = avma;
     198             :   GEN y;
     199       60886 :   if (gpwrap)
     200          84 :     factmod_init(&f, pp);
     201             :   else
     202       60802 :     ZX_factmod_init(&f, pp);
     203       60886 :   switch(lg(f))
     204             :   {
     205          14 :     case 2: pari_err_ROOTS0("rootmod");
     206          35 :     case 3: avma = av; return cgetg(1,t_COL);
     207             :   }
     208       60837 :   if (typ(f) == t_VECSMALL)
     209             :   {
     210       58040 :     ulong p = pp[2];
     211       58040 :     if (p == 2)
     212        6377 :       y = Flx_root_mod_2(f);
     213             :     else
     214             :     {
     215       51663 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     216       51663 :       y = Roots(f, p);
     217             :     }
     218       58033 :     y = Flc_to_ZC(y);
     219             :   }
     220             :   else
     221        2797 :     y = FpX_roots_i(f, pp);
     222       60823 :   if (gpwrap) y = FpC_to_mod(y, pp);
     223       60823 :   return gerepileupto(av, y);
     224             : }
     225             : /* assume that f is a ZX an pp a prime */
     226             : GEN
     227       60802 : FpX_roots(GEN f, GEN pp)
     228       60802 : { 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           0 :   return NULL; /* not reached */
     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        6482 : FpX_quad_root(GEN x, GEN p, int unknown)
     257             : {
     258        6482 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     259             : 
     260        6482 :   if (equaliu(p, 2)) {
     261           0 :     if (!signe(b)) return c;
     262           0 :     return signe(c)? NULL: gen_1;
     263             :   }
     264        6482 :   D = subii(sqri(b), shifti(c,2));
     265        6482 :   D = remii(D,p);
     266        6482 :   if (unknown && kronecker(D,p) == -1) return NULL;
     267             : 
     268        6474 :   s = Fp_sqrt(D,p);
     269             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     270        6474 :   if (!s) return NULL;
     271        6464 :   return Fp_halve(Fp_sub(s,b, p), p);
     272             : }
     273             : static GEN
     274        2871 : FpX_otherroot(GEN x, GEN r, GEN p)
     275        2871 : { 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    17286352 : Fl_disc_bc(ulong b, ulong c, ulong p)
     280    17286352 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     281             : /* p > 2 */
     282             : static ulong
     283    17217146 : Flx_quad_root(GEN x, ulong p, int unknown)
     284             : {
     285    17217146 :   ulong s, b = x[3], c = x[2];
     286    17217146 :   ulong D = Fl_disc_bc(b, c, p);
     287    17207650 :   if (unknown && krouu(D,p) == -1) return p;
     288    11891412 :   s = Fl_sqrt(D,p);
     289    11896115 :   if (s==~0UL) return p;
     290    11896090 :   return Fl_halve(Fl_sub(s,b, p), p);
     291             : }
     292             : static ulong
     293    11461068 : Flx_otherroot(GEN x, ulong r, ulong p)
     294    11461068 : { 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     4552720 : split_init(struct split_t *S, long max)
     302             : {
     303     4552720 :   S->todo = vectrunc_init(max);
     304     4552727 :   S->done = vectrunc_init(max);
     305     4552446 : }
     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     4746814 : 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     4746791 : split_moveto_done(struct split_t *S, long i, GEN t)
     323             : {
     324     4746791 :   long n = lg(S->todo)-1;
     325     4746791 :   vectrunc_append(S->done, t);
     326     4747226 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     327     4747226 :   setlg(S->todo, n);
     328             : 
     329     4747271 : }
     330             : /* append t to done */
     331             : static void
     332      321323 : split_add_done(struct split_t *S, GEN t)
     333      321323 : { vectrunc_append(S->done, t); }
     334             : /* split todo[i] into a and b */
     335             : static void
     336      277756 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     337             : {
     338      277756 :   gel(S->todo, i) = a;
     339      277756 :   split_add(S, b);
     340      277760 : }
     341             : /* split todo[i] into a and b, moved to done */
     342             : static void
     343      302416 : split_done(struct split_t *S, long i, GEN a, GEN b)
     344             : {
     345      302416 :   split_moveto_done(S, i, a);
     346      302417 :   split_add_done(S, b);
     347      302420 : }
     348             : 
     349             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     350             : static GEN
     351        2797 : FpX_roots_i(GEN f, GEN p)
     352             : {
     353             :   GEN pol, pol0, a, q;
     354             :   struct split_t S;
     355             : 
     356        2797 :   split_init(&S, lg(f)-1);
     357        2797 :   settyp(S.done, t_COL);
     358        2797 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     359        2797 :   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        1603 :       GEN s, r = FpX_quad_root(f, p, 1);
     365        1603 :       if (r) {
     366        1603 :         split_add_done(&S, r);
     367        1603 :         s = FpX_otherroot(f,r, p);
     368             :         /* f not known to be square free yet */
     369        1603 :         if (!equalii(r, s)) split_add_done(&S, s);
     370             :       }
     371        1603 :       return sort(S.done);
     372             :     }
     373             :   }
     374             : 
     375        1173 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     376        1173 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     377        1173 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     378        1173 :   a = FpX_gcd(f,a, p);
     379        1173 :   if (!degpol(a)) return ZC_copy(S.done);
     380        1173 :   split_add(&S, FpX_normalize(a,p));
     381             : 
     382        1173 :   q = shifti(p,-1);
     383        1173 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     384        1173 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     385        2444 :   for (pol0[2] = 1;; pol0[2]++)
     386             :   {
     387        2444 :     long j, l = lg(S.todo);
     388        2444 :     if (l == 1) return sort(S.done);
     389        1278 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     390        2647 :     for (j = 1; j < l; j++)
     391             :     {
     392        1376 :       GEN c = gel(S.todo,j);
     393        1376 :       switch(degpol(c))
     394             :       { /* convert linear and quadratics to roots, try to split the rest */
     395             :         case 1:
     396          63 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     397          63 :           j--; l--; break;
     398             :         case 2: {
     399        1208 :           GEN r = FpX_quad_root(c, p, 0), s;
     400        1208 :           if (!r) pari_err_PRIME("polrootsmod",p);
     401        1201 :           s = FpX_otherroot(c,r, p);
     402        1201 :           split_done(&S, j, r, s);
     403        1201 :           j--; l--; break;
     404             :         }
     405             :         default: {
     406             :           /* b = pol^(p-1)/2 - 1 */
     407         105 :           GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
     408             :           long db;
     409         105 :           b = FpX_gcd(c,b, p); db = degpol(b);
     410         105 :           if (db && db < degpol(c))
     411             :           {
     412          98 :             b = FpX_normalize(b, p);
     413          98 :             c = FpX_div(c,b, p);
     414          98 :             split_todo(&S, j, b, c);
     415             :           }
     416             :         }
     417             :       }
     418             :     }
     419        1271 :   }
     420             : }
     421             : 
     422             : /* Assume f is normalized */
     423             : static ulong
     424      226208 : Flx_cubic_root(GEN ff, ulong p)
     425             : {
     426      226208 :   GEN f = Flx_normalize(ff,p);
     427      226208 :   ulong pi = get_Fl_red(p);
     428      226209 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     429      226209 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     430      226209 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     431      226208 :   ulong B = Fl_addmul_pre(t, Fl_sub(Fl_double(t2, p), b, p), c, p, pi);
     432      226209 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     433      226209 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     434      226209 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     435      226209 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     436      226209 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     437      226209 :   if (s!=~0UL)
     438             :   {
     439      130946 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     440      130946 :     if (p%3==2) /* 1 solutions */
     441       26276 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     442             :     else
     443             :     {
     444      104670 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     445      104670 :       if (vS1==~0UL) return p; /*0 solutions*/
     446             :       /*3 solutions*/
     447             :     }
     448       82258 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     449       82258 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     450             :   }
     451             :   else
     452             :   {
     453       95263 :     pari_sp av = avma;
     454       95263 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     455       95263 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     456             :     ulong Sa;
     457       95262 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     458       92356 :     Sa = vS1[1];
     459       92356 :     if (p%3==1) /*1 solutions*/
     460             :     {
     461       58826 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     462       58826 :       if (Fa!=P)
     463       38738 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     464             :     }
     465       92356 :     avma = av;
     466       92356 :     return Fl_sub(Fl_double(Sa,p),t,p);
     467             :   }
     468             : }
     469             : 
     470             : /* assume p > 2 prime */
     471             : static ulong
     472     1418307 : Flx_oneroot_i(GEN f, ulong p, long fl)
     473             : {
     474             :   GEN pol, a;
     475             :   ulong q;
     476             :   long da;
     477             : 
     478     1418307 :   if (Flx_val(f)) return 0;
     479     1418133 :   switch(degpol(f))
     480             :   {
     481        7099 :     case 1: return Fl_neg(f[2], p);
     482      358126 :     case 2: return Flx_quad_root(f, p, 1);
     483      207874 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     484             :   }
     485             : 
     486      845003 :   if (!fl)
     487             :   {
     488      817652 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     489      817652 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     490      817652 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     491      817652 :     a = Flx_gcd(f,a, p);
     492       27351 :   } else a = f;
     493      845003 :   da = degpol(a);
     494      844998 :   if (!da) return p;
     495      383771 :   a = Flx_normalize(a,p);
     496             : 
     497      383796 :   q = p >> 1;
     498      383796 :   pol = polx_Flx(f[1]);
     499      555234 :   for(pol[2] = 1;; pol[2]++)
     500             :   {
     501      555234 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     502      555233 :     switch(da)
     503             :     {
     504      137652 :       case 1: return Fl_neg(a[2], p);
     505      227875 :       case 2: return Flx_quad_root(a, p, 0);
     506       18334 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     507             :       default: {
     508      171372 :         GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
     509             :         long db;
     510      171502 :         b = Flx_gcd(a,b, p); db = degpol(b);
     511      171492 :         if (db && db < da)
     512             :         {
     513      157467 :           b = Flx_normalize(b, p);
     514      157476 :           if (db <= (da >> 1)) {
     515      103134 :             a = b;
     516      103134 :             da = db;
     517             :           } else {
     518       54342 :             a = Flx_div(a,b, p);
     519       54340 :             da -= db;
     520             :           }
     521             :         }
     522             :       }
     523             :     }
     524      171499 :   }
     525             : }
     526             : 
     527             : /* assume p > 2 prime */
     528             : static GEN
     529        3649 : FpX_oneroot_i(GEN f, GEN p)
     530             : {
     531             :   GEN pol, pol0, a, q;
     532             :   long da;
     533             : 
     534        3649 :   if (ZX_val(f)) return gen_0;
     535        3642 :   switch(degpol(f))
     536             :   {
     537         140 :     case 1: return subii(p, gel(f,2));
     538        3432 :     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     1383292 : Flx_oneroot(GEN f, ulong p)
     581             : {
     582     1383292 :   pari_sp av = avma;
     583             :   ulong r;
     584     1383292 :   switch(lg(f))
     585             :   {
     586           0 :     case 2: return 0;
     587           1 :     case 3: avma = av; return p;
     588             :   }
     589     1383291 :   if (p == 2) return Flx_oneroot_mod_2(f);
     590     1383291 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     591     1383291 :   avma = av; return r;
     592             : }
     593             : 
     594             : ulong
     595       27666 : Flx_oneroot_split(GEN f, ulong p)
     596             : {
     597       27666 :   pari_sp av = avma;
     598             :   ulong r;
     599       27666 :   switch(lg(f))
     600             :   {
     601           0 :     case 2: return 0;
     602           0 :     case 3: avma = av; return p;
     603             :   }
     604       27666 :   if (p == 2) return Flx_oneroot_mod_2(f);
     605       27666 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     606       27742 :   avma = av; return r;
     607             : }
     608             : 
     609             : /* assume that p is prime */
     610             : GEN
     611       11004 : FpX_oneroot(GEN f, GEN pp) {
     612       11004 :   pari_sp av = avma;
     613       11004 :   ZX_factmod_init(&f, pp);
     614       11004 :   switch(lg(f))
     615             :   {
     616           0 :     case 2: avma = av; return gen_0;
     617           0 :     case 3: avma = av; return NULL;
     618             :   }
     619       11004 :   if (typ(f) == t_VECSMALL)
     620             :   {
     621        7355 :     ulong r, p = pp[2];
     622        7355 :     if (p == 2)
     623          14 :       r = Flx_oneroot_mod_2(f);
     624             :     else
     625        7341 :       r = Flx_oneroot_i(f, p, 0);
     626        7355 :     avma = av;
     627        7355 :     return (r == p)? NULL: utoi(r);
     628             :   }
     629        3649 :   f = FpX_oneroot_i(f, pp);
     630        3649 :   if (!f) { avma = av; return NULL; }
     631        3649 :   return gerepileuptoint(av, f);
     632             : }
     633             : 
     634             : /*******************************************************************/
     635             : /*                                                                 */
     636             : /*                     FACTORISATION MODULO p                      */
     637             : /*                                                                 */
     638             : /*******************************************************************/
     639             : 
     640             : /* Functions giving information on the factorisation. */
     641             : 
     642             : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
     643             : GEN
     644          55 : FpX_Berlekamp_ker(GEN u, GEN p)
     645             : {
     646          55 :   pari_sp ltop=avma;
     647          55 :   long j,N = degpol(u);
     648          55 :   GEN Q  = FpX_matFrobenius(u, p);
     649        2139 :   for (j=1; j<=N; j++)
     650        2084 :     gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
     651          55 :   return gerepileupto(ltop, FpM_ker(Q,p));
     652             : }
     653             : 
     654             : GEN
     655       22554 : F2x_Berlekamp_ker(GEN u)
     656             : {
     657       22554 :   pari_sp ltop=avma;
     658       22554 :   long j,N = F2x_degree(u);
     659             :   GEN Q;
     660             :   pari_timer T;
     661       22554 :   timer_start(&T);
     662       22554 :   Q = F2x_matFrobenius(u);
     663      121730 :   for (j=1; j<=N; j++)
     664       99176 :     F2m_flip(Q,j,j);
     665       22554 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     666       22554 :   Q = F2m_ker_sp(Q,0);
     667       22554 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     668       22554 :   return gerepileupto(ltop,Q);
     669             : }
     670             : 
     671             : GEN
     672      207101 : Flx_Berlekamp_ker(GEN u, ulong l)
     673             : {
     674      207101 :   pari_sp ltop=avma;
     675      207101 :   long j,N = degpol(u);
     676             :   GEN Q;
     677             :   pari_timer T;
     678      207101 :   timer_start(&T);
     679      207101 :   Q  = Flx_matFrobenius(u, l);
     680      973450 :   for (j=1; j<=N; j++)
     681      766349 :     coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
     682      207101 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     683      207101 :   Q = Flm_ker_sp(Q,l,0);
     684      207101 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     685      207101 :   return gerepileupto(ltop,Q);
     686             : }
     687             : 
     688             : /* product of terms of degree 1 in factorization of f */
     689             : GEN
     690      140469 : FpX_split_part(GEN f, GEN p)
     691             : {
     692      140469 :   long n = degpol(f);
     693      140469 :   GEN z, X = pol_x(varn(f));
     694      140469 :   if (n <= 1) return f;
     695      140441 :   f = FpX_red(f, p);
     696      140441 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     697      140441 :   return FpX_gcd(z,f,p);
     698             : }
     699             : 
     700             : /* Compute the number of roots in Fp without counting multiplicity
     701             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     702             : long
     703       98539 : FpX_nbroots(GEN f, GEN p)
     704             : {
     705       98539 :   pari_sp av = avma;
     706       98539 :   GEN z = FpX_split_part(f, p);
     707       98539 :   avma = av; return degpol(z);
     708             : }
     709             : 
     710             : int
     711           0 : FpX_is_totally_split(GEN f, GEN p)
     712             : {
     713           0 :   long n=degpol(f);
     714           0 :   pari_sp av = avma;
     715           0 :   if (n <= 1) return 1;
     716           0 :   if (cmpui(n, p) > 0) return 0;
     717           0 :   f = FpX_red(f, p);
     718           0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     719             : }
     720             : 
     721             : long
     722     2548963 : Flx_nbroots(GEN f, ulong p)
     723             : {
     724     2548963 :   long n = degpol(f);
     725     2548963 :   pari_sp av = avma;
     726             :   GEN z;
     727     2548963 :   if (n <= 1) return n;
     728     2548788 :   if (n == 2)
     729             :   {
     730             :     ulong D;
     731        6811 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     732        5978 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     733        5978 :     return 1 + krouu(D,p);
     734             :   }
     735     2541977 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     736     2541977 :   z = Flx_gcd(z, f, p);
     737     2541977 :   avma = av; return degpol(z);
     738             : }
     739             : 
     740             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     741             : static GEN
     742        4102 : FpX_ddf(GEN T, GEN XP, GEN p)
     743             : {
     744        4102 :   pari_sp av = avma;
     745             :   GEN b, g, h, F, f, Tr, xq;
     746             :   long i, j, n, v;
     747             :   long B, l, m;
     748             :   pari_timer ti;
     749        4102 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     750        4102 :   if (n == 0) return cgetg(1, t_VEC);
     751        4102 :   if (n == 1) return mkvec(get_FpX_mod(T));
     752        4095 :   B = n/2;
     753        4095 :   l = usqrt(B);
     754        4095 :   m = (B+l-1)/l;
     755        4095 :   T = FpX_get_red(T, p);
     756        4095 :   b = cgetg(l+2, t_VEC);
     757        4095 :   gel(b, 1) = pol_x(v);
     758        4095 :   gel(b, 2) = XP;
     759        4095 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     760        4095 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     761        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
     762       10661 :   for (i = 3; i <= l+1; i++)
     763        6566 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     764        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
     765        4095 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     766        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
     767        4095 :   g = cgetg(m+1, t_VEC);
     768        4095 :   gel(g, 1) = gel(xq, 2);
     769       14931 :   for(i = 2; i <= m; i++)
     770       10836 :     gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     771        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
     772        4095 :   h = cgetg(m+1, t_VEC);
     773       19026 :   for (j = 1; j <= m; j++)
     774             :   {
     775       14931 :     pari_sp av = avma;
     776       14931 :     GEN gj = gel(g, j);
     777       14931 :     GEN e = FpX_sub(gj, gel(b, 1), p);
     778       42938 :     for (i = 2; i <= l; i++)
     779       28007 :       e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
     780       14931 :     gel(h, j) = gerepileupto(av, e);
     781             :   }
     782        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
     783        4095 :   Tr = get_FpX_mod(T);
     784        4095 :   F = cgetg(m+1, t_VEC);
     785       19026 :   for (j = 1; j <= m; j++)
     786             :   {
     787       14931 :     gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
     788       14931 :     Tr = FpX_div(Tr, gel(F,j), p);
     789             :   }
     790        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
     791        4095 :   f = const_vec(n, pol_1(v));
     792       19026 :   for (j = 1; j <= m; j++)
     793             :   {
     794       14931 :     GEN e = gel(F, j);
     795       16919 :     for (i=l-1; i >= 0; i--)
     796             :     {
     797       16919 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     798       16919 :       if (degpol(u))
     799             :       {
     800        2709 :         gel(f, l*j-i) = u;
     801        2709 :         e = FpX_div(e, u, p);
     802             :       }
     803       16919 :       if (!degpol(e)) break;
     804             :     }
     805             :   }
     806        4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
     807        4095 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     808        4095 :   return gerepilecopy(av, f);
     809             : }
     810             : 
     811             : static void
     812           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     813             : {
     814           0 :   long n = degpol(Tp), r = n/d;
     815             :   GEN T, f, ff;
     816             :   GEN p2;
     817           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     818           0 :   p2 = shifti(p,-1);
     819           0 :   T = FpX_get_red(Tp, p);
     820           0 :   XP = FpX_rem(XP, T, p);
     821             :   while (1)
     822             :   {
     823           0 :     pari_sp btop = avma;
     824             :     long i;
     825           0 :     GEN g = random_FpX(n, varn(Tp), p);
     826           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     827           0 :     if (signe(t) == 0) continue;
     828           0 :     for(i=1; i<=10; i++)
     829             :     {
     830           0 :       pari_sp btop2 = avma;
     831           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     832           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     833           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     834           0 :       avma = btop2;
     835             :     }
     836           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     837           0 :     avma = btop;
     838           0 :   }
     839           0 :   f = FpX_normalize(f, p);
     840           0 :   ff = FpX_div(Tp, f ,p);
     841           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     842           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     843             : }
     844             : 
     845             : static void
     846         210 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     847             : {
     848             :   pari_sp av;
     849         210 :   GEN Tp = get_FpX_mod(T);
     850         210 :   long n = degpol(hp), vT = varn(Tp);
     851             :   GEN u1, u2, f1, f2;
     852             :   GEN R, h;
     853         210 :   h = FpX_get_red(hp, p);
     854         210 :   t = FpX_rem(t, T, p);
     855         210 :   av = avma;
     856             :   do
     857             :   {
     858         343 :     avma = av;
     859         343 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     860         343 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     861         343 :   } while (degpol(u1)==0 || degpol(u1)==n);
     862         210 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     863         210 :   f1 = FpX_normalize(f1, p);
     864         210 :   u2 = FpX_div(hp, u1, p);
     865         210 :   f2 = FpX_div(Tp, f1, p);
     866         210 :   if (degpol(u1)==1)
     867         105 :     gel(V, idx) = f1;
     868             :   else
     869         105 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     870         210 :   idx += degpol(f1)/d;
     871         210 :   if (degpol(u2)==1)
     872         119 :     gel(V, idx) = f2;
     873             :   else
     874          91 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     875         210 : }
     876             : 
     877             : static void
     878          14 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     879             : {
     880          14 :   long n = degpol(Tp), r = n/d, vT = varn(Tp);
     881             :   GEN T, h, t;
     882             :   pari_timer ti;
     883          28 :   if (r==1) { gel(V, idx) = Tp; return; }
     884          14 :   T = FpX_get_red(Tp, p);
     885          14 :   XP = FpX_rem(XP, T, p);
     886          14 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     887             :   do
     888             :   {
     889          14 :     GEN g = random_FpX(n, vT, p);
     890          14 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     891          14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     892          14 :     h = FpXQ_minpoly(t, T, p);
     893          14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     894          14 :   } while (degpol(h) != r);
     895          14 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     896             : }
     897             : 
     898             : static GEN
     899          14 : FpX_factor_Shoup(GEN T, GEN p)
     900             : {
     901          14 :   long i, n, s = 0;
     902             :   GEN XP, D, V;
     903          14 :   long e = expi(p);
     904             :   pari_timer ti;
     905          14 :   n = get_FpX_degree(T);
     906          14 :   T = FpX_get_red(T, p);
     907          14 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     908          14 :   XP = FpX_Frobenius(T, p);
     909          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     910          14 :   D = FpX_ddf(T, XP, p);
     911          14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     912         301 :   for (i = 1; i <= n; i++)
     913         287 :     s += degpol(gel(D,i))/i;
     914          14 :   V = cgetg(s+1, t_COL);
     915         301 :   for (i = 1, s = 1; i <= n; i++)
     916             :   {
     917         287 :     GEN Di = gel(D,i);
     918         287 :     long ni = degpol(Di), ri = ni/i;
     919         287 :     if (ni == 0) continue;
     920          21 :     Di = FpX_normalize(Di, p);
     921          21 :     if (ni == i) { gel(V, s++) = Di; continue; }
     922          14 :     if (ri <= e*expu(e))
     923          14 :       FpX_edf(Di, XP, i, p, V, s);
     924             :     else
     925           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     926          14 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     927          14 :     s += ri;
     928             :   }
     929          14 :   return V;
     930             : }
     931             : 
     932             : static GEN
     933           0 : FpX_simplefact_Shoup(GEN T, GEN p)
     934             : {
     935           0 :   long i, n, s = 0, j = 1, k;
     936             :   GEN XP, D, V;
     937             :   pari_timer ti;
     938           0 :   n = get_FpX_degree(T);
     939           0 :   T = FpX_get_red(T, p);
     940           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     941           0 :   XP = FpX_Frobenius(T, p);
     942           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     943           0 :   D = FpX_ddf(T, XP, p);
     944           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     945           0 :   for (i = 1; i <= n; i++)
     946           0 :     s += degpol(gel(D,i))/i;
     947           0 :   V = cgetg(s+1, t_VEC);
     948           0 :   for (i = 1; i <= n; i++)
     949             :   {
     950           0 :     long ni = degpol(gel(D,i)), ri = ni/i;
     951           0 :     if (ni == 0) continue;
     952           0 :     for (k = 1; k <= ri; k++)
     953           0 :       gel(V, j++) = utoi(i);
     954             :   }
     955           0 :   return V;
     956             : }
     957             : 
     958             : /* Yun algorithm: Assume p > degpol(T) */
     959             : 
     960             : static GEN
     961           7 : FpX_factor_Yun(GEN T, GEN p)
     962             : {
     963           7 :   pari_sp av = avma;
     964           7 :   long n = degpol(T);
     965           7 :   long i = 1;
     966           7 :   GEN d = FpX_deriv(T, p);
     967             :   GEN a, b, c;
     968           7 :   GEN V = cgetg(n+1,t_VEC);
     969           7 :   a = FpX_gcd(T, d, p);
     970           7 :   if (degpol(a) == 0) return mkvec(T);
     971           7 :   b = FpX_div(T, a, p);
     972             :   do
     973             :   {
     974          14 :     c = FpX_div(d, a, p);
     975          14 :     d = FpX_sub(c, FpX_deriv(b, p), p);
     976          14 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
     977          14 :     gel(V, i++) = a;
     978          14 :     b = FpX_div(b, a, p);
     979          14 :   } while (degpol(b));
     980           7 :   setlg(V, i);
     981           7 :   return gerepilecopy(av, V);
     982             : }
     983             : 
     984             : static GEN
     985           7 : FpX_factor_Cantor(GEN T, GEN p)
     986             : {
     987             :   GEN E, F, M, V;
     988             :   long i, j, s, l;
     989           7 :   V = FpX_factor_Yun(T, p);
     990           7 :   l = lg(V);
     991          21 :   for (s=0, i=1; i < l; i++)
     992          14 :     s += !!degpol(gel(V,i));
     993           7 :   F = cgetg(s+1, t_VEC);
     994           7 :   E = cgetg(s+1, t_VEC);
     995          21 :   for (i=1, j=1; i < l; i++)
     996          14 :     if (degpol(gel(V,i)))
     997             :     {
     998          14 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
     999          14 :       gel(F, j) = Fj;
    1000          14 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1001          14 :       j++;
    1002             :     }
    1003           7 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1004           7 :   return sort_factor_pol(M, cmpii);
    1005             : }
    1006             : 
    1007             : static GEN
    1008           0 : FpX_simplefact_Cantor(GEN T, GEN p)
    1009             : {
    1010             :   GEN E, F, M, V;
    1011             :   long i, j, s, l;
    1012           0 :   V = FpX_factor_Yun(get_FpX_mod(T), p);
    1013           0 :   l = lg(V);
    1014           0 :   for (s=0, i=1; i < l; i++)
    1015           0 :     s += !!degpol(gel(V,i));
    1016           0 :   F = cgetg(s+1, t_VEC);
    1017           0 :   E = cgetg(s+1, t_VEC);
    1018           0 :   for (i=1, j=1; i < l; i++)
    1019           0 :     if (degpol(gel(V,i)))
    1020             :     {
    1021           0 :       GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
    1022           0 :       gel(F, j) = Fj;
    1023           0 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1024           0 :       j++;
    1025             :     }
    1026           0 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1027           0 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1028             : }
    1029             : 
    1030             : static int
    1031           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1032             : {
    1033           0 :   pari_sp av = avma;
    1034             :   pari_timer ti;
    1035             :   long n, d;
    1036           0 :   GEN T = get_FpX_mod(Tp);
    1037           0 :   GEN dT = FpX_deriv(T, p);
    1038             :   GEN XP, D;
    1039           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1040           0 :   n = get_FpX_degree(T);
    1041           0 :   T = FpX_get_red(Tp, p);
    1042           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1043           0 :   XP = FpX_Frobenius(T, p);
    1044           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1045           0 :   D = FpX_ddf(T, XP, p);
    1046           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1047           0 :   d = degpol(gel(D, n));
    1048           0 :   avma = av; return d==n;
    1049             : }
    1050             : 
    1051             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1052             : 
    1053             : /*Assume that p is large and odd*/
    1054             : static GEN
    1055          35 : FpX_factcantor_i(GEN f, GEN pp, long flag)
    1056             : {
    1057          35 :   long d = degpol(f);
    1058          35 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1059           7 :   switch(flag)
    1060             :   {
    1061           7 :     default: return FpX_factor_Cantor(f, pp);
    1062           0 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1063           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1064             :   }
    1065             : }
    1066             : 
    1067             : long
    1068        4088 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1069             : {
    1070        4088 :   pari_sp av = avma;
    1071        4088 :   GEN ddf = FpX_ddf(T, XP, p);
    1072        4088 :   long l = lg(ddf), i, s=0;
    1073       86716 :   for(i = 1; i < l; i++)
    1074       82628 :     s += degpol(gel(ddf,i))/i;
    1075        4088 :   avma = av; return s;
    1076             : }
    1077             : 
    1078             : long
    1079          28 : FpX_nbfact(GEN T, GEN p)
    1080             : {
    1081          28 :   pari_sp av = avma;
    1082          28 :   GEN XP = FpX_Frobenius(T, p);
    1083          28 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1084          28 :   avma = av; return n;
    1085             : }
    1086             : 
    1087             : /* p > 2 */
    1088             : static GEN
    1089           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1090             : {
    1091           7 :   switch(d)
    1092             :   {
    1093             :     case -1:
    1094           0 :     case 0: return NULL;
    1095           0 :     case 1: return gen_1;
    1096             :   }
    1097           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1098             : }
    1099             : /* p > 2 */
    1100             : static GEN
    1101          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1102             : {
    1103          14 :   switch(d)
    1104             :   {
    1105           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1106           0 :     case 0: return trivial_fact();
    1107           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1108             :   }
    1109          14 :   switch(FpX_quad_factortype(f, p)) {
    1110           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1111           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1112           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1113             :   }
    1114             : }
    1115             : 
    1116             : GEN
    1117          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1118             : GEN
    1119         763 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1120             : 
    1121             : /* Mod(0,p) * x, where x is f's main variable */
    1122             : static GEN
    1123          14 : Mod0pX(GEN f, GEN p)
    1124          14 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
    1125             : static GEN
    1126          14 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
    1127             : 
    1128             : /* not gerepile safe */
    1129             : static GEN
    1130          74 : FpX_factor_2(GEN f, GEN p, long d)
    1131             : {
    1132             :   GEN r, s, R, S;
    1133             :   long v;
    1134             :   int sgn;
    1135          74 :   switch(d)
    1136             :   {
    1137           0 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1138           2 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1139          10 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1140             :   }
    1141          62 :   r = FpX_quad_root(f, p, 1);
    1142          62 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1143          53 :   v = varn(f);
    1144          53 :   s = FpX_otherroot(f, r, p);
    1145          53 :   if (signe(r)) r = subii(p, r);
    1146          53 :   if (signe(s)) s = subii(p, s);
    1147          53 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1148          53 :   R = deg1pol_shallow(gen_1, r, v);
    1149          53 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1150          46 :   S = deg1pol_shallow(gen_1, s, v);
    1151          46 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1152             : }
    1153             : static GEN
    1154          95 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1155             : {
    1156          95 :   switch(flag) {
    1157           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1158          14 :     case 1: return FpX_degfact_2(f, p, d);
    1159          74 :     default: return FpX_factor_2(f, p, d);
    1160             :   }
    1161             : }
    1162             : 
    1163             : static int
    1164       28651 : F2x_quad_factortype(GEN x)
    1165       28651 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1166             : 
    1167             : static GEN
    1168           7 : F2x_is_irred_2(GEN f, long d)
    1169           7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1170             : 
    1171             : static GEN
    1172         854 : F2x_degfact_2(GEN f, long d)
    1173             : {
    1174         854 :   if (!d) return trivial_fact();
    1175         854 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1176         735 :   switch(F2x_quad_factortype(f)) {
    1177         133 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1178         147 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1179         455 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1180             :   }
    1181             : }
    1182             : 
    1183             : static GEN
    1184       73262 : F2x_factor_2(GEN f, long d)
    1185             : {
    1186       73262 :   long v = f[1];
    1187       73262 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1188       73262 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1189       65737 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1190       20083 :   switch(F2x_quad_factortype(f))
    1191             :   {
    1192        3955 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1193       13559 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1194        2569 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1195             :   }
    1196             : }
    1197             : static GEN
    1198       74123 : F2x_factor_deg2(GEN f, long d, long flag)
    1199             : {
    1200       74123 :   switch(flag) {
    1201           7 :     case 2: return F2x_is_irred_2(f, d);
    1202         854 :     case 1: return F2x_degfact_2(f, d);
    1203       73262 :     default: return F2x_factor_2(f, d);
    1204             :   }
    1205             : }
    1206             : 
    1207             : static void
    1208          19 : split_squares(struct split_t *S, GEN g, ulong p)
    1209             : {
    1210          19 :   ulong q = p >> 1;
    1211          19 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1212          19 :   if (degpol(a) < 0)
    1213             :   {
    1214             :     ulong i;
    1215           6 :     split_add_done(S, (GEN)1);
    1216           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1217             :   } else {
    1218          13 :     (void)Flx_valrem(a, &a);
    1219          13 :     if (degpol(a))
    1220             :     {
    1221          13 :       a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
    1222          13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1223             :     }
    1224             :   }
    1225          19 : }
    1226             : static void
    1227          19 : split_nonsquares(struct split_t *S, GEN g, ulong p)
    1228             : {
    1229          19 :   ulong q = p >> 1;
    1230          19 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1231          19 :   if (degpol(a) < 0)
    1232             :   {
    1233           6 :     ulong i, z = Fl_nonsquare(p);
    1234           6 :     split_add_done(S, (GEN)z);
    1235           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1236             :   } else {
    1237          13 :     (void)Flx_valrem(a, &a);
    1238          13 :     if (degpol(a))
    1239             :     {
    1240          13 :       a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
    1241          13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1242             :     }
    1243             :   }
    1244          19 : }
    1245             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1246             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1247             : static int
    1248     4549543 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1249             : {
    1250     4549543 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1251     4549582 :   long d = degpol(g);
    1252     4549561 :   if (d < 0) return 0;
    1253     4549494 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1254     4549681 :   if (!d) return 1;
    1255     4534141 :   if (p <= 1.4 * (ulong)d) {
    1256             :     /* small p; split further using x^((p-1)/2) +/- 1.
    1257             :      * 30% degree drop makes the extra gcd worth it. */
    1258          19 :     split_squares(S, g, p);
    1259          19 :     split_nonsquares(S, g, p);
    1260             :   } else { /* large p; use x^(p-1) - 1 directly */
    1261     4534122 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1262     4533811 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1263     4533811 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1264     4533598 :     g = Flx_gcd(g,a, p);
    1265     4533509 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1266             :   }
    1267     4534229 :   return 1;
    1268             : }
    1269             : 
    1270             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1271             : static GEN
    1272    20601399 : Flx_roots_i(GEN f, ulong p)
    1273             : {
    1274             :   GEN pol, g;
    1275    20601399 :   long v = Flx_valrem(f, &g);
    1276             :   ulong q;
    1277             :   struct split_t S;
    1278             : 
    1279             :   /* optimization: test for small degree first */
    1280    20606178 :   switch(degpol(g))
    1281             :   {
    1282             :     case 1: {
    1283       23548 :       ulong r = p - g[2];
    1284       23548 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1285             :     }
    1286             :     case 2: {
    1287    16032127 :       ulong r = Flx_quad_root(g, p, 1), s;
    1288    16045699 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1289    10964489 :       s = Flx_otherroot(g,r, p);
    1290    10969838 :       if (r < s)
    1291     2737781 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1292     8232057 :       else if (r > s)
    1293     8231896 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1294             :       else
    1295         161 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1296             :     }
    1297             :   }
    1298     4549947 :   q = p >> 1;
    1299     4549947 :   split_init(&S, lg(f)-1);
    1300     4549640 :   settyp(S.done, t_VECSMALL);
    1301     4549640 :   if (v) split_add_done(&S, (GEN)0);
    1302     4549640 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1303          35 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1304     4549448 :   pol = polx_Flx(f[1]);
    1305     9365798 :   for (pol[2]=1; ; pol[2]++)
    1306             :   {
    1307     9365798 :     long j, l = lg(S.todo);
    1308     9365798 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1309     4816187 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1310     9921850 :     for (j = 1; j < l; j++)
    1311             :     {
    1312     5105947 :       GEN c = gel(S.todo,j);
    1313     5105947 :       switch(degpol(c))
    1314             :       {
    1315             :         case 1:
    1316     4444348 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1317     4444649 :           j--; l--; break;
    1318             :         case 2: {
    1319      301221 :           ulong r = Flx_quad_root(c, p, 0), s;
    1320      301225 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1321      301218 :           s = Flx_otherroot(c,r, p);
    1322      301215 :           split_done(&S, j, (GEN)r, (GEN)s);
    1323      301220 :           j--; l--; break;
    1324             :         }
    1325             :         default: {
    1326      360357 :           GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
    1327             :           long db;
    1328      360345 :           b = Flx_gcd(c,b, p); db = degpol(b);
    1329      360355 :           if (db && db < degpol(c))
    1330             :           {
    1331      277662 :             b = Flx_normalize(b, p);
    1332      277664 :             c = Flx_div(c,b, p);
    1333      277659 :             split_todo(&S, j, b, c);
    1334             :           }
    1335             :         }
    1336             :       }
    1337             :     }
    1338     4815903 :   }
    1339             : }
    1340             : 
    1341             : GEN
    1342    20548922 : Flx_roots(GEN f, ulong p)
    1343             : {
    1344    20548922 :   pari_sp av = avma;
    1345    20548922 :   switch(lg(f))
    1346             :   {
    1347           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1348           0 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1349             :   }
    1350    20559007 :   if (p == 2) return Flx_root_mod_2(f);
    1351    20559007 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1352             : }
    1353             : 
    1354             : /* assume x reduced mod p, monic. */
    1355             : static int
    1356       75747 : Flx_quad_factortype(GEN x, ulong p)
    1357             : {
    1358       75747 :   ulong b = x[3], c = x[2];
    1359       75747 :   return krouu(Fl_disc_bc(b, c, p), p);
    1360             : }
    1361             : static GEN
    1362           0 : Flx_is_irred_2(GEN f, ulong p, long d)
    1363             : {
    1364           0 :   if (!d) return NULL;
    1365           0 :   if (d == 1) return gen_1;
    1366           0 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1367             : }
    1368             : static GEN
    1369       77959 : Flx_degfact_2(GEN f, ulong p, long d)
    1370             : {
    1371       77959 :   if (!d) return trivial_fact();
    1372       77959 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1373       75747 :   switch(Flx_quad_factortype(f, p)) {
    1374       35469 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1375       39340 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1376         938 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1377             :   }
    1378             : }
    1379             : /* p > 2 */
    1380             : static GEN
    1381      340684 : Flx_factor_2(GEN f, ulong p, long d)
    1382             : {
    1383             :   ulong r, s;
    1384             :   GEN R,S;
    1385      340684 :   long v = f[1];
    1386      340684 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1387      340684 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1388      327015 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1389      229550 :   r = Flx_quad_root(f, p, 1);
    1390      229550 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1391      136564 :   s = Flx_otherroot(f, r, p);
    1392      136564 :   r = Fl_neg(r, p);
    1393      136564 :   s = Fl_neg(s, p);
    1394      136564 :   if (s < r) lswap(s,r);
    1395      136564 :   R = mkvecsmall3(v,r,1);
    1396      136564 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1397      122081 :   S = mkvecsmall3(v,s,1);
    1398      122081 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1399             : }
    1400             : static GEN
    1401      418643 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1402             : {
    1403      418643 :   switch(flag) {
    1404           0 :     case 2: return Flx_is_irred_2(f, p, d);
    1405       77959 :     case 1: return Flx_degfact_2(f, p, d);
    1406      340684 :     default: return Flx_factor_2(f, p, d);
    1407             :   }
    1408             : }
    1409             : 
    1410             : void
    1411       14945 : F2xV_to_FlxV_inplace(GEN v)
    1412             : {
    1413             :   long i;
    1414       14945 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1415       14945 : }
    1416             : void
    1417      601707 : FlxV_to_ZXV_inplace(GEN v)
    1418             : {
    1419             :   long i;
    1420      601707 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1421      601707 : }
    1422             : void
    1423      119936 : F2xV_to_ZXV_inplace(GEN v)
    1424             : {
    1425             :   long i;
    1426      119936 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1427      119936 : }
    1428             : 
    1429             : /* Adapted from Shoup NTL */
    1430             : static GEN
    1431       11186 : F2x_factor_squarefree(GEN f)
    1432             : {
    1433       11186 :   pari_sp av = avma;
    1434             :   GEN r, t, v, tv;
    1435       11186 :   long q, n = F2x_degree(f);
    1436       11186 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1437       15274 :   for(q = 1;;q *= 2)
    1438             :   {
    1439       15274 :     r = F2x_gcd(f, F2x_deriv(f));
    1440       15274 :     if (F2x_degree(r) == 0)
    1441             :     {
    1442       10486 :       gel(u, q) = f;
    1443       10486 :       break;
    1444             :     }
    1445        4788 :     t = F2x_div(f, r);
    1446        4788 :     if (F2x_degree(t) > 0)
    1447             :     {
    1448             :       long j;
    1449        5733 :       for(j = 1;;j++)
    1450             :       {
    1451        5733 :         v = F2x_gcd(r, t);
    1452        5733 :         tv = F2x_div(t, v);
    1453        5733 :         if (F2x_degree(tv) > 0)
    1454        4207 :           gel(u, j*q) = tv;
    1455        5733 :         if (F2x_degree(v) <= 0) break;
    1456        2142 :         r = F2x_div(r, v);
    1457        2142 :         t = v;
    1458        2142 :       }
    1459        3591 :       if (F2x_degree(r) == 0) break;
    1460             :     }
    1461        4088 :     f = F2x_sqrt(r);
    1462        4088 :   }
    1463       11186 :   return gerepilecopy(av, u);
    1464             : }
    1465             : 
    1466             : static GEN
    1467       14791 : F2x_ddf_simple(GEN T, GEN XP)
    1468             : {
    1469       14791 :   pari_sp av = avma, av2;
    1470             :   GEN f, z, Tr, X;
    1471       14791 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1472       14791 :   if (n == 0) return cgetg(1, t_VEC);
    1473       14791 :   if (n == 1) return mkvec(T);
    1474       11235 :   z = XP; Tr = T; X = polx_F2x(v);
    1475       11235 :   f = const_vec(n, pol1_F2x(v));
    1476       11235 :   av2 = avma;
    1477       84028 :   for (j = 1; j <= B; j++)
    1478             :   {
    1479       77560 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1480       77560 :     if (F2x_degree(u))
    1481             :     {
    1482       20776 :       gel(f, j) = u;
    1483       20776 :       Tr = F2x_div(Tr, u);
    1484       20776 :       av2 = avma;
    1485       56784 :     } else z = gerepileuptoleaf(av2, z);
    1486       77560 :     if (!F2x_degree(Tr)) break;
    1487       72793 :     z = F2xq_sqr(z, Tr);
    1488             :   }
    1489       11235 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1490       11235 :   return gerepilecopy(av, f);
    1491             : }
    1492             : 
    1493             : static GEN
    1494        4424 : F2xq_frobtrace(GEN a, long d, GEN T)
    1495             : {
    1496        4424 :   pari_sp av = avma;
    1497             :   long i;
    1498        4424 :   GEN x = a;
    1499       27951 :   for(i=1; i<d; i++)
    1500             :   {
    1501       23527 :     x = F2x_add(a, F2xq_sqr(x,T));
    1502       23527 :     if (gc_needed(av, 2))
    1503           0 :       x = gerepileuptoleaf(av, x);
    1504             :   }
    1505        4424 :   return x;
    1506             : }
    1507             : 
    1508             : static void
    1509        6671 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1510             : {
    1511        6671 :   long n = F2x_degree(Tp), r = n/d;
    1512             :   GEN T, f, ff;
    1513       13342 :   if (r==1) { gel(V, idx) = Tp; return; }
    1514        2345 :   T = Tp;
    1515        2345 :   XP = F2x_rem(XP, T);
    1516             :   while (1)
    1517             :   {
    1518        4424 :     pari_sp btop = avma;
    1519             :     long df;
    1520        4424 :     GEN g = random_F2x(n, Tp[1]);
    1521        4424 :     GEN t = F2xq_frobtrace(g, d, T);
    1522        4424 :     if (lgpol(t) == 0) continue;
    1523        3409 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1524        3409 :     if (df > 0 && df < n) break;
    1525        1064 :     avma = btop;
    1526        2079 :   }
    1527        2345 :   ff = F2x_div(Tp, f);
    1528        2345 :   F2x_edf_simple(f, XP, d, V, idx);
    1529        2345 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1530             : }
    1531             : 
    1532             : static GEN
    1533       13034 : F2x_factor_Shoup(GEN T)
    1534             : {
    1535       13034 :   long i, n, s = 0;
    1536             :   GEN XP, D, V;
    1537             :   pari_timer ti;
    1538       13034 :   n = F2x_degree(T);
    1539       13034 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1540       13034 :   XP = F2x_Frobenius(T);
    1541       13034 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1542       13034 :   D = F2x_ddf_simple(T, XP);
    1543       13034 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1544      183393 :   for (i = 1; i <= n; i++)
    1545      170359 :     s += F2x_degree(gel(D,i))/i;
    1546       13034 :   V = cgetg(s+1, t_COL);
    1547      183393 :   for (i = 1, s = 1; i <= n; i++)
    1548             :   {
    1549      170359 :     GEN Di = gel(D,i);
    1550      170359 :     long ni = F2x_degree(Di), ri = ni/i;
    1551      170359 :     if (ni == 0) continue;
    1552       28497 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1553        1981 :     F2x_edf_simple(Di, XP, i, V, s);
    1554        1981 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1555        1981 :     s += ri;
    1556             :   }
    1557       13034 :   return V;
    1558             : }
    1559             : 
    1560             : static GEN
    1561        1659 : F2x_simplefact_Shoup(GEN T)
    1562             : {
    1563        1659 :   long i, n, s = 0, j = 1, k;
    1564             :   GEN XP, D, V;
    1565             :   pari_timer ti;
    1566        1659 :   n = F2x_degree(T);
    1567        1659 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1568        1659 :   XP = F2x_Frobenius(T);
    1569        1659 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1570        1659 :   D = F2x_ddf_simple(T, XP);
    1571        1659 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1572        8470 :   for (i = 1; i <= n; i++)
    1573        6811 :     s += F2x_degree(gel(D,i))/i;
    1574        1659 :   V = cgetg(s+1, t_VECSMALL);
    1575        8470 :   for (i = 1; i <= n; i++)
    1576             :   {
    1577        6811 :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1578        6811 :     if (ni == 0) continue;
    1579        4046 :     for (k = 1; k <= ri; k++)
    1580        2072 :       V[j++] = i;
    1581             :   }
    1582        1659 :   return V;
    1583             : }
    1584             : 
    1585             : static GEN
    1586        9569 : F2x_factor_Cantor(GEN T)
    1587             : {
    1588             :   GEN E, F, M, V;
    1589             :   long i, j, s, l;
    1590        9569 :   V = F2x_factor_squarefree(T);
    1591        9569 :   l = lg(V);
    1592      198184 :   for (s=0, i=1; i < l; i++)
    1593      188615 :     s += !!F2x_degree(gel(V,i));
    1594        9569 :   F = cgetg(s+1, t_VEC);
    1595        9569 :   E = cgetg(s+1, t_VEC);
    1596      198184 :   for (i=1, j=1; i < l; i++)
    1597      188615 :     if (F2x_degree(gel(V,i)))
    1598             :     {
    1599       13034 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1600       13034 :       gel(F, j) = Fj;
    1601       13034 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1602       13034 :       j++;
    1603             :     }
    1604        9569 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1605        9569 :   return sort_factor_pol(M, cmpGuGu);
    1606             : }
    1607             : 
    1608             : static GEN
    1609        1617 : F2x_simplefact_Cantor(GEN T)
    1610             : {
    1611             :   GEN E, F, M, V;
    1612             :   long i, j, s, l;
    1613        1617 :   V = F2x_factor_squarefree(T);
    1614        1617 :   l = lg(V);
    1615       11298 :   for (s=0, i=1; i < l; i++)
    1616        9681 :     s += !!F2x_degree(gel(V,i));
    1617        1617 :   F = cgetg(s+1, t_VEC);
    1618        1617 :   E = cgetg(s+1, t_VEC);
    1619       11298 :   for (i=1, j=1; i < l; i++)
    1620        9681 :     if (F2x_degree(gel(V,i)))
    1621             :     {
    1622        1659 :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1623        1659 :       gel(F, j) = Fj;
    1624        1659 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1625        1659 :       j++;
    1626             :     }
    1627        1617 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1628        1617 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1629             : }
    1630             : 
    1631             : static int
    1632         140 : F2x_isirred_Cantor(GEN T)
    1633             : {
    1634         140 :   pari_sp av = avma;
    1635             :   pari_timer ti;
    1636             :   long n, d;
    1637         140 :   GEN dT = F2x_deriv(T);
    1638             :   GEN XP, D;
    1639         140 :   if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
    1640          98 :   n = F2x_degree(T);
    1641          98 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1642          98 :   XP = F2x_Frobenius(T);
    1643          98 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1644          98 :   D = F2x_ddf_simple(T, XP);
    1645          98 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1646          98 :   d = F2x_degree(gel(D, n));
    1647          98 :   avma = av; return d==n;
    1648             : }
    1649             : 
    1650             : static GEN
    1651       14161 : F2x_factcantor_i(GEN f, long flag)
    1652             : {
    1653       14161 :   long d = F2x_degree(f);
    1654       14161 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1655       11326 :   switch(flag)
    1656             :   {
    1657        9569 :     default: return F2x_factor_Cantor(f);
    1658        1617 :     case 1: return F2x_simplefact_Cantor(f);
    1659         140 :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1660             :   }
    1661             : }
    1662             : 
    1663             : GEN
    1664        7504 : F2x_factcantor(GEN f, long flag)
    1665             : {
    1666        7504 :   pari_sp av = avma;
    1667        7504 :   GEN z = F2x_factcantor_i(f, flag);
    1668        7504 :   if (flag == 2) { avma = av; return z; }
    1669        7504 :   return gerepilecopy(av, z);
    1670             : }
    1671             : 
    1672             : GEN
    1673           0 : F2x_degfact(GEN f)
    1674             : {
    1675           0 :   pari_sp av = avma;
    1676           0 :   GEN z = F2x_factcantor_i(f, 1);
    1677           0 :   return gerepilecopy(av, z);
    1678             : }
    1679             : 
    1680             : int
    1681         140 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f, 2); }
    1682             : 
    1683             : /* Adapted from Shoup NTL */
    1684             : static GEN
    1685      311851 : Flx_factor_squarefree(GEN f, ulong p)
    1686             : {
    1687      311851 :   pari_sp av = avma;
    1688             :   GEN r, t, v, tv;
    1689      311851 :   long q, n = degpol(f);
    1690      311851 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1691      312236 :   for(q = 1;;q *= p)
    1692             :   {
    1693      312236 :     r = Flx_gcd(f, Flx_deriv(f, p), p);
    1694      312236 :     if (degpol(r) == 0)
    1695             :     {
    1696      307364 :       gel(u, q) = f;
    1697      307364 :       break;
    1698             :     }
    1699        4872 :     t = Flx_div(f, r, p);
    1700        4872 :     if (degpol(t) > 0)
    1701             :     {
    1702             :       long j;
    1703       10598 :       for(j = 1;;j++)
    1704             :       {
    1705       10598 :         v = Flx_gcd(r, t, p);
    1706       10598 :         tv = Flx_div(t, v, p);
    1707       10598 :         if (degpol(tv) > 0)
    1708        8414 :           gel(u, j*q) = tv;
    1709       10598 :         if (degpol(v) <= 0) break;
    1710        5929 :         r = Flx_div(r, v, p);
    1711        5929 :         t = v;
    1712        5929 :       }
    1713        4669 :       if (degpol(r) == 0) break;
    1714             :     }
    1715         385 :     f = Flx_deflate(r, p);
    1716         385 :   }
    1717      311851 :   return gerepilecopy(av, u);
    1718             : }
    1719             : 
    1720             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1721             : static GEN
    1722      423778 : Flx_ddf(GEN T, GEN XP, ulong p)
    1723             : {
    1724      423778 :   pari_sp av = avma;
    1725             :   GEN b, g, h, F, f, Tr, xq;
    1726             :   long i, j, n, v, bo, ro;
    1727             :   long B, l, m;
    1728             :   pari_timer ti;
    1729      423778 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1730      423778 :   if (n == 0) return cgetg(1, t_VEC);
    1731      423778 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1732      416729 :   B = n/2;
    1733      416729 :   l = usqrt(B);
    1734      416729 :   m = (B+l-1)/l;
    1735      416729 :   T = Flx_get_red(T, p);
    1736      416729 :   b = cgetg(l+2, t_VEC);
    1737      416729 :   gel(b, 1) = polx_Flx(v);
    1738      416729 :   gel(b, 2) = XP;
    1739      416729 :   bo = brent_kung_optpow(n, l-1, 1);
    1740      416729 :   ro = (bo-1) + (l-1)*((n-1)/bo);
    1741      416729 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1742      416729 :   if (expu(p) <= ro)
    1743      177856 :     for (i = 3; i <= l+1; i++)
    1744      109620 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1745             :   else
    1746             :   {
    1747      348493 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1748      348493 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
    1749      388827 :     for (i = 3; i <= l+1; i++)
    1750       40334 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1751             :   }
    1752      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
    1753      416729 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1754      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
    1755      416729 :   g = cgetg(m+1, t_VEC);
    1756      416729 :   gel(g, 1) = gel(xq, 2);
    1757      923200 :   for(i = 2; i <= m; i++)
    1758      506471 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1759      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
    1760      416729 :   h = cgetg(m+1, t_VEC);
    1761     1339929 :   for (j = 1; j <= m; j++)
    1762             :   {
    1763      923200 :     pari_sp av = avma;
    1764      923200 :     GEN gj = gel(g, j);
    1765      923200 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1766     1507917 :     for (i = 2; i <= l; i++)
    1767      584717 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1768      923200 :     gel(h, j) = gerepileupto(av, e);
    1769             :   }
    1770      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
    1771      416729 :   Tr = get_Flx_mod(T);
    1772      416729 :   F = cgetg(m+1, t_VEC);
    1773     1339929 :   for (j = 1; j <= m; j++)
    1774             :   {
    1775      923200 :     gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
    1776      923200 :     Tr = Flx_div(Tr, gel(F,j), p);
    1777             :   }
    1778      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
    1779      416729 :   f = const_vec(n, pol1_Flx(v));
    1780     1339929 :   for (j = 1; j <= m; j++)
    1781             :   {
    1782      923200 :     GEN e = gel(F, j);
    1783     1028844 :     for (i=l-1; i >= 0; i--)
    1784             :     {
    1785     1028844 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    1786     1028844 :       if (degpol(u))
    1787             :       {
    1788      372030 :         gel(f, l*j-i) = u;
    1789      372030 :         e = Flx_div(e, u, p);
    1790             :       }
    1791     1028844 :       if (!degpol(e)) break;
    1792             :     }
    1793             :   }
    1794      416729 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
    1795      416729 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    1796      416729 :   return gerepilecopy(av, f);
    1797             : }
    1798             : 
    1799             : static void
    1800       45682 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1801             : {
    1802       45682 :   long n = degpol(Tp), r = n/d;
    1803             :   GEN T, f, ff;
    1804             :   ulong p2;
    1805       91364 :   if (r==1) { gel(V, idx) = Tp; return; }
    1806       20776 :   p2 = p>>1;
    1807       20776 :   T = Flx_get_red(Tp, p);
    1808       20776 :   XP = Flx_rem(XP, T, p);
    1809             :   while (1)
    1810             :   {
    1811       22414 :     pari_sp btop = avma;
    1812             :     long i;
    1813       22414 :     GEN g = random_Flx(n, Tp[1], p);
    1814       22414 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1815       22414 :     if (lgpol(t) == 0) continue;
    1816       44429 :     for(i=1; i<=10; i++)
    1817             :     {
    1818       43197 :       pari_sp btop2 = avma;
    1819       43197 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    1820       43197 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    1821       43197 :       if (degpol(f) > 0 && degpol(f) < n) break;
    1822       22421 :       avma = btop2;
    1823             :     }
    1824       22008 :     if (degpol(f) > 0 && degpol(f) < n) break;
    1825        1232 :     avma = btop;
    1826        1638 :   }
    1827       20776 :   f = Flx_normalize(f, p);
    1828       20776 :   ff = Flx_div(Tp, f ,p);
    1829       20776 :   Flx_edf_simple(f, XP, d, p, V, idx);
    1830       20776 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    1831             : }
    1832             : static void
    1833             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    1834             : 
    1835             : static void
    1836       33299 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    1837             : {
    1838             :   pari_sp av;
    1839       33299 :   GEN Tp = get_Flx_mod(T);
    1840       33299 :   long n = degpol(hp), vT = Tp[1];
    1841             :   GEN u1, u2, f1, f2;
    1842       33299 :   ulong p2 = p>>1;
    1843             :   GEN R, h;
    1844       33299 :   h = Flx_get_red(hp, p);
    1845       33299 :   t = Flx_rem(t, T, p);
    1846       33299 :   av = avma;
    1847             :   do
    1848             :   {
    1849       55972 :     avma = av;
    1850       55972 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    1851       55972 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    1852       55972 :   } while (degpol(u1)==0 || degpol(u1)==n);
    1853       33299 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    1854       33299 :   f1 = Flx_normalize(f1, p);
    1855       33299 :   u2 = Flx_div(hp, u1, p);
    1856       33299 :   f2 = Flx_div(Tp, f1, p);
    1857       33299 :   if (degpol(u1)==1)
    1858             :   {
    1859       26194 :     if (degpol(f1)==d)
    1860       25732 :       gel(V, idx) = f1;
    1861             :     else
    1862         462 :       Flx_edf(f1, XP, d, p, V, idx);
    1863             :   }
    1864             :   else
    1865        7105 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    1866       33299 :   idx += degpol(f1)/d;
    1867       33299 :   if (degpol(u2)==1)
    1868             :   {
    1869       25956 :     if (degpol(f2)==d)
    1870       25508 :       gel(V, idx) = f2;
    1871             :     else
    1872         448 :       Flx_edf(f2, XP, d, p, V, idx);
    1873             :   }
    1874             :   else
    1875        7343 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    1876       33299 : }
    1877             : 
    1878             : static void
    1879       18851 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1880             : {
    1881       18851 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    1882             :   GEN T, h, t;
    1883             :   pari_timer ti;
    1884       37702 :   if (r==1) { gel(V, idx) = Tp; return; }
    1885       18851 :   T = Flx_get_red(Tp, p);
    1886       18851 :   XP = Flx_rem(XP, T, p);
    1887       18851 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1888             :   do
    1889             :   {
    1890       20650 :     GEN g = random_Flx(n, vT, p);
    1891       20650 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1892       20650 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    1893       20650 :     h = Flxq_minpoly(t, T, p);
    1894       20650 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    1895       20650 :   } while (degpol(h) <= 1);
    1896       18851 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    1897             : }
    1898             : 
    1899             : static GEN
    1900       31479 : Flx_factor_Shoup(GEN T, ulong p)
    1901             : {
    1902       31479 :   long i, n, s = 0;
    1903             :   GEN XP, D, V;
    1904       31479 :   long e = expu(p);
    1905             :   pari_timer ti;
    1906       31479 :   n = get_Flx_degree(T);
    1907       31479 :   T = Flx_get_red(T, p);
    1908       31479 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1909       31479 :   XP = Flx_Frobenius(T, p);
    1910       31479 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1911       31479 :   D = Flx_ddf(T, XP, p);
    1912       31479 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1913      344855 :   for (i = 1; i <= n; i++)
    1914      313376 :     s += degpol(gel(D,i))/i;
    1915       31479 :   V = cgetg(s+1, t_COL);
    1916      344855 :   for (i = 1, s = 1; i <= n; i++)
    1917             :   {
    1918      313376 :     GEN Di = gel(D,i);
    1919      313376 :     long ni = degpol(Di), ri = ni/i;
    1920      313376 :     if (ni == 0) continue;
    1921       47544 :     Di = Flx_normalize(Di, p);
    1922       47544 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1923       22071 :     if (ri <= e*expu(e))
    1924       17941 :       Flx_edf(Di, XP, i, p, V, s);
    1925             :     else
    1926        4130 :       Flx_edf_simple(Di, XP, i, p, V, s);
    1927       22071 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    1928       22071 :     s += ri;
    1929             :   }
    1930       31479 :   return V;
    1931             : }
    1932             : 
    1933             : static GEN
    1934      284299 : Flx_simplefact_Shoup(GEN T, ulong p)
    1935             : {
    1936      284299 :   long i, n, s = 0, j = 1, k;
    1937             :   GEN XP, D, V;
    1938             :   pari_timer ti;
    1939      284299 :   n = get_Flx_degree(T);
    1940      284299 :   T = Flx_get_red(T, p);
    1941      284299 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1942      284299 :   XP = Flx_Frobenius(T, p);
    1943      284299 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1944      284299 :   D = Flx_ddf(T, XP, p);
    1945      284299 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1946     1802161 :   for (i = 1; i <= n; i++)
    1947     1517862 :     s += degpol(gel(D,i))/i;
    1948      284299 :   V = cgetg(s+1, t_VECSMALL);
    1949     1802161 :   for (i = 1; i <= n; i++)
    1950             :   {
    1951     1517862 :     long ni = degpol(gel(D,i)), ri = ni/i;
    1952     1517862 :     if (ni == 0) continue;
    1953     1054376 :     for (k = 1; k <= ri; k++)
    1954      689494 :       V[j++] = i;
    1955             :   }
    1956      284299 :   return V;
    1957             : }
    1958             : 
    1959             : static GEN
    1960       28441 : Flx_factor_Cantor(GEN T, ulong p)
    1961             : {
    1962             :   GEN E, F, M, V;
    1963             :   long i, j, s, l;
    1964       28441 :   V = Flx_factor_squarefree(T, p);
    1965       28441 :   l = lg(V);
    1966      375991 :   for (s=0, i=1; i < l; i++)
    1967      347550 :     s += !!degpol(gel(V,i));
    1968       28441 :   F = cgetg(s+1, t_VEC);
    1969       28441 :   E = cgetg(s+1, t_VEC);
    1970      375991 :   for (i=1, j=1; i < l; i++)
    1971      347550 :     if (degpol(gel(V,i)))
    1972             :     {
    1973       31479 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    1974       31479 :       gel(F, j) = Fj;
    1975       31479 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1976       31479 :       j++;
    1977             :     }
    1978       28441 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1979       28441 :   return sort_factor_pol(M, cmpGuGu);
    1980             : }
    1981             : 
    1982             : static GEN
    1983      283410 : Flx_simplefact_Cantor(GEN T, ulong p)
    1984             : {
    1985             :   GEN E, F, M, V;
    1986             :   long i, j, s, l;
    1987      283410 :   V = Flx_factor_squarefree(get_Flx_mod(T), p);
    1988      283410 :   l = lg(V);
    1989     2088371 :   for (s=0, i=1; i < l; i++)
    1990     1804961 :     s += !!degpol(gel(V,i));
    1991      283410 :   F = cgetg(s+1, t_VEC);
    1992      283410 :   E = cgetg(s+1, t_VEC);
    1993     2088371 :   for (i=1, j=1; i < l; i++)
    1994     1804961 :     if (degpol(gel(V,i)))
    1995             :     {
    1996      284299 :       GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
    1997      284299 :       gel(F, j) = Fj;
    1998      284299 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1999      284299 :       j++;
    2000             :     }
    2001      283410 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    2002      283410 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    2003             : }
    2004             : 
    2005             : static int
    2006         581 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2007             : {
    2008         581 :   pari_sp av = avma;
    2009             :   pari_timer ti;
    2010             :   long n, d;
    2011         581 :   GEN T = get_Flx_mod(Tp);
    2012         581 :   GEN dT = Flx_deriv(T, p);
    2013             :   GEN XP, D;
    2014         581 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    2015         441 :   n = get_Flx_degree(T);
    2016         441 :   T = Flx_get_red(Tp, p);
    2017         441 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2018         441 :   XP = Flx_Frobenius(T, p);
    2019         441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2020         441 :   D = Flx_ddf(T, XP, p);
    2021         441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2022         441 :   d = degpol(gel(D, n));
    2023         441 :   avma = av; return d==n;
    2024             : }
    2025             : 
    2026             : static GEN
    2027      398392 : Flx_factcantor_i(GEN f, ulong pp, long flag)
    2028             : {
    2029             :   long d;
    2030      398392 :   if (pp==2) { /*We need to handle 2 specially */
    2031        6482 :     GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2032        6482 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2033        6482 :     return F;
    2034             :   }
    2035      391910 :   d = degpol(f);
    2036      391910 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2037      312432 :   switch(flag)
    2038             :   {
    2039       28441 :     default: return Flx_factor_Cantor(f, pp);
    2040      283410 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2041         581 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2042             :   }
    2043             : }
    2044             : 
    2045             : GEN
    2046        7252 : Flx_factcantor(GEN f, ulong p, long flag)
    2047             : {
    2048        7252 :   pari_sp av = avma;
    2049        7252 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
    2050        7252 :   if (flag == 2) { avma = av; return z; }
    2051        7252 :   return gerepilecopy(av, z);
    2052             : }
    2053             : 
    2054             : GEN
    2055      363826 : Flx_degfact(GEN f, ulong p)
    2056             : {
    2057      363826 :   pari_sp av = avma;
    2058      363826 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
    2059      363826 :   return gerepilecopy(av, z);
    2060             : }
    2061             : 
    2062             : /* T must be squarefree mod p*/
    2063             : GEN
    2064       55640 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2065             : {
    2066             :   GEN XP, D;
    2067             :   pari_timer ti;
    2068       55640 :   long i, s, n = get_Flx_degree(T);
    2069       55640 :   GEN V = const_vecsmall(n, 0);
    2070       55640 :   pari_sp av = avma;
    2071       55640 :   T = Flx_get_red(T, p);
    2072       55640 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2073       55640 :   XP = Flx_Frobenius(T, p);
    2074       55640 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2075       55640 :   D = Flx_ddf(T, XP, p);
    2076       55640 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2077      730788 :   for (i = 1, s = 0; i <= n; i++)
    2078             :   {
    2079      675148 :     V[i] = degpol(gel(D,i))/i;
    2080      675148 :     s += V[i];
    2081             :   }
    2082       55640 :   *nb = s;
    2083       55640 :   avma = av; return V;
    2084             : }
    2085             : 
    2086             : long
    2087       51919 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2088             : {
    2089       51919 :   pari_sp av = avma;
    2090       51919 :   GEN ddf = Flx_ddf(T, XP, p);
    2091       51919 :   long l = lg(ddf), i, s=0;
    2092      612556 :   for(i = 1; i < l; i++)
    2093      560637 :     s += degpol(gel(ddf,i))/i;
    2094       51919 :   avma = av; return s;
    2095             : }
    2096             : 
    2097             : /* T must be squarefree mod p*/
    2098             : long
    2099       51919 : Flx_nbfact(GEN T, ulong p)
    2100             : {
    2101       51919 :   pari_sp av = avma;
    2102       51919 :   GEN XP = Flx_Frobenius(T, p);
    2103       51919 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2104       51919 :   avma = av; return n;
    2105             : }
    2106             : 
    2107             : int
    2108         581 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
    2109             : 
    2110             : /* factor f (FpX or Flx) mod pp.
    2111             :  * flag = 1: return the degrees, not the factors
    2112             :  * flag = 2: return NULL if f is not irreducible.
    2113             :  * Not gerepile-safe */
    2114             : static GEN
    2115          98 : factcantor_i(GEN f, GEN pp, long flag)
    2116             : {
    2117          98 :   if (typ(f) == t_VECSMALL)
    2118             :   { /* lgefint(pp) = 3 */
    2119             :     GEN F;
    2120          63 :     ulong p = pp[2];
    2121          63 :     if (p==2) {
    2122          35 :       F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2123          35 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2124             :     } else {
    2125          28 :       F = Flx_factcantor_i(f,p,flag);
    2126          28 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2127             :     }
    2128          63 :     return F;
    2129             :   }
    2130          35 :   return FpX_factcantor_i(f, pp, flag);
    2131             : }
    2132             : GEN
    2133           0 : FpX_factcantor(GEN f, GEN pp, long flag)
    2134             : {
    2135           0 :   pari_sp av = avma;
    2136             :   GEN z;
    2137           0 :   ZX_factmod_init(&f,pp);
    2138           0 :   z = factcantor_i(f,pp,flag);
    2139           0 :   if (flag == 2) { avma = av; return z; }
    2140           0 :   return gerepilecopy(av, z);
    2141             : }
    2142             : 
    2143             : static GEN
    2144         182 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
    2145             : {
    2146         182 :   pari_sp av = avma;
    2147             :   long j, nbfact;
    2148             :   GEN z, y, t, E, u, v;
    2149             : 
    2150         182 :   factmod_init(&f, p);
    2151         182 :   switch(lg(f))
    2152             :   {
    2153          14 :     case 3: avma = av; return trivial_fact();
    2154          14 :     case 2: return gerepileupto(av, zero_fact_intmod(f, p));
    2155             :   }
    2156         154 :   z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
    2157         154 :   y = cgetg(3, t_MAT); nbfact = lg(t);
    2158         154 :   u = cgetg(nbfact,t_COL); gel(y,1) = u;
    2159         154 :   v = cgetg(nbfact,t_COL); gel(y,2) = v;
    2160         154 :   if (flag)
    2161          84 :     for (j=1; j<nbfact; j++)
    2162             :     {
    2163          56 :       gel(u,j) = utoi(uel(t,j));
    2164          56 :       gel(v,j) = utoi(uel(E,j));
    2165             :     }
    2166             :   else
    2167        7896 :     for (j=1; j<nbfact; j++)
    2168             :     {
    2169        7770 :       gel(u,j) = FpX_to_mod(gel(t,j), p);
    2170        7770 :       gel(v,j) = utoi(uel(E,j));
    2171             :     }
    2172         154 :   return gerepileupto(av, y);
    2173             : }
    2174             : GEN
    2175          98 : factcantor0(GEN f, GEN p, long flag)
    2176          98 : { return factmod_aux(f, p, &factcantor_i, flag); }
    2177             : GEN
    2178          84 : factmod(GEN f, GEN p)
    2179          84 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
    2180             : 
    2181             : /* Use this function when you think f is reducible, and that there are lots of
    2182             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2183             : int
    2184          14 : FpX_is_irred(GEN f, GEN p) {
    2185          14 :   ZX_factmod_init(&f,p);
    2186          14 :   return !!factcantor_i(f,p,2);
    2187             : }
    2188             : GEN
    2189           0 : FpX_degfact(GEN f, GEN p) {
    2190           0 :   pari_sp av = avma;
    2191             :   GEN z;
    2192           0 :   ZX_factmod_init(&f,p);
    2193           0 :   z = factcantor_i(f,p,1);
    2194           0 :   return gerepilecopy(av, z);
    2195             : }
    2196             : GEN
    2197          70 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
    2198             : GEN
    2199          28 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
    2200             : 
    2201             : /* set x <-- x + c*y mod p */
    2202             : /* x is not required to be normalized.*/
    2203             : static void
    2204      602467 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2205             : {
    2206             :   long i, lx, ly;
    2207      602467 :   ulong *x=(ulong *)gx;
    2208      602467 :   ulong *y=(ulong *)gy;
    2209     1204934 :   if (!c) return;
    2210      582547 :   lx = lg(gx);
    2211      582547 :   ly = lg(gy);
    2212      582547 :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2213      582547 :   if (SMALL_ULONG(p))
    2214      579928 :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2215             :   else
    2216        2619 :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2217             : }
    2218             : 
    2219             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    2220             : /* assume x1 != 0 */
    2221             : static GEN
    2222      119546 : deg1_Flx(ulong x1, ulong x0, ulong sv)
    2223             : {
    2224      119546 :   return mkvecsmall3(sv, x0, x1);
    2225             : }
    2226             : 
    2227             : long
    2228       58070 : F2x_split_Berlekamp(GEN *t)
    2229             : {
    2230       58070 :   GEN u = *t, a, b, vker;
    2231       58070 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    2232             : 
    2233       58070 :   if (du == 1) return 1;
    2234       30373 :   if (du == 2)
    2235             :   {
    2236        7819 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    2237             :     {
    2238           0 :       t[0] = mkvecsmall2(sv, 2);
    2239           0 :       t[1] = mkvecsmall2(sv, 3);
    2240           0 :       return 2;
    2241             :     }
    2242        7819 :     return 1;
    2243             :   }
    2244             : 
    2245       22554 :   vker = F2x_Berlekamp_ker(u);
    2246       22554 :   lb = lgcols(vker);
    2247       22554 :   d = lg(vker)-1;
    2248       22554 :   ir = 0;
    2249             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2250       51387 :   for (L=1; L<d; )
    2251             :   {
    2252             :     GEN pol;
    2253        6279 :     if (d == 2)
    2254        5992 :       pol = F2v_to_F2x(gel(vker,2), sv);
    2255             :     else
    2256             :     {
    2257         287 :       GEN v = zero_zv(lb);
    2258         287 :       v[1] = du;
    2259         287 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    2260        1113 :       for (i=2; i<=d; i++)
    2261         826 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    2262         287 :       pol = F2v_to_F2x(v, sv);
    2263             :     }
    2264       12943 :     for (i=ir; i<L && L<d; i++)
    2265             :     {
    2266        6664 :       a = t[i]; la = F2x_degree(a);
    2267        6664 :       if (la == 1) { set_irred(i); }
    2268        6643 :       else if (la == 2)
    2269             :       {
    2270           7 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    2271             :         {
    2272           0 :           t[i] = mkvecsmall2(sv, 2);
    2273           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    2274             :         }
    2275           7 :         set_irred(i);
    2276             :       }
    2277             :       else
    2278             :       {
    2279        6636 :         pari_sp av = avma;
    2280             :         long lb;
    2281        6636 :         b = F2x_rem(pol, a);
    2282        6636 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    2283        6209 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    2284        6209 :         if (lb && lb < la)
    2285             :         {
    2286        6209 :           t[L] = F2x_div(a,b);
    2287        6209 :           t[i]= b; L++;
    2288             :         }
    2289           0 :         else avma = av;
    2290             :       }
    2291             :     }
    2292             :   }
    2293       22554 :   return d;
    2294             : }
    2295             : 
    2296             : /* p != 2 */
    2297             : long
    2298      289943 : Flx_split_Berlekamp(GEN *t, ulong p)
    2299             : {
    2300      289943 :   GEN u = *t, a,b,vker;
    2301      289943 :   long d, i, ir, L, la, lb, sv = u[1];
    2302      289943 :   long l = lg(u);
    2303             :   ulong po2;
    2304             : 
    2305      289943 :   if (p == 2)
    2306             :   {
    2307           0 :     *t = Flx_to_F2x(*t);
    2308           0 :     d = F2x_split_Berlekamp(t);
    2309           0 :     for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
    2310           0 :     return d;
    2311             :   }
    2312      289943 :   la = degpol(u);
    2313      289943 :   if (la == 1) return 1;
    2314      219699 :   if (la == 2)
    2315             :   {
    2316       12598 :     ulong r = Flx_quad_root(u,p,1);
    2317       12598 :     if (r != p)
    2318             :     {
    2319        5663 :       t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
    2320        5663 :       t[1] = deg1_Flx(1, p - r, sv);
    2321        5663 :       return 2;
    2322             :     }
    2323        6935 :     return 1;
    2324             :   }
    2325             : 
    2326      207101 :   vker = Flx_Berlekamp_ker(u,p);
    2327      207101 :   vker = Flm_to_FlxV(vker, sv);
    2328      207101 :   d = lg(vker)-1;
    2329      207101 :   po2 = p >> 1; /* (p-1) / 2 */
    2330      207101 :   ir = 0;
    2331             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2332      718417 :   for (L=1; L<d; )
    2333             :   {
    2334      304215 :     GEN pol = zero_zv(l-2);
    2335      304215 :     pol[1] = sv;
    2336      304215 :     pol[2] = random_Fl(p); /*Assume vker[1]=1*/
    2337      906682 :     for (i=2; i<=d; i++)
    2338      602467 :       Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
    2339      304215 :     (void)Flx_renormalize(pol,l-1);
    2340             : 
    2341      732769 :     for (i=ir; i<L && L<d; i++)
    2342             :     {
    2343      428554 :       a = t[i]; la = degpol(a);
    2344      428554 :       if (la == 1) { set_irred(i); }
    2345      386790 :       else if (la == 2)
    2346             :       {
    2347       66269 :         ulong r = Flx_quad_root(a,p,1);
    2348       66269 :         if (r != p)
    2349             :         {
    2350       54110 :           t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
    2351       54110 :           t[L] = deg1_Flx(1, p - r, sv); L++;
    2352             :         }
    2353       66269 :         set_irred(i);
    2354             :       }
    2355             :       else
    2356             :       {
    2357      320521 :         pari_sp av = avma;
    2358      320521 :         b = Flx_rem(pol, a, p);
    2359      320521 :         if (degpol(b) <= 0) { avma=av; continue; }
    2360      301834 :         b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
    2361      301834 :         b = Flx_gcd(a,b, p); lb = degpol(b);
    2362      301834 :         if (lb && lb < la)
    2363             :         {
    2364      173988 :           b = Flx_normalize(b, p);
    2365      173988 :           t[L] = Flx_div(a,b,p);
    2366      173988 :           t[i]= b; L++;
    2367             :         }
    2368      127846 :         else avma = av;
    2369             :       }
    2370             :     }
    2371             :   }
    2372      207101 :   return d;
    2373             : }
    2374             : 
    2375             : long
    2376         692 : FpX_split_Berlekamp(GEN *t, GEN p)
    2377             : {
    2378         692 :   GEN u = *t, a,b,po2,vker;
    2379         692 :   long d, i, ir, L, la, lb, vu = varn(u);
    2380         692 :   if (lgefint(p) == 3)
    2381             :   {
    2382         630 :     ulong up = p[2];
    2383         630 :     if (up == 2)
    2384             :     {
    2385           0 :       *t = ZX_to_F2x(*t);
    2386           0 :       d = F2x_split_Berlekamp(t);
    2387           0 :       for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
    2388             :     }
    2389             :     else
    2390             :     {
    2391         630 :       *t = ZX_to_Flx(*t, up);
    2392         630 :       d = Flx_split_Berlekamp(t, up);
    2393         630 :       for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
    2394             :     }
    2395         630 :     return d;
    2396             :   }
    2397          62 :   la = degpol(u);
    2398          62 :   if (la == 1) return 1;
    2399          57 :   if (la == 2)
    2400             :   {
    2401           2 :     GEN r = FpX_quad_root(u,p,1);
    2402           2 :     if (r)
    2403             :     {
    2404           0 :       t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
    2405           0 :       t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
    2406           0 :       return 2;
    2407             :     }
    2408           2 :     return 1;
    2409             :   }
    2410          55 :   vker = FpX_Berlekamp_ker(u,p);
    2411          55 :   vker = RgM_to_RgXV(vker,vu);
    2412          55 :   d = lg(vker)-1;
    2413          55 :   po2 = shifti(p, -1); /* (p-1) / 2 */
    2414          55 :   ir = 0;
    2415             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2416         162 :   for (L=1; L<d; )
    2417             :   {
    2418          52 :     GEN pol = scalar_ZX_shallow(randomi(p), vu);
    2419         208 :     for (i=2; i<=d; i++)
    2420         156 :       pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
    2421          52 :     pol = FpX_red(pol,p);
    2422         157 :     for (i=ir; i<L && L<d; i++)
    2423             :     {
    2424         105 :       a = t[i]; la = degpol(a);
    2425         105 :       if (la == 1) { set_irred(i); }
    2426          98 :       else if (la == 2)
    2427             :       {
    2428          14 :         GEN r = FpX_quad_root(a,p,1);
    2429          14 :         if (r)
    2430             :         {
    2431          14 :           t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
    2432          14 :           t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
    2433             :         }
    2434          14 :         set_irred(i);
    2435             :       }
    2436             :       else
    2437             :       {
    2438          84 :         pari_sp av = avma;
    2439          84 :         b = FpX_rem(pol, a, p);
    2440          84 :         if (degpol(b) <= 0) { avma=av; continue; }
    2441          63 :         b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
    2442          63 :         b = FpX_gcd(a,b, p); lb = degpol(b);
    2443          63 :         if (lb && lb < la)
    2444             :         {
    2445          43 :           b = FpX_normalize(b, p);
    2446          43 :           t[L] = FpX_div(a,b,p);
    2447          43 :           t[i]= b; L++;
    2448             :         }
    2449          20 :         else avma = av;
    2450             :       }
    2451             :     }
    2452             :   }
    2453          55 :   return d;
    2454             : }
    2455             : 
    2456             : static GEN
    2457      130856 : F2x_Berlekamp_i(GEN f, long flag)
    2458             : {
    2459      130856 :   long e, nbfact, val, d = F2x_degree(f);
    2460             :   ulong k, j;
    2461             :   GEN y, E, f2, g1, t;
    2462             : 
    2463      130856 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    2464             : 
    2465             :   /* to hold factors and exponents */
    2466       59568 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2467       59568 :   E = cgetg(d+1,t_VECSMALL);
    2468       59568 :   val = F2x_valrem(f, &f);
    2469       59568 :   e = nbfact = 1;
    2470       59568 :   if (val) {
    2471       18590 :     if (flag == 2 && val > 1) return NULL;
    2472       18590 :     if (flag == 1)
    2473           0 :       t[1] = 1;
    2474             :     else
    2475       18590 :       gel(t,1) = polx_F2x(f[1]);
    2476       18590 :     E[1] = val; nbfact++;
    2477             :   }
    2478             : 
    2479             :   for(;;)
    2480             :   {
    2481      101372 :     f2 = F2x_gcd(f,F2x_deriv(f));
    2482      101372 :     if (flag == 2 && F2x_degree(f2)) return NULL;
    2483      101372 :     g1 = F2x_degree(f2)? F2x_div(f,f2): f; /* squarefree */
    2484      101372 :     k = 0;
    2485      275444 :     while (F2x_degree(g1)>0)
    2486             :     {
    2487             :       GEN u;
    2488       72700 :       k++; if (k%2 == 0) { k++; f2 = F2x_div(f2,g1); }
    2489       72700 :       u = g1; /* deg(u) > 0 */
    2490       72700 :       if (!F2x_degree(f2)) g1 = pol1_F2x(0); /* only its degree (= 0) matters */
    2491             :       else
    2492             :       {
    2493             :         long dg1;
    2494       17570 :         g1 = F2x_gcd(f2,g1);
    2495       17570 :         dg1 = F2x_degree(g1);
    2496       17570 :         if (dg1)
    2497             :         {
    2498       15757 :           f2= F2x_div(f2,g1);
    2499       15757 :           if (F2x_degree(u) == dg1) continue;
    2500        1127 :           u = F2x_div( u,g1);
    2501             :         }
    2502             :       }
    2503             :       /* u is square-free (product of irred. of multiplicity e * k) */
    2504       58070 :       gel(t,nbfact) = u;
    2505       58070 :       d = F2x_split_Berlekamp(&gel(t,nbfact));
    2506       58070 :       if (flag == 2 && d != 1) return NULL;
    2507       58070 :       if (flag == 1)
    2508           0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = F2x_degree(gel(t,nbfact+j));
    2509       58070 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2510       58070 :       nbfact += d;
    2511             :     }
    2512      101372 :     j = F2x_degree(f2); if (!j) break;
    2513       41804 :     e *= 2; f = F2x_deflate(f2, 2);
    2514       41804 :   }
    2515       59568 :   if (flag == 2) return gen_1; /* irreducible */
    2516       59568 :   setlg(t, nbfact);
    2517       59568 :   setlg(E, nbfact); y = mkvec2(t,E);
    2518       59568 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2519       59568 :               : sort_factor_pol(y, cmpGuGu);
    2520             : }
    2521             : 
    2522             : static GEN
    2523      633669 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
    2524             : {
    2525      633669 :   long e, nbfact, val, d = degpol(f);
    2526             :   ulong k, j;
    2527             :   GEN y, E, f2, g1, t;
    2528             : 
    2529      633669 :   if (p == 2)
    2530             :   {
    2531       10927 :     GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
    2532       10927 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2533       10927 :     return F;
    2534             :   }
    2535      622742 :   if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
    2536             : 
    2537             :   /* to hold factors and exponents */
    2538      283577 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2539      283577 :   E = cgetg(d+1,t_VECSMALL);
    2540      283577 :   val = Flx_valrem(f, &f);
    2541      283577 :   e = nbfact = 1;
    2542      283577 :   if (val) {
    2543       12496 :     if (flag == 2 && val > 1) return NULL;
    2544       12496 :     if (flag == 1)
    2545           0 :       t[1] = 1;
    2546             :     else
    2547       12496 :       gel(t,1) = polx_Flx(f[1]);
    2548       12496 :     E[1] = val; nbfact++;
    2549             :   }
    2550             : 
    2551             :   for(;;)
    2552             :   {
    2553      335013 :     f2 = Flx_gcd(f,Flx_deriv(f,p), p);
    2554      335013 :     if (flag == 2 && degpol(f2)) return NULL;
    2555      335013 :     g1 = degpol(f2)? Flx_div(f,f2,p): f; /* squarefree */
    2556      335013 :     k = 0;
    2557     1002577 :     while (degpol(g1)>0)
    2558             :     {
    2559             :       GEN u;
    2560      332551 :       k++; if (k%p == 0) { k++; f2 = Flx_div(f2,g1,p); }
    2561      332551 :       u = g1; /* deg(u) > 0 */
    2562      332551 :       if (!degpol(f2)) g1 = pol1_Flx(0); /* only its degree (= 0) matters */
    2563             :       else
    2564             :       {
    2565       51332 :         g1 = Flx_gcd(f2,g1, p);
    2566       51332 :         if (degpol(g1))
    2567             :         {
    2568       51185 :           f2= Flx_div(f2,g1,p);
    2569       51185 :           if (lg(u) == lg(g1)) continue;
    2570        7947 :           u = Flx_div( u,g1,p);
    2571             :         }
    2572             :       }
    2573             :       /* u is square-free (product of irred. of multiplicity e * k) */
    2574      289313 :       u = Flx_normalize(u,p);
    2575      289313 :       gel(t,nbfact) = u;
    2576      289313 :       d = Flx_split_Berlekamp(&gel(t,nbfact), p);
    2577      289313 :       if (flag == 2 && d != 1) return NULL;
    2578      289313 :       if (flag == 1)
    2579           0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2580      289313 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2581      289313 :       nbfact += d;
    2582             :     }
    2583      335013 :     if (!p) break;
    2584      335013 :     j = degpol(f2); if (!j) break;
    2585       51436 :     if (j % p) pari_err_PRIME("factmod",utoi(p));
    2586       51436 :     e *= p; f = Flx_deflate(f2, p);
    2587       51436 :   }
    2588      283577 :   if (flag == 2) return gen_1; /* irreducible */
    2589      283577 :   setlg(t, nbfact);
    2590      283577 :   setlg(E, nbfact); y = mkvec2(t,E);
    2591      283577 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2592      283577 :               : sort_factor_pol(y, cmpGuGu);
    2593             : }
    2594             : 
    2595             : /* f an FpX or an Flx */
    2596             : static GEN
    2597      719325 : FpX_Berlekamp_i(GEN f, GEN pp, long flag)
    2598             : {
    2599      719325 :   long e, nbfact, val, d = degpol(f);
    2600             :   ulong k, j;
    2601             :   GEN y, E, f2, g1, t;
    2602             : 
    2603      719325 :   if (typ(f) == t_VECSMALL)
    2604             :   {/* lgefint(pp) == 3 */
    2605      719200 :     ulong p = pp[2];
    2606             :     GEN F;
    2607      719200 :     if (p == 2) {
    2608      119915 :       F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
    2609      119915 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2610             :     } else {
    2611      599285 :       F = Flx_Berlekamp_i(f, p, flag);
    2612      599285 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2613             :     }
    2614      719200 :     return F;
    2615             :   }
    2616             :   /* p is large (and odd) */
    2617         125 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    2618             : 
    2619             :   /* to hold factors and exponents */
    2620          58 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2621          58 :   E = cgetg(d+1,t_VECSMALL);
    2622          58 :   val = ZX_valrem(f, &f);
    2623          58 :   e = nbfact = 1;
    2624          58 :   if (val) {
    2625           3 :     if (flag == 2 && val > 1) return NULL;
    2626           3 :     if (flag == 1)
    2627           0 :       t[1] = 1;
    2628             :     else
    2629           3 :       gel(t,1) = pol_x(varn(f));
    2630           3 :     E[1] = val; nbfact++;
    2631             :   }
    2632             : 
    2633          58 :   f2 = FpX_gcd(f,ZX_deriv(f), pp);
    2634          58 :   if (flag == 2 && degpol(f2)) return NULL;
    2635          58 :   g1 = degpol(f2)? FpX_div(f,f2,pp): f; /* squarefree */
    2636          58 :   k = 0;
    2637         179 :   while (degpol(g1)>0)
    2638             :   {
    2639             :     GEN u;
    2640          63 :     k++;
    2641          63 :     u = g1; /* deg(u) > 0 */
    2642          63 :     if (!degpol(f2)) g1 = pol_1(0); /* only its degree (= 0) matters */
    2643             :     else
    2644             :     {
    2645           6 :       g1 = FpX_gcd(f2,g1, pp);
    2646           6 :       if (degpol(g1))
    2647             :       {
    2648           6 :         f2= FpX_div(f2,g1,pp);
    2649           6 :         if (lg(u) == lg(g1)) continue;
    2650           5 :         u = FpX_div( u,g1,pp);
    2651             :       }
    2652             :     }
    2653             :     /* u is square-free (product of irred. of multiplicity e * k) */
    2654          62 :     u = FpX_normalize(u,pp);
    2655          62 :     gel(t,nbfact) = u;
    2656          62 :     d = FpX_split_Berlekamp(&gel(t,nbfact), pp);
    2657          62 :     if (flag == 2 && d != 1) return NULL;
    2658          62 :     if (flag == 1)
    2659           0 :       for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2660          62 :     for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2661          62 :     nbfact += d;
    2662             :   }
    2663          58 :   if (flag == 2) return gen_1; /* irreducible */
    2664          58 :   setlg(t, nbfact);
    2665          58 :   setlg(E, nbfact); y = mkvec2(t,E);
    2666          58 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2667          58 :               : sort_factor_pol(y, cmpii);
    2668             : }
    2669             : GEN
    2670      719255 : FpX_factor(GEN f, GEN p)
    2671             : {
    2672      719255 :   pari_sp av = avma;
    2673      719255 :   ZX_factmod_init(&f, p);
    2674      719255 :   return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
    2675             : }
    2676             : GEN
    2677       61089 : Flx_factor(GEN f, ulong p)
    2678             : {
    2679       61089 :   pari_sp av = avma;
    2680       61089 :   GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
    2681       61089 :   return gerepilecopy(av, F);
    2682             : }
    2683             : GEN
    2684          14 : F2x_factor(GEN f)
    2685             : {
    2686          14 :   pari_sp av = avma;
    2687          14 :   return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
    2688             : }
    2689             : 
    2690             : GEN
    2691         105 : factormod0(GEN f, GEN p, long flag)
    2692             : {
    2693         105 :   switch(flag)
    2694             :   {
    2695          77 :     case 0: return factmod(f,p);
    2696          28 :     case 1: return simplefactmod(f,p);
    2697           0 :     default: pari_err_FLAG("factormod");
    2698             :   }
    2699           0 :   return NULL; /* not reached */
    2700             : }

Generated by: LCOV version 1.11