Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 20459-9710128) Lines: 1405 1600 87.8 %
Date: 2017-03-30 05:32:39 Functions: 126 137 92.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /***********************************************************************/
      18             : /**                                                                   **/
      19             : /**               Factorisation over finite field                     **/
      20             : /**                                                                   **/
      21             : /***********************************************************************/
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : /* Check types and replace F by a monic normalized FpX having the same roots
      29             :  * Don't bother to make constant polynomials monic */
      30             : static void
      31         228 : factmod_init(GEN *F, GEN p)
      32             : {
      33         228 :   if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
      34         228 :   if (signe(p) < 0) pari_err_PRIME("factmod",p);
      35         228 :   if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
      36         228 :   if (lgefint(p) == 3)
      37             :   {
      38         138 :     ulong pp = p[2];
      39         138 :     if (pp < 2) pari_err_PRIME("factmod", p);
      40         138 :     *F = RgX_to_Flx(*F, pp);
      41         138 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42             :   }
      43             :   else
      44             :   {
      45          90 :     *F = RgX_to_FpX(*F, p);
      46          90 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      47             :   }
      48         228 : }
      49             : /* as above, assume p prime and *F a ZX */
      50             : static void
      51      718859 : ZX_factmod_init(GEN *F, GEN p)
      52             : {
      53      718859 :   if (lgefint(p) == 3)
      54             :   {
      55      713459 :     ulong pp = p[2];
      56      713459 :     *F = ZX_to_Flx(*F, pp);
      57      713459 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      58             :   }
      59             :   else
      60             :   {
      61        5400 :     *F = FpX_red(*F, p);
      62        5400 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      63             :   }
      64      718859 : }
      65             : 
      66             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      67             : static GEN
      68          36 : all_roots_mod_p(ulong p, int not_0)
      69             : {
      70             :   GEN r;
      71             :   ulong i;
      72          36 :   if (not_0) {
      73          24 :     r = cgetg(p, t_VECSMALL);
      74          24 :     for (i = 1; i < p; i++) r[i] = i;
      75             :   } else {
      76          12 :     r = cgetg(p+1, t_VECSMALL);
      77          12 :     for (i = 0; i < p; i++) r[i+1] = i;
      78             :   }
      79          36 :   return r;
      80             : }
      81             : 
      82             : /* X^n - 1 */
      83             : static GEN
      84          12 : Flx_Xnm1(long sv, long n, ulong p)
      85             : {
      86          12 :   GEN t = cgetg(n+3, t_VECSMALL);
      87             :   long i;
      88          12 :   t[1] = sv;
      89          12 :   t[2] = p - 1;
      90          12 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      91          12 :   t[i] = 1; return t;
      92             : }
      93             : /* X^n + 1 */
      94             : static GEN
      95           6 : Flx_Xn1(long sv, long n, ulong p)
      96             : {
      97           6 :   GEN t = cgetg(n+3, t_VECSMALL);
      98             :   long i;
      99             :   (void) p;
     100           6 :   t[1] = sv;
     101           6 :   t[2] = 1;
     102           6 :   for (i = 3; i <= n+1; i++) t[i] = 0;
     103           6 :   t[i] = 1; return t;
     104             : }
     105             : 
     106             : static ulong
     107           6 : Fl_nonsquare(ulong p)
     108             : {
     109           6 :   long k = 2;
     110           6 :   for (;; k++)
     111             :   {
     112          12 :     long i = krouu(k, p);
     113          12 :     if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
     114          18 :     if (i < 0) return k;
     115           6 :   }
     116             : }
     117             : 
     118             : /* f monic Flx, f(0) != 0. Return a monic squarefree g with the same
     119             :  * roots as f */
     120             : static GEN
     121          12 : Flx_cut_out_roots(GEN f, ulong p)
     122             : {
     123          12 :   GEN g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
     124          12 :   if (g != f && degpol(g) >= 0) {
     125           6 :     (void)Flx_valrem(g, &g); /* reduction may introduce 0 root */
     126           6 :     g = Flx_gcd(g, Flx_Xnm1(g[1], p-1, p), p);
     127           6 :     g = Flx_normalize(g, p);
     128             :   }
     129          12 :   return g;
     130             : }
     131             : 
     132             : /* by checking f(0..p-1) */
     133             : GEN
     134          12 : Flx_roots_naive(GEN f, ulong p)
     135             : {
     136          12 :   long d, n = 0;
     137          12 :   ulong s = 1UL, r;
     138          12 :   GEN q, y = cgetg(degpol(f) + 1, t_VECSMALL);
     139          12 :   pari_sp av2, av = avma;
     140             : 
     141          12 :   if (Flx_valrem(f, &f)) y[++n] = 0;
     142          12 :   f = Flx_cut_out_roots(f, p);
     143          12 :   d = degpol(f);
     144          12 :   if (d < 0) return all_roots_mod_p(p, n == 0);
     145          12 :   av2 = avma;
     146         264 :   while (d > 1) /* d = current degree of f */
     147             :   {
     148         246 :     q = Flx_div_by_X_x(f, s, p, &r); /* TODO: FFT-type multi-evaluation */
     149         246 :     if (r) avma = av2; else { y[++n] = s; d--; f = q; av2 = avma; }
     150         246 :     if (++s == p) break;
     151             :   }
     152          12 :   if (d == 1)
     153             :   { /* -f[2]/f[3], root of deg 1 polynomial */
     154           6 :     r = Fl_mul(p - Fl_inv(f[3], p), f[2], p);
     155           6 :     if (r >= s) y[++n] = r; /* otherwise double root */
     156             :   }
     157          12 :   avma = av; fixlg(y, n+1); return y;
     158             : }
     159             : static GEN
     160        5934 : Flx_root_mod_2(GEN f)
     161             : {
     162        5934 :   int z1, z0 = !(f[2] & 1);
     163             :   long i,n;
     164             :   GEN y;
     165             : 
     166        5934 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     167        5934 :   z1 = n & 1;
     168        5934 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     169        5934 :   if (z0) y[i++] = 0;
     170        5934 :   if (z1) y[i  ] = 1;
     171        5934 :   return y;
     172             : }
     173             : static ulong
     174          12 : Flx_oneroot_mod_2(GEN f)
     175             : {
     176             :   long i,n;
     177          12 :   if (!(f[2] & 1)) return 0;
     178          12 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     179          12 :   if (n & 1) return 1;
     180           6 :   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     2006116 : 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       55032 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
     196             : {
     197       55032 :   pari_sp av = avma;
     198             :   GEN y;
     199       55032 :   if (gpwrap)
     200          72 :     factmod_init(&f, pp);
     201             :   else
     202       54960 :     ZX_factmod_init(&f, pp);
     203       55032 :   switch(lg(f))
     204             :   {
     205          12 :     case 2: pari_err_ROOTS0("rootmod");
     206          42 :     case 3: avma = av; return cgetg(1,t_COL);
     207             :   }
     208       54978 :   if (typ(f) == t_VECSMALL)
     209             :   {
     210       52650 :     ulong p = pp[2];
     211       52650 :     if (p == 2)
     212        5934 :       y = Flx_root_mod_2(f);
     213             :     else
     214             :     {
     215       46716 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     216       46716 :       y = Roots(f, p);
     217             :     }
     218       52644 :     y = Flc_to_ZC(y);
     219             :   }
     220             :   else
     221        2328 :     y = FpX_roots_i(f, pp);
     222       54966 :   if (gpwrap) y = FpC_to_mod(y, pp);
     223       54966 :   return gerepileupto(av, y);
     224             : }
     225             : /* assume that f is a ZX an pp a prime */
     226             : GEN
     227       54960 : FpX_roots(GEN f, GEN pp)
     228       54960 : { return rootmod_aux(f, pp, Flx_roots_i, 0); }
     229             : /* no assumptions on f and pp */
     230             : GEN
     231          12 : rootmod2(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_naive, 1); }
     232             : GEN
     233          60 : rootmod(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_i, 1); }
     234             : GEN
     235          72 : rootmod0(GEN f, GEN p, long flag)
     236             : {
     237          72 :   switch(flag)
     238             :   {
     239          60 :     case 0: return rootmod(f,p);
     240          12 :     case 1: return rootmod2(f,p);
     241           0 :     default: pari_err_FLAG("polrootsmod");
     242             :   }
     243             :   return NULL; /* LCOV_EXCL_LINE */
     244             : }
     245             : 
     246             : /* assume x reduced mod p > 2, monic. */
     247             : static int
     248          18 : FpX_quad_factortype(GEN x, GEN p)
     249             : {
     250          18 :   GEN b = gel(x,3), c = gel(x,2);
     251          18 :   GEN D = subii(sqri(b), shifti(c,2));
     252          18 :   return kronecker(D,p);
     253             : }
     254             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     255             : GEN
     256        5394 : FpX_quad_root(GEN x, GEN p, int unknown)
     257             : {
     258        5394 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     259             : 
     260        5394 :   if (absequaliu(p, 2)) {
     261           0 :     if (!signe(b)) return c;
     262           0 :     return signe(c)? NULL: gen_1;
     263             :   }
     264        5394 :   D = subii(sqri(b), shifti(c,2));
     265        5394 :   D = remii(D,p);
     266        5394 :   if (unknown && kronecker(D,p) == -1) return NULL;
     267             : 
     268        5376 :   s = Fp_sqrt(D,p);
     269             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     270        5376 :   if (!s) return NULL;
     271        5370 :   return Fp_halve(Fp_sub(s,b, p), p);
     272             : }
     273             : static GEN
     274        2376 : FpX_otherroot(GEN x, GEN r, GEN p)
     275        2376 : { 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    14087183 : Fl_disc_bc(ulong b, ulong c, ulong p)
     280    14087183 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     281             : /* p > 2 */
     282             : static ulong
     283    13989518 : Flx_quad_root(GEN x, ulong p, int unknown)
     284             : {
     285    13989518 :   ulong s, b = x[3], c = x[2];
     286    13989518 :   ulong D = Fl_disc_bc(b, c, p);
     287    13994993 :   if (unknown && krouu(D,p) == -1) return p;
     288     9648371 :   s = Fl_sqrt(D,p);
     289     9644275 :   if (s==~0UL) return p;
     290     9644251 :   return Fl_halve(Fl_sub(s,b, p), p);
     291             : }
     292             : static ulong
     293     9310396 : Flx_otherroot(GEN x, ulong r, ulong p)
     294     9310396 : { 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     3651384 : split_init(struct split_t *S, long max)
     302             : {
     303     3651384 :   S->todo = vectrunc_init(max);
     304     3651580 :   S->done = vectrunc_init(max);
     305     3650861 : }
     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     3814229 : 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     3814520 : split_moveto_done(struct split_t *S, long i, GEN t)
     323             : {
     324     3814520 :   long n = lg(S->todo)-1;
     325     3814520 :   vectrunc_append(S->done, t);
     326     3814919 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     327     3814919 :   setlg(S->todo, n);
     328             : 
     329     3814978 : }
     330             : /* append t to done */
     331             : static void
     332      271034 : split_add_done(struct split_t *S, GEN t)
     333      271034 : { vectrunc_append(S->done, t); }
     334             : /* split todo[i] into a and b */
     335             : static void
     336      234362 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     337             : {
     338      234362 :   gel(S->todo, i) = a;
     339      234362 :   split_add(S, b);
     340      234365 : }
     341             : /* split todo[i] into a and b, moved to done */
     342             : static void
     343      254909 : split_done(struct split_t *S, long i, GEN a, GEN b)
     344             : {
     345      254909 :   split_moveto_done(S, i, a);
     346      254913 :   split_add_done(S, b);
     347      254913 : }
     348             : 
     349             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     350             : static GEN
     351        2328 : FpX_roots_i(GEN f, GEN p)
     352             : {
     353             :   GEN pol, pol0, a, q;
     354             :   struct split_t S;
     355             : 
     356        2328 :   split_init(&S, lg(f)-1);
     357        2328 :   settyp(S.done, t_COL);
     358        2328 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     359        2328 :   switch(degpol(f))
     360             :   {
     361           6 :     case 0: return ZC_copy(S.done);
     362          12 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     363             :     case 2: {
     364        1320 :       GEN s, r = FpX_quad_root(f, p, 1);
     365        1320 :       if (r) {
     366        1320 :         split_add_done(&S, r);
     367        1320 :         s = FpX_otherroot(f,r, p);
     368             :         /* f not known to be square free yet */
     369        1320 :         if (!equalii(r, s)) split_add_done(&S, s);
     370             :       }
     371        1320 :       return sort(S.done);
     372             :     }
     373             :   }
     374             : 
     375         990 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     376         990 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     377         990 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     378         990 :   a = FpX_gcd(f,a, p);
     379         990 :   if (!degpol(a)) return ZC_copy(S.done);
     380         990 :   split_add(&S, FpX_normalize(a,p));
     381             : 
     382         990 :   q = shifti(p,-1);
     383         990 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     384         990 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     385        2064 :   for (pol0[2] = 1;; pol0[2]++)
     386             :   {
     387        2064 :     long j, l = lg(S.todo);
     388        2064 :     if (l == 1) return sort(S.done);
     389        1080 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     390        2238 :     for (j = 1; j < l; j++)
     391             :     {
     392        1164 :       GEN c = gel(S.todo,j);
     393        1164 :       switch(degpol(c))
     394             :       { /* convert linear and quadratics to roots, try to split the rest */
     395             :         case 1:
     396          54 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     397          54 :           j--; l--; break;
     398             :         case 2: {
     399        1020 :           GEN r = FpX_quad_root(c, p, 0), s;
     400        1020 :           if (!r) pari_err_PRIME("polrootsmod",p);
     401        1014 :           s = FpX_otherroot(c,r, p);
     402        1014 :           split_done(&S, j, r, s);
     403        1014 :           j--; l--; break;
     404             :         }
     405             :         default: {
     406             :           /* b = pol^(p-1)/2 - 1 */
     407          90 :           GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
     408             :           long db;
     409          90 :           b = FpX_gcd(c,b, p); db = degpol(b);
     410          90 :           if (db && db < degpol(c))
     411             :           {
     412          84 :             b = FpX_normalize(b, p);
     413          84 :             c = FpX_div(c,b, p);
     414          84 :             split_todo(&S, j, b, c);
     415             :           }
     416             :         }
     417             :       }
     418             :     }
     419        1074 :   }
     420             : }
     421             : 
     422             : /* Assume f is normalized */
     423             : static ulong
     424      182220 : Flx_cubic_root(GEN ff, ulong p)
     425             : {
     426      182220 :   GEN f = Flx_normalize(ff,p);
     427      182220 :   ulong pi = get_Fl_red(p);
     428      182220 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     429      182220 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     430      182219 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     431      182218 :   ulong B = Fl_addmul_pre(t, Fl_sub(Fl_double(t2, p), b, p), c, p, pi);
     432      182220 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     433      182221 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     434      182221 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     435      182221 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     436      182218 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     437      182220 :   if (s!=~0UL)
     438             :   {
     439      105815 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     440      105815 :     if (p%3==2) /* 1 solutions */
     441       21372 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     442             :     else
     443             :     {
     444       84443 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     445       84443 :       if (vS1==~0UL) return p; /*0 solutions*/
     446             :       /*3 solutions*/
     447             :     }
     448       66771 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     449       66771 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     450             :   }
     451             :   else
     452             :   {
     453       76405 :     pari_sp av = avma;
     454       76405 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     455       76404 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     456             :     ulong Sa;
     457       76406 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     458       74185 :     Sa = vS1[1];
     459       74185 :     if (p%3==1) /*1 solutions*/
     460             :     {
     461       46855 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     462       46855 :       if (Fa!=P)
     463       30817 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     464             :     }
     465       74185 :     avma = av;
     466       74185 :     return Fl_sub(Fl_double(Sa,p),t,p);
     467             :   }
     468             : }
     469             : 
     470             : /* assume p > 2 prime */
     471             : static ulong
     472     1099517 : Flx_oneroot_i(GEN f, ulong p, long fl)
     473             : {
     474             :   GEN pol, a;
     475             :   ulong q;
     476             :   long da;
     477             : 
     478     1099517 :   if (Flx_val(f)) return 0;
     479     1099384 :   switch(degpol(f))
     480             :   {
     481        6243 :     case 1: return Fl_neg(f[2], p);
     482      290715 :     case 2: return Flx_quad_root(f, p, 1);
     483      167499 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     484             :   }
     485             : 
     486      634900 :   if (!fl)
     487             :   {
     488      611774 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     489      611774 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     490      611774 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     491      611774 :     a = Flx_gcd(f,a, p);
     492       23126 :   } else a = f;
     493      634900 :   da = degpol(a);
     494      634899 :   if (!da) return p;
     495      291223 :   a = Flx_normalize(a,p);
     496             : 
     497      291243 :   q = p >> 1;
     498      291243 :   pol = polx_Flx(f[1]);
     499      425968 :   for(pol[2] = 1;; pol[2]++)
     500             :   {
     501      425968 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     502      425948 :     switch(da)
     503             :     {
     504      107050 :       case 1: return Fl_neg(a[2], p);
     505      169539 :       case 2: return Flx_quad_root(a, p, 0);
     506       14719 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     507             :       default: {
     508      134640 :         GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
     509             :         long db;
     510      134738 :         b = Flx_gcd(a,b, p); db = degpol(b);
     511      134730 :         if (db && db < da)
     512             :         {
     513      124106 :           b = Flx_normalize(b, p);
     514      124109 :           if (db <= (da >> 1)) {
     515       80637 :             a = b;
     516       80637 :             da = db;
     517             :           } else {
     518       43472 :             a = Flx_div(a,b, p);
     519       43468 :             da -= db;
     520             :           }
     521             :         }
     522             :       }
     523             :     }
     524      134729 :   }
     525             : }
     526             : 
     527             : /* assume p > 2 prime */
     528             : static GEN
     529        3018 : FpX_oneroot_i(GEN f, GEN p)
     530             : {
     531             :   GEN pol, pol0, a, q;
     532             :   long da;
     533             : 
     534        3018 :   if (ZX_val(f)) return gen_0;
     535        3012 :   switch(degpol(f))
     536             :   {
     537          96 :     case 1: return subii(p, gel(f,2));
     538        2856 :     case 2: return FpX_quad_root(f, p, 1);
     539             :   }
     540             : 
     541          60 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     542          60 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     543          60 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     544          60 :   a = FpX_gcd(f,a, p);
     545          60 :   da = degpol(a);
     546          60 :   if (!da) return NULL;
     547          60 :   a = FpX_normalize(a,p);
     548             : 
     549          60 :   q = shifti(p,-1);
     550          60 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     551          60 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     552         192 :   for (pol0[2]=1; ; pol0[2]++)
     553             :   {
     554         192 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     555         192 :     switch(da)
     556             :     {
     557          36 :       case 1: return subii(p, gel(a,2));
     558          24 :       case 2: return FpX_quad_root(a, p, 0);
     559             :       default: {
     560         132 :         GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, a,p), gen_1, p);
     561             :         long db;
     562         132 :         b = FpX_gcd(a,b, p); db = degpol(b);
     563         132 :         if (db && db < da)
     564             :         {
     565         126 :           b = FpX_normalize(b, p);
     566         126 :           if (db <= (da >> 1)) {
     567          90 :             a = b;
     568          90 :             da = db;
     569             :           } else {
     570          36 :             a = FpX_div(a,b, p);
     571          36 :             da -= db;
     572             :           }
     573             :         }
     574             :       }
     575             :     }
     576         132 :   }
     577             : }
     578             : 
     579             : ulong
     580     1069916 : Flx_oneroot(GEN f, ulong p)
     581             : {
     582     1069916 :   pari_sp av = avma;
     583             :   ulong r;
     584     1069916 :   switch(lg(f))
     585             :   {
     586           0 :     case 2: return 0;
     587           1 :     case 3: avma = av; return p;
     588             :   }
     589     1069915 :   if (p == 2) return Flx_oneroot_mod_2(f);
     590     1069915 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     591     1069915 :   avma = av; return r;
     592             : }
     593             : 
     594             : ulong
     595       23371 : Flx_oneroot_split(GEN f, ulong p)
     596             : {
     597       23371 :   pari_sp av = avma;
     598             :   ulong r;
     599       23371 :   switch(lg(f))
     600             :   {
     601           0 :     case 2: return 0;
     602           0 :     case 3: avma = av; return p;
     603             :   }
     604       23371 :   if (p == 2) return Flx_oneroot_mod_2(f);
     605       23371 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     606       23459 :   avma = av; return r;
     607             : }
     608             : 
     609             : /* assume that p is prime */
     610             : GEN
     611        9252 : FpX_oneroot(GEN f, GEN pp) {
     612        9252 :   pari_sp av = avma;
     613        9252 :   ZX_factmod_init(&f, pp);
     614        9252 :   switch(lg(f))
     615             :   {
     616           0 :     case 2: avma = av; return gen_0;
     617           0 :     case 3: avma = av; return NULL;
     618             :   }
     619        9252 :   if (typ(f) == t_VECSMALL)
     620             :   {
     621        6234 :     ulong r, p = pp[2];
     622        6234 :     if (p == 2)
     623          12 :       r = Flx_oneroot_mod_2(f);
     624             :     else
     625        6222 :       r = Flx_oneroot_i(f, p, 0);
     626        6234 :     avma = av;
     627        6234 :     return (r == p)? NULL: utoi(r);
     628             :   }
     629        3018 :   f = FpX_oneroot_i(f, pp);
     630        3018 :   if (!f) { avma = av; return NULL; }
     631        3018 :   return gerepileuptoint(av, f);
     632             : }
     633             : 
     634             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     635             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     636             : /* returned in n                                                          */
     637             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     638             : /* twice as big and decrement until it divides p-1.                       */
     639             : static GEN
     640           0 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     641             : {
     642           0 :    pari_sp ltop = avma;
     643             :    GEN pm, factn, power, base, zeta;
     644             :    long n;
     645             : 
     646           0 :    pm = subis (p, 1ul);
     647           0 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     648           0 :    factn = Z_factor(stoi(n));
     649           0 :    power = diviuexact (pm, n);
     650           0 :    base = gen_1;
     651             :    do {
     652           0 :       base = addis (base, 1l);
     653           0 :       zeta = Fp_pow (base, power, p);
     654             :    }
     655           0 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     656           0 :    *pt_n = n;
     657           0 :    return gerepileuptoint (ltop, zeta);
     658             : }
     659             : 
     660             : GEN
     661           0 : FpX_oneroot_split(GEN fact, GEN p)
     662             : {
     663           0 :   pari_sp av = avma;
     664             :   long n, deg_f, i, dmin;
     665             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     666           0 :   fact = FpX_normalize(fact, p);
     667           0 :   deg_f = degpol(fact);
     668           0 :   if (deg_f<=2) return FpX_oneroot(fact, p);
     669           0 :   minfactor = fact; /* factor of minimal degree found so far */
     670           0 :   dmin = degpol(minfactor);
     671           0 :   prim = good_root_of_unity(p, deg_f, 1, &n);
     672           0 :   expo = diviuexact(subiu(p, 1), n);
     673           0 :   xplusa = pol_x(varn(fact));
     674           0 :   zeta = gen_1;
     675           0 :   while (dmin != 1)
     676             :   {
     677             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     678             :     /* zeta varies over the roots of unity in F_p                         */
     679           0 :     fact = minfactor; deg_f = dmin;
     680             :     /* update X+a, avoid a=0 */
     681           0 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     682           0 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     683           0 :     for (i = 0; i < n; i++)
     684             :     {
     685           0 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     686           0 :       long dtmp = degpol(tmp);
     687           0 :       if (dtmp > 0 && dtmp < deg_f)
     688             :       {
     689           0 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     690           0 :         if (dtmp < dmin)
     691             :         {
     692           0 :           minfactor = FpX_normalize (tmp, p);
     693           0 :           dmin = dtmp;
     694           0 :           if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
     695             :             /* stop early to avoid too many gcds */
     696             :             break;
     697             :         }
     698             :       }
     699           0 :       zeta = Fp_mul (zeta, prim, p);
     700             :     }
     701             :   }
     702           0 :   return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                                                                 */
     707             : /*                     FACTORISATION MODULO p                      */
     708             : /*                                                                 */
     709             : /*******************************************************************/
     710             : 
     711             : /* Functions giving information on the factorisation. */
     712             : 
     713             : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
     714             : static GEN
     715          36 : FpX_Berlekamp_ker(GEN u, GEN p)
     716             : {
     717          36 :   pari_sp ltop=avma;
     718          36 :   long j,N = degpol(u);
     719          36 :   GEN Q  = FpX_matFrobenius(u, p);
     720        1722 :   for (j=1; j<=N; j++)
     721        1686 :     gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
     722          36 :   return gerepileupto(ltop, FpM_ker(Q,p));
     723             : }
     724             : 
     725             : static GEN
     726       14724 : F2x_Berlekamp_ker(GEN u)
     727             : {
     728       14724 :   pari_sp ltop=avma;
     729       14724 :   long j,N = F2x_degree(u);
     730             :   GEN Q;
     731             :   pari_timer T;
     732       14724 :   timer_start(&T);
     733       14724 :   Q = F2x_matFrobenius(u);
     734       80172 :   for (j=1; j<=N; j++)
     735       65448 :     F2m_flip(Q,j,j);
     736       14724 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     737       14724 :   Q = F2m_ker_sp(Q,0);
     738       14724 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     739       14724 :   return gerepileupto(ltop,Q);
     740             : }
     741             : 
     742             : static GEN
     743      217725 : Flx_Berlekamp_ker(GEN u, ulong l)
     744             : {
     745      217725 :   pari_sp ltop=avma;
     746      217725 :   long j,N = degpol(u);
     747             :   GEN Q;
     748             :   pari_timer T;
     749      217725 :   timer_start(&T);
     750      217725 :   Q  = Flx_matFrobenius(u, l);
     751     1019994 :   for (j=1; j<=N; j++)
     752      802269 :     coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
     753      217725 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     754      217725 :   Q = Flm_ker_sp(Q,l,0);
     755      217725 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     756      217725 :   return gerepileupto(ltop,Q);
     757             : }
     758             : 
     759             : /* product of terms of degree 1 in factorization of f */
     760             : GEN
     761      120174 : FpX_split_part(GEN f, GEN p)
     762             : {
     763      120174 :   long n = degpol(f);
     764      120174 :   GEN z, X = pol_x(varn(f));
     765      120174 :   if (n <= 1) return f;
     766      120138 :   f = FpX_red(f, p);
     767      120138 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     768      120138 :   return FpX_gcd(z,f,p);
     769             : }
     770             : 
     771             : /* Compute the number of roots in Fp without counting multiplicity
     772             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     773             : long
     774       84444 : FpX_nbroots(GEN f, GEN p)
     775             : {
     776       84444 :   pari_sp av = avma;
     777       84444 :   GEN z = FpX_split_part(f, p);
     778       84444 :   avma = av; return degpol(z);
     779             : }
     780             : 
     781             : int
     782           0 : FpX_is_totally_split(GEN f, GEN p)
     783             : {
     784           0 :   long n=degpol(f);
     785           0 :   pari_sp av = avma;
     786           0 :   if (n <= 1) return 1;
     787           0 :   if (abscmpui(n, p) > 0) return 0;
     788           0 :   f = FpX_red(f, p);
     789           0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     790             : }
     791             : 
     792             : long
     793     2047538 : Flx_nbroots(GEN f, ulong p)
     794             : {
     795     2047538 :   long n = degpol(f);
     796     2047538 :   pari_sp av = avma;
     797             :   GEN z;
     798     2047538 :   if (n <= 1) return n;
     799     2047346 :   if (n == 2)
     800             :   {
     801             :     ulong D;
     802        5436 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     803        4692 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     804        4692 :     return 1 + krouu(D,p);
     805             :   }
     806     2041910 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     807     2041910 :   z = Flx_gcd(z, f, p);
     808     2041910 :   avma = av; return degpol(z);
     809             : }
     810             : 
     811             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     812             : static GEN
     813        3516 : FpX_ddf(GEN T, GEN XP, GEN p)
     814             : {
     815        3516 :   pari_sp av = avma;
     816             :   GEN b, g, h, F, f, Tr, xq;
     817             :   long i, j, n, v;
     818             :   long B, l, m;
     819             :   pari_timer ti;
     820        3516 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     821        3516 :   if (n == 0) return cgetg(1, t_VEC);
     822        3516 :   if (n == 1) return mkvec(get_FpX_mod(T));
     823        3510 :   B = n/2;
     824        3510 :   l = usqrt(B);
     825        3510 :   m = (B+l-1)/l;
     826        3510 :   T = FpX_get_red(T, p);
     827        3510 :   b = cgetg(l+2, t_VEC);
     828        3510 :   gel(b, 1) = pol_x(v);
     829        3510 :   gel(b, 2) = XP;
     830        3510 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     831        3510 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     832        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
     833        9138 :   for (i = 3; i <= l+1; i++)
     834        5628 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     835        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
     836        3510 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     837        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
     838        3510 :   g = cgetg(m+1, t_VEC);
     839        3510 :   gel(g, 1) = gel(xq, 2);
     840       12798 :   for(i = 2; i <= m; i++)
     841        9288 :     gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     842        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
     843        3510 :   h = cgetg(m+1, t_VEC);
     844       16308 :   for (j = 1; j <= m; j++)
     845             :   {
     846       12798 :     pari_sp av = avma;
     847       12798 :     GEN gj = gel(g, j);
     848       12798 :     GEN e = FpX_sub(gj, gel(b, 1), p);
     849       36804 :     for (i = 2; i <= l; i++)
     850       24006 :       e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
     851       12798 :     gel(h, j) = gerepileupto(av, e);
     852             :   }
     853        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
     854        3510 :   Tr = get_FpX_mod(T);
     855        3510 :   F = cgetg(m+1, t_VEC);
     856       16308 :   for (j = 1; j <= m; j++)
     857             :   {
     858       12798 :     gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
     859       12798 :     Tr = FpX_div(Tr, gel(F,j), p);
     860             :   }
     861        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
     862        3510 :   f = const_vec(n, pol_1(v));
     863       16308 :   for (j = 1; j <= m; j++)
     864             :   {
     865       12798 :     GEN e = gel(F, j);
     866       14502 :     for (i=l-1; i >= 0; i--)
     867             :     {
     868       14502 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     869       14502 :       if (degpol(u))
     870             :       {
     871        2322 :         gel(f, l*j-i) = u;
     872        2322 :         e = FpX_div(e, u, p);
     873             :       }
     874       14502 :       if (!degpol(e)) break;
     875             :     }
     876             :   }
     877        3510 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
     878        3510 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     879        3510 :   return gerepilecopy(av, f);
     880             : }
     881             : 
     882             : static void
     883           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     884             : {
     885           0 :   long n = degpol(Tp), r = n/d;
     886             :   GEN T, f, ff;
     887             :   GEN p2;
     888           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     889           0 :   p2 = shifti(p,-1);
     890           0 :   T = FpX_get_red(Tp, p);
     891           0 :   XP = FpX_rem(XP, T, p);
     892             :   while (1)
     893             :   {
     894           0 :     pari_sp btop = avma;
     895             :     long i;
     896           0 :     GEN g = random_FpX(n, varn(Tp), p);
     897           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     898           0 :     if (signe(t) == 0) continue;
     899           0 :     for(i=1; i<=10; i++)
     900             :     {
     901           0 :       pari_sp btop2 = avma;
     902           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     903           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     904           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     905           0 :       avma = btop2;
     906             :     }
     907           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     908           0 :     avma = btop;
     909           0 :   }
     910           0 :   f = FpX_normalize(f, p);
     911           0 :   ff = FpX_div(Tp, f ,p);
     912           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     913           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     914             : }
     915             : 
     916             : static void
     917         180 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     918             : {
     919             :   pari_sp av;
     920         180 :   GEN Tp = get_FpX_mod(T);
     921         180 :   long n = degpol(hp), vT = varn(Tp);
     922             :   GEN u1, u2, f1, f2;
     923             :   GEN R, h;
     924         180 :   h = FpX_get_red(hp, p);
     925         180 :   t = FpX_rem(t, T, p);
     926         180 :   av = avma;
     927             :   do
     928             :   {
     929         294 :     avma = av;
     930         294 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     931         294 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     932         294 :   } while (degpol(u1)==0 || degpol(u1)==n);
     933         180 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     934         180 :   f1 = FpX_normalize(f1, p);
     935         180 :   u2 = FpX_div(hp, u1, p);
     936         180 :   f2 = FpX_div(Tp, f1, p);
     937         180 :   if (degpol(u1)==1)
     938          90 :     gel(V, idx) = f1;
     939             :   else
     940          90 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     941         180 :   idx += degpol(f1)/d;
     942         180 :   if (degpol(u2)==1)
     943         102 :     gel(V, idx) = f2;
     944             :   else
     945          78 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     946         180 : }
     947             : 
     948             : static void
     949          12 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     950             : {
     951          12 :   long n = degpol(Tp), r = n/d, vT = varn(Tp);
     952             :   GEN T, h, t;
     953             :   pari_timer ti;
     954          24 :   if (r==1) { gel(V, idx) = Tp; return; }
     955          12 :   T = FpX_get_red(Tp, p);
     956          12 :   XP = FpX_rem(XP, T, p);
     957          12 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     958             :   do
     959             :   {
     960          12 :     GEN g = random_FpX(n, vT, p);
     961          12 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     962          12 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     963          12 :     h = FpXQ_minpoly(t, T, p);
     964          12 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     965          12 :   } while (degpol(h) != r);
     966          12 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     967             : }
     968             : 
     969             : static GEN
     970          12 : FpX_factor_Shoup(GEN T, GEN p)
     971             : {
     972          12 :   long i, n, s = 0;
     973             :   GEN XP, D, V;
     974          12 :   long e = expi(p);
     975             :   pari_timer ti;
     976          12 :   n = get_FpX_degree(T);
     977          12 :   T = FpX_get_red(T, p);
     978          12 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     979          12 :   XP = FpX_Frobenius(T, p);
     980          12 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     981          12 :   D = FpX_ddf(T, XP, p);
     982          12 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     983         258 :   for (i = 1; i <= n; i++)
     984         246 :     s += degpol(gel(D,i))/i;
     985          12 :   V = cgetg(s+1, t_COL);
     986         258 :   for (i = 1, s = 1; i <= n; i++)
     987             :   {
     988         246 :     GEN Di = gel(D,i);
     989         246 :     long ni = degpol(Di), ri = ni/i;
     990         246 :     if (ni == 0) continue;
     991          18 :     Di = FpX_normalize(Di, p);
     992          18 :     if (ni == i) { gel(V, s++) = Di; continue; }
     993          12 :     if (ri <= e*expu(e))
     994          12 :       FpX_edf(Di, XP, i, p, V, s);
     995             :     else
     996           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     997          12 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     998          12 :     s += ri;
     999             :   }
    1000          12 :   return V;
    1001             : }
    1002             : 
    1003             : static GEN
    1004           0 : FpX_simplefact_Shoup(GEN T, GEN p)
    1005             : {
    1006           0 :   long i, n, s = 0, j = 1, k;
    1007             :   GEN XP, D, V;
    1008             :   pari_timer ti;
    1009           0 :   n = get_FpX_degree(T);
    1010           0 :   T = FpX_get_red(T, p);
    1011           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1012           0 :   XP = FpX_Frobenius(T, p);
    1013           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1014           0 :   D = FpX_ddf(T, XP, p);
    1015           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1016           0 :   for (i = 1; i <= n; i++)
    1017           0 :     s += degpol(gel(D,i))/i;
    1018           0 :   V = cgetg(s+1, t_VEC);
    1019           0 :   for (i = 1; i <= n; i++)
    1020             :   {
    1021           0 :     long ni = degpol(gel(D,i)), ri = ni/i;
    1022           0 :     if (ni == 0) continue;
    1023           0 :     for (k = 1; k <= ri; k++)
    1024           0 :       gel(V, j++) = utoi(i);
    1025             :   }
    1026           0 :   return V;
    1027             : }
    1028             : 
    1029             : /* Yun algorithm: Assume p > degpol(T) */
    1030             : static GEN
    1031         222 : FpX_factor_Yun(GEN T, GEN p)
    1032             : {
    1033         222 :   long n = degpol(T);
    1034         222 :   long i = 1;
    1035         222 :   GEN d = FpX_deriv(T, p);
    1036             :   GEN a, b, c;
    1037         222 :   GEN V = cgetg(n+1,t_VEC);
    1038         222 :   a = FpX_gcd(T, d, p);
    1039         222 :   if (degpol(a) == 0) return mkvec(T);
    1040         132 :   b = FpX_div(T, a, p);
    1041             :   do
    1042             :   {
    1043         288 :     c = FpX_div(d, a, p);
    1044         288 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1045         288 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1046         288 :     gel(V, i++) = a;
    1047         288 :     b = FpX_div(b, a, p);
    1048         288 :   } while (degpol(b));
    1049         132 :   setlg(V, i); return V;
    1050             : }
    1051             : GEN
    1052         624 : FpX_factor_squarefree(GEN T, GEN p)
    1053             : {
    1054         624 :   if (abscmpiu(p, degpol(T)) <= 0)
    1055             :   {
    1056         444 :     ulong pp = (ulong)p[2];
    1057         444 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1058         444 :     return FlxV_to_ZXV(u);
    1059             :   }
    1060         180 :   return FpX_factor_Yun(T, p);
    1061             : }
    1062             : 
    1063             : /* F / E  a vector of factors / exponents of virtual length l
    1064             :  * (their real lg may be larger). Set their lg to j and return [F,E] */
    1065             : static GEN
    1066      329721 : FE_setlg(GEN F, GEN E, long l)
    1067             : {
    1068      329721 :   setlg(E,l);
    1069      329721 :   setlg(F,l); return mkvec2(F,E);
    1070             : }
    1071             : /* F / E  a vector of vectors of factors / exponents of virtual length l
    1072             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
    1073             : static GEN
    1074      340254 : FE_concat(GEN F, GEN E, long l)
    1075             : {
    1076      340254 :   setlg(E,l); E = shallowconcat1(E);
    1077      340252 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
    1078             : }
    1079             : 
    1080             : static GEN
    1081           6 : FpX_factor_Cantor(GEN T, GEN p)
    1082             : {
    1083           6 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1084           6 :   long i, j, l = lg(V);
    1085           6 :   F = cgetg(l, t_VEC);
    1086           6 :   E = cgetg(l, t_VEC);
    1087          18 :   for (i=1, j=1; i < l; i++)
    1088          12 :     if (degpol(gel(V,i)))
    1089             :     {
    1090          12 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1091          12 :       gel(F, j) = Fj;
    1092          12 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1093          12 :       j++;
    1094             :     }
    1095           6 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1096             : }
    1097             : 
    1098             : static GEN
    1099           0 : FpX_simplefact_Cantor(GEN T, GEN p)
    1100             : {
    1101             :   GEN E, F, V;
    1102             :   long i, j, l;
    1103           0 :   V = FpX_factor_Yun(get_FpX_mod(T), p);
    1104           0 :   l = lg(V);
    1105           0 :   E = cgetg(l, t_VEC);
    1106           0 :   F = cgetg(l, t_VEC);
    1107           0 :   for (i=1, j=1; i < l; i++)
    1108           0 :     if (degpol(gel(V,i)))
    1109             :     {
    1110           0 :       GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
    1111           0 :       gel(F, j) = Fj;
    1112           0 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1113           0 :       j++;
    1114             :     }
    1115           0 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1116             : }
    1117             : 
    1118             : static int
    1119           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1120             : {
    1121           0 :   pari_sp av = avma;
    1122             :   pari_timer ti;
    1123             :   long n, d;
    1124           0 :   GEN T = get_FpX_mod(Tp);
    1125           0 :   GEN dT = FpX_deriv(T, p);
    1126             :   GEN XP, D;
    1127           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1128           0 :   n = get_FpX_degree(T);
    1129           0 :   T = FpX_get_red(Tp, p);
    1130           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1131           0 :   XP = FpX_Frobenius(T, p);
    1132           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1133           0 :   D = FpX_ddf(T, XP, p);
    1134           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1135           0 :   d = degpol(gel(D, n));
    1136           0 :   avma = av; return d==n;
    1137             : }
    1138             : 
    1139             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1140             : 
    1141             : /*Assume that p is large and odd*/
    1142             : static GEN
    1143          30 : FpX_factcantor_i(GEN f, GEN pp, long flag)
    1144             : {
    1145          30 :   long d = degpol(f);
    1146          30 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1147           6 :   switch(flag)
    1148             :   {
    1149           6 :     default: return FpX_factor_Cantor(f, pp);
    1150           0 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1151           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1152             :   }
    1153             : }
    1154             : 
    1155             : long
    1156        3504 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1157             : {
    1158        3504 :   pari_sp av = avma;
    1159        3504 :   GEN ddf = FpX_ddf(T, XP, p);
    1160        3504 :   long l = lg(ddf), i, s=0;
    1161       74328 :   for(i = 1; i < l; i++)
    1162       70824 :     s += degpol(gel(ddf,i))/i;
    1163        3504 :   avma = av; return s;
    1164             : }
    1165             : 
    1166             : long
    1167          24 : FpX_nbfact(GEN T, GEN p)
    1168             : {
    1169          24 :   pari_sp av = avma;
    1170          24 :   GEN XP = FpX_Frobenius(T, p);
    1171          24 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1172          24 :   avma = av; return n;
    1173             : }
    1174             : 
    1175             : /* p > 2 */
    1176             : static GEN
    1177           6 : FpX_is_irred_2(GEN f, GEN p, long d)
    1178             : {
    1179           6 :   switch(d)
    1180             :   {
    1181             :     case -1:
    1182           0 :     case 0: return NULL;
    1183           0 :     case 1: return gen_1;
    1184             :   }
    1185           6 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1186             : }
    1187             : /* p > 2 */
    1188             : static GEN
    1189          12 : FpX_degfact_2(GEN f, GEN p, long d)
    1190             : {
    1191          12 :   switch(d)
    1192             :   {
    1193           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1194           0 :     case 0: return trivial_fact();
    1195           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1196             :   }
    1197          12 :   switch(FpX_quad_factortype(f, p)) {
    1198           6 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1199           6 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1200           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1201             :   }
    1202             : }
    1203             : 
    1204             : GEN
    1205          60 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1206             : GEN
    1207         696 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1208             : 
    1209             : /* Mod(0,p) * x, where x is f's main variable */
    1210             : static GEN
    1211          12 : Mod0pX(GEN f, GEN p)
    1212          12 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
    1213             : static GEN
    1214          12 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
    1215             : 
    1216             : /* not gerepile safe */
    1217             : static GEN
    1218          60 : FpX_factor_2(GEN f, GEN p, long d)
    1219             : {
    1220             :   GEN r, s, R, S;
    1221             :   long v;
    1222             :   int sgn;
    1223          60 :   switch(d)
    1224             :   {
    1225           0 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1226           0 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1227           0 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1228             :   }
    1229          60 :   r = FpX_quad_root(f, p, 1);
    1230          60 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1231          42 :   v = varn(f);
    1232          42 :   s = FpX_otherroot(f, r, p);
    1233          42 :   if (signe(r)) r = subii(p, r);
    1234          42 :   if (signe(s)) s = subii(p, s);
    1235          42 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1236          42 :   R = deg1pol_shallow(gen_1, r, v);
    1237          42 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1238          36 :   S = deg1pol_shallow(gen_1, s, v);
    1239          36 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1240             : }
    1241             : static GEN
    1242          78 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1243             : {
    1244          78 :   switch(flag) {
    1245           6 :     case 2: return FpX_is_irred_2(f, p, d);
    1246          12 :     case 1: return FpX_degfact_2(f, p, d);
    1247          60 :     default: return FpX_factor_2(f, p, d);
    1248             :   }
    1249             : }
    1250             : 
    1251             : static int
    1252       22890 : F2x_quad_factortype(GEN x)
    1253       22890 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1254             : 
    1255             : static GEN
    1256           6 : F2x_is_irred_2(GEN f, long d)
    1257           6 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1258             : 
    1259             : static GEN
    1260         858 : F2x_degfact_2(GEN f, long d)
    1261             : {
    1262         858 :   if (!d) return trivial_fact();
    1263         858 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1264         684 :   switch(F2x_quad_factortype(f)) {
    1265         126 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1266         138 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1267         420 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1268             :   }
    1269             : }
    1270             : 
    1271             : static GEN
    1272       63570 : F2x_factor_2(GEN f, long d)
    1273             : {
    1274       63570 :   long v = f[1];
    1275       63570 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1276       63570 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1277       57066 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1278       16788 :   switch(F2x_quad_factortype(f))
    1279             :   {
    1280        3372 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1281       11406 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1282        2010 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1283             :   }
    1284             : }
    1285             : static GEN
    1286       64434 : F2x_factor_deg2(GEN f, long d, long flag)
    1287             : {
    1288       64434 :   switch(flag) {
    1289           6 :     case 2: return F2x_is_irred_2(f, d);
    1290         858 :     case 1: return F2x_degfact_2(f, d);
    1291       63570 :     default: return F2x_factor_2(f, d);
    1292             :   }
    1293             : }
    1294             : 
    1295             : static void
    1296          12 : split_squares(struct split_t *S, GEN g, ulong p)
    1297             : {
    1298          12 :   ulong q = p >> 1;
    1299          12 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1300          12 :   if (degpol(a) < 0)
    1301             :   {
    1302             :     ulong i;
    1303           6 :     split_add_done(S, (GEN)1);
    1304           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1305             :   } else {
    1306           6 :     (void)Flx_valrem(a, &a);
    1307           6 :     if (degpol(a))
    1308             :     {
    1309           6 :       a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
    1310           6 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1311             :     }
    1312             :   }
    1313          12 : }
    1314             : static void
    1315          12 : split_nonsquares(struct split_t *S, GEN g, ulong p)
    1316             : {
    1317          12 :   ulong q = p >> 1;
    1318          12 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1319          12 :   if (degpol(a) < 0)
    1320             :   {
    1321           6 :     ulong i, z = Fl_nonsquare(p);
    1322           6 :     split_add_done(S, (GEN)z);
    1323           6 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1324             :   } else {
    1325           6 :     (void)Flx_valrem(a, &a);
    1326           6 :     if (degpol(a))
    1327             :     {
    1328           6 :       a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
    1329           6 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1330             :     }
    1331             :   }
    1332          12 : }
    1333             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1334             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1335             : static int
    1336     3648449 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1337             : {
    1338     3648449 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1339     3648593 :   long d = degpol(g);
    1340     3648640 :   if (d < 0) return 0;
    1341     3648505 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1342     3648869 :   if (!d) return 1;
    1343     3635531 :   if (p <= 1.4 * (ulong)d) {
    1344             :     /* small p; split further using x^((p-1)/2) +/- 1.
    1345             :      * 30% degree drop makes the extra gcd worth it. */
    1346          12 :     split_squares(S, g, p);
    1347          12 :     split_nonsquares(S, g, p);
    1348             :   } else { /* large p; use x^(p-1) - 1 directly */
    1349     3635519 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1350     3635019 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1351     3635019 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1352     3634581 :     g = Flx_gcd(g,a, p);
    1353     3634922 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1354             :   }
    1355     3635546 :   return 1;
    1356             : }
    1357             : 
    1358             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1359             : static GEN
    1360    16665793 : Flx_roots_i(GEN f, ulong p)
    1361             : {
    1362             :   GEN pol, g;
    1363    16665793 :   long v = Flx_valrem(f, &g);
    1364             :   ulong q;
    1365             :   struct split_t S;
    1366             : 
    1367             :   /* optimization: test for small degree first */
    1368    16671547 :   switch(degpol(g))
    1369             :   {
    1370             :     case 1: {
    1371       20970 :       ulong r = p - g[2];
    1372       20970 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1373             :     }
    1374             :     case 2: {
    1375    13001488 :       ulong r = Flx_quad_root(g, p, 1), s;
    1376    13014525 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1377     8876232 :       s = Flx_otherroot(g,r, p);
    1378     8882633 :       if (r < s)
    1379     2215332 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1380     6667301 :       else if (r > s)
    1381     6667163 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1382             :       else
    1383         138 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1384             :     }
    1385             :   }
    1386     3649065 :   q = p >> 1;
    1387     3649065 :   split_init(&S, lg(f)-1);
    1388     3648509 :   settyp(S.done, t_VECSMALL);
    1389     3648509 :   if (v) split_add_done(&S, (GEN)0);
    1390     3648509 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1391          36 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1392     3648545 :   pol = polx_Flx(f[1]);
    1393     7521451 :   for (pol[2]=1; ; pol[2]++)
    1394             :   {
    1395     7521451 :     long j, l = lg(S.todo);
    1396     7521451 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1397     3872663 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1398     7989705 :     for (j = 1; j < l; j++)
    1399             :     {
    1400     4117544 :       GEN c = gel(S.todo,j);
    1401     4117544 :       switch(degpol(c))
    1402             :       {
    1403             :         case 1:
    1404     3559569 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1405     3559711 :           j--; l--; break;
    1406             :         case 2: {
    1407      253907 :           ulong r = Flx_quad_root(c, p, 0), s;
    1408      253904 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1409      253898 :           s = Flx_otherroot(c,r, p);
    1410      253896 :           split_done(&S, j, (GEN)r, (GEN)s);
    1411      253898 :           j--; l--; break;
    1412             :         }
    1413             :         default: {
    1414      304014 :           GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
    1415             :           long db;
    1416      304013 :           b = Flx_gcd(c,b, p); db = degpol(b);
    1417      304017 :           if (db && db < degpol(c))
    1418             :           {
    1419      234280 :             b = Flx_normalize(b, p);
    1420      234282 :             c = Flx_div(c,b, p);
    1421      234277 :             split_todo(&S, j, b, c);
    1422             :           }
    1423             :         }
    1424             :       }
    1425             :     }
    1426     3872161 :   }
    1427             : }
    1428             : 
    1429             : GEN
    1430    16622760 : Flx_roots(GEN f, ulong p)
    1431             : {
    1432    16622760 :   pari_sp av = avma;
    1433    16622760 :   switch(lg(f))
    1434             :   {
    1435           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1436           0 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1437             :   }
    1438    16633188 :   if (p == 2) return Flx_root_mod_2(f);
    1439    16633188 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1440             : }
    1441             : 
    1442             : /* assume x reduced mod p, monic. */
    1443             : static int
    1444      100284 : Flx_quad_factortype(GEN x, ulong p)
    1445             : {
    1446      100284 :   ulong b = x[3], c = x[2];
    1447      100284 :   return krouu(Fl_disc_bc(b, c, p), p);
    1448             : }
    1449             : static GEN
    1450           0 : Flx_is_irred_2(GEN f, ulong p, long d)
    1451             : {
    1452           0 :   if (!d) return NULL;
    1453           0 :   if (d == 1) return gen_1;
    1454           0 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1455             : }
    1456             : static GEN
    1457      114978 : Flx_degfact_2(GEN f, ulong p, long d)
    1458             : {
    1459      114978 :   if (!d) return trivial_fact();
    1460      114978 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1461      100284 :   switch(Flx_quad_factortype(f, p)) {
    1462       47664 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1463       51654 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1464         966 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1465             :   }
    1466             : }
    1467             : /* p > 2 */
    1468             : static GEN
    1469      296000 : Flx_factor_2(GEN f, ulong p, long d)
    1470             : {
    1471             :   ulong r, s;
    1472             :   GEN R,S;
    1473      296000 :   long v = f[1];
    1474      296000 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1475      296000 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1476      284180 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1477      200276 :   r = Flx_quad_root(f, p, 1);
    1478      200276 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1479      119299 :   s = Flx_otherroot(f, r, p);
    1480      119299 :   r = Fl_neg(r, p);
    1481      119299 :   s = Fl_neg(s, p);
    1482      119299 :   if (s < r) lswap(s,r);
    1483      119299 :   R = mkvecsmall3(v,r,1);
    1484      119299 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1485      106429 :   S = mkvecsmall3(v,s,1);
    1486      106429 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1487             : }
    1488             : static GEN
    1489      410978 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1490             : {
    1491      410978 :   switch(flag) {
    1492           0 :     case 2: return Flx_is_irred_2(f, p, d);
    1493      114978 :     case 1: return Flx_degfact_2(f, p, d);
    1494      296000 :     default: return Flx_factor_2(f, p, d);
    1495             :   }
    1496             : }
    1497             : 
    1498             : void
    1499       14346 : F2xV_to_FlxV_inplace(GEN v)
    1500             : {
    1501             :   long i;
    1502       14346 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1503       14346 : }
    1504             : void
    1505      560015 : FlxV_to_ZXV_inplace(GEN v)
    1506             : {
    1507             :   long i;
    1508      560015 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1509      560015 : }
    1510             : void
    1511       96486 : F2xV_to_ZXV_inplace(GEN v)
    1512             : {
    1513             :   long i;
    1514       96486 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1515       96486 : }
    1516             : 
    1517             : /* Adapted from Shoup NTL */
    1518             : GEN
    1519       61664 : F2x_factor_squarefree(GEN f)
    1520             : {
    1521             :   GEN r, t, v, tv;
    1522       61664 :   long i, q, n = F2x_degree(f);
    1523       61664 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1524       98570 :   for(q = 1;;q *= 2)
    1525             :   {
    1526       98570 :     r = F2x_gcd(f, F2x_deriv(f));
    1527       98564 :     if (F2x_degree(r) == 0)
    1528             :     {
    1529       51329 :       gel(u, q) = f;
    1530       51329 :       break;
    1531             :     }
    1532       47236 :     t = F2x_div(f, r);
    1533       47235 :     if (F2x_degree(t) > 0)
    1534             :     {
    1535             :       long j;
    1536       44548 :       for(j = 1;;j++)
    1537             :       {
    1538       44548 :         v = F2x_gcd(r, t);
    1539       44548 :         tv = F2x_div(t, v);
    1540       44548 :         if (F2x_degree(tv) > 0)
    1541       16811 :           gel(u, j*q) = tv;
    1542       44547 :         if (F2x_degree(v) <= 0) break;
    1543       28898 :         r = F2x_div(r, v);
    1544       28898 :         t = v;
    1545       28898 :       }
    1546       15649 :       if (F2x_degree(r) == 0) break;
    1547             :     }
    1548       36903 :     f = F2x_sqrt(r);
    1549       36904 :   }
    1550      422601 :   for (i = n; i; i--)
    1551      418646 :     if (F2x_degree(gel(u,i))) break;
    1552       61678 :   setlg(u,i+1); return u;
    1553             : }
    1554             : 
    1555             : static GEN
    1556       22005 : F2x_ddf_simple(GEN T, GEN XP)
    1557             : {
    1558       22005 :   pari_sp av = avma, av2;
    1559             :   GEN f, z, Tr, X;
    1560       22005 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1561       22011 :   if (n == 0) return cgetg(1, t_VEC);
    1562       22011 :   if (n == 1) return mkvec(T);
    1563       16613 :   z = XP; Tr = T; X = polx_F2x(v);
    1564       16614 :   f = const_vec(n, pol1_F2x(v));
    1565       16614 :   av2 = avma;
    1566      134453 :   for (j = 1; j <= B; j++)
    1567             :   {
    1568      126249 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1569      126256 :     if (F2x_degree(u))
    1570             :     {
    1571       35886 :       gel(f, j) = u;
    1572       35886 :       Tr = F2x_div(Tr, u);
    1573       35879 :       av2 = avma;
    1574       90398 :     } else z = gerepileuptoleaf(av2, z);
    1575      126293 :     if (!F2x_degree(Tr)) break;
    1576      117868 :     z = F2xq_sqr(z, Tr);
    1577             :   }
    1578       16612 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1579       16614 :   return gerepilecopy(av, f);
    1580             : }
    1581             : 
    1582             : static GEN
    1583        7023 : F2xq_frobtrace(GEN a, long d, GEN T)
    1584             : {
    1585        7023 :   pari_sp av = avma;
    1586             :   long i;
    1587        7023 :   GEN x = a;
    1588       46609 :   for(i=1; i<d; i++)
    1589             :   {
    1590       39586 :     x = F2x_add(a, F2xq_sqr(x,T));
    1591       39586 :     if (gc_needed(av, 2))
    1592           0 :       x = gerepileuptoleaf(av, x);
    1593             :   }
    1594        7023 :   return x;
    1595             : }
    1596             : 
    1597             : static void
    1598       10753 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1599             : {
    1600       10753 :   long n = F2x_degree(Tp), r = n/d;
    1601             :   GEN T, f, ff;
    1602       21506 :   if (r==1) { gel(V, idx) = Tp; return; }
    1603        3707 :   T = Tp;
    1604        3707 :   XP = F2x_rem(XP, T);
    1605             :   while (1)
    1606             :   {
    1607        7023 :     pari_sp btop = avma;
    1608             :     long df;
    1609        7023 :     GEN g = random_F2x(n, Tp[1]);
    1610        7023 :     GEN t = F2xq_frobtrace(g, d, T);
    1611        7023 :     if (lgpol(t) == 0) continue;
    1612        5386 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1613        5386 :     if (df > 0 && df < n) break;
    1614        1679 :     avma = btop;
    1615        3316 :   }
    1616        3707 :   ff = F2x_div(Tp, f);
    1617        3707 :   F2x_edf_simple(f, XP, d, V, idx);
    1618        3707 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1619             : }
    1620             : 
    1621             : static GEN
    1622       20350 : F2x_factor_Shoup(GEN T)
    1623             : {
    1624       20350 :   long i, n, s = 0;
    1625             :   GEN XP, D, V;
    1626             :   pari_timer ti;
    1627       20350 :   n = F2x_degree(T);
    1628       20350 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1629       20350 :   XP = F2x_Frobenius(T);
    1630       20343 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1631       20343 :   D = F2x_ddf_simple(T, XP);
    1632       20350 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1633      306001 :   for (i = 1; i <= n; i++)
    1634      285652 :     s += F2x_degree(gel(D,i))/i;
    1635       20349 :   V = cgetg(s+1, t_COL);
    1636      305974 :   for (i = 1, s = 1; i <= n; i++)
    1637             :   {
    1638      285626 :     GEN Di = gel(D,i);
    1639      285626 :     long ni = F2x_degree(Di), ri = ni/i;
    1640      285614 :     if (ni == 0) continue;
    1641       47226 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1642        3339 :     F2x_edf_simple(Di, XP, i, V, s);
    1643        3339 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1644        3339 :     s += ri;
    1645             :   }
    1646       20348 :   return V;
    1647             : }
    1648             : 
    1649             : static GEN
    1650        1518 : F2x_simplefact_Shoup(GEN T)
    1651             : {
    1652        1518 :   long i, n, s = 0, j = 1, k;
    1653             :   GEN XP, D, V;
    1654             :   pari_timer ti;
    1655        1518 :   n = F2x_degree(T);
    1656        1518 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1657        1518 :   XP = F2x_Frobenius(T);
    1658        1518 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1659        1518 :   D = F2x_ddf_simple(T, XP);
    1660        1518 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1661        7698 :   for (i = 1; i <= n; i++)
    1662        6180 :     s += F2x_degree(gel(D,i))/i;
    1663        1518 :   V = cgetg(s+1, t_VECSMALL);
    1664        7698 :   for (i = 1; i <= n; i++)
    1665             :   {
    1666        6180 :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1667        6180 :     if (ni == 0) continue;
    1668        3690 :     for (k = 1; k <= ri; k++)
    1669        1890 :       V[j++] = i;
    1670             :   }
    1671        1518 :   return V;
    1672             : }
    1673             : 
    1674             : static GEN
    1675       14678 : F2x_factor_Cantor(GEN T)
    1676             : {
    1677       14678 :   GEN E, F, V = F2x_factor_squarefree(T);
    1678       14679 :   long i, j, l = lg(V);
    1679       14679 :   E = cgetg(l, t_VEC);
    1680       14680 :   F = cgetg(l, t_VEC);
    1681       38275 :   for (i=1, j=1; i < l; i++)
    1682       23596 :     if (F2x_degree(gel(V,i)))
    1683             :     {
    1684       20350 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1685       20348 :       gel(F, j) = Fj;
    1686       20348 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1687       20349 :       j++;
    1688             :     }
    1689       14679 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1690             : }
    1691             : 
    1692             : static GEN
    1693        1482 : F2x_simplefact_Cantor(GEN T)
    1694             : {
    1695        1482 :   GEN E, F, V = F2x_factor_squarefree(T);
    1696        1482 :   long i, j, l = lg(V);
    1697        1482 :   F = cgetg(l, t_VEC);
    1698        1482 :   E = cgetg(l, t_VEC);
    1699        4008 :   for (i=1, j=1; i < l; i++)
    1700        2526 :     if (F2x_degree(gel(V,i)))
    1701             :     {
    1702        1518 :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1703        1518 :       gel(F, j) = Fj;
    1704        1518 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1705        1518 :       j++;
    1706             :     }
    1707        1482 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1708             : }
    1709             : 
    1710             : static int
    1711         204 : F2x_isirred_Cantor(GEN T)
    1712             : {
    1713         204 :   pari_sp av = avma;
    1714             :   pari_timer ti;
    1715             :   long n, d;
    1716         204 :   GEN dT = F2x_deriv(T);
    1717             :   GEN XP, D;
    1718         204 :   if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
    1719         144 :   n = F2x_degree(T);
    1720         144 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1721         144 :   XP = F2x_Frobenius(T);
    1722         144 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1723         144 :   D = F2x_ddf_simple(T, XP);
    1724         144 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf");
    1725         144 :   d = F2x_degree(gel(D, n));
    1726         144 :   avma = av; return d==n;
    1727             : }
    1728             : 
    1729             : static GEN
    1730       19454 : F2x_factcantor_i(GEN f, long flag)
    1731             : {
    1732       19454 :   long d = F2x_degree(f);
    1733       19454 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1734       16364 :   switch(flag)
    1735             :   {
    1736       14678 :     default: return F2x_factor_Cantor(f);
    1737        1482 :     case 1: return F2x_simplefact_Cantor(f);
    1738         204 :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1739             :   }
    1740             : }
    1741             : 
    1742             : GEN
    1743       12908 : F2x_factcantor(GEN f, long flag)
    1744             : {
    1745       12908 :   pari_sp av = avma;
    1746       12908 :   GEN z = F2x_factcantor_i(f, flag);
    1747       12910 :   if (flag == 2) { avma = av; return z; }
    1748       12910 :   return gerepilecopy(av, z);
    1749             : }
    1750             : 
    1751             : GEN
    1752           0 : F2x_degfact(GEN f)
    1753             : {
    1754           0 :   pari_sp av = avma;
    1755           0 :   GEN z = F2x_factcantor_i(f, 1);
    1756           0 :   return gerepilecopy(av, z);
    1757             : }
    1758             : 
    1759             : int
    1760         204 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f, 2); }
    1761             : 
    1762             : /* Adapted from Shoup NTL */
    1763             : GEN
    1764      608711 : Flx_factor_squarefree(GEN f, ulong p)
    1765             : {
    1766      608711 :   long i, q, n = degpol(f);
    1767      608711 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1768      653512 :   for(q = 1;;q *= p)
    1769             :   {
    1770      653512 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1771      653511 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1772       71035 :     t = Flx_div(f, r, p);
    1773       71035 :     if (degpol(t) > 0)
    1774             :     {
    1775             :       long j;
    1776       79012 :       for(j = 1;;j++)
    1777             :       {
    1778       79012 :         v = Flx_gcd(r, t, p);
    1779       79012 :         tv = Flx_div(t, v, p);
    1780       79012 :         if (degpol(tv) > 0) gel(u, j*q) = tv;
    1781       79012 :         if (degpol(v) <= 0) break;
    1782       52479 :         r = Flx_div(r, v, p);
    1783       52479 :         t = v;
    1784       52479 :       }
    1785       26533 :       if (degpol(r) == 0) break;
    1786             :     }
    1787       44801 :     f = Flx_deflate(r, p);
    1788       44801 :   }
    1789     2863796 :   for (i = n; i; i--)
    1790     2861816 :     if (degpol(gel(u,i))) break;
    1791      608711 :   setlg(u,i+1); return u;
    1792             : }
    1793             : 
    1794             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1795             : static GEN
    1796      429156 : Flx_ddf(GEN T, GEN XP, ulong p)
    1797             : {
    1798      429156 :   pari_sp av = avma;
    1799             :   GEN b, g, h, F, f, Tr, xq;
    1800             :   long i, j, n, v, bo, ro;
    1801             :   long B, l, m;
    1802             :   pari_timer ti;
    1803      429156 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1804      429156 :   if (n == 0) return cgetg(1, t_VEC);
    1805      429156 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1806      422186 :   B = n/2;
    1807      422186 :   l = usqrt(B);
    1808      422186 :   m = (B+l-1)/l;
    1809      422186 :   T = Flx_get_red(T, p);
    1810      422186 :   b = cgetg(l+2, t_VEC);
    1811      422186 :   gel(b, 1) = polx_Flx(v);
    1812      422186 :   gel(b, 2) = XP;
    1813      422186 :   bo = brent_kung_optpow(n, l-1, 1);
    1814      422186 :   ro = (bo-1) + (l-1)*((n-1)/bo);
    1815      422186 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1816      422186 :   if (expu(p) <= ro)
    1817      155859 :     for (i = 3; i <= l+1; i++)
    1818       94835 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1819             :   else
    1820             :   {
    1821      361161 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1822      361161 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
    1823      420429 :     for (i = 3; i <= l+1; i++)
    1824       59268 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1825             :   }
    1826      422185 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
    1827      422185 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1828      422186 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
    1829      422186 :   g = cgetg(m+1, t_VEC);
    1830      422186 :   gel(g, 1) = gel(xq, 2);
    1831      932463 :   for(i = 2; i <= m; i++)
    1832      510277 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1833      422186 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
    1834      422186 :   h = cgetg(m+1, t_VEC);
    1835     1354649 :   for (j = 1; j <= m; j++)
    1836             :   {
    1837      932463 :     pari_sp av = avma;
    1838      932463 :     GEN gj = gel(g, j);
    1839      932463 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1840     1493561 :     for (i = 2; i <= l; i++)
    1841      561101 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1842      932460 :     gel(h, j) = gerepileupto(av, e);
    1843             :   }
    1844      422186 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
    1845      422186 :   Tr = get_Flx_mod(T);
    1846      422186 :   F = cgetg(m+1, t_VEC);
    1847     1354648 :   for (j = 1; j <= m; j++)
    1848             :   {
    1849      932462 :     gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
    1850      932463 :     Tr = Flx_div(Tr, gel(F,j), p);
    1851             :   }
    1852      422186 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
    1853      422186 :   f = const_vec(n, pol1_Flx(v));
    1854     1354648 :   for (j = 1; j <= m; j++)
    1855             :   {
    1856      932462 :     GEN e = gel(F, j);
    1857     1046426 :     for (i=l-1; i >= 0; i--)
    1858             :     {
    1859     1046426 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    1860     1046425 :       if (degpol(u))
    1861             :       {
    1862      380915 :         gel(f, l*j-i) = u;
    1863      380915 :         e = Flx_div(e, u, p);
    1864             :       }
    1865     1046426 :       if (!degpol(e)) break;
    1866             :     }
    1867             :   }
    1868      422186 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
    1869      422186 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    1870      422186 :   return gerepilecopy(av, f);
    1871             : }
    1872             : 
    1873             : static void
    1874       36980 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1875             : {
    1876       36980 :   long n = degpol(Tp), r = n/d;
    1877             :   GEN T, f, ff;
    1878             :   ulong p2;
    1879       73960 :   if (r==1) { gel(V, idx) = Tp; return; }
    1880       16711 :   p2 = p>>1;
    1881       16711 :   T = Flx_get_red(Tp, p);
    1882       16711 :   XP = Flx_rem(XP, T, p);
    1883             :   while (1)
    1884             :   {
    1885       18312 :     pari_sp btop = avma;
    1886             :     long i;
    1887       18312 :     GEN g = random_Flx(n, Tp[1], p);
    1888       18312 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1889       18312 :     if (lgpol(t) == 0) continue;
    1890       38502 :     for(i=1; i<=10; i++)
    1891             :     {
    1892       37283 :       pari_sp btop2 = avma;
    1893       37283 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    1894       37283 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    1895       37283 :       if (degpol(f) > 0 && degpol(f) < n) break;
    1896       20572 :       avma = btop2;
    1897             :     }
    1898       17930 :     if (degpol(f) > 0 && degpol(f) < n) break;
    1899        1219 :     avma = btop;
    1900        1601 :   }
    1901       16711 :   f = Flx_normalize(f, p);
    1902       16711 :   ff = Flx_div(Tp, f ,p);
    1903       16711 :   Flx_edf_simple(f, XP, d, p, V, idx);
    1904       16711 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    1905             : }
    1906             : static void
    1907             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    1908             : 
    1909             : static void
    1910       27716 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    1911             : {
    1912             :   pari_sp av;
    1913       27716 :   GEN Tp = get_Flx_mod(T);
    1914       27716 :   long n = degpol(hp), vT = Tp[1];
    1915             :   GEN u1, u2, f1, f2;
    1916       27716 :   ulong p2 = p>>1;
    1917             :   GEN R, h;
    1918       27716 :   h = Flx_get_red(hp, p);
    1919       27716 :   t = Flx_rem(t, T, p);
    1920       27716 :   av = avma;
    1921             :   do
    1922             :   {
    1923       46490 :     avma = av;
    1924       46490 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    1925       46490 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    1926       46490 :   } while (degpol(u1)==0 || degpol(u1)==n);
    1927       27716 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    1928       27716 :   f1 = Flx_normalize(f1, p);
    1929       27716 :   u2 = Flx_div(hp, u1, p);
    1930       27716 :   f2 = Flx_div(Tp, f1, p);
    1931       27716 :   if (degpol(u1)==1)
    1932             :   {
    1933       21776 :     if (degpol(f1)==d)
    1934       21422 :       gel(V, idx) = f1;
    1935             :     else
    1936         354 :       Flx_edf(f1, XP, d, p, V, idx);
    1937             :   }
    1938             :   else
    1939        5940 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    1940       27716 :   idx += degpol(f1)/d;
    1941       27716 :   if (degpol(u2)==1)
    1942             :   {
    1943       21776 :     if (degpol(f2)==d)
    1944       21428 :       gel(V, idx) = f2;
    1945             :     else
    1946         348 :       Flx_edf(f2, XP, d, p, V, idx);
    1947             :   }
    1948             :   else
    1949        5940 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    1950       27716 : }
    1951             : 
    1952             : static void
    1953       15836 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1954             : {
    1955       15836 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    1956             :   GEN T, h, t;
    1957             :   pari_timer ti;
    1958       31672 :   if (r==1) { gel(V, idx) = Tp; return; }
    1959       15836 :   T = Flx_get_red(Tp, p);
    1960       15836 :   XP = Flx_rem(XP, T, p);
    1961       15836 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1962             :   do
    1963             :   {
    1964       17414 :     GEN g = random_Flx(n, vT, p);
    1965       17414 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1966       17414 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    1967       17414 :     h = Flxq_minpoly(t, T, p);
    1968       17414 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    1969       17414 :   } while (degpol(h) <= 1);
    1970       15836 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    1971             : }
    1972             : 
    1973             : static GEN
    1974       26430 : Flx_factor_Shoup(GEN T, ulong p)
    1975             : {
    1976       26430 :   long i, n, s = 0;
    1977             :   GEN XP, D, V;
    1978       26430 :   long e = expu(p);
    1979             :   pari_timer ti;
    1980       26430 :   n = get_Flx_degree(T);
    1981       26430 :   T = Flx_get_red(T, p);
    1982       26430 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1983       26430 :   XP = Flx_Frobenius(T, p);
    1984       26430 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1985       26430 :   D = Flx_ddf(T, XP, p);
    1986       26430 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1987      278126 :   for (i = 1; i <= n; i++)
    1988      251696 :     s += degpol(gel(D,i))/i;
    1989       26430 :   V = cgetg(s+1, t_COL);
    1990      278125 :   for (i = 1, s = 1; i <= n; i++)
    1991             :   {
    1992      251695 :     GEN Di = gel(D,i);
    1993      251695 :     long ni = degpol(Di), ri = ni/i;
    1994      251695 :     if (ni == 0) continue;
    1995       39978 :     Di = Flx_normalize(Di, p);
    1996       39978 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1997       18692 :     if (ri <= e*expu(e))
    1998       15134 :       Flx_edf(Di, XP, i, p, V, s);
    1999             :     else
    2000        3558 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2001       18692 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2002       18692 :     s += ri;
    2003             :   }
    2004       26430 :   return V;
    2005             : }
    2006             : 
    2007             : static GEN
    2008      301158 : Flx_simplefact_Shoup(GEN T, ulong p)
    2009             : {
    2010      301158 :   long i, n, s = 0, j = 1, k;
    2011             :   GEN XP, D, V;
    2012             :   pari_timer ti;
    2013      301158 :   n = get_Flx_degree(T);
    2014      301158 :   T = Flx_get_red(T, p);
    2015      301158 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2016      301158 :   XP = Flx_Frobenius(T, p);
    2017      301158 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2018      301158 :   D = Flx_ddf(T, XP, p);
    2019      301158 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2020     1989852 :   for (i = 1; i <= n; i++)
    2021     1688694 :     s += degpol(gel(D,i))/i;
    2022      301158 :   V = cgetg(s+1, t_VECSMALL);
    2023     1989852 :   for (i = 1; i <= n; i++)
    2024             :   {
    2025     1688694 :     long ni = degpol(gel(D,i)), ri = ni/i;
    2026     1688694 :     if (ni == 0) continue;
    2027     1153806 :     for (k = 1; k <= ri; k++)
    2028      771270 :       V[j++] = i;
    2029             :   }
    2030      301158 :   return V;
    2031             : }
    2032             : 
    2033             : static GEN
    2034       23768 : Flx_factor_Cantor(GEN T, ulong p)
    2035             : {
    2036       23768 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2037       23768 :   long i, j, l = lg(V);
    2038       23768 :   F = cgetg(l, t_VEC);
    2039       23768 :   E = cgetg(l, t_VEC);
    2040       51411 :   for (i=1, j=1; i < l; i++)
    2041       27643 :     if (degpol(gel(V,i)))
    2042             :     {
    2043       26430 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2044       26430 :       gel(F, j) = Fj;
    2045       26430 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2046       26430 :       j++;
    2047             :     }
    2048       23768 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2049             : }
    2050             : 
    2051             : static GEN
    2052      300318 : Flx_simplefact_Cantor(GEN T, ulong p)
    2053             : {
    2054      300318 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2055      300318 :   long i, j, l = lg(V);
    2056      300318 :   F = cgetg(l, t_VEC);
    2057      300318 :   E = cgetg(l, t_VEC);
    2058      603192 :   for (i=1, j=1; i < l; i++)
    2059      302874 :     if (degpol(gel(V,i)))
    2060             :     {
    2061      301158 :       GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
    2062      301158 :       gel(F, j) = Fj;
    2063      301158 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2064      301158 :       j++;
    2065             :     }
    2066      300318 :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    2067             : }
    2068             : 
    2069             : static int
    2070         498 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2071             : {
    2072         498 :   pari_sp av = avma;
    2073             :   pari_timer ti;
    2074             :   long n, d;
    2075         498 :   GEN T = get_Flx_mod(Tp);
    2076         498 :   GEN dT = Flx_deriv(T, p);
    2077             :   GEN XP, D;
    2078         498 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    2079         378 :   n = get_Flx_degree(T);
    2080         378 :   T = Flx_get_red(Tp, p);
    2081         378 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2082         378 :   XP = Flx_Frobenius(T, p);
    2083         378 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2084         378 :   D = Flx_ddf(T, XP, p);
    2085         378 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2086         378 :   d = degpol(gel(D, n));
    2087         378 :   avma = av; return d==n;
    2088             : }
    2089             : 
    2090             : static GEN
    2091      447320 : Flx_factcantor_i(GEN f, ulong pp, long flag)
    2092             : {
    2093             :   long d;
    2094      447320 :   if (pp==2) { /*We need to handle 2 specially */
    2095        6312 :     GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2096        6312 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2097        6312 :     return F;
    2098             :   }
    2099      441008 :   d = degpol(f);
    2100      441008 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2101      324584 :   switch(flag)
    2102             :   {
    2103       23768 :     default: return Flx_factor_Cantor(f, pp);
    2104      300318 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2105         498 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2106             :   }
    2107             : }
    2108             : 
    2109             : GEN
    2110        6242 : Flx_factcantor(GEN f, ulong p, long flag)
    2111             : {
    2112        6242 :   pari_sp av = avma;
    2113        6242 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
    2114        6242 :   if (flag == 2) { avma = av; return z; }
    2115        6242 :   return gerepilecopy(av, z);
    2116             : }
    2117             : 
    2118             : GEN
    2119      417624 : Flx_degfact(GEN f, ulong p)
    2120             : {
    2121      417624 :   pari_sp av = avma;
    2122      417624 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
    2123      417624 :   return gerepilecopy(av, z);
    2124             : }
    2125             : 
    2126             : /* T must be squarefree mod p*/
    2127             : GEN
    2128       52716 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2129             : {
    2130             :   GEN XP, D;
    2131             :   pari_timer ti;
    2132       52716 :   long i, s, n = get_Flx_degree(T);
    2133       52716 :   GEN V = const_vecsmall(n, 0);
    2134       52716 :   pari_sp av = avma;
    2135       52716 :   T = Flx_get_red(T, p);
    2136       52716 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2137       52716 :   XP = Flx_Frobenius(T, p);
    2138       52716 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2139       52716 :   D = Flx_ddf(T, XP, p);
    2140       52716 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    2141      692058 :   for (i = 1, s = 0; i <= n; i++)
    2142             :   {
    2143      639342 :     V[i] = degpol(gel(D,i))/i;
    2144      639342 :     s += V[i];
    2145             :   }
    2146       52716 :   *nb = s;
    2147       52716 :   avma = av; return V;
    2148             : }
    2149             : 
    2150             : long
    2151       48474 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2152             : {
    2153       48474 :   pari_sp av = avma;
    2154       48474 :   GEN ddf = Flx_ddf(T, XP, p);
    2155       48474 :   long l = lg(ddf), i, s=0;
    2156      523338 :   for(i = 1; i < l; i++)
    2157      474864 :     s += degpol(gel(ddf,i))/i;
    2158       48474 :   avma = av; return s;
    2159             : }
    2160             : 
    2161             : /* T must be squarefree mod p*/
    2162             : long
    2163       48474 : Flx_nbfact(GEN T, ulong p)
    2164             : {
    2165       48474 :   pari_sp av = avma;
    2166       48474 :   GEN XP = Flx_Frobenius(T, p);
    2167       48474 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2168       48474 :   avma = av; return n;
    2169             : }
    2170             : 
    2171             : int
    2172         498 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
    2173             : 
    2174             : /* factor f (FpX or Flx) mod pp.
    2175             :  * flag = 1: return the degrees, not the factors
    2176             :  * flag = 2: return NULL if f is not irreducible.
    2177             :  * Not gerepile-safe */
    2178             : static GEN
    2179          84 : factcantor_i(GEN f, GEN pp, long flag)
    2180             : {
    2181          84 :   if (typ(f) == t_VECSMALL)
    2182             :   { /* lgefint(pp) = 3 */
    2183             :     GEN F;
    2184          54 :     ulong p = pp[2];
    2185          54 :     if (p==2) {
    2186          30 :       F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    2187          30 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2188             :     } else {
    2189          24 :       F = Flx_factcantor_i(f,p,flag);
    2190          24 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2191             :     }
    2192          54 :     return F;
    2193             :   }
    2194          30 :   return FpX_factcantor_i(f, pp, flag);
    2195             : }
    2196             : GEN
    2197           0 : FpX_factcantor(GEN f, GEN pp, long flag)
    2198             : {
    2199           0 :   pari_sp av = avma;
    2200             :   GEN z;
    2201           0 :   ZX_factmod_init(&f,pp);
    2202           0 :   z = factcantor_i(f,pp,flag);
    2203           0 :   if (flag == 2) { avma = av; return z; }
    2204           0 :   return gerepilecopy(av, z);
    2205             : }
    2206             : 
    2207             : static GEN
    2208         156 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
    2209             : {
    2210         156 :   pari_sp av = avma;
    2211             :   long j, lfact;
    2212             :   GEN z, t, E, y, u, v;
    2213             : 
    2214         156 :   factmod_init(&f, p);
    2215         156 :   switch(lg(f))
    2216             :   {
    2217          12 :     case 3: avma = av; return trivial_fact();
    2218          12 :     case 2: return gerepileupto(av, zero_fact_intmod(f, p));
    2219             :   }
    2220         132 :   z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
    2221         132 :   lfact = lg(t); y = cgetg(3, t_MAT);
    2222         132 :   gel(y,1) = u = cgetg(lfact,t_COL);
    2223         132 :   gel(y,2) = v = cgetg(lfact,t_COL);
    2224         132 :   if (flag)
    2225          72 :     for (j=1; j<lfact; j++)
    2226             :     {
    2227          48 :       gel(u,j) = utoi(uel(t,j));
    2228          48 :       gel(v,j) = utoi(uel(E,j));
    2229             :     }
    2230             :   else
    2231        6768 :     for (j=1; j<lfact; j++)
    2232             :     {
    2233        6660 :       gel(u,j) = FpX_to_mod(gel(t,j), p);
    2234        6660 :       gel(v,j) = utoi(uel(E,j));
    2235             :     }
    2236         132 :   return gerepileupto(av, y);
    2237             : }
    2238             : GEN
    2239          84 : factcantor0(GEN f, GEN p, long flag)
    2240          84 : { return factmod_aux(f, p, &factcantor_i, flag); }
    2241             : GEN
    2242          72 : factmod(GEN f, GEN p)
    2243          72 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
    2244             : 
    2245             : /* Use this function when you think f is reducible, and that there are lots of
    2246             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2247             : int
    2248          12 : FpX_is_irred(GEN f, GEN p) {
    2249          12 :   ZX_factmod_init(&f,p);
    2250          12 :   return !!factcantor_i(f,p,2);
    2251             : }
    2252             : GEN
    2253           0 : FpX_degfact(GEN f, GEN p) {
    2254           0 :   pari_sp av = avma;
    2255             :   GEN z;
    2256           0 :   ZX_factmod_init(&f,p);
    2257           0 :   z = factcantor_i(f,p,1);
    2258           0 :   return gerepilecopy(av, z);
    2259             : }
    2260             : GEN
    2261          60 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
    2262             : GEN
    2263          24 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
    2264             : 
    2265             : /* set x <-- x + c*y mod p */
    2266             : /* x is not required to be normalized.*/
    2267             : static void
    2268      641220 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2269             : {
    2270             :   long i, lx, ly;
    2271      641220 :   ulong *x=(ulong *)gx;
    2272      641220 :   ulong *y=(ulong *)gy;
    2273     1282440 :   if (!c) return;
    2274      621846 :   lx = lg(gx);
    2275      621846 :   ly = lg(gy);
    2276      621846 :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2277      621846 :   if (SMALL_ULONG(p))
    2278      621336 :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2279             :   else
    2280         510 :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2281             : }
    2282             : 
    2283             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    2284             : /* assume x1 != 0 */
    2285             : static GEN
    2286      123228 : deg1_Flx(ulong x1, ulong x0, ulong sv)
    2287             : {
    2288      123228 :   return mkvecsmall3(sv, x0, x1);
    2289             : }
    2290             : 
    2291             : static long
    2292       42336 : F2x_split_Berlekamp(GEN *t)
    2293             : {
    2294       42336 :   GEN u = *t, a, b, vker;
    2295       42336 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    2296             : 
    2297       42336 :   if (du == 1) return 1;
    2298       20118 :   if (du == 2)
    2299             :   {
    2300        5394 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    2301             :     {
    2302           0 :       t[0] = mkvecsmall2(sv, 2);
    2303           0 :       t[1] = mkvecsmall2(sv, 3);
    2304           0 :       return 2;
    2305             :     }
    2306        5394 :     return 1;
    2307             :   }
    2308             : 
    2309       14724 :   vker = F2x_Berlekamp_ker(u);
    2310       14724 :   lb = lgcols(vker);
    2311       14724 :   d = lg(vker)-1;
    2312       14724 :   ir = 0;
    2313             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2314       34488 :   for (L=1; L<d; )
    2315             :   {
    2316             :     GEN pol;
    2317        5040 :     if (d == 2)
    2318        4782 :       pol = F2v_to_F2x(gel(vker,2), sv);
    2319             :     else
    2320             :     {
    2321         258 :       GEN v = zero_zv(lb);
    2322         258 :       v[1] = du;
    2323         258 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    2324         990 :       for (i=2; i<=d; i++)
    2325         732 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    2326         258 :       pol = F2v_to_F2x(v, sv);
    2327             :     }
    2328       10410 :     for (i=ir; i<L && L<d; i++)
    2329             :     {
    2330        5370 :       a = t[i]; la = F2x_degree(a);
    2331        5370 :       if (la == 1) { set_irred(i); }
    2332        5358 :       else if (la == 2)
    2333             :       {
    2334          18 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    2335             :         {
    2336           0 :           t[i] = mkvecsmall2(sv, 2);
    2337           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    2338             :         }
    2339          18 :         set_irred(i);
    2340             :       }
    2341             :       else
    2342             :       {
    2343        5340 :         pari_sp av = avma;
    2344             :         long lb;
    2345        5340 :         b = F2x_rem(pol, a);
    2346        5340 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    2347        4968 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    2348        4968 :         if (lb && lb < la)
    2349             :         {
    2350        4968 :           t[L] = F2x_div(a,b);
    2351        4968 :           t[i]= b; L++;
    2352             :         }
    2353           0 :         else avma = av;
    2354             :       }
    2355             :     }
    2356             :   }
    2357       14724 :   return d;
    2358             : }
    2359             : 
    2360             : /* p != 2 */
    2361             : static long
    2362      289683 : Flx_split_Berlekamp(GEN *t, ulong p)
    2363             : {
    2364      289683 :   GEN u = *t, a,b,vker;
    2365      289683 :   long d, i, ir, L, la, lb, sv = u[1];
    2366      289683 :   long l = lg(u);
    2367             :   ulong po2;
    2368             : 
    2369      289683 :   if (p == 2)
    2370             :   {
    2371           0 :     *t = Flx_to_F2x(*t);
    2372           0 :     d = F2x_split_Berlekamp(t);
    2373           0 :     for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
    2374           0 :     return d;
    2375             :   }
    2376      289683 :   la = degpol(u);
    2377      289683 :   if (la == 1) return 1;
    2378      228999 :   if (la == 2)
    2379             :   {
    2380       11274 :     ulong r = Flx_quad_root(u,p,1);
    2381       11274 :     if (r != p)
    2382             :     {
    2383        5406 :       t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
    2384        5406 :       t[1] = deg1_Flx(1, p - r, sv);
    2385        5406 :       return 2;
    2386             :     }
    2387        5868 :     return 1;
    2388             :   }
    2389             : 
    2390      217725 :   vker = Flx_Berlekamp_ker(u,p);
    2391      217725 :   vker = Flm_to_FlxV(vker, sv);
    2392      217725 :   d = lg(vker)-1;
    2393      217725 :   po2 = p >> 1; /* (p-1) / 2 */
    2394      217725 :   ir = 0;
    2395             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2396      759720 :   for (L=1; L<d; )
    2397             :   {
    2398      324270 :     GEN pol = zero_zv(l-2);
    2399      324270 :     pol[1] = sv;
    2400      324270 :     pol[2] = random_Fl(p); /*Assume vker[1]=1*/
    2401      965490 :     for (i=2; i<=d; i++)
    2402      641220 :       Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
    2403      324270 :     (void)Flx_renormalize(pol,l-1);
    2404             : 
    2405      779400 :     for (i=ir; i<L && L<d; i++)
    2406             :     {
    2407      455130 :       a = t[i]; la = degpol(a);
    2408      455130 :       if (la == 1) { set_irred(i); }
    2409      411624 :       else if (la == 2)
    2410             :       {
    2411       71304 :         ulong r = Flx_quad_root(a,p,1);
    2412       71304 :         if (r != p)
    2413             :         {
    2414       56208 :           t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
    2415       56208 :           t[L] = deg1_Flx(1, p - r, sv); L++;
    2416             :         }
    2417       71304 :         set_irred(i);
    2418             :       }
    2419             :       else
    2420             :       {
    2421      340320 :         pari_sp av = avma;
    2422      340320 :         b = Flx_rem(pol, a, p);
    2423      340320 :         if (degpol(b) <= 0) { avma=av; continue; }
    2424      325038 :         b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
    2425      325038 :         b = Flx_gcd(a,b, p); lb = degpol(b);
    2426      325038 :         if (lb && lb < la)
    2427             :         {
    2428      185145 :           b = Flx_normalize(b, p);
    2429      185145 :           t[L] = Flx_div(a,b,p);
    2430      185145 :           t[i]= b; L++;
    2431             :         }
    2432      139893 :         else avma = av;
    2433             :       }
    2434             :     }
    2435             :   }
    2436      217725 :   return d;
    2437             : }
    2438             : 
    2439             : static long
    2440          36 : FpX_split_Berlekamp(GEN *t, GEN p)
    2441             : {
    2442          36 :   GEN u = *t, a,b,po2,vker;
    2443          36 :   long d, i, ir, L, la, lb, vu = varn(u);
    2444          36 :   if (lgefint(p) == 3)
    2445             :   {
    2446           0 :     ulong up = p[2];
    2447           0 :     if (up == 2)
    2448             :     {
    2449           0 :       *t = ZX_to_F2x(*t);
    2450           0 :       d = F2x_split_Berlekamp(t);
    2451           0 :       for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
    2452             :     }
    2453             :     else
    2454             :     {
    2455           0 :       *t = ZX_to_Flx(*t, up);
    2456           0 :       d = Flx_split_Berlekamp(t, up);
    2457           0 :       for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
    2458             :     }
    2459           0 :     return d;
    2460             :   }
    2461          36 :   la = degpol(u);
    2462          36 :   if (la == 1) return 1;
    2463          36 :   if (la == 2)
    2464             :   {
    2465           0 :     GEN r = FpX_quad_root(u,p,1);
    2466           0 :     if (r)
    2467             :     {
    2468           0 :       t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
    2469           0 :       t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
    2470           0 :       return 2;
    2471             :     }
    2472           0 :     return 1;
    2473             :   }
    2474          36 :   vker = FpX_Berlekamp_ker(u,p);
    2475          36 :   vker = RgM_to_RgXV(vker,vu);
    2476          36 :   d = lg(vker)-1;
    2477          36 :   po2 = shifti(p, -1); /* (p-1) / 2 */
    2478          36 :   ir = 0;
    2479             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2480         102 :   for (L=1; L<d; )
    2481             :   {
    2482          30 :     GEN pol = scalar_ZX_shallow(randomi(p), vu);
    2483          84 :     for (i=2; i<=d; i++)
    2484          54 :       pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
    2485          30 :     pol = FpX_red(pol,p);
    2486          72 :     for (i=ir; i<L && L<d; i++)
    2487             :     {
    2488          42 :       a = t[i]; la = degpol(a);
    2489          42 :       if (la == 1) { set_irred(i); }
    2490          42 :       else if (la == 2)
    2491             :       {
    2492           0 :         GEN r = FpX_quad_root(a,p,1);
    2493           0 :         if (r)
    2494             :         {
    2495           0 :           t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
    2496           0 :           t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
    2497             :         }
    2498           0 :         set_irred(i);
    2499             :       }
    2500             :       else
    2501             :       {
    2502          42 :         pari_sp av = avma;
    2503          42 :         b = FpX_rem(pol, a, p);
    2504          42 :         if (degpol(b) <= 0) { avma=av; continue; }
    2505          42 :         b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
    2506          42 :         b = FpX_gcd(a,b, p); lb = degpol(b);
    2507          42 :         if (lb && lb < la)
    2508             :         {
    2509          24 :           b = FpX_normalize(b, p);
    2510          24 :           t[L] = FpX_div(a,b,p);
    2511          24 :           t[i]= b; L++;
    2512             :         }
    2513          18 :         else avma = av;
    2514             :       }
    2515             :     }
    2516             :   }
    2517          36 :   return d;
    2518             : }
    2519             : 
    2520             : static GEN
    2521      106848 : F2x_Berlekamp_i(GEN f, long flag)
    2522             : {
    2523      106848 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    2524             :   GEN y, E, t, V;
    2525             : 
    2526      106848 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    2527             : 
    2528       45504 :   val = F2x_valrem(f, &f);
    2529       45504 :   if (flag == 2 && val > 1) return NULL;
    2530       45504 :   V = F2x_factor_squarefree(f); lV = lg(V);
    2531       45504 :   if (flag == 2 && lV > 2) return NULL;
    2532             : 
    2533             :   /* to hold factors and exponents */
    2534       45504 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2535       45504 :   E = cgetg(d+1,t_VECSMALL);
    2536       45504 :   lfact = 1;
    2537       45504 :   if (val) {
    2538       14634 :     if (flag == 1)
    2539           0 :       t[1] = 1;
    2540             :     else
    2541       14634 :       gel(t,1) = polx_F2x(f[1]);
    2542       14634 :     E[1] = val; lfact++;
    2543             :   }
    2544             : 
    2545      180552 :   for (k=1; k<lV; k++)
    2546             :   {
    2547      135048 :     if (F2x_degree(gel(V, k))==0) continue;
    2548       42336 :     gel(t,lfact) = gel(V, k);
    2549       42336 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    2550       42336 :     if (flag == 2 && d != 1) return NULL;
    2551       42336 :     if (flag == 1)
    2552           0 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    2553       42336 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2554       42336 :     lfact += d;
    2555             :   }
    2556       45504 :   if (flag == 2) return gen_1; /* irreducible */
    2557       45504 :   y = FE_setlg(t,E, lfact);
    2558       45504 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2559       45504 :               : sort_factor_pol(y, cmpGuGu);
    2560             : }
    2561             : 
    2562             : static GEN
    2563      589103 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
    2564             : {
    2565      589103 :   long lfact, val, d = degpol(f), j, k, lV;
    2566             :   GEN y, E, t, V;
    2567             : 
    2568      589103 :   if (p == 2)
    2569             :   {
    2570       10368 :     GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
    2571       10368 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2572       10368 :     return F;
    2573             :   }
    2574      578735 :   if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
    2575      284181 :   val = Flx_valrem(f, &f);
    2576      284181 :   if (flag == 2 && val > 1) return NULL;
    2577      284181 :   V = Flx_factor_squarefree(f, p); lV = lg(V);
    2578      284181 :   if (flag == 2 && lV > 2) return NULL;
    2579             : 
    2580             :   /* to hold factors and exponents */
    2581      284181 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2582      284181 :   E = cgetg(d+1,t_VECSMALL);
    2583      284181 :   lfact = 1;
    2584      284181 :   if (val) {
    2585       11448 :     if (flag == 1)
    2586           0 :       t[1] = 1;
    2587             :     else
    2588       11448 :       gel(t,1) = polx_Flx(f[1]);
    2589       11448 :     E[1] = val; lfact++;
    2590             :   }
    2591             : 
    2592      711636 :   for (k=1; k<lV; k++)
    2593             :   {
    2594      427455 :     if (degpol(gel(V, k))==0) continue;
    2595      289683 :     gel(t,lfact) = Flx_normalize(gel(V, k), p);
    2596      289683 :     d = Flx_split_Berlekamp(&gel(t,lfact), p);
    2597      289683 :     if (flag == 2 && d != 1) return NULL;
    2598      289683 :     if (flag == 1)
    2599           0 :       for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
    2600      289683 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2601      289683 :     lfact += d;
    2602             :   }
    2603      284181 :   if (flag == 2) return gen_1; /* irreducible */
    2604      284181 :   y = FE_setlg(t,E,lfact);
    2605      284181 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2606      284181 :               : sort_factor_pol(y, cmpGuGu);
    2607             : }
    2608             : 
    2609             : /* f an FpX or an Flx */
    2610             : static GEN
    2611      654695 : FpX_Berlekamp_i(GEN f, GEN p, long flag)
    2612             : {
    2613      654695 :   long lfact, val, d = degpol(f), j, k, lV;
    2614             :   GEN y, E, t ,V;
    2615             : 
    2616      654695 :   if (typ(f) == t_VECSMALL)
    2617             :   {/* lgefint(p) == 3 */
    2618      654605 :     ulong pp = p[2];
    2619             :     GEN F;
    2620      654605 :     if (pp == 2) {
    2621       96468 :       F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
    2622       96468 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2623             :     } else {
    2624      558137 :       F = Flx_Berlekamp_i(f, pp, flag);
    2625      558137 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2626             :     }
    2627      654605 :     return F;
    2628             :   }
    2629             :   /* p is large (and odd) */
    2630          90 :   if (d <= 2) return FpX_factor_deg2(f, p, d, flag);
    2631          36 :   val = ZX_valrem(f, &f);
    2632          36 :   if (flag == 2 && val > 1) return NULL;
    2633          36 :   V = FpX_factor_Yun(f, p); lV = lg(V);
    2634          36 :   if (flag == 2 && lg(V) > 2) return NULL;
    2635             : 
    2636             :   /* to hold factors and exponents */
    2637          36 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2638          36 :   E = cgetg(d+1,t_VECSMALL);
    2639          36 :   lfact = 1;
    2640          36 :   if (val) {
    2641           0 :     if (flag == 1)
    2642           0 :       t[1] = 1;
    2643             :     else
    2644           0 :       gel(t,1) = pol_x(varn(f));
    2645           0 :     E[1] = val; lfact++;
    2646             :   }
    2647             : 
    2648          72 :   for (k=1; k<lV; k++)
    2649             :   {
    2650          36 :     if (degpol(gel(V,k))==0) continue;
    2651          36 :     gel(t,lfact) = FpX_normalize(gel(V, k), p);
    2652          36 :     d = FpX_split_Berlekamp(&gel(t,lfact), p);
    2653          36 :     if (flag == 2 && d != 1) return NULL;
    2654          36 :     if (flag == 1)
    2655           0 :       for (j=0; j<d; j++) t[lfact+j] = degpol(gel(t,lfact+j));
    2656          36 :     for (j=0; j<d; j++) E[lfact+j] = k;
    2657          36 :     lfact += d;
    2658             :   }
    2659          36 :   if (flag == 2) return gen_1; /* irreducible */
    2660          36 :   y = FE_setlg(t,E, lfact);
    2661          36 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2662          36 :               : sort_factor_pol(y, cmpii);
    2663             : }
    2664             : GEN
    2665      654635 : FpX_factor(GEN f, GEN p)
    2666             : {
    2667      654635 :   pari_sp av = avma;
    2668      654635 :   ZX_factmod_init(&f, p);
    2669      654635 :   return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
    2670             : }
    2671             : GEN
    2672       53898 : Flx_factor(GEN f, ulong p)
    2673             : {
    2674       53898 :   pari_sp av = avma;
    2675       53898 :   GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
    2676       53898 :   return gerepilecopy(av, F);
    2677             : }
    2678             : GEN
    2679          12 : F2x_factor(GEN f)
    2680             : {
    2681          12 :   pari_sp av = avma;
    2682          12 :   return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
    2683             : }
    2684             : 
    2685             : GEN
    2686          90 : factormod0(GEN f, GEN p, long flag)
    2687             : {
    2688          90 :   switch(flag)
    2689             :   {
    2690          66 :     case 0: return factmod(f,p);
    2691          24 :     case 1: return simplefactmod(f,p);
    2692           0 :     default: pari_err_FLAG("factormod");
    2693             :   }
    2694             :   return NULL; /* LCOV_EXCL_LINE */
    2695             : }

Generated by: LCOV version 1.11