Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 18906-a425091) Lines: 2154 2325 92.6 %
Date: 2016-05-05 Functions: 181 194 93.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1101 1422 77.4 %

           Branch data     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                 :        224 : factmod_init(GEN *F, GEN p)
      32                 :            : {
      33         [ -  + ]:        224 :   if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
      34         [ -  + ]:        224 :   if (signe(p) < 0) pari_err_PRIME("factmod",p);
      35         [ -  + ]:        224 :   if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
      36         [ +  + ]:        224 :   if (lgefint(p) == 3)
      37                 :            :   {
      38                 :        126 :     ulong pp = p[2];
      39         [ -  + ]:        126 :     if (pp < 2) pari_err_PRIME("factmod", p);
      40                 :        126 :     *F = RgX_to_Flx(*F, pp);
      41         [ +  + ]:        126 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42                 :            :   }
      43                 :            :   else
      44                 :            :   {
      45                 :         98 :     *F = RgX_to_FpX(*F, p);
      46         [ +  + ]:         98 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      47                 :            :   }
      48                 :        224 : }
      49                 :            : /* as above, assume p prime and *F a ZX */
      50                 :            : static void
      51                 :     607682 : ZX_factmod_init(GEN *F, GEN p)
      52                 :            : {
      53         [ +  + ]:     607682 :   if (lgefint(p) == 3)
      54                 :            :   {
      55                 :     601153 :     ulong pp = p[2];
      56                 :     601153 :     *F = ZX_to_Flx(*F, pp);
      57         [ +  + ]:     601153 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      58                 :            :   }
      59                 :            :   else
      60                 :            :   {
      61                 :       6529 :     *F = FpX_red(*F, p);
      62         [ +  + ]:       6529 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      63                 :            :   }
      64                 :     607682 : }
      65                 :            : 
      66                 :            : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      67                 :            : static GEN
      68                 :         35 : all_roots_mod_p(ulong p, int not_0)
      69                 :            : {
      70                 :            :   GEN r;
      71                 :            :   ulong i;
      72         [ +  + ]:         35 :   if (not_0) {
      73                 :         21 :     r = cgetg(p, t_VECSMALL);
      74         [ +  + ]:         63 :     for (i = 1; i < p; i++) r[i] = i;
      75                 :            :   } else {
      76                 :         14 :     r = cgetg(p+1, t_VECSMALL);
      77         [ +  + ]:         84 :     for (i = 0; i < p; i++) r[i+1] = i;
      78                 :            :   }
      79                 :         35 :   return r;
      80                 :            : }
      81                 :            : 
      82                 :            : /* X^n - 1 */
      83                 :            : static GEN
      84                 :         20 : Flx_Xnm1(long sv, long n, ulong p)
      85                 :            : {
      86                 :         20 :   GEN t = cgetg(n+3, t_VECSMALL);
      87                 :            :   long i;
      88                 :         20 :   t[1] = sv;
      89                 :         20 :   t[2] = p - 1;
      90         [ +  + ]:         81 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      91                 :         20 :   t[i] = 1; return t;
      92                 :            : }
      93                 :            : /* X^n + 1 */
      94                 :            : static GEN
      95                 :         13 : Flx_Xn1(long sv, long n, ulong p)
      96                 :            : {
      97                 :         13 :   GEN t = cgetg(n+3, t_VECSMALL);
      98                 :            :   long i;
      99                 :            :   (void) p;
     100                 :         13 :   t[1] = sv;
     101                 :         13 :   t[2] = 1;
     102         [ +  + ]:         53 :   for (i = 3; i <= n+1; i++) t[i] = 0;
     103                 :         13 :   t[i] = 1; return t;
     104                 :            : }
     105                 :            : 
     106                 :            : static ulong
     107                 :          6 : Fl_nonsquare(ulong p)
     108                 :            : {
     109                 :          6 :   long k = 2;
     110                 :          6 :   for (;; k++)
     111                 :            :   {
     112                 :         12 :     long i = krouu(k, p);
     113         [ -  + ]:         12 :     if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
     114         [ +  + ]:         12 :     if (i < 0) return k;
     115                 :          6 :   }
     116                 :            : }
     117                 :            : 
     118                 :            : /* f monic Flx, f(0) != 0. Return a monic squarefree g with the same
     119                 :            :  * roots as f */
     120                 :            : static GEN
     121                 :         14 : Flx_cut_out_roots(GEN f, ulong p)
     122                 :            : {
     123                 :         14 :   GEN g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
     124 [ +  + ][ +  - ]:         14 :   if (g != f && degpol(g) >= 0) {
     125                 :          7 :     (void)Flx_valrem(g, &g); /* reduction may introduce 0 root */
     126                 :          7 :     g = Flx_gcd(g, Flx_Xnm1(g[1], p-1, p), p);
     127                 :          7 :     g = Flx_normalize(g, p);
     128                 :            :   }
     129                 :         14 :   return g;
     130                 :            : }
     131                 :            : 
     132                 :            : /* by checking f(0..p-1) */
     133                 :            : GEN
     134                 :         14 : Flx_roots_naive(GEN f, ulong p)
     135                 :            : {
     136                 :         14 :   long d, n = 0;
     137                 :         14 :   ulong s = 1UL, r;
     138                 :         14 :   GEN q, y = cgetg(degpol(f) + 1, t_VECSMALL);
     139                 :         14 :   pari_sp av2, av = avma;
     140                 :            : 
     141         [ -  + ]:         14 :   if (Flx_valrem(f, &f)) y[++n] = 0;
     142                 :         14 :   f = Flx_cut_out_roots(f, p);
     143                 :         14 :   d = degpol(f);
     144         [ -  + ]:         14 :   if (d < 0) return all_roots_mod_p(p, n == 0);
     145                 :         14 :   av2 = avma;
     146         [ +  + ]:        294 :   while (d > 1) /* d = current degree of f */
     147                 :            :   {
     148                 :        287 :     q = Flx_div_by_X_x(f, s, p, &r); /* TODO: FFT-type multi-evaluation */
     149         [ +  + ]:        287 :     if (r) avma = av2; else { y[++n] = s; d--; f = q; av2 = avma; }
     150         [ +  + ]:        287 :     if (++s == p) break;
     151                 :            :   }
     152         [ +  + ]:         14 :   if (d == 1)
     153                 :            :   { /* -f[2]/f[3], root of deg 1 polynomial */
     154                 :          7 :     r = Fl_mul(p - Fl_inv(f[3], p), f[2], p);
     155         [ +  - ]:          7 :     if (r >= s) y[++n] = r; /* otherwise double root */
     156                 :            :   }
     157                 :         14 :   avma = av; fixlg(y, n+1); return y;
     158                 :            : }
     159                 :            : static GEN
     160                 :       6328 : Flx_root_mod_2(GEN f)
     161                 :            : {
     162                 :       6328 :   int z1, z0 = !(f[2] & 1);
     163                 :            :   long i,n;
     164                 :            :   GEN y;
     165                 :            : 
     166         [ +  + ]:      24472 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     167                 :       6328 :   z1 = n & 1;
     168                 :       6328 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     169         [ +  + ]:       6328 :   if (z0) y[i++] = 0;
     170         [ +  + ]:       6328 :   if (z1) y[i  ] = 1;
     171                 :       6328 :   return y;
     172                 :            : }
     173                 :            : static ulong
     174                 :         14 : Flx_oneroot_mod_2(GEN f)
     175                 :            : {
     176                 :            :   long i,n;
     177         [ -  + ]:         14 :   if (!(f[2] & 1)) return 0;
     178         [ +  + ]:         56 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     179         [ +  + ]:         14 :   if (n & 1) return 1;
     180                 :         14 :   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         [ +  + ]:    1816350 : 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                 :      60809 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
     196                 :            : {
     197                 :      60809 :   pari_sp av = avma;
     198                 :            :   GEN y;
     199         [ +  + ]:      60809 :   if (gpwrap)
     200                 :         70 :     factmod_init(&f, pp);
     201                 :            :   else
     202                 :      60739 :     ZX_factmod_init(&f, pp);
     203      [ +  +  + ]:      60809 :   switch(lg(f))
     204                 :            :   {
     205                 :         14 :     case 2: pari_err_ROOTS0("rootmod");
     206                 :         35 :     case 3: avma = av; return cgetg(1,t_COL);
     207                 :            :   }
     208         [ +  + ]:      60760 :   if (typ(f) == t_VECSMALL)
     209                 :            :   {
     210                 :      57970 :     ulong p = pp[2];
     211         [ +  + ]:      57970 :     if (p == 2)
     212                 :       6328 :       y = Flx_root_mod_2(f);
     213                 :            :     else
     214                 :            :     {
     215         [ -  + ]:      51642 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     216                 :      51642 :       y = Roots(f, p);
     217                 :            :     }
     218                 :      57970 :     y = Flc_to_ZC(y);
     219                 :            :   }
     220                 :            :   else
     221                 :       2790 :     y = FpX_roots_i(f, pp);
     222         [ +  + ]:      60760 :   if (gpwrap) y = FpC_to_mod(y, pp);
     223                 :      60795 :   return gerepileupto(av, y);
     224                 :            : }
     225                 :            : /* assume that f is a ZX an pp a prime */
     226                 :            : GEN
     227                 :      60739 : FpX_roots(GEN f, GEN pp)
     228                 :      60739 : { return rootmod_aux(f, pp, Flx_roots_i, 0); }
     229                 :            : /* no assumptions on f and pp */
     230                 :            : GEN
     231                 :         14 : rootmod2(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_naive, 1); }
     232                 :            : GEN
     233                 :         56 : rootmod(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_i, 1); }
     234                 :            : GEN
     235                 :         70 : rootmod0(GEN f, GEN p, long flag)
     236                 :            : {
     237      [ +  +  - ]:         70 :   switch(flag)
     238                 :            :   {
     239                 :         56 :     case 0: return rootmod(f,p);
     240                 :         14 :     case 1: return rootmod2(f,p);
     241                 :          0 :     default: pari_err_FLAG("polrootsmod");
     242                 :            :   }
     243                 :         56 :   return NULL; /* not reached */
     244                 :            : }
     245                 :            : 
     246                 :            : /* assume x reduced mod p > 2, monic. */
     247                 :            : static int
     248                 :         21 : FpX_quad_factortype(GEN x, GEN p)
     249                 :            : {
     250                 :         21 :   GEN b = gel(x,3), c = gel(x,2);
     251                 :         21 :   GEN D = subii(sqri(b), shifti(c,2));
     252                 :         21 :   return kronecker(D,p);
     253                 :            : }
     254                 :            : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     255                 :            : GEN
     256                 :       6475 : FpX_quad_root(GEN x, GEN p, int unknown)
     257                 :            : {
     258                 :       6475 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     259                 :            : 
     260         [ -  + ]:       6475 :   if (equaliu(p, 2)) {
     261         [ #  # ]:          0 :     if (!signe(b)) return c;
     262         [ #  # ]:          0 :     return signe(c)? NULL: gen_1;
     263                 :            :   }
     264                 :       6475 :   D = subii(sqri(b), shifti(c,2));
     265                 :       6475 :   D = remii(D,p);
     266 [ +  + ][ +  + ]:       6475 :   if (unknown && kronecker(D,p) == -1) return NULL;
     267                 :            : 
     268                 :       6467 :   s = Fp_sqrt(D,p);
     269                 :            :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     270         [ +  + ]:       6467 :   if (!s) return NULL;
     271                 :       6475 :   return Fp_halve(Fp_sub(s,b, p), p);
     272                 :            : }
     273                 :            : static GEN
     274                 :       2871 : FpX_otherroot(GEN x, GEN r, GEN p)
     275                 :       2871 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     276                 :            : 
     277                 :            : /* disc(x^2+bx+c) = b^2 - 4c */
     278                 :            : static ulong
     279                 :    5974397 : Fl_disc_bc(ulong b, ulong c, ulong p)
     280                 :    5974397 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     281                 :            : /* p > 2 */
     282                 :            : static ulong
     283                 :    5900241 : Flx_quad_root(GEN x, ulong p, int unknown)
     284                 :            : {
     285                 :    5900241 :   ulong s, b = x[3], c = x[2];
     286                 :    5900241 :   ulong D = Fl_disc_bc(b, c, p);
     287 [ +  + ][ +  + ]:    5900649 :   if (unknown && krouu(D,p) == -1) return p;
     288                 :    5786942 :   s = Fl_sqrt(D,p);
     289         [ +  + ]:    5776113 :   if (s==~0UL) return p;
     290                 :    5889266 :   return Fl_halve(Fl_sub(s,b, p), p);
     291                 :            : }
     292                 :            : static ulong
     293                 :    5249914 : Flx_otherroot(GEN x, ulong r, ulong p)
     294                 :    5249914 : { 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                 :     614047 : split_init(struct split_t *S, long max)
     302                 :            : {
     303                 :     614047 :   S->todo = vectrunc_init(max);
     304                 :     614048 :   S->done = vectrunc_init(max);
     305                 :     614053 : }
     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                 :     671915 : 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                 :     671906 : split_moveto_done(struct split_t *S, long i, GEN t)
     323                 :            : {
     324                 :     671906 :   long n = lg(S->todo)-1;
     325                 :     671906 :   vectrunc_append(S->done, t);
     326         [ +  - ]:     671901 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     327                 :     671901 :   setlg(S->todo, n);
     328                 :            : 
     329                 :     671902 : }
     330                 :            : /* append t to done */
     331                 :            : static void
     332                 :      99933 : split_add_done(struct split_t *S, GEN t)
     333                 :      99933 : { vectrunc_append(S->done, t); }
     334                 :            : /* split todo[i] into a and b */
     335                 :            : static void
     336                 :      75053 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     337                 :            : {
     338                 :      75053 :   gel(S->todo, i) = a;
     339                 :      75053 :   split_add(S, b);
     340                 :      75054 : }
     341                 :            : /* split todo[i] into a and b, moved to done */
     342                 :            : static void
     343                 :      81033 : split_done(struct split_t *S, long i, GEN a, GEN b)
     344                 :            : {
     345                 :      81033 :   split_moveto_done(S, i, a);
     346                 :      81031 :   split_add_done(S, b);
     347                 :      81032 : }
     348                 :            : 
     349                 :            : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     350                 :            : static GEN
     351                 :       2790 : FpX_roots_i(GEN f, GEN p)
     352                 :            : {
     353                 :            :   GEN pol, pol0, a, q;
     354                 :            :   struct split_t S;
     355                 :            : 
     356                 :       2790 :   split_init(&S, lg(f)-1);
     357                 :       2790 :   settyp(S.done, t_COL);
     358         [ +  + ]:       2790 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     359   [ +  +  +  + ]:       2790 :   switch(degpol(f))
     360                 :            :   {
     361                 :          7 :     case 0: return ZC_copy(S.done);
     362                 :         14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     363                 :            :     case 2: {
     364                 :       1603 :       GEN s, r = FpX_quad_root(f, p, 1);
     365         [ +  - ]:       1603 :       if (r) {
     366                 :       1603 :         split_add_done(&S, r);
     367                 :       1603 :         s = FpX_otherroot(f,r, p);
     368                 :            :         /* f not known to be square free yet */
     369         [ +  - ]:       1603 :         if (!equalii(r, s)) split_add_done(&S, s);
     370                 :            :       }
     371                 :       1603 :       return sort(S.done);
     372                 :            :     }
     373                 :            :   }
     374                 :            : 
     375                 :       1166 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     376         [ -  + ]:       1166 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     377                 :       1166 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     378                 :       1166 :   a = FpX_gcd(f,a, p);
     379         [ -  + ]:       1166 :   if (!degpol(a)) return ZC_copy(S.done);
     380                 :       1166 :   split_add(&S, FpX_normalize(a,p));
     381                 :            : 
     382                 :       1166 :   q = shifti(p,-1);
     383                 :       1166 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     384                 :       1166 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     385                 :       1166 :   for (pol0[2] = 1;; pol0[2]++)
     386                 :            :   {
     387                 :       2430 :     long j, l = lg(S.todo);
     388         [ +  + ]:       2430 :     if (l == 1) return sort(S.done);
     389 [ -  + ][ #  # ]:       1264 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     390         [ +  + ]:       2619 :     for (j = 1; j < l; j++)
     391                 :            :     {
     392                 :       1355 :       GEN c = gel(S.todo,j);
     393      [ +  +  + ]:       1355 :       switch(degpol(c))
     394                 :            :       { /* convert linear and quadratics to roots, try to split the rest */
     395                 :            :         case 1:
     396                 :         56 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     397                 :         56 :           j--; l--; break;
     398                 :            :         case 2: {
     399                 :       1201 :           GEN r = FpX_quad_root(c, p, 0), s = FpX_otherroot(c,r, p);
     400                 :       1201 :           split_done(&S, j, r, s);
     401                 :       1201 :           j--; l--; break;
     402                 :            :         }
     403                 :            :         default: {
     404                 :            :           /* b = pol^(p-1)/2 - 1 */
     405                 :         98 :           GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
     406                 :            :           long db;
     407                 :         98 :           b = FpX_gcd(c,b, p); db = degpol(b);
     408 [ +  - ][ +  + ]:         98 :           if (db && db < degpol(c))
     409                 :            :           {
     410                 :         91 :             b = FpX_normalize(b, p);
     411                 :         91 :             c = FpX_div(c,b, p);
     412                 :         91 :             split_todo(&S, j, b, c);
     413                 :            :           }
     414                 :            :         }
     415                 :            :       }
     416                 :            :     }
     417                 :       4054 :   }
     418                 :            : }
     419                 :            : 
     420                 :            : /* Assume f is normalized */
     421                 :            : static ulong
     422                 :      53915 : Flx_cubic_root(GEN ff, ulong p)
     423                 :            : {
     424                 :      53915 :   GEN f = Flx_normalize(ff,p);
     425                 :      53915 :   ulong pi = get_Fl_red(p);
     426         [ +  + ]:      53915 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     427                 :      53915 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     428                 :      53915 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     429                 :      53913 :   ulong B = Fl_addmul_pre(t, Fl_sub(Fl_double(t2, p), b, p), c, p, pi);
     430                 :      53915 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     431                 :      53915 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     432                 :      53915 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     433                 :      53913 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     434                 :      53915 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     435         [ +  + ]:      53915 :   if (s!=~0UL)
     436                 :            :   {
     437         [ +  + ]:      33392 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     438         [ +  + ]:      33392 :     if (p%3==2) /* 1 solutions */
     439                 :       5166 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     440                 :            :     else
     441                 :            :     {
     442                 :      28226 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     443         [ +  + ]:      28226 :       if (vS1==~0UL) return p; /*0 solutions*/
     444                 :            :       /*3 solutions*/
     445                 :            :     }
     446         [ +  + ]:      22680 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     447                 :      22680 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     448                 :            :   }
     449                 :            :   else
     450                 :            :   {
     451                 :      20523 :     pari_sp av = avma;
     452                 :      20523 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     453                 :      20523 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     454                 :            :     ulong Sa;
     455         [ +  + ]:      20523 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     456                 :      19583 :     Sa = vS1[1];
     457         [ +  + ]:      19583 :     if (p%3==1) /*1 solutions*/
     458                 :            :     {
     459                 :      14132 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     460         [ +  + ]:      14132 :       if (Fa!=P)
     461                 :       9542 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     462                 :            :     }
     463                 :      19583 :     avma = av;
     464                 :      53915 :     return Fl_sub(Fl_double(Sa,p),t,p);
     465                 :            :   }
     466                 :            : }
     467                 :            : 
     468                 :            : /* assume p > 2 prime */
     469                 :            : static ulong
     470                 :     694288 : Flx_oneroot_i(GEN f, ulong p, long fl)
     471                 :            : {
     472                 :            :   GEN pol, a;
     473                 :            :   ulong q;
     474                 :            :   long da;
     475                 :            : 
     476         [ +  + ]:     694288 :   if (Flx_val(f)) return 0;
     477   [ +  +  +  + ]:     694123 :   switch(degpol(f))
     478                 :            :   {
     479                 :        190 :     case 1: return Fl_neg(f[2], p);
     480                 :     540348 :     case 2: return Flx_quad_root(f, p, 1);
     481         [ +  - ]:      50102 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     482                 :            :   }
     483                 :            : 
     484         [ +  + ]:     103474 :   if (!fl)
     485                 :            :   {
     486                 :      93113 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     487         [ -  + ]:      93113 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     488                 :      93113 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     489                 :      93113 :     a = Flx_gcd(f,a, p);
     490                 :      10361 :   } else a = f;
     491                 :     103474 :   da = degpol(a);
     492         [ +  + ]:     103478 :   if (!da) return p;
     493                 :      53583 :   a = Flx_normalize(a,p);
     494                 :            : 
     495                 :      53595 :   q = p >> 1;
     496                 :      53595 :   pol = polx_Flx(f[1]);
     497                 :      53580 :   for(pol[2] = 1;; pol[2]++)
     498                 :            :   {
     499 [ -  + ][ #  # ]:      90759 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     500   [ +  +  +  + ]:      90757 :     switch(da)
     501                 :            :     {
     502                 :      22191 :       case 1: return Fl_neg(a[2], p);
     503                 :      27600 :       case 2: return Flx_quad_root(a, p, 0);
     504         [ +  - ]:       3813 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     505                 :            :       default: {
     506                 :      37153 :         GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
     507                 :            :         long db;
     508                 :      37175 :         b = Flx_gcd(a,b, p); db = degpol(b);
     509 [ +  + ][ +  + ]:      37178 :         if (db && db < da)
     510                 :            :         {
     511                 :      35395 :           b = Flx_normalize(b, p);
     512         [ +  + ]:      35395 :           if (db <= (da >> 1)) {
     513                 :      21897 :             a = b;
     514                 :      21897 :             da = db;
     515                 :            :           } else {
     516                 :      13498 :             a = Flx_div(a,b, p);
     517                 :      13499 :             da -= db;
     518                 :            :           }
     519                 :            :         }
     520                 :            :       }
     521                 :            :     }
     522                 :     731491 :   }
     523                 :            : }
     524                 :            : 
     525                 :            : /* assume p > 2 prime */
     526                 :            : static GEN
     527                 :       3649 : FpX_oneroot_i(GEN f, GEN p)
     528                 :            : {
     529                 :            :   GEN pol, pol0, a, q;
     530                 :            :   long da;
     531                 :            : 
     532         [ +  + ]:       3649 :   if (ZX_val(f)) return gen_0;
     533      [ +  +  + ]:       3642 :   switch(degpol(f))
     534                 :            :   {
     535                 :        140 :     case 1: return subii(p, gel(f,2));
     536                 :       3432 :     case 2: return FpX_quad_root(f, p, 1);
     537                 :            :   }
     538                 :            : 
     539                 :         70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     540         [ -  + ]:         70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     541                 :         70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     542                 :         70 :   a = FpX_gcd(f,a, p);
     543                 :         70 :   da = degpol(a);
     544         [ -  + ]:         70 :   if (!da) return NULL;
     545                 :         70 :   a = FpX_normalize(a,p);
     546                 :            : 
     547                 :         70 :   q = shifti(p,-1);
     548                 :         70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     549                 :         70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     550                 :         70 :   for (pol0[2]=1; ; pol0[2]++)
     551                 :            :   {
     552 [ -  + ][ #  # ]:        224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     553      [ +  +  + ]:        224 :     switch(da)
     554                 :            :     {
     555                 :         42 :       case 1: return subii(p, gel(a,2));
     556                 :         28 :       case 2: return FpX_quad_root(a, p, 0);
     557                 :            :       default: {
     558                 :        154 :         GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, a,p), gen_1, p);
     559                 :            :         long db;
     560                 :        154 :         b = FpX_gcd(a,b, p); db = degpol(b);
     561 [ +  + ][ +  - ]:        154 :         if (db && db < da)
     562                 :            :         {
     563                 :        147 :           b = FpX_normalize(b, p);
     564         [ +  + ]:        147 :           if (db <= (da >> 1)) {
     565                 :        105 :             a = b;
     566                 :        105 :             da = db;
     567                 :            :           } else {
     568                 :         42 :             a = FpX_div(a,b, p);
     569                 :         42 :             da -= db;
     570                 :            :           }
     571                 :            :         }
     572                 :            :       }
     573                 :            :     }
     574                 :       3803 :   }
     575                 :            : }
     576                 :            : 
     577                 :            : ulong
     578                 :     676500 : Flx_oneroot(GEN f, ulong p)
     579                 :            : {
     580                 :     676500 :   pari_sp av = avma;
     581                 :            :   ulong r;
     582      [ -  -  + ]:     676500 :   switch(lg(f))
     583                 :            :   {
     584                 :          0 :     case 2: return 0;
     585                 :          0 :     case 3: avma = av; return p;
     586                 :            :   }
     587         [ -  + ]:     676500 :   if (p == 2) return Flx_oneroot_mod_2(f);
     588                 :     676500 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     589                 :     676500 :   avma = av; return r;
     590                 :            : }
     591                 :            : 
     592                 :            : ulong
     593                 :      10447 : Flx_oneroot_split(GEN f, ulong p)
     594                 :            : {
     595                 :      10447 :   pari_sp av = avma;
     596                 :            :   ulong r;
     597      [ -  -  + ]:      10447 :   switch(lg(f))
     598                 :            :   {
     599                 :          0 :     case 2: return 0;
     600                 :          0 :     case 3: avma = av; return p;
     601                 :            :   }
     602         [ -  + ]:      10447 :   if (p == 2) return Flx_oneroot_mod_2(f);
     603                 :      10447 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     604                 :      10471 :   avma = av; return r;
     605                 :            : }
     606                 :            : 
     607                 :            : /* assume that p is prime */
     608                 :            : GEN
     609                 :      11004 : FpX_oneroot(GEN f, GEN pp) {
     610                 :      11004 :   pari_sp av = avma;
     611                 :      11004 :   ZX_factmod_init(&f, pp);
     612      [ -  -  + ]:      11004 :   switch(lg(f))
     613                 :            :   {
     614                 :          0 :     case 2: avma = av; return gen_0;
     615                 :          0 :     case 3: avma = av; return NULL;
     616                 :            :   }
     617         [ +  + ]:      11004 :   if (typ(f) == t_VECSMALL)
     618                 :            :   {
     619                 :       7355 :     ulong r, p = pp[2];
     620         [ +  + ]:       7355 :     if (p == 2)
     621                 :         14 :       r = Flx_oneroot_mod_2(f);
     622                 :            :     else
     623                 :       7341 :       r = Flx_oneroot_i(f, p, 0);
     624                 :       7355 :     avma = av;
     625         [ +  + ]:       7355 :     return (r == p)? NULL: utoi(r);
     626                 :            :   }
     627                 :       3649 :   f = FpX_oneroot_i(f, pp);
     628         [ -  + ]:       3649 :   if (!f) { avma = av; return NULL; }
     629                 :      11004 :   return gerepileuptoint(av, f);
     630                 :            : }
     631                 :            : 
     632                 :            : /*******************************************************************/
     633                 :            : /*                                                                 */
     634                 :            : /*                     FACTORISATION MODULO p                      */
     635                 :            : /*                                                                 */
     636                 :            : /*******************************************************************/
     637                 :            : 
     638                 :            : /* Functions giving information on the factorisation. */
     639                 :            : 
     640                 :            : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
     641                 :            : GEN
     642                 :         55 : FpX_Berlekamp_ker(GEN u, GEN p)
     643                 :            : {
     644                 :         55 :   pari_sp ltop=avma;
     645                 :         55 :   long j,N = degpol(u);
     646                 :         55 :   GEN Q  = FpX_matFrobenius(u, p);
     647         [ +  + ]:       2139 :   for (j=1; j<=N; j++)
     648                 :       2084 :     gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
     649                 :         55 :   return gerepileupto(ltop, FpM_ker(Q,p));
     650                 :            : }
     651                 :            : 
     652                 :            : GEN
     653                 :      20825 : F2x_Berlekamp_ker(GEN u)
     654                 :            : {
     655                 :      20825 :   pari_sp ltop=avma;
     656                 :      20825 :   long j,N = F2x_degree(u);
     657                 :            :   GEN Q;
     658                 :            :   pari_timer T;
     659                 :      20825 :   timer_start(&T);
     660                 :      20825 :   Q = F2x_matFrobenius(u);
     661         [ +  + ]:     107884 :   for (j=1; j<=N; j++)
     662                 :      87059 :     F2m_flip(Q,j,j);
     663         [ -  + ]:      20825 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     664                 :      20825 :   Q = F2m_ker_sp(Q,0);
     665         [ -  + ]:      20825 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     666                 :      20825 :   return gerepileupto(ltop,Q);
     667                 :            : }
     668                 :            : 
     669                 :            : GEN
     670                 :     205526 : Flx_Berlekamp_ker(GEN u, ulong l)
     671                 :            : {
     672                 :     205526 :   pari_sp ltop=avma;
     673                 :     205526 :   long j,N = degpol(u);
     674                 :            :   GEN Q;
     675                 :            :   pari_timer T;
     676                 :     205526 :   timer_start(&T);
     677                 :     205526 :   Q  = Flx_matFrobenius(u, l);
     678         [ +  + ]:     967164 :   for (j=1; j<=N; j++)
     679                 :     761638 :     coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
     680         [ -  + ]:     205526 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     681                 :     205526 :   Q = Flm_ker_sp(Q,l,0);
     682         [ -  + ]:     205526 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     683                 :     205526 :   return gerepileupto(ltop,Q);
     684                 :            : }
     685                 :            : 
     686                 :            : /* product of terms of degree 1 in factorization of f */
     687                 :            : GEN
     688                 :      60368 : FpX_split_part(GEN f, GEN p)
     689                 :            : {
     690                 :      60368 :   long n = degpol(f);
     691                 :      60368 :   GEN z, X = pol_x(varn(f));
     692         [ +  + ]:      60368 :   if (n <= 1) return f;
     693                 :      60340 :   f = FpX_red(f, p);
     694                 :      60340 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     695                 :      60368 :   return FpX_gcd(z,f,p);
     696                 :            : }
     697                 :            : 
     698                 :            : /* Compute the number of roots in Fp without counting multiplicity
     699                 :            :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     700                 :            : long
     701                 :      18452 : FpX_nbroots(GEN f, GEN p)
     702                 :            : {
     703                 :      18452 :   pari_sp av = avma;
     704                 :      18452 :   GEN z = FpX_split_part(f, p);
     705                 :      18452 :   avma = av; return degpol(z);
     706                 :            : }
     707                 :            : 
     708                 :            : int
     709                 :          0 : FpX_is_totally_split(GEN f, GEN p)
     710                 :            : {
     711                 :          0 :   long n=degpol(f);
     712                 :          0 :   pari_sp av = avma;
     713         [ #  # ]:          0 :   if (n <= 1) return 1;
     714         [ #  # ]:          0 :   if (cmpui(n, p) > 0) return 0;
     715                 :          0 :   f = FpX_red(f, p);
     716                 :          0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     717                 :            : }
     718                 :            : 
     719                 :            : long
     720                 :    8118919 : Flx_nbroots(GEN f, ulong p)
     721                 :            : {
     722                 :    8118919 :   long n = degpol(f);
     723                 :    8116919 :   pari_sp av = avma;
     724                 :            :   GEN z;
     725         [ +  + ]:    8116919 :   if (n <= 1) return n;
     726         [ +  + ]:    8116744 :   if (n == 2)
     727                 :            :   {
     728                 :            :     ulong D;
     729         [ +  + ]:    4894192 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     730                 :    4893359 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     731                 :    4899093 :     return 1 + krouu(D,p);
     732                 :            :   }
     733                 :    3222552 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     734                 :    3222531 :   z = Flx_gcd(z, f, p);
     735                 :    8121543 :   avma = av; return degpol(z);
     736                 :            : }
     737                 :            : 
     738                 :            : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     739                 :            : static GEN
     740                 :       4102 : FpX_ddf(GEN T, GEN XP, GEN p)
     741                 :            : {
     742                 :       4102 :   pari_sp av = avma;
     743                 :            :   GEN b, g, h, F, f, Tr, xq;
     744                 :            :   long i, j, n, v;
     745                 :            :   long B, l, m;
     746                 :            :   pari_timer ti;
     747                 :       4102 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     748         [ -  + ]:       4102 :   if (n == 0) return cgetg(1, t_VEC);
     749         [ +  + ]:       4102 :   if (n == 1) return mkvec(get_FpX_mod(T));
     750                 :       4095 :   B = n/2;
     751                 :       4095 :   l = usqrt(B);
     752                 :       4095 :   m = (B+l-1)/l;
     753                 :       4095 :   T = FpX_get_red(T, p);
     754                 :       4095 :   b = cgetg(l+2, t_VEC);
     755                 :       4095 :   gel(b, 1) = pol_x(v);
     756                 :       4095 :   gel(b, 2) = XP;
     757         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     758                 :       4095 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     759         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
     760         [ +  + ]:      10661 :   for (i = 3; i <= l+1; i++)
     761                 :       6566 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     762         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
     763                 :       4095 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     764         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
     765                 :       4095 :   g = cgetg(m+1, t_VEC);
     766                 :       4095 :   gel(g, 1) = gel(xq, 2);
     767         [ +  + ]:      14931 :   for(i = 2; i <= m; i++)
     768                 :      10836 :     gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     769         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
     770                 :       4095 :   h = cgetg(m+1, t_VEC);
     771         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     772                 :            :   {
     773                 :      14931 :     pari_sp av = avma;
     774                 :      14931 :     GEN gj = gel(g, j);
     775                 :      14931 :     GEN e = FpX_sub(gj, gel(b, 1), p);
     776         [ +  + ]:      42938 :     for (i = 2; i <= l; i++)
     777                 :      28007 :       e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
     778                 :      14931 :     gel(h, j) = gerepileupto(av, e);
     779                 :            :   }
     780         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
     781                 :       4095 :   Tr = get_FpX_mod(T);
     782                 :       4095 :   F = cgetg(m+1, t_VEC);
     783         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     784                 :            :   {
     785                 :      14931 :     gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
     786                 :      14931 :     Tr = FpX_div(Tr, gel(F,j), p);
     787                 :            :   }
     788         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
     789                 :       4095 :   f = const_vec(n, pol_1(v));
     790         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     791                 :            :   {
     792                 :      14931 :     GEN e = gel(F, j);
     793         [ +  - ]:      16919 :     for (i=l-1; i >= 0; i--)
     794                 :            :     {
     795                 :      16919 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     796         [ +  + ]:      16919 :       if (degpol(u))
     797                 :            :       {
     798                 :       2709 :         gel(f, l*j-i) = u;
     799                 :       2709 :         e = FpX_div(e, u, p);
     800                 :            :       }
     801         [ +  + ]:      16919 :       if (!degpol(e)) break;
     802                 :            :     }
     803                 :            :   }
     804         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
     805         [ +  + ]:       4095 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     806                 :       4102 :   return gerepilecopy(av, f);
     807                 :            : }
     808                 :            : 
     809                 :            : static void
     810                 :          0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     811                 :            : {
     812                 :          0 :   long n = degpol(Tp), r = n/d;
     813                 :            :   GEN T, f, ff;
     814                 :            :   GEN p2;
     815         [ #  # ]:          0 :   if (r==1) { gel(V, idx) = Tp; return; }
     816                 :          0 :   p2 = shifti(p,-1);
     817                 :          0 :   T = FpX_get_red(Tp, p);
     818                 :          0 :   XP = FpX_rem(XP, T, p);
     819                 :            :   while (1)
     820                 :            :   {
     821                 :          0 :     pari_sp btop = avma;
     822                 :            :     long i;
     823                 :          0 :     GEN g = random_FpX(n, varn(Tp), p);
     824                 :          0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     825         [ #  # ]:          0 :     if (signe(t) == 0) continue;
     826         [ #  # ]:          0 :     for(i=1; i<=10; i++)
     827                 :            :     {
     828                 :          0 :       pari_sp btop2 = avma;
     829                 :          0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     830                 :          0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     831 [ #  # ][ #  # ]:          0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     832                 :          0 :       avma = btop2;
     833                 :            :     }
     834 [ #  # ][ #  # ]:          0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     835                 :          0 :     avma = btop;
     836                 :          0 :   }
     837                 :          0 :   f = FpX_normalize(f, p);
     838                 :          0 :   ff = FpX_div(Tp, f ,p);
     839                 :          0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     840                 :          0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     841                 :            : }
     842                 :            : 
     843                 :            : static void
     844                 :        210 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     845                 :            : {
     846                 :            :   pari_sp av;
     847                 :        210 :   GEN Tp = get_FpX_mod(T);
     848                 :        210 :   long n = degpol(hp), vT = varn(Tp);
     849                 :            :   GEN u1, u2, f1, f2;
     850                 :            :   GEN R, h;
     851                 :        210 :   h = FpX_get_red(hp, p);
     852                 :        210 :   t = FpX_rem(t, T, p);
     853                 :        210 :   av = avma;
     854                 :            :   do
     855                 :            :   {
     856                 :        238 :     avma = av;
     857                 :        238 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     858                 :        238 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     859 [ +  + ][ +  + ]:        238 :   } while (degpol(u1)==0 || degpol(u1)==n);
     860                 :        210 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     861                 :        210 :   f1 = FpX_normalize(f1, p);
     862                 :        210 :   u2 = FpX_div(hp, u1, p);
     863                 :        210 :   f2 = FpX_div(Tp, f1, p);
     864         [ +  + ]:        210 :   if (degpol(u1)==1)
     865                 :        126 :     gel(V, idx) = f1;
     866                 :            :   else
     867                 :         84 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     868                 :        210 :   idx += degpol(f1)/d;
     869         [ +  + ]:        210 :   if (degpol(u2)==1)
     870                 :         98 :     gel(V, idx) = f2;
     871                 :            :   else
     872                 :        112 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     873                 :        210 : }
     874                 :            : 
     875                 :            : static void
     876                 :         14 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     877                 :            : {
     878                 :         14 :   long n = degpol(Tp), r = n/d, vT = varn(Tp);
     879                 :            :   GEN T, h, t;
     880                 :            :   pari_timer ti;
     881         [ -  + ]:         28 :   if (r==1) { gel(V, idx) = Tp; return; }
     882                 :         14 :   T = FpX_get_red(Tp, p);
     883                 :         14 :   XP = FpX_rem(XP, T, p);
     884         [ -  + ]:         14 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     885                 :            :   do
     886                 :            :   {
     887                 :         14 :     GEN g = random_FpX(n, vT, p);
     888                 :         14 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     889         [ -  + ]:         14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     890                 :         14 :     h = FpXQ_minpoly(t, T, p);
     891         [ -  + ]:         14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     892         [ -  + ]:         14 :   } while (degpol(h) != r);
     893                 :         14 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     894                 :            : }
     895                 :            : 
     896                 :            : static GEN
     897                 :         14 : FpX_factor_Shoup(GEN T, GEN p)
     898                 :            : {
     899                 :         14 :   long i, n, s = 0;
     900                 :            :   GEN XP, D, V;
     901                 :         14 :   long e = expi(p);
     902                 :            :   pari_timer ti;
     903                 :         14 :   n = get_FpX_degree(T);
     904                 :         14 :   T = FpX_get_red(T, p);
     905         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     906                 :         14 :   XP = FpX_Frobenius(T, p);
     907         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     908                 :         14 :   D = FpX_ddf(T, XP, p);
     909         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     910         [ +  + ]:        301 :   for (i = 1; i <= n; i++)
     911                 :        287 :     s += degpol(gel(D,i))/i;
     912                 :         14 :   V = cgetg(s+1, t_COL);
     913         [ +  + ]:        301 :   for (i = 1, s = 1; i <= n; i++)
     914                 :            :   {
     915                 :        287 :     GEN Di = gel(D,i);
     916                 :        287 :     long ni = degpol(Di), ri = ni/i;
     917         [ +  + ]:        287 :     if (ni == 0) continue;
     918                 :         21 :     Di = FpX_normalize(Di, p);
     919         [ +  + ]:         21 :     if (ni == i) { gel(V, s++) = Di; continue; }
     920         [ +  - ]:         14 :     if (ri <= e*expu(e))
     921                 :         14 :       FpX_edf(Di, XP, i, p, V, s);
     922                 :            :     else
     923                 :          0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     924         [ -  + ]:         14 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     925                 :         14 :     s += ri;
     926                 :            :   }
     927                 :         14 :   return V;
     928                 :            : }
     929                 :            : 
     930                 :            : static GEN
     931                 :          0 : FpX_simplefact_Shoup(GEN T, GEN p)
     932                 :            : {
     933                 :          0 :   long i, n, s = 0, j = 1, k;
     934                 :            :   GEN XP, D, V;
     935                 :            :   pari_timer ti;
     936                 :          0 :   n = get_FpX_degree(T);
     937                 :          0 :   T = FpX_get_red(T, p);
     938         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     939                 :          0 :   XP = FpX_Frobenius(T, p);
     940         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     941                 :          0 :   D = FpX_ddf(T, XP, p);
     942         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     943         [ #  # ]:          0 :   for (i = 1; i <= n; i++)
     944                 :          0 :     s += degpol(gel(D,i))/i;
     945                 :          0 :   V = cgetg(s+1, t_VEC);
     946         [ #  # ]:          0 :   for (i = 1; i <= n; i++)
     947                 :            :   {
     948                 :          0 :     long ni = degpol(gel(D,i)), ri = ni/i;
     949         [ #  # ]:          0 :     if (ni == 0) continue;
     950         [ #  # ]:          0 :     for (k = 1; k <= ri; k++)
     951                 :          0 :       gel(V, j++) = utoi(i);
     952                 :            :   }
     953                 :          0 :   return V;
     954                 :            : }
     955                 :            : 
     956                 :            : /* Yun algorithm: Assume p > degpol(T) */
     957                 :            : 
     958                 :            : static GEN
     959                 :          7 : FpX_factor_Yun(GEN T, GEN p)
     960                 :            : {
     961                 :          7 :   pari_sp av = avma;
     962                 :          7 :   long n = degpol(T);
     963                 :          7 :   long i = 1;
     964                 :          7 :   GEN d = FpX_deriv(T, p);
     965                 :            :   GEN a, b, c;
     966                 :          7 :   GEN V = cgetg(n+1,t_VEC);
     967                 :          7 :   a = FpX_gcd(T, d, p);
     968         [ -  + ]:          7 :   if (degpol(a) == 0) return mkvec(T);
     969                 :          7 :   b = FpX_div(T, a, p);
     970                 :            :   do
     971                 :            :   {
     972                 :         14 :     c = FpX_div(d, a, p);
     973                 :         14 :     d = FpX_sub(c, FpX_deriv(b, p), p);
     974                 :         14 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
     975                 :         14 :     gel(V, i++) = a;
     976                 :         14 :     b = FpX_div(b, a, p);
     977         [ +  + ]:         14 :   } while (degpol(b));
     978                 :          7 :   setlg(V, i);
     979                 :          7 :   return gerepilecopy(av, V);
     980                 :            : }
     981                 :            : 
     982                 :            : static GEN
     983                 :          7 : FpX_factor_Cantor(GEN T, GEN p)
     984                 :            : {
     985                 :            :   GEN E, F, M, V;
     986                 :            :   long i, j, s, l;
     987                 :          7 :   V = FpX_factor_Yun(T, p);
     988                 :          7 :   l = lg(V);
     989         [ +  + ]:         21 :   for (s=0, i=1; i < l; i++)
     990                 :         14 :     s += !!degpol(gel(V,i));
     991                 :          7 :   F = cgetg(s+1, t_VEC);
     992                 :          7 :   E = cgetg(s+1, t_VEC);
     993         [ +  + ]:         21 :   for (i=1, j=1; i < l; i++)
     994         [ +  - ]:         14 :     if (degpol(gel(V,i)))
     995                 :            :     {
     996                 :         14 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
     997                 :         14 :       gel(F, j) = Fj;
     998                 :         14 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
     999                 :         14 :       j++;
    1000                 :            :     }
    1001                 :          7 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1002                 :          7 :   return sort_factor_pol(M, cmpii);
    1003                 :            : }
    1004                 :            : 
    1005                 :            : static GEN
    1006                 :          0 : FpX_simplefact_Cantor(GEN T, GEN p)
    1007                 :            : {
    1008                 :            :   GEN E, F, M, V;
    1009                 :            :   long i, j, s, l;
    1010                 :          0 :   V = FpX_factor_Yun(get_FpX_mod(T), p);
    1011                 :          0 :   l = lg(V);
    1012         [ #  # ]:          0 :   for (s=0, i=1; i < l; i++)
    1013                 :          0 :     s += !!degpol(gel(V,i));
    1014                 :          0 :   F = cgetg(s+1, t_VEC);
    1015                 :          0 :   E = cgetg(s+1, t_VEC);
    1016         [ #  # ]:          0 :   for (i=1, j=1; i < l; i++)
    1017         [ #  # ]:          0 :     if (degpol(gel(V,i)))
    1018                 :            :     {
    1019                 :          0 :       GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
    1020                 :          0 :       gel(F, j) = Fj;
    1021                 :          0 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1022                 :          0 :       j++;
    1023                 :            :     }
    1024                 :          0 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1025                 :          0 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1026                 :            : }
    1027                 :            : 
    1028                 :            : static int
    1029                 :          0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1030                 :            : {
    1031                 :          0 :   pari_sp av = avma;
    1032                 :            :   pari_timer ti;
    1033                 :            :   long n, d;
    1034                 :          0 :   GEN T = get_FpX_mod(Tp);
    1035                 :          0 :   GEN dT = FpX_deriv(T, p);
    1036                 :            :   GEN XP, D;
    1037         [ #  # ]:          0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1038                 :          0 :   n = get_FpX_degree(T);
    1039                 :          0 :   T = FpX_get_red(Tp, p);
    1040         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1041                 :          0 :   XP = FpX_Frobenius(T, p);
    1042         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1043                 :          0 :   D = FpX_ddf(T, XP, p);
    1044         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1045                 :          0 :   d = degpol(gel(D, n));
    1046                 :          0 :   avma = av; return d==n;
    1047                 :            : }
    1048                 :            : 
    1049                 :            : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1050                 :            : 
    1051                 :            : /*Assume that p is large and odd*/
    1052                 :            : static GEN
    1053                 :         35 : FpX_factcantor_i(GEN f, GEN pp, long flag)
    1054                 :            : {
    1055                 :         35 :   long d = degpol(f);
    1056         [ +  + ]:         35 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1057      [ +  -  - ]:          7 :   switch(flag)
    1058                 :            :   {
    1059                 :          7 :     default: return FpX_factor_Cantor(f, pp);
    1060                 :          0 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1061         [ #  # ]:         35 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1062                 :            :   }
    1063                 :            : }
    1064                 :            : 
    1065                 :            : long
    1066                 :       4088 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1067                 :            : {
    1068                 :       4088 :   pari_sp av = avma;
    1069                 :       4088 :   GEN ddf = FpX_ddf(T, XP, p);
    1070                 :       4088 :   long l = lg(ddf), i, s=0;
    1071         [ +  + ]:      86716 :   for(i = 1; i < l; i++)
    1072                 :      82628 :     s += degpol(gel(ddf,i))/i;
    1073                 :       4088 :   avma = av; return s;
    1074                 :            : }
    1075                 :            : 
    1076                 :            : long
    1077                 :         28 : FpX_nbfact(GEN T, GEN p)
    1078                 :            : {
    1079                 :         28 :   pari_sp av = avma;
    1080                 :         28 :   GEN XP = FpX_Frobenius(T, p);
    1081                 :         28 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1082                 :         28 :   avma = av; return n;
    1083                 :            : }
    1084                 :            : 
    1085                 :            : /* INPUT:
    1086                 :            :  *  m integer (converted to polynomial w in Z[X] by stoFpX)
    1087                 :            :  *  p prime; q = (p^d-1) / 2^r
    1088                 :            :  *  t[0] polynomial of degree k*d product of k irreducible factors of degree d
    1089                 :            :  *  t[0] is expected to be normalized (leading coeff = 1)
    1090                 :            :  * OUTPUT:
    1091                 :            :  *  t[0],t[1]...t[k-1] the k factors, normalized */
    1092                 :            : 
    1093                 :            : static void
    1094                 :      21308 : F2x_split(ulong m, GEN *t, long d)
    1095                 :            : {
    1096                 :            :   long l, v, dv;
    1097                 :            :   pari_sp av0, av;
    1098                 :            :   GEN w;
    1099                 :            : 
    1100         [ +  + ]:      23653 :   dv = F2x_degree(*t); if (dv==d) return;
    1101                 :       2345 :   v=(*t)[1]; av0=avma;
    1102                 :       2345 :   for(av=avma;;avma=av)
    1103                 :            :   {
    1104                 :       3948 :     GEN w0 = w = F2xq_powu(polx_F2x(v), m-1, *t); m += 2;
    1105         [ +  + ]:      27006 :     for (l=1; l<d; l++) w = F2x_add(w0, F2xq_sqr(w, *t));
    1106                 :       3948 :     w = F2x_gcd(*t,w);
    1107 [ +  + ][ +  + ]:       3948 :     l = F2x_degree(w); if (l && l!=dv) break;
    1108                 :       1603 :   }
    1109                 :       2345 :   w = gerepileupto(av0, w);
    1110                 :       2345 :   l /= d; t[l]=F2x_div(*t,w); *t=w;
    1111                 :       2345 :   F2x_split(m,t+l,d);
    1112                 :       2345 :   F2x_split(m,t,  d);
    1113                 :            : }
    1114                 :            : 
    1115                 :            : /* p > 2 */
    1116                 :            : static GEN
    1117                 :          7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1118                 :            : {
    1119      [ -  -  + ]:          7 :   switch(d)
    1120                 :            :   {
    1121                 :            :     case -1:
    1122                 :          0 :     case 0: return NULL;
    1123                 :          0 :     case 1: return gen_1;
    1124                 :            :   }
    1125         [ -  + ]:          7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1126                 :            : }
    1127                 :            : /* p > 2 */
    1128                 :            : static GEN
    1129                 :         14 : FpX_degfact_2(GEN f, GEN p, long d)
    1130                 :            : {
    1131   [ -  -  -  + ]:         14 :   switch(d)
    1132                 :            :   {
    1133                 :          0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1134                 :          0 :     case 0: return trivial_fact();
    1135                 :          0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1136                 :            :   }
    1137      [ +  +  - ]:         14 :   switch(FpX_quad_factortype(f, p)) {
    1138                 :          7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1139                 :          7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1140                 :         14 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1141                 :            :   }
    1142                 :            : }
    1143                 :            : 
    1144                 :            : GEN
    1145                 :         70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1146                 :            : GEN
    1147                 :        742 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1148                 :            : 
    1149                 :            : /* Mod(0,p) * x, where x is f's main variable */
    1150                 :            : static GEN
    1151                 :         14 : Mod0pX(GEN f, GEN p)
    1152                 :         14 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
    1153                 :            : static GEN
    1154                 :         14 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
    1155                 :            : 
    1156                 :            : /* not gerepile safe */
    1157                 :            : static GEN
    1158                 :         74 : FpX_factor_2(GEN f, GEN p, long d)
    1159                 :            : {
    1160                 :            :   GEN r, s, R, S;
    1161                 :            :   long v;
    1162                 :            :   int sgn;
    1163   [ -  +  +  + ]:         74 :   switch(d)
    1164                 :            :   {
    1165                 :          0 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1166                 :          2 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1167                 :         10 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1168                 :            :   }
    1169                 :         62 :   r = FpX_quad_root(f, p, 1);
    1170         [ +  + ]:         62 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1171                 :         53 :   v = varn(f);
    1172                 :         53 :   s = FpX_otherroot(f, r, p);
    1173         [ +  - ]:         53 :   if (signe(r)) r = subii(p, r);
    1174         [ +  - ]:         53 :   if (signe(s)) s = subii(p, s);
    1175         [ +  + ]:         53 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1176                 :         53 :   R = deg1pol_shallow(gen_1, r, v);
    1177         [ +  + ]:         53 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1178                 :         46 :   S = deg1pol_shallow(gen_1, s, v);
    1179                 :         74 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1180                 :            : }
    1181                 :            : static GEN
    1182                 :         95 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1183                 :            : {
    1184      [ +  +  + ]:         95 :   switch(flag) {
    1185                 :          7 :     case 2: return FpX_is_irred_2(f, p, d);
    1186                 :         14 :     case 1: return FpX_degfact_2(f, p, d);
    1187                 :         95 :     default: return FpX_factor_2(f, p, d);
    1188                 :            :   }
    1189                 :            : }
    1190                 :            : 
    1191                 :            : static int
    1192                 :      27790 : F2x_quad_factortype(GEN x)
    1193         [ +  + ]:      27790 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1194                 :            : 
    1195                 :            : static GEN
    1196                 :          7 : F2x_is_irred_2(GEN f, long d)
    1197 [ +  - ][ +  - ]:          7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
                 [ +  - ]
    1198                 :            : 
    1199                 :            : static GEN
    1200                 :        847 : F2x_degfact_2(GEN f, long d)
    1201                 :            : {
    1202         [ -  + ]:        847 :   if (!d) return trivial_fact();
    1203         [ +  + ]:        847 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1204      [ +  +  + ]:        728 :   switch(F2x_quad_factortype(f)) {
    1205                 :        133 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1206                 :        154 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1207                 :        847 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1208                 :            :   }
    1209                 :            : }
    1210                 :            : 
    1211                 :            : static GEN
    1212                 :      46102 : F2x_factor_2(GEN f, long d)
    1213                 :            : {
    1214                 :      46102 :   long v = f[1];
    1215         [ -  + ]:      46102 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1216         [ +  + ]:      46102 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1217         [ +  + ]:      38577 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1218      [ +  +  + ]:      19796 :   switch(F2x_quad_factortype(f))
    1219                 :            :   {
    1220                 :       3822 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1221                 :      13447 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1222                 :      46102 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1223                 :            :   }
    1224                 :            : }
    1225                 :            : static GEN
    1226                 :      46956 : F2x_factor_deg2(GEN f, long d, long flag)
    1227                 :            : {
    1228      [ +  +  + ]:      46956 :   switch(flag) {
    1229                 :          7 :     case 2: return F2x_is_irred_2(f, d);
    1230                 :        847 :     case 1: return F2x_degfact_2(f, d);
    1231                 :      46956 :     default: return F2x_factor_2(f, d);
    1232                 :            :   }
    1233                 :            : }
    1234                 :            : 
    1235                 :            : static void
    1236                 :         19 : split_squares(struct split_t *S, GEN g, ulong p)
    1237                 :            : {
    1238                 :         19 :   ulong q = p >> 1;
    1239                 :         19 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1240         [ +  + ]:         19 :   if (degpol(a) < 0)
    1241                 :            :   {
    1242                 :            :     ulong i;
    1243                 :          6 :     split_add_done(S, (GEN)1);
    1244         [ +  + ]:         18 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1245                 :            :   } else {
    1246                 :         13 :     (void)Flx_valrem(a, &a);
    1247         [ +  - ]:         13 :     if (degpol(a))
    1248                 :            :     {
    1249                 :         13 :       a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
    1250         [ +  - ]:         13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1251                 :            :     }
    1252                 :            :   }
    1253                 :         19 : }
    1254                 :            : static void
    1255                 :         19 : split_nonsquares(struct split_t *S, GEN g, ulong p)
    1256                 :            : {
    1257                 :         19 :   ulong q = p >> 1;
    1258                 :         19 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1259         [ +  + ]:         19 :   if (degpol(a) < 0)
    1260                 :            :   {
    1261                 :          6 :     ulong i, z = Fl_nonsquare(p);
    1262                 :          6 :     split_add_done(S, (GEN)z);
    1263         [ +  + ]:         18 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1264                 :            :   } else {
    1265                 :         13 :     (void)Flx_valrem(a, &a);
    1266         [ +  - ]:         13 :     if (degpol(a))
    1267                 :            :     {
    1268                 :         13 :       a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
    1269         [ +  - ]:         13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1270                 :            :     }
    1271                 :            :   }
    1272                 :         19 : }
    1273                 :            : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1274                 :            :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1275                 :            : static int
    1276                 :     611259 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1277                 :            : {
    1278                 :     611259 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1279                 :     611257 :   long d = degpol(g);
    1280         [ +  + ]:     611258 :   if (d < 0) return 0;
    1281         [ -  + ]:     611223 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1282         [ +  + ]:     611225 :   if (!d) return 1;
    1283         [ +  + ]:     595678 :   if (p <= 1.4 * (ulong)d) {
    1284                 :            :     /* small p; split further using x^((p-1)/2) +/- 1.
    1285                 :            :      * 30% degree drop makes the extra gcd worth it. */
    1286                 :         19 :     split_squares(S, g, p);
    1287                 :         19 :     split_nonsquares(S, g, p);
    1288                 :            :   } else { /* large p; use x^(p-1) - 1 directly */
    1289                 :     595659 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1290         [ -  + ]:     595674 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1291                 :     595674 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1292                 :     595659 :     g = Flx_gcd(g,a, p);
    1293         [ +  - ]:     595652 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1294                 :            :   }
    1295                 :     611268 :   return 1;
    1296                 :            : }
    1297                 :            : 
    1298                 :            : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1299                 :            : static GEN
    1300                 :    5646784 : Flx_roots_i(GEN f, ulong p)
    1301                 :            : {
    1302                 :            :   GEN pol, g;
    1303                 :    5646784 :   long v = Flx_valrem(f, &g);
    1304                 :            :   ulong q;
    1305                 :            :   struct split_t S;
    1306                 :            : 
    1307                 :            :   /* optimization: test for small degree first */
    1308      [ +  +  + ]:    5646965 :   switch(degpol(g))
    1309                 :            :   {
    1310                 :            :     case 1: {
    1311                 :      23555 :       ulong r = p - g[2];
    1312         [ +  + ]:      23555 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1313                 :            :     }
    1314                 :            :     case 2: {
    1315                 :    5012224 :       ulong r = Flx_quad_root(g, p, 1), s;
    1316 [ +  + ][ +  + ]:    5006822 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1317                 :    5006800 :       s = Flx_otherroot(g,r, p);
    1318         [ +  + ]:    5011731 :       if (r < s)
    1319         [ +  + ]:    1251968 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1320         [ +  + ]:    3759763 :       else if (r > s)
    1321         [ +  + ]:    3759700 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1322                 :            :       else
    1323         [ -  + ]:         63 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1324                 :            :     }
    1325                 :            :   }
    1326                 :     611262 :   q = p >> 1;
    1327                 :     611262 :   split_init(&S, lg(f)-1);
    1328                 :     611257 :   settyp(S.done, t_VECSMALL);
    1329         [ +  + ]:     611257 :   if (v) split_add_done(&S, (GEN)0);
    1330         [ +  + ]:     611257 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1331                 :         35 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1332                 :     611231 :   pol = polx_Flx(f[1]);
    1333                 :     611237 :   for (pol[2]=1; ; pol[2]++)
    1334                 :            :   {
    1335                 :    1296133 :     long j, l = lg(S.todo);
    1336         [ +  + ]:    1296133 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1337 [ -  + ][ #  # ]:     684910 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1338         [ +  + ]:    1452472 :     for (j = 1; j < l; j++)
    1339                 :            :     {
    1340                 :     767576 :       GEN c = gel(S.todo,j);
    1341      [ +  +  + ]:     767576 :       switch(degpol(c))
    1342                 :            :       {
    1343                 :            :         case 1:
    1344                 :     590820 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1345                 :     590808 :           j--; l--; break;
    1346                 :            :         case 2: {
    1347                 :      79831 :           ulong r = Flx_quad_root(c, p, 0), s = Flx_otherroot(c,r, p);
    1348                 :      79832 :           split_done(&S, j, (GEN)r, (GEN)s);
    1349                 :      79830 :           j--; l--; break;
    1350                 :            :         }
    1351                 :            :         default: {
    1352                 :      96921 :           GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
    1353                 :            :           long db;
    1354                 :      96917 :           b = Flx_gcd(c,b, p); db = degpol(b);
    1355 [ +  + ][ +  + ]:      96923 :           if (db && db < degpol(c))
    1356                 :            :           {
    1357                 :      74965 :             b = Flx_normalize(b, p);
    1358                 :      74964 :             c = Flx_div(c,b, p);
    1359                 :      74962 :             split_todo(&S, j, b, c);
    1360                 :            :           }
    1361                 :            :         }
    1362                 :            :       }
    1363                 :            :     }
    1364                 :    6327042 :   }
    1365                 :            : }
    1366                 :            : 
    1367                 :            : GEN
    1368                 :    5596428 : Flx_roots(GEN f, ulong p)
    1369                 :            : {
    1370                 :    5596428 :   pari_sp av = avma;
    1371      [ -  +  + ]:    5596428 :   switch(lg(f))
    1372                 :            :   {
    1373                 :          0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1374                 :        592 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1375                 :            :   }
    1376         [ -  + ]:    5595836 :   if (p == 2) return Flx_root_mod_2(f);
    1377                 :    5595836 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1378                 :            : }
    1379                 :            : 
    1380                 :            : /* assume x reduced mod p, monic. */
    1381                 :            : static int
    1382                 :      74620 : Flx_quad_factortype(GEN x, ulong p)
    1383                 :            : {
    1384                 :      74620 :   ulong b = x[3], c = x[2];
    1385                 :      74620 :   return krouu(Fl_disc_bc(b, c, p), p);
    1386                 :            : }
    1387                 :            : static GEN
    1388                 :          0 : Flx_is_irred_2(GEN f, ulong p, long d)
    1389                 :            : {
    1390         [ #  # ]:          0 :   if (!d) return NULL;
    1391         [ #  # ]:          0 :   if (d == 1) return gen_1;
    1392         [ #  # ]:          0 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1393                 :            : }
    1394                 :            : static GEN
    1395                 :      76832 : Flx_degfact_2(GEN f, ulong p, long d)
    1396                 :            : {
    1397         [ -  + ]:      76832 :   if (!d) return trivial_fact();
    1398         [ +  + ]:      76832 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1399      [ +  +  + ]:      74620 :   switch(Flx_quad_factortype(f, p)) {
    1400                 :      34916 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1401                 :      38794 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1402                 :      76832 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1403                 :            :   }
    1404                 :            : }
    1405                 :            : /* p > 2 */
    1406                 :            : static GEN
    1407                 :     253401 : Flx_factor_2(GEN f, ulong p, long d)
    1408                 :            : {
    1409                 :            :   ulong r, s;
    1410                 :            :   GEN R,S;
    1411                 :     253401 :   long v = f[1];
    1412         [ -  + ]:     253401 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1413         [ +  + ]:     253401 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1414         [ +  + ]:     239711 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1415                 :     163911 :   r = Flx_quad_root(f, p, 1);
    1416         [ +  + ]:     163911 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1417                 :     104098 :   s = Flx_otherroot(f, r, p);
    1418                 :     104098 :   r = Fl_neg(r, p);
    1419                 :     104098 :   s = Fl_neg(s, p);
    1420         [ +  + ]:     104098 :   if (s < r) lswap(s,r);
    1421                 :     104098 :   R = mkvecsmall3(v,r,1);
    1422         [ +  + ]:     104098 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1423                 :      89783 :   S = mkvecsmall3(v,s,1);
    1424                 :     253401 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1425                 :            : }
    1426                 :            : static GEN
    1427                 :     330233 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1428                 :            : {
    1429      [ -  +  + ]:     330233 :   switch(flag) {
    1430                 :          0 :     case 2: return Flx_is_irred_2(f, p, d);
    1431                 :      76832 :     case 1: return Flx_degfact_2(f, p, d);
    1432                 :     330233 :     default: return Flx_factor_2(f, p, d);
    1433                 :            :   }
    1434                 :            : }
    1435                 :            : 
    1436                 :            : static GEN
    1437                 :      14168 : F2x_factcantor_i(GEN f, long flag)
    1438                 :            : {
    1439                 :      14168 :   long j, e, nbfact, d = F2x_degree(f);
    1440                 :            :   ulong k;
    1441                 :            :   GEN X, E, f2, g, g1, u, v, y, t;
    1442                 :            : 
    1443         [ +  + ]:      14168 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    1444                 :            :   /* to hold factors and exponents */
    1445         [ +  + ]:      11354 :   t = flag ? cgetg(d+1,t_VECSMALL): cgetg(d+1,t_VEC);
    1446                 :      11354 :   E = cgetg(d+1, t_VECSMALL);
    1447                 :      11354 :   X = polx_F2x(f[1]);
    1448                 :      11354 :   e = nbfact = 1;
    1449                 :            :   for(;;)
    1450                 :            :   {
    1451                 :      15428 :     f2 = F2x_gcd(f,F2x_deriv(f));
    1452 [ +  + ][ +  + ]:      15428 :     if (flag == 2 && F2x_degree(f2) > 0) return NULL;
    1453                 :      15386 :     g1 = F2x_div(f,f2);
    1454                 :      15386 :     k = 0;
    1455         [ +  + ]:      30562 :     while (F2x_degree(g1)>0)
    1456                 :            :     {
    1457                 :            :       pari_sp av;
    1458                 :            :       long du, dg, dg1;
    1459         [ +  + ]:      15267 :       k++; if (k%2==0) { k++; f2 = F2x_div(f2,g1); }
    1460                 :      15267 :       u = g1; g1 = F2x_gcd(f2,g1);
    1461                 :      15267 :       du = F2x_degree(u);
    1462                 :      15267 :       dg1 = F2x_degree(g1);
    1463         [ +  + ]:      15267 :       if (dg1>0)
    1464                 :            :       {
    1465                 :       1071 :         f2= F2x_div(f2,g1);
    1466         [ +  + ]:       1071 :         if (du == dg1) continue;
    1467                 :        616 :         u = F2x_div( u,g1);
    1468                 :        616 :         du -= dg1;
    1469                 :            :       }
    1470                 :            :       /* here u is square-free (product of irred. of multiplicity e * k) */
    1471                 :      14812 :       v = X;
    1472                 :      14812 :       av = avma;
    1473         [ +  + ]:      73150 :       for (d=1; d <= du>>1; d++)
    1474                 :            :       {
    1475                 :      58429 :         v = F2xq_sqr(v, u);
    1476                 :      58429 :         g = F2x_gcd(F2x_add(v, X), u);
    1477                 :      58429 :         dg = F2x_degree(g);
    1478         [ +  + ]:      58429 :         if (dg <= 0) {avma = (pari_sp)v; v = gerepileuptoleaf(av,v); continue;}
    1479                 :            :         /* g is a product of irred. pols, all of which have degree d */
    1480                 :      17143 :         j = nbfact+dg/d;
    1481         [ +  + ]:      17143 :         if (flag)
    1482                 :            :         {
    1483         [ +  + ]:        525 :           if (flag == 2) return NULL;
    1484         [ +  + ]:        987 :           for ( ; nbfact<j; nbfact++) { t[nbfact]=d; E[nbfact]=e*k; }
    1485                 :            :         }
    1486                 :            :         else
    1487                 :            :         {
    1488                 :      16618 :           gel(t,nbfact) = g;
    1489                 :      16618 :           F2x_split(2,&gel(t,nbfact),d);
    1490         [ +  + ]:      35581 :           for (; nbfact<j; nbfact++) E[nbfact]=e*k;
    1491                 :            :         }
    1492                 :      17052 :         du -= dg;
    1493                 :      17052 :         u = F2x_div(u,g);
    1494                 :      17052 :         v = F2x_rem(v,u);
    1495                 :      17052 :         av = avma;
    1496                 :            :       }
    1497         [ +  + ]:      14721 :       if (du)
    1498                 :            :       {
    1499         [ +  + ]:      13384 :         if (flag) t[nbfact] = du;
    1500                 :      11816 :         else  gel(t,nbfact) = u;
    1501                 :      13384 :         E[nbfact++]=e*k;
    1502                 :            :       }
    1503                 :            :     }
    1504         [ +  + ]:      15295 :     if (F2x_degree(f2)==0) break;
    1505                 :       4074 :     e <<=1; f = F2x_sqrt(f2);
    1506                 :       4074 :   }
    1507         [ +  + ]:      11221 :   if (flag == 2) return gen_1; /* irreducible */
    1508                 :      11214 :   setlg(t, nbfact);
    1509                 :      11214 :   setlg(E, nbfact); y = mkvec2(t, E);
    1510                 :      14168 :   return flag ? sort_factor(y, (void*)cmpGuGu, cmp_nodata)
    1511         [ +  + ]:      11214 :               : sort_factor_pol(y, cmpGuGu);
    1512                 :            : }
    1513                 :            : GEN
    1514                 :       7504 : F2x_factcantor(GEN f, long flag)
    1515                 :            : {
    1516                 :       7504 :   pari_sp av = avma;
    1517                 :       7504 :   GEN z = F2x_factcantor_i(f, flag);
    1518         [ -  + ]:       7504 :   if (flag == 2) { avma = av; return z; }
    1519                 :       7504 :   return gerepilecopy(av, z);
    1520                 :            : }
    1521                 :            : 
    1522                 :            : int
    1523                 :        140 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f,2); }
    1524                 :            : 
    1525                 :            : void
    1526                 :      14987 : F2xV_to_FlxV_inplace(GEN v)
    1527                 :            : {
    1528                 :            :   long i;
    1529         [ +  + ]:      32774 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1530                 :      14987 : }
    1531                 :            : void
    1532                 :     459761 : FlxV_to_ZXV_inplace(GEN v)
    1533                 :            : {
    1534                 :            :   long i;
    1535         [ +  + ]:    1230164 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1536                 :     459761 : }
    1537                 :            : void
    1538                 :      78510 : F2xV_to_ZXV_inplace(GEN v)
    1539                 :            : {
    1540                 :            :   long i;
    1541         [ +  + ]:     164389 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1542                 :      78510 : }
    1543                 :            : 
    1544                 :            : /* Adapted from Shoup NTL */
    1545                 :            : static GEN
    1546                 :     311872 : Flx_factor_squarefree(GEN f, ulong p)
    1547                 :            : {
    1548                 :     311872 :   pari_sp av = avma;
    1549                 :            :   GEN r, t, v, tv;
    1550                 :     311872 :   long q, n = degpol(f);
    1551                 :     311872 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1552                 :     311872 :   for(q = 1;;q *= p)
    1553                 :            :   {
    1554                 :     312250 :     r = Flx_gcd(f, Flx_deriv(f, p), p);
    1555                 :     312250 :     t = Flx_div(f, r, p);
    1556         [ +  + ]:     312250 :     if (degpol(t) > 0)
    1557                 :            :     {
    1558                 :            :       long j;
    1559                 :     312047 :       for(j = 1;;j++)
    1560                 :            :       {
    1561                 :     317983 :         v = Flx_gcd(r, t, p);
    1562                 :     317983 :         tv = Flx_div(t, v, p);
    1563         [ +  + ]:     317983 :         if (degpol(tv) > 0)
    1564                 :     315799 :           gel(u, j*q) = tv;
    1565         [ +  + ]:     317983 :         if (degpol(v) <= 0) break;
    1566                 :       5936 :         r = Flx_div(r, v, p);
    1567                 :       5936 :         t = v;
    1568                 :       5936 :       }
    1569         [ +  + ]:     312047 :       if (degpol(r) == 0) break;
    1570                 :            :     }
    1571                 :        378 :     f = Flx_deflate(r, p);
    1572                 :        378 :   }
    1573                 :     311872 :   return gerepilecopy(av, u);
    1574                 :            : }
    1575                 :            : 
    1576                 :            : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1577                 :            : static GEN
    1578                 :     424786 : Flx_ddf(GEN T, GEN XP, ulong p)
    1579                 :            : {
    1580                 :     424786 :   pari_sp av = avma;
    1581                 :            :   GEN b, g, h, F, f, Tr, xq;
    1582                 :            :   long i, j, n, v, bo, ro;
    1583                 :            :   long B, l, m;
    1584                 :            :   pari_timer ti;
    1585                 :     424786 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1586         [ -  + ]:     424786 :   if (n == 0) return cgetg(1, t_VEC);
    1587         [ +  + ]:     424786 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1588                 :     417786 :   B = n/2;
    1589                 :     417786 :   l = usqrt(B);
    1590                 :     417786 :   m = (B+l-1)/l;
    1591                 :     417786 :   T = Flx_get_red(T, p);
    1592                 :     417786 :   b = cgetg(l+2, t_VEC);
    1593                 :     417786 :   gel(b, 1) = polx_Flx(v);
    1594                 :     417786 :   gel(b, 2) = XP;
    1595                 :     417786 :   bo = brent_kung_optpow(n, l-1, 1);
    1596                 :     417786 :   ro = (bo-1) + (l-1)*((n-1)/bo);
    1597         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1598         [ +  + ]:     417786 :   if (expu(p) <= ro)
    1599         [ +  + ]:     177674 :     for (i = 3; i <= l+1; i++)
    1600                 :     109529 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1601                 :            :   else
    1602                 :            :   {
    1603                 :     349641 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1604         [ -  + ]:     349641 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
    1605         [ +  + ]:     389954 :     for (i = 3; i <= l+1; i++)
    1606                 :      40313 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1607                 :            :   }
    1608         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
    1609                 :     417786 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1610         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
    1611                 :     417786 :   g = cgetg(m+1, t_VEC);
    1612                 :     417786 :   gel(g, 1) = gel(xq, 2);
    1613         [ +  + ]:     925160 :   for(i = 2; i <= m; i++)
    1614                 :     507374 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1615         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
    1616                 :     417786 :   h = cgetg(m+1, t_VEC);
    1617         [ +  + ]:    1342946 :   for (j = 1; j <= m; j++)
    1618                 :            :   {
    1619                 :     925160 :     pari_sp av = avma;
    1620                 :     925160 :     GEN gj = gel(g, j);
    1621                 :     925160 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1622         [ +  + ]:    1509380 :     for (i = 2; i <= l; i++)
    1623                 :     584220 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1624                 :     925160 :     gel(h, j) = gerepileupto(av, e);
    1625                 :            :   }
    1626         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
    1627                 :     417786 :   Tr = get_Flx_mod(T);
    1628                 :     417786 :   F = cgetg(m+1, t_VEC);
    1629         [ +  + ]:    1342946 :   for (j = 1; j <= m; j++)
    1630                 :            :   {
    1631                 :     925160 :     gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
    1632                 :     925160 :     Tr = Flx_div(Tr, gel(F,j), p);
    1633                 :            :   }
    1634         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
    1635                 :     417786 :   f = const_vec(n, pol1_Flx(v));
    1636         [ +  + ]:    1342946 :   for (j = 1; j <= m; j++)
    1637                 :            :   {
    1638                 :     925160 :     GEN e = gel(F, j);
    1639         [ +  - ]:    1030622 :     for (i=l-1; i >= 0; i--)
    1640                 :            :     {
    1641                 :    1030622 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    1642         [ +  + ]:    1030622 :       if (degpol(u))
    1643                 :            :       {
    1644                 :     372695 :         gel(f, l*j-i) = u;
    1645                 :     372695 :         e = Flx_div(e, u, p);
    1646                 :            :       }
    1647         [ +  + ]:    1030622 :       if (!degpol(e)) break;
    1648                 :            :     }
    1649                 :            :   }
    1650         [ -  + ]:     417786 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
    1651         [ +  + ]:     417786 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    1652                 :     424786 :   return gerepilecopy(av, f);
    1653                 :            : }
    1654                 :            : 
    1655                 :            : static void
    1656                 :      45640 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1657                 :            : {
    1658                 :      45640 :   long n = degpol(Tp), r = n/d;
    1659                 :            :   GEN T, f, ff;
    1660                 :            :   ulong p2;
    1661         [ +  + ]:      66402 :   if (r==1) { gel(V, idx) = Tp; return; }
    1662                 :      20762 :   p2 = p>>1;
    1663                 :      20762 :   T = Flx_get_red(Tp, p);
    1664                 :      20762 :   XP = Flx_rem(XP, T, p);
    1665                 :            :   while (1)
    1666                 :            :   {
    1667                 :      22624 :     pari_sp btop = avma;
    1668                 :            :     long i;
    1669                 :      22624 :     GEN g = random_Flx(n, Tp[1], p);
    1670                 :      22624 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1671         [ +  + ]:      22624 :     if (lgpol(t) == 0) continue;
    1672         [ +  + ]:      46284 :     for(i=1; i<=10; i++)
    1673                 :            :     {
    1674                 :      44891 :       pari_sp btop2 = avma;
    1675                 :      44891 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    1676                 :      44891 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    1677 [ +  + ][ +  + ]:      44891 :       if (degpol(f) > 0 && degpol(f) < n) break;
    1678                 :      24129 :       avma = btop2;
    1679                 :            :     }
    1680 [ +  + ][ +  + ]:      22155 :     if (degpol(f) > 0 && degpol(f) < n) break;
    1681                 :       1393 :     avma = btop;
    1682                 :       1862 :   }
    1683                 :      20762 :   f = Flx_normalize(f, p);
    1684                 :      20762 :   ff = Flx_div(Tp, f ,p);
    1685                 :      20762 :   Flx_edf_simple(f, XP, d, p, V, idx);
    1686                 :      20762 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    1687                 :            : }
    1688                 :            : static void
    1689                 :            : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    1690                 :            : 
    1691                 :            : static void
    1692                 :      33243 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    1693                 :            : {
    1694                 :            :   pari_sp av;
    1695                 :      33243 :   GEN Tp = get_Flx_mod(T);
    1696                 :      33243 :   long n = degpol(hp), vT = Tp[1];
    1697                 :            :   GEN u1, u2, f1, f2;
    1698                 :      33243 :   ulong p2 = p>>1;
    1699                 :            :   GEN R, h;
    1700                 :      33243 :   h = Flx_get_red(hp, p);
    1701                 :      33243 :   t = Flx_rem(t, T, p);
    1702                 :      33243 :   av = avma;
    1703                 :            :   do
    1704                 :            :   {
    1705                 :      56217 :     avma = av;
    1706                 :      56217 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    1707                 :      56217 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    1708 [ +  + ][ +  + ]:      56217 :   } while (degpol(u1)==0 || degpol(u1)==n);
    1709                 :      33243 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    1710                 :      33243 :   f1 = Flx_normalize(f1, p);
    1711                 :      33243 :   u2 = Flx_div(hp, u1, p);
    1712                 :      33243 :   f2 = Flx_div(Tp, f1, p);
    1713         [ +  + ]:      33243 :   if (degpol(u1)==1)
    1714                 :            :   {
    1715         [ +  + ]:      26159 :     if (degpol(f1)==d)
    1716                 :      25697 :       gel(V, idx) = f1;
    1717                 :            :     else
    1718                 :        462 :       Flx_edf(f1, XP, d, p, V, idx);
    1719                 :            :   }
    1720                 :            :   else
    1721                 :       7084 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    1722                 :      33243 :   idx += degpol(f1)/d;
    1723         [ +  + ]:      33243 :   if (degpol(u2)==1)
    1724                 :            :   {
    1725         [ +  + ]:      25851 :     if (degpol(f2)==d)
    1726                 :      25445 :       gel(V, idx) = f2;
    1727                 :            :     else
    1728                 :        406 :       Flx_edf(f2, XP, d, p, V, idx);
    1729                 :            :   }
    1730                 :            :   else
    1731                 :       7392 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    1732                 :      33243 : }
    1733                 :            : 
    1734                 :            : static void
    1735                 :      18767 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1736                 :            : {
    1737                 :      18767 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    1738                 :            :   GEN T, h, t;
    1739                 :            :   pari_timer ti;
    1740         [ -  + ]:      37534 :   if (r==1) { gel(V, idx) = Tp; return; }
    1741                 :      18767 :   T = Flx_get_red(Tp, p);
    1742                 :      18767 :   XP = Flx_rem(XP, T, p);
    1743         [ -  + ]:      18767 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1744                 :            :   do
    1745                 :            :   {
    1746                 :      20594 :     GEN g = random_Flx(n, vT, p);
    1747                 :      20594 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1748         [ -  + ]:      20594 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    1749                 :      20594 :     h = Flxq_minpoly(t, T, p);
    1750         [ -  + ]:      20594 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    1751         [ +  + ]:      20594 :   } while (degpol(h) <= 1);
    1752                 :      18767 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    1753                 :            : }
    1754                 :            : 
    1755                 :            : static GEN
    1756                 :      31465 : Flx_factor_Shoup(GEN T, ulong p)
    1757                 :            : {
    1758                 :      31465 :   long i, n, s = 0;
    1759                 :            :   GEN XP, D, V;
    1760                 :      31465 :   long e = expu(p);
    1761                 :            :   pari_timer ti;
    1762                 :      31465 :   n = get_Flx_degree(T);
    1763                 :      31465 :   T = Flx_get_red(T, p);
    1764         [ -  + ]:      31465 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1765                 :      31465 :   XP = Flx_Frobenius(T, p);
    1766         [ -  + ]:      31465 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1767                 :      31465 :   D = Flx_ddf(T, XP, p);
    1768         [ -  + ]:      31465 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1769         [ +  + ]:     344561 :   for (i = 1; i <= n; i++)
    1770                 :     313096 :     s += degpol(gel(D,i))/i;
    1771                 :      31465 :   V = cgetg(s+1, t_COL);
    1772         [ +  + ]:     344561 :   for (i = 1, s = 1; i <= n; i++)
    1773                 :            :   {
    1774                 :     313096 :     GEN Di = gel(D,i);
    1775                 :     313096 :     long ni = degpol(Di), ri = ni/i;
    1776         [ +  + ]:     313096 :     if (ni == 0) continue;
    1777                 :      47432 :     Di = Flx_normalize(Di, p);
    1778         [ +  + ]:      47432 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1779         [ +  + ]:      22015 :     if (ri <= e*expu(e))
    1780                 :      17899 :       Flx_edf(Di, XP, i, p, V, s);
    1781                 :            :     else
    1782                 :       4116 :       Flx_edf_simple(Di, XP, i, p, V, s);
    1783         [ -  + ]:      22015 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    1784                 :      22015 :     s += ri;
    1785                 :            :   }
    1786                 :      31465 :   return V;
    1787                 :            : }
    1788                 :            : 
    1789                 :            : static GEN
    1790                 :     284334 : Flx_simplefact_Shoup(GEN T, ulong p)
    1791                 :            : {
    1792                 :     284334 :   long i, n, s = 0, j = 1, k;
    1793                 :            :   GEN XP, D, V;
    1794                 :            :   pari_timer ti;
    1795                 :     284334 :   n = get_Flx_degree(T);
    1796                 :     284334 :   T = Flx_get_red(T, p);
    1797         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1798                 :     284334 :   XP = Flx_Frobenius(T, p);
    1799         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1800                 :     284334 :   D = Flx_ddf(T, XP, p);
    1801         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1802         [ +  + ]:    1802462 :   for (i = 1; i <= n; i++)
    1803                 :    1518128 :     s += degpol(gel(D,i))/i;
    1804                 :     284334 :   V = cgetg(s+1, t_VECSMALL);
    1805         [ +  + ]:    1802462 :   for (i = 1; i <= n; i++)
    1806                 :            :   {
    1807                 :    1518128 :     long ni = degpol(gel(D,i)), ri = ni/i;
    1808         [ +  + ]:    1518128 :     if (ni == 0) continue;
    1809         [ +  + ]:    1054621 :     for (k = 1; k <= ri; k++)
    1810                 :     689669 :       V[j++] = i;
    1811                 :            :   }
    1812                 :     284334 :   return V;
    1813                 :            : }
    1814                 :            : 
    1815                 :            : static GEN
    1816                 :      28427 : Flx_factor_Cantor(GEN T, ulong p)
    1817                 :            : {
    1818                 :            :   GEN E, F, M, V;
    1819                 :            :   long i, j, s, l;
    1820                 :      28427 :   V = Flx_factor_squarefree(T, p);
    1821                 :      28427 :   l = lg(V);
    1822         [ +  + ]:     375697 :   for (s=0, i=1; i < l; i++)
    1823                 :     347270 :     s += !!degpol(gel(V,i));
    1824                 :      28427 :   F = cgetg(s+1, t_VEC);
    1825                 :      28427 :   E = cgetg(s+1, t_VEC);
    1826         [ +  + ]:     375697 :   for (i=1, j=1; i < l; i++)
    1827         [ +  + ]:     347270 :     if (degpol(gel(V,i)))
    1828                 :            :     {
    1829                 :      31465 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    1830                 :      31465 :       gel(F, j) = Fj;
    1831                 :      31465 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1832                 :      31465 :       j++;
    1833                 :            :     }
    1834                 :      28427 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1835                 :      28427 :   return sort_factor_pol(M, cmpGuGu);
    1836                 :            : }
    1837                 :            : 
    1838                 :            : static GEN
    1839                 :     283445 : Flx_simplefact_Cantor(GEN T, ulong p)
    1840                 :            : {
    1841                 :            :   GEN E, F, M, V;
    1842                 :            :   long i, j, s, l;
    1843                 :     283445 :   V = Flx_factor_squarefree(get_Flx_mod(T), p);
    1844                 :     283445 :   l = lg(V);
    1845         [ +  + ]:    2088707 :   for (s=0, i=1; i < l; i++)
    1846                 :    1805262 :     s += !!degpol(gel(V,i));
    1847                 :     283445 :   F = cgetg(s+1, t_VEC);
    1848                 :     283445 :   E = cgetg(s+1, t_VEC);
    1849         [ +  + ]:    2088707 :   for (i=1, j=1; i < l; i++)
    1850         [ +  + ]:    1805262 :     if (degpol(gel(V,i)))
    1851                 :            :     {
    1852                 :     284334 :       GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
    1853                 :     284334 :       gel(F, j) = Fj;
    1854                 :     284334 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1855                 :     284334 :       j++;
    1856                 :            :     }
    1857                 :     283445 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1858                 :     283445 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1859                 :            : }
    1860                 :            : 
    1861                 :            : static int
    1862                 :        581 : Flx_isirred_Cantor(GEN Tp, ulong p)
    1863                 :            : {
    1864                 :        581 :   pari_sp av = avma;
    1865                 :            :   pari_timer ti;
    1866                 :            :   long n, d;
    1867                 :        581 :   GEN T = get_Flx_mod(Tp);
    1868                 :        581 :   GEN dT = Flx_deriv(T, p);
    1869                 :            :   GEN XP, D;
    1870         [ +  + ]:        581 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1871                 :        441 :   n = get_Flx_degree(T);
    1872                 :        441 :   T = Flx_get_red(Tp, p);
    1873         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1874                 :        441 :   XP = Flx_Frobenius(T, p);
    1875         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1876                 :        441 :   D = Flx_ddf(T, XP, p);
    1877         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1878                 :        441 :   d = degpol(gel(D, n));
    1879                 :        581 :   avma = av; return d==n;
    1880                 :            : }
    1881                 :            : 
    1882                 :            : static GEN
    1883                 :     397279 : Flx_factcantor_i(GEN f, ulong pp, long flag)
    1884                 :            : {
    1885                 :            :   long d;
    1886         [ +  + ]:     397279 :   if (pp==2) { /*We need to handle 2 specially */
    1887                 :       6503 :     GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    1888         [ +  + ]:       6503 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    1889                 :       6503 :     return F;
    1890                 :            :   }
    1891                 :     390776 :   d = degpol(f);
    1892         [ +  + ]:     390776 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    1893      [ +  +  + ]:     312453 :   switch(flag)
    1894                 :            :   {
    1895                 :      28427 :     default: return Flx_factor_Cantor(f, pp);
    1896                 :     283445 :     case 1: return Flx_simplefact_Cantor(f, pp);
    1897         [ +  + ]:     397279 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    1898                 :            :   }
    1899                 :            : }
    1900                 :            : 
    1901                 :            : GEN
    1902                 :       7266 : Flx_factcantor(GEN f, ulong p, long flag)
    1903                 :            : {
    1904                 :       7266 :   pari_sp av = avma;
    1905                 :       7266 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
    1906         [ -  + ]:       7266 :   if (flag == 2) { avma = av; return z; }
    1907                 :       7266 :   return gerepilecopy(av, z);
    1908                 :            : }
    1909                 :            : 
    1910                 :            : GEN
    1911                 :     362748 : Flx_degfact(GEN f, ulong p)
    1912                 :            : {
    1913                 :     362748 :   pari_sp av = avma;
    1914                 :     362748 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
    1915                 :     362748 :   return gerepilecopy(av, z);
    1916                 :            : }
    1917                 :            : 
    1918                 :            : /* T must be squarefree mod p*/
    1919                 :            : GEN
    1920                 :      55710 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    1921                 :            : {
    1922                 :            :   GEN XP, D;
    1923                 :            :   pari_timer ti;
    1924                 :      55710 :   long i, s, n = get_Flx_degree(T);
    1925                 :      55710 :   GEN V = const_vecsmall(n, 0);
    1926                 :      55710 :   pari_sp av = avma;
    1927                 :      55710 :   T = Flx_get_red(T, p);
    1928         [ -  + ]:      55710 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1929                 :      55710 :   XP = Flx_Frobenius(T, p);
    1930         [ -  + ]:      55710 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1931                 :      55710 :   D = Flx_ddf(T, XP, p);
    1932         [ -  + ]:      55710 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1933         [ +  + ]:     730998 :   for (i = 1, s = 0; i <= n; i++)
    1934                 :            :   {
    1935                 :     675288 :     V[i] = degpol(gel(D,i))/i;
    1936                 :     675288 :     s += V[i];
    1937                 :            :   }
    1938                 :      55710 :   *nb = s;
    1939                 :      55710 :   avma = av; return V;
    1940                 :            : }
    1941                 :            : 
    1942                 :            : long
    1943                 :      52836 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    1944                 :            : {
    1945                 :      52836 :   pari_sp av = avma;
    1946                 :      52836 :   GEN ddf = Flx_ddf(T, XP, p);
    1947                 :      52836 :   long l = lg(ddf), i, s=0;
    1948         [ +  + ]:     616672 :   for(i = 1; i < l; i++)
    1949                 :     563836 :     s += degpol(gel(ddf,i))/i;
    1950                 :      52836 :   avma = av; return s;
    1951                 :            : }
    1952                 :            : 
    1953                 :            : /* T must be squarefree mod p*/
    1954                 :            : long
    1955                 :      52836 : Flx_nbfact(GEN T, ulong p)
    1956                 :            : {
    1957                 :      52836 :   pari_sp av = avma;
    1958                 :      52836 :   GEN XP = Flx_Frobenius(T, p);
    1959                 :      52836 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    1960                 :      52836 :   avma = av; return n;
    1961                 :            : }
    1962                 :            : 
    1963                 :            : int
    1964                 :        581 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
    1965                 :            : 
    1966                 :            : /* factor f (FpX or Flx) mod pp.
    1967                 :            :  * flag = 1: return the degrees, not the factors
    1968                 :            :  * flag = 2: return NULL if f is not irreducible.
    1969                 :            :  * Not gerepile-safe */
    1970                 :            : static GEN
    1971                 :         77 : factcantor_i(GEN f, GEN pp, long flag)
    1972                 :            : {
    1973         [ +  + ]:         77 :   if (typ(f) == t_VECSMALL)
    1974                 :            :   { /* lgefint(pp) = 3 */
    1975                 :            :     GEN F;
    1976                 :         42 :     ulong p = pp[2];
    1977         [ +  + ]:         42 :     if (p==2) {
    1978                 :         21 :       F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    1979         [ +  + ]:         21 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    1980                 :            :     } else {
    1981                 :         21 :       F = Flx_factcantor_i(f,p,flag);
    1982         [ +  + ]:         21 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    1983                 :            :     }
    1984                 :         42 :     return F;
    1985                 :            :   }
    1986                 :         77 :   return FpX_factcantor_i(f, pp, flag);
    1987                 :            : }
    1988                 :            : GEN
    1989                 :          0 : FpX_factcantor(GEN f, GEN pp, long flag)
    1990                 :            : {
    1991                 :          0 :   pari_sp av = avma;
    1992                 :            :   GEN z;
    1993                 :          0 :   ZX_factmod_init(&f,pp);
    1994                 :          0 :   z = factcantor_i(f,pp,flag);
    1995         [ #  # ]:          0 :   if (flag == 2) { avma = av; return z; }
    1996                 :          0 :   return gerepilecopy(av, z);
    1997                 :            : }
    1998                 :            : 
    1999                 :            : static GEN
    2000                 :        154 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
    2001                 :            : {
    2002                 :        154 :   pari_sp av = avma;
    2003                 :            :   long j, nbfact;
    2004                 :            :   GEN z, y, t, E, u, v;
    2005                 :            : 
    2006                 :        154 :   factmod_init(&f, p);
    2007      [ +  +  + ]:        154 :   switch(lg(f))
    2008                 :            :   {
    2009                 :         14 :     case 3: avma = av; return trivial_fact();
    2010                 :         14 :     case 2: return gerepileupto(av, zero_fact_intmod(f, p));
    2011                 :            :   }
    2012                 :        126 :   z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
    2013                 :        126 :   y = cgetg(3, t_MAT); nbfact = lg(t);
    2014                 :        126 :   u = cgetg(nbfact,t_COL); gel(y,1) = u;
    2015                 :        126 :   v = cgetg(nbfact,t_COL); gel(y,2) = v;
    2016         [ +  + ]:        126 :   if (flag)
    2017         [ +  + ]:         84 :     for (j=1; j<nbfact; j++)
    2018                 :            :     {
    2019                 :         56 :       gel(u,j) = utoi(uel(t,j));
    2020                 :         56 :       gel(v,j) = utoi(uel(E,j));
    2021                 :            :     }
    2022                 :            :   else
    2023         [ +  + ]:       7686 :     for (j=1; j<nbfact; j++)
    2024                 :            :     {
    2025                 :       7588 :       gel(u,j) = FpX_to_mod(gel(t,j), p);
    2026                 :       7588 :       gel(v,j) = utoi(uel(E,j));
    2027                 :            :     }
    2028                 :        154 :   return gerepileupto(av, y);
    2029                 :            : }
    2030                 :            : GEN
    2031                 :         77 : factcantor0(GEN f, GEN p, long flag)
    2032                 :         77 : { return factmod_aux(f, p, &factcantor_i, flag); }
    2033                 :            : GEN
    2034                 :         77 : factmod(GEN f, GEN p)
    2035                 :         77 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
    2036                 :            : 
    2037                 :            : /* Use this function when you think f is reducible, and that there are lots of
    2038                 :            :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2039                 :            : int
    2040                 :         14 : FpX_is_irred(GEN f, GEN p) {
    2041                 :         14 :   ZX_factmod_init(&f,p);
    2042                 :         14 :   return !!factcantor_i(f,p,2);
    2043                 :            : }
    2044                 :            : GEN
    2045                 :          0 : FpX_degfact(GEN f, GEN p) {
    2046                 :          0 :   pari_sp av = avma;
    2047                 :            :   GEN z;
    2048                 :          0 :   ZX_factmod_init(&f,p);
    2049                 :          0 :   z = factcantor_i(f,p,1);
    2050                 :          0 :   return gerepilecopy(av, z);
    2051                 :            : }
    2052                 :            : GEN
    2053                 :         49 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
    2054                 :            : GEN
    2055                 :         28 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
    2056                 :            : 
    2057                 :            : /* set x <-- x + c*y mod p */
    2058                 :            : /* x is not required to be normalized.*/
    2059                 :            : static void
    2060                 :     600640 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2061                 :            : {
    2062                 :            :   long i, lx, ly;
    2063                 :     600640 :   ulong *x=(ulong *)gx;
    2064                 :     600640 :   ulong *y=(ulong *)gy;
    2065         [ +  + ]:    1181983 :   if (!c) return;
    2066                 :     581343 :   lx = lg(gx);
    2067                 :     581343 :   ly = lg(gy);
    2068         [ -  + ]:     581343 :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2069         [ +  + ]:     581343 :   if (SMALL_ULONG(p))
    2070         [ +  + ]:    4065064 :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2071                 :            :   else
    2072         [ +  + ]:      17095 :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2073                 :            : }
    2074                 :            : 
    2075                 :            : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    2076                 :            : /* assume x1 != 0 */
    2077                 :            : static GEN
    2078                 :     118928 : deg1_Flx(ulong x1, ulong x0, ulong sv)
    2079                 :            : {
    2080                 :     118928 :   return mkvecsmall3(sv, x0, x1);
    2081                 :            : }
    2082                 :            : 
    2083                 :            : long
    2084                 :      43832 : F2x_split_Berlekamp(GEN *t)
    2085                 :            : {
    2086                 :      43832 :   GEN u = *t, a, b, vker;
    2087                 :      43832 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    2088                 :            : 
    2089         [ +  + ]:      43832 :   if (du == 1) return 1;
    2090         [ +  + ]:      28084 :   if (du == 2)
    2091                 :            :   {
    2092         [ -  + ]:       7259 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    2093                 :            :     {
    2094                 :          0 :       t[0] = mkvecsmall2(sv, 2);
    2095                 :          0 :       t[1] = mkvecsmall2(sv, 3);
    2096                 :          0 :       return 2;
    2097                 :            :     }
    2098                 :       7259 :     return 1;
    2099                 :            :   }
    2100                 :            : 
    2101                 :      20825 :   vker = F2x_Berlekamp_ker(u);
    2102                 :      20825 :   lb = lgcols(vker);
    2103                 :      20825 :   d = lg(vker)-1;
    2104                 :      20825 :   ir = 0;
    2105                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2106         [ +  + ]:      26327 :   for (L=1; L<d; )
    2107                 :            :   {
    2108                 :            :     GEN pol;
    2109         [ +  + ]:       5502 :     if (d == 2)
    2110                 :       5250 :       pol = F2v_to_F2x(gel(vker,2), sv);
    2111                 :            :     else
    2112                 :            :     {
    2113                 :        252 :       GEN v = zero_zv(lb);
    2114                 :        252 :       v[1] = du;
    2115                 :        252 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    2116         [ +  + ]:        798 :       for (i=2; i<=d; i++)
    2117         [ +  + ]:        546 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    2118                 :        252 :       pol = F2v_to_F2x(v, sv);
    2119                 :            :     }
    2120 [ +  + ][ +  + ]:      11249 :     for (i=ir; i<L && L<d; i++)
    2121                 :            :     {
    2122                 :       5747 :       a = t[i]; la = F2x_degree(a);
    2123 [ +  + ][ +  + ]:       5747 :       if (la == 1) { set_irred(i); }
    2124         [ -  + ]:       5726 :       else if (la == 2)
    2125                 :            :       {
    2126         [ #  # ]:          0 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    2127                 :            :         {
    2128                 :          0 :           t[i] = mkvecsmall2(sv, 2);
    2129                 :          0 :           t[L] = mkvecsmall2(sv, 3); L++;
    2130                 :            :         }
    2131         [ #  # ]:          0 :         set_irred(i);
    2132                 :            :       }
    2133                 :            :       else
    2134                 :            :       {
    2135                 :       5726 :         pari_sp av = avma;
    2136                 :            :         long lb;
    2137                 :       5726 :         b = F2x_rem(pol, a);
    2138         [ +  + ]:       5726 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    2139                 :       5411 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    2140 [ +  - ][ +  - ]:       5411 :         if (lb && lb < la)
    2141                 :            :         {
    2142                 :       5411 :           t[L] = F2x_div(a,b);
    2143                 :       5411 :           t[i]= b; L++;
    2144                 :            :         }
    2145                 :          0 :         else avma = av;
    2146                 :            :       }
    2147                 :            :     }
    2148                 :            :   }
    2149                 :      43832 :   return d;
    2150                 :            : }
    2151                 :            : 
    2152                 :            : /* p != 2 */
    2153                 :            : long
    2154                 :     232760 : Flx_split_Berlekamp(GEN *t, ulong p)
    2155                 :            : {
    2156                 :     232760 :   GEN u = *t, a,b,vker;
    2157                 :     232760 :   long d, i, ir, L, la, lb, sv = u[1];
    2158                 :     232760 :   long l = lg(u);
    2159                 :            :   ulong po2;
    2160                 :            : 
    2161         [ -  + ]:     232760 :   if (p == 2)
    2162                 :            :   {
    2163                 :          0 :     *t = Flx_to_F2x(*t);
    2164                 :          0 :     d = F2x_split_Berlekamp(t);
    2165         [ #  # ]:          0 :     for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
    2166                 :          0 :     return d;
    2167                 :            :   }
    2168                 :     232760 :   la = degpol(u);
    2169         [ +  + ]:     232760 :   if (la == 1) return 1;
    2170         [ +  + ]:     217459 :   if (la == 2)
    2171                 :            :   {
    2172                 :      11933 :     ulong r = Flx_quad_root(u,p,1);
    2173         [ +  + ]:      11933 :     if (r != p)
    2174                 :            :     {
    2175                 :       5439 :       t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
    2176                 :       5439 :       t[1] = deg1_Flx(1, p - r, sv);
    2177                 :       5439 :       return 2;
    2178                 :            :     }
    2179                 :       6494 :     return 1;
    2180                 :            :   }
    2181                 :            : 
    2182                 :     205526 :   vker = Flx_Berlekamp_ker(u,p);
    2183                 :     205526 :   vker = Flm_to_FlxV(vker, sv);
    2184                 :     205526 :   d = lg(vker)-1;
    2185                 :     205526 :   po2 = p >> 1; /* (p-1) / 2 */
    2186                 :     205526 :   ir = 0;
    2187                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2188         [ +  + ]:     507738 :   for (L=1; L<d; )
    2189                 :            :   {
    2190                 :     302212 :     GEN pol = zero_zv(l-2);
    2191                 :     302212 :     pol[1] = sv;
    2192                 :     302212 :     pol[2] = random_Fl(p); /*Assume vker[1]=1*/
    2193         [ +  + ]:     902852 :     for (i=2; i<=d; i++)
    2194                 :     600640 :       Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
    2195                 :     302212 :     (void)Flx_renormalize(pol,l-1);
    2196                 :            : 
    2197 [ +  + ][ +  + ]:     728902 :     for (i=ir; i<L && L<d; i++)
    2198                 :            :     {
    2199                 :     426690 :       a = t[i]; la = degpol(a);
    2200 [ +  + ][ +  + ]:     426690 :       if (la == 1) { set_irred(i); }
    2201         [ +  + ]:     384758 :       else if (la == 2)
    2202                 :            :       {
    2203                 :      66100 :         ulong r = Flx_quad_root(a,p,1);
    2204         [ +  + ]:      66100 :         if (r != p)
    2205                 :            :         {
    2206                 :      54025 :           t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
    2207                 :      54025 :           t[L] = deg1_Flx(1, p - r, sv); L++;
    2208                 :            :         }
    2209         [ +  + ]:      66100 :         set_irred(i);
    2210                 :            :       }
    2211                 :            :       else
    2212                 :            :       {
    2213                 :     318658 :         pari_sp av = avma;
    2214                 :     318658 :         b = Flx_rem(pol, a, p);
    2215         [ +  + ]:     318658 :         if (degpol(b) <= 0) { avma=av; continue; }
    2216                 :     300481 :         b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
    2217                 :     300481 :         b = Flx_gcd(a,b, p); lb = degpol(b);
    2218 [ +  + ][ +  + ]:     300481 :         if (lb && lb < la)
    2219                 :            :         {
    2220                 :     173156 :           b = Flx_normalize(b, p);
    2221                 :     173156 :           t[L] = Flx_div(a,b,p);
    2222                 :     173156 :           t[i]= b; L++;
    2223                 :            :         }
    2224                 :     127325 :         else avma = av;
    2225                 :            :       }
    2226                 :            :     }
    2227                 :            :   }
    2228                 :     232760 :   return d;
    2229                 :            : }
    2230                 :            : 
    2231                 :            : long
    2232                 :        692 : FpX_split_Berlekamp(GEN *t, GEN p)
    2233                 :            : {
    2234                 :        692 :   GEN u = *t, a,b,po2,vker;
    2235                 :        692 :   long d, i, ir, L, la, lb, vu = varn(u);
    2236         [ +  + ]:        692 :   if (lgefint(p) == 3)
    2237                 :            :   {
    2238                 :        630 :     ulong up = p[2];
    2239         [ -  + ]:        630 :     if (up == 2)
    2240                 :            :     {
    2241                 :          0 :       *t = ZX_to_F2x(*t);
    2242                 :          0 :       d = F2x_split_Berlekamp(t);
    2243         [ #  # ]:          0 :       for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
    2244                 :            :     }
    2245                 :            :     else
    2246                 :            :     {
    2247                 :        630 :       *t = ZX_to_Flx(*t, up);
    2248                 :        630 :       d = Flx_split_Berlekamp(t, up);
    2249         [ +  + ]:       3850 :       for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
    2250                 :            :     }
    2251                 :        630 :     return d;
    2252                 :            :   }
    2253                 :         62 :   la = degpol(u);
    2254         [ +  + ]:         62 :   if (la == 1) return 1;
    2255         [ +  + ]:         57 :   if (la == 2)
    2256                 :            :   {
    2257                 :          2 :     GEN r = FpX_quad_root(u,p,1);
    2258         [ -  + ]:          2 :     if (r)
    2259                 :            :     {
    2260                 :          0 :       t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
    2261                 :          0 :       t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
    2262                 :          0 :       return 2;
    2263                 :            :     }
    2264                 :          2 :     return 1;
    2265                 :            :   }
    2266                 :         55 :   vker = FpX_Berlekamp_ker(u,p);
    2267                 :         55 :   vker = RgM_to_RgXV(vker,vu);
    2268                 :         55 :   d = lg(vker)-1;
    2269                 :         55 :   po2 = shifti(p, -1); /* (p-1) / 2 */
    2270                 :         55 :   ir = 0;
    2271                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2272         [ +  + ]:        107 :   for (L=1; L<d; )
    2273                 :            :   {
    2274                 :         52 :     GEN pol = scalar_ZX_shallow(randomi(p), vu);
    2275         [ +  + ]:        208 :     for (i=2; i<=d; i++)
    2276                 :        156 :       pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
    2277                 :         52 :     pol = FpX_red(pol,p);
    2278 [ +  + ][ +  + ]:        157 :     for (i=ir; i<L && L<d; i++)
    2279                 :            :     {
    2280                 :        105 :       a = t[i]; la = degpol(a);
    2281 [ +  + ][ +  - ]:        105 :       if (la == 1) { set_irred(i); }
    2282         [ +  + ]:         98 :       else if (la == 2)
    2283                 :            :       {
    2284                 :         14 :         GEN r = FpX_quad_root(a,p,1);
    2285         [ +  - ]:         14 :         if (r)
    2286                 :            :         {
    2287                 :         14 :           t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
    2288                 :         14 :           t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
    2289                 :            :         }
    2290         [ +  + ]:         14 :         set_irred(i);
    2291                 :            :       }
    2292                 :            :       else
    2293                 :            :       {
    2294                 :         84 :         pari_sp av = avma;
    2295                 :         84 :         b = FpX_rem(pol, a, p);
    2296         [ +  + ]:         84 :         if (degpol(b) <= 0) { avma=av; continue; }
    2297                 :         63 :         b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
    2298                 :         63 :         b = FpX_gcd(a,b, p); lb = degpol(b);
    2299 [ +  + ][ +  + ]:         63 :         if (lb && lb < la)
    2300                 :            :         {
    2301                 :         43 :           b = FpX_normalize(b, p);
    2302                 :         43 :           t[L] = FpX_div(a,b,p);
    2303                 :         43 :           t[i]= b; L++;
    2304                 :            :         }
    2305                 :         20 :         else avma = av;
    2306                 :            :       }
    2307                 :            :     }
    2308                 :            :   }
    2309                 :        692 :   return d;
    2310                 :            : }
    2311                 :            : 
    2312                 :            : long
    2313                 :       7616 : FqX_is_squarefree(GEN P, GEN T, GEN p)
    2314                 :            : {
    2315                 :       7616 :   pari_sp av = avma;
    2316                 :       7616 :   GEN z = FqX_gcd(P, FqX_deriv(P, T, p), T, p);
    2317                 :       7616 :   avma = av;
    2318                 :       7616 :   return degpol(z)==0;
    2319                 :            : }
    2320                 :            : 
    2321                 :            : static GEN
    2322                 :      89479 : F2x_Berlekamp_i(GEN f, long flag)
    2323                 :            : {
    2324                 :      89479 :   long e, nbfact, val, d = F2x_degree(f);
    2325                 :            :   ulong k, j;
    2326                 :            :   GEN y, E, f2, g1, t;
    2327                 :            : 
    2328         [ +  + ]:      89479 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    2329                 :            : 
    2330                 :            :   /* to hold factors and exponents */
    2331         [ -  + ]:      45337 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2332                 :      45337 :   E = cgetg(d+1,t_VECSMALL);
    2333                 :      45337 :   val = F2x_valrem(f, &f);
    2334                 :      45337 :   e = nbfact = 1;
    2335         [ +  + ]:      45337 :   if (val) {
    2336 [ -  + ][ #  # ]:       9343 :     if (flag == 2 && val > 1) return NULL;
    2337         [ -  + ]:       9343 :     if (flag == 1)
    2338                 :          0 :       t[1] = 1;
    2339                 :            :     else
    2340                 :       9343 :       gel(t,1) = polx_F2x(f[1]);
    2341                 :       9343 :     E[1] = val; nbfact++;
    2342                 :            :   }
    2343                 :            : 
    2344                 :            :   for(;;)
    2345                 :            :   {
    2346                 :      82108 :     f2 = F2x_gcd(f,F2x_deriv(f));
    2347 [ -  + ][ #  # ]:      82108 :     if (flag == 2 && F2x_degree(f2)) return NULL;
    2348         [ +  + ]:      82108 :     g1 = F2x_degree(f2)? F2x_div(f,f2): f; /* squarefree */
    2349                 :      82108 :     k = 0;
    2350         [ +  + ]:     137308 :     while (F2x_degree(g1)>0)
    2351                 :            :     {
    2352                 :            :       GEN u;
    2353         [ +  + ]:      55200 :       k++; if (k%2 == 0) { k++; f2 = F2x_div(f2,g1); }
    2354                 :      55200 :       u = g1; /* deg(u) > 0 */
    2355         [ +  + ]:      55200 :       if (!F2x_degree(f2)) g1 = pol1_F2x(0); /* only its degree (= 0) matters */
    2356                 :            :       else
    2357                 :            :       {
    2358                 :            :         long dg1;
    2359                 :      14301 :         g1 = F2x_gcd(f2,g1);
    2360                 :      14301 :         dg1 = F2x_degree(g1);
    2361         [ +  + ]:      14301 :         if (dg1)
    2362                 :            :         {
    2363                 :      12495 :           f2= F2x_div(f2,g1);
    2364         [ +  + ]:      12495 :           if (F2x_degree(u) == dg1) continue;
    2365                 :       1127 :           u = F2x_div( u,g1);
    2366                 :            :         }
    2367                 :            :       }
    2368                 :            :       /* u is square-free (product of irred. of multiplicity e * k) */
    2369                 :      43832 :       gel(t,nbfact) = u;
    2370                 :      43832 :       d = F2x_split_Berlekamp(&gel(t,nbfact));
    2371 [ -  + ][ #  # ]:      43832 :       if (flag == 2 && d != 1) return NULL;
    2372         [ -  + ]:      43832 :       if (flag == 1)
    2373         [ #  # ]:          0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = F2x_degree(gel(t,nbfact+j));
    2374         [ +  + ]:      93075 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2375                 :      43832 :       nbfact += d;
    2376                 :            :     }
    2377         [ +  + ]:      82108 :     j = F2x_degree(f2); if (!j) break;
    2378                 :      36771 :     e *= 2; f = F2x_deflate(f2, 2);
    2379                 :      36771 :   }
    2380         [ -  + ]:      45337 :   if (flag == 2) return gen_1; /* irreducible */
    2381                 :      45337 :   setlg(t, nbfact);
    2382                 :      45337 :   setlg(E, nbfact); y = mkvec2(t,E);
    2383                 :      89479 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2384         [ -  + ]:      45337 :               : sort_factor_pol(y, cmpGuGu);
    2385                 :            : }
    2386                 :            : 
    2387                 :            : static GEN
    2388                 :     491751 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
    2389                 :            : {
    2390                 :     491751 :   long e, nbfact, val, d = degpol(f);
    2391                 :            :   ulong k, j;
    2392                 :            :   GEN y, E, f2, g1, t;
    2393                 :            : 
    2394         [ +  + ]:     491751 :   if (p == 2)
    2395                 :            :   {
    2396                 :      10962 :     GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
    2397         [ +  - ]:      10962 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2398                 :      10962 :     return F;
    2399                 :            :   }
    2400         [ +  + ]:     480789 :   if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
    2401                 :            : 
    2402                 :            :   /* to hold factors and exponents */
    2403         [ -  + ]:     228879 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2404                 :     228879 :   E = cgetg(d+1,t_VECSMALL);
    2405                 :     228879 :   val = Flx_valrem(f, &f);
    2406                 :     228879 :   e = nbfact = 1;
    2407         [ +  + ]:     228879 :   if (val) {
    2408 [ -  + ][ #  # ]:       7505 :     if (flag == 2 && val > 1) return NULL;
    2409         [ -  + ]:       7505 :     if (flag == 1)
    2410                 :          0 :       t[1] = 1;
    2411                 :            :     else
    2412                 :       7505 :       gel(t,1) = polx_Flx(f[1]);
    2413                 :       7505 :     E[1] = val; nbfact++;
    2414                 :            :   }
    2415                 :            : 
    2416                 :            :   for(;;)
    2417                 :            :   {
    2418                 :     234500 :     f2 = Flx_gcd(f,Flx_deriv(f,p), p);
    2419 [ -  + ][ #  # ]:     234500 :     if (flag == 2 && degpol(f2)) return NULL;
    2420         [ +  + ]:     234500 :     g1 = degpol(f2)? Flx_div(f,f2,p): f; /* squarefree */
    2421                 :     234500 :     k = 0;
    2422         [ +  + ]:     507684 :     while (degpol(g1)>0)
    2423                 :            :     {
    2424                 :            :       GEN u;
    2425         [ +  + ]:     273184 :       k++; if (k%p == 0) { k++; f2 = Flx_div(f2,g1,p); }
    2426                 :     273184 :       u = g1; /* deg(u) > 0 */
    2427         [ +  + ]:     273184 :       if (!degpol(f2)) g1 = pol1_Flx(0); /* only its degree (= 0) matters */
    2428                 :            :       else
    2429                 :            :       {
    2430                 :      46838 :         g1 = Flx_gcd(f2,g1, p);
    2431         [ +  + ]:      46838 :         if (degpol(g1))
    2432                 :            :         {
    2433                 :      46691 :           f2= Flx_div(f2,g1,p);
    2434         [ +  + ]:      46691 :           if (lg(u) == lg(g1)) continue;
    2435                 :       5637 :           u = Flx_div( u,g1,p);
    2436                 :            :         }
    2437                 :            :       }
    2438                 :            :       /* u is square-free (product of irred. of multiplicity e * k) */
    2439                 :     232130 :       u = Flx_normalize(u,p);
    2440                 :     232130 :       gel(t,nbfact) = u;
    2441                 :     232130 :       d = Flx_split_Berlekamp(&gel(t,nbfact), p);
    2442 [ -  + ][ #  # ]:     232130 :       if (flag == 2 && d != 1) return NULL;
    2443         [ -  + ]:     232130 :       if (flag == 1)
    2444         [ #  # ]:          0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2445         [ +  + ]:     694290 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2446                 :     232130 :       nbfact += d;
    2447                 :            :     }
    2448         [ -  + ]:     234500 :     if (!p) break;
    2449         [ +  + ]:     234500 :     j = degpol(f2); if (!j) break;
    2450         [ -  + ]:       5621 :     if (j % p) pari_err_PRIME("factmod",utoi(p));
    2451                 :       5621 :     e *= p; f = Flx_deflate(f2, p);
    2452                 :       5621 :   }
    2453         [ -  + ]:     228879 :   if (flag == 2) return gen_1; /* irreducible */
    2454                 :     228879 :   setlg(t, nbfact);
    2455                 :     228879 :   setlg(E, nbfact); y = mkvec2(t,E);
    2456                 :     491751 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2457         [ -  + ]:     228879 :               : sort_factor_pol(y, cmpGuGu);
    2458                 :            : }
    2459                 :            : 
    2460                 :            : /* f an FpX or an Flx */
    2461                 :            : static GEN
    2462                 :     535988 : FpX_Berlekamp_i(GEN f, GEN pp, long flag)
    2463                 :            : {
    2464                 :     535988 :   long e, nbfact, val, d = degpol(f);
    2465                 :            :   ulong k, j;
    2466                 :            :   GEN y, E, f2, g1, t;
    2467                 :            : 
    2468         [ +  + ]:     535988 :   if (typ(f) == t_VECSMALL)
    2469                 :            :   {/* lgefint(pp) == 3 */
    2470                 :     535863 :     ulong p = pp[2];
    2471                 :            :     GEN F;
    2472         [ +  + ]:     535863 :     if (p == 2) {
    2473                 :      78503 :       F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
    2474         [ +  - ]:      78503 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2475                 :            :     } else {
    2476                 :     457360 :       F = Flx_Berlekamp_i(f, p, flag);
    2477         [ +  - ]:     457360 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2478                 :            :     }
    2479                 :     535863 :     return F;
    2480                 :            :   }
    2481                 :            :   /* p is large (and odd) */
    2482         [ +  + ]:        125 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    2483                 :            : 
    2484                 :            :   /* to hold factors and exponents */
    2485         [ -  + ]:         58 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2486                 :         58 :   E = cgetg(d+1,t_VECSMALL);
    2487                 :         58 :   val = ZX_valrem(f, &f);
    2488                 :         58 :   e = nbfact = 1;
    2489         [ +  + ]:         58 :   if (val) {
    2490 [ -  + ][ #  # ]:          3 :     if (flag == 2 && val > 1) return NULL;
    2491         [ -  + ]:          3 :     if (flag == 1)
    2492                 :          0 :       t[1] = 1;
    2493                 :            :     else
    2494                 :          3 :       gel(t,1) = pol_x(varn(f));
    2495                 :          3 :     E[1] = val; nbfact++;
    2496                 :            :   }
    2497                 :            : 
    2498                 :         58 :   f2 = FpX_gcd(f,ZX_deriv(f), pp);
    2499 [ -  + ][ #  # ]:         58 :   if (flag == 2 && degpol(f2)) return NULL;
    2500         [ +  + ]:         58 :   g1 = degpol(f2)? FpX_div(f,f2,pp): f; /* squarefree */
    2501                 :         58 :   k = 0;
    2502         [ +  + ]:        121 :   while (degpol(g1)>0)
    2503                 :            :   {
    2504                 :            :     GEN u;
    2505                 :         63 :     k++;
    2506                 :         63 :     u = g1; /* deg(u) > 0 */
    2507         [ +  + ]:         63 :     if (!degpol(f2)) g1 = pol_1(0); /* only its degree (= 0) matters */
    2508                 :            :     else
    2509                 :            :     {
    2510                 :          6 :       g1 = FpX_gcd(f2,g1, pp);
    2511         [ +  - ]:          6 :       if (degpol(g1))
    2512                 :            :       {
    2513                 :          6 :         f2= FpX_div(f2,g1,pp);
    2514         [ +  + ]:          6 :         if (lg(u) == lg(g1)) continue;
    2515                 :          5 :         u = FpX_div( u,g1,pp);
    2516                 :            :       }
    2517                 :            :     }
    2518                 :            :     /* u is square-free (product of irred. of multiplicity e * k) */
    2519                 :         62 :     u = FpX_normalize(u,pp);
    2520                 :         62 :     gel(t,nbfact) = u;
    2521                 :         62 :     d = FpX_split_Berlekamp(&gel(t,nbfact), pp);
    2522 [ -  + ][ #  # ]:         62 :     if (flag == 2 && d != 1) return NULL;
    2523         [ -  + ]:         62 :     if (flag == 1)
    2524         [ #  # ]:          0 :       for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2525         [ +  + ]:        181 :     for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2526                 :         62 :     nbfact += d;
    2527                 :            :   }
    2528         [ -  + ]:         58 :   if (flag == 2) return gen_1; /* irreducible */
    2529                 :         58 :   setlg(t, nbfact);
    2530                 :         58 :   setlg(E, nbfact); y = mkvec2(t,E);
    2531                 :     535988 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2532         [ -  + ]:         58 :               : sort_factor_pol(y, cmpii);
    2533                 :            : }
    2534                 :            : GEN
    2535                 :     535925 : FpX_factor(GEN f, GEN p)
    2536                 :            : {
    2537                 :     535925 :   pari_sp av = avma;
    2538                 :     535925 :   ZX_factmod_init(&f, p);
    2539                 :     535925 :   return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
    2540                 :            : }
    2541                 :            : GEN
    2542                 :      61054 : Flx_factor(GEN f, ulong p)
    2543                 :            : {
    2544                 :      61054 :   pari_sp av = avma;
    2545         [ +  + ]:      61054 :   GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
    2546                 :      61054 :   return gerepilecopy(av, F);
    2547                 :            : }
    2548                 :            : GEN
    2549                 :         14 : F2x_factor(GEN f)
    2550                 :            : {
    2551                 :         14 :   pari_sp av = avma;
    2552                 :         14 :   return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
    2553                 :            : }
    2554                 :            : 
    2555                 :            : GEN
    2556                 :         98 : factormod0(GEN f, GEN p, long flag)
    2557                 :            : {
    2558      [ +  +  - ]:         98 :   switch(flag)
    2559                 :            :   {
    2560                 :         70 :     case 0: return factmod(f,p);
    2561                 :         28 :     case 1: return simplefactmod(f,p);
    2562                 :          0 :     default: pari_err_FLAG("factormod");
    2563                 :            :   }
    2564                 :         98 :   return NULL; /* not reached */
    2565                 :            : }
    2566                 :            : 
    2567                 :            : /*******************************************************************/
    2568                 :            : /*                                                                 */
    2569                 :            : /*                     FACTORIZATION IN F_q                        */
    2570                 :            : /*                                                                 */
    2571                 :            : /*******************************************************************/
    2572                 :            : static GEN
    2573                 :       8771 : to_Fq(GEN x, GEN T, GEN p)
    2574                 :            : {
    2575                 :       8771 :   long i, lx, tx = typ(x);
    2576                 :            :   GEN y;
    2577                 :            : 
    2578         [ +  + ]:       8771 :   if (tx == t_INT)
    2579                 :       4382 :     y = mkintmod(x,p);
    2580                 :            :   else
    2581                 :            :   {
    2582         [ -  + ]:       4389 :     if (tx != t_POL) pari_err_TYPE("to_Fq",x);
    2583                 :       4389 :     lx = lg(x);
    2584                 :       4389 :     y = cgetg(lx,t_POL); y[1] = x[1];
    2585         [ +  + ]:     133518 :     for (i=2; i<lx; i++) gel(y,i) = mkintmod(gel(x,i), p);
    2586                 :            :   }
    2587                 :       8771 :   return mkpolmod(y, T);
    2588                 :            : }
    2589                 :            : 
    2590                 :            : static GEN
    2591                 :       4354 : to_Fq_pol(GEN x, GEN T, GEN p)
    2592                 :            : {
    2593                 :       4354 :   long i, lx = lg(x);
    2594         [ +  + ]:      13076 :   for (i=2; i<lx; i++) gel(x,i) = to_Fq(gel(x,i),T,p);
    2595                 :       4354 :   return x;
    2596                 :            : }
    2597                 :            : 
    2598                 :            : static GEN
    2599                 :        245 : to_Fq_fact(GEN P, GEN E, GEN T, GEN p, pari_sp av)
    2600                 :            : {
    2601                 :            :   GEN y, u, v;
    2602                 :        245 :   long j, l = lg(P), nbf = lg(P);
    2603                 :            : 
    2604                 :        245 :   u = cgetg(nbf,t_COL);
    2605                 :        245 :   v = cgetg(nbf,t_COL);
    2606         [ +  + ]:       4599 :   for (j=1; j<l; j++)
    2607                 :            :   {
    2608                 :       4354 :     gel(u,j) = simplify_shallow(gel(P,j)); /* may contain pols of degree 0 */
    2609                 :       4354 :     gel(v,j) = utoi(uel(E,j));
    2610                 :            :   }
    2611                 :        245 :   y = gerepilecopy(av, mkmat2(u, v));
    2612                 :        245 :   u = gel(y,1);
    2613                 :        245 :   p = icopy(p);
    2614                 :        245 :   T = FpX_to_mod(T, p);
    2615         [ +  + ]:       4599 :   for (j=1; j<nbf; j++) gel(u,j) = to_Fq_pol(gel(u,j), T,p);
    2616                 :        245 :   return y;
    2617                 :            : }
    2618                 :            : static GEN
    2619                 :         28 : to_FqC(GEN P, GEN T, GEN p, pari_sp av)
    2620                 :            : {
    2621                 :            :   GEN u;
    2622                 :         28 :   long j, l = lg(P), nbf = lg(P);
    2623                 :            : 
    2624                 :         28 :   u = cgetg(nbf,t_COL);
    2625         [ +  + ]:         77 :   for (j=1; j<l; j++)
    2626                 :         49 :     gel(u,j) = simplify_shallow(gel(P,j)); /* may contain pols of degree 0 */
    2627                 :         28 :   u = gerepilecopy(av, u);
    2628                 :         28 :   p = icopy(p);
    2629                 :         28 :   T = FpX_to_mod(T, p);
    2630         [ +  + ]:         77 :   for (j=1; j<nbf; j++) gel(u,j) = to_Fq(gel(u,j), T,p);
    2631                 :         28 :   return u;
    2632                 :            : }
    2633                 :            : 
    2634                 :            : GEN
    2635                 :       4802 : FlxqXQ_halfFrobenius(GEN a, GEN S, GEN T, ulong p)
    2636                 :            : {
    2637                 :       4802 :   long vT = get_Flx_var(T);
    2638                 :       4802 :   GEN xp = Flx_Frobenius(T, p);
    2639                 :       4802 :   GEN Xp = FlxqXQ_powu(polx_FlxX(get_FlxqX_var(S), vT), p, S, T, p);
    2640                 :       4802 :   GEN ap2 = FlxqXQ_powu(a,p>>1, S, T, p);
    2641                 :       4802 :   GEN V = FlxqXQV_autsum(mkvec3(xp, Xp, ap2), get_Flx_degree(T), S, T, p);
    2642                 :       4802 :   return gel(V,3);
    2643                 :            : }
    2644                 :            : 
    2645                 :            : GEN
    2646                 :       4682 : FpXQXQ_halfFrobenius(GEN a, GEN S, GEN T, GEN p)
    2647                 :            : {
    2648         [ +  + ]:       4682 :   if (lgefint(p)==3)
    2649                 :            :   {
    2650                 :       4531 :     ulong pp = p[2];
    2651                 :       4531 :     long v = get_FpX_var(T);
    2652                 :       4531 :     GEN Tp = ZXT_to_FlxT(T,pp), Sp = ZXXT_to_FlxXT(S, pp, v);
    2653                 :       4531 :     return FlxX_to_ZXX(FlxqXQ_halfFrobenius(ZXX_to_FlxX(a,pp,v),Sp,Tp,pp));
    2654                 :            :   }
    2655                 :            :   else
    2656                 :            :   {
    2657                 :            :     GEN xp, Xp, ap2, V;
    2658                 :        151 :     T = FpX_get_red(T, p);
    2659                 :        151 :     S = FpXQX_get_red(S, T, p);
    2660                 :        151 :     xp = FpX_Frobenius(T, p);
    2661                 :        151 :     Xp = FpXQXQ_pow(pol_x(get_FpXQX_var(S)), p, S, T, p);
    2662                 :        151 :     ap2 = FpXQXQ_pow(a,shifti(p,-1), S, T, p);
    2663                 :        151 :     V = FpXQXQV_autsum(mkvec3(xp,Xp,ap2), get_FpX_degree(T), S, T, p);
    2664                 :       4682 :     return gel(V,3);
    2665                 :            :   }
    2666                 :            : }
    2667                 :            : 
    2668                 :            : GEN
    2669                 :      14065 : FlxqX_Frobenius(GEN S, GEN T, ulong p)
    2670                 :            : {
    2671                 :      14065 :   pari_sp av = avma;
    2672                 :      14065 :   long n = get_Flx_degree(T), vT = get_Flx_var(T);
    2673                 :      14065 :   GEN X  = polx_FlxX(get_FlxqX_var(S), vT);
    2674                 :      14065 :   GEN xp = Flx_Frobenius(T, p);
    2675                 :      14065 :   GEN Xp = FlxqXQ_powu(X, p, S, T, p);
    2676                 :      14065 :   GEN Xq = gel(FlxqXQV_autpow(mkvec2(xp,Xp), n, S, T, p), 2);
    2677                 :      14065 :   return gerepilecopy(av, Xq);
    2678                 :            : }
    2679                 :            : 
    2680                 :            : GEN
    2681                 :        138 : FpXQX_Frobenius(GEN S, GEN T, GEN p)
    2682                 :            : {
    2683                 :        138 :   pari_sp av = avma;
    2684                 :        138 :   long n = get_FpX_degree(T);
    2685                 :        138 :   GEN X  = pol_x(get_FpXQX_var(S));
    2686                 :        138 :   GEN xp = FpX_Frobenius(T, p);
    2687                 :        138 :   GEN Xp = FpXQXQ_pow(X, p, S, T, p);
    2688                 :        138 :   GEN Xq = gel(FpXQXQV_autpow(mkvec2(xp,Xp), n, S, T, p), 2);
    2689                 :        138 :   return gerepilecopy(av, Xq);
    2690                 :            : }
    2691                 :            : 
    2692                 :            : static GEN
    2693                 :        728 : FqX_Frobenius_powers(GEN S, GEN T, GEN p)
    2694                 :            : {
    2695                 :        728 :   long N = get_FpXQX_degree(S);
    2696         [ +  - ]:        728 :   if (lgefint(p)==3)
    2697                 :            :   {
    2698                 :        728 :     ulong pp = p[2];
    2699                 :        728 :     GEN Tp = ZXT_to_FlxT(T, pp), Sp = ZXXT_to_FlxXT(S, pp, get_FpX_var(T));
    2700                 :        728 :     GEN Xq = FlxqX_Frobenius(Sp, Tp, pp);
    2701                 :        728 :     return FlxqXQ_powers(Xq, N-1, Sp, Tp, pp);
    2702                 :            :   } else
    2703                 :            :   {
    2704                 :          0 :     GEN Xq = FpXQX_Frobenius(S, T, p);
    2705                 :        728 :     return FpXQXQ_powers(Xq, N-1, S, T, p);
    2706                 :            :   }
    2707                 :            : }
    2708                 :            : 
    2709                 :            : static GEN
    2710                 :       3332 : FqX_Frobenius_eval(GEN x, GEN V, GEN S, GEN T, GEN p)
    2711                 :            : {
    2712         [ +  - ]:       3332 :   if (lgefint(p)==3)
    2713                 :            :   {
    2714                 :       3332 :     ulong pp = p[2];
    2715                 :       3332 :     long v = get_FpX_var(T);
    2716                 :       3332 :     GEN Tp = ZXT_to_FlxT(T, pp), Sp = ZXXT_to_FlxXT(S, pp, v);
    2717                 :       3332 :     GEN xp = ZXX_to_FlxX(x, pp, v);
    2718                 :       3332 :     return FlxX_to_ZXX(FlxqX_FlxqXQV_eval(xp, V, Sp, Tp, pp));
    2719                 :            :   }
    2720                 :            :   else
    2721                 :       3332 :     return FpXQX_FpXQXQV_eval(x, V, S, T, p);
    2722                 :            : }
    2723                 :            : 
    2724                 :            : static GEN
    2725                 :          8 : FpXQX_split_part(GEN f, GEN T, GEN p)
    2726                 :            : {
    2727                 :          8 :   long n = degpol(f);
    2728                 :          8 :   GEN z, X = pol_x(varn(f));
    2729         [ -  + ]:          8 :   if (n <= 1) return f;
    2730                 :          8 :   f = FpXQX_red(f, T, p);
    2731                 :          8 :   z = FpXQX_Frobenius(f, T, p);
    2732                 :          8 :   z = FpXX_sub(z, X , p);
    2733                 :          8 :   return FpXQX_gcd(z, f, T, p);
    2734                 :            : }
    2735                 :            : 
    2736                 :            : static GEN
    2737                 :      12893 : FlxqX_split_part(GEN f, GEN T, ulong p)
    2738                 :            : {
    2739                 :      12893 :   long n = degpol(f);
    2740                 :      12893 :   GEN z, Xq, X = polx_FlxX(varn(f),get_Flx_var(T));
    2741         [ -  + ]:      12893 :   if (n <= 1) return f;
    2742                 :      12893 :   f = FlxqX_red(f, T, p);
    2743                 :      12893 :   Xq = FlxqX_Frobenius(f, T, p);
    2744                 :      12893 :   z = FlxX_sub(Xq, X , p);
    2745                 :      12893 :   return FlxqX_gcd(z, f, T, p);
    2746                 :            : }
    2747                 :            : 
    2748                 :            : long
    2749                 :        301 : FpXQX_nbroots(GEN f, GEN T, GEN p)
    2750                 :            : {
    2751                 :        301 :   pari_sp av = avma;
    2752                 :            :   GEN z;
    2753         [ +  + ]:        301 :   if(lgefint(p)==3)
    2754                 :            :   {
    2755                 :        293 :     ulong pp=p[2];
    2756                 :        293 :     z = FlxqX_split_part(ZXX_to_FlxX(f,pp,varn(T)),ZXT_to_FlxT(T,pp),pp);
    2757                 :            :   }
    2758                 :            :   else
    2759                 :          8 :     z = FpXQX_split_part(f, T, p);
    2760                 :        301 :   avma = av; return degpol(z);
    2761                 :            : }
    2762                 :            : 
    2763                 :            : long
    2764                 :       1785 : FqX_nbroots(GEN f, GEN T, GEN p)
    2765         [ +  + ]:       1785 : { return T ? FpXQX_nbroots(f, T, p): FpX_nbroots(f, p); }
    2766                 :            : 
    2767                 :            : long
    2768                 :      12600 : FlxqX_nbroots(GEN f, GEN T, ulong p)
    2769                 :            : {
    2770                 :      12600 :   pari_sp av = avma;
    2771                 :      12600 :   GEN z = FlxqX_split_part(f, T, p);
    2772                 :      12600 :   avma = av; return degpol(z);
    2773                 :            : }
    2774                 :            : 
    2775                 :            : GEN
    2776                 :         48 : FlxqX_Berlekamp_ker(GEN S, GEN T, ulong p)
    2777                 :            : {
    2778                 :         48 :   pari_sp ltop=avma;
    2779                 :         48 :   long j, N = get_FlxqX_degree(S);
    2780                 :         48 :   GEN Xq = FlxqX_Frobenius(S,T,p);
    2781                 :         48 :   GEN Q  = FlxqXQ_matrix_pow(Xq,N,N,S,T,p);
    2782         [ +  + ]:        266 :   for (j=1; j<=N; j++)
    2783                 :        218 :     gcoeff(Q,j,j) = Flx_Fl_add(gcoeff(Q,j,j), p-1, p);
    2784                 :         48 :   return gerepileupto(ltop, FlxqM_ker(Q,T,p));
    2785                 :            : }
    2786                 :            : 
    2787                 :            : GEN
    2788                 :         56 : FpXQX_Berlekamp_ker(GEN S, GEN T, GEN p)
    2789                 :            : {
    2790         [ +  + ]:         56 :   if (lgefint(p)==3)
    2791                 :            :   {
    2792                 :         48 :     ulong pp=p[2];
    2793                 :         48 :     long v = get_FpX_var(T);
    2794                 :         48 :     GEN Tp = ZXT_to_FlxT(T,pp), Sp = ZXX_to_FlxX(S,pp,v);
    2795                 :         48 :     return FlxM_to_ZXM(FlxqX_Berlekamp_ker(Sp, Tp, pp));
    2796                 :            :   } else
    2797                 :            :   {
    2798                 :          8 :     pari_sp ltop=avma;
    2799                 :          8 :     long j,N = get_FpXQX_degree(S);
    2800                 :          8 :     GEN Xq = FpXQX_Frobenius(S,T,p);
    2801                 :          8 :     GEN Q  = FpXQXQ_matrix_pow(Xq,N,N,S,T,p);
    2802         [ +  + ]:         84 :     for (j=1; j<=N; j++)
    2803                 :         76 :       gcoeff(Q,j,j) = Fq_sub(gcoeff(Q,j,j), gen_1, T, p);
    2804                 :         56 :     return gerepileupto(ltop, FqM_ker(Q,T,p));
    2805                 :            :   }
    2806                 :            : }
    2807                 :            : 
    2808                 :            : long
    2809                 :          0 : FpXQX_nbfact(GEN u, GEN T, GEN p)
    2810                 :            : {
    2811                 :          0 :   pari_sp av = avma;
    2812                 :          0 :   GEN vker = FpXQX_Berlekamp_ker(u, T, p);
    2813                 :          0 :   avma = av; return lg(vker)-1;
    2814                 :            : }
    2815                 :            : 
    2816                 :            : long
    2817                 :          0 : FqX_nbfact(GEN u, GEN T, GEN p)
    2818                 :            : {
    2819         [ #  # ]:          0 :   return T ? FpX_nbfact(u, p): FpXQX_nbfact(u, T, p);
    2820                 :            : }
    2821                 :            : 
    2822                 :            : long
    2823                 :         56 : FqX_split_Berlekamp(GEN *t, GEN T, GEN p)
    2824                 :            : {
    2825                 :         56 :   GEN u = *t, a,b,vker,pol;
    2826                 :         56 :   long vu = varn(u), vT = varn(T), dT = degpol(T);
    2827                 :            :   long d, i, ir, L, la, lb;
    2828                 :         56 :   T = FpX_get_red(T, p);
    2829                 :         56 :   vker = FpXQX_Berlekamp_ker(u,T,p);
    2830                 :         56 :   vker = RgM_to_RgXV(vker,vu);
    2831                 :         56 :   d = lg(vker)-1;
    2832                 :         56 :   ir = 0;
    2833                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2834         [ +  + ]:        235 :   for (L=1; L<d; )
    2835                 :            :   {
    2836                 :        179 :     pol= scalarpol(random_FpX(dT,vT,p),vu);
    2837         [ +  + ]:        705 :     for (i=2; i<=d; i++)
    2838                 :        526 :       pol = FqX_add(pol, FqX_Fq_mul(gel(vker,i),
    2839                 :            :                                     random_FpX(dT,vT,p), T, p), T, p);
    2840                 :        179 :     pol = FpXQX_red(pol,T,p);
    2841 [ +  + ][ +  + ]:        526 :     for (i=ir; i<L && L<d; i++)
    2842                 :            :     {
    2843                 :        347 :       a = t[i]; la = degpol(a);
    2844 [ +  + ][ +  + ]:        347 :       if (la == 1) { set_irred(i); }
    2845                 :            :       else
    2846                 :            :       {
    2847                 :        305 :         pari_sp av = avma;
    2848                 :        305 :         b = FqX_rem(pol, a, T,p);
    2849         [ +  + ]:        305 :         if (degpol(b)<=0) { avma=av; continue; }
    2850                 :        242 :         b = FpXQXQ_halfFrobenius(b, a,T,p);
    2851         [ +  + ]:        242 :         if (degpol(b)<=0) { avma=av; continue; }
    2852                 :        133 :         gel(b,2) = Fq_sub(gel(b,2), gen_1,T,p);
    2853                 :        133 :         b = FqX_gcd(a,b, T,p); lb = degpol(b);
    2854 [ +  - ][ +  - ]:        133 :         if (lb && lb < la)
    2855                 :            :         {
    2856                 :        133 :           b = FpXQX_normalize(b, T,p);
    2857                 :        133 :           t[L] = FqX_div(a,b,T,p);
    2858                 :        133 :           t[i]= b; L++;
    2859                 :            :         }
    2860                 :          0 :         else avma = av;
    2861                 :            :       }
    2862                 :            :     }
    2863                 :            :   }
    2864                 :         56 :   return d;
    2865                 :            : }
    2866                 :            : 
    2867                 :            : /* split into r factors of degree d */
    2868                 :            : static void
    2869                 :       6118 : FqX_split(GEN *t, long d, GEN q, GEN S, GEN T, GEN p)
    2870                 :            : {
    2871                 :       6118 :   GEN u = *t;
    2872                 :       6118 :   long l, v, is2, cnt, dt = degpol(u), dT = degpol(T);
    2873                 :            :   pari_sp av;
    2874                 :            :   pari_timer ti;
    2875                 :            :   GEN w,w0;
    2876                 :            : 
    2877         [ +  + ]:      12236 :   if (dt == d) return;
    2878                 :       2947 :   v = varn(*t);
    2879         [ -  + ]:       2947 :   if (DEBUGLEVEL > 6) timer_start(&ti);
    2880                 :       2947 :   av = avma; is2 = equaliu(p, 2);
    2881                 :       2947 :   for(cnt = 1;;cnt++, avma = av)
    2882                 :            :   { /* splits *t with probability ~ 1 - 2^(1-r) */
    2883                 :       4410 :     w = w0 = random_FpXQX(dt,v, T,p);
    2884         [ -  + ]:       4410 :     if (degpol(w) <= 0) continue;
    2885         [ +  + ]:       6769 :     for (l=1; l<d; l++) /* sum_{0<i<d} w^(q^i), result in (F_q)^r */
    2886                 :       2359 :       w = RgX_add(w0, FqX_Frobenius_eval(w, S, u, T, p));
    2887                 :       4410 :     w = FpXQX_red(w, T,p);
    2888         [ +  + ]:       4410 :     if (is2)
    2889                 :            :     {
    2890                 :         42 :       w0 = w;
    2891         [ +  + ]:        210 :       for (l=1; l<dT; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */
    2892                 :            :       {
    2893                 :        168 :         w = FqX_rem(FqX_sqr(w,T,p), *t, T,p);
    2894                 :        168 :         w = FpXX_red(RgX_add(w0,w), p);
    2895                 :            :       }
    2896                 :            :     }
    2897                 :            :     else
    2898                 :            :     {
    2899                 :       4368 :       w = FpXQXQ_halfFrobenius(w, *t, T, p);
    2900                 :            :       /* w in {-1,0,1}^r */
    2901         [ +  + ]:       4368 :       if (degpol(w) <= 0) continue;
    2902                 :       2940 :       gel(w,2) = gadd(gel(w,2), gen_1);
    2903                 :            :     }
    2904                 :       2982 :     w = FqX_gcd(*t,w, T,p); l = degpol(w);
    2905 [ +  + ][ +  + ]:       2982 :     if (l && l != dt) break;
    2906                 :       1463 :   }
    2907                 :       2947 :   w = gerepileupto(av,FpXQX_normalize(w,T,p));
    2908         [ -  + ]:       2947 :   if (DEBUGLEVEL > 6)
    2909                 :          0 :     err_printf("[FqX_split] splitting time: %ld (%ld trials)\n",
    2910                 :            :                timer_delay(&ti),cnt);
    2911                 :       2947 :   l /= d; t[l] = FqX_div(*t,w, T,p); *t = w;
    2912                 :       2947 :   FqX_split(t+l,d,q,S,T,p);
    2913                 :       6118 :   FqX_split(t  ,d,q,S,T,p);
    2914                 :            : }
    2915                 :            : 
    2916                 :            : /*******************************************************************/
    2917                 :            : /*                                                                 */
    2918                 :            : /*                  FACTOR USING TRAGER'S TRICK                    */
    2919                 :            : /*                                                                 */
    2920                 :            : /*******************************************************************/
    2921                 :            : static GEN
    2922                 :        182 : FqX_frobinv_inplace(GEN F, GEN T, GEN p)
    2923                 :            : {
    2924         [ +  + ]:        182 :   if (T)
    2925                 :            :   {
    2926                 :        161 :     GEN frobinv = powiu(p, degpol(T)-1);
    2927                 :        161 :     long i, l = lg(F);
    2928         [ +  + ]:       2100 :     for (i=2; i<l; i++) gel(F,i) = Fq_pow(gel(F,i), frobinv, T,p);
    2929                 :            :   }
    2930                 :        182 :   return F;
    2931                 :            : }
    2932                 :            : static GEN
    2933                 :        119 : FqX_frob_deflate(GEN f, GEN T, GEN p)
    2934                 :        119 : { return FqX_frobinv_inplace(RgX_deflate(f, itos(p)), T, p); }
    2935                 :            : 
    2936                 :            : static long
    2937                 :        792 : isabsolutepol(GEN f)
    2938                 :            : {
    2939                 :        792 :   long i, l = lg(f);
    2940         [ +  + ]:       6651 :   for(i=2; i<l; i++)
    2941                 :            :   {
    2942                 :       6217 :     GEN c = gel(f,i);
    2943 [ +  + ][ +  + ]:       6217 :     if (typ(c) == t_POL && degpol(c) > 0) return 0;
    2944                 :            :   }
    2945                 :        792 :   return 1;
    2946                 :            : }
    2947                 :            : 
    2948                 :            : static void
    2949                 :        770 : add(GEN z, GEN g, long d) { vectrunc_append(z, mkvec2(utoipos(d), g)); }
    2950                 :            : /* return number of roots of u; assume deg u >= 0 */
    2951                 :            : long
    2952                 :        147 : FqX_split_deg1(GEN *pz, GEN u, GEN T, GEN p)
    2953                 :            : {
    2954                 :        147 :   long dg, N = degpol(u);
    2955                 :        147 :   GEN v, S, g, X, z = vectrunc_init(N+1);
    2956                 :            : 
    2957                 :        147 :   *pz = z;
    2958         [ -  + ]:        147 :   if (N == 0) return 0;
    2959         [ -  + ]:        147 :   if (N == 1) return 1;
    2960                 :        147 :   v = X = pol_x(varn(u));
    2961                 :        147 :   S = FqX_Frobenius_powers(u, T, p);
    2962                 :        147 :   vectrunc_append(z, S);
    2963                 :        147 :   v = FqX_Frobenius_eval(v, S, u, T, p);
    2964                 :        147 :   g = FqX_gcd(FpXX_sub(v,X,p),u, T,p);
    2965                 :        147 :   dg = degpol(g);
    2966         [ +  - ]:        147 :   if (dg > 0) add(z, FpXQX_normalize(g,T,p), dg);
    2967                 :        147 :   return dg;
    2968                 :            : }
    2969                 :            : 
    2970                 :            : /* return number of factors; z not properly initialized if deg(u) <= 1 */
    2971                 :            : long
    2972                 :        553 : FqX_split_by_degree(GEN *pz, GEN u, GEN T, GEN p)
    2973                 :            : {
    2974                 :        553 :   long nb = 0, d, dg, N = degpol(u);
    2975                 :        553 :   GEN v, S, g, X, z = vectrunc_init(N+1);
    2976                 :            : 
    2977                 :        553 :   *pz = z;
    2978         [ -  + ]:        553 :   if (N <= 1) return 1;
    2979                 :        553 :   v = X = pol_x(varn(u));
    2980                 :        553 :   S = FqX_Frobenius_powers(u, T, p);
    2981                 :        553 :   vectrunc_append(z, S);
    2982         [ +  + ]:       1330 :   for (d=1; d <= N>>1; d++)
    2983                 :            :   {
    2984                 :        777 :     v = FqX_Frobenius_eval(v, S, u, T, p);
    2985                 :        777 :     g = FqX_gcd(FpXX_sub(v,X,p),u, T,p);
    2986         [ +  + ]:        777 :     dg = degpol(g); if (dg <= 0) continue;
    2987                 :            :     /* all factors of g have degree d */
    2988                 :        623 :     add(z, FpXQX_normalize(g, T,p), dg / d); nb += dg / d;
    2989                 :        623 :     N -= dg;
    2990         [ +  + ]:        623 :     if (N)
    2991                 :            :     {
    2992                 :         70 :       u = FqX_div(u,g, T,p);
    2993                 :         70 :       v = FqX_rem(v,u, T,p);
    2994                 :            :     }
    2995                 :            :   }
    2996         [ -  + ]:        553 :   if (N) { add(z, FpXQX_normalize(u, T,p), 1); nb++; }
    2997                 :        553 :   return nb;
    2998                 :            : }
    2999                 :            : 
    3000                 :            : static GEN
    3001                 :        196 : FqX_split_equal(GEN L, GEN S, GEN T, GEN p)
    3002                 :            : {
    3003                 :        196 :   long n = itos(gel(L,1));
    3004                 :        196 :   GEN u = gel(L,2), z = cgetg(n + 1, t_COL);
    3005                 :        196 :   gel(z,1) = u;
    3006                 :        196 :   FqX_split((GEN*)(z+1), degpol(u) / n, powiu(p, degpol(T)), S, T, p);
    3007                 :        196 :   return z;
    3008                 :            : }
    3009                 :            : GEN
    3010                 :         63 : FqX_split_roots(GEN z, GEN T, GEN p, GEN pol)
    3011                 :            : {
    3012                 :         63 :   GEN S = gel(z,1), L = gel(z,2), rep = FqX_split_equal(L, S, T, p);
    3013         [ +  + ]:         63 :   if (pol) rep = shallowconcat(rep, FqX_div(pol, gel(L,2), T,p));
    3014                 :         63 :   return rep;
    3015                 :            : }
    3016                 :            : GEN
    3017                 :        119 : FqX_split_all(GEN z, GEN T, GEN p)
    3018                 :            : {
    3019                 :        119 :   GEN S = gel(z,1), rep = cgetg(1, t_VEC);
    3020                 :        119 :   long i, l = lg(z);
    3021         [ +  + ]:        252 :   for (i = 2; i < l; i++)
    3022                 :        133 :     rep = shallowconcat(rep, FqX_split_equal(gel(z,i), S, T, p));
    3023                 :        119 :   return rep;
    3024                 :            : }
    3025                 :            : 
    3026                 :            : /* not memory-clean, as Flx_factorff_i, returning only linear factors */
    3027                 :            : static GEN
    3028                 :      22288 : Flx_rootsff_i(GEN P, GEN T, ulong p)
    3029                 :            : {
    3030                 :      22288 :   GEN V, F = gel(Flx_factor(P,p), 1);
    3031                 :      22288 :   long i, lfact = 1, nmax = lgpol(P), n = lg(F), dT = get_Flx_degree(T);
    3032                 :            : 
    3033                 :      22288 :   V = cgetg(nmax,t_COL);
    3034         [ +  + ]:      47817 :   for(i=1;i<n;i++)
    3035                 :            :   {
    3036                 :      25529 :     GEN R, Fi = gel(F,i);
    3037                 :      25529 :     long di = degpol(Fi), j, r;
    3038         [ +  + ]:      25529 :     if (dT % di) continue;
    3039                 :      24010 :     R = Flx_factorff_irred(gel(F,i),T,p);
    3040                 :      24010 :     r = lg(R);
    3041         [ +  + ]:      51429 :     for (j=1; j<r; j++,lfact++)
    3042                 :      27419 :       gel(V,lfact) = Flx_neg(gmael(R,j, 2), p);
    3043                 :            :   }
    3044                 :      22288 :   setlg(V,lfact);
    3045                 :      22288 :   gen_sort_inplace(V, (void*) &cmp_Flx, &cmp_nodata, NULL);
    3046                 :      22288 :   return V;
    3047                 :            : }
    3048                 :            : GEN
    3049                 :          0 : Flx_rootsff(GEN P, GEN T, ulong p)
    3050                 :            : {
    3051                 :          0 :   pari_sp av = avma;
    3052                 :          0 :   return gerepilecopy(av, Flx_rootsff_i(P, T, p));
    3053                 :            : }
    3054                 :            : 
    3055                 :            : /* dummy implementation */
    3056                 :            : static GEN
    3057                 :      14091 : F2x_rootsff_i(GEN P, GEN T)
    3058                 :            : {
    3059                 :      14091 :   return FlxC_to_F2xC(Flx_rootsff_i(F2x_to_Flx(P), F2x_to_Flx(T), 2UL));
    3060                 :            : }
    3061                 :            : 
    3062                 :            : /* not memory-clean, as FpX_factorff_i, returning only linear factors */
    3063                 :            : static GEN
    3064                 :          7 : FpX_rootsff_i(GEN P, GEN T, GEN p)
    3065                 :            : {
    3066                 :            :   GEN V, F;
    3067                 :            :   long i, lfact, nmax, n, dT;
    3068         [ -  + ]:          7 :   if (lgefint(p)==3)
    3069                 :            :   {
    3070                 :          0 :     ulong pp = p[2];
    3071                 :          0 :     GEN V = Flx_rootsff_i(ZX_to_Flx(P,pp), ZXT_to_FlxT(T,pp), pp);
    3072                 :          0 :     return FlxC_to_ZXC(V);
    3073                 :            :   }
    3074                 :          7 :   F = gel(FpX_factor(P,p), 1);
    3075                 :          7 :   lfact = 1; nmax = lgpol(P); n = lg(F); dT = get_FpX_degree(T);
    3076                 :            : 
    3077                 :          7 :   V = cgetg(nmax,t_COL);
    3078         [ +  + ]:         14 :   for(i=1;i<n;i++)
    3079                 :            :   {
    3080                 :          7 :     GEN R, Fi = gel(F,i);
    3081                 :          7 :     long di = degpol(Fi), j, r;
    3082         [ -  + ]:          7 :     if (dT % di) continue;
    3083                 :          7 :     R = FpX_factorff_irred(gel(F,i),T,p);
    3084                 :          7 :     r = lg(R);
    3085         [ +  + ]:         14 :     for (j=1; j<r; j++,lfact++)
    3086                 :          7 :       gel(V,lfact) = Fq_to_FpXQ(Fq_neg(gmael(R,j, 2), T, p), T, p);
    3087                 :            :   }
    3088                 :          7 :   setlg(V,lfact);
    3089                 :          7 :   gen_sort_inplace(V, (void*) &cmp_RgX, &cmp_nodata, NULL);
    3090                 :          7 :   return V;
    3091                 :            : }
    3092                 :            : GEN
    3093                 :          0 : FpX_rootsff(GEN P, GEN T, GEN p)
    3094                 :            : {
    3095                 :          0 :   pari_sp av = avma;
    3096                 :          0 :   return gerepilecopy(av, FpX_rootsff_i(P, T, p));
    3097                 :            : }
    3098                 :            : 
    3099                 :            : static GEN
    3100                 :      10801 : F2xqX_quad_roots(GEN P, GEN T)
    3101                 :            : {
    3102                 :      10801 :   GEN b= gel(P,3), c = gel(P,2);
    3103         [ +  + ]:      10801 :   if (lgpol(b))
    3104                 :            :   {
    3105                 :       9814 :     GEN z, d = F2xq_div(c, F2xq_sqr(b,T),T);
    3106         [ +  + ]:       9814 :     if (F2xq_trace(d,T))
    3107                 :       1050 :       return cgetg(1, t_COL);
    3108                 :       8764 :     z = F2xq_mul(b, F2xq_Artin_Schreier(d, T), T);
    3109                 :       8764 :     return mkcol2(z, F2x_add(b, z));
    3110                 :            :   }
    3111                 :            :   else
    3112                 :      10801 :     return mkcol(F2xq_sqrt(c, T));
    3113                 :            : }
    3114                 :            : 
    3115                 :            : /* Assume p>2 and x monic */
    3116                 :            : static GEN
    3117                 :      10369 : FlxqX_quad_roots(GEN x, GEN T, ulong p)
    3118                 :            : {
    3119                 :      10369 :   GEN s, D, nb, b = gel(x,3), c = gel(x,2);
    3120                 :      10369 :   D = Flx_sub(Flxq_sqr(b,T,p), Flx_mulu(c,4,p), p);
    3121                 :      10369 :   nb = Flx_neg(b,p);
    3122         [ -  + ]:      10369 :   if (lgpol(D)==0)
    3123                 :          0 :     return mkcol(Flx_halve(nb, p));
    3124                 :      10369 :   s = Flxq_sqrt(D,T,p);
    3125         [ +  + ]:      10369 :   if (!s) return cgetg(1, t_COL);
    3126                 :       9970 :   s = Flx_halve(Flx_add(s,nb,p),p);
    3127                 :      10369 :   return mkcol2(s, Flx_sub(nb,s,p));
    3128                 :            : }
    3129                 :            : 
    3130                 :            : static GEN
    3131                 :        207 : FpXQX_quad_roots(GEN x, GEN T, GEN p)
    3132                 :            : {
    3133                 :        207 :   GEN s, D, nb, b = gel(x,3), c = gel(x,2);
    3134         [ +  + ]:        207 :   if (equaliu(p, 2))
    3135                 :            :   {
    3136                 :         49 :     GEN f2 = ZXX_to_F2xX(x, get_FpX_var(T));
    3137                 :         49 :     s = F2xqX_quad_roots(f2, ZX_to_F2x(get_FpX_mod(T)));
    3138                 :         49 :     return F2xC_to_ZXC(s);
    3139                 :            :   }
    3140                 :        158 :   D = Fq_sub(Fq_sqr(b,T,p), Fq_Fp_mul(c,utoi(4),T,p), T,p);
    3141                 :        158 :   nb = Fq_neg(b,T,p);
    3142         [ -  + ]:        158 :   if (signe(D)==0)
    3143                 :          0 :     return mkcol(Fq_to_FpXQ(Fq_halve(nb,T, p),T,p));
    3144                 :        158 :   s = Fq_sqrt(D,T,p);
    3145         [ +  + ]:        158 :   if (!s) return cgetg(1, t_COL);
    3146                 :        151 :   s = Fq_halve(Fq_add(s,nb,T,p),T, p);
    3147                 :        207 :   return mkcol2(Fq_to_FpXQ(s,T,p), Fq_to_FpXQ(Fq_sub(nb,s,T,p),T,p));
    3148                 :            : }
    3149                 :            : 
    3150                 :            : static GEN
    3151                 :       9051 : F2xqX_Frobenius_deflate(GEN S, GEN T)
    3152                 :            : {
    3153                 :       9051 :   GEN F = RgX_deflate(S, 2);
    3154                 :       9051 :   long i, l = lg(F);
    3155         [ +  + ]:      29764 :   for (i=2; i<l; i++)
    3156                 :      20713 :     gel(F,i) = F2xq_sqrt(gel(F,i), T);
    3157                 :       9051 :   return F;
    3158                 :            : }
    3159                 :            : 
    3160                 :            : static GEN
    3161                 :      14091 : F2xX_to_F2x(GEN x)
    3162                 :            : {
    3163                 :      14091 :   long l=nbits2lg(lgpol(x));
    3164                 :      14091 :   GEN z=cgetg(l,t_VECSMALL);
    3165                 :            :   long i,j,k;
    3166                 :      14091 :   z[1]=x[1];
    3167         [ +  + ]:      49329 :   for(i=2, k=1,j=BITS_IN_LONG;i<lg(x);i++,j++)
    3168                 :            :   {
    3169         [ +  + ]:      35238 :     if (j==BITS_IN_LONG)
    3170                 :            :     {
    3171                 :      14091 :       j=0; k++; z[k]=0;
    3172                 :            :     }
    3173         [ +  + ]:      35238 :     if (lgpol(gel(x,i)))
    3174                 :      25991 :       z[k]|=1UL<<j;
    3175                 :            :   }
    3176                 :      14091 :   return F2x_renormalize(z,l);
    3177                 :            : }
    3178                 :            : 
    3179                 :            : static GEN
    3180                 :     161147 : F2xqX_easyroots(GEN f, GEN T)
    3181                 :            : {
    3182         [ +  + ]:     161147 :   if (F2xY_degreex(f) <= 0) return F2x_rootsff_i(F2xX_to_F2x(f), T);
    3183         [ +  + ]:     147056 :   if (degpol(f)==1) return mkcol(constant_coeff(f));
    3184         [ +  + ]:     125328 :   if (degpol(f)==2) return F2xqX_quad_roots(f,T);
    3185                 :     161147 :   return NULL;
    3186                 :            : }
    3187                 :            : 
    3188                 :            : /* Adapted from Shoup NTL */
    3189                 :            : static GEN
    3190                 :      57029 : F2xqX_factor_squarefree(GEN f, GEN T)
    3191                 :            : {
    3192                 :      57029 :   pari_sp av = avma;
    3193                 :            :   GEN r, t, v, tv;
    3194                 :      57029 :   long q, n = degpol(f);
    3195                 :      57029 :   GEN u = const_vec(n+1, pol1_F2xX(varn(f),T[1]));
    3196                 :      57029 :   for(q = 1;;q *= 2)
    3197                 :            :   {
    3198                 :      66080 :     r = F2xqX_gcd(f, F2xX_deriv(f), T);
    3199                 :      66080 :     t = F2xqX_div(f, r, T);
    3200         [ +  + ]:      66080 :     if (degpol(t) > 0)
    3201                 :            :     {
    3202                 :            :       long j;
    3203                 :      65310 :       for(j = 1;;j++)
    3204                 :            :       {
    3205                 :      68894 :         v = F2xqX_gcd(r, t, T);
    3206                 :      68894 :         tv = F2xqX_div(t, v, T);
    3207         [ +  + ]:      68894 :         if (degpol(tv) > 0)
    3208                 :      66934 :           gel(u, j*q) = F2xqX_normalize(tv, T);
    3209         [ +  + ]:      68894 :         if (degpol(v) <= 0) break;
    3210                 :       3584 :         r = F2xqX_div(r, v, T);
    3211                 :       3584 :         t = v;
    3212                 :       3584 :       }
    3213         [ +  + ]:      65310 :       if (degpol(r) == 0) break;
    3214                 :            :     }
    3215                 :       9051 :     f = F2xqX_Frobenius_deflate(r, T);
    3216                 :       9051 :   }
    3217                 :      57029 :   return gerepilecopy(av, u);
    3218                 :            : }
    3219                 :            : 
    3220                 :            : static void
    3221                 :      34237 : F2xqX_roots_edf(GEN Sp, GEN xp, GEN Xp, GEN T, GEN V, long idx)
    3222                 :            : {
    3223                 :            :   pari_sp btop;
    3224                 :      34237 :   long n = degpol(Sp);
    3225                 :            :   GEN S, f, ff;
    3226                 :      34237 :   GEN R = F2xqX_easyroots(Sp, T);
    3227         [ +  + ]:      34237 :   if (R)
    3228                 :            :   {
    3229                 :      33012 :     long i, l = lg(R)-1;
    3230         [ +  + ]:      75978 :     for (i=0; i<l; i++)
    3231                 :      42966 :       gel(V, idx+i) = gel(R,1+i);
    3232                 :      34237 :     return;
    3233                 :            :   }
    3234                 :       1225 :   S = Sp;
    3235                 :       1225 :   Xp = F2xqX_rem(Xp, S, T);
    3236                 :       1225 :   btop = avma;
    3237                 :            :   while (1)
    3238                 :            :   {
    3239                 :       1624 :     GEN a = random_F2xqX(degpol(Sp), varn(Sp), T);
    3240                 :       1624 :     GEN R = gel(F2xqXQV_auttrace(mkvec3(xp, Xp, a), F2x_degree(T), S, T), 3);
    3241                 :       1624 :     f = F2xqX_gcd(R, Sp, T);
    3242 [ +  + ][ +  + ]:       1624 :     if (degpol(f) > 0 && degpol(f) < n) break;
    3243                 :        399 :     avma = btop;
    3244                 :        399 :   }
    3245                 :       1225 :   f = gerepileupto(btop, F2xqX_normalize(f, T));
    3246                 :       1225 :   ff = F2xqX_div(Sp, f, T);
    3247                 :       1225 :   F2xqX_roots_edf(f, xp, Xp, T, V, idx);
    3248                 :       1225 :   F2xqX_roots_edf(ff,xp, Xp, T, V, idx+degpol(f));
    3249                 :            : }
    3250                 :            : 
    3251                 :            : static GEN
    3252                 :      66934 : F2xqX_roots_ddf(GEN f, GEN xp, GEN T)
    3253                 :            : {
    3254                 :            :   GEN X, Xp, Xq, g, V;
    3255                 :      66934 :   long n, dT = F2x_degree(T);
    3256                 :      66934 :   GEN R = F2xqX_easyroots(f, T);
    3257         [ +  + ]:      66934 :   if (R) return R;
    3258                 :      56322 :   X  = pol_x(varn(f));
    3259                 :      56322 :   Xp = F2xqXQ_sqr(X, f, T);
    3260                 :      56322 :   Xq = gel(F2xqXQV_autpow(mkvec2(xp, Xp), dT, f, T), 2);
    3261                 :      56322 :   g = F2xqX_gcd(F2xX_add(Xq, X), f, T);
    3262                 :      56322 :   n = degpol(g);
    3263         [ +  + ]:      56322 :   if (n==0) return cgetg(1, t_COL);
    3264                 :      31787 :   g = F2xqX_normalize(g, T);
    3265                 :      31787 :   V = cgetg(n+1,t_COL);
    3266                 :      31787 :   F2xqX_roots_edf(g, xp, Xp, T, V, 1);
    3267                 :      66934 :   return V;
    3268                 :            : }
    3269                 :            : static GEN
    3270                 :      59983 : F2xqX_roots_i(GEN S, GEN T)
    3271                 :            : {
    3272                 :            :   GEN xp, F, M, V, R;
    3273                 :            :   long i, j, s, l;
    3274                 :      59983 :   S = F2xqX_red(S, T);
    3275         [ -  + ]:      59983 :   if (!signe(S)) pari_err_ROOTS0("F2xqX_roots");
    3276         [ +  + ]:      59983 :   if (degpol(S)==0) return cgetg(1, t_COL);
    3277                 :      59976 :   S = F2xqX_normalize(S, T);
    3278                 :      59976 :   R = F2xqX_easyroots(S, T);
    3279         [ +  + ]:      59976 :   if (R) return gen_sort(R, (void*) &cmp_Flx, &cmp_nodata);
    3280                 :      57029 :   xp = F2x_Frobenius(T);
    3281                 :      57029 :   V = F2xqX_factor_squarefree(S, T);
    3282                 :      57029 :   l = lg(V);
    3283         [ +  + ]:     483973 :   for (s=0, i=1; i < l; i++)
    3284                 :     426944 :     s += !!degpol(gel(V,i));
    3285                 :      57029 :   F = cgetg(s+1, t_VEC);
    3286         [ +  + ]:     483973 :   for (i=1, j=1; i < l; i++)
    3287         [ +  + ]:     426944 :     if (degpol(gel(V,i)))
    3288                 :      66934 :       gel(F, j++) = F2xqX_roots_ddf(gel(V,i), xp, T);
    3289                 :      57029 :   M = shallowconcat1(F);
    3290                 :      57029 :   gen_sort_inplace(M, (void*) &cmp_Flx, &cmp_nodata, NULL);
    3291                 :      59983 :   return M;
    3292                 :            : }
    3293                 :            : 
    3294                 :            : static GEN
    3295                 :       8197 : FlxX_to_Flx(GEN f)
    3296                 :            : {
    3297                 :       8197 :   long i, l = lg(f);
    3298                 :       8197 :   GEN V = cgetg(l, t_VECSMALL);
    3299                 :       8197 :   V[1] = ((ulong)f[1])&VARNBITS;
    3300         [ +  + ]:      30499 :   for(i=2; i<l; i++)
    3301         [ +  + ]:      22302 :     V[i] = lgpol(gel(f,i)) ? mael(f,i,2): 0L;
    3302                 :       8197 :   return V;
    3303                 :            : }
    3304                 :            : 
    3305                 :            : static GEN
    3306                 :     165110 : FlxqX_easyroots(GEN f, GEN T, ulong p)
    3307                 :            : {
    3308         [ +  + ]:     165110 :   if (RgXY_degreex(f) <= 0) return Flx_rootsff_i(FlxX_to_Flx(f), T, p);
    3309         [ +  + ]:     156913 :   if (degpol(f)==1) return mkcol(Flx_neg(constant_coeff(f), p));
    3310         [ +  + ]:     130994 :   if (degpol(f)==2) return FlxqX_quad_roots(f,T,p);
    3311                 :     165110 :   return NULL;
    3312                 :            : }
    3313                 :            : 
    3314                 :            : static GEN
    3315                 :        427 : FlxqX_invFrobenius(GEN xp, GEN T, ulong p)
    3316                 :            : {
    3317                 :        427 :   return Flxq_autpow(xp, get_Flx_degree(T)-1, T, p);
    3318                 :            : }
    3319                 :            : 
    3320                 :            : static GEN
    3321                 :        427 : FlxqX_Frobenius_deflate(GEN S, GEN ixp, GEN T, ulong p)
    3322                 :            : {
    3323                 :        427 :   GEN F = RgX_deflate(S, p);
    3324                 :        427 :   long i, l = lg(F);
    3325         [ -  + ]:        427 :   if (typ(ixp)==t_INT)
    3326         [ #  # ]:          0 :     for (i=2; i<l; i++)
    3327                 :          0 :       gel(F,i) = Flxq_pow(gel(F,i), ixp, T, p);
    3328                 :            :   else
    3329                 :            :   {
    3330                 :        427 :     long n = brent_kung_optpow(get_Flx_degree(T)-1, l-2, 1);
    3331                 :        427 :     GEN V = Flxq_powers(ixp, n, T, p);
    3332         [ +  + ]:       1330 :     for (i=2; i<l; i++)
    3333                 :        903 :       gel(F,i) = Flx_FlxqV_eval(gel(F,i), V, T, p);
    3334                 :            :   }
    3335                 :        427 :   return F;
    3336                 :            : }
    3337                 :            : 
    3338                 :            : /* Adapted from Shoup NTL */
    3339                 :            : static GEN
    3340                 :      58632 : FlxqX_factor_squarefree(GEN f, GEN xp, GEN T, ulong p)
    3341                 :            : {
    3342                 :      58632 :   pari_sp av = avma;
    3343                 :            :   GEN r, t, v, tv;
    3344                 :      58632 :   long q, n = degpol(f);
    3345                 :      58632 :   GEN u = const_vec(n+1, pol1_FlxX(varn(f),get_Flx_var(T)));
    3346                 :      58632 :   GEN ixp = NULL;
    3347                 :      58632 :   for(q = 1;;q *= p)
    3348                 :            :   {
    3349                 :      59059 :     r = FlxqX_gcd(f, FlxX_deriv(f, p), T, p);
    3350                 :      59059 :     t = FlxqX_div(f, r, T, p);
    3351         [ +  + ]:      59059 :     if (degpol(t) > 0)
    3352                 :            :     {
    3353                 :            :       long j;
    3354                 :      59017 :       for(j = 1;;j++)
    3355                 :            :       {
    3356                 :      63147 :         v = FlxqX_gcd(r, t, T, p);
    3357                 :      63147 :         tv = FlxqX_div(t, v, T, p);
    3358         [ +  + ]:      63147 :         if (degpol(tv) > 0)
    3359                 :      62979 :           gel(u, j*q) = FlxqX_normalize(tv, T, p);
    3360         [ +  + ]:      63147 :         if (degpol(v) <= 0) break;
    3361                 :       4130 :         r = FlxqX_div(r, v, T, p);
    3362                 :       4130 :         t = v;
    3363                 :       4130 :       }
    3364         [ +  + ]:      59017 :       if (degpol(r) == 0) break;
    3365                 :            :     }
    3366         [ +  - ]:        427 :     if (!ixp) ixp = FlxqX_invFrobenius(xp, T, p);
    3367                 :        427 :     f = FlxqX_Frobenius_deflate(r, ixp, T, p);
    3368                 :        427 :   }
    3369                 :      58632 :   return gerepilecopy(av, u);
    3370                 :            : }
    3371                 :            : 
    3372                 :            : static void
    3373                 :      41897 : FlxqX_roots_edf(GEN Sp, GEN xp, GEN Xp, GEN T, ulong p, GEN V, long idx)
    3374                 :            : {
    3375                 :            :   pari_sp btop;
    3376                 :      41897 :   long n = degpol(Sp);
    3377                 :            :   GEN S, f, ff;
    3378                 :      41897 :   long vT = get_Flx_var(T), dT = get_Flx_degree(T);
    3379                 :      41897 :   GEN R = FlxqX_easyroots(Sp, T, p);
    3380         [ +  + ]:      41897 :   if (R)
    3381                 :            :   {
    3382                 :      38137 :     long i, l = lg(R)-1;
    3383         [ +  + ]:      86766 :     for (i=0; i<l; i++)
    3384                 :      48629 :       gel(V, idx+i) = gel(R,1+i);
    3385                 :      41897 :     return;
    3386                 :            :   }
    3387                 :       3760 :   S = FlxqX_get_red(Sp, T, p);
    3388                 :       3760 :   Xp = FlxqX_rem(Xp, S, T, p);
    3389                 :       3760 :   btop = avma;
    3390                 :            :   while (1)
    3391                 :            :   {
    3392                 :       4775 :     GEN a = deg1pol(pol1_Flx(vT), random_Flx(dT, vT, p), varn(Sp));
    3393                 :       4775 :     GEN ap2 = FlxqXQ_powu(a, p>>1, S, T, p);
    3394                 :       4775 :     GEN R = gel(FlxqXQV_autsum(mkvec3(xp, Xp, ap2), get_Flx_degree(T), S, T, p), 3);
    3395                 :       4775 :     f = FlxqX_gcd(FlxX_Flx_add(R, Flx_neg(pol1_Flx(vT), p), p), Sp, T, p);
    3396 [ +  + ][ +  + ]:       4775 :     if (degpol(f) > 0 && degpol(f) < n) break;
    3397                 :       1015 :     avma = btop;
    3398                 :       1015 :   }
    3399                 :       3760 :   f = gerepileupto(btop, FlxqX_normalize(f, T, p));
    3400                 :       3760 :   ff = FlxqX_div(Sp, f, T, p);
    3401                 :       3760 :   FlxqX_roots_edf(f, xp, Xp, T, p, V, idx);
    3402                 :       3760 :   FlxqX_roots_edf(ff,xp, Xp, T, p, V, idx+degpol(f));
    3403                 :            : }
    3404                 :            : 
    3405                 :            : static GEN
    3406                 :      62979 : FlxqX_roots_ddf(GEN f, GEN xp, GEN T, ulong p)
    3407                 :            : {
    3408                 :            :   GEN X, Xp, Xq, g, V;
    3409                 :      62979 :   long n, dT = get_Flx_degree(T);
    3410                 :      62979 :   GEN R = FlxqX_easyroots(f, T, p);
    3411         [ +  + ]:      62979 :   if (R) return R;
    3412                 :      58233 :   X  = pol_x(varn(f));
    3413                 :      58233 :   Xp = FlxqXQ_powu(X, p, f, T, p);
    3414                 :      58233 :   Xq = gel(FlxqXQV_autpow(mkvec2(xp, Xp), dT, f, T, p), 2);
    3415                 :      58233 :   g = FlxqX_gcd(FlxX_sub(Xq, X, p), f, T, p);
    3416                 :      58233 :   n = degpol(g);
    3417         [ +  + ]:      58233 :   if (n==0) return cgetg(1, t_COL);
    3418                 :      34377 :   g = FlxqX_normalize(g, T, p);
    3419                 :      34377 :   V = cgetg(n+1,t_COL);
    3420                 :      34377 :   FlxqX_roots_edf(g, xp, Xp, T, p, V, 1);
    3421                 :      62979 :   return V;
    3422                 :            : }
    3423                 :            : 
    3424                 :            : /* do not handle p==2 */
    3425                 :            : static GEN
    3426                 :      60241 : FlxqX_roots_i(GEN S, GEN T, ulong p)
    3427                 :            : {
    3428                 :            :   GEN xp, F, M, V, R;
    3429                 :            :   long i, j, s, l;
    3430                 :      60241 :   S = FlxqX_red(S, T, p);
    3431         [ -  + ]:      60241 :   if (!signe(S)) pari_err_ROOTS0("FlxqX_roots");
    3432         [ +  + ]:      60241 :   if (degpol(S)==0) return cgetg(1, t_COL);
    3433                 :      60234 :   S = FlxqX_normalize(S, T, p);
    3434                 :      60234 :   R = FlxqX_easyroots(S, T, p);
    3435         [ +  + ]:      60234 :   if (R) return gen_sort(R, (void*) &cmp_Flx, &cmp_nodata);
    3436                 :      58632 :   xp = Flx_Frobenius(T, p);
    3437                 :      58632 :   V = FlxqX_factor_squarefree(S, xp, T, p);
    3438                 :      58632 :   l = lg(V);
    3439         [ +  + ]:     429562 :   for (s=0, i=1; i < l; i++)
    3440                 :     370930 :     s += !!degpol(gel(V,i));
    3441                 :      58632 :   F = cgetg(s+1, t_VEC);
    3442         [ +  + ]:     429562 :   for (i=1, j=1; i < l; i++)
    3443         [ +  + ]:     370930 :     if (degpol(gel(V,i)))
    3444                 :      62979 :       gel(F, j++) = FlxqX_roots_ddf(gel(V,i), xp, T, p);
    3445                 :      58632 :   M = shallowconcat1(F);
    3446                 :      58632 :   gen_sort_inplace(M, (void*) &cmp_Flx, &cmp_nodata, NULL);
    3447                 :      60241 :   return M;
    3448                 :            : }
    3449                 :            : 
    3450                 :            : static GEN
    3451                 :        169 : FpXQX_easyroots(GEN f, GEN T, GEN p)
    3452                 :            : {
    3453         [ +  + ]:        169 :   if (isabsolutepol(f)) return FpX_rootsff_i(simplify_shallow(f), T, p);
    3454         [ +  + ]:        162 :   if (degpol(f)==1) return mkcol(Fq_to_FpXQ(Fq_neg(constant_coeff(f),T,p),T,p));
    3455         [ +  + ]:        144 :   if (degpol(f)==2) return FpXQX_quad_roots(f,T,p);
    3456                 :        169 :   return NULL;
    3457                 :            : }
    3458                 :            : 
    3459                 :            : /* Adapted from Shoup NTL */
    3460                 :            : static GEN
    3461                 :          7 : FpXQX_factor_squarefree(GEN f, GEN T, GEN p)
    3462                 :            : {
    3463                 :          7 :   pari_sp av = avma;
    3464                 :            :   GEN r, t, v, tv;
    3465                 :          7 :   long j, n = degpol(f);
    3466                 :          7 :   GEN u = const_vec(n+1, pol_1(varn(f)));
    3467                 :          7 :   r = FpXQX_gcd(f, FpXX_deriv(f, p), T, p);
    3468                 :          7 :   t = FpXQX_div(f, r, T, p);
    3469                 :          7 :   for (j = 1;;j++)
    3470                 :            :   {
    3471                 :          7 :     v = FpXQX_gcd(r, t, T, p);
    3472                 :          7 :     tv = FpXQX_div(t, v, T, p);
    3473         [ +  - ]:          7 :     if (degpol(tv) > 0)
    3474                 :          7 :       gel(u, j) = FpXQX_normalize(tv, T, p);
    3475         [ +  - ]:          7 :     if (degpol(v) <= 0) break;
    3476                 :          0 :     r = FpXQX_div(r, v, T, p);
    3477                 :          0 :     t = v;
    3478                 :          0 :   }
    3479                 :          7 :   return gerepilecopy(av, u);
    3480                 :            : }
    3481                 :            : 
    3482                 :            : static void
    3483                 :         35 : FpXQX_roots_edf(GEN Sp, GEN xp, GEN Xp, GEN T, GEN p, GEN V, long idx)
    3484                 :            : {
    3485                 :            :   pari_sp btop;
    3486                 :         35 :   long n = degpol(Sp);
    3487                 :            :   GEN S, f, ff;
    3488                 :         35 :   long vT = get_FpX_var(T), dT = get_FpX_degree(T);
    3489                 :         35 :   GEN R = FpXQX_easyroots(Sp, T, p);
    3490         [ +  + ]:         35 :   if (R)
    3491                 :            :   {
    3492                 :         21 :     long i, l = lg(R)-1;
    3493         [ +  + ]:         49 :     for (i=0; i<l; i++)
    3494                 :         28 :       gel(V, idx+i) = gel(R,1+i);
    3495                 :         35 :     return;
    3496                 :            :   }
    3497                 :         14 :   S = FpXQX_get_red(Sp, T, p);
    3498                 :         14 :   Xp = FpXQX_rem(Xp, S, T, p);
    3499                 :         14 :   btop = avma;
    3500                 :            :   while (1)
    3501                 :            :   {
    3502                 :         14 :     GEN a = deg1pol(pol_1(vT), random_FpX(dT, vT, p), varn(Sp));
    3503                 :         14 :     GEN ap2 = FpXQXQ_pow(a, shifti(p,-1), S, T, p);
    3504                 :         14 :     GEN R = gel(FpXQXQV_autsum(mkvec3(xp, Xp, ap2), get_FpX_degree(T), S, T, p), 3);
    3505                 :         14 :     f = FpXQX_gcd(FqX_Fq_add(R, FpX_neg(pol_1(vT), p), T, p), Sp, T, p);
    3506 [ +  - ][ +  - ]:         14 :     if (degpol(f) > 0 && degpol(f) < n) break;
    3507                 :          0 :     avma = btop;
    3508                 :          0 :   }
    3509                 :         14 :   f = gerepileupto(btop, FpXQX_normalize(f, T, p));
    3510                 :         14 :   ff = FpXQX_div(Sp, f, T, p);
    3511                 :         14 :   FpXQX_roots_edf(f, xp, Xp, T, p, V, idx);
    3512                 :         14 :   FpXQX_roots_edf(ff,xp, Xp, T, p, V, idx+degpol(f));
    3513                 :            : }
    3514                 :            : 
    3515                 :            : static GEN
    3516                 :          7 : FpXQX_roots_ddf(GEN f, GEN xp, GEN T, GEN p)
    3517                 :            : {
    3518                 :            :   GEN X, Xp, Xq, g, V;
    3519                 :          7 :   long n, dT = get_FpX_degree(T);
    3520                 :          7 :   GEN R = FpXQX_easyroots(f, T, p);
    3521         [ -  + ]:          7 :   if (R) return R;
    3522                 :          7 :   X  = pol_x(varn(f));
    3523                 :          7 :   Xp = FpXQXQ_pow(X, p, f, T, p);
    3524                 :          7 :   Xq = gel(FpXQXQV_autpow(mkvec2(xp, Xp), dT, f, T, p), 2);
    3525                 :          7 :   g = FpXQX_gcd(FpXX_sub(Xq, X, p), f, T, p);
    3526                 :          7 :   n = degpol(g);
    3527         [ -  + ]:          7 :   if (n==0) return cgetg(1, t_COL);
    3528                 :          7 :   g = FpXQX_normalize(g, T, p);
    3529                 :          7 :   V = cgetg(n+1,t_COL);
    3530                 :          7 :   FpXQX_roots_edf(g, xp, Xp, T, p, V, 1);
    3531                 :          7 :   return V;
    3532                 :            : }
    3533                 :            : 
    3534                 :            : /* do not handle small p */
    3535                 :            : static GEN
    3536                 :      16677 : FpXQX_roots_i(GEN S, GEN T, GEN p)
    3537                 :            : {
    3538                 :            :   GEN xp, F, M, V, R;
    3539                 :            :   long i, j, s, l;
    3540         [ +  + ]:      16677 :   if (lgefint(p)==3)
    3541                 :            :   {
    3542                 :      16550 :     ulong pp = p[2];
    3543         [ +  + ]:      16550 :     if (pp == 2)
    3544                 :            :     {
    3545                 :       2548 :       GEN V = F2xqX_roots_i(ZXX_to_F2xX(S,get_FpX_var(T)), ZX_to_F2x(get_FpX_mod(T)));
    3546                 :       2548 :       return F2xC_to_ZXC(V);
    3547                 :            :     }
    3548                 :            :     else
    3549                 :            :     {
    3550                 :      14002 :       GEN V = FlxqX_roots_i(ZXX_to_FlxX(S,pp,get_FpX_var(T)), ZXT_to_FlxT(T,pp), pp);
    3551                 :      14002 :       return FlxC_to_ZXC(V);
    3552                 :            :     }
    3553                 :            :   }
    3554                 :        127 :   S = FpXQX_red(S, T, p);
    3555         [ -  + ]:        127 :   if (!signe(S)) pari_err_ROOTS0("FpXQX_roots");
    3556         [ -  + ]:        127 :   if (degpol(S)==0) return cgetg(1, t_COL);
    3557                 :        127 :   S = FpXQX_normalize(S, T, p);
    3558                 :        127 :   R = FpXQX_easyroots(S, T, p);
    3559         [ +  + ]:        127 :   if (R) return gen_sort(R, (void*) &cmp_RgX, &cmp_nodata);
    3560                 :          7 :   xp = FpX_Frobenius(T, p);
    3561                 :          7 :   V = FpXQX_factor_squarefree(S, T, p);
    3562                 :          7 :   l = lg(V);
    3563         [ +  + ]:         84 :   for (s=0, i=1; i < l; i++)
    3564                 :         77 :     s += !!degpol(gel(V,i));
    3565                 :          7 :   F = cgetg(s+1, t_VEC);
    3566         [ +  + ]:         84 :   for (i=1, j=1; i < l; i++)
    3567         [ +  + ]:         77 :     if (degpol(gel(V,i)))
    3568                 :          7 :       gel(F, j++) = FpXQX_roots_ddf(gel(V,i), xp, T, p);
    3569                 :          7 :   M = shallowconcat1(F);
    3570                 :          7 :   gen_sort_inplace(M, (void*) &cmp_RgX, &cmp_nodata, NULL);
    3571                 :      16677 :   return M;
    3572                 :            : }
    3573                 :            : 
    3574                 :            : GEN
    3575                 :      57435 : F2xqX_roots(GEN x, GEN T)
    3576                 :            : {
    3577                 :      57435 :   pari_sp av = avma;
    3578                 :      57435 :   return gerepilecopy(av, F2xqX_roots_i(x, T));
    3579                 :            : }
    3580                 :            : 
    3581                 :            : GEN
    3582                 :      46239 : FlxqX_roots(GEN x, GEN T, ulong p)
    3583                 :            : {
    3584                 :      46239 :   pari_sp av = avma;
    3585         [ -  + ]:      46239 :   if (p==2)
    3586                 :            :   {
    3587                 :          0 :     GEN V = F2xqX_roots_i(FlxX_to_F2xX(x), Flx_to_F2x(get_Flx_mod(T)));
    3588                 :          0 :     return gerepileupto(av, F2xC_to_FlxC(V));
    3589                 :            :   }
    3590                 :      46239 :   return gerepilecopy(av, FlxqX_roots_i(x, T, p));
    3591                 :            : }
    3592                 :            : 
    3593                 :            : GEN
    3594                 :      16649 : FpXQX_roots(GEN x, GEN T, GEN p)
    3595                 :            : {
    3596                 :      16649 :   pari_sp av = avma;
    3597                 :      16649 :   return gerepilecopy(av, FpXQX_roots_i(x, T, p));
    3598                 :            : }
    3599                 :            : 
    3600                 :            : static long
    3601                 :         28 : FqX_sqf_split(GEN *t0, GEN q, GEN T, GEN p)
    3602                 :            : {
    3603                 :         28 :   GEN *t = t0, u = *t, v, S, g, X;
    3604                 :         28 :   long d, dg, N = degpol(u);
    3605                 :            : 
    3606         [ -  + ]:         28 :   if (N == 1) return 1;
    3607                 :         28 :   v = X = pol_x(varn(u));
    3608                 :         28 :   S = FqX_Frobenius_powers(u, T, p);
    3609         [ +  + ]:         77 :   for (d=1; d <= N>>1; d++)
    3610                 :            :   {
    3611                 :         49 :     v = FqX_Frobenius_eval(v, S, u, T, p);
    3612                 :         49 :     g = FpXQX_normalize(FqX_gcd(FpXX_sub(v,X,p),u, T,p),T,p);
    3613         [ +  + ]:         49 :     dg = degpol(g); if (dg <= 0) continue;
    3614                 :            : 
    3615                 :            :     /* all factors of g have degree d */
    3616                 :         28 :     *t = g;
    3617                 :         28 :     FqX_split(t, d, q, S, T, p);
    3618                 :         28 :     t += dg / d;
    3619                 :         28 :     N -= dg;
    3620         [ -  + ]:         28 :     if (N)
    3621                 :            :     {
    3622                 :          0 :       u = FqX_div(u,g, T,p);
    3623                 :          0 :       v = FqX_rem(v,u, T,p);
    3624                 :            :     }
    3625                 :            :   }
    3626         [ -  + ]:         28 :   if (N) *t++ = u;
    3627                 :         28 :   return t - t0;
    3628                 :            : }
    3629                 :            : 
    3630                 :            : /* not memory-clean */
    3631                 :            : static GEN
    3632                 :        427 : FpX_factorff_i(GEN P, GEN T, GEN p)
    3633                 :            : {
    3634                 :        427 :   GEN V, E, F = FpX_factor(P,p);
    3635                 :        427 :   long i, lfact = 1, nmax = lgpol(P), n = lgcols(F);
    3636                 :            : 
    3637                 :        427 :   V = cgetg(nmax,t_VEC);
    3638                 :        427 :   E = cgetg(nmax,t_VECSMALL);
    3639         [ +  + ]:       1043 :   for(i=1;i<n;i++)
    3640                 :            :   {
    3641                 :        616 :     GEN R = FpX_factorff_irred(gmael(F,1,i),T,p), e = gmael(F,2,i);
    3642                 :        616 :     long j, r = lg(R);
    3643         [ +  + ]:       5474 :     for (j=1; j<r; j++,lfact++)
    3644                 :            :     {
    3645                 :       4858 :       gel(V,lfact) = gel(R,j);
    3646                 :       4858 :       gel(E,lfact) = e;
    3647                 :            :     }
    3648                 :            :   }
    3649                 :        427 :   setlg(V,lfact);
    3650                 :        427 :   setlg(E,lfact); return sort_factor_pol(mkvec2(V,E), cmp_RgX);
    3651                 :            : }
    3652                 :            : GEN
    3653                 :          0 : FpX_factorff(GEN P, GEN T, GEN p)
    3654                 :            : {
    3655                 :          0 :   pari_sp av = avma;
    3656                 :          0 :   return gerepilecopy(av, FpX_factorff_i(P, T, p));
    3657                 :            : }
    3658                 :            : 
    3659                 :            : /* assumes varncmp (varn(T), varn(f)) > 0 */
    3660                 :            : static GEN
    3661                 :        637 : FpXQX_factor_i(GEN f, GEN T, GEN p)
    3662                 :            : {
    3663                 :        637 :   long pg, j, k, e, N, lfact, pk, d = degpol(f);
    3664                 :            :   GEN E, f2, f3, df1, df2, g1, u, q, t;
    3665                 :            : 
    3666      [ +  +  + ]:        637 :   switch(d)
    3667                 :            :   {
    3668                 :          7 :     case -1: retmkmat2(mkcolcopy(f), mkvecsmall(1));
    3669                 :          7 :     case 0: return trivial_fact();
    3670                 :            :   }
    3671                 :        623 :   T = FpX_normalize(T, p);
    3672                 :        623 :   f = FpXQX_normalize(f, T, p);
    3673         [ +  + ]:        623 :   if (isabsolutepol(f)) return FpX_factorff_i(simplify_shallow(f), T, p);
    3674         [ +  + ]:        196 :   if (degpol(f)==2)
    3675                 :            :   {
    3676                 :         91 :     long v = varn(f);
    3677                 :         91 :     GEN r = FpXQX_quad_roots(f,T,p);
    3678   [ +  +  +  - ]:         91 :     switch(lg(r)-1)
    3679                 :            :     {
    3680                 :            :     case 0:
    3681                 :          7 :       return mkvec2(mkcolcopy(f), mkvecsmall(1));
    3682                 :            :     case 1:
    3683                 :         42 :       return mkvec2(mkcol(deg1pol_shallow(gen_1, Fq_neg(gel(r,1), T, p), v)),
    3684                 :            :                     mkvecsmall(2));
    3685                 :            :     case 2:
    3686                 :            :       {
    3687                 :         42 :         GEN f1 = deg1pol_shallow(gen_1, Fq_neg(gel(r,1), T, p), v);
    3688                 :         42 :         GEN f2 = deg1pol_shallow(gen_1, Fq_neg(gel(r,2), T, p), v);
    3689                 :         42 :         t = mkcol2(f1, f2); E = mkvecsmall2(1, 1);
    3690                 :         42 :         sort_factor_pol(mkvec2(t, E), cmp_RgX);
    3691                 :         42 :         return mkvec2(t, E);
    3692                 :            :       }
    3693                 :            :     }
    3694                 :            :   }
    3695                 :            : 
    3696                 :        105 :   pg = itos_or_0(p);
    3697                 :        105 :   df2  = NULL; /* gcc -Wall */
    3698                 :        105 :   t = cgetg(d+1,t_VEC);
    3699                 :        105 :   E = cgetg(d+1, t_VECSMALL);
    3700                 :            : 
    3701                 :        105 :   q = powiu(p, degpol(T));
    3702                 :        105 :   e = lfact = 1;
    3703                 :        105 :   pk = 1;
    3704                 :        105 :   f3 = NULL;
    3705                 :        105 :   df1 = FqX_deriv(f, T, p);
    3706                 :            :   for(;;)
    3707                 :            :   {
    3708                 :            :     long nb0;
    3709         [ +  + ]:        252 :     while (!signe(df1))
    3710                 :            :     { /* needs d >= p: pg = 0 can't happen  */
    3711                 :         56 :       pk *= pg; e = pk;
    3712                 :         56 :       f = FqX_frob_deflate(f, T, p);
    3713                 :         56 :       df1 = FqX_deriv(f, T, p); f3 = NULL;
    3714                 :            :     }
    3715         [ +  + ]:        196 :     f2 = f3? f3: FqX_gcd(f,df1, T,p);
    3716         [ +  + ]:        196 :     if (!degpol(f2)) u = f;
    3717                 :            :     else
    3718                 :            :     {
    3719                 :         91 :       g1 = FqX_div(f,f2, T,p);
    3720                 :         91 :       df2 = FqX_deriv(f2, T,p);
    3721         [ +  + ]:         91 :       if (gequal0(df2)) { u = g1; f3 = f2; }
    3722                 :            :       else
    3723                 :            :       {
    3724                 :         70 :         f3 = FqX_gcd(f2,df2, T,p);
    3725         [ +  + ]:         70 :         u = degpol(f3)? FqX_div(f2, f3, T,p): f2;
    3726                 :         70 :         u = FqX_div(g1, u, T,p);
    3727                 :            :       }
    3728                 :            :     }
    3729                 :            :     /* u is square-free (product of irreducibles of multiplicity e) */
    3730                 :        196 :     N = degpol(u);
    3731         [ +  + ]:        196 :     if (N) {
    3732                 :        147 :       nb0 = lfact;
    3733                 :        147 :       gel(t,lfact) = FpXQX_normalize(u, T,p);
    3734         [ +  + ]:        147 :       if (N == 1) lfact++;
    3735                 :            :       else
    3736                 :            :       {
    3737         [ +  + ]:         84 :         if (!equaliu(p,2))
    3738                 :         56 :           lfact += FqX_split_Berlekamp(&gel(t,lfact), T, p);
    3739                 :            :         else
    3740                 :         28 :           lfact += FqX_sqf_split(&gel(t,lfact), q, T, p);
    3741                 :            :       }
    3742         [ +  + ]:        455 :       for (j = nb0; j < lfact; j++) E[j] = e;
    3743                 :            :     }
    3744                 :            : 
    3745         [ +  + ]:        196 :     if (!degpol(f2)) break;
    3746                 :         91 :     f = f2; df1 = df2; e += pk;
    3747                 :         91 :   }
    3748                 :        105 :   setlg(t, lfact);
    3749                 :        105 :   setlg(E, lfact);
    3750         [ +  + ]:        413 :   for (j=1; j<lfact; j++) gel(t,j) = FpXQX_normalize(gel(t,j), T,p);
    3751                 :        105 :   (void)sort_factor_pol(mkvec2(t, E), cmp_RgX);
    3752                 :        105 :   k = 1;
    3753         [ +  + ]:        308 :   for (j = 2; j < lfact; j++)
    3754                 :            :   {
    3755         [ +  + ]:        203 :     if (RgX_equal(gel(t,j), gel(t,k)))
    3756                 :         21 :       E[k] += E[j];
    3757                 :            :     else
    3758                 :            :     { /* new factor */
    3759                 :        182 :       k++;
    3760                 :        182 :       E[k] = E[j];
    3761                 :        182 :       gel(t,k) = gel(t,j);
    3762                 :            :     }
    3763                 :            :   }
    3764                 :        105 :   setlg(t, k+1);
    3765                 :        637 :   setlg(E, k+1); return mkvec2(t, E);
    3766                 :            : }
    3767                 :            : 
    3768                 :            : long
    3769                 :        168 : FqX_ispower(GEN f, ulong k, GEN T, GEN p, GEN *pt)
    3770                 :            : {
    3771                 :        168 :   pari_sp av = avma;
    3772                 :            :   long v, w;
    3773                 :            :   ulong pp;
    3774                 :            :   GEN lc, F;
    3775                 :            : 
    3776         [ -  + ]:        168 :   if (degpol(f) % k) return 0;
    3777                 :        168 :   lc = leading_coeff(f);
    3778                 :        168 :   lc = Fq_sqrtn(lc, stoi(k), T, p, NULL);
    3779         [ -  + ]:        168 :   if (!lc) { av = avma; return 0; }
    3780                 :        168 :   pp = itou_or_0(p);
    3781                 :        168 :   f = FqX_normalize(f, T, p);
    3782         [ +  - ]:        168 :   v = pp? u_lvalrem(k,pp,&k): 0;
    3783         [ +  + ]:        168 :   if (v)
    3784                 :            :   {
    3785                 :            :     long i;
    3786                 :         63 :     w = u_lval(RgX_deflate_order(f), pp);
    3787         [ +  + ]:         63 :     if (w < v) { avma = av; return 0; }
    3788                 :            :     /* deflate as much as possible using frobenius, unless k reduced to 1 */
    3789         [ +  + ]:         56 :     if (k == 1) w = v;
    3790                 :         56 :     f = RgX_deflate(f, upowuu(pp,w));
    3791 [ +  + ][ +  + ]:        119 :     if (T) for (i = 0; i < w; i++) f = FqX_frobinv_inplace(f, T, p);
    3792                 :         56 :     w -= v;
    3793                 :            :   }
    3794                 :            :   else
    3795                 :        105 :     w = 0;
    3796                 :            :   /* k coprime to p; true f we're testing is f^(p^w) */
    3797         [ +  + ]:        161 :   if (k == 1)
    3798                 :         28 :     F = f;
    3799                 :            :   else
    3800                 :            :   {
    3801                 :        133 :     ulong pow = upowuu(pp,w);
    3802         [ +  - ]:        133 :     F = pt? pol_1(varn(f)): NULL;
    3803         [ +  + ]:        308 :     while (degpol(f) > 0)
    3804                 :            :     {
    3805                 :        196 :       GEN gk, g, df = FqX_deriv(f, T, p);
    3806                 :            :       long v;
    3807         [ +  + ]:        196 :       if (!signe(df)) { pow *= pp; f = FqX_frob_deflate(f, T, p); continue; }
    3808                 :        133 :       g = FqX_div(f, FqX_normalize(FqX_gcd(f,df,T,p),T,p), T,p);
    3809                 :            :       /* g | f is squarefree,monic; remove (g^k)^oo from f */
    3810                 :        133 :       gk = FqX_powu(g, k, T,p);
    3811                 :        133 :       v = 0;
    3812                 :        133 :       for(v = 0;; v++)
    3813                 :            :       {
    3814                 :        329 :         GEN q = FqX_divrem(f, gk, T,p, ONLY_DIVIDES);
    3815         [ +  + ]:        329 :         if (!q) break;
    3816                 :        196 :         f = q;
    3817                 :        196 :       }
    3818                 :            :       /* some factor from g remains in f ? */
    3819 [ +  + ][ -  + ]:        133 :       if (!v || degpol(FqX_gcd(f,g,T,p))) { avma = av; return 0; }
    3820         [ +  - ]:        112 :       if (F) F = FqX_mul(F, FqX_powu(g, v*pow, T,p), T,p);
    3821                 :            :     }
    3822                 :            :   }
    3823         [ +  - ]:        140 :   if (pt) *pt = gerepileupto(av, FqX_Fq_mul(F, lc, T,p));
    3824                 :        168 :   return 1;
    3825                 :            : }
    3826                 :            : 
    3827                 :            : static void
    3828                 :        273 : ffcheck(pari_sp *av, GEN *f, GEN *T, GEN p)
    3829                 :            : {
    3830                 :            :   long v;
    3831         [ -  + ]:        273 :   if (typ(*T)!=t_POL) pari_err_TYPE("factorff",*T);
    3832         [ -  + ]:        273 :   if (typ(*f)!=t_POL) pari_err_TYPE("factorff",*f);
    3833         [ -  + ]:        273 :   if (typ(p)!=t_INT) pari_err_TYPE("factorff",p);
    3834                 :        273 :   v = varn(*T);
    3835         [ -  + ]:        273 :   if (varncmp(v, varn(*f)) <= 0)
    3836                 :          0 :     pari_err_PRIORITY("factorff", *T, "<=", varn(*f));
    3837                 :        273 :   *T = RgX_to_FpX(*T, p); *av = avma;
    3838                 :        273 :   *f = RgX_to_FqX(*f, *T,p);
    3839                 :        273 : }
    3840                 :            : GEN
    3841                 :        287 : factorff(GEN f, GEN p, GEN T)
    3842                 :            : {
    3843                 :            :   pari_sp av;
    3844                 :            :   GEN z;
    3845 [ +  + ][ -  + ]:        287 :   if (!p || !T)
    3846                 :            :   {
    3847                 :            :     long pa, t;
    3848         [ -  + ]:         42 :     if (typ(f) != t_POL) pari_err_TYPE("factorff",f);
    3849                 :         42 :     T = p = NULL;
    3850                 :         42 :     t = RgX_type(f, &p, &T, &pa);
    3851         [ -  + ]:         42 :     if (t != t_FFELT) pari_err_TYPE("factorff",f);
    3852                 :         42 :     return FFX_factor(f,T);
    3853                 :            :   }
    3854                 :        245 :   ffcheck(&av, &f, &T, p); z = FpXQX_factor_i(f, T, p);
    3855                 :        287 :   return to_Fq_fact(gel(z,1),gel(z,2), T,p, av);
    3856                 :            : }
    3857                 :            : GEN
    3858                 :     103369 : polrootsff(GEN f, GEN p, GEN T)
    3859                 :            : {
    3860                 :            :   pari_sp av;
    3861                 :            :   GEN z;
    3862 [ +  + ][ -  + ]:     103369 :   if (!p || !T)
    3863                 :            :   {
    3864                 :            :     long pa, t;
    3865         [ -  + ]:     103341 :     if (typ(f) != t_POL) pari_err_TYPE("polrootsff",f);
    3866                 :     103341 :     T = p = NULL;
    3867                 :     103341 :     t = RgX_type(f, &p, &T, &pa);
    3868         [ -  + ]:     103341 :     if (t != t_FFELT) pari_err_TYPE("polrootsff",f);
    3869                 :     103341 :     return FFX_roots(f,T);
    3870                 :            :   }
    3871                 :         28 :   ffcheck(&av, &f, &T, p); z = FpXQX_roots_i(f, T, p);
    3872                 :     103369 :   return to_FqC(z, T,p, av);
    3873                 :            : }
    3874                 :            : 
    3875                 :            : /* factorization of x modulo (T,p). Assume x already reduced */
    3876                 :            : GEN
    3877                 :        392 : FpXQX_factor(GEN x, GEN T, GEN p)
    3878                 :            : {
    3879                 :        392 :   pari_sp av = avma;
    3880                 :        392 :   return gerepilecopy(av, FpXQX_factor_i(x, T, p));
    3881                 :            : }
    3882                 :            : /* See also: Isomorphisms between finite field and relative
    3883                 :            :  * factorization in polarit3.c */

Generated by: LCOV version 1.9