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 to exceed 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.12.1 lcov report (development 25819-e703fe1174) Lines: 1245 1356 91.8 %
Date: 2020-09-18 06:10:04 Functions: 107 117 91.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /***********************************************************************/
      18             : /**                                                                   **/
      19             : /**               Factorisation over finite field                     **/
      20             : /**                                                                   **/
      21             : /***********************************************************************/
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : /* Replace F by a monic normalized FpX having the same factors;
      29             :  * assume p prime and *F a ZX */
      30             : static int
      31     1047744 : ZX_factmod_init(GEN *F, GEN p)
      32             : {
      33     1047744 :   if (lgefint(p) == 3)
      34             :   {
      35     1046193 :     ulong pp = p[2];
      36     1046193 :     if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
      37      875067 :     *F = ZX_to_Flx(*F, pp);
      38      875066 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      39      875059 :     return 1;
      40             :   }
      41        1551 :   *F = FpX_red(*F, p);
      42        1551 :   if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      43        1551 :   return 2;
      44             : }
      45             : static void
      46       91685 : ZX_rootmod_init(GEN *F, GEN p)
      47             : {
      48       91685 :   if (lgefint(p) == 3)
      49             :   {
      50       83714 :     ulong pp = p[2];
      51       83714 :     *F = ZX_to_Flx(*F, pp);
      52       83714 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      53             :   }
      54             :   else
      55             :   {
      56        7971 :     *F = FpX_red(*F, p);
      57        7971 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      58             :   }
      59       91685 : }
      60             : 
      61             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      62             : static GEN
      63          42 : all_roots_mod_p(ulong p, int not_0)
      64             : {
      65             :   GEN r;
      66             :   ulong i;
      67          42 :   if (not_0) {
      68          28 :     r = cgetg(p, t_VECSMALL);
      69          84 :     for (i = 1; i < p; i++) r[i] = i;
      70             :   } else {
      71          14 :     r = cgetg(p+1, t_VECSMALL);
      72          84 :     for (i = 0; i < p; i++) r[i+1] = i;
      73             :   }
      74          42 :   return r;
      75             : }
      76             : 
      77             : /* X^n - 1 */
      78             : static GEN
      79         231 : Flx_Xnm1(long sv, long n, ulong p)
      80             : {
      81         231 :   GEN t = cgetg(n+3, t_VECSMALL);
      82             :   long i;
      83         231 :   t[1] = sv;
      84         231 :   t[2] = p - 1;
      85         665 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      86         231 :   t[i] = 1; return t;
      87             : }
      88             : /* X^n + 1 */
      89             : static GEN
      90         504 : Flx_Xn1(long sv, long n, ulong p)
      91             : {
      92         504 :   GEN t = cgetg(n+3, t_VECSMALL);
      93             :   long i;
      94             :   (void) p;
      95         504 :   t[1] = sv;
      96         504 :   t[2] = 1;
      97        1456 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      98         504 :   t[i] = 1; return t;
      99             : }
     100             : 
     101             : static GEN
     102       12538 : Flx_root_mod_2(GEN f)
     103             : {
     104       12538 :   int z1, z0 = !(f[2] & 1);
     105             :   long i,n;
     106             :   GEN y;
     107             : 
     108       48151 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     109       12538 :   z1 = n & 1;
     110       12538 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     111       12538 :   if (z0) y[i++] = 0;
     112       12538 :   if (z1) y[i  ] = 1;
     113       12538 :   return y;
     114             : }
     115             : static ulong
     116          70 : Flx_oneroot_mod_2(GEN f)
     117             : {
     118             :   long i,n;
     119          70 :   if (!(f[2] & 1)) return 0;
     120         280 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     121          70 :   if (n & 1) return 1;
     122          28 :   return 2;
     123             : }
     124             : 
     125             : static GEN FpX_roots_i(GEN f, GEN p);
     126             : static GEN Flx_roots_i(GEN f, ulong p);
     127             : 
     128             : static int
     129     5460429 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     130             : 
     131             : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
     132             :  * pp is a small prime.
     133             :  * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
     134             :  * assume that f is an FpX, pp a prime and return t_INTs */
     135             : static GEN
     136       79029 : rootmod_aux(GEN f, GEN pp)
     137             : {
     138             :   GEN y;
     139       79029 :   switch(lg(f))
     140             :   {
     141          14 :     case 2: pari_err_ROOTS0("rootmod");
     142          49 :     case 3: return cgetg(1,t_COL);
     143             :   }
     144       78966 :   if (typ(f) == t_VECSMALL)
     145             :   {
     146       75856 :     ulong p = pp[2];
     147       75856 :     if (p == 2)
     148       12538 :       y = Flx_root_mod_2(f);
     149             :     else
     150             :     {
     151       63318 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     152       63318 :       y = Flx_roots_i(f, p);
     153             :     }
     154       75849 :     y = Flc_to_ZC(y);
     155             :   }
     156             :   else
     157        3110 :     y = FpX_roots_i(f, pp);
     158       78952 :   return y;
     159             : }
     160             : /* assume that f is a ZX and p a prime */
     161             : GEN
     162       79029 : FpX_roots(GEN f, GEN p)
     163             : {
     164       79029 :   pari_sp av = avma;
     165       79029 :   GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
     166       79001 :   return gerepileupto(av, y);
     167             : }
     168             : 
     169             : /* assume x reduced mod p > 2, monic. */
     170             : static int
     171          21 : FpX_quad_factortype(GEN x, GEN p)
     172             : {
     173          21 :   GEN b = gel(x,3), c = gel(x,2);
     174          21 :   GEN D = subii(sqri(b), shifti(c,2));
     175          21 :   return kronecker(D,p);
     176             : }
     177             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     178             : static GEN
     179        7586 : FpX_quad_root(GEN x, GEN p, int unknown)
     180             : {
     181        7586 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     182             : 
     183        7586 :   if (absequaliu(p, 2)) {
     184           0 :     if (!signe(b)) return c;
     185           0 :     return signe(c)? NULL: gen_1;
     186             :   }
     187        7586 :   D = subii(sqri(b), shifti(c,2));
     188        7586 :   D = remii(D,p);
     189        7586 :   if (unknown && kronecker(D,p) == -1) return NULL;
     190             : 
     191        7052 :   s = Fp_sqrt(D,p);
     192             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     193        7052 :   if (!s) return NULL;
     194        7044 :   return Fp_halve(Fp_sub(s,b, p), p);
     195             : }
     196             : static GEN
     197        3179 : FpX_otherroot(GEN x, GEN r, GEN p)
     198        3179 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     199             : 
     200             : /* disc(x^2+bx+c) = b^2 - 4c */
     201             : static ulong
     202    21650700 : Fl_disc_bc(ulong b, ulong c, ulong p)
     203    21650700 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     204             : /* p > 2 */
     205             : static ulong
     206    20512986 : Flx_quad_root(GEN x, ulong p, int unknown)
     207             : {
     208    20512986 :   ulong s, b = x[3], c = x[2];
     209    20512986 :   ulong D = Fl_disc_bc(b, c, p);
     210    20503001 :   if (unknown && krouu(D,p) == -1) return p;
     211    13696545 :   s = Fl_sqrt(D,p);
     212    13755596 :   if (s==~0UL) return p;
     213    13755583 :   return Fl_halve(Fl_sub(s,b, p), p);
     214             : }
     215             : static ulong
     216    12318192 : Flx_otherroot(GEN x, ulong r, ulong p)
     217    12318192 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     218             : 
     219             : /* 'todo' contains the list of factors to be split.
     220             :  * 'done' the list of finished factors, no longer touched */
     221             : struct split_t { GEN todo, done; };
     222             : static void
     223     4933570 : split_init(struct split_t *S, long max)
     224             : {
     225     4933570 :   S->todo = vectrunc_init(max);
     226     4933292 :   S->done = vectrunc_init(max);
     227     4932899 : }
     228             : #if 0
     229             : /* move todo[i] to done */
     230             : static void
     231             : split_convert(struct split_t *S, long i)
     232             : {
     233             :   long n = lg(S->todo)-1;
     234             :   vectrunc_append(S->done, gel(S->todo,i));
     235             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     236             :   setlg(S->todo, n);
     237             : }
     238             : #endif
     239             : /* append t to todo */
     240             : static void
     241     5180871 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     242             : /* delete todo[i], add t to done */
     243             : static void
     244     5180565 : split_moveto_done(struct split_t *S, long i, GEN t)
     245             : {
     246     5180565 :   long n = lg(S->todo)-1;
     247     5180565 :   vectrunc_append(S->done, t);
     248     5180816 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     249     5180816 :   setlg(S->todo, n);
     250             : 
     251     5180778 : }
     252             : /* append t to done */
     253             : static void
     254      398112 : split_add_done(struct split_t *S, GEN t)
     255      398112 : { vectrunc_append(S->done, t); }
     256             : /* split todo[i] into a and b */
     257             : static void
     258      340696 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     259             : {
     260      340696 :   gel(S->todo, i) = a;
     261      340696 :   split_add(S, b);
     262      340697 : }
     263             : /* split todo[i] into a and b, moved to done */
     264             : static void
     265      377196 : split_done(struct split_t *S, long i, GEN a, GEN b)
     266             : {
     267      377196 :   split_moveto_done(S, i, a);
     268      377197 :   split_add_done(S, b);
     269      377197 : }
     270             : 
     271             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     272             : static GEN
     273        3110 : FpX_roots_i(GEN f, GEN p)
     274             : {
     275             :   GEN pol, pol0, a, q;
     276             :   struct split_t S;
     277             : 
     278        3110 :   split_init(&S, lg(f)-1);
     279        3110 :   settyp(S.done, t_COL);
     280        3110 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     281        3110 :   switch(degpol(f))
     282             :   {
     283           7 :     case 0: return ZC_copy(S.done);
     284          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     285        1744 :     case 2: {
     286        1744 :       GEN s, r = FpX_quad_root(f, p, 1);
     287        1744 :       if (r) {
     288        1744 :         split_add_done(&S, r);
     289        1744 :         s = FpX_otherroot(f,r, p);
     290             :         /* f not known to be square free yet */
     291        1744 :         if (!equalii(r, s)) split_add_done(&S, s);
     292             :       }
     293        1744 :       return sort(S.done);
     294             :     }
     295             :   }
     296             : 
     297        1345 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     298        1345 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     299        1345 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     300        1345 :   a = FpX_gcd(f,a, p);
     301        1345 :   if (!degpol(a)) return ZC_copy(S.done);
     302        1345 :   split_add(&S, FpX_normalize(a,p));
     303             : 
     304        1345 :   q = shifti(p,-1);
     305        1345 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     306        1345 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     307        1345 :   for (pol0[2] = 1;; pol0[2]++)
     308        1475 :   {
     309        2820 :     long j, l = lg(S.todo);
     310        2820 :     if (l == 1) return sort(S.done);
     311        1482 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     312        3080 :     for (j = 1; j < l; j++)
     313             :     {
     314        1605 :       GEN b, r, s, c = gel(S.todo,j);
     315        1605 :       switch(degpol(c))
     316             :       { /* convert linear and quadratics to roots, try to split the rest */
     317         196 :         case 1:
     318         196 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     319         196 :           j--; l--; break;
     320        1272 :         case 2:
     321        1272 :           r = FpX_quad_root(c, p, 0);
     322        1272 :           if (!r) pari_err_PRIME("polrootsmod",p);
     323        1265 :           s = FpX_otherroot(c,r, p);
     324        1265 :           split_done(&S, j, r, s);
     325        1265 :           j--; l--; break;
     326         137 :         default:
     327         137 :           b = FpXQ_pow(pol,q, c,p);
     328         137 :           if (degpol(b) <= 0) continue;
     329         123 :           b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
     330         123 :           if (!degpol(b)) continue;
     331         123 :           b = FpX_normalize(b, p);
     332         123 :           c = FpX_div(c,b, p);
     333         123 :           split_todo(&S, j, b, c);
     334             :       }
     335             :     }
     336             :   }
     337             : }
     338             : 
     339             : /* Assume f is normalized */
     340             : static ulong
     341      174769 : Flx_cubic_root(GEN ff, ulong p)
     342             : {
     343      174769 :   GEN f = Flx_normalize(ff,p);
     344      174769 :   ulong pi = get_Fl_red(p);
     345      174769 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     346      174769 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     347      174770 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     348      174770 :   ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
     349      174770 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     350      174770 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     351      174770 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     352      174770 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     353      174770 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     354      174770 :   if (s!=~0UL)
     355             :   {
     356      110246 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     357      110246 :     if (p%3==2) /* 1 solutions */
     358       19533 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     359             :     else
     360             :     {
     361       90713 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     362       90713 :       if (vS1==~0UL) return p; /*0 solutions*/
     363             :       /*3 solutions*/
     364             :     }
     365       79594 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     366       79594 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     367             :   }
     368             :   else
     369             :   {
     370       64524 :     pari_sp av = avma;
     371       64524 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     372       64524 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     373             :     ulong Sa;
     374       64524 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     375       64524 :     Sa = vS1[1];
     376       64524 :     if (p%3==1) /*1 solutions*/
     377             :     {
     378       24280 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     379       24280 :       if (Fa!=P)
     380       16174 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     381             :     }
     382       64524 :     set_avma(av);
     383       64524 :     return Fl_sub(Fl_double(Sa,p),t,p);
     384             :   }
     385             : }
     386             : 
     387             : /* assume p > 2 prime */
     388             : static ulong
     389     3147361 : Flx_oneroot_i(GEN f, ulong p, long fl)
     390             : {
     391             :   GEN pol, a;
     392             :   ulong q;
     393             :   long da;
     394             : 
     395     3147361 :   if (Flx_val(f)) return 0;
     396     3146567 :   switch(degpol(f))
     397             :   {
     398       12371 :     case 1: return Fl_neg(f[2], p);
     399     2696574 :     case 2: return Flx_quad_root(f, p, 1);
     400      159102 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     401             :   }
     402             : 
     403      278526 :   if (!fl)
     404             :   {
     405      246721 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     406      246721 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     407      246721 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     408      246721 :     a = Flx_gcd(f,a, p);
     409       31805 :   } else a = f;
     410      278526 :   da = degpol(a);
     411      278526 :   if (!da) return p;
     412      199036 :   a = Flx_normalize(a,p);
     413             : 
     414      199036 :   q = p >> 1;
     415      199036 :   pol = polx_Flx(f[1]);
     416      199036 :   for(pol[2] = 1;; pol[2]++)
     417             :   {
     418      300065 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     419      300064 :     switch(da)
     420             :     {
     421      124222 :       case 1: return Fl_neg(a[2], p);
     422       59140 :       case 2: return Flx_quad_root(a, p, 0);
     423       15674 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     424             :       default: {
     425      101028 :         GEN b = Flxq_powu(pol,q, a,p);
     426             :         long db;
     427      101017 :         if (degpol(b) <= 0) continue;
     428       95964 :         b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
     429       95959 :         db = degpol(b); if (!db) continue;
     430       95958 :         b = Flx_normalize(b, p);
     431       95962 :         if (db <= (da >> 1)) {
     432       58549 :           a = b;
     433       58549 :           da = db;
     434             :         } else {
     435       37413 :           a = Flx_div(a,b, p);
     436       37422 :           da -= db;
     437             :         }
     438             :       }
     439             :     }
     440             :   }
     441             : }
     442             : 
     443             : /* assume p > 2 prime */
     444             : static GEN
     445        4847 : FpX_oneroot_i(GEN f, GEN p)
     446             : {
     447             :   GEN pol, pol0, a, q;
     448             :   long da;
     449             : 
     450        4847 :   if (ZX_val(f)) return gen_0;
     451        4581 :   switch(degpol(f))
     452             :   {
     453         674 :     case 1: return subii(p, gel(f,2));
     454        3837 :     case 2: return FpX_quad_root(f, p, 1);
     455             :   }
     456             : 
     457          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     458          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     459          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     460          70 :   a = FpX_gcd(f,a, p);
     461          70 :   da = degpol(a);
     462          70 :   if (!da) return NULL;
     463          70 :   a = FpX_normalize(a,p);
     464             : 
     465          70 :   q = shifti(p,-1);
     466          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     467          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     468          70 :   for (pol0[2]=1; ; pol0[2]++)
     469             :   {
     470         224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     471         224 :     switch(da)
     472             :     {
     473          42 :       case 1: return subii(p, gel(a,2));
     474          28 :       case 2: return FpX_quad_root(a, p, 0);
     475         154 :       default: {
     476         154 :         GEN b = FpXQ_pow(pol,q, a,p);
     477             :         long db;
     478         154 :         if (degpol(b) <= 0) continue;
     479         147 :         b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
     480         147 :         db = degpol(b); if (!db) continue;
     481         147 :         b = FpX_normalize(b, p);
     482         147 :         if (db <= (da >> 1)) {
     483         105 :           a = b;
     484         105 :           da = db;
     485             :         } else {
     486          42 :           a = FpX_div(a,b, p);
     487          42 :           da -= db;
     488             :         }
     489             :       }
     490             :     }
     491             :   }
     492             : }
     493             : 
     494             : ulong
     495     3106957 : Flx_oneroot(GEN f, ulong p)
     496             : {
     497     3106957 :   pari_sp av = avma;
     498             :   ulong r;
     499     3106957 :   switch(lg(f))
     500             :   {
     501           0 :     case 2: return 0;
     502           0 :     case 3: return p;
     503             :   }
     504     3106957 :   if (p == 2) return Flx_oneroot_mod_2(f);
     505     3106957 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     506     3106957 :   return gc_ulong(av,r);
     507             : }
     508             : 
     509             : ulong
     510       32667 : Flx_oneroot_split(GEN f, ulong p)
     511             : {
     512       32667 :   pari_sp av = avma;
     513             :   ulong r;
     514       32667 :   switch(lg(f))
     515             :   {
     516           0 :     case 2: return 0;
     517           0 :     case 3: return p;
     518             :   }
     519       32667 :   if (p == 2) return Flx_oneroot_mod_2(f);
     520       32667 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     521       32667 :   return gc_ulong(av,r);
     522             : }
     523             : 
     524             : /* assume that p is prime */
     525             : GEN
     526       12656 : FpX_oneroot(GEN f, GEN pp) {
     527       12656 :   pari_sp av = avma;
     528       12656 :   ZX_rootmod_init(&f, pp);
     529       12656 :   switch(lg(f))
     530             :   {
     531           0 :     case 2: set_avma(av); return gen_0;
     532           0 :     case 3: return gc_NULL(av);
     533             :   }
     534       12656 :   if (typ(f) == t_VECSMALL)
     535             :   {
     536        7809 :     ulong r, p = pp[2];
     537        7809 :     if (p == 2)
     538          70 :       r = Flx_oneroot_mod_2(f);
     539             :     else
     540        7739 :       r = Flx_oneroot_i(f, p, 0);
     541        7809 :     set_avma(av);
     542        7809 :     return (r == p)? NULL: utoi(r);
     543             :   }
     544        4847 :   f = FpX_oneroot_i(f, pp);
     545        4847 :   if (!f) return gc_NULL(av);
     546        4847 :   return gerepileuptoint(av, f);
     547             : }
     548             : 
     549             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     550             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     551             : /* returned in n                                                          */
     552             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     553             : /* twice as big and decrement until it divides p-1.                       */
     554             : static GEN
     555         651 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     556             : {
     557         651 :    pari_sp ltop = avma;
     558             :    GEN pm, factn, power, base, zeta;
     559             :    long n;
     560             : 
     561         651 :    pm = subis (p, 1ul);
     562         917 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     563         651 :    factn = Z_factor(stoi(n));
     564         651 :    power = diviuexact (pm, n);
     565         651 :    base = gen_1;
     566             :    do {
     567         938 :       base = addis (base, 1l);
     568         938 :       zeta = Fp_pow (base, power, p);
     569             :    }
     570         938 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     571         651 :    *pt_n = n;
     572         651 :    return gerepileuptoint (ltop, zeta);
     573             : }
     574             : 
     575             : GEN
     576         924 : FpX_oneroot_split(GEN fact, GEN p)
     577             : {
     578         924 :   pari_sp av = avma;
     579             :   long n, deg_f, i, dmin;
     580             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     581         924 :   fact = FpX_normalize(fact, p);
     582         924 :   deg_f = degpol(fact);
     583         924 :   if (deg_f<=2) return FpX_oneroot(fact, p);
     584         217 :   minfactor = fact; /* factor of minimal degree found so far */
     585         217 :   dmin = degpol(minfactor);
     586         217 :   xplusa = pol_x(varn(fact));
     587         868 :   while (dmin != 1)
     588             :   {
     589             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     590             :     /* zeta varies over the roots of unity in F_p                         */
     591         651 :     fact = minfactor; deg_f = dmin;
     592         651 :     zeta = gen_1;
     593         651 :     prim = good_root_of_unity(p, deg_f, 1, &n);
     594         651 :     expo = diviuexact(subiu(p, 1), n);
     595             :     /* update X+a, avoid a=0 */
     596         651 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     597         651 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     598        1491 :     for (i = 0; i < n; i++)
     599             :     {
     600        1099 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     601        1099 :       long dtmp = degpol(tmp);
     602        1099 :       if (dtmp > 0 && dtmp < deg_f)
     603             :       {
     604         455 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     605         455 :         if (dtmp < dmin)
     606             :         {
     607         448 :           minfactor = FpX_normalize (tmp, p);
     608         448 :           dmin = dtmp;
     609         448 :           if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
     610             :             /* stop early to avoid too many gcds */
     611             :             break;
     612             :         }
     613             :       }
     614         840 :       zeta = Fp_mul (zeta, prim, p);
     615             :     }
     616             :   }
     617         217 :   return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
     618             : }
     619             : 
     620             : /*******************************************************************/
     621             : /*                                                                 */
     622             : /*                     FACTORISATION MODULO p                      */
     623             : /*                                                                 */
     624             : /*******************************************************************/
     625             : 
     626             : /* F / E  a vector of vectors of factors / exponents of virtual length l
     627             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
     628             : static GEN
     629      626573 : FE_concat(GEN F, GEN E, long l)
     630             : {
     631      626573 :   setlg(E,l); E = shallowconcat1(E);
     632      626576 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
     633             : }
     634             : 
     635             : static GEN
     636          14 : ddf_to_ddf2_i(GEN V, long fl)
     637             : {
     638             :   GEN F, D;
     639          14 :   long i, j, l = lg(V);
     640          14 :   F = cgetg(l, t_VEC);
     641          14 :   D = cgetg(l, t_VECSMALL);
     642         112 :   for (i = j = 1; i < l; i++)
     643             :   {
     644          98 :     GEN Vi = gel(V,i);
     645          98 :     if ((fl==2 && F2x_degree(Vi) == 0)
     646          98 :       ||(fl==0 && degpol(Vi) == 0)) continue;
     647          35 :     gel(F,j) = Vi;
     648          35 :     uel(D,j) = i; j++;
     649             :   }
     650          14 :   setlg(F,j);
     651          14 :   setlg(D,j); return mkvec2(F,D);
     652             : }
     653             : 
     654             : GEN
     655           7 : ddf_to_ddf2(GEN V)
     656           7 : { return ddf_to_ddf2_i(V, 0); }
     657             : 
     658             : static GEN
     659           7 : F2x_ddf_to_ddf2(GEN V)
     660           7 : { return ddf_to_ddf2_i(V, 2); }
     661             : 
     662             : GEN
     663      752703 : vddf_to_simplefact(GEN V, long d)
     664             : {
     665             :   GEN E, F;
     666      752703 :   long i, j, c, l = lg(V);
     667      752703 :   F = cgetg(d+1, t_VECSMALL);
     668      752703 :   E = cgetg(d+1, t_VECSMALL);
     669     1510390 :   for (i = c = 1; i < l; i++)
     670             :   {
     671      757687 :     GEN Vi = gel(V,i);
     672      757687 :     long l = lg(Vi);
     673     5357387 :     for (j = 1; j < l; j++)
     674             :     {
     675     4599700 :       long k, n = degpol(gel(Vi,j)) / j;
     676     6623904 :       for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
     677             :     }
     678             :   }
     679      752703 :   setlg(F,c);
     680      752703 :   setlg(E,c);
     681      752703 :   return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
     682             : }
     683             : 
     684             : /* product of terms of degree 1 in factorization of f */
     685             : GEN
     686      147403 : FpX_split_part(GEN f, GEN p)
     687             : {
     688      147403 :   long n = degpol(f);
     689      147403 :   GEN z, X = pol_x(varn(f));
     690      147403 :   if (n <= 1) return f;
     691      145906 :   f = FpX_red(f, p);
     692      145906 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     693      145906 :   return FpX_gcd(z,f,p);
     694             : }
     695             : 
     696             : /* Compute the number of roots in Fp without counting multiplicity
     697             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     698             : long
     699       99261 : FpX_nbroots(GEN f, GEN p)
     700             : {
     701       99261 :   pari_sp av = avma;
     702       99261 :   GEN z = FpX_split_part(f, p);
     703       99261 :   return gc_long(av, degpol(z));
     704             : }
     705             : 
     706             : /* 1 < deg(f) <= p */
     707             : static int
     708           0 : Flx_is_totally_split_i(GEN f, ulong p)
     709             : {
     710           0 :   GEN F = Flx_Frobenius(f, p);
     711           0 :   return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
     712             : }
     713             : int
     714           0 : Flx_is_totally_split(GEN f, ulong p)
     715             : {
     716           0 :   pari_sp av = avma;
     717           0 :   ulong n = degpol(f);
     718           0 :   if (n <= 1) return 1;
     719           0 :   if (n > p) return 0; /* includes n < 0 */
     720           0 :   return gc_bool(av, Flx_is_totally_split_i(f,p));
     721             : }
     722             : int
     723           0 : FpX_is_totally_split(GEN f, GEN p)
     724             : {
     725           0 :   pari_sp av = avma;
     726           0 :   ulong n = degpol(f);
     727             :   int u;
     728           0 :   if (n <= 1) return 1;
     729           0 :   if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
     730           0 :   if (lgefint(p) != 3)
     731           0 :     u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
     732             :   else
     733             :   {
     734           0 :     ulong pp = (ulong)p[2];
     735           0 :     u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
     736             :   }
     737           0 :   return gc_bool(av, u);
     738             : }
     739             : 
     740             : long
     741     2539456 : Flx_nbroots(GEN f, ulong p)
     742             : {
     743     2539456 :   long n = degpol(f);
     744     2539456 :   pari_sp av = avma;
     745             :   GEN z;
     746     2539456 :   if (n <= 1) return n;
     747     2528977 :   if (n == 2)
     748             :   {
     749             :     ulong D;
     750       20226 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     751       18950 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     752       18950 :     return 1 + krouu(D,p);
     753             :   }
     754     2508751 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     755     2508757 :   z = Flx_gcd(z, f, p);
     756     2508762 :   return gc_long(av, degpol(z));
     757             : }
     758             : 
     759             : long
     760        4403 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
     761             : {
     762        4403 :   pari_sp av = avma;
     763             :   GEN X, b, g, xq;
     764             :   long i, j, n, v, B, l, m;
     765             :   pari_timer ti;
     766             :   hashtable h;
     767             : 
     768        4403 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     769        4403 :   X = pol_x(v);
     770        4403 :   if (ZX_equal(X,XP)) return 1;
     771        4403 :   B = n/2;
     772        4403 :   l = usqrt(B);
     773        4403 :   m = (B+l-1)/l;
     774        4403 :   T = FpX_get_red(T, p);
     775        4403 :   hash_init_GEN(&h, l+2, ZX_equal, 1);
     776        4403 :   hash_insert_long(&h, X,  0);
     777        4403 :   hash_insert_long(&h, XP, 1);
     778        4403 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     779        4403 :   b = XP;
     780        4403 :   xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
     781        4403 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
     782       10402 :   for (i = 3; i <= l+1; i++)
     783             :   {
     784        6790 :     b = FpX_FpXQV_eval(b, xq, T, p);
     785        6790 :     if (gequalX(b)) return gc_long(av,i-1);
     786        5999 :     hash_insert_long(&h, b, i-1);
     787             :   }
     788        3612 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
     789        3612 :   g = b;
     790        3612 :   xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
     791        3612 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
     792       12341 :   for(i = 2; i <= m+1; i++)
     793             :   {
     794       10745 :     g = FpX_FpXQV_eval(g, xq, T, p);
     795       10745 :     if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
     796             :   }
     797        1596 :   return gc_long(av,n);
     798             : }
     799             : 
     800             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     801             : static GEN
     802         753 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
     803             : {
     804             :   GEN b, g, h, F, f, Tr, xq;
     805             :   long i, j, n, v, B, l, m;
     806             :   pari_timer ti;
     807             : 
     808         753 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     809         753 :   if (n == 0) return cgetg(1, t_VEC);
     810         753 :   if (n == 1) return mkvec(get_FpX_mod(T));
     811         707 :   B = n/2;
     812         707 :   l = usqrt(B);
     813         707 :   m = (B+l-1)/l;
     814         707 :   T = FpX_get_red(T, p);
     815         707 :   b = cgetg(l+2, t_VEC);
     816         707 :   gel(b, 1) = pol_x(v);
     817         707 :   gel(b, 2) = XP;
     818         707 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     819         707 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     820         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
     821         893 :   for (i = 3; i <= l+1; i++)
     822         186 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     823         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
     824         707 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     825         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
     826         707 :   g = cgetg(m+1, t_VEC);
     827         707 :   gel(g, 1) = gel(xq, 2);
     828        1548 :   for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     829         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
     830         707 :   h = cgetg(m+1, t_VEC);
     831        2255 :   for (j = 1; j <= m; j++)
     832             :   {
     833        1548 :     pari_sp av = avma;
     834        1548 :     GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
     835        2643 :     for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
     836        1548 :     gel(h,j) = gerepileupto(av, e);
     837             :   }
     838         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
     839         707 :   Tr = get_FpX_mod(T);
     840         707 :   F = cgetg(m+1, t_VEC);
     841        2255 :   for (j = 1; j <= m; j++)
     842             :   {
     843        1548 :     GEN u = FpX_gcd(Tr, gel(h,j), p);
     844        1548 :     if (degpol(u))
     845             :     {
     846         303 :       u = FpX_normalize(u, p);
     847         303 :       Tr = FpX_div(Tr, u, p);
     848             :     }
     849        1548 :     gel(F,j) = u;
     850             :   }
     851         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
     852         707 :   f = const_vec(n, pol_1(v));
     853        2255 :   for (j = 1; j <= m; j++)
     854             :   {
     855        1548 :     GEN e = gel(F, j);
     856        1583 :     for (i=l-1; i >= 0; i--)
     857             :     {
     858        1583 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     859        1583 :       if (degpol(u))
     860             :       {
     861         310 :         u = FpX_normalize(u, p);
     862         310 :         gel(f, l*j-i) = u;
     863         310 :         e = FpX_div(e, u, p);
     864             :       }
     865        1583 :       if (!degpol(e)) break;
     866             :     }
     867             :   }
     868         707 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
     869         707 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     870         707 :   return f;
     871             : }
     872             : 
     873             : static void
     874           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     875             : {
     876           0 :   long n = degpol(Tp), r = n/d, ct = 0;
     877             :   GEN T, f, ff, p2;
     878           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     879           0 :   p2 = shifti(p,-1);
     880           0 :   T = FpX_get_red(Tp, p);
     881           0 :   XP = FpX_rem(XP, T, p);
     882             :   while (1)
     883           0 :   {
     884           0 :     pari_sp btop = avma;
     885             :     long i;
     886           0 :     GEN g = random_FpX(n, varn(Tp), p);
     887           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     888           0 :     if (signe(t) == 0) continue;
     889           0 :     for(i=1; i<=10; i++)
     890             :     {
     891           0 :       pari_sp btop2 = avma;
     892           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     893           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     894           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     895           0 :       set_avma(btop2);
     896             :     }
     897           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     898           0 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
     899           0 :     set_avma(btop);
     900             :   }
     901           0 :   f = FpX_normalize(f, p);
     902           0 :   ff = FpX_div(Tp, f ,p);
     903           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     904           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     905             : }
     906             : 
     907             : static void
     908         787 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     909             : {
     910             :   pari_sp av;
     911         787 :   GEN Tp = get_FpX_mod(T);
     912         787 :   long n = degpol(hp), vT = varn(Tp), ct = 0;
     913             :   GEN u1, u2, f1, f2, R, h;
     914         787 :   h = FpX_get_red(hp, p);
     915         787 :   t = FpX_rem(t, T, p);
     916         787 :   av = avma;
     917             :   do
     918             :   {
     919        1276 :     set_avma(av);
     920        1276 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     921        1276 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     922        1276 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
     923        1276 :   } while (degpol(u1)==0 || degpol(u1)==n);
     924         787 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     925         787 :   f1 = FpX_normalize(f1, p);
     926         787 :   u2 = FpX_div(hp, u1, p);
     927         787 :   f2 = FpX_div(Tp, f1, p);
     928         787 :   if (degpol(u1)==1)
     929         529 :     gel(V, idx) = f1;
     930             :   else
     931         258 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     932         787 :   idx += degpol(f1)/d;
     933         787 :   if (degpol(u2)==1)
     934         559 :     gel(V, idx) = f2;
     935             :   else
     936         228 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     937         787 : }
     938             : 
     939             : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
     940             : static void
     941         301 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     942             : {
     943         301 :   long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
     944             :   GEN T, h, t;
     945             :   pari_timer ti;
     946             : 
     947         301 :   T = FpX_get_red(Tp, p);
     948         301 :   XP = FpX_rem(XP, T, p);
     949         301 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     950             :   do
     951             :   {
     952         301 :     GEN g = random_FpX(n, vT, p);
     953         301 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     954         301 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     955         301 :     h = FpXQ_minpoly(t, T, p);
     956         301 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     957         301 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
     958         301 :   } while (degpol(h) != r);
     959         301 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     960         301 : }
     961             : 
     962             : static GEN
     963         739 : FpX_factor_Shoup(GEN T, GEN p)
     964             : {
     965         739 :   long i, n, s = 0;
     966             :   GEN XP, D, V;
     967         739 :   long e = expi(p);
     968             :   pari_timer ti;
     969         739 :   n = get_FpX_degree(T);
     970         739 :   T = FpX_get_red(T, p);
     971         739 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     972         739 :   XP = FpX_Frobenius(T, p);
     973         739 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     974         739 :   D = FpX_ddf_Shoup(T, XP, p);
     975         739 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
     976         739 :   s = ddf_to_nbfact(D);
     977         739 :   V = cgetg(s+1, t_COL);
     978        6042 :   for (i = 1, s = 1; i <= n; i++)
     979             :   {
     980        5303 :     GEN Di = gel(D,i);
     981        5303 :     long ni = degpol(Di), ri = ni/i;
     982        5303 :     if (ni == 0) continue;
     983         748 :     Di = FpX_normalize(Di, p);
     984         748 :     if (ni == i) { gel(V, s++) = Di; continue; }
     985         301 :     if (ri <= e*expu(e))
     986         301 :       FpX_edf(Di, XP, i, p, V, s);
     987             :     else
     988           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     989         301 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     990         301 :     s += ri;
     991             :   }
     992         739 :   return V;
     993             : }
     994             : 
     995             : long
     996      716094 : ddf_to_nbfact(GEN D)
     997             : {
     998      716094 :   long l = lg(D), i, s = 0;
     999     5129854 :   for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
    1000      716101 :   return s;
    1001             : }
    1002             : 
    1003             : /* Yun algorithm: Assume p > degpol(T) */
    1004             : static GEN
    1005         756 : FpX_factor_Yun(GEN T, GEN p)
    1006             : {
    1007         756 :   long n = degpol(T), i = 1;
    1008         756 :   GEN a, b, c, d = FpX_deriv(T, p);
    1009         756 :   GEN V = cgetg(n+1,t_VEC);
    1010         756 :   a = FpX_gcd(T, d, p);
    1011         756 :   if (degpol(a) == 0) return mkvec(T);
    1012         488 :   b = FpX_div(T, a, p);
    1013             :   do
    1014             :   {
    1015        2432 :     c = FpX_div(d, a, p);
    1016        2432 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1017        2432 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1018        2432 :     gel(V, i++) = a;
    1019        2432 :     b = FpX_div(b, a, p);
    1020        2432 :   } while (degpol(b));
    1021         488 :   setlg(V, i); return V;
    1022             : }
    1023             : GEN
    1024         777 : FpX_factor_squarefree(GEN T, GEN p)
    1025             : {
    1026         777 :   if (lgefint(p)==3)
    1027             :   {
    1028         777 :     ulong pp = (ulong)p[2];
    1029         777 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1030         777 :     return FlxV_to_ZXV(u);
    1031             :   }
    1032           0 :   return FpX_factor_Yun(T, p);
    1033             : }
    1034             : 
    1035             : long
    1036         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1037             : {
    1038         168 :   pari_sp av = avma;
    1039             :   GEN lc, F;
    1040         168 :   long i, l, n = degpol(f), v = varn(f);
    1041         168 :   if (n % k) return 0;
    1042         168 :   if (lgefint(p)==3)
    1043             :   {
    1044         126 :     ulong pp = p[2];
    1045         126 :     GEN fp = ZX_to_Flx(f, pp);
    1046         126 :     if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
    1047         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
    1048         105 :     return 1;
    1049             :   }
    1050          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1051          42 :   if (!lc) { av = avma; return 0; }
    1052          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1053        1491 :   for(i=1; i <= l; i++)
    1054        1456 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1055          35 :   if (pt_r)
    1056             :   {
    1057          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1058        1484 :     for (i=l; i>=1; i--)
    1059             :     {
    1060        1449 :       if (i%k) continue;
    1061         294 :       s = FpX_mul(s, gel(F,i), p);
    1062         294 :       r = FpX_mul(r, s, p);
    1063             :     }
    1064          35 :     *pt_r = gerepileupto(av, r);
    1065           0 :   } else av = avma;
    1066          35 :   return 1;
    1067             : }
    1068             : 
    1069             : static GEN
    1070         700 : FpX_factor_Cantor(GEN T, GEN p)
    1071             : {
    1072         700 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1073         700 :   long i, j, l = lg(V);
    1074         700 :   F = cgetg(l, t_VEC);
    1075         700 :   E = cgetg(l, t_VEC);
    1076        1916 :   for (i=1, j=1; i < l; i++)
    1077        1216 :     if (degpol(gel(V,i)))
    1078             :     {
    1079         739 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1080         739 :       gel(F, j) = Fj;
    1081         739 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1082         739 :       j++;
    1083             :     }
    1084         700 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1085             : }
    1086             : 
    1087             : static GEN
    1088           0 : FpX_ddf_i(GEN T, GEN p)
    1089             : {
    1090             :   GEN XP;
    1091           0 :   T = FpX_get_red(T, p);
    1092           0 :   XP = FpX_Frobenius(T, p);
    1093           0 :   return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
    1094             : }
    1095             : 
    1096             : GEN
    1097           7 : FpX_ddf(GEN f, GEN p)
    1098             : {
    1099           7 :   pari_sp av = avma;
    1100             :   GEN F;
    1101           7 :   switch(ZX_factmod_init(&f, p))
    1102             :   {
    1103           7 :     case 0:  F = F2x_ddf(f);
    1104           7 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    1105           0 :     case 1:  F = Flx_ddf(f,p[2]);
    1106           0 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    1107           0 :     default: F = FpX_ddf_i(f,p); break;
    1108             :   }
    1109           7 :   return gerepilecopy(av, F);
    1110             : }
    1111             : 
    1112             : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
    1113             : static GEN
    1114          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1115             : {
    1116             :   GEN V;
    1117             :   long i, l;
    1118          14 :   if (lgefint(p) == 3)
    1119             :   {
    1120           0 :     ulong pp = p[2];
    1121           0 :     return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
    1122             :   }
    1123          14 :   T = FpX_get_red(T, p);
    1124          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
    1125          28 :   for (i=1; i < l; i++)
    1126          14 :     gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
    1127          14 :   return vddf_to_simplefact(V, get_FpX_degree(T));
    1128             : }
    1129             : 
    1130             : static int
    1131           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1132             : {
    1133           0 :   pari_sp av = avma;
    1134             :   pari_timer ti;
    1135             :   long n;
    1136           0 :   GEN T = get_FpX_mod(Tp);
    1137           0 :   GEN dT = FpX_deriv(T, p);
    1138             :   GEN XP, D;
    1139           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    1140           0 :   n = get_FpX_degree(T);
    1141           0 :   T = FpX_get_red(Tp, p);
    1142           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1143           0 :   XP = FpX_Frobenius(T, p);
    1144           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1145           0 :   D = FpX_ddf_Shoup(T, XP, p);
    1146           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1147           0 :   return gc_bool(av, degpol(gel(D,n)) == n);
    1148             : }
    1149             : 
    1150             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1151             : 
    1152             : /*Assume that p is large and odd*/
    1153             : static GEN
    1154        1551 : FpX_factor_i(GEN f, GEN pp, long flag)
    1155             : {
    1156        1551 :   long d = degpol(f);
    1157        1551 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1158         714 :   switch(flag)
    1159             :   {
    1160         700 :     default: return FpX_factor_Cantor(f, pp);
    1161          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1162           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1163             :   }
    1164             : }
    1165             : 
    1166             : long
    1167           0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1168             : {
    1169           0 :   pari_sp av = avma;
    1170           0 :   long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
    1171           0 :   return gc_long(av,s);
    1172             : }
    1173             : 
    1174             : long
    1175           0 : FpX_nbfact(GEN T, GEN p)
    1176             : {
    1177           0 :   pari_sp av = avma;
    1178           0 :   GEN XP = FpX_Frobenius(T, p);
    1179           0 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1180           0 :   return gc_long(av,n);
    1181             : }
    1182             : 
    1183             : /* p > 2 */
    1184             : static GEN
    1185           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1186             : {
    1187           7 :   switch(d)
    1188             :   {
    1189           0 :     case -1:
    1190           0 :     case 0: return NULL;
    1191           0 :     case 1: return gen_1;
    1192             :   }
    1193           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1194             : }
    1195             : /* p > 2 */
    1196             : static GEN
    1197          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1198             : {
    1199          14 :   switch(d)
    1200             :   {
    1201           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1202           0 :     case 0: return trivial_fact();
    1203           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1204             :   }
    1205          14 :   switch(FpX_quad_factortype(f, p)) {
    1206           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1207           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1208           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1209             :   }
    1210             : }
    1211             : 
    1212             : GEN
    1213          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1214             : GEN
    1215       55471 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1216             : 
    1217             : /* not gerepile safe */
    1218             : static GEN
    1219         816 : FpX_factor_2(GEN f, GEN p, long d)
    1220             : {
    1221             :   GEN r, s, R, S;
    1222             :   long v;
    1223             :   int sgn;
    1224         816 :   switch(d)
    1225             :   {
    1226           7 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1227          30 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1228          74 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1229             :   }
    1230         705 :   r = FpX_quad_root(f, p, 1);
    1231         705 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1232         170 :   v = varn(f);
    1233         170 :   s = FpX_otherroot(f, r, p);
    1234         170 :   if (signe(r)) r = subii(p, r);
    1235         170 :   if (signe(s)) s = subii(p, s);
    1236         170 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1237         170 :   R = deg1pol_shallow(gen_1, r, v);
    1238         170 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1239          77 :   S = deg1pol_shallow(gen_1, s, v);
    1240          77 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1241             : }
    1242             : static GEN
    1243         837 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1244             : {
    1245         837 :   switch(flag) {
    1246           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1247          14 :     case 1: return FpX_degfact_2(f, p, d);
    1248         816 :     default: return FpX_factor_2(f, p, d);
    1249             :   }
    1250             : }
    1251             : 
    1252             : static int
    1253       77703 : F2x_quad_factortype(GEN x)
    1254       77703 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1255             : 
    1256             : static GEN
    1257           7 : F2x_is_irred_2(GEN f, long d)
    1258           7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1259             : 
    1260             : static GEN
    1261       10969 : F2x_degfact_2(GEN f, long d)
    1262             : {
    1263       10969 :   if (!d) return trivial_fact();
    1264       10969 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1265       10759 :   switch(F2x_quad_factortype(f)) {
    1266        3339 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1267        3276 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1268        4144 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1269             :   }
    1270             : }
    1271             : 
    1272             : static GEN
    1273      124459 : F2x_factor_2(GEN f, long d)
    1274             : {
    1275      124459 :   long v = f[1];
    1276      124459 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1277      114303 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1278       65470 :   switch(F2x_quad_factortype(f))
    1279             :   {
    1280       15551 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1281       29338 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1282       20581 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1283             :   }
    1284             : }
    1285             : static GEN
    1286      135435 : F2x_factor_deg2(GEN f, long d, long flag)
    1287             : {
    1288      135435 :   switch(flag) {
    1289           7 :     case 2: return F2x_is_irred_2(f, d);
    1290       10969 :     case 1: return F2x_degfact_2(f, d);
    1291      124459 :     default: return F2x_factor_2(f, d);
    1292             :   }
    1293             : }
    1294             : 
    1295             : /* xt = NULL or x^(p-1)/2 mod g */
    1296             : static void
    1297        6199 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
    1298             : {
    1299        6199 :   ulong q = p >> 1;
    1300        6199 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1301        6199 :   long d = degpol(a);
    1302        6199 :   if (d < 0)
    1303             :   {
    1304             :     ulong i;
    1305         322 :     split_add_done(S, (GEN)1);
    1306         917 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1307             :   } else {
    1308        5877 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1309        5877 :     if (d)
    1310             :     {
    1311        5877 :       if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
    1312        5877 :       a = Flx_gcd(a, xt, p);
    1313        5877 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1314             :     }
    1315             :   }
    1316        6199 : }
    1317             : static void
    1318        6199 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
    1319             : {
    1320        6199 :   ulong q = p >> 1;
    1321        6199 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1322        6199 :   long d = degpol(a);
    1323        6199 :   if (d < 0)
    1324             :   {
    1325          28 :     ulong i, z = nonsquare_Fl(p);
    1326          28 :     split_add_done(S, (GEN)z);
    1327          63 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1328             :   } else {
    1329        6171 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1330        6171 :     if (d)
    1331             :     {
    1332        6150 :       if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
    1333        6150 :       a = Flx_gcd(a, xt, p);
    1334        6150 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1335             :     }
    1336             :   }
    1337        6199 : }
    1338             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1339             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1340             : static int
    1341     4929751 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1342             : {
    1343     4929751 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1344     4929665 :   long d = degpol(g);
    1345     4929582 :   if (d < 0) return 0;
    1346     4929540 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1347     4929544 :   if (!d) return 1;
    1348     4913258 :   if ((p >> 4) <= (ulong)d)
    1349             :   { /* small p; split directly using x^((p-1)/2) +/- 1 */
    1350        5646 :     GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
    1351        6199 :                                 : NULL;
    1352        6199 :     split_squares(S, g, p, xt);
    1353        6199 :     split_nonsquares(S, g, p, xt);
    1354             :   } else { /* large p; use x^(p-1) - 1 directly */
    1355     4907059 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1356     4894746 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1357     4894746 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1358     4897067 :     g = Flx_gcd(g,a, p);
    1359     4899521 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1360             :   }
    1361     4911962 :   return 1;
    1362             : }
    1363             : 
    1364             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1365             : static GEN
    1366    21982563 : Flx_roots_i(GEN f, ulong p)
    1367             : {
    1368             :   GEN pol, g;
    1369    21982563 :   long v = Flx_valrem(f, &g);
    1370             :   ulong q;
    1371             :   struct split_t S;
    1372             : 
    1373             :   /* optimization: test for small degree first */
    1374    21980846 :   switch(degpol(g))
    1375             :   {
    1376       28887 :     case 1: {
    1377       28887 :       ulong r = p - g[2];
    1378       28887 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1379             :     }
    1380    17011834 :     case 2: {
    1381    17011834 :       ulong r = Flx_quad_root(g, p, 1), s;
    1382    17029298 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1383    11620915 :       s = Flx_otherroot(g,r, p);
    1384    11692817 :       if (r < s)
    1385     2927781 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1386     8765036 :       else if (r > s)
    1387     8791389 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1388             :       else
    1389         138 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1390             :     }
    1391             :   }
    1392     4930604 :   q = p >> 1;
    1393     4930604 :   split_init(&S, lg(f)-1);
    1394     4929721 :   settyp(S.done, t_VECSMALL);
    1395     4929721 :   if (v) split_add_done(&S, (GEN)0);
    1396     4929721 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1397          42 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1398     4928130 :   pol = polx_Flx(f[1]);
    1399     4928374 :   for (pol[2]=1; ; pol[2]++)
    1400     5258597 :   {
    1401    10186971 :     long j, l = lg(S.todo);
    1402    10186971 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1403     5258385 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1404    10879352 :     for (j = 1; j < l; j++)
    1405             :     {
    1406     5620755 :       GEN b, c = gel(S.todo,j);
    1407             :       ulong r, s;
    1408     5620755 :       switch(degpol(c))
    1409             :       {
    1410     4803088 :         case 1:
    1411     4803088 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1412     4803382 :           j--; l--; break;
    1413      375918 :         case 2:
    1414      375918 :           r = Flx_quad_root(c, p, 0);
    1415      375938 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1416      375931 :           s = Flx_otherroot(c,r, p);
    1417      375931 :           split_done(&S, j, (GEN)r, (GEN)s);
    1418      375931 :           j--; l--; break;
    1419      441626 :         default:
    1420      441626 :           b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
    1421      441719 :           if (degpol(b) <= 0) continue;
    1422      340688 :           b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
    1423      340704 :           if (!degpol(b)) continue;
    1424      340566 :           b = Flx_normalize(b, p);
    1425      340574 :           c = Flx_div(c,b, p);
    1426      340569 :           split_todo(&S, j, b, c);
    1427             :       }
    1428             :     }
    1429             :   }
    1430             : }
    1431             : 
    1432             : GEN
    1433    21926558 : Flx_roots(GEN f, ulong p)
    1434             : {
    1435    21926558 :   pari_sp av = avma;
    1436    21926558 :   switch(lg(f))
    1437             :   {
    1438           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1439           0 :     case 3: set_avma(av); return cgetg(1, t_VECSMALL);
    1440             :   }
    1441    21927691 :   if (p == 2) return Flx_root_mod_2(f);
    1442    21927691 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1443             : }
    1444             : 
    1445             : /* assume x reduced mod p, monic. */
    1446             : static int
    1447     1143513 : Flx_quad_factortype(GEN x, ulong p)
    1448             : {
    1449     1143513 :   ulong b = x[3], c = x[2];
    1450     1143513 :   return krouu(Fl_disc_bc(b, c, p), p);
    1451             : }
    1452             : static GEN
    1453          56 : Flx_is_irred_2(GEN f, ulong p, long d)
    1454             : {
    1455          56 :   if (!d) return NULL;
    1456          56 :   if (d == 1) return gen_1;
    1457          56 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1458             : }
    1459             : static GEN
    1460     1170211 : Flx_degfact_2(GEN f, ulong p, long d)
    1461             : {
    1462     1170211 :   if (!d) return trivial_fact();
    1463     1170211 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1464     1143457 :   switch(Flx_quad_factortype(f, p)) {
    1465      545426 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1466      584612 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1467       13419 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1468             :   }
    1469             : }
    1470             : /* p > 2 */
    1471             : static GEN
    1472      502576 : Flx_factor_2(GEN f, ulong p, long d)
    1473             : {
    1474             :   ulong r, s;
    1475             :   GEN R,S;
    1476      502576 :   long v = f[1];
    1477      502576 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1478      484606 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1479      376969 :   r = Flx_quad_root(f, p, 1);
    1480      376990 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1481      251320 :   s = Flx_otherroot(f, r, p);
    1482      251323 :   r = Fl_neg(r, p);
    1483      251323 :   s = Fl_neg(s, p);
    1484      251323 :   if (s < r) lswap(s,r);
    1485      251323 :   R = mkvecsmall3(v,r,1);
    1486      251323 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1487      216353 :   S = mkvecsmall3(v,s,1);
    1488      216351 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1489             : }
    1490             : static GEN
    1491     1672844 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1492             : {
    1493     1672844 :   switch(flag) {
    1494          56 :     case 2: return Flx_is_irred_2(f, p, d);
    1495     1170211 :     case 1: return Flx_degfact_2(f, p, d);
    1496      502577 :     default: return Flx_factor_2(f, p, d);
    1497             :   }
    1498             : }
    1499             : 
    1500             : static GEN
    1501       11049 : F2x_Berlekamp_ker(GEN u)
    1502             : {
    1503       11049 :   pari_sp ltop=avma;
    1504       11049 :   long j,N = F2x_degree(u);
    1505             :   GEN Q;
    1506             :   pari_timer T;
    1507       11049 :   timer_start(&T);
    1508       11049 :   Q = F2x_matFrobenius(u);
    1509      254364 :   for (j=1; j<=N; j++)
    1510      243315 :     F2m_flip(Q,j,j);
    1511       11049 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
    1512       11049 :   Q = F2m_ker_sp(Q,0);
    1513       11049 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
    1514       11049 :   return gerepileupto(ltop,Q);
    1515             : }
    1516             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    1517             : static long
    1518       15543 : F2x_split_Berlekamp(GEN *t)
    1519             : {
    1520       15543 :   GEN u = *t, a, b, vker;
    1521       15543 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    1522             : 
    1523       15543 :   if (du == 1) return 1;
    1524       11803 :   if (du == 2)
    1525             :   {
    1526         755 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    1527             :     {
    1528           0 :       t[0] = mkvecsmall2(sv, 2);
    1529           0 :       t[1] = mkvecsmall2(sv, 3);
    1530           0 :       return 2;
    1531             :     }
    1532         755 :     return 1;
    1533             :   }
    1534             : 
    1535       11048 :   vker = F2x_Berlekamp_ker(u);
    1536       11049 :   lb = lgcols(vker);
    1537       11062 :   d = lg(vker)-1;
    1538       11062 :   ir = 0;
    1539             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    1540       43079 :   for (L=1; L<d; )
    1541             :   {
    1542             :     GEN pol;
    1543       32030 :     if (d == 2)
    1544        1906 :       pol = F2v_to_F2x(gel(vker,2), sv);
    1545             :     else
    1546             :     {
    1547       30124 :       GEN v = zero_zv(lb);
    1548       30127 :       v[1] = du;
    1549       30127 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    1550      112096 :       for (i=2; i<=d; i++)
    1551       81963 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    1552       30133 :       pol = F2v_to_F2x(v, sv);
    1553             :     }
    1554       98873 :     for (i=ir; i<L && L<d; i++)
    1555             :     {
    1556       66856 :       a = t[i]; la = F2x_degree(a);
    1557       66856 :       if (la == 1) { set_irred(i); }
    1558       66673 :       else if (la == 2)
    1559             :       {
    1560         712 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    1561             :         {
    1562           0 :           t[i] = mkvecsmall2(sv, 2);
    1563           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    1564             :         }
    1565         712 :         set_irred(i);
    1566             :       }
    1567             :       else
    1568             :       {
    1569       65961 :         pari_sp av = avma;
    1570             :         long lb;
    1571       65961 :         b = F2x_rem(pol, a);
    1572       65955 :         if (F2x_degree(b) <= 0) { set_avma(av); continue; }
    1573       22010 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    1574       22012 :         if (lb && lb < la)
    1575             :         {
    1576       22012 :           t[L] = F2x_div(a,b);
    1577       22011 :           t[i]= b; L++;
    1578             :         }
    1579           0 :         else set_avma(av);
    1580             :       }
    1581             :     }
    1582             :   }
    1583       11049 :   return d;
    1584             : }
    1585             : /* assume deg f > 2 */
    1586             : static GEN
    1587       12781 : F2x_Berlekamp_i(GEN f, long flag)
    1588             : {
    1589       12781 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    1590             :   GEN y, E, t, V;
    1591             : 
    1592       12781 :   val = F2x_valrem(f, &f);
    1593       12781 :   if (flag == 2 && val) return NULL;
    1594       12767 :   V = F2x_factor_squarefree(f); lV = lg(V);
    1595       12766 :   if (flag == 2 && lV > 2) return NULL;
    1596             : 
    1597             :   /* to hold factors and exponents */
    1598       12696 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    1599       12697 :   E = cgetg(d+1,t_VECSMALL);
    1600       12697 :   lfact = 1;
    1601       12697 :   if (val) {
    1602        3643 :     if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
    1603        3643 :     E[1] = val; lfact++;
    1604             :   }
    1605             : 
    1606       64381 :   for (k=1; k<lV; k++)
    1607             :   {
    1608       51840 :     if (F2x_degree(gel(V, k))==0) continue;
    1609       15543 :     gel(t,lfact) = gel(V, k);
    1610       15543 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    1611       15542 :     if (flag == 2 && d != 1) return NULL;
    1612       15388 :     if (flag == 1)
    1613        4361 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    1614       52375 :     for (j=0; j<d; j++) E[lfact+j] = k;
    1615       15388 :     lfact += d;
    1616             :   }
    1617       12541 :   if (flag == 2) return gen_1; /* irreducible */
    1618       12527 :   setlg(t, lfact);
    1619       12527 :   setlg(E, lfact); y = mkvec2(t,E);
    1620        2247 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    1621       14775 :               : sort_factor_pol(y, cmpGuGu);
    1622             : }
    1623             : 
    1624             : /* Adapted from Shoup NTL */
    1625             : GEN
    1626       85575 : F2x_factor_squarefree(GEN f)
    1627             : {
    1628             :   GEN r, t, v, tv;
    1629       85575 :   long i, q, n = F2x_degree(f);
    1630       85573 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1631       85576 :   for(q = 1;;q *= 2)
    1632             :   {
    1633      154930 :     r = F2x_gcd(f, F2x_deriv(f));
    1634      154924 :     if (F2x_degree(r) == 0)
    1635             :     {
    1636       70220 :       gel(u, q) = f;
    1637       70220 :       break;
    1638             :     }
    1639       84708 :     t = F2x_div(f, r);
    1640       84708 :     if (F2x_degree(t) > 0)
    1641             :     {
    1642             :       long j;
    1643       34476 :       for(j = 1;;j++)
    1644             :       {
    1645       76632 :         v = F2x_gcd(r, t);
    1646       76631 :         tv = F2x_div(t, v);
    1647       76630 :         if (F2x_degree(tv) > 0)
    1648       36019 :           gel(u, j*q) = tv;
    1649       76631 :         if (F2x_degree(v) <= 0) break;
    1650       42156 :         r = F2x_div(r, v);
    1651       42156 :         t = v;
    1652             :       }
    1653       34476 :       if (F2x_degree(r) == 0) break;
    1654             :     }
    1655       69353 :     f = F2x_sqrt(r);
    1656             :   }
    1657      601987 :   for (i = n; i; i--)
    1658      601615 :     if (F2x_degree(gel(u,i))) break;
    1659       85572 :   setlg(u,i+1); return u;
    1660             : }
    1661             : 
    1662             : static GEN
    1663       90191 : F2x_ddf_simple(GEN T, GEN XP)
    1664             : {
    1665       90191 :   pari_sp av = avma, av2;
    1666             :   GEN f, z, Tr, X;
    1667       90191 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1668       90191 :   if (n == 0) return cgetg(1, t_VEC);
    1669       90191 :   if (n == 1) return mkvec(T);
    1670       43349 :   z = XP; Tr = T; X = polx_F2x(v);
    1671       43351 :   f = const_vec(n, pol1_F2x(v));
    1672       43351 :   av2 = avma;
    1673      148040 :   for (j = 1; j <= B; j++)
    1674             :   {
    1675      110236 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1676      110246 :     if (F2x_degree(u))
    1677             :     {
    1678       24164 :       gel(f, j) = u;
    1679       24164 :       Tr = F2x_div(Tr, u);
    1680       24166 :       av2 = avma;
    1681       86106 :     } else z = gerepileuptoleaf(av2, z);
    1682      110297 :     if (!F2x_degree(Tr)) break;
    1683      104715 :     z = F2xq_sqr(z, Tr);
    1684             :   }
    1685       43341 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1686       43350 :   return gerepilecopy(av, f);
    1687             : }
    1688             : 
    1689             : GEN
    1690           7 : F2x_ddf(GEN T)
    1691             : {
    1692             :   GEN XP;
    1693           7 :   T = F2x_get_red(T);
    1694           7 :   XP = F2x_Frobenius(T);
    1695           7 :   return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
    1696             : }
    1697             : 
    1698             : static GEN
    1699        8447 : F2xq_frobtrace(GEN a, long d, GEN T)
    1700             : {
    1701        8447 :   pari_sp av = avma;
    1702             :   long i;
    1703        8447 :   GEN x = a;
    1704       33901 :   for(i=1; i<d; i++)
    1705             :   {
    1706       25455 :     x = F2x_add(a, F2xq_sqr(x,T));
    1707       25454 :     if (gc_needed(av, 2))
    1708           0 :       x = gerepileuptoleaf(av, x);
    1709             :   }
    1710        8446 :   return x;
    1711             : }
    1712             : 
    1713             : static void
    1714       12846 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1715             : {
    1716       12846 :   long n = F2x_degree(Tp), r = n/d;
    1717             :   GEN T, f, ff;
    1718       12846 :   if (r==1) { gel(V, idx) = Tp; return; }
    1719        4364 :   T = Tp;
    1720        4364 :   XP = F2x_rem(XP, T);
    1721             :   while (1)
    1722        4083 :   {
    1723        8447 :     pari_sp btop = avma;
    1724             :     long df;
    1725        8447 :     GEN g = random_F2x(n, Tp[1]);
    1726        8447 :     GEN t = F2xq_frobtrace(g, d, T);
    1727        8447 :     if (lgpol(t) == 0) continue;
    1728        6349 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1729        6349 :     if (df > 0 && df < n) break;
    1730        1985 :     set_avma(btop);
    1731             :   }
    1732        4364 :   ff = F2x_div(Tp, f);
    1733        4364 :   F2x_edf_simple(f, XP, d, V, idx);
    1734        4364 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1735             : }
    1736             : 
    1737             : static GEN
    1738       90185 : F2x_factor_Shoup(GEN T)
    1739             : {
    1740       90185 :   long i, n, s = 0;
    1741             :   GEN XP, D, V;
    1742             :   pari_timer ti;
    1743       90185 :   n = F2x_degree(T);
    1744       90185 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1745       90185 :   XP = F2x_Frobenius(T);
    1746       90184 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1747       90184 :   D = F2x_ddf_simple(T, XP);
    1748       90184 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1749      384211 :   for (i = 1; i <= n; i++)
    1750      294033 :     s += F2x_degree(gel(D,i))/i;
    1751       90178 :   V = cgetg(s+1, t_COL);
    1752      384257 :   for (i = 1, s = 1; i <= n; i++)
    1753             :   {
    1754      294036 :     GEN Di = gel(D,i);
    1755      294036 :     long ni = F2x_degree(Di), ri = ni/i;
    1756      294021 :     if (ni == 0) continue;
    1757      108757 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1758        4068 :     F2x_edf_simple(Di, XP, i, V, s);
    1759        4118 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1760        4118 :     s += ri;
    1761             :   }
    1762       90221 :   return V;
    1763             : }
    1764             : 
    1765             : static GEN
    1766       72808 : F2x_factor_Cantor(GEN T)
    1767             : {
    1768       72808 :   GEN E, F, V = F2x_factor_squarefree(T);
    1769       72806 :   long i, j, l = lg(V);
    1770       72806 :   E = cgetg(l, t_VEC);
    1771       72809 :   F = cgetg(l, t_VEC);
    1772      292156 :   for (i=1, j=1; i < l; i++)
    1773      219347 :     if (F2x_degree(gel(V,i)))
    1774             :     {
    1775       90185 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1776       90185 :       gel(F, j) = Fj;
    1777       90185 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1778       90185 :       j++;
    1779             :     }
    1780       72809 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1781             : }
    1782             : 
    1783             : #if 0
    1784             : static GEN
    1785             : F2x_simplefact_Shoup(GEN T)
    1786             : {
    1787             :   long i, n, s = 0, j = 1, k;
    1788             :   GEN XP, D, V;
    1789             :   pari_timer ti;
    1790             :   n = F2x_degree(T);
    1791             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1792             :   XP = F2x_Frobenius(T);
    1793             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1794             :   D = F2x_ddf_simple(T, XP);
    1795             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1796             :   for (i = 1; i <= n; i++)
    1797             :     s += F2x_degree(gel(D,i))/i;
    1798             :   V = cgetg(s+1, t_VECSMALL);
    1799             :   for (i = 1; i <= n; i++)
    1800             :   {
    1801             :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1802             :     if (ni == 0) continue;
    1803             :     for (k = 1; k <= ri; k++)
    1804             :       V[j++] = i;
    1805             :   }
    1806             :   return V;
    1807             : }
    1808             : static GEN
    1809             : F2x_simplefact_Cantor(GEN T)
    1810             : {
    1811             :   GEN E, F, V = F2x_factor_squarefree(T);
    1812             :   long i, j, l = lg(V);
    1813             :   F = cgetg(l, t_VEC);
    1814             :   E = cgetg(l, t_VEC);
    1815             :   for (i=1, j=1; i < l; i++)
    1816             :     if (F2x_degree(gel(V,i)))
    1817             :     {
    1818             :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1819             :       gel(F, j) = Fj;
    1820             :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1821             :       j++;
    1822             :     }
    1823             :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1824             : }
    1825             : static int
    1826             : F2x_isirred_Cantor(GEN T)
    1827             : {
    1828             :   pari_sp av = avma;
    1829             :   pari_timer ti;
    1830             :   long n;
    1831             :   GEN dT = F2x_deriv(T);
    1832             :   GEN XP, D;
    1833             :   if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
    1834             :   n = F2x_degree(T);
    1835             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1836             :   XP = F2x_Frobenius(T);
    1837             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1838             :   D = F2x_ddf_simple(T, XP);
    1839             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1840             :   return gc_bool(av, F2x_degree(gel(D,n)) == n);
    1841             : }
    1842             : #endif
    1843             : 
    1844             : /* driver for Cantor factorization, assume deg f > 2; not competitive for
    1845             :  * flag != 0, or as deg f increases */
    1846             : static GEN
    1847       72808 : F2x_Cantor_i(GEN f, long flag)
    1848             : {
    1849             :   switch(flag)
    1850             :   {
    1851       72808 :     default: return F2x_factor_Cantor(f);
    1852             : #if 0
    1853             :     case 1: return F2x_simplefact_Cantor(f);
    1854             :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1855             : #endif
    1856             :   }
    1857             : }
    1858             : static GEN
    1859      221023 : F2x_factor_i(GEN f, long flag)
    1860             : {
    1861      221023 :   long d = F2x_degree(f);
    1862      221024 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1863       83090 :   return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
    1864      168680 :                                : F2x_Berlekamp_i(f, flag);
    1865             : }
    1866             : 
    1867             : GEN
    1868           0 : F2x_degfact(GEN f)
    1869             : {
    1870           0 :   pari_sp av = avma;
    1871           0 :   GEN z = F2x_factor_i(f, 1);
    1872           0 :   return gerepilecopy(av, z);
    1873             : }
    1874             : 
    1875             : int
    1876         238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
    1877             : 
    1878             : /* Adapted from Shoup NTL */
    1879             : GEN
    1880     1309674 : Flx_factor_squarefree(GEN f, ulong p)
    1881             : {
    1882     1309674 :   long i, q, n = degpol(f);
    1883     1309670 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1884     1309675 :   for(q = 1;;q *= p)
    1885       55336 :   {
    1886     1365011 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1887     1364991 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1888      113856 :     t = Flx_div(f, r, p);
    1889      113856 :     if (degpol(t) > 0)
    1890             :     {
    1891             :       long j;
    1892       59611 :       for(j = 1;;j++)
    1893             :       {
    1894      204575 :         v = Flx_gcd(r, t, p);
    1895      204575 :         tv = Flx_div(t, v, p);
    1896      204575 :         if (degpol(tv) > 0)
    1897       87773 :           gel(u, j*q) = Flx_normalize(tv, p);
    1898      204575 :         if (degpol(v) <= 0) break;
    1899      144964 :         r = Flx_div(r, v, p);
    1900      144964 :         t = v;
    1901             :       }
    1902       59611 :       if (degpol(r) == 0) break;
    1903             :     }
    1904       55336 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1905             :   }
    1906     7937224 :   for (i = n; i; i--)
    1907     7937240 :     if (degpol(gel(u,i))) break;
    1908     1309654 :   setlg(u,i+1); return u;
    1909             : }
    1910             : 
    1911             : long
    1912        3276 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1913             : {
    1914        3276 :   pari_sp av = avma;
    1915             :   ulong lc;
    1916             :   GEN F;
    1917        3276 :   long i, n = degpol(f), v = f[1], l;
    1918        3276 :   if (n % k) return 0;
    1919        3276 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    1920        3276 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    1921        3276 :   F = Flx_factor_squarefree(f, p); l = lg(F)-1;
    1922       37492 :   for (i = 1; i <= l; i++)
    1923       34237 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1924        3255 :   if (pt_r)
    1925             :   {
    1926        3255 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    1927       37471 :     for(i = l; i >= 1; i--)
    1928             :     {
    1929       34216 :       if (i%k) continue;
    1930       12110 :       s = Flx_mul(s, gel(F,i), p);
    1931       12110 :       r = Flx_mul(r, s, p);
    1932             :     }
    1933        3255 :     *pt_r = gerepileuptoleaf(av, r);
    1934           0 :   } else set_avma(av);
    1935        3255 :   return 1;
    1936             : }
    1937             : 
    1938             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1939             : static GEN
    1940     1764157 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
    1941             : {
    1942     1764157 :   pari_sp av = avma;
    1943             :   GEN b, g, h, F, f, Tr, xq;
    1944             :   long i, j, n, v, bo, ro;
    1945             :   long B, l, m;
    1946             :   pari_timer ti;
    1947     1764157 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1948     1764159 :   if (n == 0) return cgetg(1, t_VEC);
    1949     1760554 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1950     1655392 :   B = n/2;
    1951     1655392 :   l = usqrt(B);
    1952     1655390 :   m = (B+l-1)/l;
    1953     1655390 :   T = Flx_get_red(T, p);
    1954     1655388 :   b = cgetg(l+2, t_VEC);
    1955     1655389 :   gel(b, 1) = polx_Flx(v);
    1956     1655396 :   gel(b, 2) = XP;
    1957     1655396 :   bo = brent_kung_optpow(n, l-1, 1);
    1958     1655407 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    1959     1655407 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1960     1655407 :   if (expu(p) <= ro)
    1961      311879 :     for (i = 3; i <= l+1; i++)
    1962      176131 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1963             :   else
    1964             :   {
    1965     1519659 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1966     1519646 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
    1967     1792677 :     for (i = 3; i <= l+1; i++)
    1968      273027 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1969             :   }
    1970     1655398 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
    1971     1655398 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1972     1655383 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
    1973     1655383 :   g = cgetg(m+1, t_VEC);
    1974     1655388 :   gel(g, 1) = gel(xq, 2);
    1975     3389838 :   for(i = 2; i <= m; i++)
    1976     1734471 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1977     1655367 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
    1978     1655367 :   h = cgetg(m+1, t_VEC);
    1979     5045243 :   for (j = 1; j <= m; j++)
    1980             :   {
    1981     3389876 :     pari_sp av = avma;
    1982     3389876 :     GEN gj = gel(g, j);
    1983     3389876 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1984     5015073 :     for (i = 2; i <= l; i++)
    1985     1625308 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1986     3389765 :     gel(h, j) = gerepileupto(av, e);
    1987             :   }
    1988     1655367 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
    1989     1655367 :   Tr = get_Flx_mod(T);
    1990     1655390 :   F = cgetg(m+1, t_VEC);
    1991     5045226 :   for (j = 1; j <= m; j++)
    1992             :   {
    1993     3389836 :     GEN u = Flx_gcd(Tr, gel(h, j), p);
    1994     3389822 :     if (degpol(u))
    1995             :     {
    1996     1316502 :       u = Flx_normalize(u, p);
    1997     1316502 :       Tr = Flx_div(Tr, u, p);
    1998             :     }
    1999     3389798 :     gel(F, j) = u;
    2000             :   }
    2001     1655390 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
    2002     1655390 :   f = const_vec(n, pol1_Flx(v));
    2003     5045206 :   for (j = 1; j <= m; j++)
    2004             :   {
    2005     3389814 :     GEN e = gel(F, j);
    2006     3713564 :     for (i=l-1; i >= 0; i--)
    2007             :     {
    2008     3713574 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    2009     3713537 :       if (degpol(u))
    2010             :       {
    2011     1347790 :         gel(f, l*j-i) = u;
    2012     1347790 :         e = Flx_div(e, u, p);
    2013             :       }
    2014     3713519 :       if (!degpol(e)) break;
    2015             :     }
    2016             :   }
    2017     1655392 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
    2018     1655392 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    2019     1655393 :   return gerepilecopy(av, f);
    2020             : }
    2021             : 
    2022             : static void
    2023      186633 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2024             : {
    2025      186633 :   long n = degpol(Tp), r = n/d;
    2026             :   GEN T, f, ff;
    2027             :   ulong p2;
    2028      186633 :   if (r==1) { gel(V, idx) = Tp; return; }
    2029       82109 :   p2 = p>>1;
    2030       82109 :   T = Flx_get_red(Tp, p);
    2031       82109 :   XP = Flx_rem(XP, T, p);
    2032             :   while (1)
    2033        8065 :   {
    2034       90174 :     pari_sp btop = avma;
    2035             :     long i;
    2036       90174 :     GEN g = random_Flx(n, Tp[1], p);
    2037       90174 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2038       90174 :     if (lgpol(t) == 0) continue;
    2039      188839 :     for(i=1; i<=10; i++)
    2040             :     {
    2041      182738 :       pari_sp btop2 = avma;
    2042      182738 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    2043      182734 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    2044      182737 :       if (degpol(f) > 0 && degpol(f) < n) break;
    2045      100628 :       set_avma(btop2);
    2046             :     }
    2047       88210 :     if (degpol(f) > 0 && degpol(f) < n) break;
    2048        6101 :     set_avma(btop);
    2049             :   }
    2050       82109 :   f = Flx_normalize(f, p);
    2051       82109 :   ff = Flx_div(Tp, f ,p);
    2052       82109 :   Flx_edf_simple(f, XP, d, p, V, idx);
    2053       82109 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    2054             : }
    2055             : static void
    2056             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    2057             : 
    2058             : static void
    2059      604479 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    2060             : {
    2061             :   pari_sp av;
    2062      604479 :   GEN Tp = get_Flx_mod(T);
    2063      604478 :   long n = degpol(hp), vT = Tp[1];
    2064             :   GEN u1, u2, f1, f2;
    2065      604481 :   ulong p2 = p>>1;
    2066             :   GEN R, h;
    2067      604481 :   h = Flx_get_red(hp, p);
    2068      604455 :   t = Flx_rem(t, T, p);
    2069      604385 :   av = avma;
    2070             :   do
    2071             :   {
    2072      990329 :     set_avma(av);
    2073      990332 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    2074      990153 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    2075      990346 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2076      604434 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    2077      604409 :   f1 = Flx_normalize(f1, p);
    2078      604447 :   u2 = Flx_div(hp, u1, p);
    2079      604465 :   f2 = Flx_div(Tp, f1, p);
    2080      604457 :   if (degpol(u1)==1)
    2081             :   {
    2082      434047 :     if (degpol(f1)==d)
    2083      427696 :       gel(V, idx) = f1;
    2084             :     else
    2085        6339 :       Flx_edf(f1, XP, d, p, V, idx);
    2086             :   }
    2087             :   else
    2088      170466 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    2089      604498 :   idx += degpol(f1)/d;
    2090      604478 :   if (degpol(u2)==1)
    2091             :   {
    2092      432546 :     if (degpol(f2)==d)
    2093      426589 :       gel(V, idx) = f2;
    2094             :     else
    2095        5957 :       Flx_edf(f2, XP, d, p, V, idx);
    2096             :   }
    2097             :   else
    2098      171973 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    2099      604530 : }
    2100             : 
    2101             : static void
    2102      262113 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2103             : {
    2104      262113 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2105             :   GEN T, h, t;
    2106             :   pari_timer ti;
    2107      262111 :   if (r==1) { gel(V, idx) = Tp; return; }
    2108      262111 :   T = Flx_get_red(Tp, p);
    2109      262103 :   XP = Flx_rem(XP, T, p);
    2110      262096 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2111             :   do
    2112             :   {
    2113      268165 :     GEN g = random_Flx(n, vT, p);
    2114      268181 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2115      268187 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2116      268187 :     h = Flxq_minpoly(t, T, p);
    2117      268170 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2118      268170 :   } while (degpol(h) <= 1);
    2119      262104 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    2120             : }
    2121             : 
    2122             : static GEN
    2123      580462 : Flx_factor_Shoup(GEN T, ulong p)
    2124             : {
    2125      580462 :   long i, n, s = 0;
    2126             :   GEN XP, D, V;
    2127      580462 :   long e = expu(p);
    2128             :   pari_timer ti;
    2129      580462 :   n = get_Flx_degree(T);
    2130      580462 :   T = Flx_get_red(T, p);
    2131      580458 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2132      580458 :   XP = Flx_Frobenius(T, p);
    2133      580450 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2134      580450 :   D = Flx_ddf_Shoup(T, XP, p);
    2135      580481 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2136      580481 :   s = ddf_to_nbfact(D);
    2137      580474 :   V = cgetg(s+1, t_COL);
    2138     3811225 :   for (i = 1, s = 1; i <= n; i++)
    2139             :   {
    2140     3230757 :     GEN Di = gel(D,i);
    2141     3230757 :     long ni = degpol(Di), ri = ni/i;
    2142     3230760 :     if (ni == 0) continue;
    2143      733484 :     Di = Flx_normalize(Di, p);
    2144      733505 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2145      272219 :     if (ri <= e*expu(e))
    2146      249807 :       Flx_edf(Di, XP, i, p, V, s);
    2147             :     else
    2148       22413 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2149      272218 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2150      272218 :     s += ri;
    2151             :   }
    2152      580468 :   return V;
    2153             : }
    2154             : 
    2155             : static GEN
    2156      553059 : Flx_factor_Cantor(GEN T, ulong p)
    2157             : {
    2158      553059 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2159      553052 :   long i, j, l = lg(V);
    2160      553052 :   F = cgetg(l, t_VEC);
    2161      553055 :   E = cgetg(l, t_VEC);
    2162     1344193 :   for (i=1, j=1; i < l; i++)
    2163      791129 :     if (degpol(gel(V,i)))
    2164             :     {
    2165      580464 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2166      580468 :       gel(F, j) = Fj;
    2167      580468 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2168      580471 :       j++;
    2169             :     }
    2170      553064 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2171             : }
    2172             : 
    2173             : GEN
    2174           0 : Flx_ddf(GEN T, ulong p)
    2175             : {
    2176             :   GEN XP;
    2177           0 :   T = Flx_get_red(T, p);
    2178           0 :   XP = Flx_Frobenius(T, p);
    2179           0 :   return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
    2180             : }
    2181             : 
    2182             : static GEN
    2183      752563 : Flx_simplefact_Cantor(GEN T, ulong p)
    2184             : {
    2185             :   GEN V;
    2186             :   long i, l;
    2187      752563 :   T = Flx_get_red(T, p);
    2188      752563 :   V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
    2189     1510061 :   for (i=1; i < l; i++)
    2190      757498 :     gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
    2191      752563 :   return vddf_to_simplefact(V, get_Flx_degree(T));
    2192             : }
    2193             : 
    2194             : static int
    2195        1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2196             : {
    2197        1078 :   pari_sp av = avma;
    2198             :   pari_timer ti;
    2199             :   long n;
    2200        1078 :   GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
    2201        1078 :   if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    2202         791 :   n = get_Flx_degree(T);
    2203         791 :   T = Flx_get_red(Tp, p);
    2204         791 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2205         791 :   XP = Flx_Frobenius(T, p);
    2206         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2207         791 :   D = Flx_ddf_Shoup(T, XP, p);
    2208         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2209         791 :   return gc_bool(av, degpol(gel(D,n)) == n);
    2210             : }
    2211             : 
    2212             : /* f monic */
    2213             : static GEN
    2214     3013917 : Flx_factor_i(GEN f, ulong pp, long flag)
    2215             : {
    2216             :   long d;
    2217     3013917 :   if (pp==2) { /*We need to handle 2 specially */
    2218       34377 :     GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
    2219       34377 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2220       34377 :     return F;
    2221             :   }
    2222     2979540 :   d = degpol(f);
    2223     2979543 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2224     1306699 :   switch(flag)
    2225             :   {
    2226      553058 :     default: return Flx_factor_Cantor(f, pp);
    2227      752563 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2228        1078 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2229             :   }
    2230             : }
    2231             : 
    2232             : GEN
    2233     1934156 : Flx_degfact(GEN f, ulong p)
    2234             : {
    2235     1934156 :   pari_sp av = avma;
    2236     1934156 :   GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
    2237     1934156 :   return gerepilecopy(av, z);
    2238             : }
    2239             : 
    2240             : /* T must be squarefree mod p*/
    2241             : GEN
    2242      291130 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2243             : {
    2244             :   GEN XP, D;
    2245             :   pari_timer ti;
    2246      291130 :   long i, s, n = get_Flx_degree(T);
    2247      291130 :   GEN V = const_vecsmall(n, 0);
    2248      291130 :   pari_sp av = avma;
    2249      291130 :   T = Flx_get_red(T, p);
    2250      291130 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2251      291130 :   XP = Flx_Frobenius(T, p);
    2252      291130 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2253      291130 :   D = Flx_ddf_Shoup(T, XP, p);
    2254      291130 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2255     1733676 :   for (i = 1, s = 0; i <= n; i++)
    2256             :   {
    2257     1442546 :     V[i] = degpol(gel(D,i))/i;
    2258     1442546 :     s += V[i];
    2259             :   }
    2260      291130 :   *nb = s;
    2261      291130 :   set_avma(av); return V;
    2262             : }
    2263             : 
    2264             : long
    2265      134284 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2266             : {
    2267      134284 :   pari_sp av = avma;
    2268      134284 :   long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
    2269      134287 :   return gc_long(av,s);
    2270             : }
    2271             : 
    2272             : /* T must be squarefree mod p*/
    2273             : long
    2274      134281 : Flx_nbfact(GEN T, ulong p)
    2275             : {
    2276      134281 :   pari_sp av = avma;
    2277      134281 :   GEN XP = Flx_Frobenius(T, p);
    2278      134284 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2279      134287 :   return gc_long(av,n);
    2280             : }
    2281             : 
    2282             : int
    2283        1057 : Flx_is_irred(GEN f, ulong p)
    2284             : {
    2285        1057 :   pari_sp av = avma;
    2286        1057 :   f = Flx_normalize(f,p);
    2287        1057 :   return gc_bool(av, !!Flx_factor_i(f,p,2));
    2288             : }
    2289             : 
    2290             : /* Use this function when you think f is reducible, and that there are lots of
    2291             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2292             : int
    2293         105 : FpX_is_irred(GEN f, GEN p)
    2294             : {
    2295         105 :   pari_sp av = avma;
    2296             :   int z;
    2297         105 :   switch(ZX_factmod_init(&f,p))
    2298             :   {
    2299          21 :     case 0:  z = !!F2x_factor_i(f,2); break;
    2300          77 :     case 1:  z = !!Flx_factor_i(f,p[2],2); break;
    2301           7 :     default: z = !!FpX_factor_i(f,p,2); break;
    2302             :   }
    2303         105 :   return gc_bool(av,z);
    2304             : }
    2305             : GEN
    2306        1862 : FpX_degfact(GEN f, GEN p) {
    2307        1862 :   pari_sp av = avma;
    2308             :   GEN F;
    2309        1862 :   switch(ZX_factmod_init(&f,p))
    2310             :   {
    2311           7 :     case 0:  F = F2x_factor_i(f,1); break;
    2312        1827 :     case 1:  F = Flx_factor_i(f,p[2],1); break;
    2313          28 :     default: F = FpX_factor_i(f,p,1); break;
    2314             :   }
    2315        1862 :   return gerepilecopy(av, F);
    2316             : }
    2317             : 
    2318             : #if 0
    2319             : /* set x <-- x + c*y mod p */
    2320             : /* x is not required to be normalized.*/
    2321             : static void
    2322             : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2323             : {
    2324             :   long i, lx, ly;
    2325             :   ulong *x=(ulong *)gx;
    2326             :   ulong *y=(ulong *)gy;
    2327             :   if (!c) return;
    2328             :   lx = lg(gx);
    2329             :   ly = lg(gy);
    2330             :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2331             :   if (SMALL_ULONG(p))
    2332             :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2333             :   else
    2334             :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2335             : }
    2336             : #endif
    2337             : 
    2338             : GEN
    2339     1045772 : FpX_factor(GEN f, GEN p)
    2340             : {
    2341     1045772 :   pari_sp av = avma;
    2342             :   GEN F;
    2343     1045772 :   switch(ZX_factmod_init(&f, p))
    2344             :   {
    2345      171091 :     case 0:  F = F2x_factor_i(f,0);
    2346      171091 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    2347      873151 :     case 1:  F = Flx_factor_i(f,p[2],0);
    2348      873182 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    2349        1516 :     default: F = FpX_factor_i(f,p,0); break;
    2350             :   }
    2351     1045783 :   return gerepilecopy(av, F);
    2352             : }
    2353             : 
    2354             : GEN
    2355      203650 : Flx_factor(GEN f, ulong p)
    2356             : {
    2357      203650 :   pari_sp av = avma;
    2358      203650 :   return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
    2359             : }
    2360             : GEN
    2361       15289 : F2x_factor(GEN f)
    2362             : {
    2363       15289 :   pari_sp av = avma;
    2364       15289 :   return gerepilecopy(av, F2x_factor_i(f,0));
    2365             : }

Generated by: LCOV version 1.13