Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.11.0 lcov report (development 22860-5579deb0b) Lines: 1252 1351 92.7 %
Date: 2018-07-18 05:36:42 Functions: 111 119 93.3 %
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      874607 : ZX_factmod_init(GEN *F, GEN p)
      32             : {
      33      874607 :   if (lgefint(p) == 3)
      34             :   {
      35      873222 :     ulong pp = p[2];
      36      873222 :     if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
      37      728520 :     *F = ZX_to_Flx(*F, pp);
      38      728520 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      39      728520 :     return 1;
      40             :   }
      41        1385 :   *F = FpX_red(*F, p);
      42        1385 :   if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      43        1385 :   return 2;
      44             : }
      45             : static void
      46       83252 : ZX_rootmod_init(GEN *F, GEN p)
      47             : {
      48       83252 :   if (lgefint(p) == 3)
      49             :   {
      50       75280 :     ulong pp = p[2];
      51       75280 :     *F = ZX_to_Flx(*F, pp);
      52       75280 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      53             :   }
      54             :   else
      55             :   {
      56        7972 :     *F = FpX_red(*F, p);
      57        7972 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      58             :   }
      59       83252 : }
      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          28 :     for (i = 1; i < p; i++) r[i] = i;
      70             :   } else {
      71          14 :     r = cgetg(p+1, t_VECSMALL);
      72          14 :     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         126 : Flx_Xnm1(long sv, long n, ulong p)
      80             : {
      81         126 :   GEN t = cgetg(n+3, t_VECSMALL);
      82             :   long i;
      83         126 :   t[1] = sv;
      84         126 :   t[2] = p - 1;
      85         126 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      86         126 :   t[i] = 1; return t;
      87             : }
      88             : /* X^n + 1 */
      89             : static GEN
      90         140 : Flx_Xn1(long sv, long n, ulong p)
      91             : {
      92         140 :   GEN t = cgetg(n+3, t_VECSMALL);
      93             :   long i;
      94             :   (void) p;
      95         140 :   t[1] = sv;
      96         140 :   t[2] = 1;
      97         140 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      98         140 :   t[i] = 1; return t;
      99             : }
     100             : 
     101             : static ulong
     102          28 : Fl_nonsquare(ulong p)
     103             : {
     104          28 :   long k = 2;
     105           7 :   for (;; k++)
     106           7 :   {
     107          35 :     long i = krouu(k, p);
     108          35 :     if (!i) pari_err_PRIME("Fl_nonsquare",utoipos(p));
     109          63 :     if (i < 0) return k;
     110             :   }
     111             : }
     112             : 
     113             : static GEN
     114        8881 : Flx_root_mod_2(GEN f)
     115             : {
     116        8881 :   int z1, z0 = !(f[2] & 1);
     117             :   long i,n;
     118             :   GEN y;
     119             : 
     120        8881 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     121        8881 :   z1 = n & 1;
     122        8881 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     123        8881 :   if (z0) y[i++] = 0;
     124        8881 :   if (z1) y[i  ] = 1;
     125        8881 :   return y;
     126             : }
     127             : static ulong
     128          56 : Flx_oneroot_mod_2(GEN f)
     129             : {
     130             :   long i,n;
     131          56 :   if (!(f[2] & 1)) return 0;
     132          56 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     133          56 :   if (n & 1) return 1;
     134          28 :   return 2;
     135             : }
     136             : 
     137             : static GEN FpX_roots_i(GEN f, GEN p);
     138             : static GEN Flx_roots_i(GEN f, ulong p);
     139             : 
     140             : static int
     141     3212314 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     142             : 
     143             : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
     144             :  * pp is a small prime.
     145             :  * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
     146             :  * assume that f is an FpX, pp a prime and return t_INTs */
     147             : static GEN
     148       70939 : rootmod_aux(GEN f, GEN pp)
     149             : {
     150             :   GEN y;
     151       70939 :   switch(lg(f))
     152             :   {
     153          14 :     case 2: pari_err_ROOTS0("rootmod");
     154          49 :     case 3: return cgetg(1,t_COL);
     155             :   }
     156       70876 :   if (typ(f) == t_VECSMALL)
     157             :   {
     158       67766 :     ulong p = pp[2];
     159       67766 :     if (p == 2)
     160        8881 :       y = Flx_root_mod_2(f);
     161             :     else
     162             :     {
     163       58885 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     164       58885 :       y = Flx_roots_i(f, p);
     165             :     }
     166       67759 :     y = Flc_to_ZC(y);
     167             :   }
     168             :   else
     169        3110 :     y = FpX_roots_i(f, pp);
     170       70862 :   return y;
     171             : }
     172             : /* assume that f is a ZX and p a prime */
     173             : GEN
     174       70939 : FpX_roots(GEN f, GEN p)
     175             : {
     176       70939 :   pari_sp av = avma;
     177       70939 :   GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
     178       70911 :   return gerepileupto(av, y);
     179             : }
     180             : 
     181             : /* assume x reduced mod p > 2, monic. */
     182             : static int
     183          21 : FpX_quad_factortype(GEN x, GEN p)
     184             : {
     185          21 :   GEN b = gel(x,3), c = gel(x,2);
     186          21 :   GEN D = subii(sqri(b), shifti(c,2));
     187          21 :   return kronecker(D,p);
     188             : }
     189             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     190             : static GEN
     191        7532 : FpX_quad_root(GEN x, GEN p, int unknown)
     192             : {
     193        7532 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     194             : 
     195        7532 :   if (absequaliu(p, 2)) {
     196           0 :     if (!signe(b)) return c;
     197           0 :     return signe(c)? NULL: gen_1;
     198             :   }
     199        7532 :   D = subii(sqri(b), shifti(c,2));
     200        7532 :   D = remii(D,p);
     201        7532 :   if (unknown && kronecker(D,p) == -1) return NULL;
     202             : 
     203        7019 :   s = Fp_sqrt(D,p);
     204             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     205        7019 :   if (!s) return NULL;
     206        7011 :   return Fp_halve(Fp_sub(s,b, p), p);
     207             : }
     208             : static GEN
     209        3146 : FpX_otherroot(GEN x, GEN r, GEN p)
     210        3146 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     211             : 
     212             : /* disc(x^2+bx+c) = b^2 - 4c */
     213             : static ulong
     214    21200559 : Fl_disc_bc(ulong b, ulong c, ulong p)
     215    21200559 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     216             : /* p > 2 */
     217             : static ulong
     218    20466578 : Flx_quad_root(GEN x, ulong p, int unknown)
     219             : {
     220    20466578 :   ulong s, b = x[3], c = x[2];
     221    20466578 :   ulong D = Fl_disc_bc(b, c, p);
     222    20462537 :   if (unknown && krouu(D,p) == -1) return p;
     223    13692784 :   s = Fl_sqrt(D,p);
     224    13686700 :   if (s==~0UL) return p;
     225    13686687 :   return Fl_halve(Fl_sub(s,b, p), p);
     226             : }
     227             : static ulong
     228    12262873 : Flx_otherroot(GEN x, ulong r, ulong p)
     229    12262873 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     230             : 
     231             : 
     232             : /* 'todo' contains the list of factors to be split.
     233             :  * 'done' the list of finished factors, no longer touched */
     234             : struct split_t { GEN todo, done; };
     235             : static void
     236     4930559 : split_init(struct split_t *S, long max)
     237             : {
     238     4930559 :   S->todo = vectrunc_init(max);
     239     4930555 :   S->done = vectrunc_init(max);
     240     4930409 : }
     241             : #if 0
     242             : /* move todo[i] to done */
     243             : static void
     244             : split_convert(struct split_t *S, long i)
     245             : {
     246             :   long n = lg(S->todo)-1;
     247             :   vectrunc_append(S->done, gel(S->todo,i));
     248             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     249             :   setlg(S->todo, n);
     250             : }
     251             : #endif
     252             : /* append t to todo */
     253             : static void
     254     5177019 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     255             : /* delete todo[i], add t to done */
     256             : static void
     257     5177052 : split_moveto_done(struct split_t *S, long i, GEN t)
     258             : {
     259     5177052 :   long n = lg(S->todo)-1;
     260     5177052 :   vectrunc_append(S->done, t);
     261     5177023 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     262     5177023 :   setlg(S->todo, n);
     263             : 
     264     5177561 : }
     265             : /* append t to done */
     266             : static void
     267      391697 : split_add_done(struct split_t *S, GEN t)
     268      391697 : { vectrunc_append(S->done, t); }
     269             : /* split todo[i] into a and b */
     270             : static void
     271      337162 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     272             : {
     273      337162 :   gel(S->todo, i) = a;
     274      337162 :   split_add(S, b);
     275      337161 : }
     276             : /* split todo[i] into a and b, moved to done */
     277             : static void
     278      372319 : split_done(struct split_t *S, long i, GEN a, GEN b)
     279             : {
     280      372319 :   split_moveto_done(S, i, a);
     281      372321 :   split_add_done(S, b);
     282      372322 : }
     283             : 
     284             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     285             : static GEN
     286        3110 : FpX_roots_i(GEN f, GEN p)
     287             : {
     288             :   GEN pol, pol0, a, q;
     289             :   struct split_t S;
     290             : 
     291        3110 :   split_init(&S, lg(f)-1);
     292        3110 :   settyp(S.done, t_COL);
     293        3110 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     294        3110 :   switch(degpol(f))
     295             :   {
     296           7 :     case 0: return ZC_copy(S.done);
     297          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     298             :     case 2: {
     299        1743 :       GEN s, r = FpX_quad_root(f, p, 1);
     300        1743 :       if (r) {
     301        1743 :         split_add_done(&S, r);
     302        1743 :         s = FpX_otherroot(f,r, p);
     303             :         /* f not known to be square free yet */
     304        1743 :         if (!equalii(r, s)) split_add_done(&S, s);
     305             :       }
     306        1743 :       return sort(S.done);
     307             :     }
     308             :   }
     309             : 
     310        1346 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     311        1346 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     312        1346 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     313        1346 :   a = FpX_gcd(f,a, p);
     314        1346 :   if (!degpol(a)) return ZC_copy(S.done);
     315        1346 :   split_add(&S, FpX_normalize(a,p));
     316             : 
     317        1346 :   q = shifti(p,-1);
     318        1346 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     319        1346 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     320        2819 :   for (pol0[2] = 1;; pol0[2]++)
     321        1473 :   {
     322        2819 :     long j, l = lg(S.todo);
     323        2819 :     if (l == 1) return sort(S.done);
     324        1480 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     325        3074 :     for (j = 1; j < l; j++)
     326             :     {
     327        1601 :       GEN b, r, s, c = gel(S.todo,j);
     328        1601 :       switch(degpol(c))
     329             :       { /* convert linear and quadratics to roots, try to split the rest */
     330             :         case 1:
     331         195 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     332         195 :           j--; l--; break;
     333             :         case 2:
     334        1272 :           r = FpX_quad_root(c, p, 0);
     335        1272 :           if (!r) pari_err_PRIME("polrootsmod",p);
     336        1265 :           s = FpX_otherroot(c,r, p);
     337        1265 :           split_done(&S, j, r, s);
     338        1265 :           j--; l--; break;
     339             :         default:
     340         134 :           b = FpXQ_pow(pol,q, c,p);
     341         134 :           if (degpol(b) <= 0) continue;
     342         121 :           b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
     343         121 :           if (!degpol(b)) continue;
     344         121 :           b = FpX_normalize(b, p);
     345         121 :           c = FpX_div(c,b, p);
     346         121 :           split_todo(&S, j, b, c);
     347             :       }
     348             :     }
     349             :   }
     350             : }
     351             : 
     352             : /* Assume f is normalized */
     353             : static ulong
     354      174315 : Flx_cubic_root(GEN ff, ulong p)
     355             : {
     356      174315 :   GEN f = Flx_normalize(ff,p);
     357      174316 :   ulong pi = get_Fl_red(p);
     358      174317 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     359      174317 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     360      174318 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     361      174316 :   ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
     362      174318 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     363      174318 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     364      174317 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     365      174317 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     366      174315 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     367      174317 :   if (s!=~0UL)
     368             :   {
     369      109726 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     370      109726 :     if (p%3==2) /* 1 solutions */
     371       19356 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     372             :     else
     373             :     {
     374       90370 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     375       90370 :       if (vS1==~0UL) return p; /*0 solutions*/
     376             :       /*3 solutions*/
     377             :     }
     378       79212 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     379       79211 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     380             :   }
     381             :   else
     382             :   {
     383       64591 :     pari_sp av = avma;
     384       64591 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     385       64591 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     386             :     ulong Sa;
     387       64592 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     388       64592 :     Sa = vS1[1];
     389       64592 :     if (p%3==1) /*1 solutions*/
     390             :     {
     391       24086 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     392       24086 :       if (Fa!=P)
     393       16026 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     394             :     }
     395       64592 :     avma = av;
     396       64592 :     return Fl_sub(Fl_double(Sa,p),t,p);
     397             :   }
     398             : }
     399             : 
     400             : /* assume p > 2 prime */
     401             : static ulong
     402     3133256 : Flx_oneroot_i(GEN f, ulong p, long fl)
     403             : {
     404             :   GEN pol, a;
     405             :   ulong q;
     406             :   long da;
     407             : 
     408     3133256 :   if (Flx_val(f)) return 0;
     409     3132534 :   switch(degpol(f))
     410             :   {
     411       12188 :     case 1: return Fl_neg(f[2], p);
     412     2683232 :     case 2: return Flx_quad_root(f, p, 1);
     413      158508 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     414             :   }
     415             : 
     416      278605 :   if (!fl)
     417             :   {
     418      247223 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     419      247223 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     420      247223 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     421      247223 :     a = Flx_gcd(f,a, p);
     422       31382 :   } else a = f;
     423      278605 :   da = degpol(a);
     424      278598 :   if (!da) return p;
     425      199088 :   a = Flx_normalize(a,p);
     426             : 
     427      199099 :   q = p >> 1;
     428      199099 :   pol = polx_Flx(f[1]);
     429      299221 :   for(pol[2] = 1;; pol[2]++)
     430             :   {
     431      399345 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     432      299233 :     switch(da)
     433             :     {
     434      124300 :       case 1: return Fl_neg(a[2], p);
     435       59023 :       case 2: return Flx_quad_root(a, p, 0);
     436       15814 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     437             :       default: {
     438      100096 :         GEN b = Flxq_powu(pol,q, a,p);
     439             :         long db;
     440      100120 :         if (degpol(b) <= 0) continue;
     441       95038 :         b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
     442       95052 :         db = degpol(b); if (!db) continue;
     443       95045 :         b = Flx_normalize(b, p);
     444       95053 :         if (db <= (da >> 1)) {
     445       57911 :           a = b;
     446       57911 :           da = db;
     447             :         } else {
     448       37142 :           a = Flx_div(a,b, p);
     449       37139 :           da -= db;
     450             :         }
     451             :       }
     452             :     }
     453             :   }
     454             : }
     455             : 
     456             : /* assume p > 2 prime */
     457             : static GEN
     458        4848 : FpX_oneroot_i(GEN f, GEN p)
     459             : {
     460             :   GEN pol, pol0, a, q;
     461             :   long da;
     462             : 
     463        4848 :   if (ZX_val(f)) return gen_0;
     464        4582 :   switch(degpol(f))
     465             :   {
     466         675 :     case 1: return subii(p, gel(f,2));
     467        3837 :     case 2: return FpX_quad_root(f, p, 1);
     468             :   }
     469             : 
     470          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     471          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     472          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     473          70 :   a = FpX_gcd(f,a, p);
     474          70 :   da = degpol(a);
     475          70 :   if (!da) return NULL;
     476          70 :   a = FpX_normalize(a,p);
     477             : 
     478          70 :   q = shifti(p,-1);
     479          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     480          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     481         224 :   for (pol0[2]=1; ; pol0[2]++)
     482             :   {
     483         378 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     484         224 :     switch(da)
     485             :     {
     486          42 :       case 1: return subii(p, gel(a,2));
     487          28 :       case 2: return FpX_quad_root(a, p, 0);
     488             :       default: {
     489         154 :         GEN b = FpXQ_pow(pol,q, a,p);
     490             :         long db;
     491         154 :         if (degpol(b) <= 0) continue;
     492         147 :         b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
     493         147 :         db = degpol(b); if (!db) continue;
     494         147 :         b = FpX_normalize(b, p);
     495         147 :         if (db <= (da >> 1)) {
     496         105 :           a = b;
     497         105 :           da = db;
     498             :         } else {
     499          42 :           a = FpX_div(a,b, p);
     500          42 :           da -= db;
     501             :         }
     502             :       }
     503             :     }
     504             :   }
     505             : }
     506             : 
     507             : ulong
     508     3093602 : Flx_oneroot(GEN f, ulong p)
     509             : {
     510     3093602 :   pari_sp av = avma;
     511             :   ulong r;
     512     3093602 :   switch(lg(f))
     513             :   {
     514           0 :     case 2: return 0;
     515           0 :     case 3: avma = av; return p;
     516             :   }
     517     3093602 :   if (p == 2) return Flx_oneroot_mod_2(f);
     518     3093602 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     519     3093602 :   avma = av; return r;
     520             : }
     521             : 
     522             : ulong
     523       32245 : Flx_oneroot_split(GEN f, ulong p)
     524             : {
     525       32245 :   pari_sp av = avma;
     526             :   ulong r;
     527       32245 :   switch(lg(f))
     528             :   {
     529           0 :     case 2: return 0;
     530           0 :     case 3: avma = av; return p;
     531             :   }
     532       32245 :   if (p == 2) return Flx_oneroot_mod_2(f);
     533       32245 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     534       32283 :   avma = av; return r;
     535             : }
     536             : 
     537             : /* assume that p is prime */
     538             : GEN
     539       12313 : FpX_oneroot(GEN f, GEN pp) {
     540       12313 :   pari_sp av = avma;
     541       12313 :   ZX_rootmod_init(&f, pp);
     542       12313 :   switch(lg(f))
     543             :   {
     544           0 :     case 2: avma = av; return gen_0;
     545           0 :     case 3: avma = av; return NULL;
     546             :   }
     547       12313 :   if (typ(f) == t_VECSMALL)
     548             :   {
     549        7465 :     ulong r, p = pp[2];
     550        7465 :     if (p == 2)
     551          56 :       r = Flx_oneroot_mod_2(f);
     552             :     else
     553        7409 :       r = Flx_oneroot_i(f, p, 0);
     554        7465 :     avma = av;
     555        7465 :     return (r == p)? NULL: utoi(r);
     556             :   }
     557        4848 :   f = FpX_oneroot_i(f, pp);
     558        4848 :   if (!f) { avma = av; return NULL; }
     559        4848 :   return gerepileuptoint(av, f);
     560             : }
     561             : 
     562             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     563             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     564             : /* returned in n                                                          */
     565             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     566             : /* twice as big and decrement until it divides p-1.                       */
     567             : static GEN
     568         217 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     569             : {
     570         217 :    pari_sp ltop = avma;
     571             :    GEN pm, factn, power, base, zeta;
     572             :    long n;
     573             : 
     574         217 :    pm = subis (p, 1ul);
     575         217 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     576         217 :    factn = Z_factor(stoi(n));
     577         217 :    power = diviuexact (pm, n);
     578         217 :    base = gen_1;
     579             :    do {
     580         378 :       base = addis (base, 1l);
     581         378 :       zeta = Fp_pow (base, power, p);
     582             :    }
     583         378 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     584         217 :    *pt_n = n;
     585         217 :    return gerepileuptoint (ltop, zeta);
     586             : }
     587             : 
     588             : GEN
     589         924 : FpX_oneroot_split(GEN fact, GEN p)
     590             : {
     591         924 :   pari_sp av = avma;
     592             :   long n, deg_f, i, dmin;
     593             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     594         924 :   fact = FpX_normalize(fact, p);
     595         924 :   deg_f = degpol(fact);
     596         924 :   if (deg_f<=2) return FpX_oneroot(fact, p);
     597         217 :   minfactor = fact; /* factor of minimal degree found so far */
     598         217 :   dmin = degpol(minfactor);
     599         217 :   prim = good_root_of_unity(p, deg_f, 1, &n);
     600         217 :   expo = diviuexact(subiu(p, 1), n);
     601         217 :   xplusa = pol_x(varn(fact));
     602         217 :   zeta = gen_1;
     603        1071 :   while (dmin != 1)
     604             :   {
     605             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     606             :     /* zeta varies over the roots of unity in F_p                         */
     607         637 :     fact = minfactor; deg_f = dmin;
     608             :     /* update X+a, avoid a=0 */
     609         637 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     610         637 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     611        1260 :     for (i = 0; i < n; i++)
     612             :     {
     613         980 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     614         980 :       long dtmp = degpol(tmp);
     615         980 :       if (dtmp > 0 && dtmp < deg_f)
     616             :       {
     617         427 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     618         427 :         if (dtmp < dmin)
     619             :         {
     620         427 :           minfactor = FpX_normalize (tmp, p);
     621         427 :           dmin = dtmp;
     622         427 :           if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
     623             :             /* stop early to avoid too many gcds */
     624             :             break;
     625             :         }
     626             :       }
     627         623 :       zeta = Fp_mul (zeta, prim, p);
     628             :     }
     629             :   }
     630         217 :   return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
     631             : }
     632             : 
     633             : /*******************************************************************/
     634             : /*                                                                 */
     635             : /*                     FACTORISATION MODULO p                      */
     636             : /*                                                                 */
     637             : /*******************************************************************/
     638             : 
     639             : /* F / E  a vector of vectors of factors / exponents of virtual length l
     640             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
     641             : static GEN
     642      482087 : FE_concat(GEN F, GEN E, long l)
     643             : {
     644      482087 :   setlg(E,l); E = shallowconcat1(E);
     645      482087 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
     646             : }
     647             : 
     648             : static GEN
     649          14 : ddf_to_ddf2_i(GEN V, long fl)
     650             : {
     651             :   GEN F, D;
     652          14 :   long i, j, l = lg(V);
     653          14 :   F = cgetg(l, t_VEC);
     654          14 :   D = cgetg(l, t_VECSMALL);
     655         112 :   for (i = j = 1; i < l; i++)
     656             :   {
     657          98 :     GEN Vi = gel(V,i);
     658          98 :     if ((fl==2 && F2x_degree(Vi) == 0)
     659          84 :       ||(fl==0 && degpol(Vi) == 0)) continue;
     660          35 :     gel(F,j) = Vi;
     661          35 :     uel(D,j) = i; j++;
     662             :   }
     663          14 :   setlg(F,j);
     664          14 :   setlg(D,j); return mkvec2(F,D);
     665             : }
     666             : 
     667             : GEN
     668           7 : ddf_to_ddf2(GEN V)
     669           7 : { return ddf_to_ddf2_i(V, 0); }
     670             : 
     671             : static GEN
     672           7 : F2x_ddf_to_ddf2(GEN V)
     673           7 : { return ddf_to_ddf2_i(V, 2); }
     674             : 
     675             : GEN
     676      354779 : vddf_to_simplefact(GEN V, long d)
     677             : {
     678             :   GEN E, F;
     679      354779 :   long i, j, c, l = lg(V);
     680      354779 :   F = cgetg(d+1, t_VECSMALL);
     681      354779 :   E = cgetg(d+1, t_VECSMALL);
     682      712611 :   for (i = c = 1; i < l; i++)
     683             :   {
     684      357832 :     GEN Vi = gel(V,i);
     685      357832 :     long l = lg(Vi);
     686     2463663 :     for (j = 1; j < l; j++)
     687             :     {
     688     2105831 :       long k, n = degpol(gel(Vi,j)) / j;
     689     2105831 :       for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
     690             :     }
     691             :   }
     692      354779 :   setlg(F,c);
     693      354779 :   setlg(E,c);
     694      354779 :   return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
     695             : }
     696             : 
     697             : /* product of terms of degree 1 in factorization of f */
     698             : GEN
     699      142823 : FpX_split_part(GEN f, GEN p)
     700             : {
     701      142823 :   long n = degpol(f);
     702      142823 :   GEN z, X = pol_x(varn(f));
     703      142823 :   if (n <= 1) return f;
     704      142508 :   f = FpX_red(f, p);
     705      142508 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     706      142508 :   return FpX_gcd(z,f,p);
     707             : }
     708             : 
     709             : /* Compute the number of roots in Fp without counting multiplicity
     710             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     711             : long
     712       99148 : FpX_nbroots(GEN f, GEN p)
     713             : {
     714       99148 :   pari_sp av = avma;
     715       99148 :   GEN z = FpX_split_part(f, p);
     716       99148 :   avma = av; return degpol(z);
     717             : }
     718             : 
     719             : int
     720           0 : FpX_is_totally_split(GEN f, GEN p)
     721             : {
     722           0 :   long n=degpol(f);
     723           0 :   pari_sp av = avma;
     724           0 :   if (n <= 1) return 1;
     725           0 :   if (abscmpui(n, p) > 0) return 0;
     726           0 :   f = FpX_red(f, p);
     727           0 :   avma = av; return gequalX(FpX_Frobenius(f, p));
     728             : }
     729             : 
     730             : long
     731     2406352 : Flx_nbroots(GEN f, ulong p)
     732             : {
     733     2406352 :   long n = degpol(f);
     734     2406352 :   pari_sp av = avma;
     735             :   GEN z;
     736     2406352 :   if (n <= 1) return n;
     737     2404147 :   if (n == 2)
     738             :   {
     739             :     ulong D;
     740       11844 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     741       10962 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     742       10962 :     return 1 + krouu(D,p);
     743             :   }
     744     2392303 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     745     2392303 :   z = Flx_gcd(z, f, p);
     746     2392303 :   avma = av; return degpol(z);
     747             : }
     748             : 
     749             : long
     750        4389 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
     751             : {
     752        4389 :   pari_sp av = avma;
     753             :   GEN X, b, g, xq;
     754             :   long i, j, n, v, B, l, m;
     755             :   pari_timer ti;
     756             :   hashtable h;
     757             : 
     758        4389 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     759        4389 :   X = pol_x(v);
     760        4389 :   if (ZX_equal(X,XP)) return 1;
     761        4389 :   B = n/2;
     762        4389 :   l = usqrt(B);
     763        4389 :   m = (B+l-1)/l;
     764        4389 :   T = FpX_get_red(T, p);
     765        4389 :   hash_init_GEN(&h, l+2, ZX_equal, 1);
     766        4389 :   hash_insert_long(&h, X,  0);
     767        4389 :   hash_insert_long(&h, XP, 1);
     768        4389 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     769        4389 :   b = XP;
     770        4389 :   xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
     771        4389 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
     772       10367 :   for (i = 3; i <= l+1; i++)
     773             :   {
     774        6776 :     b = FpX_FpXQV_eval(b, xq, T, p);
     775        6776 :     if (gequalX(b)) { avma = av; return i-1; }
     776        5978 :     hash_insert_long(&h, b, i-1);
     777             :   }
     778        3591 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
     779        3591 :   g = b;
     780        3591 :   xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
     781        3591 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
     782       12278 :   for(i = 2; i <= m+1; i++)
     783             :   {
     784       10696 :     g = FpX_FpXQV_eval(g, xq, T, p);
     785       10696 :     if (hash_haskey_long(&h, g, &j)) { avma=av; return l*i-j; }
     786             :   }
     787        1582 :   avma = av; return n;
     788             : }
     789             : 
     790             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     791             : static GEN
     792         669 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
     793             : {
     794             :   GEN b, g, h, F, f, Tr, xq;
     795             :   long i, j, n, v, B, l, m;
     796             :   pari_timer ti;
     797             : 
     798         669 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     799         669 :   if (n == 0) return cgetg(1, t_VEC);
     800         669 :   if (n == 1) return mkvec(get_FpX_mod(T));
     801         623 :   B = n/2;
     802         623 :   l = usqrt(B);
     803         623 :   m = (B+l-1)/l;
     804         623 :   T = FpX_get_red(T, p);
     805         623 :   b = cgetg(l+2, t_VEC);
     806         623 :   gel(b, 1) = pol_x(v);
     807         623 :   gel(b, 2) = XP;
     808         623 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     809         623 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     810         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
     811         811 :   for (i = 3; i <= l+1; i++)
     812         188 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     813         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
     814         623 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     815         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
     816         623 :   g = cgetg(m+1, t_VEC);
     817         623 :   gel(g, 1) = gel(xq, 2);
     818         623 :   for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     819         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
     820         623 :   h = cgetg(m+1, t_VEC);
     821        1945 :   for (j = 1; j <= m; j++)
     822             :   {
     823        1322 :     pari_sp av = avma;
     824        1322 :     GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
     825        1322 :     for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
     826        1322 :     gel(h,j) = gerepileupto(av, e);
     827             :   }
     828         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
     829         623 :   Tr = get_FpX_mod(T);
     830         623 :   F = cgetg(m+1, t_VEC);
     831        1945 :   for (j = 1; j <= m; j++)
     832             :   {
     833        1322 :     GEN u = FpX_gcd(Tr, gel(h,j), p);
     834        1322 :     if (degpol(u))
     835             :     {
     836         226 :       u = FpX_normalize(u, p);
     837         226 :       Tr = FpX_div(Tr, u, p);
     838             :     }
     839        1322 :     gel(F,j) = u;
     840             :   }
     841         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
     842         623 :   f = const_vec(n, pol_1(v));
     843        1945 :   for (j = 1; j <= m; j++)
     844             :   {
     845        1322 :     GEN e = gel(F, j);
     846        1357 :     for (i=l-1; i >= 0; i--)
     847             :     {
     848        1357 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     849        1357 :       if (degpol(u))
     850             :       {
     851         233 :         u = FpX_normalize(u, p);
     852         233 :         gel(f, l*j-i) = u;
     853         233 :         e = FpX_div(e, u, p);
     854             :       }
     855        1357 :       if (!degpol(e)) break;
     856             :     }
     857             :   }
     858         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
     859         623 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     860         623 :   return f;
     861             : }
     862             : 
     863             : static void
     864           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     865             : {
     866           0 :   long n = degpol(Tp), r = n/d, ct = 0;
     867             :   GEN T, f, ff, p2;
     868           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     869           0 :   p2 = shifti(p,-1);
     870           0 :   T = FpX_get_red(Tp, p);
     871           0 :   XP = FpX_rem(XP, T, p);
     872             :   while (1)
     873           0 :   {
     874           0 :     pari_sp btop = avma;
     875             :     long i;
     876           0 :     GEN g = random_FpX(n, varn(Tp), p);
     877           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     878           0 :     if (signe(t) == 0) continue;
     879           0 :     for(i=1; i<=10; i++)
     880             :     {
     881           0 :       pari_sp btop2 = avma;
     882           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     883           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     884           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     885           0 :       avma = btop2;
     886             :     }
     887           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     888           0 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
     889           0 :     avma = btop;
     890             :   }
     891           0 :   f = FpX_normalize(f, p);
     892           0 :   ff = FpX_div(Tp, f ,p);
     893           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     894           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     895             : }
     896             : 
     897             : static void
     898         462 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     899             : {
     900             :   pari_sp av;
     901         462 :   GEN Tp = get_FpX_mod(T);
     902         462 :   long n = degpol(hp), vT = varn(Tp), ct = 0;
     903             :   GEN u1, u2, f1, f2, R, h;
     904         462 :   h = FpX_get_red(hp, p);
     905         462 :   t = FpX_rem(t, T, p);
     906         462 :   av = avma;
     907             :   do
     908             :   {
     909         761 :     avma = av;
     910         761 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     911         761 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     912         761 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
     913         761 :   } while (degpol(u1)==0 || degpol(u1)==n);
     914         462 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     915         462 :   f1 = FpX_normalize(f1, p);
     916         462 :   u2 = FpX_div(hp, u1, p);
     917         462 :   f2 = FpX_div(Tp, f1, p);
     918         462 :   if (degpol(u1)==1)
     919         336 :     gel(V, idx) = f1;
     920             :   else
     921         126 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     922         462 :   idx += degpol(f1)/d;
     923         462 :   if (degpol(u2)==1)
     924         350 :     gel(V, idx) = f2;
     925             :   else
     926         112 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     927         462 : }
     928             : 
     929             : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
     930             : static void
     931         224 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     932             : {
     933         224 :   long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
     934             :   GEN T, h, t;
     935             :   pari_timer ti;
     936             : 
     937         224 :   T = FpX_get_red(Tp, p);
     938         224 :   XP = FpX_rem(XP, T, p);
     939         224 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     940             :   do
     941             :   {
     942         224 :     GEN g = random_FpX(n, vT, p);
     943         224 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     944         224 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     945         224 :     h = FpXQ_minpoly(t, T, p);
     946         224 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     947         224 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
     948         224 :   } while (degpol(h) != r);
     949         224 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     950         224 : }
     951             : 
     952             : static GEN
     953         655 : FpX_factor_Shoup(GEN T, GEN p)
     954             : {
     955         655 :   long i, n, s = 0;
     956             :   GEN XP, D, V;
     957         655 :   long e = expi(p);
     958             :   pari_timer ti;
     959         655 :   n = get_FpX_degree(T);
     960         655 :   T = FpX_get_red(T, p);
     961         655 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     962         655 :   XP = FpX_Frobenius(T, p);
     963         655 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     964         655 :   D = FpX_ddf_Shoup(T, XP, p);
     965         655 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
     966         655 :   s = ddf_to_nbfact(D);
     967         655 :   V = cgetg(s+1, t_COL);
     968        5507 :   for (i = 1, s = 1; i <= n; i++)
     969             :   {
     970        4852 :     GEN Di = gel(D,i);
     971        4852 :     long ni = degpol(Di), ri = ni/i;
     972        4852 :     if (ni == 0) continue;
     973         664 :     Di = FpX_normalize(Di, p);
     974         664 :     if (ni == i) { gel(V, s++) = Di; continue; }
     975         224 :     if (ri <= e*expu(e))
     976         224 :       FpX_edf(Di, XP, i, p, V, s);
     977             :     else
     978           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     979         224 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     980         224 :     s += ri;
     981             :   }
     982         655 :   return V;
     983             : }
     984             : 
     985             : long
     986      536358 : ddf_to_nbfact(GEN D)
     987             : {
     988      536358 :   long l = lg(D), i, s = 0;
     989      536358 :   for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
     990      536358 :   return s;
     991             : }
     992             : 
     993             : /* Yun algorithm: Assume p > degpol(T) */
     994             : static GEN
     995         672 : FpX_factor_Yun(GEN T, GEN p)
     996             : {
     997         672 :   long n = degpol(T), i = 1;
     998         672 :   GEN a, b, c, d = FpX_deriv(T, p);
     999         672 :   GEN V = cgetg(n+1,t_VEC);
    1000         672 :   a = FpX_gcd(T, d, p);
    1001         672 :   if (degpol(a) == 0) return mkvec(T);
    1002         488 :   b = FpX_div(T, a, p);
    1003             :   do
    1004             :   {
    1005        2432 :     c = FpX_div(d, a, p);
    1006        2432 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1007        2432 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1008        2432 :     gel(V, i++) = a;
    1009        2432 :     b = FpX_div(b, a, p);
    1010        2432 :   } while (degpol(b));
    1011         488 :   setlg(V, i); return V;
    1012             : }
    1013             : GEN
    1014         777 : FpX_factor_squarefree(GEN T, GEN p)
    1015             : {
    1016         777 :   if (lgefint(p)==3)
    1017             :   {
    1018         777 :     ulong pp = (ulong)p[2];
    1019         777 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1020         777 :     return FlxV_to_ZXV(u);
    1021             :   }
    1022           0 :   return FpX_factor_Yun(T, p);
    1023             : }
    1024             : 
    1025             : long
    1026         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1027             : {
    1028         168 :   pari_sp av = avma;
    1029             :   GEN lc, F;
    1030         168 :   long i, l, n = degpol(f), v = varn(f);
    1031         168 :   if (n % k) return 0;
    1032         168 :   if (lgefint(p)==3)
    1033             :   {
    1034         126 :     ulong pp = p[2];
    1035         126 :     GEN fp = ZX_to_Flx(f, pp);
    1036         126 :     if (!Flx_ispower(fp, k, pp, pt_r)) { avma = av; return 0; }
    1037         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else avma = av;
    1038         105 :     return 1;
    1039             :   }
    1040          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1041          42 :   if (!lc) { av = avma; return 0; }
    1042          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1043        1491 :   for(i=1; i <= l; i++)
    1044        1456 :     if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
    1045          35 :   if (pt_r)
    1046             :   {
    1047          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1048        1484 :     for (i=l; i>=1; i--)
    1049             :     {
    1050        1449 :       if (i%k) continue;
    1051         294 :       s = FpX_mul(s, gel(F,i), p);
    1052         294 :       r = FpX_mul(r, s, p);
    1053             :     }
    1054          35 :     *pt_r = gerepileupto(av, r);
    1055           0 :   } else av = avma;
    1056          35 :   return 1;
    1057             : }
    1058             : 
    1059             : static GEN
    1060         616 : FpX_factor_Cantor(GEN T, GEN p)
    1061             : {
    1062         616 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1063         616 :   long i, j, l = lg(V);
    1064         616 :   F = cgetg(l, t_VEC);
    1065         616 :   E = cgetg(l, t_VEC);
    1066        1748 :   for (i=1, j=1; i < l; i++)
    1067        1132 :     if (degpol(gel(V,i)))
    1068             :     {
    1069         655 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1070         655 :       gel(F, j) = Fj;
    1071         655 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1072         655 :       j++;
    1073             :     }
    1074         616 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1075             : }
    1076             : 
    1077             : static GEN
    1078           0 : FpX_ddf_i(GEN T, GEN p)
    1079             : {
    1080             :   GEN XP;
    1081           0 :   T = FpX_get_red(T, p);
    1082           0 :   XP = FpX_Frobenius(T, p);
    1083           0 :   return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
    1084             : }
    1085             : 
    1086             : GEN
    1087           7 : FpX_ddf(GEN f, GEN p)
    1088             : {
    1089           7 :   pari_sp av = avma;
    1090             :   GEN F;
    1091           7 :   switch(ZX_factmod_init(&f, p))
    1092             :   {
    1093           7 :     case 0:  F = F2x_ddf(f);
    1094           7 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    1095           0 :     case 1:  F = Flx_ddf(f,p[2]);
    1096           0 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    1097           0 :     default: F = FpX_ddf_i(f,p); break;
    1098             :   }
    1099           7 :   return gerepilecopy(av, F);
    1100             : }
    1101             : 
    1102             : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
    1103             : static GEN
    1104          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1105             : {
    1106             :   GEN V, XP;
    1107             :   long i, l;
    1108          14 :   if (lgefint(p) == 3)
    1109             :   {
    1110           0 :     ulong pp = p[2];
    1111           0 :     return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
    1112             :   }
    1113          14 :   T = FpX_get_red(T, p);
    1114          14 :   XP = FpX_Frobenius(T, p);
    1115          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
    1116          14 :   for (i=1; i < l; i++) gel(V,i) = FpX_ddf_Shoup(gel(V,i), XP, p);
    1117          14 :   return vddf_to_simplefact(V, get_FpX_degree(T));
    1118             : }
    1119             : 
    1120             : static int
    1121           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1122             : {
    1123           0 :   pari_sp av = avma;
    1124             :   pari_timer ti;
    1125             :   long n, d;
    1126           0 :   GEN T = get_FpX_mod(Tp);
    1127           0 :   GEN dT = FpX_deriv(T, p);
    1128             :   GEN XP, D;
    1129           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    1130           0 :   n = get_FpX_degree(T);
    1131           0 :   T = FpX_get_red(Tp, p);
    1132           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1133           0 :   XP = FpX_Frobenius(T, p);
    1134           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1135           0 :   D = FpX_ddf_Shoup(T, XP, p);
    1136           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1137           0 :   d = degpol(gel(D, n));
    1138           0 :   avma = av; return d==n;
    1139             : }
    1140             : 
    1141             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1142             : 
    1143             : /*Assume that p is large and odd*/
    1144             : static GEN
    1145        1385 : FpX_factor_i(GEN f, GEN pp, long flag)
    1146             : {
    1147        1385 :   long d = degpol(f);
    1148        1385 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1149         630 :   switch(flag)
    1150             :   {
    1151         616 :     default: return FpX_factor_Cantor(f, pp);
    1152          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1153           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1154             :   }
    1155             : }
    1156             : 
    1157             : long
    1158           0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1159             : {
    1160           0 :   pari_sp av = avma;
    1161           0 :   long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
    1162           0 :   avma = av; return s;
    1163             : }
    1164             : 
    1165             : long
    1166           0 : FpX_nbfact(GEN T, GEN p)
    1167             : {
    1168           0 :   pari_sp av = avma;
    1169           0 :   GEN XP = FpX_Frobenius(T, p);
    1170           0 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1171           0 :   avma = av; return n;
    1172             : }
    1173             : 
    1174             : /* p > 2 */
    1175             : static GEN
    1176           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1177             : {
    1178           7 :   switch(d)
    1179             :   {
    1180             :     case -1:
    1181           0 :     case 0: return NULL;
    1182           0 :     case 1: return gen_1;
    1183             :   }
    1184           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1185             : }
    1186             : /* p > 2 */
    1187             : static GEN
    1188          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1189             : {
    1190          14 :   switch(d)
    1191             :   {
    1192           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1193           0 :     case 0: return trivial_fact();
    1194           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1195             :   }
    1196          14 :   switch(FpX_quad_factortype(f, p)) {
    1197           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1198           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1199           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1200             :   }
    1201             : }
    1202             : 
    1203             : GEN
    1204          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1205             : GEN
    1206       67720 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1207             : 
    1208             : /* not gerepile safe */
    1209             : static GEN
    1210         734 : FpX_factor_2(GEN f, GEN p, long d)
    1211             : {
    1212             :   GEN r, s, R, S;
    1213             :   long v;
    1214             :   int sgn;
    1215         734 :   switch(d)
    1216             :   {
    1217           7 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1218          30 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1219          45 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1220             :   }
    1221         652 :   r = FpX_quad_root(f, p, 1);
    1222         652 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1223         138 :   v = varn(f);
    1224         138 :   s = FpX_otherroot(f, r, p);
    1225         138 :   if (signe(r)) r = subii(p, r);
    1226         138 :   if (signe(s)) s = subii(p, s);
    1227         138 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1228         138 :   R = deg1pol_shallow(gen_1, r, v);
    1229         138 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1230          47 :   S = deg1pol_shallow(gen_1, s, v);
    1231          47 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1232             : }
    1233             : static GEN
    1234         755 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1235             : {
    1236         755 :   switch(flag) {
    1237           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1238          14 :     case 1: return FpX_degfact_2(f, p, d);
    1239         734 :     default: return FpX_factor_2(f, p, d);
    1240             :   }
    1241             : }
    1242             : 
    1243             : static int
    1244       58734 : F2x_quad_factortype(GEN x)
    1245       58734 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1246             : 
    1247             : static GEN
    1248           7 : F2x_is_irred_2(GEN f, long d)
    1249           7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1250             : 
    1251             : static GEN
    1252        7140 : F2x_degfact_2(GEN f, long d)
    1253             : {
    1254        7140 :   if (!d) return trivial_fact();
    1255        7140 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1256        6937 :   switch(F2x_quad_factortype(f)) {
    1257        2268 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1258        2219 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1259        2450 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1260             :   }
    1261             : }
    1262             : 
    1263             : static GEN
    1264      105616 : F2x_factor_2(GEN f, long d)
    1265             : {
    1266      105616 :   long v = f[1];
    1267      105616 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1268       97734 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1269       50442 :   switch(F2x_quad_factortype(f))
    1270             :   {
    1271       12208 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1272       23464 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1273       14770 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1274             :   }
    1275             : }
    1276             : static GEN
    1277      112763 : F2x_factor_deg2(GEN f, long d, long flag)
    1278             : {
    1279      112763 :   switch(flag) {
    1280           7 :     case 2: return F2x_is_irred_2(f, d);
    1281        7140 :     case 1: return F2x_degfact_2(f, d);
    1282      105616 :     default: return F2x_factor_2(f, d);
    1283             :   }
    1284             : }
    1285             : 
    1286             : /* xt = NULL or x^(p-1)/2 mod g */
    1287             : static void
    1288        4809 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
    1289             : {
    1290        4809 :   ulong q = p >> 1;
    1291        4809 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1292        4809 :   long d = degpol(a);
    1293        4809 :   if (d < 0)
    1294             :   {
    1295             :     ulong i;
    1296          42 :     split_add_done(S, (GEN)1);
    1297          42 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1298             :   } else {
    1299        4767 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1300        4767 :     if (d)
    1301             :     {
    1302        4767 :       if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
    1303        4767 :       a = Flx_gcd(a, xt, p);
    1304        4767 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1305             :     }
    1306             :   }
    1307        4809 : }
    1308             : static void
    1309        4809 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
    1310             : {
    1311        4809 :   ulong q = p >> 1;
    1312        4809 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1313        4809 :   long d = degpol(a);
    1314        4809 :   if (d < 0)
    1315             :   {
    1316          28 :     ulong i, z = Fl_nonsquare(p);
    1317          28 :     split_add_done(S, (GEN)z);
    1318          28 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1319             :   } else {
    1320        4781 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1321        4781 :     if (d)
    1322             :     {
    1323        4781 :       if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
    1324        4781 :       a = Flx_gcd(a, xt, p);
    1325        4781 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1326             :     }
    1327             :   }
    1328        4809 : }
    1329             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1330             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1331             : static int
    1332     4926904 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1333             : {
    1334     4926904 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1335     4926690 :   long d = degpol(g);
    1336     4926603 :   if (d < 0) return 0;
    1337     4926571 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1338     4926607 :   if (!d) return 1;
    1339     4910983 :   if ((p >> 4) <= (ulong)d)
    1340             :   { /* small p; split directly using x^((p-1)/2) +/- 1 */
    1341       14259 :     GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
    1342        9450 :                                 : NULL;
    1343        4809 :     split_squares(S, g, p, xt);
    1344        4809 :     split_nonsquares(S, g, p, xt);
    1345             :   } else { /* large p; use x^(p-1) - 1 directly */
    1346     4906174 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1347     4903474 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1348     4903474 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1349     4903747 :     g = Flx_gcd(g,a, p);
    1350     4904510 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1351             :   }
    1352     4912230 :   return 1;
    1353             : }
    1354             : 
    1355             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1356             : static GEN
    1357    21999947 : Flx_roots_i(GEN f, ulong p)
    1358             : {
    1359             :   GEN pol, g;
    1360    21999947 :   long v = Flx_valrem(f, &g);
    1361             :   ulong q;
    1362             :   struct split_t S;
    1363             : 
    1364             :   /* optimization: test for small degree first */
    1365    21999750 :   switch(degpol(g))
    1366             :   {
    1367             :     case 1: {
    1368       25515 :       ulong r = p - g[2];
    1369       25515 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1370             :     }
    1371             :     case 2: {
    1372    17040574 :       ulong r = Flx_quad_root(g, p, 1), s;
    1373    17067628 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1374    11655793 :       s = Flx_otherroot(g,r, p);
    1375    11703986 :       if (r < s)
    1376     2921934 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1377     8782052 :       else if (r > s)
    1378     8781891 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1379             :       else
    1380         161 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1381             :     }
    1382             :   }
    1383     4927525 :   q = p >> 1;
    1384     4927525 :   split_init(&S, lg(f)-1);
    1385     4926933 :   settyp(S.done, t_VECSMALL);
    1386     4926933 :   if (v) split_add_done(&S, (GEN)0);
    1387     4926933 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1388          42 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1389     4927694 :   pol = polx_Flx(f[1]);
    1390    10186195 :   for (pol[2]=1; ; pol[2]++)
    1391     5258047 :   {
    1392    10186195 :     long j, l = lg(S.todo);
    1393    10186195 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1394     5257915 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1395    10872105 :     for (j = 1; j < l; j++)
    1396             :     {
    1397     5614058 :       GEN b, c = gel(S.todo,j);
    1398             :       ulong r, s;
    1399     5614058 :       switch(degpol(c))
    1400             :       {
    1401             :         case 1:
    1402     4804548 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1403     4805068 :           j--; l--; break;
    1404             :         case 2:
    1405      371057 :           r = Flx_quad_root(c, p, 0);
    1406      371062 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1407      371055 :           s = Flx_otherroot(c,r, p);
    1408      371055 :           split_done(&S, j, (GEN)r, (GEN)s);
    1409      371057 :           j--; l--; break;
    1410             :         default:
    1411      438427 :           b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
    1412      438418 :           if (degpol(b) <= 0) continue;
    1413      337174 :           b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
    1414      337177 :           if (!degpol(b)) continue;
    1415      337043 :           b = Flx_normalize(b, p);
    1416      337045 :           c = Flx_div(c,b, p);
    1417      337040 :           split_todo(&S, j, b, c);
    1418             :       }
    1419             :     }
    1420             :   }
    1421             : }
    1422             : 
    1423             : GEN
    1424    21956158 : Flx_roots(GEN f, ulong p)
    1425             : {
    1426    21956158 :   pari_sp av = avma;
    1427    21956158 :   switch(lg(f))
    1428             :   {
    1429           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1430           0 :     case 3: avma = av; return cgetg(1, t_VECSMALL);
    1431             :   }
    1432    21957323 :   if (p == 2) return Flx_root_mod_2(f);
    1433    21957323 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1434             : }
    1435             : 
    1436             : /* assume x reduced mod p, monic. */
    1437             : static int
    1438      734573 : Flx_quad_factortype(GEN x, ulong p)
    1439             : {
    1440      734573 :   ulong b = x[3], c = x[2];
    1441      734573 :   return krouu(Fl_disc_bc(b, c, p), p);
    1442             : }
    1443             : static GEN
    1444          56 : Flx_is_irred_2(GEN f, ulong p, long d)
    1445             : {
    1446          56 :   if (!d) return NULL;
    1447          56 :   if (d == 1) return gen_1;
    1448          56 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1449             : }
    1450             : static GEN
    1451      751660 : Flx_degfact_2(GEN f, ulong p, long d)
    1452             : {
    1453      751660 :   if (!d) return trivial_fact();
    1454      751660 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1455      734517 :   switch(Flx_quad_factortype(f, p)) {
    1456      349076 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1457      376978 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1458        8463 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1459             :   }
    1460             : }
    1461             : /* p > 2 */
    1462             : static GEN
    1463      432348 : Flx_factor_2(GEN f, ulong p, long d)
    1464             : {
    1465             :   ulong r, s;
    1466             :   GEN R,S;
    1467      432348 :   long v = f[1];
    1468      432348 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1469      418175 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1470      318902 :   r = Flx_quad_root(f, p, 1);
    1471      318902 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1472      210978 :   s = Flx_otherroot(f, r, p);
    1473      210978 :   r = Fl_neg(r, p);
    1474      210978 :   s = Fl_neg(s, p);
    1475      210978 :   if (s < r) lswap(s,r);
    1476      210978 :   R = mkvecsmall3(v,r,1);
    1477      210978 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1478      183208 :   S = mkvecsmall3(v,s,1);
    1479      183208 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1480             : }
    1481             : static GEN
    1482     1184064 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1483             : {
    1484     1184064 :   switch(flag) {
    1485          56 :     case 2: return Flx_is_irred_2(f, p, d);
    1486      751660 :     case 1: return Flx_degfact_2(f, p, d);
    1487      432348 :     default: return Flx_factor_2(f, p, d);
    1488             :   }
    1489             : }
    1490             : 
    1491             : void
    1492       19957 : F2xV_to_FlxV_inplace(GEN v)
    1493             : {
    1494             :   long i;
    1495       19957 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1496       19957 : }
    1497             : void
    1498      730859 : FlxV_to_ZXV_inplace(GEN v)
    1499             : {
    1500             :   long i;
    1501      730859 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1502      730859 : }
    1503             : void
    1504      144681 : F2xV_to_ZXV_inplace(GEN v)
    1505             : {
    1506             :   long i;
    1507      144681 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1508      144681 : }
    1509             : 
    1510             : static GEN
    1511       10608 : F2x_Berlekamp_ker(GEN u)
    1512             : {
    1513       10608 :   pari_sp ltop=avma;
    1514       10608 :   long j,N = F2x_degree(u);
    1515             :   GEN Q;
    1516             :   pari_timer T;
    1517       10607 :   timer_start(&T);
    1518       10609 :   Q = F2x_matFrobenius(u);
    1519      249234 :   for (j=1; j<=N; j++)
    1520      238625 :     F2m_flip(Q,j,j);
    1521       10609 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
    1522       10609 :   Q = F2m_ker_sp(Q,0);
    1523       10609 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
    1524       10609 :   return gerepileupto(ltop,Q);
    1525             : }
    1526             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    1527             : static long
    1528       14567 : F2x_split_Berlekamp(GEN *t)
    1529             : {
    1530       14567 :   GEN u = *t, a, b, vker;
    1531       14567 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    1532             : 
    1533       14570 :   if (du == 1) return 1;
    1534       11212 :   if (du == 2)
    1535             :   {
    1536         604 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    1537             :     {
    1538           0 :       t[0] = mkvecsmall2(sv, 2);
    1539           0 :       t[1] = mkvecsmall2(sv, 3);
    1540           0 :       return 2;
    1541             :     }
    1542         604 :     return 1;
    1543             :   }
    1544             : 
    1545       10608 :   vker = F2x_Berlekamp_ker(u);
    1546       10609 :   lb = lgcols(vker);
    1547       10628 :   d = lg(vker)-1;
    1548       10628 :   ir = 0;
    1549             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    1550       52986 :   for (L=1; L<d; )
    1551             :   {
    1552             :     GEN pol;
    1553       31751 :     if (d == 2)
    1554        1744 :       pol = F2v_to_F2x(gel(vker,2), sv);
    1555             :     else
    1556             :     {
    1557       30007 :       GEN v = zero_zv(lb);
    1558       30014 :       v[1] = du;
    1559       30014 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    1560      112207 :       for (i=2; i<=d; i++)
    1561       82192 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    1562       30015 :       pol = F2v_to_F2x(v, sv);
    1563             :     }
    1564       98654 :     for (i=ir; i<L && L<d; i++)
    1565             :     {
    1566       66924 :       a = t[i]; la = F2x_degree(a);
    1567       66897 :       if (la == 1) { set_irred(i); }
    1568       66735 :       else if (la == 2)
    1569             :       {
    1570         744 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    1571             :         {
    1572           0 :           t[i] = mkvecsmall2(sv, 2);
    1573           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    1574             :         }
    1575         744 :         set_irred(i);
    1576             :       }
    1577             :       else
    1578             :       {
    1579       65991 :         pari_sp av = avma;
    1580             :         long lb;
    1581       65991 :         b = F2x_rem(pol, a);
    1582       65977 :         if (F2x_degree(b) <= 0) { avma=av; continue; }
    1583       21694 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    1584       21701 :         if (lb && lb < la)
    1585             :         {
    1586       21701 :           t[L] = F2x_div(a,b);
    1587       21694 :           t[i]= b; L++;
    1588             :         }
    1589           0 :         else avma = av;
    1590             :       }
    1591             :     }
    1592             :   }
    1593       10607 :   return d;
    1594             : }
    1595             : /* assume deg f > 2 */
    1596             : static GEN
    1597       11699 : F2x_Berlekamp_i(GEN f, long flag)
    1598             : {
    1599       11699 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    1600             :   GEN y, E, t, V;
    1601             : 
    1602       11700 :   val = F2x_valrem(f, &f);
    1603       11700 :   if (flag == 2 && val) return NULL;
    1604       11693 :   V = F2x_factor_squarefree(f); lV = lg(V);
    1605       11693 :   if (flag == 2 && lV > 2) return NULL;
    1606             : 
    1607             :   /* to hold factors and exponents */
    1608       11623 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    1609       11623 :   E = cgetg(d+1,t_VECSMALL);
    1610       11624 :   lfact = 1;
    1611       11624 :   if (val) {
    1612        3482 :     if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
    1613        3482 :     E[1] = val; lfact++;
    1614             :   }
    1615             : 
    1616       49054 :   for (k=1; k<lV; k++)
    1617             :   {
    1618       37586 :     if (F2x_degree(gel(V, k))==0) continue;
    1619       14568 :     gel(t,lfact) = gel(V, k);
    1620       14568 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    1621       14568 :     if (flag == 2 && d != 1) return NULL;
    1622       14414 :     if (flag == 1)
    1623        1533 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    1624       14414 :     for (j=0; j<d; j++) E[lfact+j] = k;
    1625       14414 :     lfact += d;
    1626             :   }
    1627       11468 :   if (flag == 2) return gen_1; /* irreducible */
    1628       11454 :   setlg(t, lfact);
    1629       11454 :   setlg(E, lfact); y = mkvec2(t,E);
    1630             :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    1631       11454 :               : sort_factor_pol(y, cmpGuGu);
    1632             : }
    1633             : 
    1634             : /* Adapted from Shoup NTL */
    1635             : GEN
    1636       76228 : F2x_factor_squarefree(GEN f)
    1637             : {
    1638             :   GEN r, t, v, tv;
    1639       76228 :   long i, q, n = F2x_degree(f);
    1640       76229 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1641      135235 :   for(q = 1;;q *= 2)
    1642             :   {
    1643      194241 :     r = F2x_gcd(f, F2x_deriv(f));
    1644      135228 :     if (F2x_degree(r) == 0)
    1645             :     {
    1646       62418 :       gel(u, q) = f;
    1647       62418 :       break;
    1648             :     }
    1649       72814 :     t = F2x_div(f, r);
    1650       72814 :     if (F2x_degree(t) > 0)
    1651             :     {
    1652             :       long j;
    1653       69267 :       for(j = 1;;j++)
    1654             :       {
    1655      107319 :         v = F2x_gcd(r, t);
    1656       69266 :         tv = F2x_div(t, v);
    1657       69266 :         if (F2x_degree(tv) > 0)
    1658       32511 :           gel(u, j*q) = tv;
    1659       69266 :         if (F2x_degree(v) <= 0) break;
    1660       38052 :         r = F2x_div(r, v);
    1661       38052 :         t = v;
    1662             :       }
    1663       31214 :       if (F2x_degree(r) == 0) break;
    1664             :     }
    1665       59005 :     f = F2x_sqrt(r);
    1666             :   }
    1667      560996 :   for (i = n; i; i--)
    1668      560683 :     if (F2x_degree(gel(u,i))) break;
    1669       76277 :   setlg(u,i+1); return u;
    1670             : }
    1671             : 
    1672             : static GEN
    1673       79965 : F2x_ddf_simple(GEN T, GEN XP)
    1674             : {
    1675       79965 :   pari_sp av = avma, av2;
    1676             :   GEN f, z, Tr, X;
    1677       79965 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1678       79967 :   if (n == 0) return cgetg(1, t_VEC);
    1679       79967 :   if (n == 1) return mkvec(T);
    1680       37482 :   z = XP; Tr = T; X = polx_F2x(v);
    1681       37484 :   f = const_vec(n, pol1_F2x(v));
    1682       37484 :   av2 = avma;
    1683      133515 :   for (j = 1; j <= B; j++)
    1684             :   {
    1685      100568 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1686      100589 :     if (F2x_degree(u))
    1687             :     {
    1688       22044 :       gel(f, j) = u;
    1689       22044 :       Tr = F2x_div(Tr, u);
    1690       22049 :       av2 = avma;
    1691       78552 :     } else z = gerepileuptoleaf(av2, z);
    1692      100615 :     if (!F2x_degree(Tr)) break;
    1693       96067 :     z = F2xq_sqr(z, Tr);
    1694             :   }
    1695       37475 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1696       37484 :   return gerepilecopy(av, f);
    1697             : }
    1698             : 
    1699             : GEN
    1700           7 : F2x_ddf(GEN T)
    1701             : {
    1702             :   GEN XP;
    1703           7 :   T = F2x_get_red(T);
    1704           7 :   XP = F2x_Frobenius(T);
    1705           7 :   return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
    1706             : }
    1707             : 
    1708             : static GEN
    1709        6574 : F2xq_frobtrace(GEN a, long d, GEN T)
    1710             : {
    1711        6574 :   pari_sp av = avma;
    1712             :   long i;
    1713        6574 :   GEN x = a;
    1714       27106 :   for(i=1; i<d; i++)
    1715             :   {
    1716       20530 :     x = F2x_add(a, F2xq_sqr(x,T));
    1717       20532 :     if (gc_needed(av, 2))
    1718           0 :       x = gerepileuptoleaf(av, x);
    1719             :   }
    1720        6576 :   return x;
    1721             : }
    1722             : 
    1723             : static void
    1724        9782 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1725             : {
    1726        9782 :   long n = F2x_degree(Tp), r = n/d;
    1727             :   GEN T, f, ff;
    1728        9782 :   if (r==1) { gel(V, idx) = Tp; return; }
    1729        3333 :   T = Tp;
    1730        3333 :   XP = F2x_rem(XP, T);
    1731             :   while (1)
    1732        3243 :   {
    1733        6576 :     pari_sp btop = avma;
    1734             :     long df;
    1735        6576 :     GEN g = random_F2x(n, Tp[1]);
    1736        6576 :     GEN t = F2xq_frobtrace(g, d, T);
    1737        6576 :     if (lgpol(t) == 0) continue;
    1738        4883 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1739        4883 :     if (df > 0 && df < n) break;
    1740        1550 :     avma = btop;
    1741             :   }
    1742        3333 :   ff = F2x_div(Tp, f);
    1743        3333 :   F2x_edf_simple(f, XP, d, V, idx);
    1744        3333 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1745             : }
    1746             : 
    1747             : static GEN
    1748       79963 : F2x_factor_Shoup(GEN T)
    1749             : {
    1750       79963 :   long i, n, s = 0;
    1751             :   GEN XP, D, V;
    1752             :   pari_timer ti;
    1753       79963 :   n = F2x_degree(T);
    1754       79962 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1755       79962 :   XP = F2x_Frobenius(T);
    1756       79958 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1757       79958 :   D = F2x_ddf_simple(T, XP);
    1758       79963 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1759      348114 :   for (i = 1; i <= n; i++)
    1760      268152 :     s += F2x_degree(gel(D,i))/i;
    1761       79962 :   V = cgetg(s+1, t_COL);
    1762      348083 :   for (i = 1, s = 1; i <= n; i++)
    1763             :   {
    1764      268122 :     GEN Di = gel(D,i);
    1765      268122 :     long ni = F2x_degree(Di), ri = ni/i;
    1766      268115 :     if (ni == 0) continue;
    1767       97473 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1768        3116 :     F2x_edf_simple(Di, XP, i, V, s);
    1769        3116 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1770        3116 :     s += ri;
    1771             :   }
    1772       79961 :   return V;
    1773             : }
    1774             : 
    1775             : static GEN
    1776       64535 : F2x_factor_Cantor(GEN T)
    1777             : {
    1778       64535 :   GEN E, F, V = F2x_factor_squarefree(T);
    1779       64535 :   long i, j, l = lg(V);
    1780       64535 :   E = cgetg(l, t_VEC);
    1781       64536 :   F = cgetg(l, t_VEC);
    1782      256622 :   for (i=1, j=1; i < l; i++)
    1783      192086 :     if (F2x_degree(gel(V,i)))
    1784             :     {
    1785       79963 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1786       79961 :       gel(F, j) = Fj;
    1787       79961 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1788       79963 :       j++;
    1789             :     }
    1790       64536 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1791             : }
    1792             : 
    1793             : #if 0
    1794             : static GEN
    1795             : F2x_simplefact_Shoup(GEN T)
    1796             : {
    1797             :   long i, n, s = 0, j = 1, k;
    1798             :   GEN XP, D, V;
    1799             :   pari_timer ti;
    1800             :   n = F2x_degree(T);
    1801             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1802             :   XP = F2x_Frobenius(T);
    1803             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1804             :   D = F2x_ddf_simple(T, XP);
    1805             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1806             :   for (i = 1; i <= n; i++)
    1807             :     s += F2x_degree(gel(D,i))/i;
    1808             :   V = cgetg(s+1, t_VECSMALL);
    1809             :   for (i = 1; i <= n; i++)
    1810             :   {
    1811             :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1812             :     if (ni == 0) continue;
    1813             :     for (k = 1; k <= ri; k++)
    1814             :       V[j++] = i;
    1815             :   }
    1816             :   return V;
    1817             : }
    1818             : static GEN
    1819             : F2x_simplefact_Cantor(GEN T)
    1820             : {
    1821             :   GEN E, F, V = F2x_factor_squarefree(T);
    1822             :   long i, j, l = lg(V);
    1823             :   F = cgetg(l, t_VEC);
    1824             :   E = cgetg(l, t_VEC);
    1825             :   for (i=1, j=1; i < l; i++)
    1826             :     if (F2x_degree(gel(V,i)))
    1827             :     {
    1828             :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1829             :       gel(F, j) = Fj;
    1830             :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1831             :       j++;
    1832             :     }
    1833             :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1834             : }
    1835             : static int
    1836             : F2x_isirred_Cantor(GEN T)
    1837             : {
    1838             :   pari_sp av = avma;
    1839             :   pari_timer ti;
    1840             :   long n, d;
    1841             :   GEN dT = F2x_deriv(T);
    1842             :   GEN XP, D;
    1843             :   if (F2x_degree(F2x_gcd(T, dT)) != 0) { avma = av; return 0; }
    1844             :   n = F2x_degree(T);
    1845             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1846             :   XP = F2x_Frobenius(T);
    1847             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1848             :   D = F2x_ddf_simple(T, XP);
    1849             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1850             :   d = F2x_degree(gel(D, n));
    1851             :   avma = av; return d==n;
    1852             : }
    1853             : #endif
    1854             : 
    1855             : /* driver for Cantor factorization, assume deg f > 2; not competitive for
    1856             :  * flag != 0, or as deg f increases */
    1857             : static GEN
    1858       64536 : F2x_Cantor_i(GEN f, long flag)
    1859             : {
    1860             :   switch(flag)
    1861             :   {
    1862       64536 :     default: return F2x_factor_Cantor(f);
    1863             : #if 0
    1864             :     case 1: return F2x_simplefact_Cantor(f);
    1865             :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1866             : #endif
    1867             :   }
    1868             : }
    1869             : static GEN
    1870      188999 : F2x_factor_i(GEN f, long flag)
    1871             : {
    1872      188999 :   long d = F2x_degree(f);
    1873      188998 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1874       74219 :   return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
    1875      140771 :                                : F2x_Berlekamp_i(f, flag);
    1876             : }
    1877             : 
    1878             : GEN
    1879           0 : F2x_degfact(GEN f)
    1880             : {
    1881           0 :   pari_sp av = avma;
    1882           0 :   GEN z = F2x_factor_i(f, 1);
    1883           0 :   return gerepilecopy(av, z);
    1884             : }
    1885             : 
    1886             : int
    1887         238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
    1888             : 
    1889             : /* Adapted from Shoup NTL */
    1890             : GEN
    1891      772475 : Flx_factor_squarefree(GEN f, ulong p)
    1892             : {
    1893      772475 :   long i, q, n = degpol(f);
    1894      772476 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1895      826792 :   for(q = 1;;q *= p)
    1896       54316 :   {
    1897      826792 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1898      826793 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1899       99718 :     t = Flx_div(f, r, p);
    1900       99718 :     if (degpol(t) > 0)
    1901             :     {
    1902             :       long j;
    1903      146861 :       for(j = 1;;j++)
    1904             :       {
    1905      247535 :         v = Flx_gcd(r, t, p);
    1906      146861 :         tv = Flx_div(t, v, p);
    1907      146860 :         if (degpol(tv) > 0)
    1908       70504 :           gel(u, j*q) = Flx_normalize(tv, p);
    1909      146860 :         if (degpol(v) <= 0) break;
    1910      100674 :         r = Flx_div(r, v, p);
    1911      100674 :         t = v;
    1912             :       }
    1913       46187 :       if (degpol(r) == 0) break;
    1914             :     }
    1915       54316 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1916             :   }
    1917     4216880 :   for (i = n; i; i--)
    1918     4216880 :     if (degpol(gel(u,i))) break;
    1919      772477 :   setlg(u,i+1); return u;
    1920             : }
    1921             : 
    1922             : long
    1923         126 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1924             : {
    1925         126 :   pari_sp av = avma;
    1926             :   ulong lc;
    1927             :   GEN F;
    1928         126 :   long i, n = degpol(f), v = f[1], l;
    1929         126 :   if (n % k) return 0;
    1930         126 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    1931         126 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    1932         126 :   F = Flx_factor_squarefree(f, p); l = lg(F)-1;
    1933        6405 :   for (i = 1; i <= l; i++)
    1934        6300 :     if (i%k && degpol(gel(F,i))) { avma = av; return 0; }
    1935         105 :   if (pt_r)
    1936             :   {
    1937         105 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    1938        6384 :     for(i = l; i >= 1; i--)
    1939             :     {
    1940        6279 :       if (i%k) continue;
    1941        1274 :       s = Flx_mul(s, gel(F,i), p);
    1942        1274 :       r = Flx_mul(r, s, p);
    1943             :     }
    1944         105 :     *pt_r = gerepileuptoleaf(av, r);
    1945           0 :   } else av = avma;
    1946         105 :   return 1;
    1947             : }
    1948             : 
    1949             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1950             : static GEN
    1951     1085569 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
    1952             : {
    1953     1085569 :   pari_sp av = avma;
    1954             :   GEN b, g, h, F, f, Tr, xq;
    1955             :   long i, j, n, v, bo, ro;
    1956             :   long B, l, m;
    1957             :   pari_timer ti;
    1958     1085569 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1959     1085569 :   if (n == 0) return cgetg(1, t_VEC);
    1960     1083554 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1961      986536 :   B = n/2;
    1962      986536 :   l = usqrt(B);
    1963      986536 :   m = (B+l-1)/l;
    1964      986536 :   T = Flx_get_red(T, p);
    1965      986535 :   b = cgetg(l+2, t_VEC);
    1966      986535 :   gel(b, 1) = polx_Flx(v);
    1967      986535 :   gel(b, 2) = XP;
    1968      986535 :   bo = brent_kung_optpow(n, l-1, 1);
    1969      986535 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    1970      986535 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1971      986535 :   if (expu(p) <= ro)
    1972      237627 :     for (i = 3; i <= l+1; i++)
    1973      133915 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1974             :   else
    1975             :   {
    1976      882823 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1977      882824 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
    1978     1027391 :     for (i = 3; i <= l+1; i++)
    1979      144567 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1980             :   }
    1981      986536 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
    1982      986536 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1983      986536 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
    1984      986536 :   g = cgetg(m+1, t_VEC);
    1985      986535 :   gel(g, 1) = gel(xq, 2);
    1986     1847628 :   for(i = 2; i <= m; i++)
    1987      861093 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1988      986535 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
    1989      986535 :   h = cgetg(m+1, t_VEC);
    1990     2834165 :   for (j = 1; j <= m; j++)
    1991             :   {
    1992     1847629 :     pari_sp av = avma;
    1993     1847629 :     GEN gj = gel(g, j);
    1994     1847629 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1995     2883889 :     for (i = 2; i <= l; i++)
    1996     1036269 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1997     1847620 :     gel(h, j) = gerepileupto(av, e);
    1998             :   }
    1999      986536 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
    2000      986536 :   Tr = get_Flx_mod(T);
    2001      986536 :   F = cgetg(m+1, t_VEC);
    2002     2834160 :   for (j = 1; j <= m; j++)
    2003             :   {
    2004     1847625 :     GEN u = Flx_gcd(Tr, gel(h, j), p);
    2005     1847626 :     if (degpol(u))
    2006             :     {
    2007      762574 :       u = Flx_normalize(u, p);
    2008      762574 :       Tr = Flx_div(Tr, u, p);
    2009             :     }
    2010     1847623 :     gel(F, j) = u;
    2011             :   }
    2012      986535 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
    2013      986535 :   f = const_vec(n, pol1_Flx(v));
    2014     2834159 :   for (j = 1; j <= m; j++)
    2015             :   {
    2016     1847623 :     GEN e = gel(F, j);
    2017     2038234 :     for (i=l-1; i >= 0; i--)
    2018             :     {
    2019     2038230 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    2020     2038233 :       if (degpol(u))
    2021             :       {
    2022      783756 :         gel(f, l*j-i) = u;
    2023      783756 :         e = Flx_div(e, u, p);
    2024             :       }
    2025     2038227 :       if (!degpol(e)) break;
    2026             :     }
    2027             :   }
    2028      986536 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
    2029      986536 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    2030      986536 :   return gerepilecopy(av, f);
    2031             : }
    2032             : 
    2033             : static void
    2034      128014 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2035             : {
    2036      128014 :   long n = degpol(Tp), r = n/d;
    2037             :   GEN T, f, ff;
    2038             :   ulong p2;
    2039      128014 :   if (r==1) { gel(V, idx) = Tp; return; }
    2040       55422 :   p2 = p>>1;
    2041       55422 :   T = Flx_get_red(Tp, p);
    2042       55422 :   XP = Flx_rem(XP, T, p);
    2043             :   while (1)
    2044        6253 :   {
    2045       61675 :     pari_sp btop = avma;
    2046             :     long i;
    2047       61675 :     GEN g = random_Flx(n, Tp[1], p);
    2048       61675 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2049       61674 :     if (lgpol(t) == 0) continue;
    2050      134633 :     for(i=1; i<=10; i++)
    2051             :     {
    2052      129968 :       pari_sp btop2 = avma;
    2053      129968 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    2054      129967 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    2055      129963 :       if (degpol(f) > 0 && degpol(f) < n) break;
    2056       74547 :       avma = btop2;
    2057             :     }
    2058       60087 :     if (degpol(f) > 0 && degpol(f) < n) break;
    2059        4665 :     avma = btop;
    2060             :   }
    2061       55422 :   f = Flx_normalize(f, p);
    2062       55422 :   ff = Flx_div(Tp, f ,p);
    2063       55422 :   Flx_edf_simple(f, XP, d, p, V, idx);
    2064       55422 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    2065             : }
    2066             : static void
    2067             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    2068             : 
    2069             : static void
    2070      348902 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    2071             : {
    2072             :   pari_sp av;
    2073      348902 :   GEN Tp = get_Flx_mod(T);
    2074      348902 :   long n = degpol(hp), vT = Tp[1];
    2075             :   GEN u1, u2, f1, f2;
    2076      348902 :   ulong p2 = p>>1;
    2077             :   GEN R, h;
    2078      348902 :   h = Flx_get_red(hp, p);
    2079      348902 :   t = Flx_rem(t, T, p);
    2080      348902 :   av = avma;
    2081             :   do
    2082             :   {
    2083      572502 :     avma = av;
    2084      572502 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    2085      572502 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    2086      572502 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2087      348902 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    2088      348902 :   f1 = Flx_normalize(f1, p);
    2089      348902 :   u2 = Flx_div(hp, u1, p);
    2090      348902 :   f2 = Flx_div(Tp, f1, p);
    2091      348902 :   if (degpol(u1)==1)
    2092             :   {
    2093      252665 :     if (degpol(f1)==d)
    2094      248470 :       gel(V, idx) = f1;
    2095             :     else
    2096        4195 :       Flx_edf(f1, XP, d, p, V, idx);
    2097             :   }
    2098             :   else
    2099       96237 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    2100      348902 :   idx += degpol(f1)/d;
    2101      348902 :   if (degpol(u2)==1)
    2102             :   {
    2103      250323 :     if (degpol(f2)==d)
    2104      246468 :       gel(V, idx) = f2;
    2105             :     else
    2106        3855 :       Flx_edf(f2, XP, d, p, V, idx);
    2107             :   }
    2108             :   else
    2109       98579 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    2110      348902 : }
    2111             : 
    2112             : static void
    2113      154086 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2114             : {
    2115      154086 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2116             :   GEN T, h, t;
    2117             :   pari_timer ti;
    2118      154086 :   if (r==1) { gel(V, idx) = Tp; return; }
    2119      154086 :   T = Flx_get_red(Tp, p);
    2120      154086 :   XP = Flx_rem(XP, T, p);
    2121      154086 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2122             :   do
    2123             :   {
    2124      158993 :     GEN g = random_Flx(n, vT, p);
    2125      158993 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2126      158993 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2127      158993 :     h = Flxq_minpoly(t, T, p);
    2128      158993 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2129      158993 :   } while (degpol(h) <= 1);
    2130      154086 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    2131             : }
    2132             : 
    2133             : static GEN
    2134      440537 : Flx_factor_Shoup(GEN T, ulong p)
    2135             : {
    2136      440537 :   long i, n, s = 0;
    2137             :   GEN XP, D, V;
    2138      440537 :   long e = expu(p);
    2139             :   pari_timer ti;
    2140      440537 :   n = get_Flx_degree(T);
    2141      440537 :   T = Flx_get_red(T, p);
    2142      440537 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2143      440537 :   XP = Flx_Frobenius(T, p);
    2144      440538 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2145      440538 :   D = Flx_ddf_Shoup(T, XP, p);
    2146      440538 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2147      440538 :   s = ddf_to_nbfact(D);
    2148      440538 :   V = cgetg(s+1, t_COL);
    2149     2481985 :   for (i = 1, s = 1; i <= n; i++)
    2150             :   {
    2151     2041447 :     GEN Di = gel(D,i);
    2152     2041447 :     long ni = degpol(Di), ri = ni/i;
    2153     2041447 :     if (ni == 0) continue;
    2154      555913 :     Di = Flx_normalize(Di, p);
    2155      555913 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2156      163206 :     if (ri <= e*expu(e))
    2157      146036 :       Flx_edf(Di, XP, i, p, V, s);
    2158             :     else
    2159       17170 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2160      163206 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2161      163206 :     s += ri;
    2162             :   }
    2163      440538 :   return V;
    2164             : }
    2165             : 
    2166             : static GEN
    2167      416934 : Flx_factor_Cantor(GEN T, ulong p)
    2168             : {
    2169      416934 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2170      416935 :   long i, j, l = lg(V);
    2171      416935 :   F = cgetg(l, t_VEC);
    2172      416935 :   E = cgetg(l, t_VEC);
    2173     1051205 :   for (i=1, j=1; i < l; i++)
    2174      634270 :     if (degpol(gel(V,i)))
    2175             :     {
    2176      440537 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2177      440538 :       gel(F, j) = Fj;
    2178      440538 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2179      440538 :       j++;
    2180             :     }
    2181      416935 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2182             : }
    2183             : 
    2184             : GEN
    2185           0 : Flx_ddf(GEN T, ulong p)
    2186             : {
    2187             :   GEN XP;
    2188           0 :   T = Flx_get_red(T, p);
    2189           0 :   XP = Flx_Frobenius(T, p);
    2190           0 :   return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
    2191             : }
    2192             : 
    2193             : static GEN
    2194      354639 : Flx_simplefact_Cantor(GEN T, ulong p)
    2195             : {
    2196             :   GEN XP, V;
    2197             :   long i, l;
    2198      354639 :   T = Flx_get_red(T, p);
    2199      354639 :   XP = Flx_Frobenius(T, p);
    2200      354639 :   V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
    2201      354639 :   for (i=1; i < l; i++) gel(V,i) = Flx_ddf_Shoup(gel(V,i), XP, p);
    2202      354639 :   return vddf_to_simplefact(V, get_Flx_degree(T));
    2203             : }
    2204             : 
    2205             : static int
    2206        1057 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2207             : {
    2208        1057 :   pari_sp av = avma;
    2209             :   pari_timer ti;
    2210             :   long n, d;
    2211        1057 :   GEN T = get_Flx_mod(Tp);
    2212        1057 :   GEN dT = Flx_deriv(T, p);
    2213             :   GEN XP, D;
    2214        1057 :   if (degpol(Flx_gcd(T, dT, p)) != 0) { avma = av; return 0; }
    2215         770 :   n = get_Flx_degree(T);
    2216         770 :   T = Flx_get_red(Tp, p);
    2217         770 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2218         770 :   XP = Flx_Frobenius(T, p);
    2219         770 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2220         770 :   D = Flx_ddf_Shoup(T, XP, p);
    2221         770 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2222         770 :   d = degpol(gel(D, n));
    2223         770 :   avma = av; return d==n;
    2224             : }
    2225             : 
    2226             : /* f monic */
    2227             : static GEN
    2228     1985555 : Flx_factor_i(GEN f, ulong pp, long flag)
    2229             : {
    2230             :   long d;
    2231     1985555 :   if (pp==2) { /*We need to handle 2 specially */
    2232       28861 :     GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
    2233       28861 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2234       28861 :     return F;
    2235             :   }
    2236     1956694 :   d = degpol(f);
    2237     1956694 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2238      772630 :   switch(flag)
    2239             :   {
    2240      416934 :     default: return Flx_factor_Cantor(f, pp);
    2241      354639 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2242        1057 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2243             :   }
    2244             : }
    2245             : 
    2246             : GEN
    2247     1115196 : Flx_degfact(GEN f, ulong p)
    2248             : {
    2249     1115196 :   pari_sp av = avma;
    2250     1115196 :   GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
    2251     1115196 :   return gerepilecopy(av, z);
    2252             : }
    2253             : 
    2254             : /* T must be squarefree mod p*/
    2255             : GEN
    2256      191992 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2257             : {
    2258             :   GEN XP, D;
    2259             :   pari_timer ti;
    2260      191992 :   long i, s, n = get_Flx_degree(T);
    2261      191992 :   GEN V = const_vecsmall(n, 0);
    2262      191992 :   pari_sp av = avma;
    2263      191992 :   T = Flx_get_red(T, p);
    2264      191992 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2265      191992 :   XP = Flx_Frobenius(T, p);
    2266      191992 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2267      191992 :   D = Flx_ddf_Shoup(T, XP, p);
    2268      191992 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2269     1228714 :   for (i = 1, s = 0; i <= n; i++)
    2270             :   {
    2271     1036722 :     V[i] = degpol(gel(D,i))/i;
    2272     1036722 :     s += V[i];
    2273             :   }
    2274      191992 :   *nb = s;
    2275      191992 :   avma = av; return V;
    2276             : }
    2277             : 
    2278             : long
    2279       94626 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2280             : {
    2281       94626 :   pari_sp av = avma;
    2282       94626 :   long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
    2283       94626 :   avma = av; return s;
    2284             : }
    2285             : 
    2286             : /* T must be squarefree mod p*/
    2287             : long
    2288       94626 : Flx_nbfact(GEN T, ulong p)
    2289             : {
    2290       94626 :   pari_sp av = avma;
    2291       94626 :   GEN XP = Flx_Frobenius(T, p);
    2292       94626 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2293       94626 :   avma = av; return n;
    2294             : }
    2295             : 
    2296             : int
    2297        1057 : Flx_is_irred(GEN f, ulong p)
    2298             : {
    2299        1057 :   pari_sp av = avma;
    2300        1057 :   int z = !!Flx_factor_i(Flx_normalize(f,p),p,2);
    2301        1057 :   avma = av; return z;
    2302             : }
    2303             : 
    2304             : /* Use this function when you think f is reducible, and that there are lots of
    2305             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2306             : int
    2307          77 : FpX_is_irred(GEN f, GEN p)
    2308             : {
    2309          77 :   pari_sp av = avma;
    2310             :   int z;
    2311          77 :   switch(ZX_factmod_init(&f,p))
    2312             :   {
    2313          14 :     case 0:  z = !!F2x_factor_i(f,2); break;
    2314          56 :     case 1:  z = !!Flx_factor_i(f,p[2],2); break;
    2315           7 :     default: z = !!FpX_factor_i(f,p,2); break;
    2316             :   }
    2317          77 :   avma = av; return z;
    2318             : }
    2319             : GEN
    2320          42 : FpX_degfact(GEN f, GEN p) {
    2321          42 :   pari_sp av = avma;
    2322             :   GEN F;
    2323          42 :   switch(ZX_factmod_init(&f,p))
    2324             :   {
    2325           7 :     case 0:  F = F2x_factor_i(f,1); break;
    2326           7 :     case 1:  F = Flx_factor_i(f,p[2],1); break;
    2327          28 :     default: F = FpX_factor_i(f,p,1); break;
    2328             :   }
    2329          42 :   return gerepilecopy(av, F);
    2330             : }
    2331             : 
    2332             : #if 0
    2333             : /* set x <-- x + c*y mod p */
    2334             : /* x is not required to be normalized.*/
    2335             : static void
    2336             : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2337             : {
    2338             :   long i, lx, ly;
    2339             :   ulong *x=(ulong *)gx;
    2340             :   ulong *y=(ulong *)gy;
    2341             :   if (!c) return;
    2342             :   lx = lg(gx);
    2343             :   ly = lg(gy);
    2344             :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2345             :   if (SMALL_ULONG(p))
    2346             :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2347             :   else
    2348             :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2349             : }
    2350             : #endif
    2351             : 
    2352             : GEN
    2353      874481 : FpX_factor(GEN f, GEN p)
    2354             : {
    2355      874481 :   pari_sp av = avma;
    2356             :   GEN F;
    2357      874481 :   switch(ZX_factmod_init(&f, p))
    2358             :   {
    2359      144674 :     case 0:  F = F2x_factor_i(f,0);
    2360      144674 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    2361      728457 :     case 1:  F = Flx_factor_i(f,p[2],0);
    2362      728457 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    2363        1350 :     default: F = FpX_factor_i(f,p,0); break;
    2364             :   }
    2365      874481 :   return gerepilecopy(av, F);
    2366             : }
    2367             : 
    2368             : GEN
    2369      140781 : Flx_factor(GEN f, ulong p)
    2370             : {
    2371      140781 :   pari_sp av = avma;
    2372      140781 :   return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
    2373             : }
    2374             : GEN
    2375       15205 : F2x_factor(GEN f)
    2376             : {
    2377       15205 :   pari_sp av = avma;
    2378       15205 :   return gerepilecopy(av, F2x_factor_i(f,0));
    2379             : }

Generated by: LCOV version 1.13