Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 1294 1399 92.5 %
Date: 2022-07-03 07:33:15 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_factormod
      19             : 
      20             : /***********************************************************************/
      21             : /**                                                                   **/
      22             : /**               Factorisation over finite field                     **/
      23             : /**                                                                   **/
      24             : /***********************************************************************/
      25             : 
      26             : /*******************************************************************/
      27             : /*                                                                 */
      28             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      29             : /*                                                                 */
      30             : /*******************************************************************/
      31             : /* Replace F by a monic normalized FpX having the same factors;
      32             :  * assume p prime and *F a ZX */
      33             : static int
      34     3459405 : ZX_factmod_init(GEN *F, GEN p)
      35             : {
      36     3459405 :   if (lgefint(p) == 3)
      37             :   {
      38     3457715 :     ulong pp = p[2];
      39     3457715 :     if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
      40     2079705 :     *F = ZX_to_Flx(*F, pp);
      41     2079795 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42     2079786 :     return 1;
      43             :   }
      44        1690 :   *F = FpX_red(*F, p);
      45        1695 :   if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      46        1695 :   return 2;
      47             : }
      48             : static void
      49      367625 : ZX_rootmod_init(GEN *F, GEN p)
      50             : {
      51      367625 :   if (lgefint(p) == 3)
      52             :   {
      53      359229 :     ulong pp = p[2];
      54      359229 :     *F = ZX_to_Flx(*F, pp);
      55      359225 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      56             :   }
      57             :   else
      58             :   {
      59        8396 :     *F = FpX_red(*F, p);
      60        8395 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      61             :   }
      62      367621 : }
      63             : 
      64             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      65             : static GEN
      66         434 : all_roots_mod_p(ulong p, int not_0)
      67             : {
      68             :   GEN r;
      69             :   ulong i;
      70         434 :   if (not_0) {
      71         329 :     r = cgetg(p, t_VECSMALL);
      72         987 :     for (i = 1; i < p; i++) r[i] = i;
      73             :   } else {
      74         105 :     r = cgetg(p+1, t_VECSMALL);
      75         448 :     for (i = 0; i < p; i++) r[i+1] = i;
      76             :   }
      77         434 :   return r;
      78             : }
      79             : 
      80             : /* X^n - 1 */
      81             : static GEN
      82        3188 : Flx_Xnm1(long sv, long n, ulong p)
      83             : {
      84        3188 :   GEN t = cgetg(n+3, t_VECSMALL);
      85             :   long i;
      86        3188 :   t[1] = sv;
      87        3188 :   t[2] = p - 1;
      88        9824 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      89        3188 :   t[i] = 1; return t;
      90             : }
      91             : /* X^n + 1 */
      92             : static GEN
      93        3135 : Flx_Xn1(long sv, long n, ulong p)
      94             : {
      95        3135 :   GEN t = cgetg(n+3, t_VECSMALL);
      96             :   long i;
      97             :   (void) p;
      98        3135 :   t[1] = sv;
      99        3135 :   t[2] = 1;
     100        9718 :   for (i = 3; i <= n+1; i++) t[i] = 0;
     101        3135 :   t[i] = 1; return t;
     102             : }
     103             : 
     104             : static GEN
     105       89430 : Flx_root_mod_2(GEN f)
     106             : {
     107       89430 :   int z1, z0 = !(f[2] & 1);
     108             :   long i,n;
     109             :   GEN y;
     110             : 
     111      343641 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     112       89430 :   z1 = n & 1;
     113       89430 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     114       89428 :   if (z0) y[i++] = 0;
     115       89428 :   if (z1) y[i  ] = 1;
     116       89428 :   return y;
     117             : }
     118             : static ulong
     119          91 : Flx_oneroot_mod_2(GEN f)
     120             : {
     121             :   long i,n;
     122          91 :   if (!(f[2] & 1)) return 0;
     123         364 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     124          91 :   if (n & 1) return 1;
     125          28 :   return 2;
     126             : }
     127             : 
     128             : static GEN FpX_roots_i(GEN f, GEN p);
     129             : static GEN Flx_roots_i(GEN f, ulong p);
     130             : 
     131             : static int
     132    16049042 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     133             : 
     134             : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
     135             :  * pp is a small prime.
     136             :  * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
     137             :  * assume that f is an FpX, pp a prime and return t_INTs */
     138             : static GEN
     139      354476 : rootmod_aux(GEN f, GEN pp)
     140             : {
     141             :   GEN y;
     142      354476 :   switch(lg(f))
     143             :   {
     144          14 :     case 2: pari_err_ROOTS0("rootmod");
     145       46221 :     case 3: return cgetg(1,t_COL);
     146             :   }
     147      308241 :   if (typ(f) == t_VECSMALL)
     148             :   {
     149      305105 :     ulong p = pp[2];
     150      305105 :     if (p == 2)
     151       89423 :       y = Flx_root_mod_2(f);
     152             :     else
     153             :     {
     154      215682 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     155      215682 :       y = Flx_roots_i(f, p);
     156             :     }
     157      305099 :     y = Flc_to_ZC(y);
     158             :   }
     159             :   else
     160        3136 :     y = FpX_roots_i(f, pp);
     161      308229 :   return y;
     162             : }
     163             : /* assume that f is a ZX and p a prime */
     164             : GEN
     165      354479 : FpX_roots(GEN f, GEN p)
     166             : {
     167      354479 :   pari_sp av = avma;
     168      354479 :   GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
     169      354450 :   return gerepileupto(av, y);
     170             : }
     171             : 
     172             : /* assume x reduced mod p > 2, monic. */
     173             : static int
     174          21 : FpX_quad_factortype(GEN x, GEN p)
     175             : {
     176          21 :   GEN b = gel(x,3), c = gel(x,2);
     177          21 :   GEN D = subii(sqri(b), shifti(c,2));
     178          21 :   return kronecker(D,p);
     179             : }
     180             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     181             : static GEN
     182        7706 : FpX_quad_root(GEN x, GEN p, int unknown)
     183             : {
     184        7706 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     185             : 
     186        7706 :   if (absequaliu(p, 2)) {
     187           0 :     if (!signe(b)) return c;
     188           0 :     return signe(c)? NULL: gen_1;
     189             :   }
     190        7706 :   D = subii(sqri(b), shifti(c,2));
     191        7706 :   D = remii(D,p);
     192        7706 :   if (unknown && kronecker(D,p) == -1) return NULL;
     193             : 
     194        7172 :   s = Fp_sqrt(D,p);
     195             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     196        7172 :   if (!s) return NULL;
     197        7164 :   return Fp_halve(Fp_sub(s,b, p), p);
     198             : }
     199             : static GEN
     200        3222 : FpX_otherroot(GEN x, GEN r, GEN p)
     201        3222 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     202             : 
     203             : /* disc(x^2+bx+c) = b^2 - 4c */
     204             : static ulong
     205    25286197 : Fl_disc_bc(ulong b, ulong c, ulong p)
     206    25286197 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     207             : /* p > 2 */
     208             : static ulong
     209    22821997 : Flx_quad_root(GEN x, ulong p, int unknown)
     210             : {
     211    22821997 :   ulong s, b = x[3], c = x[2];
     212    22821997 :   ulong D = Fl_disc_bc(b, c, p);
     213    22797285 :   if (unknown && krouu(D,p) == -1) return p;
     214    15257236 :   s = Fl_sqrt(D,p);
     215    15314611 :   if (s==~0UL) return p;
     216    15314598 :   return Fl_halve(Fl_sub(s,b, p), p);
     217             : }
     218             : static ulong
     219    13832601 : Flx_otherroot(GEN x, ulong r, ulong p)
     220    13832601 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     221             : 
     222             : /* 'todo' contains the list of factors to be split.
     223             :  * 'done' the list of finished factors, no longer touched */
     224             : struct split_t { GEN todo, done; };
     225             : static void
     226     5022311 : split_init(struct split_t *S, long max)
     227             : {
     228     5022311 :   S->todo = vectrunc_init(max);
     229     5021857 :   S->done = vectrunc_init(max);
     230     5021446 : }
     231             : #if 0
     232             : /* move todo[i] to done */
     233             : static void
     234             : split_convert(struct split_t *S, long i)
     235             : {
     236             :   long n = lg(S->todo)-1;
     237             :   vectrunc_append(S->done, gel(S->todo,i));
     238             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     239             :   setlg(S->todo, n);
     240             : }
     241             : #endif
     242             : /* append t to todo */
     243             : static void
     244     5336575 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     245             : /* delete todo[i], add t to done */
     246             : static void
     247     5336826 : split_moveto_done(struct split_t *S, long i, GEN t)
     248             : {
     249     5336826 :   long n = lg(S->todo)-1;
     250     5336826 :   vectrunc_append(S->done, t);
     251     5337304 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     252     5337304 :   setlg(S->todo, n);
     253             : 
     254     5337193 : }
     255             : /* append t to done */
     256             : static void
     257      484551 : split_add_done(struct split_t *S, GEN t)
     258      484551 : { vectrunc_append(S->done, t); }
     259             : /* split todo[i] into a and b */
     260             : static void
     261      404803 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     262             : {
     263      404803 :   gel(S->todo, i) = a;
     264      404803 :   split_add(S, b);
     265      404804 : }
     266             : /* split todo[i] into a and b, moved to done */
     267             : static void
     268      449281 : split_done(struct split_t *S, long i, GEN a, GEN b)
     269             : {
     270      449281 :   split_moveto_done(S, i, a);
     271      449287 :   split_add_done(S, b);
     272      449287 : }
     273             : 
     274             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     275             : static GEN
     276        3136 : FpX_roots_i(GEN f, GEN p)
     277             : {
     278             :   GEN pol, pol0, a, q;
     279             :   struct split_t S;
     280             : 
     281        3136 :   split_init(&S, lg(f)-1);
     282        3136 :   settyp(S.done, t_COL);
     283        3136 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     284        3136 :   switch(degpol(f))
     285             :   {
     286           7 :     case 0: return ZC_copy(S.done);
     287          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     288        1744 :     case 2: {
     289        1744 :       GEN s, r = FpX_quad_root(f, p, 1);
     290        1744 :       if (r) {
     291        1744 :         split_add_done(&S, r);
     292        1744 :         s = FpX_otherroot(f,r, p);
     293             :         /* f not known to be square free yet */
     294        1744 :         if (!equalii(r, s)) split_add_done(&S, s);
     295             :       }
     296        1744 :       return sort(S.done);
     297             :     }
     298             :   }
     299             : 
     300        1371 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     301        1371 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     302        1371 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     303        1371 :   a = FpX_gcd(f,a, p);
     304        1371 :   if (!degpol(a)) return ZC_copy(S.done);
     305        1371 :   split_add(&S, FpX_normalize(a,p));
     306             : 
     307        1371 :   q = shifti(p,-1);
     308        1371 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     309        1371 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     310        1371 :   for (pol0[2] = 1;; pol0[2]++)
     311        1501 :   {
     312        2872 :     long j, l = lg(S.todo);
     313        2872 :     if (l == 1) return sort(S.done);
     314        1508 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     315        3132 :     for (j = 1; j < l; j++)
     316             :     {
     317        1631 :       GEN b, r, s, c = gel(S.todo,j);
     318        1631 :       switch(degpol(c))
     319             :       { /* convert linear and quadratics to roots, try to split the rest */
     320         196 :         case 1:
     321         196 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     322         196 :           j--; l--; break;
     323        1298 :         case 2:
     324        1298 :           r = FpX_quad_root(c, p, 0);
     325        1298 :           if (!r) pari_err_PRIME("polrootsmod",p);
     326        1291 :           s = FpX_otherroot(c,r, p);
     327        1291 :           split_done(&S, j, r, s);
     328        1291 :           j--; l--; break;
     329         137 :         default:
     330         137 :           b = FpXQ_pow(pol,q, c,p);
     331         137 :           if (degpol(b) <= 0) continue;
     332         123 :           b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
     333         123 :           if (!degpol(b)) continue;
     334         123 :           b = FpX_normalize(b, p);
     335         123 :           c = FpX_div(c,b, p);
     336         123 :           split_todo(&S, j, b, c);
     337             :       }
     338             :     }
     339             :   }
     340             : }
     341             : 
     342             : /* Assume f is normalized */
     343             : static ulong
     344      180228 : Flx_cubic_root(GEN ff, ulong p)
     345             : {
     346      180228 :   GEN f = Flx_normalize(ff,p);
     347      180228 :   ulong pi = get_Fl_red(p);
     348      180228 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     349      180228 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     350      180229 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     351      180228 :   ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
     352      180229 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     353      180229 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     354      180229 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     355      180229 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     356      180229 :   if (krouu(D,p) >= 0)
     357             :   {
     358      115363 :     ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     359      115363 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     360      115363 :     if (p%3==2) /* 1 solutions */
     361       20567 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     362             :     else
     363             :     {
     364       94796 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     365       94794 :       if (vS1==~0UL) return p; /*0 solutions*/
     366             :       /*3 solutions*/
     367             :     }
     368       83875 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     369       83877 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     370             :   }
     371             :   else
     372             :   {
     373       64866 :     pari_sp av = avma;
     374       64866 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     375       64866 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     376             :     ulong Sa;
     377       64866 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     378       64866 :     Sa = vS1[1];
     379       64866 :     if (p%3==1) /*1 solutions*/
     380             :     {
     381       24235 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     382       24235 :       if (Fa!=P)
     383       15952 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     384             :     }
     385       64866 :     set_avma(av);
     386       64866 :     return Fl_sub(Fl_double(Sa,p),t,p);
     387             :   }
     388             : }
     389             : 
     390             : /* Assume f is normalized */
     391             : static GEN
     392         126 : FpX_cubic_root(GEN ff, GEN p)
     393             : {
     394         126 :   GEN f = FpX_normalize(ff,p);
     395         126 :   GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
     396         126 :   ulong pm3 = umodiu(p,3);
     397          28 :   GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
     398         126 :                   : diviuexact(addiu(p,1),3);
     399         126 :   GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
     400         126 :   GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
     401         126 :   GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
     402         126 :   GEN A3 =  Fp_mul(A, p3, p);
     403         126 :   GEN A32 = Fp_sqr(A3, p), A33 = Fp_mul(A3, A32, p);
     404         126 :   GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
     405         126 :   GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
     406         126 :   if (kronecker(D,p) >= 0)
     407             :   {
     408          28 :     GEN s = Fp_sqrt(D, p), vS1, vS2;
     409          28 :     GEN S1 = S==s ? S: Fp_halve(Fp_sub(S, s, p), p);
     410          28 :     if (pm3 == 2) /* 1 solutions */
     411           0 :       vS1 = Fp_pow(S1, diviuexact(addiu(shifti(p, 1), 1), 3), p);
     412             :     else
     413             :     {
     414          28 :       vS1 = Fp_sqrtn(S1, utoi(3), p, NULL);
     415          28 :       if (!vS1) return p; /*0 solutions*/
     416             :       /*3 solutions*/
     417             :     }
     418          28 :     vS2 = P? Fp_mul(P, Fp_inv(vS1, p), p): 0;
     419          28 :     return Fp_sub(Fp_add(vS1,vS2, p), t, p);
     420             :   }
     421             :   else
     422             :   {
     423          98 :     pari_sp av = avma;
     424          98 :     GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
     425          98 :     GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
     426          98 :     GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
     427             :     GEN Sa;
     428          98 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     429          98 :     Sa = gel(vS1,2);
     430          98 :     if (pm3 == 1) /*1 solutions*/
     431             :     {
     432           0 :       GEN Fa = FpXQ_norm(vS1, T, p);
     433           0 :       if (!equalii(Fa,P))
     434           0 :         Sa = Fp_mul(Sa, Fp_div(Fa, P, p),p);
     435             :     }
     436          98 :     set_avma(av);
     437          98 :     return Fp_sub(shifti(Sa,1),t,p);
     438             :   }
     439             : }
     440             : 
     441             : /* assume p > 2 prime */
     442             : static ulong
     443     3248591 : Flx_oneroot_i(GEN f, ulong p, long fl)
     444             : {
     445             :   GEN pol, a;
     446             :   ulong q;
     447             :   long da;
     448             : 
     449     3248591 :   if (Flx_val(f)) return 0;
     450     3247807 :   switch(degpol(f))
     451             :   {
     452       13014 :     case 1: return Fl_neg(f[2], p);
     453     2795644 :     case 2: return Flx_quad_root(f, p, 1);
     454      164308 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     455             :   }
     456             : 
     457      274841 :   if (!fl)
     458             :   {
     459      242297 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     460      242297 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     461      242297 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     462      242297 :     a = Flx_gcd(f,a, p);
     463       32544 :   } else a = f;
     464      274841 :   da = degpol(a);
     465      274838 :   if (!da) return p;
     466      196505 :   a = Flx_normalize(a,p);
     467             : 
     468      196505 :   q = p >> 1;
     469      196505 :   pol = polx_Flx(f[1]);
     470      196504 :   for(pol[2] = 1;; pol[2]++)
     471             :   {
     472      299873 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     473      299874 :     switch(da)
     474             :     {
     475      123298 :       case 1: return Fl_neg(a[2], p);
     476       57293 :       case 2: return Flx_quad_root(a, p, 0);
     477       15920 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     478             :       default: {
     479      103363 :         GEN b = Flxq_powu(pol,q, a,p);
     480             :         long db;
     481      103347 :         if (degpol(b) <= 0) continue;
     482       98391 :         b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
     483       98395 :         db = degpol(b); if (!db) continue;
     484       98394 :         b = Flx_normalize(b, p);
     485       98402 :         if (db <= (da >> 1)) {
     486       60388 :           a = b;
     487       60388 :           da = db;
     488             :         } else {
     489       38014 :           a = Flx_div(a,b, p);
     490       38021 :           da -= db;
     491             :         }
     492             :       }
     493             :     }
     494             :   }
     495             : }
     496             : 
     497             : /* assume p > 3 prime */
     498             : static GEN
     499        5245 : FpX_oneroot_i(GEN f, GEN p)
     500             : {
     501             :   GEN pol, pol0, a, q;
     502             :   long da;
     503             : 
     504        5245 :   if (ZX_val(f)) return gen_0;
     505        4923 :   switch(degpol(f))
     506             :   {
     507         814 :     case 1: return subii(p, gel(f,2));
     508        3914 :     case 2: return FpX_quad_root(f, p, 1);
     509         126 :     case 3: return FpX_cubic_root(f, p);
     510             :   }
     511             : 
     512          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     513          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     514          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     515          70 :   a = FpX_gcd(f,a, p);
     516          70 :   da = degpol(a);
     517          70 :   if (!da) return NULL;
     518          70 :   a = FpX_normalize(a,p);
     519             : 
     520          70 :   q = shifti(p,-1);
     521          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     522          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     523          70 :   for (pol0[2]=1; ; pol0[2]++)
     524             :   {
     525         224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     526         224 :     switch(da)
     527             :     {
     528          42 :       case 1: return subii(p, gel(a,2));
     529          28 :       case 2: return FpX_quad_root(a, p, 0);
     530         154 :       default: {
     531         154 :         GEN b = FpXQ_pow(pol,q, a,p);
     532             :         long db;
     533         154 :         if (degpol(b) <= 0) continue;
     534         147 :         b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
     535         147 :         db = degpol(b); if (!db) continue;
     536         147 :         b = FpX_normalize(b, p);
     537         147 :         if (db <= (da >> 1)) {
     538         105 :           a = b;
     539         105 :           da = db;
     540             :         } else {
     541          42 :           a = FpX_div(a,b, p);
     542          42 :           da -= db;
     543             :         }
     544             :       }
     545             :     }
     546             :   }
     547             : }
     548             : 
     549             : ulong
     550     3206776 : Flx_oneroot(GEN f, ulong p)
     551             : {
     552     3206776 :   pari_sp av = avma;
     553             :   ulong r;
     554     3206776 :   switch(lg(f))
     555             :   {
     556           0 :     case 2: return 0;
     557           0 :     case 3: return p;
     558             :   }
     559     3206776 :   if (p == 2) return Flx_oneroot_mod_2(f);
     560     3206776 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     561     3206776 :   return gc_ulong(av,r);
     562             : }
     563             : 
     564             : ulong
     565       34006 : Flx_oneroot_split(GEN f, ulong p)
     566             : {
     567       34006 :   pari_sp av = avma;
     568             :   ulong r;
     569       34006 :   switch(lg(f))
     570             :   {
     571           0 :     case 2: return 0;
     572           0 :     case 3: return p;
     573             :   }
     574       34006 :   if (p == 2) return Flx_oneroot_mod_2(f);
     575       34006 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     576       34010 :   return gc_ulong(av,r);
     577             : }
     578             : 
     579             : /* assume that p is prime */
     580             : GEN
     581       13146 : FpX_oneroot(GEN f, GEN pp) {
     582       13146 :   pari_sp av = avma;
     583       13146 :   ZX_rootmod_init(&f, pp);
     584       13145 :   switch(lg(f))
     585             :   {
     586           0 :     case 2: set_avma(av); return gen_0;
     587           0 :     case 3: return gc_NULL(av);
     588             :   }
     589       13145 :   if (typ(f) == t_VECSMALL)
     590             :   {
     591        7900 :     ulong r, p = pp[2];
     592        7900 :     if (p == 2)
     593          91 :       r = Flx_oneroot_mod_2(f);
     594             :     else
     595        7809 :       r = Flx_oneroot_i(f, p, 0);
     596        7900 :     set_avma(av);
     597        7900 :     return (r == p)? NULL: utoi(r);
     598             :   }
     599        5245 :   f = FpX_oneroot_i(f, pp);
     600        5246 :   if (!f) return gc_NULL(av);
     601        5246 :   return gerepileuptoint(av, f);
     602             : }
     603             : 
     604             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     605             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     606             : /* returned in n                                                          */
     607             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     608             : /* twice as big and decrement until it divides p-1.                       */
     609             : static GEN
     610         154 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     611             : {
     612         154 :    pari_sp ltop = avma;
     613             :    GEN pm, factn, power, base, zeta;
     614             :    long n;
     615             : 
     616         154 :    pm = subis (p, 1ul);
     617         336 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     618         154 :    factn = Z_factor(stoi(n));
     619         154 :    power = diviuexact (pm, n);
     620         154 :    base = gen_1;
     621             :    do {
     622         259 :       base = addis (base, 1l);
     623         259 :       zeta = Fp_pow (base, power, p);
     624             :    }
     625         259 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     626         154 :    *pt_n = n;
     627         154 :    return gerepileuptoint (ltop, zeta);
     628             : }
     629             : 
     630             : GEN
     631        1092 : FpX_oneroot_split(GEN fact, GEN p)
     632             : {
     633        1092 :   pari_sp av = avma;
     634             :   long n, deg_f, i, dmin;
     635             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     636        1092 :   fact = FpX_normalize(fact, p);
     637        1092 :   deg_f = degpol(fact);
     638        1092 :   if (deg_f <= 3) return FpX_oneroot(fact, p);
     639         133 :   minfactor = fact; /* factor of minimal degree found so far */
     640         133 :   dmin = degpol(minfactor);
     641         133 :   xplusa = pol_x(varn(fact));
     642         287 :   while (dmin > 3)
     643             :   {
     644             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     645             :     /* zeta varies over the roots of unity in F_p                         */
     646         154 :     fact = minfactor; deg_f = dmin;
     647         154 :     zeta = gen_1;
     648         154 :     prim = good_root_of_unity(p, deg_f, 1, &n);
     649         154 :     expo = diviuexact(subiu(p, 1), n);
     650             :     /* update X+a, avoid a=0 */
     651         154 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     652         154 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     653         329 :     for (i = 0; i < n; i++)
     654             :     {
     655         252 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     656         252 :       long dtmp = degpol(tmp);
     657         252 :       if (dtmp > 0 && dtmp < deg_f)
     658             :       {
     659         161 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     660         161 :         if (dtmp < dmin)
     661             :         {
     662         154 :           minfactor = FpX_normalize (tmp, p);
     663         154 :           dmin = dtmp;
     664         154 :           if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
     665             :             /* stop early to avoid too many gcds */
     666             :             break;
     667             :         }
     668             :       }
     669         175 :       zeta = Fp_mul (zeta, prim, p);
     670             :     }
     671             :   }
     672         133 :   return gerepileuptoint(av, FpX_oneroot(minfactor, p));
     673             : }
     674             : 
     675             : /*******************************************************************/
     676             : /*                                                                 */
     677             : /*                     FACTORISATION MODULO p                      */
     678             : /*                                                                 */
     679             : /*******************************************************************/
     680             : 
     681             : /* F / E  a vector of vectors of factors / exponents of virtual length l
     682             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
     683             : static GEN
     684     1587632 : FE_concat(GEN F, GEN E, long l)
     685             : {
     686     1587632 :   setlg(E,l); E = shallowconcat1(E);
     687     1587628 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
     688             : }
     689             : 
     690             : static GEN
     691          14 : ddf_to_ddf2_i(GEN V, long fl)
     692             : {
     693             :   GEN F, D;
     694          14 :   long i, j, l = lg(V);
     695          14 :   F = cgetg(l, t_VEC);
     696          14 :   D = cgetg(l, t_VECSMALL);
     697         112 :   for (i = j = 1; i < l; i++)
     698             :   {
     699          98 :     GEN Vi = gel(V,i);
     700          98 :     if ((fl==2 && F2x_degree(Vi) == 0)
     701          98 :       ||(fl==0 && degpol(Vi) == 0)) continue;
     702          35 :     gel(F,j) = Vi;
     703          35 :     uel(D,j) = i; j++;
     704             :   }
     705          14 :   setlg(F,j);
     706          14 :   setlg(D,j); return mkvec2(F,D);
     707             : }
     708             : 
     709             : GEN
     710           7 : ddf_to_ddf2(GEN V)
     711           7 : { return ddf_to_ddf2_i(V, 0); }
     712             : 
     713             : static GEN
     714           7 : F2x_ddf_to_ddf2(GEN V)
     715           7 : { return ddf_to_ddf2_i(V, 2); }
     716             : 
     717             : GEN
     718     5368745 : vddf_to_simplefact(GEN V, long d)
     719             : {
     720             :   GEN E, F;
     721     5368745 :   long i, j, c, l = lg(V);
     722     5368745 :   F = cgetg(d+1, t_VECSMALL);
     723     5368092 :   E = cgetg(d+1, t_VECSMALL);
     724    10821622 :   for (i = c = 1; i < l; i++)
     725             :   {
     726     5453037 :     GEN Vi = gel(V,i);
     727     5453037 :     long l = lg(Vi);
     728    27713413 :     for (j = 1; j < l; j++)
     729             :     {
     730    22259413 :       long k, n = degpol(gel(Vi,j)) / j;
     731    33619515 :       for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
     732             :     }
     733             :   }
     734     5368585 :   setlg(F,c);
     735     5368656 :   setlg(E,c);
     736     5368399 :   return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
     737             : }
     738             : 
     739             : /* product of terms of degree 1 in factorization of f */
     740             : GEN
     741      182695 : FpX_split_part(GEN f, GEN p)
     742             : {
     743      182695 :   long n = degpol(f);
     744      182695 :   GEN z, X = pol_x(varn(f));
     745      182695 :   if (n <= 1) return f;
     746      180762 :   f = FpX_red(f, p);
     747      180761 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     748      180762 :   return FpX_gcd(z,f,p);
     749             : }
     750             : 
     751             : /* Compute the number of roots in Fp without counting multiplicity
     752             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     753             : long
     754       99275 : FpX_nbroots(GEN f, GEN p)
     755             : {
     756       99275 :   pari_sp av = avma;
     757       99275 :   GEN z = FpX_split_part(f, p);
     758       99275 :   return gc_long(av, degpol(z));
     759             : }
     760             : 
     761             : /* 1 < deg(f) <= p */
     762             : static int
     763       80569 : Flx_is_totally_split_i(GEN f, ulong p)
     764             : {
     765       80569 :   GEN F = Flx_Frobenius(f, p);
     766       80574 :   return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
     767             : }
     768             : int
     769       80576 : Flx_is_totally_split(GEN f, ulong p)
     770             : {
     771       80576 :   pari_sp av = avma;
     772       80576 :   ulong n = degpol(f);
     773       80576 :   if (n <= 1) return 1;
     774       80569 :   if (n > p) return 0; /* includes n < 0 */
     775       80569 :   return gc_bool(av, Flx_is_totally_split_i(f,p));
     776             : }
     777             : int
     778           0 : FpX_is_totally_split(GEN f, GEN p)
     779             : {
     780           0 :   pari_sp av = avma;
     781           0 :   ulong n = degpol(f);
     782             :   int u;
     783           0 :   if (n <= 1) return 1;
     784           0 :   if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
     785           0 :   if (lgefint(p) != 3)
     786           0 :     u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
     787             :   else
     788             :   {
     789           0 :     ulong pp = (ulong)p[2];
     790           0 :     u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
     791             :   }
     792           0 :   return gc_bool(av, u);
     793             : }
     794             : 
     795             : long
     796     2702744 : Flx_nbroots(GEN f, ulong p)
     797             : {
     798     2702744 :   long n = degpol(f);
     799     2702745 :   pari_sp av = avma;
     800             :   GEN z;
     801     2702745 :   if (n <= 1) return n;
     802     2689263 :   if (n == 2)
     803             :   {
     804             :     ulong D;
     805       62461 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     806       61467 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     807       61467 :     return 1 + krouu(D,p);
     808             :   }
     809     2626802 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     810     2626799 :   z = Flx_gcd(z, f, p);
     811     2626806 :   return gc_long(av, degpol(z));
     812             : }
     813             : 
     814             : long
     815        4403 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
     816             : {
     817        4403 :   pari_sp av = avma;
     818             :   GEN X, b, g, xq;
     819             :   long i, j, n, v, B, l, m;
     820             :   pari_timer ti;
     821             :   hashtable h;
     822             : 
     823        4403 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     824        4403 :   X = pol_x(v);
     825        4403 :   if (ZX_equal(X,XP)) return 1;
     826        4403 :   B = n/2;
     827        4403 :   l = usqrt(B);
     828        4403 :   m = (B+l-1)/l;
     829        4403 :   T = FpX_get_red(T, p);
     830        4403 :   hash_init_GEN(&h, l+2, ZX_equal, 1);
     831        4403 :   hash_insert_long(&h, X,  0);
     832        4403 :   hash_insert_long(&h, XP, 1);
     833        4403 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     834        4403 :   b = XP;
     835        4403 :   xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
     836        4403 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
     837       10402 :   for (i = 3; i <= l+1; i++)
     838             :   {
     839        6790 :     b = FpX_FpXQV_eval(b, xq, T, p);
     840        6790 :     if (gequalX(b)) return gc_long(av,i-1);
     841        5999 :     hash_insert_long(&h, b, i-1);
     842             :   }
     843        3612 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
     844        3612 :   g = b;
     845        3612 :   xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
     846        3612 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
     847       12341 :   for(i = 2; i <= m+1; i++)
     848             :   {
     849       10745 :     g = FpX_FpXQV_eval(g, xq, T, p);
     850       10745 :     if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
     851             :   }
     852        1596 :   return gc_long(av,n);
     853             : }
     854             : 
     855             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     856             : static GEN
     857         911 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
     858             : {
     859             :   GEN b, g, h, F, f, Tr, xq;
     860             :   long i, j, n, v, B, l, m;
     861             :   pari_timer ti;
     862             : 
     863         911 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     864         911 :   if (n == 0) return cgetg(1, t_VEC);
     865         911 :   if (n == 1) return mkvec(get_FpX_mod(T));
     866         805 :   B = n/2;
     867         805 :   l = usqrt(B);
     868         805 :   m = (B+l-1)/l;
     869         805 :   T = FpX_get_red(T, p);
     870         805 :   b = cgetg(l+2, t_VEC);
     871         805 :   gel(b, 1) = pol_x(v);
     872         805 :   gel(b, 2) = XP;
     873         805 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     874         805 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     875         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
     876        1065 :   for (i = 3; i <= l+1; i++)
     877         260 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     878         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
     879         805 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     880         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
     881         805 :   g = cgetg(m+1, t_VEC);
     882         805 :   gel(g, 1) = gel(xq, 2);
     883        1804 :   for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     884         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
     885         805 :   h = cgetg(m+1, t_VEC);
     886        2609 :   for (j = 1; j <= m; j++)
     887             :   {
     888        1804 :     pari_sp av = avma;
     889        1804 :     GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
     890        3175 :     for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
     891        1804 :     gel(h,j) = gerepileupto(av, e);
     892             :   }
     893         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
     894         805 :   Tr = get_FpX_mod(T);
     895         805 :   F = cgetg(m+1, t_VEC);
     896        2609 :   for (j = 1; j <= m; j++)
     897             :   {
     898        1804 :     GEN u = FpX_gcd(Tr, gel(h,j), p);
     899        1804 :     if (degpol(u))
     900             :     {
     901         409 :       u = FpX_normalize(u, p);
     902         409 :       Tr = FpX_div(Tr, u, p);
     903             :     }
     904        1804 :     gel(F,j) = u;
     905             :   }
     906         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
     907         805 :   f = const_vec(n, pol_1(v));
     908        2609 :   for (j = 1; j <= m; j++)
     909             :   {
     910        1804 :     GEN e = gel(F, j);
     911        1861 :     for (i=l-1; i >= 0; i--)
     912             :     {
     913        1861 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     914        1861 :       if (degpol(u))
     915             :       {
     916         417 :         u = FpX_normalize(u, p);
     917         417 :         gel(f, l*j-i) = u;
     918         417 :         e = FpX_div(e, u, p);
     919             :       }
     920        1861 :       if (!degpol(e)) break;
     921             :     }
     922             :   }
     923         805 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
     924         805 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     925         805 :   return f;
     926             : }
     927             : 
     928             : static void
     929           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     930             : {
     931           0 :   long n = degpol(Tp), r = n/d, ct = 0;
     932             :   GEN T, f, ff, p2;
     933           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     934           0 :   p2 = shifti(p,-1);
     935           0 :   T = FpX_get_red(Tp, p);
     936           0 :   XP = FpX_rem(XP, T, p);
     937             :   while (1)
     938           0 :   {
     939           0 :     pari_sp btop = avma;
     940             :     long i;
     941           0 :     GEN g = random_FpX(n, varn(Tp), p);
     942           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     943           0 :     if (signe(t) == 0) continue;
     944           0 :     for(i=1; i<=10; i++)
     945             :     {
     946           0 :       pari_sp btop2 = avma;
     947           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     948           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     949           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     950           0 :       set_avma(btop2);
     951             :     }
     952           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     953           0 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
     954           0 :     set_avma(btop);
     955             :   }
     956           0 :   f = FpX_normalize(f, p);
     957           0 :   ff = FpX_div(Tp, f ,p);
     958           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     959           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     960             : }
     961             : 
     962             : static void
     963        1007 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     964             : {
     965             :   pari_sp av;
     966        1007 :   GEN Tp = get_FpX_mod(T);
     967        1007 :   long n = degpol(hp), vT = varn(Tp), ct = 0;
     968             :   GEN u1, u2, f1, f2, R, h;
     969        1007 :   h = FpX_get_red(hp, p);
     970        1007 :   t = FpX_rem(t, T, p);
     971        1007 :   av = avma;
     972             :   do
     973             :   {
     974        1682 :     set_avma(av);
     975        1682 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     976        1682 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     977        1682 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
     978        1682 :   } while (degpol(u1)==0 || degpol(u1)==n);
     979        1007 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     980        1007 :   f1 = FpX_normalize(f1, p);
     981        1007 :   u2 = FpX_div(hp, u1, p);
     982        1007 :   f2 = FpX_div(Tp, f1, p);
     983        1007 :   if (degpol(u1)==1)
     984         669 :     gel(V, idx) = f1;
     985             :   else
     986         338 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     987        1007 :   idx += degpol(f1)/d;
     988        1007 :   if (degpol(u2)==1)
     989         730 :     gel(V, idx) = f2;
     990             :   else
     991         277 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     992        1007 : }
     993             : 
     994             : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
     995             : static void
     996         392 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     997             : {
     998         392 :   long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
     999             :   GEN T, h, t;
    1000             :   pari_timer ti;
    1001             : 
    1002         392 :   T = FpX_get_red(Tp, p);
    1003         392 :   XP = FpX_rem(XP, T, p);
    1004         392 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1005             :   do
    1006             :   {
    1007         392 :     GEN g = random_FpX(n, vT, p);
    1008         392 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
    1009         392 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
    1010         392 :     h = FpXQ_minpoly(t, T, p);
    1011         392 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
    1012         392 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
    1013         392 :   } while (degpol(h) != r);
    1014         392 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
    1015         392 : }
    1016             : 
    1017             : static GEN
    1018         897 : FpX_factor_Shoup(GEN T, GEN p)
    1019             : {
    1020         897 :   long i, n, s = 0;
    1021             :   GEN XP, D, V;
    1022         897 :   long e = expi(p);
    1023             :   pari_timer ti;
    1024         897 :   n = get_FpX_degree(T);
    1025         897 :   T = FpX_get_red(T, p);
    1026         897 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1027         897 :   XP = FpX_Frobenius(T, p);
    1028         897 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1029         897 :   D = FpX_ddf_Shoup(T, XP, p);
    1030         897 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1031         897 :   s = ddf_to_nbfact(D);
    1032         897 :   V = cgetg(s+1, t_COL);
    1033        7279 :   for (i = 1, s = 1; i <= n; i++)
    1034             :   {
    1035        6382 :     GEN Di = gel(D,i);
    1036        6382 :     long ni = degpol(Di), ri = ni/i;
    1037        6382 :     if (ni == 0) continue;
    1038         930 :     Di = FpX_normalize(Di, p);
    1039         930 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1040         392 :     if (ri <= e*expu(e))
    1041         392 :       FpX_edf(Di, XP, i, p, V, s);
    1042             :     else
    1043           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
    1044         392 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
    1045         392 :     s += ri;
    1046             :   }
    1047         897 :   return V;
    1048             : }
    1049             : 
    1050             : long
    1051     2047063 : ddf_to_nbfact(GEN D)
    1052             : {
    1053     2047063 :   long l = lg(D), i, s = 0;
    1054    13733030 :   for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
    1055     2047064 :   return s;
    1056             : }
    1057             : 
    1058             : /* Yun algorithm: Assume p > degpol(T) */
    1059             : static GEN
    1060        1595 : FpX_factor_Yun(GEN T, GEN p)
    1061             : {
    1062        1595 :   long n = degpol(T), i = 1;
    1063        1595 :   GEN a, b, c, d = FpX_deriv(T, p);
    1064        1595 :   GEN V = cgetg(n+1,t_VEC);
    1065        1595 :   a = FpX_gcd(T, d, p);
    1066        1595 :   if (degpol(a) == 0) return mkvec(T);
    1067         535 :   b = FpX_div(T, a, p);
    1068             :   do
    1069             :   {
    1070        2526 :     c = FpX_div(d, a, p);
    1071        2526 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1072        2526 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1073        2526 :     gel(V, i++) = a;
    1074        2526 :     b = FpX_div(b, a, p);
    1075        2526 :   } while (degpol(b));
    1076         535 :   setlg(V, i); return V;
    1077             : }
    1078             : GEN
    1079      333676 : FpX_factor_squarefree(GEN T, GEN p)
    1080             : {
    1081      333676 :   if (lgefint(p)==3)
    1082             :   {
    1083      332960 :     ulong pp = (ulong)p[2];
    1084      332960 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1085      332960 :     return FlxV_to_ZXV(u);
    1086             :   }
    1087         716 :   return FpX_factor_Yun(T, p);
    1088             : }
    1089             : 
    1090             : GEN
    1091      218470 : FpX_roots_mult(GEN T, long n, GEN p)
    1092             : {
    1093      218470 :   GEN V = FpX_factor_squarefree(T,p), W;
    1094      218470 :   long l = lg(V), i;
    1095      218470 :   if (l<=n) return cgetg(1,t_COL);
    1096       35455 :   W = cgetg(l-n+1,t_VEC);
    1097      116473 :   for (i = n; i < l; i++)
    1098       81018 :     gel(W,i-n+1) = FpX_roots(gel(V,i), p);
    1099       35455 :   return shallowconcat1(W);
    1100             : }
    1101             : 
    1102             : long
    1103         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1104             : {
    1105         168 :   pari_sp av = avma;
    1106             :   GEN lc, F;
    1107         168 :   long i, l, n = degpol(f), v = varn(f);
    1108         168 :   if (n % k) return 0;
    1109         168 :   if (lgefint(p)==3)
    1110             :   {
    1111         126 :     ulong pp = p[2];
    1112         126 :     GEN fp = ZX_to_Flx(f, pp);
    1113         126 :     if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
    1114         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
    1115         105 :     return 1;
    1116             :   }
    1117          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1118          42 :   if (!lc) { av = avma; return 0; }
    1119          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1120        1491 :   for(i=1; i <= l; i++)
    1121        1456 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1122          35 :   if (pt_r)
    1123             :   {
    1124          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1125        1484 :     for (i=l; i>=1; i--)
    1126             :     {
    1127        1449 :       if (i%k) continue;
    1128         294 :       s = FpX_mul(s, gel(F,i), p);
    1129         294 :       r = FpX_mul(r, s, p);
    1130             :     }
    1131          35 :     *pt_r = gerepileupto(av, r);
    1132           0 :   } else av = avma;
    1133          35 :   return 1;
    1134             : }
    1135             : 
    1136             : static GEN
    1137         823 : FpX_factor_Cantor(GEN T, GEN p)
    1138             : {
    1139         823 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1140         823 :   long i, j, l = lg(V);
    1141         823 :   F = cgetg(l, t_VEC);
    1142         823 :   E = cgetg(l, t_VEC);
    1143        2197 :   for (i=1, j=1; i < l; i++)
    1144        1374 :     if (degpol(gel(V,i)))
    1145             :     {
    1146         897 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1147         897 :       gel(F, j) = Fj;
    1148         897 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1149         897 :       j++;
    1150             :     }
    1151         823 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1152             : }
    1153             : 
    1154             : static GEN
    1155           0 : FpX_ddf_i(GEN T, GEN p)
    1156             : {
    1157             :   GEN XP;
    1158           0 :   T = FpX_get_red(T, p);
    1159           0 :   XP = FpX_Frobenius(T, p);
    1160           0 :   return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
    1161             : }
    1162             : 
    1163             : GEN
    1164           7 : FpX_ddf(GEN f, GEN p)
    1165             : {
    1166           7 :   pari_sp av = avma;
    1167             :   GEN F;
    1168           7 :   switch(ZX_factmod_init(&f, p))
    1169             :   {
    1170           7 :     case 0:  F = F2x_ddf(f);
    1171           7 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    1172           0 :     case 1:  F = Flx_ddf(f,p[2]);
    1173           0 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    1174           0 :     default: F = FpX_ddf_i(f,p); break;
    1175             :   }
    1176           7 :   return gerepilecopy(av, F);
    1177             : }
    1178             : 
    1179             : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
    1180             : static GEN
    1181          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1182             : {
    1183             :   GEN V;
    1184             :   long i, l;
    1185          14 :   if (lgefint(p) == 3)
    1186             :   {
    1187           0 :     ulong pp = p[2];
    1188           0 :     return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
    1189             :   }
    1190          14 :   T = FpX_get_red(T, p);
    1191          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
    1192          28 :   for (i=1; i < l; i++)
    1193          14 :     gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
    1194          14 :   return vddf_to_simplefact(V, get_FpX_degree(T));
    1195             : }
    1196             : 
    1197             : static int
    1198           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1199             : {
    1200           0 :   pari_sp av = avma;
    1201             :   pari_timer ti;
    1202             :   long n;
    1203           0 :   GEN T = get_FpX_mod(Tp);
    1204           0 :   GEN dT = FpX_deriv(T, p);
    1205             :   GEN XP, D;
    1206           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    1207           0 :   n = get_FpX_degree(T);
    1208           0 :   T = FpX_get_red(Tp, p);
    1209           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1210           0 :   XP = FpX_Frobenius(T, p);
    1211           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1212           0 :   D = FpX_ddf_Shoup(T, XP, p);
    1213           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1214           0 :   return gc_bool(av, degpol(gel(D,n)) == n);
    1215             : }
    1216             : 
    1217             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1218             : 
    1219             : /*Assume that p is large and odd*/
    1220             : static GEN
    1221        1695 : FpX_factor_i(GEN f, GEN pp, long flag)
    1222             : {
    1223        1695 :   long d = degpol(f);
    1224        1695 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1225         837 :   switch(flag)
    1226             :   {
    1227         823 :     default: return FpX_factor_Cantor(f, pp);
    1228          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1229           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1230             :   }
    1231             : }
    1232             : 
    1233             : long
    1234           0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1235             : {
    1236           0 :   pari_sp av = avma;
    1237           0 :   long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
    1238           0 :   return gc_long(av,s);
    1239             : }
    1240             : 
    1241             : long
    1242           0 : FpX_nbfact(GEN T, GEN p)
    1243             : {
    1244           0 :   pari_sp av = avma;
    1245           0 :   GEN XP = FpX_Frobenius(T, p);
    1246           0 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1247           0 :   return gc_long(av,n);
    1248             : }
    1249             : 
    1250             : /* p > 2 */
    1251             : static GEN
    1252           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1253             : {
    1254           7 :   switch(d)
    1255             :   {
    1256           0 :     case -1:
    1257           0 :     case 0: return NULL;
    1258           0 :     case 1: return gen_1;
    1259             :   }
    1260           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1261             : }
    1262             : /* p > 2 */
    1263             : static GEN
    1264          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1265             : {
    1266          14 :   switch(d)
    1267             :   {
    1268           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1269           0 :     case 0: return trivial_fact();
    1270           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1271             :   }
    1272          14 :   switch(FpX_quad_factortype(f, p)) {
    1273           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1274           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1275           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1276             :   }
    1277             : }
    1278             : 
    1279             : GEN
    1280          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1281             : GEN
    1282      164419 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1283             : 
    1284             : /* not gerepile safe */
    1285             : static GEN
    1286         837 : FpX_factor_2(GEN f, GEN p, long d)
    1287             : {
    1288             :   GEN r, s, R, S;
    1289             :   long v;
    1290             :   int sgn;
    1291         837 :   switch(d)
    1292             :   {
    1293           7 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1294          30 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1295          78 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1296             :   }
    1297         722 :   r = FpX_quad_root(f, p, 1);
    1298         722 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1299         187 :   v = varn(f);
    1300         187 :   s = FpX_otherroot(f, r, p);
    1301         187 :   if (signe(r)) r = subii(p, r);
    1302         187 :   if (signe(s)) s = subii(p, s);
    1303         187 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1304         187 :   R = deg1pol_shallow(gen_1, r, v);
    1305         187 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1306          93 :   S = deg1pol_shallow(gen_1, s, v);
    1307          93 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1308             : }
    1309             : static GEN
    1310         858 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1311             : {
    1312         858 :   switch(flag) {
    1313           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1314          14 :     case 1: return FpX_degfact_2(f, p, d);
    1315         837 :     default: return FpX_factor_2(f, p, d);
    1316             :   }
    1317             : }
    1318             : 
    1319             : static int
    1320      450366 : F2x_quad_factortype(GEN x)
    1321      450366 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1322             : 
    1323             : static GEN
    1324          14 : F2x_is_irred_2(GEN f, long d)
    1325          14 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1326             : 
    1327             : static GEN
    1328       19593 : F2x_degfact_2(GEN f, long d)
    1329             : {
    1330       19593 :   if (!d) return trivial_fact();
    1331       19593 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1332       19404 :   switch(F2x_quad_factortype(f)) {
    1333        5964 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1334        6202 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1335        7238 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1336             :   }
    1337             : }
    1338             : 
    1339             : static GEN
    1340     1165906 : F2x_factor_2(GEN f, long d)
    1341             : {
    1342     1165906 :   long v = f[1];
    1343     1165906 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1344      909186 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1345      426161 :   switch(F2x_quad_factortype(f))
    1346             :   {
    1347       72507 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1348      326235 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1349       27438 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1350             :   }
    1351             : }
    1352             : static GEN
    1353     1185515 : F2x_factor_deg2(GEN f, long d, long flag)
    1354             : {
    1355     1185515 :   switch(flag) {
    1356          14 :     case 2: return F2x_is_irred_2(f, d);
    1357       19593 :     case 1: return F2x_degfact_2(f, d);
    1358     1165908 :     default: return F2x_factor_2(f, d);
    1359             :   }
    1360             : }
    1361             : 
    1362             : /* xt = NULL or x^(p-1)/2 mod g */
    1363             : static void
    1364       27725 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
    1365             : {
    1366       27725 :   ulong q = p >> 1;
    1367       27725 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1368       27725 :   long d = degpol(a);
    1369       27725 :   if (d < 0)
    1370             :   {
    1371             :     ulong i;
    1372         611 :     split_add_done(S, (GEN)1);
    1373        2454 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1374             :   } else {
    1375       27114 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1376       27114 :     if (d)
    1377             :     {
    1378       27079 :       if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
    1379       27079 :       a = Flx_gcd(a, xt, p);
    1380       27079 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1381             :     }
    1382             :   }
    1383       27725 : }
    1384             : static void
    1385       27725 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
    1386             : {
    1387       27725 :   ulong q = p >> 1;
    1388       27725 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1389       27725 :   long d = degpol(a);
    1390       27725 :   if (d < 0)
    1391             :   {
    1392         167 :     ulong i, z = nonsquare_Fl(p);
    1393         167 :     split_add_done(S, (GEN)z);
    1394         306 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1395             :   } else {
    1396       27558 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1397       27558 :     if (d)
    1398             :     {
    1399       27026 :       if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
    1400       27026 :       a = Flx_gcd(a, xt, p);
    1401       27025 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1402             :     }
    1403             :   }
    1404       27725 : }
    1405             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1406             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1407             : static int
    1408     5018248 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1409             : {
    1410     5018248 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1411     5018060 :   long d = degpol(g);
    1412     5017955 :   if (d < 0) return 0;
    1413     5017521 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1414     5017533 :   if (!d) return 1;
    1415     4989018 :   if ((p >> 4) <= (ulong)d)
    1416             :   { /* small p; split directly using x^((p-1)/2) +/- 1 */
    1417       23891 :     GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
    1418       27725 :                                 : NULL;
    1419       27725 :     split_squares(S, g, p, xt);
    1420       27725 :     split_nonsquares(S, g, p, xt);
    1421             :   } else { /* large p; use x^(p-1) - 1 directly */
    1422     4961293 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1423     4947902 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1424     4947902 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1425     4949997 :     g = Flx_gcd(g,a, p);
    1426     4954920 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1427             :   }
    1428     4987814 :   return 1;
    1429             : }
    1430             : 
    1431             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1432             : static GEN
    1433    23733161 : Flx_roots_i(GEN f, ulong p)
    1434             : {
    1435             :   GEN pol, g;
    1436    23733161 :   long v = Flx_valrem(f, &g);
    1437             :   ulong q;
    1438             :   struct split_t S;
    1439             : 
    1440             :   /* optimization: test for small degree first */
    1441    23744146 :   switch(degpol(g))
    1442             :   {
    1443       99469 :     case 1: {
    1444       99469 :       ulong r = p - g[2];
    1445       99469 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1446             :     }
    1447    18614932 :     case 2: {
    1448    18614932 :       ulong r = Flx_quad_root(g, p, 1), s;
    1449    18604954 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1450    12685477 :       s = Flx_otherroot(g,r, p);
    1451    12771476 :       if (r < s)
    1452     3196349 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1453     9575127 :       else if (r > s)
    1454     9598952 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1455             :       else
    1456        4554 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1457             :     }
    1458             :   }
    1459     5019275 :   q = p >> 1;
    1460     5019275 :   split_init(&S, lg(f)-1);
    1461     5018253 :   settyp(S.done, t_VECSMALL);
    1462     5018253 :   if (v) split_add_done(&S, (GEN)0);
    1463     5018253 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1464         434 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1465     5016206 :   pol = polx_Flx(f[1]);
    1466     5017061 :   for (pol[2]=1; ; pol[2]++)
    1467     5371602 :   {
    1468    10388663 :     long j, l = lg(S.todo);
    1469    10388663 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1470     5371324 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1471    11224779 :     for (j = 1; j < l; j++)
    1472             :     {
    1473     5853177 :       GEN b, c = gel(S.todo,j);
    1474             :       ulong r, s;
    1475     5853177 :       switch(degpol(c))
    1476             :       {
    1477     4887282 :         case 1:
    1478     4887282 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1479     4887712 :           j--; l--; break;
    1480      447988 :         case 2:
    1481      447988 :           r = Flx_quad_root(c, p, 0);
    1482      447995 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1483      447988 :           s = Flx_otherroot(c,r, p);
    1484      447989 :           split_done(&S, j, (GEN)r, (GEN)s);
    1485      447996 :           j--; l--; break;
    1486      517799 :         default:
    1487      517799 :           b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
    1488      517895 :           if (degpol(b) <= 0) continue;
    1489      405132 :           b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
    1490      405148 :           if (!degpol(b)) continue;
    1491      404683 :           b = Flx_normalize(b, p);
    1492      404692 :           c = Flx_div(c,b, p);
    1493      404681 :           split_todo(&S, j, b, c);
    1494             :       }
    1495             :     }
    1496             :   }
    1497             : }
    1498             : 
    1499             : GEN
    1500    23526617 : Flx_roots(GEN f, ulong p)
    1501             : {
    1502    23526617 :   pari_sp av = avma;
    1503    23526617 :   switch(lg(f))
    1504             :   {
    1505           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1506           0 :     case 3: set_avma(av); return cgetg(1, t_VECSMALL);
    1507             :   }
    1508    23527619 :   if (p == 2) return Flx_root_mod_2(f);
    1509    23527612 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1510             : }
    1511             : 
    1512             : /* assume x reduced mod p, monic. */
    1513             : static int
    1514     2467329 : Flx_quad_factortype(GEN x, ulong p)
    1515             : {
    1516     2467329 :   ulong b = x[3], c = x[2];
    1517     2467329 :   return krouu(Fl_disc_bc(b, c, p), p);
    1518             : }
    1519             : static GEN
    1520          56 : Flx_is_irred_2(GEN f, ulong p, long d)
    1521             : {
    1522          56 :   if (!d) return NULL;
    1523          56 :   if (d == 1) return gen_1;
    1524          56 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1525             : }
    1526             : static GEN
    1527     2491533 : Flx_degfact_2(GEN f, ulong p, long d)
    1528             : {
    1529     2491533 :   if (!d) return trivial_fact();
    1530     2491533 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1531     2467439 :   switch(Flx_quad_factortype(f, p)) {
    1532     1171057 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1533     1269130 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1534       27829 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1535             :   }
    1536             : }
    1537             : /* p > 2 */
    1538             : static GEN
    1539     1221206 : Flx_factor_2(GEN f, ulong p, long d)
    1540             : {
    1541             :   ulong r, s;
    1542             :   GEN R,S;
    1543     1221206 :   long v = f[1];
    1544     1221206 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1545     1144764 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1546      916507 :   r = Flx_quad_root(f, p, 1);
    1547      916537 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1548      615314 :   s = Flx_otherroot(f, r, p);
    1549      615328 :   r = Fl_neg(r, p);
    1550      615326 :   s = Fl_neg(s, p);
    1551      615324 :   if (s < r) lswap(s,r);
    1552      615324 :   R = mkvecsmall3(v,r,1);
    1553      615323 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1554      506424 :   S = mkvecsmall3(v,s,1);
    1555      506428 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1556             : }
    1557             : static GEN
    1558     3712084 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1559             : {
    1560     3712084 :   switch(flag) {
    1561          56 :     case 2: return Flx_is_irred_2(f, p, d);
    1562     2491619 :     case 1: return Flx_degfact_2(f, p, d);
    1563     1220409 :     default: return Flx_factor_2(f, p, d);
    1564             :   }
    1565             : }
    1566             : 
    1567             : static GEN
    1568       23359 : F2x_Berlekamp_ker(GEN u)
    1569             : {
    1570       23359 :   pari_sp ltop=avma;
    1571       23359 :   long j,N = F2x_degree(u);
    1572             :   GEN Q;
    1573             :   pari_timer T;
    1574       23356 :   timer_start(&T);
    1575       23359 :   Q = F2x_matFrobenius(u);
    1576      322302 :   for (j=1; j<=N; j++)
    1577      298943 :     F2m_flip(Q,j,j);
    1578       23359 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
    1579       23359 :   Q = F2m_ker_sp(Q,0);
    1580       23359 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
    1581       23359 :   return gerepileupto(ltop,Q);
    1582             : }
    1583             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    1584             : static long
    1585       35947 : F2x_split_Berlekamp(GEN *t)
    1586             : {
    1587       35947 :   GEN u = *t, a, b, vker;
    1588       35947 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    1589             : 
    1590       35947 :   if (du == 1) return 1;
    1591       27425 :   if (du == 2)
    1592             :   {
    1593        4066 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    1594             :     {
    1595           0 :       t[0] = mkvecsmall2(sv, 2);
    1596           0 :       t[1] = mkvecsmall2(sv, 3);
    1597           0 :       return 2;
    1598             :     }
    1599        4066 :     return 1;
    1600             :   }
    1601             : 
    1602       23359 :   vker = F2x_Berlekamp_ker(u);
    1603       23359 :   lb = lgcols(vker);
    1604       23377 :   d = lg(vker)-1;
    1605       23377 :   ir = 0;
    1606             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    1607       60811 :   for (L=1; L<d; )
    1608             :   {
    1609             :     GEN pol;
    1610       37454 :     if (d == 2)
    1611        5616 :       pol = F2v_to_F2x(gel(vker,2), sv);
    1612             :     else
    1613             :     {
    1614       31838 :       GEN v = zero_zv(lb);
    1615       31838 :       v[1] = du;
    1616       31838 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    1617      119792 :       for (i=2; i<=d; i++)
    1618       87953 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    1619       31839 :       pol = F2v_to_F2x(v, sv);
    1620             :     }
    1621      112614 :     for (i=ir; i<L && L<d; i++)
    1622             :     {
    1623       75180 :       a = t[i]; la = F2x_degree(a);
    1624       75174 :       if (la == 1) { set_irred(i); }
    1625       74996 :       else if (la == 2)
    1626             :       {
    1627         709 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    1628             :         {
    1629           0 :           t[i] = mkvecsmall2(sv, 2);
    1630           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    1631             :         }
    1632         709 :         set_irred(i);
    1633             :       }
    1634             :       else
    1635             :       {
    1636       74287 :         pari_sp av = avma;
    1637             :         long lb;
    1638       74287 :         b = F2x_rem(pol, a);
    1639       74278 :         if (F2x_degree(b) <= 0) { set_avma(av); continue; }
    1640       27053 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    1641       27060 :         if (lb && lb < la)
    1642             :         {
    1643       27060 :           t[L] = F2x_div(a,b);
    1644       27058 :           t[i]= b; L++;
    1645             :         }
    1646           0 :         else set_avma(av);
    1647             :       }
    1648             :     }
    1649             :   }
    1650       23357 :   return d;
    1651             : }
    1652             : /* assume deg f > 2 */
    1653             : static GEN
    1654       35306 : F2x_Berlekamp_i(GEN f, long flag)
    1655             : {
    1656       35306 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    1657             :   GEN y, E, t, V;
    1658             : 
    1659       35306 :   val = F2x_valrem(f, &f);
    1660       35306 :   if (flag == 2 && val) return NULL;
    1661       35292 :   V = F2x_factor_squarefree(f); lV = lg(V);
    1662       35291 :   if (flag == 2 && lV > 2) return NULL;
    1663             : 
    1664             :   /* to hold factors and exponents */
    1665       35221 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    1666       35223 :   E = cgetg(d+1,t_VECSMALL);
    1667       35223 :   lfact = 1;
    1668       35223 :   if (val) {
    1669       11272 :     if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
    1670       11272 :     E[1] = val; lfact++;
    1671             :   }
    1672             : 
    1673      116286 :   for (k=1; k<lV; k++)
    1674             :   {
    1675       81219 :     if (F2x_degree(gel(V, k))==0) continue;
    1676       35947 :     gel(t,lfact) = gel(V, k);
    1677       35947 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    1678       35945 :     if (flag == 2 && d != 1) return NULL;
    1679       35791 :     if (flag == 1)
    1680       48066 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    1681       98231 :     for (j=0; j<d; j++) E[lfact+j] = k;
    1682       35791 :     lfact += d;
    1683             :   }
    1684       35067 :   if (flag == 2) return gen_1; /* irreducible */
    1685       35053 :   setlg(t, lfact);
    1686       35054 :   setlg(E, lfact); y = mkvec2(t,E);
    1687       24352 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    1688       59406 :               : sort_factor_pol(y, cmpGuGu);
    1689             : }
    1690             : 
    1691             : /* Adapted from Shoup NTL */
    1692             : GEN
    1693      281678 : F2x_factor_squarefree(GEN f)
    1694             : {
    1695             :   GEN r, t, v, tv;
    1696      281678 :   long i, q, n = F2x_degree(f);
    1697      281677 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1698      281682 :   for(q = 1;;q *= 2)
    1699             :   {
    1700      486209 :     r = F2x_gcd(f, F2x_deriv(f));
    1701      486209 :     if (F2x_degree(r) == 0)
    1702             :     {
    1703      232676 :       gel(u, q) = f;
    1704      232676 :       break;
    1705             :     }
    1706      253534 :     t = F2x_div(f, r);
    1707      253531 :     if (F2x_degree(t) > 0)
    1708             :     {
    1709             :       long j;
    1710       95499 :       for(j = 1;;j++)
    1711             :       {
    1712      209289 :         v = F2x_gcd(r, t);
    1713      209288 :         tv = F2x_div(t, v);
    1714      209288 :         if (F2x_degree(tv) > 0)
    1715       97353 :           gel(u, j*q) = tv;
    1716      209288 :         if (F2x_degree(v) <= 0) break;
    1717      113789 :         r = F2x_div(r, v);
    1718      113790 :         t = v;
    1719             :       }
    1720       95499 :       if (F2x_degree(r) == 0) break;
    1721             :     }
    1722      204529 :     f = F2x_sqrt(r);
    1723             :   }
    1724     1136043 :   for (i = n; i; i--)
    1725     1133429 :     if (F2x_degree(gel(u,i))) break;
    1726      281668 :   setlg(u,i+1); return u;
    1727             : }
    1728             : 
    1729             : static GEN
    1730      291326 : F2x_ddf_simple(GEN T, GEN XP)
    1731             : {
    1732      291326 :   pari_sp av = avma, av2;
    1733             :   GEN f, z, Tr, X;
    1734      291326 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1735      291326 :   if (n == 0) return cgetg(1, t_VEC);
    1736      291326 :   if (n == 1) return mkvec(T);
    1737      135248 :   z = XP; Tr = T; X = polx_F2x(v);
    1738      135248 :   f = const_vec(n, pol1_F2x(v));
    1739      135249 :   av2 = avma;
    1740      355299 :   for (j = 1; j <= B; j++)
    1741             :   {
    1742      232370 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1743      232387 :     if (F2x_degree(u))
    1744             :     {
    1745       49140 :       gel(f, j) = u;
    1746       49140 :       Tr = F2x_div(Tr, u);
    1747       49135 :       av2 = avma;
    1748      183261 :     } else z = gerepileuptoleaf(av2, z);
    1749      232435 :     if (!F2x_degree(Tr)) break;
    1750      220066 :     z = F2xq_sqr(z, Tr);
    1751             :   }
    1752      135240 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1753      135244 :   return gerepilecopy(av, f);
    1754             : }
    1755             : 
    1756             : GEN
    1757           7 : F2x_ddf(GEN T)
    1758             : {
    1759             :   GEN XP;
    1760           7 :   T = F2x_get_red(T);
    1761           7 :   XP = F2x_Frobenius(T);
    1762           7 :   return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
    1763             : }
    1764             : 
    1765             : static GEN
    1766       22165 : F2xq_frobtrace(GEN a, long d, GEN T)
    1767             : {
    1768       22165 :   pari_sp av = avma;
    1769             :   long i;
    1770       22165 :   GEN x = a;
    1771       63818 :   for(i=1; i<d; i++)
    1772             :   {
    1773       41653 :     x = F2x_add(a, F2xq_sqr(x,T));
    1774       41653 :     if (gc_needed(av, 2))
    1775           0 :       x = gerepileuptoleaf(av, x);
    1776             :   }
    1777       22165 :   return x;
    1778             : }
    1779             : 
    1780             : static void
    1781       33112 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1782             : {
    1783       33112 :   long n = F2x_degree(Tp), r = n/d;
    1784             :   GEN T, f, ff;
    1785       33112 :   if (r==1) { gel(V, idx) = Tp; return; }
    1786       11133 :   T = Tp;
    1787       11133 :   XP = F2x_rem(XP, T);
    1788             :   while (1)
    1789       11032 :   {
    1790       22165 :     pari_sp btop = avma;
    1791             :     long df;
    1792       22165 :     GEN g = random_F2x(n, Tp[1]);
    1793       22165 :     GEN t = F2xq_frobtrace(g, d, T);
    1794       22165 :     if (lgpol(t) == 0) continue;
    1795       16570 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1796       16570 :     if (df > 0 && df < n) break;
    1797        5437 :     set_avma(btop);
    1798             :   }
    1799       11133 :   ff = F2x_div(Tp, f);
    1800       11133 :   F2x_edf_simple(f, XP, d, V, idx);
    1801       11133 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1802             : }
    1803             : 
    1804             : static GEN
    1805      291317 : F2x_factor_Shoup(GEN T)
    1806             : {
    1807      291317 :   long i, n, s = 0;
    1808             :   GEN XP, D, V;
    1809             :   pari_timer ti;
    1810      291317 :   n = F2x_degree(T);
    1811      291319 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1812      291319 :   XP = F2x_Frobenius(T);
    1813      291319 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1814      291319 :   D = F2x_ddf_simple(T, XP);
    1815      291321 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1816      977056 :   for (i = 1; i <= n; i++)
    1817      685751 :     s += F2x_degree(gel(D,i))/i;
    1818      291305 :   V = cgetg(s+1, t_COL);
    1819      977107 :   for (i = 1, s = 1; i <= n; i++)
    1820             :   {
    1821      685742 :     GEN Di = gel(D,i);
    1822      685742 :     long ni = F2x_degree(Di), ri = ni/i;
    1823      685729 :     if (ni == 0) continue;
    1824      328077 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1825       10789 :     F2x_edf_simple(Di, XP, i, V, s);
    1826       10846 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1827       10846 :     s += ri;
    1828             :   }
    1829      291365 :   return V;
    1830             : }
    1831             : 
    1832             : static GEN
    1833      246387 : F2x_factor_Cantor(GEN T)
    1834             : {
    1835      246387 :   GEN E, F, V = F2x_factor_squarefree(T);
    1836      246390 :   long i, j, l = lg(V);
    1837      246390 :   E = cgetg(l, t_VEC);
    1838      246388 :   F = cgetg(l, t_VEC);
    1839      881271 :   for (i=1, j=1; i < l; i++)
    1840      634880 :     if (F2x_degree(gel(V,i)))
    1841             :     {
    1842      291318 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1843      291321 :       gel(F, j) = Fj;
    1844      291321 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1845      291320 :       j++;
    1846             :     }
    1847      246391 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1848             : }
    1849             : 
    1850             : #if 0
    1851             : static GEN
    1852             : F2x_simplefact_Shoup(GEN T)
    1853             : {
    1854             :   long i, n, s = 0, j = 1, k;
    1855             :   GEN XP, D, V;
    1856             :   pari_timer ti;
    1857             :   n = F2x_degree(T);
    1858             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1859             :   XP = F2x_Frobenius(T);
    1860             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1861             :   D = F2x_ddf_simple(T, XP);
    1862             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1863             :   for (i = 1; i <= n; i++)
    1864             :     s += F2x_degree(gel(D,i))/i;
    1865             :   V = cgetg(s+1, t_VECSMALL);
    1866             :   for (i = 1; i <= n; i++)
    1867             :   {
    1868             :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1869             :     if (ni == 0) continue;
    1870             :     for (k = 1; k <= ri; k++)
    1871             :       V[j++] = i;
    1872             :   }
    1873             :   return V;
    1874             : }
    1875             : static GEN
    1876             : F2x_simplefact_Cantor(GEN T)
    1877             : {
    1878             :   GEN E, F, V = F2x_factor_squarefree(T);
    1879             :   long i, j, l = lg(V);
    1880             :   F = cgetg(l, t_VEC);
    1881             :   E = cgetg(l, t_VEC);
    1882             :   for (i=1, j=1; i < l; i++)
    1883             :     if (F2x_degree(gel(V,i)))
    1884             :     {
    1885             :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1886             :       gel(F, j) = Fj;
    1887             :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1888             :       j++;
    1889             :     }
    1890             :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1891             : }
    1892             : static int
    1893             : F2x_isirred_Cantor(GEN T)
    1894             : {
    1895             :   pari_sp av = avma;
    1896             :   pari_timer ti;
    1897             :   long n;
    1898             :   GEN dT = F2x_deriv(T);
    1899             :   GEN XP, D;
    1900             :   if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
    1901             :   n = F2x_degree(T);
    1902             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1903             :   XP = F2x_Frobenius(T);
    1904             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1905             :   D = F2x_ddf_simple(T, XP);
    1906             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1907             :   return gc_bool(av, F2x_degree(gel(D,n)) == n);
    1908             : }
    1909             : #endif
    1910             : 
    1911             : /* driver for Cantor factorization, assume deg f > 2; not competitive for
    1912             :  * flag != 0, or as deg f increases */
    1913             : static GEN
    1914      246387 : F2x_Cantor_i(GEN f, long flag)
    1915             : {
    1916             :   switch(flag)
    1917             :   {
    1918      246387 :     default: return F2x_factor_Cantor(f);
    1919             : #if 0
    1920             :     case 1: return F2x_simplefact_Cantor(f);
    1921             :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1922             : #endif
    1923             :   }
    1924             : }
    1925             : static GEN
    1926     1467174 : F2x_factor_i(GEN f, long flag)
    1927             : {
    1928     1467174 :   long d = F2x_degree(f);
    1929     1467205 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1930      257088 :   return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
    1931      538772 :                                : F2x_Berlekamp_i(f, flag);
    1932             : }
    1933             : 
    1934             : GEN
    1935           0 : F2x_degfact(GEN f)
    1936             : {
    1937           0 :   pari_sp av = avma;
    1938           0 :   GEN z = F2x_factor_i(f, 1);
    1939           0 :   return gerepilecopy(av, z);
    1940             : }
    1941             : 
    1942             : int
    1943         238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
    1944             : 
    1945             : /* Adapted from Shoup NTL */
    1946             : GEN
    1947     7040909 : Flx_factor_squarefree(GEN f, ulong p)
    1948             : {
    1949     7040909 :   long i, q, n = degpol(f);
    1950     7040330 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1951     7043777 :   for(q = 1;;q *= p)
    1952       86646 :   {
    1953     7130423 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1954     7124275 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1955      565534 :     t = Flx_div(f, r, p);
    1956      565531 :     if (degpol(t) > 0)
    1957             :     {
    1958             :       long j;
    1959      485820 :       for(j = 1;;j++)
    1960             :       {
    1961     1195589 :         v = Flx_gcd(r, t, p);
    1962     1195572 :         tv = Flx_div(t, v, p);
    1963     1195566 :         if (degpol(tv) > 0)
    1964      728406 :           gel(u, j*q) = Flx_normalize(tv, p);
    1965     1195560 :         if (degpol(v) <= 0) break;
    1966      709744 :         r = Flx_div(r, v, p);
    1967      709769 :         t = v;
    1968             :       }
    1969      485821 :       if (degpol(r) == 0) break;
    1970             :     }
    1971       86644 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1972             :   }
    1973    30217568 :   for (i = n; i; i--)
    1974    30217862 :     if (degpol(gel(u,i))) break;
    1975     7038604 :   setlg(u,i+1); return u;
    1976             : }
    1977             : 
    1978             : long
    1979        3451 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1980             : {
    1981        3451 :   pari_sp av = avma;
    1982             :   ulong lc;
    1983             :   GEN F;
    1984        3451 :   long i, n = degpol(f), v = f[1], l;
    1985        3451 :   if (n % k) return 0;
    1986        3451 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    1987        3451 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    1988        3451 :   F = Flx_factor_squarefree(f, p); l = lg(F)-1;
    1989       38325 :   for (i = 1; i <= l; i++)
    1990       34895 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1991        3430 :   if (pt_r)
    1992             :   {
    1993        3430 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    1994       38304 :     for(i = l; i >= 1; i--)
    1995             :     {
    1996       34874 :       if (i%k) continue;
    1997       12453 :       s = Flx_mul(s, gel(F,i), p);
    1998       12453 :       r = Flx_mul(r, s, p);
    1999             :     }
    2000        3430 :     *pt_r = gerepileuptoleaf(av, r);
    2001           0 :   } else set_avma(av);
    2002        3430 :   return 1;
    2003             : }
    2004             : 
    2005             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    2006             : static GEN
    2007     8937695 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
    2008             : {
    2009     8937695 :   pari_sp av = avma;
    2010             :   GEN b, g, h, F, f, Tr, xq;
    2011             :   long i, j, n, v, bo, ro;
    2012             :   long B, l, m;
    2013             :   pari_timer ti;
    2014     8937695 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    2015     8939491 :   if (n == 0) return cgetg(1, t_VEC);
    2016     8894917 :   if (n == 1) return mkvec(get_Flx_mod(T));
    2017     8598992 :   B = n/2;
    2018     8598992 :   l = usqrt(B);
    2019     8598831 :   m = (B+l-1)/l;
    2020     8598831 :   T = Flx_get_red(T, p);
    2021     8596863 :   b = cgetg(l+2, t_VEC);
    2022     8596567 :   gel(b, 1) = polx_Flx(v);
    2023     8597352 :   gel(b, 2) = XP;
    2024     8597352 :   bo = brent_kung_optpow(n, l-1, 1);
    2025     8598788 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    2026     8598788 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2027     8598788 :   if (expu(p) <= ro)
    2028      877553 :     for (i = 3; i <= l+1; i++)
    2029      486713 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    2030             :   else
    2031             :   {
    2032     8208437 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    2033     8207461 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
    2034     8936954 :     for (i = 3; i <= l+1; i++)
    2035      729380 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    2036             :   }
    2037     8598414 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
    2038     8598414 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    2039     8596851 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
    2040     8596851 :   g = cgetg(m+1, t_VEC);
    2041     8597676 :   gel(g, 1) = gel(xq, 2);
    2042    13958064 :   for(i = 2; i <= m; i++)
    2043     5360380 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    2044     8597684 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
    2045     8597684 :   h = cgetg(m+1, t_VEC);
    2046    22557049 :   for (j = 1; j <= m; j++)
    2047             :   {
    2048    13959519 :     pari_sp av = avma;
    2049    13959519 :     GEN gj = gel(g, j);
    2050    13959519 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    2051    17903308 :     for (i = 2; i <= l; i++)
    2052     3953940 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    2053    13949368 :     gel(h, j) = gerepileupto(av, e);
    2054             :   }
    2055     8597530 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
    2056     8597530 :   Tr = get_Flx_mod(T);
    2057     8598622 :   F = cgetg(m+1, t_VEC);
    2058    22552908 :   for (j = 1; j <= m; j++)
    2059             :   {
    2060    13954287 :     GEN u = Flx_gcd(Tr, gel(h, j), p);
    2061    13956354 :     if (degpol(u))
    2062             :     {
    2063     6359042 :       u = Flx_normalize(u, p);
    2064     6358878 :       Tr = Flx_div(Tr, u, p);
    2065             :     }
    2066    13953111 :     gel(F, j) = u;
    2067             :   }
    2068     8598621 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
    2069     8598621 :   f = const_vec(n, pol1_Flx(v));
    2070    22559273 :   for (j = 1; j <= m; j++)
    2071             :   {
    2072    13960855 :     GEN e = gel(F, j);
    2073    14845075 :     for (i=l-1; i >= 0; i--)
    2074             :     {
    2075    14844367 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    2076    14841411 :       if (degpol(u))
    2077             :       {
    2078     6455262 :         gel(f, l*j-i) = u;
    2079     6455262 :         e = Flx_div(e, u, p);
    2080             :       }
    2081    14838716 :       if (!degpol(e)) break;
    2082             :     }
    2083             :   }
    2084     8598418 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
    2085     8598418 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    2086     8598798 :   return gerepilecopy(av, f);
    2087             : }
    2088             : 
    2089             : static void
    2090      441517 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2091             : {
    2092      441517 :   long n = degpol(Tp), r = n/d;
    2093             :   GEN T, f, ff;
    2094             :   ulong p2;
    2095      441517 :   if (r==1) { gel(V, idx) = Tp; return; }
    2096      190217 :   p2 = p>>1;
    2097      190217 :   T = Flx_get_red(Tp, p);
    2098      190219 :   XP = Flx_rem(XP, T, p);
    2099             :   while (1)
    2100       20691 :   {
    2101      210911 :     pari_sp btop = avma;
    2102             :     long i;
    2103      210911 :     GEN g = random_Flx(n, Tp[1], p);
    2104      210911 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2105      210911 :     if (lgpol(t) == 0) continue;
    2106      458467 :     for(i=1; i<=10; i++)
    2107             :     {
    2108      442869 :       pari_sp btop2 = avma;
    2109      442869 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    2110      442854 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    2111      442869 :       if (degpol(f) > 0 && degpol(f) < n) break;
    2112      252649 :       set_avma(btop2);
    2113             :     }
    2114      205817 :     if (degpol(f) > 0 && degpol(f) < n) break;
    2115       15597 :     set_avma(btop);
    2116             :   }
    2117      190219 :   f = Flx_normalize(f, p);
    2118      190219 :   ff = Flx_div(Tp, f ,p);
    2119      190219 :   Flx_edf_simple(f, XP, d, p, V, idx);
    2120      190220 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    2121             : }
    2122             : static void
    2123             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    2124             : 
    2125             : static void
    2126     1229908 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    2127             : {
    2128             :   pari_sp av;
    2129     1229908 :   GEN Tp = get_Flx_mod(T);
    2130     1229909 :   long n = degpol(hp), vT = Tp[1];
    2131             :   GEN u1, u2, f1, f2;
    2132     1229915 :   ulong p2 = p>>1;
    2133             :   GEN R, h;
    2134     1229915 :   h = Flx_get_red(hp, p);
    2135     1229896 :   t = Flx_rem(t, T, p);
    2136     1229865 :   av = avma;
    2137             :   do
    2138             :   {
    2139     2052614 :     set_avma(av);
    2140     2052652 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    2141     2052468 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    2142     2052615 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2143     1229884 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    2144     1229882 :   f1 = Flx_normalize(f1, p);
    2145     1229894 :   u2 = Flx_div(hp, u1, p);
    2146     1229889 :   f2 = Flx_div(Tp, f1, p);
    2147     1229915 :   if (degpol(u1)==1)
    2148             :   {
    2149      926432 :     if (degpol(f1)==d)
    2150      912086 :       gel(V, idx) = f1;
    2151             :     else
    2152       14333 :       Flx_edf(f1, XP, d, p, V, idx);
    2153             :   }
    2154             :   else
    2155      303528 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    2156     1229949 :   idx += degpol(f1)/d;
    2157     1229923 :   if (degpol(u2)==1)
    2158             :   {
    2159      921165 :     if (degpol(f2)==d)
    2160      907528 :       gel(V, idx) = f2;
    2161             :     else
    2162       13636 :       Flx_edf(f2, XP, d, p, V, idx);
    2163             :   }
    2164             :   else
    2165      308794 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    2166     1229969 : }
    2167             : 
    2168             : static void
    2169      617660 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2170             : {
    2171      617660 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2172             :   GEN T, h, t;
    2173             :   pari_timer ti;
    2174      617660 :   if (r==1) { gel(V, idx) = Tp; return; }
    2175      617660 :   T = Flx_get_red(Tp, p);
    2176      617658 :   XP = Flx_rem(XP, T, p);
    2177      617659 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2178             :   do
    2179             :   {
    2180      635838 :     GEN g = random_Flx(n, vT, p);
    2181      635835 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2182      635854 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2183      635854 :     h = Flxq_minpoly(t, T, p);
    2184      635827 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2185      635827 :   } while (degpol(h) <= 1);
    2186      617649 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    2187             : }
    2188             : 
    2189             : static GEN
    2190     1410116 : Flx_factor_Shoup(GEN T, ulong p)
    2191             : {
    2192     1410116 :   long i, n, s = 0;
    2193             :   GEN XP, D, V;
    2194     1410116 :   long e = expu(p);
    2195             :   pari_timer ti;
    2196     1410114 :   n = get_Flx_degree(T);
    2197     1410114 :   T = Flx_get_red(T, p);
    2198     1410105 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2199     1410105 :   XP = Flx_Frobenius(T, p);
    2200     1410095 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2201     1410095 :   D = Flx_ddf_Shoup(T, XP, p);
    2202     1410145 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2203     1410145 :   s = ddf_to_nbfact(D);
    2204     1410130 :   V = cgetg(s+1, t_COL);
    2205     8191211 :   for (i = 1, s = 1; i <= n; i++)
    2206             :   {
    2207     6781088 :     GEN Di = gel(D,i);
    2208     6781088 :     long ni = degpol(Di), ri = ni/i;
    2209     6781084 :     if (ni == 0) continue;
    2210     1863281 :     Di = Flx_normalize(Di, p);
    2211     1863320 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2212      650754 :     if (ri <= e*expu(e))
    2213      589684 :       Flx_edf(Di, XP, i, p, V, s);
    2214             :     else
    2215       61078 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2216      650760 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2217      650760 :     s += ri;
    2218             :   }
    2219     1410123 :   return V;
    2220             : }
    2221             : 
    2222             : static GEN
    2223     1340418 : Flx_factor_Cantor(GEN T, ulong p)
    2224             : {
    2225     1340418 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2226     1340418 :   long i, j, l = lg(V);
    2227     1340418 :   F = cgetg(l, t_VEC);
    2228     1340416 :   E = cgetg(l, t_VEC);
    2229     3132598 :   for (i=1, j=1; i < l; i++)
    2230     1792177 :     if (degpol(gel(V,i)))
    2231             :     {
    2232     1410116 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2233     1410121 :       gel(F, j) = Fj;
    2234     1410121 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2235     1410117 :       j++;
    2236             :     }
    2237     1340421 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2238             : }
    2239             : 
    2240             : GEN
    2241           0 : Flx_ddf(GEN T, ulong p)
    2242             : {
    2243             :   GEN XP;
    2244           0 :   T = Flx_get_red(T, p);
    2245           0 :   XP = Flx_Frobenius(T, p);
    2246           0 :   return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
    2247             : }
    2248             : 
    2249             : static GEN
    2250     5367384 : Flx_simplefact_Cantor(GEN T, ulong p)
    2251             : {
    2252             :   GEN V;
    2253             :   long i, l;
    2254     5367384 :   T = Flx_get_red(T, p);
    2255     5365297 :   V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
    2256    10818934 :   for (i=1; i < l; i++)
    2257     5448370 :     gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
    2258     5370564 :   return vddf_to_simplefact(V, get_Flx_degree(T));
    2259             : }
    2260             : 
    2261             : static int
    2262        1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2263             : {
    2264        1078 :   pari_sp av = avma;
    2265             :   pari_timer ti;
    2266             :   long n;
    2267        1078 :   GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
    2268        1078 :   if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    2269         791 :   n = get_Flx_degree(T);
    2270         791 :   T = Flx_get_red(Tp, p);
    2271         791 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2272         791 :   XP = Flx_Frobenius(T, p);
    2273         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2274         791 :   D = Flx_ddf_Shoup(T, XP, p);
    2275         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2276         791 :   return gc_bool(av, degpol(gel(D,n)) == n);
    2277             : }
    2278             : 
    2279             : /* f monic */
    2280             : static GEN
    2281    10493415 : Flx_factor_i(GEN f, ulong pp, long flag)
    2282             : {
    2283             :   long d;
    2284    10493415 :   if (pp==2) { /*We need to handle 2 specially */
    2285       73672 :     GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
    2286       73671 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2287       73671 :     return F;
    2288             :   }
    2289    10419743 :   d = degpol(f);
    2290    10420343 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2291     6708352 :   switch(flag)
    2292             :   {
    2293     1340418 :     default: return Flx_factor_Cantor(f, pp);
    2294     5366856 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2295        1078 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2296             :   }
    2297             : }
    2298             : 
    2299             : GEN
    2300     7838588 : Flx_degfact(GEN f, ulong p)
    2301             : {
    2302     7838588 :   pari_sp av = avma;
    2303     7838588 :   GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
    2304     7839084 :   return gerepilecopy(av, z);
    2305             : }
    2306             : 
    2307             : /* T must be squarefree mod p*/
    2308             : GEN
    2309     1444644 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2310             : {
    2311             :   GEN XP, D;
    2312             :   pari_timer ti;
    2313     1444644 :   long i, s, n = get_Flx_degree(T);
    2314     1444592 :   GEN V = const_vecsmall(n, 0);
    2315     1444577 :   pari_sp av = avma;
    2316     1444577 :   T = Flx_get_red(T, p);
    2317     1444545 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2318     1444545 :   XP = Flx_Frobenius(T, p);
    2319     1444569 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2320     1444569 :   D = Flx_ddf_Shoup(T, XP, p);
    2321     1444922 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2322     7424834 :   for (i = 1, s = 0; i <= n; i++)
    2323             :   {
    2324     5979967 :     V[i] = degpol(gel(D,i))/i;
    2325     5979913 :     s += V[i];
    2326             :   }
    2327     1444867 :   *nb = s;
    2328     1444867 :   set_avma(av); return V;
    2329             : }
    2330             : 
    2331             : long
    2332      635405 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2333             : {
    2334      635405 :   pari_sp av = avma;
    2335      635405 :   long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
    2336      635412 :   return gc_long(av,s);
    2337             : }
    2338             : 
    2339             : /* T must be squarefree mod p*/
    2340             : long
    2341      635402 : Flx_nbfact(GEN T, ulong p)
    2342             : {
    2343      635402 :   pari_sp av = avma;
    2344      635402 :   GEN XP = Flx_Frobenius(T, p);
    2345      635404 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2346      635412 :   return gc_long(av,n);
    2347             : }
    2348             : 
    2349             : int
    2350        1057 : Flx_is_irred(GEN f, ulong p)
    2351             : {
    2352        1057 :   pari_sp av = avma;
    2353        1057 :   f = Flx_normalize(f,p);
    2354        1057 :   return gc_bool(av, !!Flx_factor_i(f,p,2));
    2355             : }
    2356             : 
    2357             : /* Use this function when you think f is reducible, and that there are lots of
    2358             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2359             : int
    2360         112 : FpX_is_irred(GEN f, GEN p)
    2361             : {
    2362         112 :   pari_sp av = avma;
    2363             :   int z;
    2364         112 :   switch(ZX_factmod_init(&f,p))
    2365             :   {
    2366          28 :     case 0:  z = !!F2x_factor_i(f,2); break;
    2367          77 :     case 1:  z = !!Flx_factor_i(f,p[2],2); break;
    2368           7 :     default: z = !!FpX_factor_i(f,p,2); break;
    2369             :   }
    2370         112 :   return gc_bool(av,z);
    2371             : }
    2372             : GEN
    2373       63287 : FpX_degfact(GEN f, GEN p) {
    2374       63287 :   pari_sp av = avma;
    2375             :   GEN F;
    2376       63287 :   switch(ZX_factmod_init(&f,p))
    2377             :   {
    2378         567 :     case 0:  F = F2x_factor_i(f,1); break;
    2379       62692 :     case 1:  F = Flx_factor_i(f,p[2],1); break;
    2380          28 :     default: F = FpX_factor_i(f,p,1); break;
    2381             :   }
    2382       63287 :   return gerepilecopy(av, F);
    2383             : }
    2384             : 
    2385             : #if 0
    2386             : /* set x <-- x + c*y mod p */
    2387             : /* x is not required to be normalized.*/
    2388             : static void
    2389             : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2390             : {
    2391             :   long i, lx, ly;
    2392             :   ulong *x=(ulong *)gx;
    2393             :   ulong *y=(ulong *)gy;
    2394             :   if (!c) return;
    2395             :   lx = lg(gx);
    2396             :   ly = lg(gy);
    2397             :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2398             :   if (SMALL_ULONG(p))
    2399             :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2400             :   else
    2401             :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2402             : }
    2403             : #endif
    2404             : 
    2405             : GEN
    2406     3395988 : FpX_factor(GEN f, GEN p)
    2407             : {
    2408     3395988 :   pari_sp av = avma;
    2409             :   GEN F;
    2410     3395988 :   switch(ZX_factmod_init(&f, p))
    2411             :   {
    2412     1377389 :     case 0:  F = F2x_factor_i(f,0);
    2413     1377428 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    2414     2017015 :     case 1:  F = Flx_factor_i(f,p[2],0);
    2415     2017021 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    2416        1616 :     default: F = FpX_factor_i(f,p,0); break;
    2417             :   }
    2418     3396135 :   return gerepilecopy(av, F);
    2419             : }
    2420             : 
    2421             : GEN
    2422      574920 : Flx_factor(GEN f, ulong p)
    2423             : {
    2424      574920 :   pari_sp av = avma;
    2425      574920 :   return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
    2426             : }
    2427             : GEN
    2428       15285 : F2x_factor(GEN f)
    2429             : {
    2430       15285 :   pari_sp av = avma;
    2431       15285 :   return gerepilecopy(av, F2x_factor_i(f,0));
    2432             : }

Generated by: LCOV version 1.13