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 19070-36a960b) Lines: 1384 1518 91.2 %
Date: 2016-06-30 Functions: 115 123 93.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 784 1045 75.0 %

           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                 :        266 : factmod_init(GEN *F, GEN p)
      32                 :            : {
      33         [ -  + ]:        266 :   if (typ(p)!=t_INT) pari_err_TYPE("factmod",p);
      34         [ -  + ]:        266 :   if (signe(p) < 0) pari_err_PRIME("factmod",p);
      35         [ -  + ]:        266 :   if (typ(*F)!=t_POL) pari_err_TYPE("factmod",*F);
      36         [ +  + ]:        266 :   if (lgefint(p) == 3)
      37                 :            :   {
      38                 :        161 :     ulong pp = p[2];
      39         [ -  + ]:        161 :     if (pp < 2) pari_err_PRIME("factmod", p);
      40                 :        161 :     *F = RgX_to_Flx(*F, pp);
      41         [ +  + ]:        161 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42                 :            :   }
      43                 :            :   else
      44                 :            :   {
      45                 :        105 :     *F = RgX_to_FpX(*F, p);
      46         [ +  + ]:        105 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      47                 :            :   }
      48                 :        266 : }
      49                 :            : /* as above, assume p prime and *F a ZX */
      50                 :            : static void
      51                 :     612239 : ZX_factmod_init(GEN *F, GEN p)
      52                 :            : {
      53         [ +  + ]:     612239 :   if (lgefint(p) == 3)
      54                 :            :   {
      55                 :     605710 :     ulong pp = p[2];
      56                 :     605710 :     *F = ZX_to_Flx(*F, pp);
      57         [ +  + ]:     605710 :     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                 :     612239 : }
      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                 :       6363 : Flx_root_mod_2(GEN f)
     161                 :            : {
     162                 :       6363 :   int z1, z0 = !(f[2] & 1);
     163                 :            :   long i,n;
     164                 :            :   GEN y;
     165                 :            : 
     166         [ +  + ]:      24612 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     167                 :       6363 :   z1 = n & 1;
     168                 :       6363 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     169         [ +  + ]:       6363 :   if (z0) y[i++] = 0;
     170         [ +  + ]:       6363 :   if (z1) y[i  ] = 1;
     171                 :       6363 :   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         [ +  + ]:    1816409 : 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                 :      60872 : rootmod_aux(GEN f, GEN pp, GEN (*Roots)(GEN,ulong), int gpwrap)
     196                 :            : {
     197                 :      60872 :   pari_sp av = avma;
     198                 :            :   GEN y;
     199         [ +  + ]:      60872 :   if (gpwrap)
     200                 :         84 :     factmod_init(&f, pp);
     201                 :            :   else
     202                 :      60788 :     ZX_factmod_init(&f, pp);
     203      [ +  +  + ]:      60872 :   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         [ +  + ]:      60823 :   if (typ(f) == t_VECSMALL)
     209                 :            :   {
     210                 :      58026 :     ulong p = pp[2];
     211         [ +  + ]:      58026 :     if (p == 2)
     212                 :       6363 :       y = Flx_root_mod_2(f);
     213                 :            :     else
     214                 :            :     {
     215         [ -  + ]:      51663 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     216                 :      51663 :       y = Roots(f, p);
     217                 :            :     }
     218                 :      58019 :     y = Flc_to_ZC(y);
     219                 :            :   }
     220                 :            :   else
     221                 :       2797 :     y = FpX_roots_i(f, pp);
     222         [ +  + ]:      60809 :   if (gpwrap) y = FpC_to_mod(y, pp);
     223                 :      60844 :   return gerepileupto(av, y);
     224                 :            : }
     225                 :            : /* assume that f is a ZX an pp a prime */
     226                 :            : GEN
     227                 :      60788 : FpX_roots(GEN f, GEN pp)
     228                 :      60788 : { return rootmod_aux(f, pp, Flx_roots_i, 0); }
     229                 :            : /* no assumptions on f and pp */
     230                 :            : GEN
     231                 :         14 : rootmod2(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_naive, 1); }
     232                 :            : GEN
     233                 :         70 : rootmod(GEN f, GEN pp) { return rootmod_aux(f, pp, &Flx_roots_i, 1); }
     234                 :            : GEN
     235                 :         84 : rootmod0(GEN f, GEN p, long flag)
     236                 :            : {
     237      [ +  +  - ]:         84 :   switch(flag)
     238                 :            :   {
     239                 :         70 :     case 0: return rootmod(f,p);
     240                 :         14 :     case 1: return rootmod2(f,p);
     241                 :          0 :     default: pari_err_FLAG("polrootsmod");
     242                 :            :   }
     243                 :         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                 :       6482 : FpX_quad_root(GEN x, GEN p, int unknown)
     257                 :            : {
     258                 :       6482 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     259                 :            : 
     260         [ -  + ]:       6482 :   if (equaliu(p, 2)) {
     261         [ #  # ]:          0 :     if (!signe(b)) return c;
     262         [ #  # ]:          0 :     return signe(c)? NULL: gen_1;
     263                 :            :   }
     264                 :       6482 :   D = subii(sqri(b), shifti(c,2));
     265                 :       6482 :   D = remii(D,p);
     266 [ +  + ][ +  + ]:       6482 :   if (unknown && kronecker(D,p) == -1) return NULL;
     267                 :            : 
     268                 :       6474 :   s = Fp_sqrt(D,p);
     269                 :            :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     270         [ +  + ]:       6474 :   if (!s) return NULL;
     271                 :       6482 :   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                 :    5976852 : Fl_disc_bc(ulong b, ulong c, ulong p)
     280                 :    5976852 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     281                 :            : /* p > 2 */
     282                 :            : static ulong
     283                 :    5902338 : Flx_quad_root(GEN x, ulong p, int unknown)
     284                 :            : {
     285                 :    5902338 :   ulong s, b = x[3], c = x[2];
     286                 :    5902338 :   ulong D = Fl_disc_bc(b, c, p);
     287 [ +  + ][ +  + ]:    5902809 :   if (unknown && krouu(D,p) == -1) return p;
     288                 :    5789952 :   s = Fl_sqrt(D,p);
     289         [ +  + ]:    5782015 :   if (s==~0UL) return p;
     290                 :    5897520 :   return Fl_halve(Fl_sub(s,b, p), p);
     291                 :            : }
     292                 :            : static ulong
     293                 :    5251968 : Flx_otherroot(GEN x, ulong r, ulong p)
     294                 :    5251968 : { 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                 :     614061 : split_init(struct split_t *S, long max)
     302                 :            : {
     303                 :     614061 :   S->todo = vectrunc_init(max);
     304                 :     614052 :   S->done = vectrunc_init(max);
     305                 :     614060 : }
     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                 :     671921 : 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                 :     671916 : split_moveto_done(struct split_t *S, long i, GEN t)
     323                 :            : {
     324                 :     671916 :   long n = lg(S->todo)-1;
     325                 :     671916 :   vectrunc_append(S->done, t);
     326         [ +  - ]:     671911 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     327                 :     671911 :   setlg(S->todo, n);
     328                 :            : 
     329                 :     671911 : }
     330                 :            : /* append t to done */
     331                 :            : static void
     332                 :      99925 : split_add_done(struct split_t *S, GEN t)
     333                 :      99925 : { vectrunc_append(S->done, t); }
     334                 :            : /* split todo[i] into a and b */
     335                 :            : static void
     336                 :      75062 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     337                 :            : {
     338                 :      75062 :   gel(S->todo, i) = a;
     339                 :      75062 :   split_add(S, b);
     340                 :      75064 : }
     341                 :            : /* split todo[i] into a and b, moved to done */
     342                 :            : static void
     343                 :      81029 : split_done(struct split_t *S, long i, GEN a, GEN b)
     344                 :            : {
     345                 :      81029 :   split_moveto_done(S, i, a);
     346                 :      81030 :   split_add_done(S, b);
     347                 :      81029 : }
     348                 :            : 
     349                 :            : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     350                 :            : static GEN
     351                 :       2797 : FpX_roots_i(GEN f, GEN p)
     352                 :            : {
     353                 :            :   GEN pol, pol0, a, q;
     354                 :            :   struct split_t S;
     355                 :            : 
     356                 :       2797 :   split_init(&S, lg(f)-1);
     357                 :       2797 :   settyp(S.done, t_COL);
     358         [ +  + ]:       2797 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     359   [ +  +  +  + ]:       2797 :   switch(degpol(f))
     360                 :            :   {
     361                 :          7 :     case 0: return ZC_copy(S.done);
     362                 :         14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     363                 :            :     case 2: {
     364                 :       1603 :       GEN s, r = FpX_quad_root(f, p, 1);
     365         [ +  - ]:       1603 :       if (r) {
     366                 :       1603 :         split_add_done(&S, r);
     367                 :       1603 :         s = FpX_otherroot(f,r, p);
     368                 :            :         /* f not known to be square free yet */
     369         [ +  - ]:       1603 :         if (!equalii(r, s)) split_add_done(&S, s);
     370                 :            :       }
     371                 :       1603 :       return sort(S.done);
     372                 :            :     }
     373                 :            :   }
     374                 :            : 
     375                 :       1173 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     376         [ -  + ]:       1173 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     377                 :       1173 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     378                 :       1173 :   a = FpX_gcd(f,a, p);
     379         [ -  + ]:       1173 :   if (!degpol(a)) return ZC_copy(S.done);
     380                 :       1173 :   split_add(&S, FpX_normalize(a,p));
     381                 :            : 
     382                 :       1173 :   q = shifti(p,-1);
     383                 :       1173 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     384                 :       1173 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     385                 :       1173 :   for (pol0[2] = 1;; pol0[2]++)
     386                 :            :   {
     387                 :       2444 :     long j, l = lg(S.todo);
     388         [ +  + ]:       2444 :     if (l == 1) return sort(S.done);
     389 [ -  + ][ #  # ]:       1278 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     390         [ +  + ]:       2647 :     for (j = 1; j < l; j++)
     391                 :            :     {
     392                 :       1376 :       GEN c = gel(S.todo,j);
     393      [ +  +  + ]:       1376 :       switch(degpol(c))
     394                 :            :       { /* convert linear and quadratics to roots, try to split the rest */
     395                 :            :         case 1:
     396                 :         63 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     397                 :         63 :           j--; l--; break;
     398                 :            :         case 2: {
     399                 :       1208 :           GEN r = FpX_quad_root(c, p, 0), s;
     400         [ +  + ]:       1208 :           if (!r) pari_err_PRIME("polrootsmod",p);
     401                 :       1201 :           s = FpX_otherroot(c,r, p);
     402                 :       1201 :           split_done(&S, j, r, s);
     403                 :       1201 :           j--; l--; break;
     404                 :            :         }
     405                 :            :         default: {
     406                 :            :           /* b = pol^(p-1)/2 - 1 */
     407                 :        105 :           GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, c,p), gen_1, p);
     408                 :            :           long db;
     409                 :        105 :           b = FpX_gcd(c,b, p); db = degpol(b);
     410 [ +  - ][ +  + ]:        105 :           if (db && db < degpol(c))
     411                 :            :           {
     412                 :         98 :             b = FpX_normalize(b, p);
     413                 :         98 :             c = FpX_div(c,b, p);
     414                 :         98 :             split_todo(&S, j, b, c);
     415                 :            :           }
     416                 :            :         }
     417                 :            :       }
     418                 :            :     }
     419                 :       4061 :   }
     420                 :            : }
     421                 :            : 
     422                 :            : /* Assume f is normalized */
     423                 :            : static ulong
     424                 :      53915 : Flx_cubic_root(GEN ff, ulong p)
     425                 :            : {
     426                 :      53915 :   GEN f = Flx_normalize(ff,p);
     427                 :      53915 :   ulong pi = get_Fl_red(p);
     428         [ +  + ]:      53915 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     429                 :      53915 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     430                 :      53915 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     431                 :      53915 :   ulong B = Fl_addmul_pre(t, Fl_sub(Fl_double(t2, p), b, p), c, p, pi);
     432                 :      53915 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     433                 :      53915 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     434                 :      53915 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     435                 :      53915 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     436                 :      53915 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     437         [ +  + ]:      53915 :   if (s!=~0UL)
     438                 :            :   {
     439         [ +  + ]:      33392 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     440         [ +  + ]:      33392 :     if (p%3==2) /* 1 solutions */
     441                 :       5166 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     442                 :            :     else
     443                 :            :     {
     444                 :      28226 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     445         [ +  + ]:      28226 :       if (vS1==~0UL) return p; /*0 solutions*/
     446                 :            :       /*3 solutions*/
     447                 :            :     }
     448         [ +  + ]:      22680 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     449                 :      22680 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     450                 :            :   }
     451                 :            :   else
     452                 :            :   {
     453                 :      20523 :     pari_sp av = avma;
     454                 :      20523 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     455                 :      20523 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     456                 :            :     ulong Sa;
     457         [ +  + ]:      20523 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     458                 :      19583 :     Sa = vS1[1];
     459         [ +  + ]:      19583 :     if (p%3==1) /*1 solutions*/
     460                 :            :     {
     461                 :      14132 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     462         [ +  + ]:      14132 :       if (Fa!=P)
     463                 :       9542 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     464                 :            :     }
     465                 :      19583 :     avma = av;
     466                 :      53915 :     return Fl_sub(Fl_double(Sa,p),t,p);
     467                 :            :   }
     468                 :            : }
     469                 :            : 
     470                 :            : /* assume p > 2 prime */
     471                 :            : static ulong
     472                 :     694298 : Flx_oneroot_i(GEN f, ulong p, long fl)
     473                 :            : {
     474                 :            :   GEN pol, a;
     475                 :            :   ulong q;
     476                 :            :   long da;
     477                 :            : 
     478         [ +  + ]:     694298 :   if (Flx_val(f)) return 0;
     479   [ +  +  +  + ]:     694130 :   switch(degpol(f))
     480                 :            :   {
     481                 :        190 :     case 1: return Fl_neg(f[2], p);
     482                 :     540348 :     case 2: return Flx_quad_root(f, p, 1);
     483         [ +  - ]:      50102 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     484                 :            :   }
     485                 :            : 
     486         [ +  + ]:     103484 :   if (!fl)
     487                 :            :   {
     488                 :      93113 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     489         [ -  + ]:      93113 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     490                 :      93113 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     491                 :      93113 :     a = Flx_gcd(f,a, p);
     492                 :      10371 :   } else a = f;
     493                 :     103484 :   da = degpol(a);
     494         [ +  + ]:     103486 :   if (!da) return p;
     495                 :      53591 :   a = Flx_normalize(a,p);
     496                 :            : 
     497                 :      53596 :   q = p >> 1;
     498                 :      53596 :   pol = polx_Flx(f[1]);
     499                 :      53580 :   for(pol[2] = 1;; pol[2]++)
     500                 :            :   {
     501 [ -  + ][ #  # ]:      90756 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     502   [ +  +  +  + ]:      90759 :     switch(da)
     503                 :            :     {
     504                 :      22191 :       case 1: return Fl_neg(a[2], p);
     505                 :      27600 :       case 2: return Flx_quad_root(a, p, 0);
     506         [ +  - ]:       3813 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     507                 :            :       default: {
     508                 :      37155 :         GEN b = Flx_Fl_add(Flxq_powu(pol,q, a,p), p-1, p);
     509                 :            :         long db;
     510                 :      37178 :         b = Flx_gcd(a,b, p); db = degpol(b);
     511 [ +  + ][ +  + ]:      37177 :         if (db && db < da)
     512                 :            :         {
     513                 :      35395 :           b = Flx_normalize(b, p);
     514         [ +  + ]:      35394 :           if (db <= (da >> 1)) {
     515                 :      21895 :             a = b;
     516                 :      21895 :             da = db;
     517                 :            :           } else {
     518                 :      13499 :             a = Flx_div(a,b, p);
     519                 :      13499 :             da -= db;
     520                 :            :           }
     521                 :            :         }
     522                 :            :       }
     523                 :            :     }
     524                 :     731488 :   }
     525                 :            : }
     526                 :            : 
     527                 :            : /* assume p > 2 prime */
     528                 :            : static GEN
     529                 :       3649 : FpX_oneroot_i(GEN f, GEN p)
     530                 :            : {
     531                 :            :   GEN pol, pol0, a, q;
     532                 :            :   long da;
     533                 :            : 
     534         [ +  + ]:       3649 :   if (ZX_val(f)) return gen_0;
     535      [ +  +  + ]:       3642 :   switch(degpol(f))
     536                 :            :   {
     537                 :        140 :     case 1: return subii(p, gel(f,2));
     538                 :       3432 :     case 2: return FpX_quad_root(f, p, 1);
     539                 :            :   }
     540                 :            : 
     541                 :         70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     542         [ -  + ]:         70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     543                 :         70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     544                 :         70 :   a = FpX_gcd(f,a, p);
     545                 :         70 :   da = degpol(a);
     546         [ -  + ]:         70 :   if (!da) return NULL;
     547                 :         70 :   a = FpX_normalize(a,p);
     548                 :            : 
     549                 :         70 :   q = shifti(p,-1);
     550                 :         70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     551                 :         70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     552                 :         70 :   for (pol0[2]=1; ; pol0[2]++)
     553                 :            :   {
     554 [ -  + ][ #  # ]:        224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     555      [ +  +  + ]:        224 :     switch(da)
     556                 :            :     {
     557                 :         42 :       case 1: return subii(p, gel(a,2));
     558                 :         28 :       case 2: return FpX_quad_root(a, p, 0);
     559                 :            :       default: {
     560                 :        154 :         GEN b = FpX_Fp_sub_shallow(FpXQ_pow(pol,q, a,p), gen_1, p);
     561                 :            :         long db;
     562                 :        154 :         b = FpX_gcd(a,b, p); db = degpol(b);
     563 [ +  + ][ +  - ]:        154 :         if (db && db < da)
     564                 :            :         {
     565                 :        147 :           b = FpX_normalize(b, p);
     566         [ +  + ]:        147 :           if (db <= (da >> 1)) {
     567                 :        105 :             a = b;
     568                 :        105 :             da = db;
     569                 :            :           } else {
     570                 :         42 :             a = FpX_div(a,b, p);
     571                 :         42 :             da -= db;
     572                 :            :           }
     573                 :            :         }
     574                 :            :       }
     575                 :            :     }
     576                 :       3803 :   }
     577                 :            : }
     578                 :            : 
     579                 :            : ulong
     580                 :     676500 : Flx_oneroot(GEN f, ulong p)
     581                 :            : {
     582                 :     676500 :   pari_sp av = avma;
     583                 :            :   ulong r;
     584      [ -  -  + ]:     676500 :   switch(lg(f))
     585                 :            :   {
     586                 :          0 :     case 2: return 0;
     587                 :          0 :     case 3: avma = av; return p;
     588                 :            :   }
     589         [ -  + ]:     676500 :   if (p == 2) return Flx_oneroot_mod_2(f);
     590                 :     676500 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     591                 :     676500 :   avma = av; return r;
     592                 :            : }
     593                 :            : 
     594                 :            : ulong
     595                 :      10454 : Flx_oneroot_split(GEN f, ulong p)
     596                 :            : {
     597                 :      10454 :   pari_sp av = avma;
     598                 :            :   ulong r;
     599      [ -  -  + ]:      10454 :   switch(lg(f))
     600                 :            :   {
     601                 :          0 :     case 2: return 0;
     602                 :          0 :     case 3: avma = av; return p;
     603                 :            :   }
     604         [ -  + ]:      10454 :   if (p == 2) return Flx_oneroot_mod_2(f);
     605                 :      10454 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     606                 :      10471 :   avma = av; return r;
     607                 :            : }
     608                 :            : 
     609                 :            : /* assume that p is prime */
     610                 :            : GEN
     611                 :      11004 : FpX_oneroot(GEN f, GEN pp) {
     612                 :      11004 :   pari_sp av = avma;
     613                 :      11004 :   ZX_factmod_init(&f, pp);
     614      [ -  -  + ]:      11004 :   switch(lg(f))
     615                 :            :   {
     616                 :          0 :     case 2: avma = av; return gen_0;
     617                 :          0 :     case 3: avma = av; return NULL;
     618                 :            :   }
     619         [ +  + ]:      11004 :   if (typ(f) == t_VECSMALL)
     620                 :            :   {
     621                 :       7355 :     ulong r, p = pp[2];
     622         [ +  + ]:       7355 :     if (p == 2)
     623                 :         14 :       r = Flx_oneroot_mod_2(f);
     624                 :            :     else
     625                 :       7341 :       r = Flx_oneroot_i(f, p, 0);
     626                 :       7355 :     avma = av;
     627         [ +  + ]:       7355 :     return (r == p)? NULL: utoi(r);
     628                 :            :   }
     629                 :       3649 :   f = FpX_oneroot_i(f, pp);
     630         [ -  + ]:       3649 :   if (!f) { avma = av; return NULL; }
     631                 :      11004 :   return gerepileuptoint(av, f);
     632                 :            : }
     633                 :            : 
     634                 :            : /*******************************************************************/
     635                 :            : /*                                                                 */
     636                 :            : /*                     FACTORISATION MODULO p                      */
     637                 :            : /*                                                                 */
     638                 :            : /*******************************************************************/
     639                 :            : 
     640                 :            : /* Functions giving information on the factorisation. */
     641                 :            : 
     642                 :            : /* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */
     643                 :            : GEN
     644                 :         55 : FpX_Berlekamp_ker(GEN u, GEN p)
     645                 :            : {
     646                 :         55 :   pari_sp ltop=avma;
     647                 :         55 :   long j,N = degpol(u);
     648                 :         55 :   GEN Q  = FpX_matFrobenius(u, p);
     649         [ +  + ]:       2139 :   for (j=1; j<=N; j++)
     650                 :       2084 :     gcoeff(Q,j,j) = Fp_sub(gcoeff(Q,j,j), gen_1, p);
     651                 :         55 :   return gerepileupto(ltop, FpM_ker(Q,p));
     652                 :            : }
     653                 :            : 
     654                 :            : GEN
     655                 :      20853 : F2x_Berlekamp_ker(GEN u)
     656                 :            : {
     657                 :      20853 :   pari_sp ltop=avma;
     658                 :      20853 :   long j,N = F2x_degree(u);
     659                 :            :   GEN Q;
     660                 :            :   pari_timer T;
     661                 :      20853 :   timer_start(&T);
     662                 :      20853 :   Q = F2x_matFrobenius(u);
     663         [ +  + ]:     114954 :   for (j=1; j<=N; j++)
     664                 :      94101 :     F2m_flip(Q,j,j);
     665         [ -  + ]:      20853 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     666                 :      20853 :   Q = F2m_ker_sp(Q,0);
     667         [ -  + ]:      20853 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     668                 :      20853 :   return gerepileupto(ltop,Q);
     669                 :            : }
     670                 :            : 
     671                 :            : GEN
     672                 :     205561 : Flx_Berlekamp_ker(GEN u, ulong l)
     673                 :            : {
     674                 :     205561 :   pari_sp ltop=avma;
     675                 :     205561 :   long j,N = degpol(u);
     676                 :            :   GEN Q;
     677                 :            :   pari_timer T;
     678                 :     205561 :   timer_start(&T);
     679                 :     205561 :   Q  = Flx_matFrobenius(u, l);
     680         [ +  + ]:     967409 :   for (j=1; j<=N; j++)
     681                 :     761848 :     coeff(Q,j,j) = Fl_sub(coeff(Q,j,j),1,l);
     682         [ -  + ]:     205561 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
     683                 :     205561 :   Q = Flm_ker_sp(Q,l,0);
     684         [ -  + ]:     205561 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
     685                 :     205561 :   return gerepileupto(ltop,Q);
     686                 :            : }
     687                 :            : 
     688                 :            : /* product of terms of degree 1 in factorization of f */
     689                 :            : GEN
     690                 :      60389 : FpX_split_part(GEN f, GEN p)
     691                 :            : {
     692                 :      60389 :   long n = degpol(f);
     693                 :      60389 :   GEN z, X = pol_x(varn(f));
     694         [ +  + ]:      60389 :   if (n <= 1) return f;
     695                 :      60361 :   f = FpX_red(f, p);
     696                 :      60361 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     697                 :      60389 :   return FpX_gcd(z,f,p);
     698                 :            : }
     699                 :            : 
     700                 :            : /* Compute the number of roots in Fp without counting multiplicity
     701                 :            :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     702                 :            : long
     703                 :      18459 : FpX_nbroots(GEN f, GEN p)
     704                 :            : {
     705                 :      18459 :   pari_sp av = avma;
     706                 :      18459 :   GEN z = FpX_split_part(f, p);
     707                 :      18459 :   avma = av; return degpol(z);
     708                 :            : }
     709                 :            : 
     710                 :            : int
     711                 :          0 : FpX_is_totally_split(GEN f, GEN p)
     712                 :            : {
     713                 :          0 :   long n=degpol(f);
     714                 :          0 :   pari_sp av = avma;
     715         [ #  # ]:          0 :   if (n <= 1) return 1;
     716         [ #  # ]:          0 :   if (cmpui(n, p) > 0) return 0;
     717                 :          0 :   f = FpX_red(f, p);
     718                 :          0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     719                 :            : }
     720                 :            : 
     721                 :            : long
     722                 :    8114928 : Flx_nbroots(GEN f, ulong p)
     723                 :            : {
     724                 :    8114928 :   long n = degpol(f);
     725                 :    8113110 :   pari_sp av = avma;
     726                 :            :   GEN z;
     727         [ +  + ]:    8113110 :   if (n <= 1) return n;
     728         [ +  + ]:    8112935 :   if (n == 2)
     729                 :            :   {
     730                 :            :     ulong D;
     731         [ +  + ]:    4890065 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     732                 :    4889232 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     733                 :    4897221 :     return 1 + krouu(D,p);
     734                 :            :   }
     735                 :    3222870 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     736                 :    3222846 :   z = Flx_gcd(z, f, p);
     737                 :    8123174 :   avma = av; return degpol(z);
     738                 :            : }
     739                 :            : 
     740                 :            : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     741                 :            : static GEN
     742                 :       4102 : FpX_ddf(GEN T, GEN XP, GEN p)
     743                 :            : {
     744                 :       4102 :   pari_sp av = avma;
     745                 :            :   GEN b, g, h, F, f, Tr, xq;
     746                 :            :   long i, j, n, v;
     747                 :            :   long B, l, m;
     748                 :            :   pari_timer ti;
     749                 :       4102 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     750         [ -  + ]:       4102 :   if (n == 0) return cgetg(1, t_VEC);
     751         [ +  + ]:       4102 :   if (n == 1) return mkvec(get_FpX_mod(T));
     752                 :       4095 :   B = n/2;
     753                 :       4095 :   l = usqrt(B);
     754                 :       4095 :   m = (B+l-1)/l;
     755                 :       4095 :   T = FpX_get_red(T, p);
     756                 :       4095 :   b = cgetg(l+2, t_VEC);
     757                 :       4095 :   gel(b, 1) = pol_x(v);
     758                 :       4095 :   gel(b, 2) = XP;
     759         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     760                 :       4095 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     761         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq baby");
     762         [ +  + ]:      10661 :   for (i = 3; i <= l+1; i++)
     763                 :       6566 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     764         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: baby");
     765                 :       4095 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     766         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: xq giant");
     767                 :       4095 :   g = cgetg(m+1, t_VEC);
     768                 :       4095 :   gel(g, 1) = gel(xq, 2);
     769         [ +  + ]:      14931 :   for(i = 2; i <= m; i++)
     770                 :      10836 :     gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     771         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: giant");
     772                 :       4095 :   h = cgetg(m+1, t_VEC);
     773         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     774                 :            :   {
     775                 :      14931 :     pari_sp av = avma;
     776                 :      14931 :     GEN gj = gel(g, j);
     777                 :      14931 :     GEN e = FpX_sub(gj, gel(b, 1), p);
     778         [ +  + ]:      42938 :     for (i = 2; i <= l; i++)
     779                 :      28007 :       e = FpXQ_mul(e, FpX_sub(gj, gel(b, i), p), T, p);
     780                 :      14931 :     gel(h, j) = gerepileupto(av, e);
     781                 :            :   }
     782         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: diff");
     783                 :       4095 :   Tr = get_FpX_mod(T);
     784                 :       4095 :   F = cgetg(m+1, t_VEC);
     785         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     786                 :            :   {
     787                 :      14931 :     gel(F, j) = FpX_gcd(Tr, gel(h, j), p);
     788                 :      14931 :     Tr = FpX_div(Tr, gel(F,j), p);
     789                 :            :   }
     790         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: F");
     791                 :       4095 :   f = const_vec(n, pol_1(v));
     792         [ +  + ]:      19026 :   for (j = 1; j <= m; j++)
     793                 :            :   {
     794                 :      14931 :     GEN e = gel(F, j);
     795         [ +  - ]:      16919 :     for (i=l-1; i >= 0; i--)
     796                 :            :     {
     797                 :      16919 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     798         [ +  + ]:      16919 :       if (degpol(u))
     799                 :            :       {
     800                 :       2709 :         gel(f, l*j-i) = u;
     801                 :       2709 :         e = FpX_div(e, u, p);
     802                 :            :       }
     803         [ +  + ]:      16919 :       if (!degpol(e)) break;
     804                 :            :     }
     805                 :            :   }
     806         [ -  + ]:       4095 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf: f");
     807         [ +  + ]:       4095 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     808                 :       4102 :   return gerepilecopy(av, f);
     809                 :            : }
     810                 :            : 
     811                 :            : static void
     812                 :          0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     813                 :            : {
     814                 :          0 :   long n = degpol(Tp), r = n/d;
     815                 :            :   GEN T, f, ff;
     816                 :            :   GEN p2;
     817         [ #  # ]:          0 :   if (r==1) { gel(V, idx) = Tp; return; }
     818                 :          0 :   p2 = shifti(p,-1);
     819                 :          0 :   T = FpX_get_red(Tp, p);
     820                 :          0 :   XP = FpX_rem(XP, T, p);
     821                 :            :   while (1)
     822                 :            :   {
     823                 :          0 :     pari_sp btop = avma;
     824                 :            :     long i;
     825                 :          0 :     GEN g = random_FpX(n, varn(Tp), p);
     826                 :          0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     827         [ #  # ]:          0 :     if (signe(t) == 0) continue;
     828         [ #  # ]:          0 :     for(i=1; i<=10; i++)
     829                 :            :     {
     830                 :          0 :       pari_sp btop2 = avma;
     831                 :          0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     832                 :          0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     833 [ #  # ][ #  # ]:          0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     834                 :          0 :       avma = btop2;
     835                 :            :     }
     836 [ #  # ][ #  # ]:          0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     837                 :          0 :     avma = btop;
     838                 :          0 :   }
     839                 :          0 :   f = FpX_normalize(f, p);
     840                 :          0 :   ff = FpX_div(Tp, f ,p);
     841                 :          0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     842                 :          0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     843                 :            : }
     844                 :            : 
     845                 :            : static void
     846                 :        210 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     847                 :            : {
     848                 :            :   pari_sp av;
     849                 :        210 :   GEN Tp = get_FpX_mod(T);
     850                 :        210 :   long n = degpol(hp), vT = varn(Tp);
     851                 :            :   GEN u1, u2, f1, f2;
     852                 :            :   GEN R, h;
     853                 :        210 :   h = FpX_get_red(hp, p);
     854                 :        210 :   t = FpX_rem(t, T, p);
     855                 :        210 :   av = avma;
     856                 :            :   do
     857                 :            :   {
     858                 :        343 :     avma = av;
     859                 :        343 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     860                 :        343 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     861 [ +  + ][ +  + ]:        343 :   } while (degpol(u1)==0 || degpol(u1)==n);
     862                 :        210 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     863                 :        210 :   f1 = FpX_normalize(f1, p);
     864                 :        210 :   u2 = FpX_div(hp, u1, p);
     865                 :        210 :   f2 = FpX_div(Tp, f1, p);
     866         [ +  + ]:        210 :   if (degpol(u1)==1)
     867                 :        105 :     gel(V, idx) = f1;
     868                 :            :   else
     869                 :        105 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     870                 :        210 :   idx += degpol(f1)/d;
     871         [ +  + ]:        210 :   if (degpol(u2)==1)
     872                 :        119 :     gel(V, idx) = f2;
     873                 :            :   else
     874                 :         91 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     875                 :        210 : }
     876                 :            : 
     877                 :            : static void
     878                 :         14 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     879                 :            : {
     880                 :         14 :   long n = degpol(Tp), r = n/d, vT = varn(Tp);
     881                 :            :   GEN T, h, t;
     882                 :            :   pari_timer ti;
     883         [ -  + ]:         28 :   if (r==1) { gel(V, idx) = Tp; return; }
     884                 :         14 :   T = FpX_get_red(Tp, p);
     885                 :         14 :   XP = FpX_rem(XP, T, p);
     886         [ -  + ]:         14 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     887                 :            :   do
     888                 :            :   {
     889                 :         14 :     GEN g = random_FpX(n, vT, p);
     890                 :         14 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     891         [ -  + ]:         14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     892                 :         14 :     h = FpXQ_minpoly(t, T, p);
     893         [ -  + ]:         14 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     894         [ -  + ]:         14 :   } while (degpol(h) != r);
     895                 :         14 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     896                 :            : }
     897                 :            : 
     898                 :            : static GEN
     899                 :         14 : FpX_factor_Shoup(GEN T, GEN p)
     900                 :            : {
     901                 :         14 :   long i, n, s = 0;
     902                 :            :   GEN XP, D, V;
     903                 :         14 :   long e = expi(p);
     904                 :            :   pari_timer ti;
     905                 :         14 :   n = get_FpX_degree(T);
     906                 :         14 :   T = FpX_get_red(T, p);
     907         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     908                 :         14 :   XP = FpX_Frobenius(T, p);
     909         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     910                 :         14 :   D = FpX_ddf(T, XP, p);
     911         [ -  + ]:         14 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     912         [ +  + ]:        301 :   for (i = 1; i <= n; i++)
     913                 :        287 :     s += degpol(gel(D,i))/i;
     914                 :         14 :   V = cgetg(s+1, t_COL);
     915         [ +  + ]:        301 :   for (i = 1, s = 1; i <= n; i++)
     916                 :            :   {
     917                 :        287 :     GEN Di = gel(D,i);
     918                 :        287 :     long ni = degpol(Di), ri = ni/i;
     919         [ +  + ]:        287 :     if (ni == 0) continue;
     920                 :         21 :     Di = FpX_normalize(Di, p);
     921         [ +  + ]:         21 :     if (ni == i) { gel(V, s++) = Di; continue; }
     922         [ +  - ]:         14 :     if (ri <= e*expu(e))
     923                 :         14 :       FpX_edf(Di, XP, i, p, V, s);
     924                 :            :     else
     925                 :          0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     926         [ -  + ]:         14 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     927                 :         14 :     s += ri;
     928                 :            :   }
     929                 :         14 :   return V;
     930                 :            : }
     931                 :            : 
     932                 :            : static GEN
     933                 :          0 : FpX_simplefact_Shoup(GEN T, GEN p)
     934                 :            : {
     935                 :          0 :   long i, n, s = 0, j = 1, k;
     936                 :            :   GEN XP, D, V;
     937                 :            :   pari_timer ti;
     938                 :          0 :   n = get_FpX_degree(T);
     939                 :          0 :   T = FpX_get_red(T, p);
     940         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     941                 :          0 :   XP = FpX_Frobenius(T, p);
     942         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     943                 :          0 :   D = FpX_ddf(T, XP, p);
     944         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
     945         [ #  # ]:          0 :   for (i = 1; i <= n; i++)
     946                 :          0 :     s += degpol(gel(D,i))/i;
     947                 :          0 :   V = cgetg(s+1, t_VEC);
     948         [ #  # ]:          0 :   for (i = 1; i <= n; i++)
     949                 :            :   {
     950                 :          0 :     long ni = degpol(gel(D,i)), ri = ni/i;
     951         [ #  # ]:          0 :     if (ni == 0) continue;
     952         [ #  # ]:          0 :     for (k = 1; k <= ri; k++)
     953                 :          0 :       gel(V, j++) = utoi(i);
     954                 :            :   }
     955                 :          0 :   return V;
     956                 :            : }
     957                 :            : 
     958                 :            : /* Yun algorithm: Assume p > degpol(T) */
     959                 :            : 
     960                 :            : static GEN
     961                 :          7 : FpX_factor_Yun(GEN T, GEN p)
     962                 :            : {
     963                 :          7 :   pari_sp av = avma;
     964                 :          7 :   long n = degpol(T);
     965                 :          7 :   long i = 1;
     966                 :          7 :   GEN d = FpX_deriv(T, p);
     967                 :            :   GEN a, b, c;
     968                 :          7 :   GEN V = cgetg(n+1,t_VEC);
     969                 :          7 :   a = FpX_gcd(T, d, p);
     970         [ -  + ]:          7 :   if (degpol(a) == 0) return mkvec(T);
     971                 :          7 :   b = FpX_div(T, a, p);
     972                 :            :   do
     973                 :            :   {
     974                 :         14 :     c = FpX_div(d, a, p);
     975                 :         14 :     d = FpX_sub(c, FpX_deriv(b, p), p);
     976                 :         14 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
     977                 :         14 :     gel(V, i++) = a;
     978                 :         14 :     b = FpX_div(b, a, p);
     979         [ +  + ]:         14 :   } while (degpol(b));
     980                 :          7 :   setlg(V, i);
     981                 :          7 :   return gerepilecopy(av, V);
     982                 :            : }
     983                 :            : 
     984                 :            : static GEN
     985                 :          7 : FpX_factor_Cantor(GEN T, GEN p)
     986                 :            : {
     987                 :            :   GEN E, F, M, V;
     988                 :            :   long i, j, s, l;
     989                 :          7 :   V = FpX_factor_Yun(T, p);
     990                 :          7 :   l = lg(V);
     991         [ +  + ]:         21 :   for (s=0, i=1; i < l; i++)
     992                 :         14 :     s += !!degpol(gel(V,i));
     993                 :          7 :   F = cgetg(s+1, t_VEC);
     994                 :          7 :   E = cgetg(s+1, t_VEC);
     995         [ +  + ]:         21 :   for (i=1, j=1; i < l; i++)
     996         [ +  - ]:         14 :     if (degpol(gel(V,i)))
     997                 :            :     {
     998                 :         14 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
     999                 :         14 :       gel(F, j) = Fj;
    1000                 :         14 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1001                 :         14 :       j++;
    1002                 :            :     }
    1003                 :          7 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1004                 :          7 :   return sort_factor_pol(M, cmpii);
    1005                 :            : }
    1006                 :            : 
    1007                 :            : static GEN
    1008                 :          0 : FpX_simplefact_Cantor(GEN T, GEN p)
    1009                 :            : {
    1010                 :            :   GEN E, F, M, V;
    1011                 :            :   long i, j, s, l;
    1012                 :          0 :   V = FpX_factor_Yun(get_FpX_mod(T), p);
    1013                 :          0 :   l = lg(V);
    1014         [ #  # ]:          0 :   for (s=0, i=1; i < l; i++)
    1015                 :          0 :     s += !!degpol(gel(V,i));
    1016                 :          0 :   F = cgetg(s+1, t_VEC);
    1017                 :          0 :   E = cgetg(s+1, t_VEC);
    1018         [ #  # ]:          0 :   for (i=1, j=1; i < l; i++)
    1019         [ #  # ]:          0 :     if (degpol(gel(V,i)))
    1020                 :            :     {
    1021                 :          0 :       GEN Fj = FpX_simplefact_Shoup(gel(V,i), p);
    1022                 :          0 :       gel(F, j) = Fj;
    1023                 :          0 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1024                 :          0 :       j++;
    1025                 :            :     }
    1026                 :          0 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1027                 :          0 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1028                 :            : }
    1029                 :            : 
    1030                 :            : static int
    1031                 :          0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1032                 :            : {
    1033                 :          0 :   pari_sp av = avma;
    1034                 :            :   pari_timer ti;
    1035                 :            :   long n, d;
    1036                 :          0 :   GEN T = get_FpX_mod(Tp);
    1037                 :          0 :   GEN dT = FpX_deriv(T, p);
    1038                 :            :   GEN XP, D;
    1039         [ #  # ]:          0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1040                 :          0 :   n = get_FpX_degree(T);
    1041                 :          0 :   T = FpX_get_red(Tp, p);
    1042         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1043                 :          0 :   XP = FpX_Frobenius(T, p);
    1044         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1045                 :          0 :   D = FpX_ddf(T, XP, p);
    1046         [ #  # ]:          0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf");
    1047                 :          0 :   d = degpol(gel(D, n));
    1048                 :          0 :   avma = av; return d==n;
    1049                 :            : }
    1050                 :            : 
    1051                 :            : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1052                 :            : 
    1053                 :            : /*Assume that p is large and odd*/
    1054                 :            : static GEN
    1055                 :         35 : FpX_factcantor_i(GEN f, GEN pp, long flag)
    1056                 :            : {
    1057                 :         35 :   long d = degpol(f);
    1058         [ +  + ]:         35 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1059      [ +  -  - ]:          7 :   switch(flag)
    1060                 :            :   {
    1061                 :          7 :     default: return FpX_factor_Cantor(f, pp);
    1062                 :          0 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1063         [ #  # ]:         35 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1064                 :            :   }
    1065                 :            : }
    1066                 :            : 
    1067                 :            : long
    1068                 :       4088 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1069                 :            : {
    1070                 :       4088 :   pari_sp av = avma;
    1071                 :       4088 :   GEN ddf = FpX_ddf(T, XP, p);
    1072                 :       4088 :   long l = lg(ddf), i, s=0;
    1073         [ +  + ]:      86716 :   for(i = 1; i < l; i++)
    1074                 :      82628 :     s += degpol(gel(ddf,i))/i;
    1075                 :       4088 :   avma = av; return s;
    1076                 :            : }
    1077                 :            : 
    1078                 :            : long
    1079                 :         28 : FpX_nbfact(GEN T, GEN p)
    1080                 :            : {
    1081                 :         28 :   pari_sp av = avma;
    1082                 :         28 :   GEN XP = FpX_Frobenius(T, p);
    1083                 :         28 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1084                 :         28 :   avma = av; return n;
    1085                 :            : }
    1086                 :            : 
    1087                 :            : /* INPUT:
    1088                 :            :  *  m integer (converted to polynomial w in Z[X] by stoFpX)
    1089                 :            :  *  p prime; q = (p^d-1) / 2^r
    1090                 :            :  *  t[0] polynomial of degree k*d product of k irreducible factors of degree d
    1091                 :            :  *  t[0] is expected to be normalized (leading coeff = 1)
    1092                 :            :  * OUTPUT:
    1093                 :            :  *  t[0],t[1]...t[k-1] the k factors, normalized */
    1094                 :            : 
    1095                 :            : static void
    1096                 :      21385 : F2x_split(ulong m, GEN *t, long d)
    1097                 :            : {
    1098                 :            :   long l, v, dv;
    1099                 :            :   pari_sp av0, av;
    1100                 :            :   GEN w;
    1101                 :            : 
    1102         [ +  + ]:      23737 :   dv = F2x_degree(*t); if (dv==d) return;
    1103                 :       2352 :   v=(*t)[1]; av0=avma;
    1104                 :       2352 :   for(av=avma;;avma=av)
    1105                 :            :   {
    1106                 :       3969 :     GEN w0 = w = F2xq_powu(polx_F2x(v), m-1, *t); m += 2;
    1107         [ +  + ]:      27321 :     for (l=1; l<d; l++) w = F2x_add(w0, F2xq_sqr(w, *t));
    1108                 :       3969 :     w = F2x_gcd(*t,w);
    1109 [ +  + ][ +  + ]:       3969 :     l = F2x_degree(w); if (l && l!=dv) break;
    1110                 :       1617 :   }
    1111                 :       2352 :   w = gerepileupto(av0, w);
    1112                 :       2352 :   l /= d; t[l]=F2x_div(*t,w); *t=w;
    1113                 :       2352 :   F2x_split(m,t+l,d);
    1114                 :       2352 :   F2x_split(m,t,  d);
    1115                 :            : }
    1116                 :            : 
    1117                 :            : /* p > 2 */
    1118                 :            : static GEN
    1119                 :          7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1120                 :            : {
    1121      [ -  -  + ]:          7 :   switch(d)
    1122                 :            :   {
    1123                 :            :     case -1:
    1124                 :          0 :     case 0: return NULL;
    1125                 :          0 :     case 1: return gen_1;
    1126                 :            :   }
    1127         [ -  + ]:          7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1128                 :            : }
    1129                 :            : /* p > 2 */
    1130                 :            : static GEN
    1131                 :         14 : FpX_degfact_2(GEN f, GEN p, long d)
    1132                 :            : {
    1133   [ -  -  -  + ]:         14 :   switch(d)
    1134                 :            :   {
    1135                 :          0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1136                 :          0 :     case 0: return trivial_fact();
    1137                 :          0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1138                 :            :   }
    1139      [ +  +  - ]:         14 :   switch(FpX_quad_factortype(f, p)) {
    1140                 :          7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1141                 :          7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1142                 :         14 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1143                 :            :   }
    1144                 :            : }
    1145                 :            : 
    1146                 :            : GEN
    1147                 :         70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1148                 :            : GEN
    1149                 :        763 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1150                 :            : 
    1151                 :            : /* Mod(0,p) * x, where x is f's main variable */
    1152                 :            : static GEN
    1153                 :         14 : Mod0pX(GEN f, GEN p)
    1154                 :         14 : { return scalarpol(mkintmod(gen_0, p), varn(f)); }
    1155                 :            : static GEN
    1156                 :         14 : zero_fact_intmod(GEN f, GEN p) { return prime_fact(Mod0pX(f,p)); }
    1157                 :            : 
    1158                 :            : /* not gerepile safe */
    1159                 :            : static GEN
    1160                 :         74 : FpX_factor_2(GEN f, GEN p, long d)
    1161                 :            : {
    1162                 :            :   GEN r, s, R, S;
    1163                 :            :   long v;
    1164                 :            :   int sgn;
    1165   [ -  +  +  + ]:         74 :   switch(d)
    1166                 :            :   {
    1167                 :          0 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1168                 :          2 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1169                 :         10 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1170                 :            :   }
    1171                 :         62 :   r = FpX_quad_root(f, p, 1);
    1172         [ +  + ]:         62 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1173                 :         53 :   v = varn(f);
    1174                 :         53 :   s = FpX_otherroot(f, r, p);
    1175         [ +  - ]:         53 :   if (signe(r)) r = subii(p, r);
    1176         [ +  - ]:         53 :   if (signe(s)) s = subii(p, s);
    1177         [ +  + ]:         53 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1178                 :         53 :   R = deg1pol_shallow(gen_1, r, v);
    1179         [ +  + ]:         53 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1180                 :         46 :   S = deg1pol_shallow(gen_1, s, v);
    1181                 :         74 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1182                 :            : }
    1183                 :            : static GEN
    1184                 :         95 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1185                 :            : {
    1186      [ +  +  + ]:         95 :   switch(flag) {
    1187                 :          7 :     case 2: return FpX_is_irred_2(f, p, d);
    1188                 :         14 :     case 1: return FpX_degfact_2(f, p, d);
    1189                 :         95 :     default: return FpX_factor_2(f, p, d);
    1190                 :            :   }
    1191                 :            : }
    1192                 :            : 
    1193                 :            : static int
    1194                 :      27832 : F2x_quad_factortype(GEN x)
    1195         [ +  + ]:      27832 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1196                 :            : 
    1197                 :            : static GEN
    1198                 :          7 : F2x_is_irred_2(GEN f, long d)
    1199 [ +  - ][ +  - ]:          7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
                 [ +  - ]
    1200                 :            : 
    1201                 :            : static GEN
    1202                 :        861 : F2x_degfact_2(GEN f, long d)
    1203                 :            : {
    1204         [ -  + ]:        861 :   if (!d) return trivial_fact();
    1205         [ +  + ]:        861 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1206      [ +  +  + ]:        742 :   switch(F2x_quad_factortype(f)) {
    1207                 :        133 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1208                 :        161 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1209                 :        861 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1210                 :            :   }
    1211                 :            : }
    1212                 :            : 
    1213                 :            : static GEN
    1214                 :      46123 : F2x_factor_2(GEN f, long d)
    1215                 :            : {
    1216                 :      46123 :   long v = f[1];
    1217         [ -  + ]:      46123 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1218         [ +  + ]:      46123 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1219         [ +  + ]:      38598 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1220      [ +  +  + ]:      19817 :   switch(F2x_quad_factortype(f))
    1221                 :            :   {
    1222                 :       3836 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1223                 :      13454 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1224                 :      46123 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1225                 :            :   }
    1226                 :            : }
    1227                 :            : static GEN
    1228                 :      46991 : F2x_factor_deg2(GEN f, long d, long flag)
    1229                 :            : {
    1230      [ +  +  + ]:      46991 :   switch(flag) {
    1231                 :          7 :     case 2: return F2x_is_irred_2(f, d);
    1232                 :        861 :     case 1: return F2x_degfact_2(f, d);
    1233                 :      46991 :     default: return F2x_factor_2(f, d);
    1234                 :            :   }
    1235                 :            : }
    1236                 :            : 
    1237                 :            : static void
    1238                 :         19 : split_squares(struct split_t *S, GEN g, ulong p)
    1239                 :            : {
    1240                 :         19 :   ulong q = p >> 1;
    1241                 :         19 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1242         [ +  + ]:         19 :   if (degpol(a) < 0)
    1243                 :            :   {
    1244                 :            :     ulong i;
    1245                 :          6 :     split_add_done(S, (GEN)1);
    1246         [ +  + ]:         18 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1247                 :            :   } else {
    1248                 :         13 :     (void)Flx_valrem(a, &a);
    1249         [ +  - ]:         13 :     if (degpol(a))
    1250                 :            :     {
    1251                 :         13 :       a = Flx_gcd(a, Flx_Xnm1(g[1], q, p), p);
    1252         [ +  - ]:         13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1253                 :            :     }
    1254                 :            :   }
    1255                 :         19 : }
    1256                 :            : static void
    1257                 :         19 : split_nonsquares(struct split_t *S, GEN g, ulong p)
    1258                 :            : {
    1259                 :         19 :   ulong q = p >> 1;
    1260                 :         19 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1261         [ +  + ]:         19 :   if (degpol(a) < 0)
    1262                 :            :   {
    1263                 :          6 :     ulong i, z = Fl_nonsquare(p);
    1264                 :          6 :     split_add_done(S, (GEN)z);
    1265         [ +  + ]:         18 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1266                 :            :   } else {
    1267                 :         13 :     (void)Flx_valrem(a, &a);
    1268         [ +  - ]:         13 :     if (degpol(a))
    1269                 :            :     {
    1270                 :         13 :       a = Flx_gcd(a, Flx_Xn1(g[1], q, p), p);
    1271         [ +  - ]:         13 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1272                 :            :     }
    1273                 :            :   }
    1274                 :         19 : }
    1275                 :            : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1276                 :            :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1277                 :            : static int
    1278                 :     611259 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1279                 :            : {
    1280                 :     611259 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1281                 :     611264 :   long d = degpol(g);
    1282         [ +  + ]:     611259 :   if (d < 0) return 0;
    1283         [ -  + ]:     611224 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1284         [ +  + ]:     611223 :   if (!d) return 1;
    1285         [ +  + ]:     595683 :   if (p <= 1.4 * (ulong)d) {
    1286                 :            :     /* small p; split further using x^((p-1)/2) +/- 1.
    1287                 :            :      * 30% degree drop makes the extra gcd worth it. */
    1288                 :         19 :     split_squares(S, g, p);
    1289                 :         19 :     split_nonsquares(S, g, p);
    1290                 :            :   } else { /* large p; use x^(p-1) - 1 directly */
    1291                 :     595664 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1292         [ -  + ]:     595669 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1293                 :     595669 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1294                 :     595660 :     g = Flx_gcd(g,a, p);
    1295         [ +  - ]:     595661 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1296                 :            :   }
    1297                 :     611254 :   return 1;
    1298                 :            : }
    1299                 :            : 
    1300                 :            : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1301                 :            : static GEN
    1302                 :    5645298 : Flx_roots_i(GEN f, ulong p)
    1303                 :            : {
    1304                 :            :   GEN pol, g;
    1305                 :    5645298 :   long v = Flx_valrem(f, &g);
    1306                 :            :   ulong q;
    1307                 :            :   struct split_t S;
    1308                 :            : 
    1309                 :            :   /* optimization: test for small degree first */
    1310      [ +  +  + ]:    5645081 :   switch(degpol(g))
    1311                 :            :   {
    1312                 :            :     case 1: {
    1313                 :      23555 :       ulong r = p - g[2];
    1314         [ +  + ]:      23555 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1315                 :            :     }
    1316                 :            :     case 2: {
    1317                 :    5010250 :       ulong r = Flx_quad_root(g, p, 1), s;
    1318 [ +  + ][ +  + ]:    5006727 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1319                 :    5006712 :       s = Flx_otherroot(g,r, p);
    1320         [ +  + ]:    5010653 :       if (r < s)
    1321         [ +  + ]:    1252362 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1322         [ +  + ]:    3758291 :       else if (r > s)
    1323         [ +  + ]:    3758228 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1324                 :            :       else
    1325         [ -  + ]:         63 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1326                 :            :     }
    1327                 :            :   }
    1328                 :     611258 :   q = p >> 1;
    1329                 :     611258 :   split_init(&S, lg(f)-1);
    1330                 :     611260 :   settyp(S.done, t_VECSMALL);
    1331         [ +  + ]:     611260 :   if (v) split_add_done(&S, (GEN)0);
    1332         [ +  + ]:     611260 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1333                 :         35 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1334                 :     611222 :   pol = polx_Flx(f[1]);
    1335                 :     611236 :   for (pol[2]=1; ; pol[2]++)
    1336                 :            :   {
    1337                 :    1296130 :     long j, l = lg(S.todo);
    1338         [ +  + ]:    1296130 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1339 [ -  + ][ #  # ]:     684917 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1340         [ +  + ]:    1452479 :     for (j = 1; j < l; j++)
    1341                 :            :     {
    1342                 :     767585 :       GEN c = gel(S.todo,j);
    1343      [ +  +  + ]:     767585 :       switch(degpol(c))
    1344                 :            :       {
    1345                 :            :         case 1:
    1346                 :     590827 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1347                 :     590814 :           j--; l--; break;
    1348                 :            :         case 2: {
    1349                 :      79838 :           ulong r = Flx_quad_root(c, p, 0), s;
    1350         [ +  + ]:      79836 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1351                 :      79829 :           s = Flx_otherroot(c,r, p);
    1352                 :      79829 :           split_done(&S, j, (GEN)r, (GEN)s);
    1353                 :      79828 :           j--; l--; break;
    1354                 :            :         }
    1355                 :            :         default: {
    1356                 :      96921 :           GEN b = Flx_Fl_add(Flxq_powu(pol,q, c,p), p-1, p); /* pol^(p-1)/2 */
    1357                 :            :           long db;
    1358                 :      96921 :           b = Flx_gcd(c,b, p); db = degpol(b);
    1359 [ +  + ][ +  + ]:      96922 :           if (db && db < degpol(c))
    1360                 :            :           {
    1361                 :      74967 :             b = Flx_normalize(b, p);
    1362                 :      74966 :             c = Flx_div(c,b, p);
    1363                 :      74964 :             split_todo(&S, j, b, c);
    1364                 :            :           }
    1365                 :            :         }
    1366                 :            :       }
    1367                 :            :     }
    1368                 :    6324364 :   }
    1369                 :            : }
    1370                 :            : 
    1371                 :            : GEN
    1372                 :    5595308 : Flx_roots(GEN f, ulong p)
    1373                 :            : {
    1374                 :    5595308 :   pari_sp av = avma;
    1375      [ -  +  + ]:    5595308 :   switch(lg(f))
    1376                 :            :   {
    1377                 :          0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1378                 :        546 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1379                 :            :   }
    1380         [ -  + ]:    5594762 :   if (p == 2) return Flx_root_mod_2(f);
    1381                 :    5594762 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1382                 :            : }
    1383                 :            : 
    1384                 :            : /* assume x reduced mod p, monic. */
    1385                 :            : static int
    1386                 :      75488 : Flx_quad_factortype(GEN x, ulong p)
    1387                 :            : {
    1388                 :      75488 :   ulong b = x[3], c = x[2];
    1389                 :      75488 :   return krouu(Fl_disc_bc(b, c, p), p);
    1390                 :            : }
    1391                 :            : static GEN
    1392                 :          0 : Flx_is_irred_2(GEN f, ulong p, long d)
    1393                 :            : {
    1394         [ #  # ]:          0 :   if (!d) return NULL;
    1395         [ #  # ]:          0 :   if (d == 1) return gen_1;
    1396         [ #  # ]:          0 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1397                 :            : }
    1398                 :            : static GEN
    1399                 :      77700 : Flx_degfact_2(GEN f, ulong p, long d)
    1400                 :            : {
    1401         [ -  + ]:      77700 :   if (!d) return trivial_fact();
    1402         [ +  + ]:      77700 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1403      [ +  +  + ]:      75488 :   switch(Flx_quad_factortype(f, p)) {
    1404                 :      35322 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1405                 :      39249 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1406                 :      77700 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1407                 :            :   }
    1408                 :            : }
    1409                 :            : /* p > 2 */
    1410                 :            : static GEN
    1411                 :     258035 : Flx_factor_2(GEN f, ulong p, long d)
    1412                 :            : {
    1413                 :            :   ulong r, s;
    1414                 :            :   GEN R,S;
    1415                 :     258035 :   long v = f[1];
    1416         [ -  + ]:     258035 :   if (d < 0) pari_err(e_ROOTS0,"Flx_factor_2");
    1417         [ +  + ]:     258035 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1418         [ +  + ]:     244338 :   if (d == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1419                 :     168531 :   r = Flx_quad_root(f, p, 1);
    1420         [ +  + ]:     168531 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1421                 :     106387 :   s = Flx_otherroot(f, r, p);
    1422                 :     106387 :   r = Fl_neg(r, p);
    1423                 :     106387 :   s = Fl_neg(s, p);
    1424         [ +  + ]:     106387 :   if (s < r) lswap(s,r);
    1425                 :     106387 :   R = mkvecsmall3(v,r,1);
    1426         [ +  + ]:     106387 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1427                 :      92058 :   S = mkvecsmall3(v,s,1);
    1428                 :     258035 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1429                 :            : }
    1430                 :            : static GEN
    1431                 :     335735 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1432                 :            : {
    1433      [ -  +  + ]:     335735 :   switch(flag) {
    1434                 :          0 :     case 2: return Flx_is_irred_2(f, p, d);
    1435                 :      77700 :     case 1: return Flx_degfact_2(f, p, d);
    1436                 :     335735 :     default: return Flx_factor_2(f, p, d);
    1437                 :            :   }
    1438                 :            : }
    1439                 :            : 
    1440                 :            : static GEN
    1441                 :      14196 : F2x_factcantor_i(GEN f, long flag)
    1442                 :            : {
    1443                 :      14196 :   long j, e, nbfact, d = F2x_degree(f);
    1444                 :            :   ulong k;
    1445                 :            :   GEN X, E, f2, g, g1, u, v, y, t;
    1446                 :            : 
    1447         [ +  + ]:      14196 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    1448                 :            :   /* to hold factors and exponents */
    1449         [ +  + ]:      11368 :   t = flag ? cgetg(d+1,t_VECSMALL): cgetg(d+1,t_VEC);
    1450                 :      11368 :   E = cgetg(d+1, t_VECSMALL);
    1451                 :      11368 :   X = polx_F2x(f[1]);
    1452                 :      11368 :   e = nbfact = 1;
    1453                 :            :   for(;;)
    1454                 :            :   {
    1455                 :      15456 :     f2 = F2x_gcd(f,F2x_deriv(f));
    1456 [ +  + ][ +  + ]:      15456 :     if (flag == 2 && F2x_degree(f2) > 0) return NULL;
    1457                 :      15414 :     g1 = F2x_div(f,f2);
    1458                 :      15414 :     k = 0;
    1459         [ +  + ]:      30611 :     while (F2x_degree(g1)>0)
    1460                 :            :     {
    1461                 :            :       pari_sp av;
    1462                 :            :       long du, dg, dg1;
    1463         [ +  + ]:      15288 :       k++; if (k%2==0) { k++; f2 = F2x_div(f2,g1); }
    1464                 :      15288 :       u = g1; g1 = F2x_gcd(f2,g1);
    1465                 :      15288 :       du = F2x_degree(u);
    1466                 :      15288 :       dg1 = F2x_degree(g1);
    1467         [ +  + ]:      15288 :       if (dg1>0)
    1468                 :            :       {
    1469                 :       1071 :         f2= F2x_div(f2,g1);
    1470         [ +  + ]:       1071 :         if (du == dg1) continue;
    1471                 :        616 :         u = F2x_div( u,g1);
    1472                 :        616 :         du -= dg1;
    1473                 :            :       }
    1474                 :            :       /* here u is square-free (product of irred. of multiplicity e * k) */
    1475                 :      14833 :       v = X;
    1476                 :      14833 :       av = avma;
    1477         [ +  + ]:      75873 :       for (d=1; d <= du>>1; d++)
    1478                 :            :       {
    1479                 :      61131 :         v = F2xq_sqr(v, u);
    1480                 :      61131 :         g = F2x_gcd(F2x_add(v, X), u);
    1481                 :      61131 :         dg = F2x_degree(g);
    1482         [ +  + ]:      61131 :         if (dg <= 0) {avma = (pari_sp)v; v = gerepileuptoleaf(av,v); continue;}
    1483                 :            :         /* g is a product of irred. pols, all of which have degree d */
    1484                 :      17206 :         j = nbfact+dg/d;
    1485         [ +  + ]:      17206 :         if (flag)
    1486                 :            :         {
    1487         [ +  + ]:        525 :           if (flag == 2) return NULL;
    1488         [ +  + ]:        987 :           for ( ; nbfact<j; nbfact++) { t[nbfact]=d; E[nbfact]=e*k; }
    1489                 :            :         }
    1490                 :            :         else
    1491                 :            :         {
    1492                 :      16681 :           gel(t,nbfact) = g;
    1493                 :      16681 :           F2x_split(2,&gel(t,nbfact),d);
    1494         [ +  + ]:      35714 :           for (; nbfact<j; nbfact++) E[nbfact]=e*k;
    1495                 :            :         }
    1496                 :      17115 :         du -= dg;
    1497                 :      17115 :         u = F2x_div(u,g);
    1498                 :      17115 :         v = F2x_rem(v,u);
    1499                 :      17115 :         av = avma;
    1500                 :            :       }
    1501         [ +  + ]:      14742 :       if (du)
    1502                 :            :       {
    1503         [ +  + ]:      13405 :         if (flag) t[nbfact] = du;
    1504                 :      11837 :         else  gel(t,nbfact) = u;
    1505                 :      13405 :         E[nbfact++]=e*k;
    1506                 :            :       }
    1507                 :            :     }
    1508         [ +  + ]:      15323 :     if (F2x_degree(f2)==0) break;
    1509                 :       4088 :     e <<=1; f = F2x_sqrt(f2);
    1510                 :       4088 :   }
    1511         [ +  + ]:      11235 :   if (flag == 2) return gen_1; /* irreducible */
    1512                 :      11228 :   setlg(t, nbfact);
    1513                 :      11228 :   setlg(E, nbfact); y = mkvec2(t, E);
    1514                 :      14196 :   return flag ? sort_factor(y, (void*)cmpGuGu, cmp_nodata)
    1515         [ +  + ]:      11228 :               : sort_factor_pol(y, cmpGuGu);
    1516                 :            : }
    1517                 :            : GEN
    1518                 :       7504 : F2x_factcantor(GEN f, long flag)
    1519                 :            : {
    1520                 :       7504 :   pari_sp av = avma;
    1521                 :       7504 :   GEN z = F2x_factcantor_i(f, flag);
    1522         [ -  + ]:       7504 :   if (flag == 2) { avma = av; return z; }
    1523                 :       7504 :   return gerepilecopy(av, z);
    1524                 :            : }
    1525                 :            : 
    1526                 :            : int
    1527                 :        140 : F2x_is_irred(GEN f) { return !!F2x_factcantor_i(f,2); }
    1528                 :            : 
    1529                 :            : void
    1530                 :      14987 : F2xV_to_FlxV_inplace(GEN v)
    1531                 :            : {
    1532                 :            :   long i;
    1533         [ +  + ]:      32774 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1534                 :      14987 : }
    1535                 :            : void
    1536                 :     464248 : FlxV_to_ZXV_inplace(GEN v)
    1537                 :            : {
    1538                 :            :   long i;
    1539         [ +  + ]:    1241504 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1540                 :     464248 : }
    1541                 :            : void
    1542                 :      78573 : F2xV_to_ZXV_inplace(GEN v)
    1543                 :            : {
    1544                 :            :   long i;
    1545         [ +  + ]:     164676 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1546                 :      78573 : }
    1547                 :            : 
    1548                 :            : /* Adapted from Shoup NTL */
    1549                 :            : static GEN
    1550                 :     311928 : Flx_factor_squarefree(GEN f, ulong p)
    1551                 :            : {
    1552                 :     311928 :   pari_sp av = avma;
    1553                 :            :   GEN r, t, v, tv;
    1554                 :     311928 :   long q, n = degpol(f);
    1555                 :     311928 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1556                 :     311928 :   for(q = 1;;q *= p)
    1557                 :            :   {
    1558                 :     312306 :     r = Flx_gcd(f, Flx_deriv(f, p), p);
    1559                 :     312306 :     t = Flx_div(f, r, p);
    1560         [ +  + ]:     312306 :     if (degpol(t) > 0)
    1561                 :            :     {
    1562                 :            :       long j;
    1563                 :     312103 :       for(j = 1;;j++)
    1564                 :            :       {
    1565                 :     318039 :         v = Flx_gcd(r, t, p);
    1566                 :     318039 :         tv = Flx_div(t, v, p);
    1567         [ +  + ]:     318039 :         if (degpol(tv) > 0)
    1568                 :     315855 :           gel(u, j*q) = tv;
    1569         [ +  + ]:     318039 :         if (degpol(v) <= 0) break;
    1570                 :       5936 :         r = Flx_div(r, v, p);
    1571                 :       5936 :         t = v;
    1572                 :       5936 :       }
    1573         [ +  + ]:     312103 :       if (degpol(r) == 0) break;
    1574                 :            :     }
    1575                 :        378 :     f = Flx_deflate(r, p);
    1576                 :        378 :   }
    1577                 :     311928 :   return gerepilecopy(av, u);
    1578                 :            : }
    1579                 :            : 
    1580                 :            : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1581                 :            : static GEN
    1582                 :     424149 : Flx_ddf(GEN T, GEN XP, ulong p)
    1583                 :            : {
    1584                 :     424149 :   pari_sp av = avma;
    1585                 :            :   GEN b, g, h, F, f, Tr, xq;
    1586                 :            :   long i, j, n, v, bo, ro;
    1587                 :            :   long B, l, m;
    1588                 :            :   pari_timer ti;
    1589                 :     424149 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1590         [ -  + ]:     424149 :   if (n == 0) return cgetg(1, t_VEC);
    1591         [ +  + ]:     424149 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1592                 :     417135 :   B = n/2;
    1593                 :     417135 :   l = usqrt(B);
    1594                 :     417135 :   m = (B+l-1)/l;
    1595                 :     417135 :   T = Flx_get_red(T, p);
    1596                 :     417135 :   b = cgetg(l+2, t_VEC);
    1597                 :     417135 :   gel(b, 1) = polx_Flx(v);
    1598                 :     417135 :   gel(b, 2) = XP;
    1599                 :     417135 :   bo = brent_kung_optpow(n, l-1, 1);
    1600                 :     417135 :   ro = (bo-1) + (l-1)*((n-1)/bo);
    1601         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1602         [ +  + ]:     417135 :   if (expu(p) <= ro)
    1603         [ +  + ]:     178010 :     for (i = 3; i <= l+1; i++)
    1604                 :     109697 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1605                 :            :   else
    1606                 :            :   {
    1607                 :     348822 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1608         [ -  + ]:     348822 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq baby");
    1609         [ +  + ]:     389135 :     for (i = 3; i <= l+1; i++)
    1610                 :      40313 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1611                 :            :   }
    1612         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: baby");
    1613                 :     417135 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1614         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: xq giant");
    1615                 :     417135 :   g = cgetg(m+1, t_VEC);
    1616                 :     417135 :   gel(g, 1) = gel(xq, 2);
    1617         [ +  + ]:     924026 :   for(i = 2; i <= m; i++)
    1618                 :     506891 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1619         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: giant");
    1620                 :     417135 :   h = cgetg(m+1, t_VEC);
    1621         [ +  + ]:    1341161 :   for (j = 1; j <= m; j++)
    1622                 :            :   {
    1623                 :     924026 :     pari_sp av = avma;
    1624                 :     924026 :     GEN gj = gel(g, j);
    1625                 :     924026 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1626         [ +  + ]:    1508918 :     for (i = 2; i <= l; i++)
    1627                 :     584892 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1628                 :     924026 :     gel(h, j) = gerepileupto(av, e);
    1629                 :            :   }
    1630         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: diff");
    1631                 :     417135 :   Tr = get_Flx_mod(T);
    1632                 :     417135 :   F = cgetg(m+1, t_VEC);
    1633         [ +  + ]:    1341161 :   for (j = 1; j <= m; j++)
    1634                 :            :   {
    1635                 :     924026 :     gel(F, j) = Flx_gcd(Tr, gel(h, j), p);
    1636                 :     924026 :     Tr = Flx_div(Tr, gel(F,j), p);
    1637                 :            :   }
    1638         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: F");
    1639                 :     417135 :   f = const_vec(n, pol1_Flx(v));
    1640         [ +  + ]:    1341161 :   for (j = 1; j <= m; j++)
    1641                 :            :   {
    1642                 :     924026 :     GEN e = gel(F, j);
    1643         [ +  - ]:    1029719 :     for (i=l-1; i >= 0; i--)
    1644                 :            :     {
    1645                 :    1029719 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    1646         [ +  + ]:    1029719 :       if (degpol(u))
    1647                 :            :       {
    1648                 :     372457 :         gel(f, l*j-i) = u;
    1649                 :     372457 :         e = Flx_div(e, u, p);
    1650                 :            :       }
    1651         [ +  + ]:    1029719 :       if (!degpol(e)) break;
    1652                 :            :     }
    1653                 :            :   }
    1654         [ -  + ]:     417135 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf: f");
    1655         [ +  + ]:     417135 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    1656                 :     424149 :   return gerepilecopy(av, f);
    1657                 :            : }
    1658                 :            : 
    1659                 :            : static void
    1660                 :      45640 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1661                 :            : {
    1662                 :      45640 :   long n = degpol(Tp), r = n/d;
    1663                 :            :   GEN T, f, ff;
    1664                 :            :   ulong p2;
    1665         [ +  + ]:      66402 :   if (r==1) { gel(V, idx) = Tp; return; }
    1666                 :      20762 :   p2 = p>>1;
    1667                 :      20762 :   T = Flx_get_red(Tp, p);
    1668                 :      20762 :   XP = Flx_rem(XP, T, p);
    1669                 :            :   while (1)
    1670                 :            :   {
    1671                 :      22645 :     pari_sp btop = avma;
    1672                 :            :     long i;
    1673                 :      22645 :     GEN g = random_Flx(n, Tp[1], p);
    1674                 :      22645 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1675         [ +  + ]:      22645 :     if (lgpol(t) == 0) continue;
    1676         [ +  + ]:      46452 :     for(i=1; i<=10; i++)
    1677                 :            :     {
    1678                 :      45045 :       pari_sp btop2 = avma;
    1679                 :      45045 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    1680                 :      45045 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    1681 [ +  + ][ +  + ]:      45045 :       if (degpol(f) > 0 && degpol(f) < n) break;
    1682                 :      24283 :       avma = btop2;
    1683                 :            :     }
    1684 [ +  + ][ +  + ]:      22169 :     if (degpol(f) > 0 && degpol(f) < n) break;
    1685                 :       1407 :     avma = btop;
    1686                 :       1883 :   }
    1687                 :      20762 :   f = Flx_normalize(f, p);
    1688                 :      20762 :   ff = Flx_div(Tp, f ,p);
    1689                 :      20762 :   Flx_edf_simple(f, XP, d, p, V, idx);
    1690                 :      20762 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    1691                 :            : }
    1692                 :            : static void
    1693                 :            : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    1694                 :            : 
    1695                 :            : static void
    1696                 :      33278 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    1697                 :            : {
    1698                 :            :   pari_sp av;
    1699                 :      33278 :   GEN Tp = get_Flx_mod(T);
    1700                 :      33278 :   long n = degpol(hp), vT = Tp[1];
    1701                 :            :   GEN u1, u2, f1, f2;
    1702                 :      33278 :   ulong p2 = p>>1;
    1703                 :            :   GEN R, h;
    1704                 :      33278 :   h = Flx_get_red(hp, p);
    1705                 :      33278 :   t = Flx_rem(t, T, p);
    1706                 :      33278 :   av = avma;
    1707                 :            :   do
    1708                 :            :   {
    1709                 :      56301 :     avma = av;
    1710                 :      56301 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    1711                 :      56301 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    1712 [ +  + ][ +  + ]:      56301 :   } while (degpol(u1)==0 || degpol(u1)==n);
    1713                 :      33278 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    1714                 :      33278 :   f1 = Flx_normalize(f1, p);
    1715                 :      33278 :   u2 = Flx_div(hp, u1, p);
    1716                 :      33278 :   f2 = Flx_div(Tp, f1, p);
    1717         [ +  + ]:      33278 :   if (degpol(u1)==1)
    1718                 :            :   {
    1719         [ +  + ]:      26180 :     if (degpol(f1)==d)
    1720                 :      25718 :       gel(V, idx) = f1;
    1721                 :            :     else
    1722                 :        462 :       Flx_edf(f1, XP, d, p, V, idx);
    1723                 :            :   }
    1724                 :            :   else
    1725                 :       7098 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    1726                 :      33278 :   idx += degpol(f1)/d;
    1727         [ +  + ]:      33278 :   if (degpol(u2)==1)
    1728                 :            :   {
    1729         [ +  + ]:      25893 :     if (degpol(f2)==d)
    1730                 :      25487 :       gel(V, idx) = f2;
    1731                 :            :     else
    1732                 :        406 :       Flx_edf(f2, XP, d, p, V, idx);
    1733                 :            :   }
    1734                 :            :   else
    1735                 :       7385 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    1736                 :      33278 : }
    1737                 :            : 
    1738                 :            : static void
    1739                 :      18795 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    1740                 :            : {
    1741                 :      18795 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    1742                 :            :   GEN T, h, t;
    1743                 :            :   pari_timer ti;
    1744         [ -  + ]:      37590 :   if (r==1) { gel(V, idx) = Tp; return; }
    1745                 :      18795 :   T = Flx_get_red(Tp, p);
    1746                 :      18795 :   XP = Flx_rem(XP, T, p);
    1747         [ -  + ]:      18795 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1748                 :            :   do
    1749                 :            :   {
    1750                 :      20629 :     GEN g = random_Flx(n, vT, p);
    1751                 :      20629 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    1752         [ -  + ]:      20629 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    1753                 :      20629 :     h = Flxq_minpoly(t, T, p);
    1754         [ -  + ]:      20629 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    1755         [ +  + ]:      20629 :   } while (degpol(h) <= 1);
    1756                 :      18795 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    1757                 :            : }
    1758                 :            : 
    1759                 :            : static GEN
    1760                 :      31521 : Flx_factor_Shoup(GEN T, ulong p)
    1761                 :            : {
    1762                 :      31521 :   long i, n, s = 0;
    1763                 :            :   GEN XP, D, V;
    1764                 :      31521 :   long e = expu(p);
    1765                 :            :   pari_timer ti;
    1766                 :      31521 :   n = get_Flx_degree(T);
    1767                 :      31521 :   T = Flx_get_red(T, p);
    1768         [ -  + ]:      31521 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1769                 :      31521 :   XP = Flx_Frobenius(T, p);
    1770         [ -  + ]:      31521 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1771                 :      31521 :   D = Flx_ddf(T, XP, p);
    1772         [ -  + ]:      31521 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1773         [ +  + ]:     345114 :   for (i = 1; i <= n; i++)
    1774                 :     313593 :     s += degpol(gel(D,i))/i;
    1775                 :      31521 :   V = cgetg(s+1, t_COL);
    1776         [ +  + ]:     345114 :   for (i = 1, s = 1; i <= n; i++)
    1777                 :            :   {
    1778                 :     313593 :     GEN Di = gel(D,i);
    1779                 :     313593 :     long ni = degpol(Di), ri = ni/i;
    1780         [ +  + ]:     313593 :     if (ni == 0) continue;
    1781                 :      47586 :     Di = Flx_normalize(Di, p);
    1782         [ +  + ]:      47586 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1783         [ +  + ]:      22043 :     if (ri <= e*expu(e))
    1784                 :      17927 :       Flx_edf(Di, XP, i, p, V, s);
    1785                 :            :     else
    1786                 :       4116 :       Flx_edf_simple(Di, XP, i, p, V, s);
    1787         [ -  + ]:      22043 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    1788                 :      22043 :     s += ri;
    1789                 :            :   }
    1790                 :      31521 :   return V;
    1791                 :            : }
    1792                 :            : 
    1793                 :            : static GEN
    1794                 :     284334 : Flx_simplefact_Shoup(GEN T, ulong p)
    1795                 :            : {
    1796                 :     284334 :   long i, n, s = 0, j = 1, k;
    1797                 :            :   GEN XP, D, V;
    1798                 :            :   pari_timer ti;
    1799                 :     284334 :   n = get_Flx_degree(T);
    1800                 :     284334 :   T = Flx_get_red(T, p);
    1801         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1802                 :     284334 :   XP = Flx_Frobenius(T, p);
    1803         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1804                 :     284334 :   D = Flx_ddf(T, XP, p);
    1805         [ -  + ]:     284334 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1806         [ +  + ]:    1802462 :   for (i = 1; i <= n; i++)
    1807                 :    1518128 :     s += degpol(gel(D,i))/i;
    1808                 :     284334 :   V = cgetg(s+1, t_VECSMALL);
    1809         [ +  + ]:    1802462 :   for (i = 1; i <= n; i++)
    1810                 :            :   {
    1811                 :    1518128 :     long ni = degpol(gel(D,i)), ri = ni/i;
    1812         [ +  + ]:    1518128 :     if (ni == 0) continue;
    1813         [ +  + ]:    1054621 :     for (k = 1; k <= ri; k++)
    1814                 :     689669 :       V[j++] = i;
    1815                 :            :   }
    1816                 :     284334 :   return V;
    1817                 :            : }
    1818                 :            : 
    1819                 :            : static GEN
    1820                 :      28483 : Flx_factor_Cantor(GEN T, ulong p)
    1821                 :            : {
    1822                 :            :   GEN E, F, M, V;
    1823                 :            :   long i, j, s, l;
    1824                 :      28483 :   V = Flx_factor_squarefree(T, p);
    1825                 :      28483 :   l = lg(V);
    1826         [ +  + ]:     376306 :   for (s=0, i=1; i < l; i++)
    1827                 :     347823 :     s += !!degpol(gel(V,i));
    1828                 :      28483 :   F = cgetg(s+1, t_VEC);
    1829                 :      28483 :   E = cgetg(s+1, t_VEC);
    1830         [ +  + ]:     376306 :   for (i=1, j=1; i < l; i++)
    1831         [ +  + ]:     347823 :     if (degpol(gel(V,i)))
    1832                 :            :     {
    1833                 :      31521 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    1834                 :      31521 :       gel(F, j) = Fj;
    1835                 :      31521 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1836                 :      31521 :       j++;
    1837                 :            :     }
    1838                 :      28483 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1839                 :      28483 :   return sort_factor_pol(M, cmpGuGu);
    1840                 :            : }
    1841                 :            : 
    1842                 :            : static GEN
    1843                 :     283445 : Flx_simplefact_Cantor(GEN T, ulong p)
    1844                 :            : {
    1845                 :            :   GEN E, F, M, V;
    1846                 :            :   long i, j, s, l;
    1847                 :     283445 :   V = Flx_factor_squarefree(get_Flx_mod(T), p);
    1848                 :     283445 :   l = lg(V);
    1849         [ +  + ]:    2088707 :   for (s=0, i=1; i < l; i++)
    1850                 :    1805262 :     s += !!degpol(gel(V,i));
    1851                 :     283445 :   F = cgetg(s+1, t_VEC);
    1852                 :     283445 :   E = cgetg(s+1, t_VEC);
    1853         [ +  + ]:    2088707 :   for (i=1, j=1; i < l; i++)
    1854         [ +  + ]:    1805262 :     if (degpol(gel(V,i)))
    1855                 :            :     {
    1856                 :     284334 :       GEN Fj = Flx_simplefact_Shoup(gel(V,i), p);
    1857                 :     284334 :       gel(F, j) = Fj;
    1858                 :     284334 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1859                 :     284334 :       j++;
    1860                 :            :     }
    1861                 :     283445 :   M = mkvec2(shallowconcat1(F), shallowconcat1(E));
    1862                 :     283445 :   return sort_factor(M, (void*)&cmpGuGu, cmp_nodata);
    1863                 :            : }
    1864                 :            : 
    1865                 :            : static int
    1866                 :        581 : Flx_isirred_Cantor(GEN Tp, ulong p)
    1867                 :            : {
    1868                 :        581 :   pari_sp av = avma;
    1869                 :            :   pari_timer ti;
    1870                 :            :   long n, d;
    1871                 :        581 :   GEN T = get_Flx_mod(Tp);
    1872                 :        581 :   GEN dT = Flx_deriv(T, p);
    1873                 :            :   GEN XP, D;
    1874         [ +  + ]:        581 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1875                 :        441 :   n = get_Flx_degree(T);
    1876                 :        441 :   T = Flx_get_red(Tp, p);
    1877         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1878                 :        441 :   XP = Flx_Frobenius(T, p);
    1879         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1880                 :        441 :   D = Flx_ddf(T, XP, p);
    1881         [ -  + ]:        441 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1882                 :        441 :   d = degpol(gel(D, n));
    1883                 :        581 :   avma = av; return d==n;
    1884                 :            : }
    1885                 :            : 
    1886                 :            : static GEN
    1887                 :     398231 : Flx_factcantor_i(GEN f, ulong pp, long flag)
    1888                 :            : {
    1889                 :            :   long d;
    1890         [ +  + ]:     398231 :   if (pp==2) { /*We need to handle 2 specially */
    1891                 :       6517 :     GEN F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    1892         [ +  + ]:       6517 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    1893                 :       6517 :     return F;
    1894                 :            :   }
    1895                 :     391714 :   d = degpol(f);
    1896         [ +  + ]:     391714 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    1897      [ +  +  + ]:     312509 :   switch(flag)
    1898                 :            :   {
    1899                 :      28483 :     default: return Flx_factor_Cantor(f, pp);
    1900                 :     283445 :     case 1: return Flx_simplefact_Cantor(f, pp);
    1901         [ +  + ]:     398231 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    1902                 :            :   }
    1903                 :            : }
    1904                 :            : 
    1905                 :            : GEN
    1906                 :       7266 : Flx_factcantor(GEN f, ulong p, long flag)
    1907                 :            : {
    1908                 :       7266 :   pari_sp av = avma;
    1909                 :       7266 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,flag);
    1910         [ -  + ]:       7266 :   if (flag == 2) { avma = av; return z; }
    1911                 :       7266 :   return gerepilecopy(av, z);
    1912                 :            : }
    1913                 :            : 
    1914                 :            : GEN
    1915                 :     363630 : Flx_degfact(GEN f, ulong p)
    1916                 :            : {
    1917                 :     363630 :   pari_sp av = avma;
    1918                 :     363630 :   GEN z = Flx_factcantor_i(Flx_normalize(f,p),p,1);
    1919                 :     363630 :   return gerepilecopy(av, z);
    1920                 :            : }
    1921                 :            : 
    1922                 :            : /* T must be squarefree mod p*/
    1923                 :            : GEN
    1924                 :      55836 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    1925                 :            : {
    1926                 :            :   GEN XP, D;
    1927                 :            :   pari_timer ti;
    1928                 :      55836 :   long i, s, n = get_Flx_degree(T);
    1929                 :      55836 :   GEN V = const_vecsmall(n, 0);
    1930                 :      55836 :   pari_sp av = avma;
    1931                 :      55836 :   T = Flx_get_red(T, p);
    1932         [ -  + ]:      55836 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1933                 :      55836 :   XP = Flx_Frobenius(T, p);
    1934         [ -  + ]:      55836 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    1935                 :      55836 :   D = Flx_ddf(T, XP, p);
    1936         [ -  + ]:      55836 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf");
    1937         [ +  + ]:     731376 :   for (i = 1, s = 0; i <= n; i++)
    1938                 :            :   {
    1939                 :     675540 :     V[i] = degpol(gel(D,i))/i;
    1940                 :     675540 :     s += V[i];
    1941                 :            :   }
    1942                 :      55836 :   *nb = s;
    1943                 :      55836 :   avma = av; return V;
    1944                 :            : }
    1945                 :            : 
    1946                 :            : long
    1947                 :      52017 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    1948                 :            : {
    1949                 :      52017 :   pari_sp av = avma;
    1950                 :      52017 :   GEN ddf = Flx_ddf(T, XP, p);
    1951                 :      52017 :   long l = lg(ddf), i, s=0;
    1952         [ +  + ]:     613683 :   for(i = 1; i < l; i++)
    1953                 :     561666 :     s += degpol(gel(ddf,i))/i;
    1954                 :      52017 :   avma = av; return s;
    1955                 :            : }
    1956                 :            : 
    1957                 :            : /* T must be squarefree mod p*/
    1958                 :            : long
    1959                 :      52017 : Flx_nbfact(GEN T, ulong p)
    1960                 :            : {
    1961                 :      52017 :   pari_sp av = avma;
    1962                 :      52017 :   GEN XP = Flx_Frobenius(T, p);
    1963                 :      52017 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    1964                 :      52017 :   avma = av; return n;
    1965                 :            : }
    1966                 :            : 
    1967                 :            : int
    1968                 :        581 : Flx_is_irred(GEN f, ulong p) { return !!Flx_factcantor_i(f,p,2); }
    1969                 :            : 
    1970                 :            : /* factor f (FpX or Flx) mod pp.
    1971                 :            :  * flag = 1: return the degrees, not the factors
    1972                 :            :  * flag = 2: return NULL if f is not irreducible.
    1973                 :            :  * Not gerepile-safe */
    1974                 :            : static GEN
    1975                 :         98 : factcantor_i(GEN f, GEN pp, long flag)
    1976                 :            : {
    1977         [ +  + ]:         98 :   if (typ(f) == t_VECSMALL)
    1978                 :            :   { /* lgefint(pp) = 3 */
    1979                 :            :     GEN F;
    1980                 :         63 :     ulong p = pp[2];
    1981         [ +  + ]:         63 :     if (p==2) {
    1982                 :         35 :       F = F2x_factcantor_i(Flx_to_F2x(f),flag);
    1983         [ +  + ]:         35 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    1984                 :            :     } else {
    1985                 :         28 :       F = Flx_factcantor_i(f,p,flag);
    1986         [ +  + ]:         28 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    1987                 :            :     }
    1988                 :         63 :     return F;
    1989                 :            :   }
    1990                 :         98 :   return FpX_factcantor_i(f, pp, flag);
    1991                 :            : }
    1992                 :            : GEN
    1993                 :          0 : FpX_factcantor(GEN f, GEN pp, long flag)
    1994                 :            : {
    1995                 :          0 :   pari_sp av = avma;
    1996                 :            :   GEN z;
    1997                 :          0 :   ZX_factmod_init(&f,pp);
    1998                 :          0 :   z = factcantor_i(f,pp,flag);
    1999         [ #  # ]:          0 :   if (flag == 2) { avma = av; return z; }
    2000                 :          0 :   return gerepilecopy(av, z);
    2001                 :            : }
    2002                 :            : 
    2003                 :            : static GEN
    2004                 :        182 : factmod_aux(GEN f, GEN p, GEN (*Factor)(GEN,GEN,long), long flag)
    2005                 :            : {
    2006                 :        182 :   pari_sp av = avma;
    2007                 :            :   long j, nbfact;
    2008                 :            :   GEN z, y, t, E, u, v;
    2009                 :            : 
    2010                 :        182 :   factmod_init(&f, p);
    2011      [ +  +  + ]:        182 :   switch(lg(f))
    2012                 :            :   {
    2013                 :         14 :     case 3: avma = av; return trivial_fact();
    2014                 :         14 :     case 2: return gerepileupto(av, zero_fact_intmod(f, p));
    2015                 :            :   }
    2016                 :        154 :   z = Factor(f,p,flag); t = gel(z,1); E = gel(z,2);
    2017                 :        154 :   y = cgetg(3, t_MAT); nbfact = lg(t);
    2018                 :        154 :   u = cgetg(nbfact,t_COL); gel(y,1) = u;
    2019                 :        154 :   v = cgetg(nbfact,t_COL); gel(y,2) = v;
    2020         [ +  + ]:        154 :   if (flag)
    2021         [ +  + ]:         84 :     for (j=1; j<nbfact; j++)
    2022                 :            :     {
    2023                 :         56 :       gel(u,j) = utoi(uel(t,j));
    2024                 :         56 :       gel(v,j) = utoi(uel(E,j));
    2025                 :            :     }
    2026                 :            :   else
    2027         [ +  + ]:       7896 :     for (j=1; j<nbfact; j++)
    2028                 :            :     {
    2029                 :       7770 :       gel(u,j) = FpX_to_mod(gel(t,j), p);
    2030                 :       7770 :       gel(v,j) = utoi(uel(E,j));
    2031                 :            :     }
    2032                 :        182 :   return gerepileupto(av, y);
    2033                 :            : }
    2034                 :            : GEN
    2035                 :         98 : factcantor0(GEN f, GEN p, long flag)
    2036                 :         98 : { return factmod_aux(f, p, &factcantor_i, flag); }
    2037                 :            : GEN
    2038                 :         84 : factmod(GEN f, GEN p)
    2039                 :         84 : { return factmod_aux(f, p, &FpX_Berlekamp_i, 0); }
    2040                 :            : 
    2041                 :            : /* Use this function when you think f is reducible, and that there are lots of
    2042                 :            :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2043                 :            : int
    2044                 :         14 : FpX_is_irred(GEN f, GEN p) {
    2045                 :         14 :   ZX_factmod_init(&f,p);
    2046                 :         14 :   return !!factcantor_i(f,p,2);
    2047                 :            : }
    2048                 :            : GEN
    2049                 :          0 : FpX_degfact(GEN f, GEN p) {
    2050                 :          0 :   pari_sp av = avma;
    2051                 :            :   GEN z;
    2052                 :          0 :   ZX_factmod_init(&f,p);
    2053                 :          0 :   z = factcantor_i(f,p,1);
    2054                 :          0 :   return gerepilecopy(av, z);
    2055                 :            : }
    2056                 :            : GEN
    2057                 :         70 : factcantor(GEN f, GEN p) { return factcantor0(f,p,0); }
    2058                 :            : GEN
    2059                 :         28 : simplefactmod(GEN f, GEN p) { return factcantor0(f,p,1); }
    2060                 :            : 
    2061                 :            : /* set x <-- x + c*y mod p */
    2062                 :            : /* x is not required to be normalized.*/
    2063                 :            : static void
    2064                 :     601004 : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2065                 :            : {
    2066                 :            :   long i, lx, ly;
    2067                 :     601004 :   ulong *x=(ulong *)gx;
    2068                 :     601004 :   ulong *y=(ulong *)gy;
    2069         [ +  + ]:    1182662 :   if (!c) return;
    2070                 :     581658 :   lx = lg(gx);
    2071                 :     581658 :   ly = lg(gy);
    2072         [ -  + ]:     581658 :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2073         [ +  + ]:     581658 :   if (SMALL_ULONG(p))
    2074         [ +  + ]:    4066975 :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2075                 :            :   else
    2076         [ +  + ]:      17095 :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2077                 :            : }
    2078                 :            : 
    2079                 :            : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    2080                 :            : /* assume x1 != 0 */
    2081                 :            : static GEN
    2082                 :     118970 : deg1_Flx(ulong x1, ulong x0, ulong sv)
    2083                 :            : {
    2084                 :     118970 :   return mkvecsmall3(sv, x0, x1);
    2085                 :            : }
    2086                 :            : 
    2087                 :            : long
    2088                 :      43867 : F2x_split_Berlekamp(GEN *t)
    2089                 :            : {
    2090                 :      43867 :   GEN u = *t, a, b, vker;
    2091                 :      43867 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    2092                 :            : 
    2093         [ +  + ]:      43867 :   if (du == 1) return 1;
    2094         [ +  + ]:      28112 :   if (du == 2)
    2095                 :            :   {
    2096         [ -  + ]:       7259 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    2097                 :            :     {
    2098                 :          0 :       t[0] = mkvecsmall2(sv, 2);
    2099                 :          0 :       t[1] = mkvecsmall2(sv, 3);
    2100                 :          0 :       return 2;
    2101                 :            :     }
    2102                 :       7259 :     return 1;
    2103                 :            :   }
    2104                 :            : 
    2105                 :      20853 :   vker = F2x_Berlekamp_ker(u);
    2106                 :      20853 :   lb = lgcols(vker);
    2107                 :      20853 :   d = lg(vker)-1;
    2108                 :      20853 :   ir = 0;
    2109                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2110         [ +  + ]:      26411 :   for (L=1; L<d; )
    2111                 :            :   {
    2112                 :            :     GEN pol;
    2113         [ +  + ]:       5558 :     if (d == 2)
    2114                 :       5271 :       pol = F2v_to_F2x(gel(vker,2), sv);
    2115                 :            :     else
    2116                 :            :     {
    2117                 :        287 :       GEN v = zero_zv(lb);
    2118                 :        287 :       v[1] = du;
    2119                 :        287 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    2120         [ +  + ]:       1113 :       for (i=2; i<=d; i++)
    2121         [ +  + ]:        826 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    2122                 :        287 :       pol = F2v_to_F2x(v, sv);
    2123                 :            :     }
    2124 [ +  + ][ +  + ]:      11508 :     for (i=ir; i<L && L<d; i++)
    2125                 :            :     {
    2126                 :       5950 :       a = t[i]; la = F2x_degree(a);
    2127 [ +  + ][ +  + ]:       5950 :       if (la == 1) { set_irred(i); }
    2128         [ +  + ]:       5929 :       else if (la == 2)
    2129                 :            :       {
    2130         [ -  + ]:          7 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    2131                 :            :         {
    2132                 :          0 :           t[i] = mkvecsmall2(sv, 2);
    2133                 :          0 :           t[L] = mkvecsmall2(sv, 3); L++;
    2134                 :            :         }
    2135         [ +  - ]:          7 :         set_irred(i);
    2136                 :            :       }
    2137                 :            :       else
    2138                 :            :       {
    2139                 :       5922 :         pari_sp av = avma;
    2140                 :            :         long lb;
    2141                 :       5922 :         b = F2x_rem(pol, a);
    2142         [ +  + ]:       5922 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    2143                 :       5488 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    2144 [ +  - ][ +  - ]:       5488 :         if (lb && lb < la)
    2145                 :            :         {
    2146                 :       5488 :           t[L] = F2x_div(a,b);
    2147                 :       5488 :           t[i]= b; L++;
    2148                 :            :         }
    2149                 :          0 :         else avma = av;
    2150                 :            :       }
    2151                 :            :     }
    2152                 :            :   }
    2153                 :      43867 :   return d;
    2154                 :            : }
    2155                 :            : 
    2156                 :            : /* p != 2 */
    2157                 :            : long
    2158                 :     232802 : Flx_split_Berlekamp(GEN *t, ulong p)
    2159                 :            : {
    2160                 :     232802 :   GEN u = *t, a,b,vker;
    2161                 :     232802 :   long d, i, ir, L, la, lb, sv = u[1];
    2162                 :     232802 :   long l = lg(u);
    2163                 :            :   ulong po2;
    2164                 :            : 
    2165         [ -  + ]:     232802 :   if (p == 2)
    2166                 :            :   {
    2167                 :          0 :     *t = Flx_to_F2x(*t);
    2168                 :          0 :     d = F2x_split_Berlekamp(t);
    2169         [ #  # ]:          0 :     for (i = 1; i <= d; i++) t[i] = F2x_to_Flx(t[i]);
    2170                 :          0 :     return d;
    2171                 :            :   }
    2172                 :     232802 :   la = degpol(u);
    2173         [ +  + ]:     232802 :   if (la == 1) return 1;
    2174         [ +  + ]:     217494 :   if (la == 2)
    2175                 :            :   {
    2176                 :      11933 :     ulong r = Flx_quad_root(u,p,1);
    2177         [ +  + ]:      11933 :     if (r != p)
    2178                 :            :     {
    2179                 :       5439 :       t[0] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(u,r,p);
    2180                 :       5439 :       t[1] = deg1_Flx(1, p - r, sv);
    2181                 :       5439 :       return 2;
    2182                 :            :     }
    2183                 :       6494 :     return 1;
    2184                 :            :   }
    2185                 :            : 
    2186                 :     205561 :   vker = Flx_Berlekamp_ker(u,p);
    2187                 :     205561 :   vker = Flm_to_FlxV(vker, sv);
    2188                 :     205561 :   d = lg(vker)-1;
    2189                 :     205561 :   po2 = p >> 1; /* (p-1) / 2 */
    2190                 :     205561 :   ir = 0;
    2191                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2192         [ +  + ]:     507899 :   for (L=1; L<d; )
    2193                 :            :   {
    2194                 :     302338 :     GEN pol = zero_zv(l-2);
    2195                 :     302338 :     pol[1] = sv;
    2196                 :     302338 :     pol[2] = random_Fl(p); /*Assume vker[1]=1*/
    2197         [ +  + ]:     903342 :     for (i=2; i<=d; i++)
    2198                 :     601004 :       Flx_addmul_inplace(pol, gel(vker,i), random_Fl(p), p);
    2199                 :     302338 :     (void)Flx_renormalize(pol,l-1);
    2200                 :            : 
    2201 [ +  + ][ +  + ]:     729245 :     for (i=ir; i<L && L<d; i++)
    2202                 :            :     {
    2203                 :     426907 :       a = t[i]; la = degpol(a);
    2204 [ +  + ][ +  + ]:     426907 :       if (la == 1) { set_irred(i); }
    2205         [ +  + ]:     384961 :       else if (la == 2)
    2206                 :            :       {
    2207                 :      66149 :         ulong r = Flx_quad_root(a,p,1);
    2208         [ +  + ]:      66149 :         if (r != p)
    2209                 :            :         {
    2210                 :      54046 :           t[i] = deg1_Flx(1, p - r, sv); r = Flx_otherroot(a,r,p);
    2211                 :      54046 :           t[L] = deg1_Flx(1, p - r, sv); L++;
    2212                 :            :         }
    2213         [ +  + ]:      66149 :         set_irred(i);
    2214                 :            :       }
    2215                 :            :       else
    2216                 :            :       {
    2217                 :     318812 :         pari_sp av = avma;
    2218                 :     318812 :         b = Flx_rem(pol, a, p);
    2219         [ +  + ]:     318812 :         if (degpol(b) <= 0) { avma=av; continue; }
    2220                 :     300607 :         b = Flx_Fl_add(Flxq_powu(b,po2, a,p), p-1, p);
    2221                 :     300607 :         b = Flx_gcd(a,b, p); lb = degpol(b);
    2222 [ +  + ][ +  + ]:     300607 :         if (lb && lb < la)
    2223                 :            :         {
    2224                 :     173240 :           b = Flx_normalize(b, p);
    2225                 :     173240 :           t[L] = Flx_div(a,b,p);
    2226                 :     173240 :           t[i]= b; L++;
    2227                 :            :         }
    2228                 :     127367 :         else avma = av;
    2229                 :            :       }
    2230                 :            :     }
    2231                 :            :   }
    2232                 :     232802 :   return d;
    2233                 :            : }
    2234                 :            : 
    2235                 :            : long
    2236                 :        692 : FpX_split_Berlekamp(GEN *t, GEN p)
    2237                 :            : {
    2238                 :        692 :   GEN u = *t, a,b,po2,vker;
    2239                 :        692 :   long d, i, ir, L, la, lb, vu = varn(u);
    2240         [ +  + ]:        692 :   if (lgefint(p) == 3)
    2241                 :            :   {
    2242                 :        630 :     ulong up = p[2];
    2243         [ -  + ]:        630 :     if (up == 2)
    2244                 :            :     {
    2245                 :          0 :       *t = ZX_to_F2x(*t);
    2246                 :          0 :       d = F2x_split_Berlekamp(t);
    2247         [ #  # ]:          0 :       for (i = 0; i < d; i++) t[i] = F2x_to_ZX(t[i]);
    2248                 :            :     }
    2249                 :            :     else
    2250                 :            :     {
    2251                 :        630 :       *t = ZX_to_Flx(*t, up);
    2252                 :        630 :       d = Flx_split_Berlekamp(t, up);
    2253         [ +  + ]:       3850 :       for (i = 0; i < d; i++) t[i] = Flx_to_ZX(t[i]);
    2254                 :            :     }
    2255                 :        630 :     return d;
    2256                 :            :   }
    2257                 :         62 :   la = degpol(u);
    2258         [ +  + ]:         62 :   if (la == 1) return 1;
    2259         [ +  + ]:         57 :   if (la == 2)
    2260                 :            :   {
    2261                 :          2 :     GEN r = FpX_quad_root(u,p,1);
    2262         [ -  + ]:          2 :     if (r)
    2263                 :            :     {
    2264                 :          0 :       t[0] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(u,r,p);
    2265                 :          0 :       t[1] = deg1pol_shallow(gen_1, subii(p,r), vu);
    2266                 :          0 :       return 2;
    2267                 :            :     }
    2268                 :          2 :     return 1;
    2269                 :            :   }
    2270                 :         55 :   vker = FpX_Berlekamp_ker(u,p);
    2271                 :         55 :   vker = RgM_to_RgXV(vker,vu);
    2272                 :         55 :   d = lg(vker)-1;
    2273                 :         55 :   po2 = shifti(p, -1); /* (p-1) / 2 */
    2274                 :         55 :   ir = 0;
    2275                 :            :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    2276         [ +  + ]:        107 :   for (L=1; L<d; )
    2277                 :            :   {
    2278                 :         52 :     GEN pol = scalar_ZX_shallow(randomi(p), vu);
    2279         [ +  + ]:        208 :     for (i=2; i<=d; i++)
    2280                 :        156 :       pol = ZX_add(pol, ZX_Z_mul(gel(vker,i), randomi(p)));
    2281                 :         52 :     pol = FpX_red(pol,p);
    2282 [ +  + ][ +  + ]:        157 :     for (i=ir; i<L && L<d; i++)
    2283                 :            :     {
    2284                 :        105 :       a = t[i]; la = degpol(a);
    2285 [ +  + ][ +  - ]:        105 :       if (la == 1) { set_irred(i); }
    2286         [ +  + ]:         98 :       else if (la == 2)
    2287                 :            :       {
    2288                 :         14 :         GEN r = FpX_quad_root(a,p,1);
    2289         [ +  - ]:         14 :         if (r)
    2290                 :            :         {
    2291                 :         14 :           t[i] = deg1pol_shallow(gen_1, subii(p,r), vu); r = FpX_otherroot(a,r,p);
    2292                 :         14 :           t[L] = deg1pol_shallow(gen_1, subii(p,r), vu); L++;
    2293                 :            :         }
    2294         [ +  + ]:         14 :         set_irred(i);
    2295                 :            :       }
    2296                 :            :       else
    2297                 :            :       {
    2298                 :         84 :         pari_sp av = avma;
    2299                 :         84 :         b = FpX_rem(pol, a, p);
    2300         [ +  + ]:         84 :         if (degpol(b) <= 0) { avma=av; continue; }
    2301                 :         63 :         b = FpX_Fp_sub_shallow(FpXQ_pow(b,po2, a,p), gen_1, p);
    2302                 :         63 :         b = FpX_gcd(a,b, p); lb = degpol(b);
    2303 [ +  + ][ +  + ]:         63 :         if (lb && lb < la)
    2304                 :            :         {
    2305                 :         43 :           b = FpX_normalize(b, p);
    2306                 :         43 :           t[L] = FpX_div(a,b,p);
    2307                 :         43 :           t[i]= b; L++;
    2308                 :            :         }
    2309                 :         20 :         else avma = av;
    2310                 :            :       }
    2311                 :            :     }
    2312                 :            :   }
    2313                 :        692 :   return d;
    2314                 :            : }
    2315                 :            : 
    2316                 :            : static GEN
    2317                 :      89528 : F2x_Berlekamp_i(GEN f, long flag)
    2318                 :            : {
    2319                 :      89528 :   long e, nbfact, val, d = F2x_degree(f);
    2320                 :            :   ulong k, j;
    2321                 :            :   GEN y, E, f2, g1, t;
    2322                 :            : 
    2323         [ +  + ]:      89528 :   if (d <= 2) return F2x_factor_deg2(f, d, flag);
    2324                 :            : 
    2325                 :            :   /* to hold factors and exponents */
    2326         [ -  + ]:      45365 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2327                 :      45365 :   E = cgetg(d+1,t_VECSMALL);
    2328                 :      45365 :   val = F2x_valrem(f, &f);
    2329                 :      45365 :   e = nbfact = 1;
    2330         [ +  + ]:      45365 :   if (val) {
    2331 [ -  + ][ #  # ]:       9343 :     if (flag == 2 && val > 1) return NULL;
    2332         [ -  + ]:       9343 :     if (flag == 1)
    2333                 :          0 :       t[1] = 1;
    2334                 :            :     else
    2335                 :       9343 :       gel(t,1) = polx_F2x(f[1]);
    2336                 :       9343 :     E[1] = val; nbfact++;
    2337                 :            :   }
    2338                 :            : 
    2339                 :            :   for(;;)
    2340                 :            :   {
    2341                 :      82171 :     f2 = F2x_gcd(f,F2x_deriv(f));
    2342 [ -  + ][ #  # ]:      82171 :     if (flag == 2 && F2x_degree(f2)) return NULL;
    2343         [ +  + ]:      82171 :     g1 = F2x_degree(f2)? F2x_div(f,f2): f; /* squarefree */
    2344                 :      82171 :     k = 0;
    2345         [ +  + ]:     137406 :     while (F2x_degree(g1)>0)
    2346                 :            :     {
    2347                 :            :       GEN u;
    2348         [ +  + ]:      55235 :       k++; if (k%2 == 0) { k++; f2 = F2x_div(f2,g1); }
    2349                 :      55235 :       u = g1; /* deg(u) > 0 */
    2350         [ +  + ]:      55235 :       if (!F2x_degree(f2)) g1 = pol1_F2x(0); /* only its degree (= 0) matters */
    2351                 :            :       else
    2352                 :            :       {
    2353                 :            :         long dg1;
    2354                 :      14308 :         g1 = F2x_gcd(f2,g1);
    2355                 :      14308 :         dg1 = F2x_degree(g1);
    2356         [ +  + ]:      14308 :         if (dg1)
    2357                 :            :         {
    2358                 :      12495 :           f2= F2x_div(f2,g1);
    2359         [ +  + ]:      12495 :           if (F2x_degree(u) == dg1) continue;
    2360                 :       1127 :           u = F2x_div( u,g1);
    2361                 :            :         }
    2362                 :            :       }
    2363                 :            :       /* u is square-free (product of irred. of multiplicity e * k) */
    2364                 :      43867 :       gel(t,nbfact) = u;
    2365                 :      43867 :       d = F2x_split_Berlekamp(&gel(t,nbfact));
    2366 [ -  + ][ #  # ]:      43867 :       if (flag == 2 && d != 1) return NULL;
    2367         [ -  + ]:      43867 :       if (flag == 1)
    2368         [ #  # ]:          0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = F2x_degree(gel(t,nbfact+j));
    2369         [ +  + ]:      93222 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2370                 :      43867 :       nbfact += d;
    2371                 :            :     }
    2372         [ +  + ]:      82171 :     j = F2x_degree(f2); if (!j) break;
    2373                 :      36806 :     e *= 2; f = F2x_deflate(f2, 2);
    2374                 :      36806 :   }
    2375         [ -  + ]:      45365 :   if (flag == 2) return gen_1; /* irreducible */
    2376                 :      45365 :   setlg(t, nbfact);
    2377                 :      45365 :   setlg(E, nbfact); y = mkvec2(t,E);
    2378                 :      89528 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2379         [ -  + ]:      45365 :               : sort_factor_pol(y, cmpGuGu);
    2380                 :            : }
    2381                 :            : 
    2382                 :            : static GEN
    2383                 :     496238 : Flx_Berlekamp_i(GEN f, ulong p, long flag)
    2384                 :            : {
    2385                 :     496238 :   long e, nbfact, val, d = degpol(f);
    2386                 :            :   ulong k, j;
    2387                 :            :   GEN y, E, f2, g1, t;
    2388                 :            : 
    2389         [ +  + ]:     496238 :   if (p == 2)
    2390                 :            :   {
    2391                 :      10962 :     GEN F = F2x_Berlekamp_i(Flx_to_F2x(f),flag);
    2392         [ +  - ]:      10962 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2393                 :      10962 :     return F;
    2394                 :            :   }
    2395         [ +  + ]:     485276 :   if (d <= 2) return Flx_factor_deg2(f,p,d,flag);
    2396                 :            : 
    2397                 :            :   /* to hold factors and exponents */
    2398         [ -  + ]:     228746 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2399                 :     228746 :   E = cgetg(d+1,t_VECSMALL);
    2400                 :     228746 :   val = Flx_valrem(f, &f);
    2401                 :     228746 :   e = nbfact = 1;
    2402         [ +  + ]:     228746 :   if (val) {
    2403 [ -  + ][ #  # ]:       7330 :     if (flag == 2 && val > 1) return NULL;
    2404         [ -  + ]:       7330 :     if (flag == 1)
    2405                 :          0 :       t[1] = 1;
    2406                 :            :     else
    2407                 :       7330 :       gel(t,1) = polx_Flx(f[1]);
    2408                 :       7330 :     E[1] = val; nbfact++;
    2409                 :            :   }
    2410                 :            : 
    2411                 :            :   for(;;)
    2412                 :            :   {
    2413                 :     234374 :     f2 = Flx_gcd(f,Flx_deriv(f,p), p);
    2414 [ -  + ][ #  # ]:     234374 :     if (flag == 2 && degpol(f2)) return NULL;
    2415         [ +  + ]:     234374 :     g1 = degpol(f2)? Flx_div(f,f2,p): f; /* squarefree */
    2416                 :     234374 :     k = 0;
    2417         [ +  + ]:     507600 :     while (degpol(g1)>0)
    2418                 :            :     {
    2419                 :            :       GEN u;
    2420         [ +  + ]:     273226 :       k++; if (k%p == 0) { k++; f2 = Flx_div(f2,g1,p); }
    2421                 :     273226 :       u = g1; /* deg(u) > 0 */
    2422         [ +  + ]:     273226 :       if (!degpol(f2)) g1 = pol1_Flx(0); /* only its degree (= 0) matters */
    2423                 :            :       else
    2424                 :            :       {
    2425                 :      46838 :         g1 = Flx_gcd(f2,g1, p);
    2426         [ +  + ]:      46838 :         if (degpol(g1))
    2427                 :            :         {
    2428                 :      46691 :           f2= Flx_div(f2,g1,p);
    2429         [ +  + ]:      46691 :           if (lg(u) == lg(g1)) continue;
    2430                 :       5637 :           u = Flx_div( u,g1,p);
    2431                 :            :         }
    2432                 :            :       }
    2433                 :            :       /* u is square-free (product of irred. of multiplicity e * k) */
    2434                 :     232172 :       u = Flx_normalize(u,p);
    2435                 :     232172 :       gel(t,nbfact) = u;
    2436                 :     232172 :       d = Flx_split_Berlekamp(&gel(t,nbfact), p);
    2437 [ -  + ][ #  # ]:     232172 :       if (flag == 2 && d != 1) return NULL;
    2438         [ -  + ]:     232172 :       if (flag == 1)
    2439         [ #  # ]:          0 :         for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2440         [ +  + ]:     694479 :       for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2441                 :     232172 :       nbfact += d;
    2442                 :            :     }
    2443         [ -  + ]:     234374 :     if (!p) break;
    2444         [ +  + ]:     234374 :     j = degpol(f2); if (!j) break;
    2445         [ -  + ]:       5628 :     if (j % p) pari_err_PRIME("factmod",utoi(p));
    2446                 :       5628 :     e *= p; f = Flx_deflate(f2, p);
    2447                 :       5628 :   }
    2448         [ -  + ]:     228746 :   if (flag == 2) return gen_1; /* irreducible */
    2449                 :     228746 :   setlg(t, nbfact);
    2450                 :     228746 :   setlg(E, nbfact); y = mkvec2(t,E);
    2451                 :     496238 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2452         [ -  + ]:     228746 :               : sort_factor_pol(y, cmpGuGu);
    2453                 :            : }
    2454                 :            : 
    2455                 :            : /* f an FpX or an Flx */
    2456                 :            : static GEN
    2457                 :     540503 : FpX_Berlekamp_i(GEN f, GEN pp, long flag)
    2458                 :            : {
    2459                 :     540503 :   long e, nbfact, val, d = degpol(f);
    2460                 :            :   ulong k, j;
    2461                 :            :   GEN y, E, f2, g1, t;
    2462                 :            : 
    2463         [ +  + ]:     540503 :   if (typ(f) == t_VECSMALL)
    2464                 :            :   {/* lgefint(pp) == 3 */
    2465                 :     540378 :     ulong p = pp[2];
    2466                 :            :     GEN F;
    2467         [ +  + ]:     540378 :     if (p == 2) {
    2468                 :      78552 :       F = F2x_Berlekamp_i(Flx_to_F2x(f), flag);
    2469         [ +  - ]:      78552 :       if (flag==0) F2xV_to_ZXV_inplace(gel(F,1));
    2470                 :            :     } else {
    2471                 :     461826 :       F = Flx_Berlekamp_i(f, p, flag);
    2472         [ +  - ]:     461826 :       if (flag==0) FlxV_to_ZXV_inplace(gel(F,1));
    2473                 :            :     }
    2474                 :     540378 :     return F;
    2475                 :            :   }
    2476                 :            :   /* p is large (and odd) */
    2477         [ +  + ]:        125 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    2478                 :            : 
    2479                 :            :   /* to hold factors and exponents */
    2480         [ -  + ]:         58 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    2481                 :         58 :   E = cgetg(d+1,t_VECSMALL);
    2482                 :         58 :   val = ZX_valrem(f, &f);
    2483                 :         58 :   e = nbfact = 1;
    2484         [ +  + ]:         58 :   if (val) {
    2485 [ -  + ][ #  # ]:          3 :     if (flag == 2 && val > 1) return NULL;
    2486         [ -  + ]:          3 :     if (flag == 1)
    2487                 :          0 :       t[1] = 1;
    2488                 :            :     else
    2489                 :          3 :       gel(t,1) = pol_x(varn(f));
    2490                 :          3 :     E[1] = val; nbfact++;
    2491                 :            :   }
    2492                 :            : 
    2493                 :         58 :   f2 = FpX_gcd(f,ZX_deriv(f), pp);
    2494 [ -  + ][ #  # ]:         58 :   if (flag == 2 && degpol(f2)) return NULL;
    2495         [ +  + ]:         58 :   g1 = degpol(f2)? FpX_div(f,f2,pp): f; /* squarefree */
    2496                 :         58 :   k = 0;
    2497         [ +  + ]:        121 :   while (degpol(g1)>0)
    2498                 :            :   {
    2499                 :            :     GEN u;
    2500                 :         63 :     k++;
    2501                 :         63 :     u = g1; /* deg(u) > 0 */
    2502         [ +  + ]:         63 :     if (!degpol(f2)) g1 = pol_1(0); /* only its degree (= 0) matters */
    2503                 :            :     else
    2504                 :            :     {
    2505                 :          6 :       g1 = FpX_gcd(f2,g1, pp);
    2506         [ +  - ]:          6 :       if (degpol(g1))
    2507                 :            :       {
    2508                 :          6 :         f2= FpX_div(f2,g1,pp);
    2509         [ +  + ]:          6 :         if (lg(u) == lg(g1)) continue;
    2510                 :          5 :         u = FpX_div( u,g1,pp);
    2511                 :            :       }
    2512                 :            :     }
    2513                 :            :     /* u is square-free (product of irred. of multiplicity e * k) */
    2514                 :         62 :     u = FpX_normalize(u,pp);
    2515                 :         62 :     gel(t,nbfact) = u;
    2516                 :         62 :     d = FpX_split_Berlekamp(&gel(t,nbfact), pp);
    2517 [ -  + ][ #  # ]:         62 :     if (flag == 2 && d != 1) return NULL;
    2518         [ -  + ]:         62 :     if (flag == 1)
    2519         [ #  # ]:          0 :       for (j=0; j<(ulong)d; j++) t[nbfact+j] = degpol(gel(t,nbfact+j));
    2520         [ +  + ]:        181 :     for (j=0; j<(ulong)d; j++) E[nbfact+j] = e*k;
    2521                 :         62 :     nbfact += d;
    2522                 :            :   }
    2523         [ -  + ]:         58 :   if (flag == 2) return gen_1; /* irreducible */
    2524                 :         58 :   setlg(t, nbfact);
    2525                 :         58 :   setlg(E, nbfact); y = mkvec2(t,E);
    2526                 :     540503 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    2527         [ -  + ]:         58 :               : sort_factor_pol(y, cmpii);
    2528                 :            : }
    2529                 :            : GEN
    2530                 :     540433 : FpX_factor(GEN f, GEN p)
    2531                 :            : {
    2532                 :     540433 :   pari_sp av = avma;
    2533                 :     540433 :   ZX_factmod_init(&f, p);
    2534                 :     540433 :   return gerepilecopy(av, FpX_Berlekamp_i(f, p, 0));
    2535                 :            : }
    2536                 :            : GEN
    2537                 :      61138 : Flx_factor(GEN f, ulong p)
    2538                 :            : {
    2539                 :      61138 :   pari_sp av = avma;
    2540         [ +  + ]:      61138 :   GEN F = (degpol(f)>log2(p))? Flx_factcantor_i(f,p,0): Flx_Berlekamp_i(f,p,0);
    2541                 :      61138 :   return gerepilecopy(av, F);
    2542                 :            : }
    2543                 :            : GEN
    2544                 :         14 : F2x_factor(GEN f)
    2545                 :            : {
    2546                 :         14 :   pari_sp av = avma;
    2547                 :         14 :   return gerepilecopy(av, F2x_Berlekamp_i(f, 0));
    2548                 :            : }
    2549                 :            : 
    2550                 :            : GEN
    2551                 :        105 : factormod0(GEN f, GEN p, long flag)
    2552                 :            : {
    2553      [ +  +  - ]:        105 :   switch(flag)
    2554                 :            :   {
    2555                 :         77 :     case 0: return factmod(f,p);
    2556                 :         28 :     case 1: return simplefactmod(f,p);
    2557                 :          0 :     default: pari_err_FLAG("factormod");
    2558                 :            :   }
    2559                 :        105 :   return NULL; /* not reached */
    2560                 :            : }

Generated by: LCOV version 1.9