Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpX_factor.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23219-221e4f55f) Lines: 1243 1341 92.7 %
Date: 2018-11-14 05:40:44 Functions: 110 118 93.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /***********************************************************************/
      18             : /**                                                                   **/
      19             : /**               Factorisation over finite field                     **/
      20             : /**                                                                   **/
      21             : /***********************************************************************/
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : /* Replace F by a monic normalized FpX having the same factors;
      29             :  * assume p prime and *F a ZX */
      30             : static int
      31      921876 : ZX_factmod_init(GEN *F, GEN p)
      32             : {
      33      921876 :   if (lgefint(p) == 3)
      34             :   {
      35      920491 :     ulong pp = p[2];
      36      920491 :     if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
      37      762371 :     *F = ZX_to_Flx(*F, pp);
      38      762378 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      39      762374 :     return 1;
      40             :   }
      41        1385 :   *F = FpX_red(*F, p);
      42        1385 :   if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      43        1385 :   return 2;
      44             : }
      45             : static void
      46       83672 : ZX_rootmod_init(GEN *F, GEN p)
      47             : {
      48       83672 :   if (lgefint(p) == 3)
      49             :   {
      50       75700 :     ulong pp = p[2];
      51       75700 :     *F = ZX_to_Flx(*F, pp);
      52       75700 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      53             :   }
      54             :   else
      55             :   {
      56        7972 :     *F = FpX_red(*F, p);
      57        7972 :     if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      58             :   }
      59       83672 : }
      60             : 
      61             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      62             : static GEN
      63          42 : all_roots_mod_p(ulong p, int not_0)
      64             : {
      65             :   GEN r;
      66             :   ulong i;
      67          42 :   if (not_0) {
      68          28 :     r = cgetg(p, t_VECSMALL);
      69          28 :     for (i = 1; i < p; i++) r[i] = i;
      70             :   } else {
      71          14 :     r = cgetg(p+1, t_VECSMALL);
      72          14 :     for (i = 0; i < p; i++) r[i+1] = i;
      73             :   }
      74          42 :   return r;
      75             : }
      76             : 
      77             : /* X^n - 1 */
      78             : static GEN
      79         126 : Flx_Xnm1(long sv, long n, ulong p)
      80             : {
      81         126 :   GEN t = cgetg(n+3, t_VECSMALL);
      82             :   long i;
      83         126 :   t[1] = sv;
      84         126 :   t[2] = p - 1;
      85         126 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      86         126 :   t[i] = 1; return t;
      87             : }
      88             : /* X^n + 1 */
      89             : static GEN
      90         140 : Flx_Xn1(long sv, long n, ulong p)
      91             : {
      92         140 :   GEN t = cgetg(n+3, t_VECSMALL);
      93             :   long i;
      94             :   (void) p;
      95         140 :   t[1] = sv;
      96         140 :   t[2] = 1;
      97         140 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      98         140 :   t[i] = 1; return t;
      99             : }
     100             : 
     101             : static GEN
     102        9049 : Flx_root_mod_2(GEN f)
     103             : {
     104        9049 :   int z1, z0 = !(f[2] & 1);
     105             :   long i,n;
     106             :   GEN y;
     107             : 
     108        9049 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     109        9049 :   z1 = n & 1;
     110        9049 :   y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
     111        9049 :   if (z0) y[i++] = 0;
     112        9049 :   if (z1) y[i  ] = 1;
     113        9049 :   return y;
     114             : }
     115             : static ulong
     116          56 : Flx_oneroot_mod_2(GEN f)
     117             : {
     118             :   long i,n;
     119          56 :   if (!(f[2] & 1)) return 0;
     120          56 :   for (i=2, n=1; i < lg(f); i++) n += f[i];
     121          56 :   if (n & 1) return 1;
     122          28 :   return 2;
     123             : }
     124             : 
     125             : static GEN FpX_roots_i(GEN f, GEN p);
     126             : static GEN Flx_roots_i(GEN f, ulong p);
     127             : 
     128             : static int
     129     3266337 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     130             : 
     131             : /* Generic driver to computes the roots of f modulo pp, using 'Roots' when
     132             :  * pp is a small prime.
     133             :  * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
     134             :  * assume that f is an FpX, pp a prime and return t_INTs */
     135             : static GEN
     136       71359 : rootmod_aux(GEN f, GEN pp)
     137             : {
     138             :   GEN y;
     139       71359 :   switch(lg(f))
     140             :   {
     141          14 :     case 2: pari_err_ROOTS0("rootmod");
     142          49 :     case 3: return cgetg(1,t_COL);
     143             :   }
     144       71296 :   if (typ(f) == t_VECSMALL)
     145             :   {
     146       68186 :     ulong p = pp[2];
     147       68186 :     if (p == 2)
     148        9049 :       y = Flx_root_mod_2(f);
     149             :     else
     150             :     {
     151       59137 :       if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
     152       59137 :       y = Flx_roots_i(f, p);
     153             :     }
     154       68179 :     y = Flc_to_ZC(y);
     155             :   }
     156             :   else
     157        3110 :     y = FpX_roots_i(f, pp);
     158       71282 :   return y;
     159             : }
     160             : /* assume that f is a ZX and p a prime */
     161             : GEN
     162       71359 : FpX_roots(GEN f, GEN p)
     163             : {
     164       71359 :   pari_sp av = avma;
     165       71359 :   GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
     166       71331 :   return gerepileupto(av, y);
     167             : }
     168             : 
     169             : /* assume x reduced mod p > 2, monic. */
     170             : static int
     171          21 : FpX_quad_factortype(GEN x, GEN p)
     172             : {
     173          21 :   GEN b = gel(x,3), c = gel(x,2);
     174          21 :   GEN D = subii(sqri(b), shifti(c,2));
     175          21 :   return kronecker(D,p);
     176             : }
     177             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     178             : static GEN
     179        7532 : FpX_quad_root(GEN x, GEN p, int unknown)
     180             : {
     181        7532 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     182             : 
     183        7532 :   if (absequaliu(p, 2)) {
     184           0 :     if (!signe(b)) return c;
     185           0 :     return signe(c)? NULL: gen_1;
     186             :   }
     187        7532 :   D = subii(sqri(b), shifti(c,2));
     188        7532 :   D = remii(D,p);
     189        7532 :   if (unknown && kronecker(D,p) == -1) return NULL;
     190             : 
     191        7019 :   s = Fp_sqrt(D,p);
     192             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     193        7019 :   if (!s) return NULL;
     194        7011 :   return Fp_halve(Fp_sub(s,b, p), p);
     195             : }
     196             : static GEN
     197        3146 : FpX_otherroot(GEN x, GEN r, GEN p)
     198        3146 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     199             : 
     200             : /* disc(x^2+bx+c) = b^2 - 4c */
     201             : static ulong
     202    21564138 : Fl_disc_bc(ulong b, ulong c, ulong p)
     203    21564138 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     204             : /* p > 2 */
     205             : static ulong
     206    20487670 : Flx_quad_root(GEN x, ulong p, int unknown)
     207             : {
     208    20487670 :   ulong s, b = x[3], c = x[2];
     209    20487670 :   ulong D = Fl_disc_bc(b, c, p);
     210    20477001 :   if (unknown && krouu(D,p) == -1) return p;
     211    13693153 :   s = Fl_sqrt(D,p);
     212    13725522 :   if (s==~0UL) return p;
     213    13725509 :   return Fl_halve(Fl_sub(s,b, p), p);
     214             : }
     215             : static ulong
     216    12283012 : Flx_otherroot(GEN x, ulong r, ulong p)
     217    12283012 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     218             : 
     219             : 
     220             : /* 'todo' contains the list of factors to be split.
     221             :  * 'done' the list of finished factors, no longer touched */
     222             : struct split_t { GEN todo, done; };
     223             : static void
     224     4932663 : split_init(struct split_t *S, long max)
     225             : {
     226     4932663 :   S->todo = vectrunc_init(max);
     227     4932569 :   S->done = vectrunc_init(max);
     228     4932389 : }
     229             : #if 0
     230             : /* move todo[i] to done */
     231             : static void
     232             : split_convert(struct split_t *S, long i)
     233             : {
     234             :   long n = lg(S->todo)-1;
     235             :   vectrunc_append(S->done, gel(S->todo,i));
     236             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     237             :   setlg(S->todo, n);
     238             : }
     239             : #endif
     240             : /* append t to todo */
     241             : static void
     242     5176958 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     243             : /* delete todo[i], add t to done */
     244             : static void
     245     5177006 : split_moveto_done(struct split_t *S, long i, GEN t)
     246             : {
     247     5177006 :   long n = lg(S->todo)-1;
     248     5177006 :   vectrunc_append(S->done, t);
     249     5177335 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     250     5177335 :   setlg(S->todo, n);
     251             : 
     252     5177391 : }
     253             : /* append t to done */
     254             : static void
     255      391977 : split_add_done(struct split_t *S, GEN t)
     256      391977 : { vectrunc_append(S->done, t); }
     257             : /* split todo[i] into a and b */
     258             : static void
     259      337395 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     260             : {
     261      337395 :   gel(S->todo, i) = a;
     262      337395 :   split_add(S, b);
     263      337399 : }
     264             : /* split todo[i] into a and b, moved to done */
     265             : static void
     266      372603 : split_done(struct split_t *S, long i, GEN a, GEN b)
     267             : {
     268      372603 :   split_moveto_done(S, i, a);
     269      372601 :   split_add_done(S, b);
     270      372601 : }
     271             : 
     272             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
     273             : static GEN
     274        3110 : FpX_roots_i(GEN f, GEN p)
     275             : {
     276             :   GEN pol, pol0, a, q;
     277             :   struct split_t S;
     278             : 
     279        3110 :   split_init(&S, lg(f)-1);
     280        3110 :   settyp(S.done, t_COL);
     281        3110 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     282        3110 :   switch(degpol(f))
     283             :   {
     284           7 :     case 0: return ZC_copy(S.done);
     285          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     286             :     case 2: {
     287        1743 :       GEN s, r = FpX_quad_root(f, p, 1);
     288        1743 :       if (r) {
     289        1743 :         split_add_done(&S, r);
     290        1743 :         s = FpX_otherroot(f,r, p);
     291             :         /* f not known to be square free yet */
     292        1743 :         if (!equalii(r, s)) split_add_done(&S, s);
     293             :       }
     294        1743 :       return sort(S.done);
     295             :     }
     296             :   }
     297             : 
     298        1346 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     299        1346 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     300        1346 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     301        1346 :   a = FpX_gcd(f,a, p);
     302        1346 :   if (!degpol(a)) return ZC_copy(S.done);
     303        1346 :   split_add(&S, FpX_normalize(a,p));
     304             : 
     305        1346 :   q = shifti(p,-1);
     306        1346 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     307        1346 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     308        2819 :   for (pol0[2] = 1;; pol0[2]++)
     309        1473 :   {
     310        2819 :     long j, l = lg(S.todo);
     311        2819 :     if (l == 1) return sort(S.done);
     312        1480 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     313        3074 :     for (j = 1; j < l; j++)
     314             :     {
     315        1601 :       GEN b, r, s, c = gel(S.todo,j);
     316        1601 :       switch(degpol(c))
     317             :       { /* convert linear and quadratics to roots, try to split the rest */
     318             :         case 1:
     319         195 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     320         195 :           j--; l--; break;
     321             :         case 2:
     322        1272 :           r = FpX_quad_root(c, p, 0);
     323        1272 :           if (!r) pari_err_PRIME("polrootsmod",p);
     324        1265 :           s = FpX_otherroot(c,r, p);
     325        1265 :           split_done(&S, j, r, s);
     326        1265 :           j--; l--; break;
     327             :         default:
     328         134 :           b = FpXQ_pow(pol,q, c,p);
     329         134 :           if (degpol(b) <= 0) continue;
     330         121 :           b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
     331         121 :           if (!degpol(b)) continue;
     332         121 :           b = FpX_normalize(b, p);
     333         121 :           c = FpX_div(c,b, p);
     334         121 :           split_todo(&S, j, b, c);
     335             :       }
     336             :     }
     337             :   }
     338             : }
     339             : 
     340             : /* Assume f is normalized */
     341             : static ulong
     342      174062 : Flx_cubic_root(GEN ff, ulong p)
     343             : {
     344      174062 :   GEN f = Flx_normalize(ff,p);
     345      174064 :   ulong pi = get_Fl_red(p);
     346      174064 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     347      174064 :   ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
     348      174066 :   ulong A = Fl_sub(b, Fl_triple(t2, p), p);
     349      174066 :   ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
     350      174065 :   ulong A3 =  Fl_mul_pre(A, p3, p, pi);
     351      174065 :   ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
     352      174065 :   ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
     353      174064 :   ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
     354      174066 :   ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     355      174066 :   if (s!=~0UL)
     356             :   {
     357      109744 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     358      109744 :     if (p%3==2) /* 1 solutions */
     359       19354 :       vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
     360             :     else
     361             :     {
     362       90390 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     363       90387 :       if (vS1==~0UL) return p; /*0 solutions*/
     364             :       /*3 solutions*/
     365             :     }
     366       79209 :     vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
     367       79211 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     368             :   }
     369             :   else
     370             :   {
     371       64322 :     pari_sp av = avma;
     372       64322 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
     373       64322 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     374             :     ulong Sa;
     375       64322 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     376       64322 :     Sa = vS1[1];
     377       64322 :     if (p%3==1) /*1 solutions*/
     378             :     {
     379       24051 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     380       24051 :       if (Fa!=P)
     381       15987 :         Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     382             :     }
     383       64322 :     set_avma(av);
     384       64322 :     return Fl_sub(Fl_double(Sa,p),t,p);
     385             :   }
     386             : }
     387             : 
     388             : /* assume p > 2 prime */
     389             : static ulong
     390     3127934 : Flx_oneroot_i(GEN f, ulong p, long fl)
     391             : {
     392             :   GEN pol, a;
     393             :   ulong q;
     394             :   long da;
     395             : 
     396     3127934 :   if (Flx_val(f)) return 0;
     397     3127225 :   switch(degpol(f))
     398             :   {
     399       12187 :     case 1: return Fl_neg(f[2], p);
     400     2680824 :     case 2: return Flx_quad_root(f, p, 1);
     401      158507 :     case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
     402             :   }
     403             : 
     404      275705 :   if (!fl)
     405             :   {
     406      244316 :     a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
     407      244316 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     408      244316 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     409      244316 :     a = Flx_gcd(f,a, p);
     410       31389 :   } else a = f;
     411      275705 :   da = degpol(a);
     412      275704 :   if (!da) return p;
     413      197465 :   a = Flx_normalize(a,p);
     414             : 
     415      197467 :   q = p >> 1;
     416      197467 :   pol = polx_Flx(f[1]);
     417      297309 :   for(pol[2] = 1;; pol[2]++)
     418             :   {
     419      397155 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     420      297313 :     switch(da)
     421             :     {
     422      123317 :       case 1: return Fl_neg(a[2], p);
     423       58622 :       case 2: return Flx_quad_root(a, p, 0);
     424       15561 :       case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
     425             :       default: {
     426       99813 :         GEN b = Flxq_powu(pol,q, a,p);
     427             :         long db;
     428       99849 :         if (degpol(b) <= 0) continue;
     429       94843 :         b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
     430       94836 :         db = degpol(b); if (!db) continue;
     431       94830 :         b = Flx_normalize(b, p);
     432       94842 :         if (db <= (da >> 1)) {
     433       57795 :           a = b;
     434       57795 :           da = db;
     435             :         } else {
     436       37047 :           a = Flx_div(a,b, p);
     437       37043 :           da -= db;
     438             :         }
     439             :       }
     440             :     }
     441             :   }
     442             : }
     443             : 
     444             : /* assume p > 2 prime */
     445             : static GEN
     446        4848 : FpX_oneroot_i(GEN f, GEN p)
     447             : {
     448             :   GEN pol, pol0, a, q;
     449             :   long da;
     450             : 
     451        4848 :   if (ZX_val(f)) return gen_0;
     452        4582 :   switch(degpol(f))
     453             :   {
     454         675 :     case 1: return subii(p, gel(f,2));
     455        3837 :     case 2: return FpX_quad_root(f, p, 1);
     456             :   }
     457             : 
     458          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     459          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     460          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     461          70 :   a = FpX_gcd(f,a, p);
     462          70 :   da = degpol(a);
     463          70 :   if (!da) return NULL;
     464          70 :   a = FpX_normalize(a,p);
     465             : 
     466          70 :   q = shifti(p,-1);
     467          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     468          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     469         224 :   for (pol0[2]=1; ; pol0[2]++)
     470             :   {
     471         378 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     472         224 :     switch(da)
     473             :     {
     474          42 :       case 1: return subii(p, gel(a,2));
     475          28 :       case 2: return FpX_quad_root(a, p, 0);
     476             :       default: {
     477         154 :         GEN b = FpXQ_pow(pol,q, a,p);
     478             :         long db;
     479         154 :         if (degpol(b) <= 0) continue;
     480         147 :         b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
     481         147 :         db = degpol(b); if (!db) continue;
     482         147 :         b = FpX_normalize(b, p);
     483         147 :         if (db <= (da >> 1)) {
     484         105 :           a = b;
     485         105 :           da = db;
     486             :         } else {
     487          42 :           a = FpX_div(a,b, p);
     488          42 :           da -= db;
     489             :         }
     490             :       }
     491             :     }
     492             :   }
     493             : }
     494             : 
     495             : ulong
     496     3088281 : Flx_oneroot(GEN f, ulong p)
     497             : {
     498     3088281 :   pari_sp av = avma;
     499             :   ulong r;
     500     3088281 :   switch(lg(f))
     501             :   {
     502           0 :     case 2: return 0;
     503           0 :     case 3: return p;
     504             :   }
     505     3088281 :   if (p == 2) return Flx_oneroot_mod_2(f);
     506     3088281 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
     507     3088281 :   return gc_ulong(av,r);
     508             : }
     509             : 
     510             : ulong
     511       32248 : Flx_oneroot_split(GEN f, ulong p)
     512             : {
     513       32248 :   pari_sp av = avma;
     514             :   ulong r;
     515       32248 :   switch(lg(f))
     516             :   {
     517           0 :     case 2: return 0;
     518           0 :     case 3: return p;
     519             :   }
     520       32248 :   if (p == 2) return Flx_oneroot_mod_2(f);
     521       32248 :   r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
     522       32282 :   return gc_ulong(av,r);
     523             : }
     524             : 
     525             : /* assume that p is prime */
     526             : GEN
     527       12313 : FpX_oneroot(GEN f, GEN pp) {
     528       12313 :   pari_sp av = avma;
     529       12313 :   ZX_rootmod_init(&f, pp);
     530       12313 :   switch(lg(f))
     531             :   {
     532           0 :     case 2: set_avma(av); return gen_0;
     533           0 :     case 3: return gc_NULL(av);
     534             :   }
     535       12313 :   if (typ(f) == t_VECSMALL)
     536             :   {
     537        7465 :     ulong r, p = pp[2];
     538        7465 :     if (p == 2)
     539          56 :       r = Flx_oneroot_mod_2(f);
     540             :     else
     541        7409 :       r = Flx_oneroot_i(f, p, 0);
     542        7465 :     set_avma(av);
     543        7465 :     return (r == p)? NULL: utoi(r);
     544             :   }
     545        4848 :   f = FpX_oneroot_i(f, pp);
     546        4848 :   if (!f) return gc_NULL(av);
     547        4848 :   return gerepileuptoint(av, f);
     548             : }
     549             : 
     550             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     551             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     552             : /* returned in n                                                          */
     553             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     554             : /* twice as big and decrement until it divides p-1.                       */
     555             : static GEN
     556         217 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     557             : {
     558         217 :    pari_sp ltop = avma;
     559             :    GEN pm, factn, power, base, zeta;
     560             :    long n;
     561             : 
     562         217 :    pm = subis (p, 1ul);
     563         217 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     564         217 :    factn = Z_factor(stoi(n));
     565         217 :    power = diviuexact (pm, n);
     566         217 :    base = gen_1;
     567             :    do {
     568         378 :       base = addis (base, 1l);
     569         378 :       zeta = Fp_pow (base, power, p);
     570             :    }
     571         378 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     572         217 :    *pt_n = n;
     573         217 :    return gerepileuptoint (ltop, zeta);
     574             : }
     575             : 
     576             : GEN
     577         924 : FpX_oneroot_split(GEN fact, GEN p)
     578             : {
     579         924 :   pari_sp av = avma;
     580             :   long n, deg_f, i, dmin;
     581             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     582         924 :   fact = FpX_normalize(fact, p);
     583         924 :   deg_f = degpol(fact);
     584         924 :   if (deg_f<=2) return FpX_oneroot(fact, p);
     585         217 :   minfactor = fact; /* factor of minimal degree found so far */
     586         217 :   dmin = degpol(minfactor);
     587         217 :   prim = good_root_of_unity(p, deg_f, 1, &n);
     588         217 :   expo = diviuexact(subiu(p, 1), n);
     589         217 :   xplusa = pol_x(varn(fact));
     590         217 :   zeta = gen_1;
     591        1071 :   while (dmin != 1)
     592             :   {
     593             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     594             :     /* zeta varies over the roots of unity in F_p                         */
     595         637 :     fact = minfactor; deg_f = dmin;
     596             :     /* update X+a, avoid a=0 */
     597         637 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     598         637 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     599        1260 :     for (i = 0; i < n; i++)
     600             :     {
     601         980 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     602         980 :       long dtmp = degpol(tmp);
     603         980 :       if (dtmp > 0 && dtmp < deg_f)
     604             :       {
     605         427 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     606         427 :         if (dtmp < dmin)
     607             :         {
     608         427 :           minfactor = FpX_normalize (tmp, p);
     609         427 :           dmin = dtmp;
     610         427 :           if (dmin == 1 || dmin <= deg_f / (n / 2) + 1)
     611             :             /* stop early to avoid too many gcds */
     612             :             break;
     613             :         }
     614             :       }
     615         623 :       zeta = Fp_mul (zeta, prim, p);
     616             :     }
     617             :   }
     618         217 :   return gerepileuptoint(av, Fp_neg(gel(minfactor,2), p));
     619             : }
     620             : 
     621             : /*******************************************************************/
     622             : /*                                                                 */
     623             : /*                     FACTORISATION MODULO p                      */
     624             : /*                                                                 */
     625             : /*******************************************************************/
     626             : 
     627             : /* F / E  a vector of vectors of factors / exponents of virtual length l
     628             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
     629             : static GEN
     630      490076 : FE_concat(GEN F, GEN E, long l)
     631             : {
     632      490076 :   setlg(E,l); E = shallowconcat1(E);
     633      490075 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
     634             : }
     635             : 
     636             : static GEN
     637          14 : ddf_to_ddf2_i(GEN V, long fl)
     638             : {
     639             :   GEN F, D;
     640          14 :   long i, j, l = lg(V);
     641          14 :   F = cgetg(l, t_VEC);
     642          14 :   D = cgetg(l, t_VECSMALL);
     643         112 :   for (i = j = 1; i < l; i++)
     644             :   {
     645          98 :     GEN Vi = gel(V,i);
     646          98 :     if ((fl==2 && F2x_degree(Vi) == 0)
     647          84 :       ||(fl==0 && degpol(Vi) == 0)) continue;
     648          35 :     gel(F,j) = Vi;
     649          35 :     uel(D,j) = i; j++;
     650             :   }
     651          14 :   setlg(F,j);
     652          14 :   setlg(D,j); return mkvec2(F,D);
     653             : }
     654             : 
     655             : GEN
     656           7 : ddf_to_ddf2(GEN V)
     657           7 : { return ddf_to_ddf2_i(V, 0); }
     658             : 
     659             : static GEN
     660           7 : F2x_ddf_to_ddf2(GEN V)
     661           7 : { return ddf_to_ddf2_i(V, 2); }
     662             : 
     663             : GEN
     664      366805 : vddf_to_simplefact(GEN V, long d)
     665             : {
     666             :   GEN E, F;
     667      366805 :   long i, j, c, l = lg(V);
     668      366805 :   F = cgetg(d+1, t_VECSMALL);
     669      366805 :   E = cgetg(d+1, t_VECSMALL);
     670      736838 :   for (i = c = 1; i < l; i++)
     671             :   {
     672      370033 :     GEN Vi = gel(V,i);
     673      370033 :     long l = lg(Vi);
     674     2528616 :     for (j = 1; j < l; j++)
     675             :     {
     676     2158583 :       long k, n = degpol(gel(Vi,j)) / j;
     677     2158583 :       for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
     678             :     }
     679             :   }
     680      366805 :   setlg(F,c);
     681      366805 :   setlg(E,c);
     682      366805 :   return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
     683             : }
     684             : 
     685             : /* product of terms of degree 1 in factorization of f */
     686             : GEN
     687      142964 : FpX_split_part(GEN f, GEN p)
     688             : {
     689      142964 :   long n = degpol(f);
     690      142964 :   GEN z, X = pol_x(varn(f));
     691      142964 :   if (n <= 1) return f;
     692      142649 :   f = FpX_red(f, p);
     693      142649 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     694      142649 :   return FpX_gcd(z,f,p);
     695             : }
     696             : 
     697             : /* Compute the number of roots in Fp without counting multiplicity
     698             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     699             : long
     700       99149 : FpX_nbroots(GEN f, GEN p)
     701             : {
     702       99149 :   pari_sp av = avma;
     703       99149 :   GEN z = FpX_split_part(f, p);
     704       99149 :   return gc_long(av, degpol(z));
     705             : }
     706             : 
     707             : int
     708           0 : FpX_is_totally_split(GEN f, GEN p)
     709             : {
     710           0 :   long n=degpol(f);
     711           0 :   pari_sp av = avma;
     712           0 :   if (n <= 1) return 1;
     713           0 :   if (abscmpui(n, p) > 0) return 0;
     714           0 :   f = FpX_red(f, p);
     715           0 :   return gc_bool(av, gequalX(FpX_Frobenius(f, p)));
     716             : }
     717             : 
     718             : long
     719     2408412 : Flx_nbroots(GEN f, ulong p)
     720             : {
     721     2408412 :   long n = degpol(f);
     722     2408411 :   pari_sp av = avma;
     723             :   GEN z;
     724     2408411 :   if (n <= 1) return n;
     725     2406206 :   if (n == 2)
     726             :   {
     727             :     ulong D;
     728       11998 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     729       11116 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     730       11116 :     return 1 + krouu(D,p);
     731             :   }
     732     2394208 :   z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
     733     2394208 :   z = Flx_gcd(z, f, p);
     734     2394214 :   return gc_long(av, degpol(z));
     735             : }
     736             : 
     737             : long
     738        4389 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
     739             : {
     740        4389 :   pari_sp av = avma;
     741             :   GEN X, b, g, xq;
     742             :   long i, j, n, v, B, l, m;
     743             :   pari_timer ti;
     744             :   hashtable h;
     745             : 
     746        4389 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     747        4389 :   X = pol_x(v);
     748        4389 :   if (ZX_equal(X,XP)) return 1;
     749        4389 :   B = n/2;
     750        4389 :   l = usqrt(B);
     751        4389 :   m = (B+l-1)/l;
     752        4389 :   T = FpX_get_red(T, p);
     753        4389 :   hash_init_GEN(&h, l+2, ZX_equal, 1);
     754        4389 :   hash_insert_long(&h, X,  0);
     755        4389 :   hash_insert_long(&h, XP, 1);
     756        4389 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     757        4389 :   b = XP;
     758        4389 :   xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
     759        4389 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
     760       10367 :   for (i = 3; i <= l+1; i++)
     761             :   {
     762        6776 :     b = FpX_FpXQV_eval(b, xq, T, p);
     763        6776 :     if (gequalX(b)) return gc_long(av,i-1);
     764        5978 :     hash_insert_long(&h, b, i-1);
     765             :   }
     766        3591 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
     767        3591 :   g = b;
     768        3591 :   xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
     769        3591 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
     770       12278 :   for(i = 2; i <= m+1; i++)
     771             :   {
     772       10696 :     g = FpX_FpXQV_eval(g, xq, T, p);
     773       10696 :     if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
     774             :   }
     775        1582 :   return gc_long(av,n);
     776             : }
     777             : 
     778             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     779             : static GEN
     780         669 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
     781             : {
     782             :   GEN b, g, h, F, f, Tr, xq;
     783             :   long i, j, n, v, B, l, m;
     784             :   pari_timer ti;
     785             : 
     786         669 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     787         669 :   if (n == 0) return cgetg(1, t_VEC);
     788         669 :   if (n == 1) return mkvec(get_FpX_mod(T));
     789         623 :   B = n/2;
     790         623 :   l = usqrt(B);
     791         623 :   m = (B+l-1)/l;
     792         623 :   T = FpX_get_red(T, p);
     793         623 :   b = cgetg(l+2, t_VEC);
     794         623 :   gel(b, 1) = pol_x(v);
     795         623 :   gel(b, 2) = XP;
     796         623 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     797         623 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     798         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
     799         811 :   for (i = 3; i <= l+1; i++)
     800         188 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     801         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
     802         623 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     803         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
     804         623 :   g = cgetg(m+1, t_VEC);
     805         623 :   gel(g, 1) = gel(xq, 2);
     806         623 :   for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     807         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
     808         623 :   h = cgetg(m+1, t_VEC);
     809        1945 :   for (j = 1; j <= m; j++)
     810             :   {
     811        1322 :     pari_sp av = avma;
     812        1322 :     GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
     813        1322 :     for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
     814        1322 :     gel(h,j) = gerepileupto(av, e);
     815             :   }
     816         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
     817         623 :   Tr = get_FpX_mod(T);
     818         623 :   F = cgetg(m+1, t_VEC);
     819        1945 :   for (j = 1; j <= m; j++)
     820             :   {
     821        1322 :     GEN u = FpX_gcd(Tr, gel(h,j), p);
     822        1322 :     if (degpol(u))
     823             :     {
     824         226 :       u = FpX_normalize(u, p);
     825         226 :       Tr = FpX_div(Tr, u, p);
     826             :     }
     827        1322 :     gel(F,j) = u;
     828             :   }
     829         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
     830         623 :   f = const_vec(n, pol_1(v));
     831        1945 :   for (j = 1; j <= m; j++)
     832             :   {
     833        1322 :     GEN e = gel(F, j);
     834        1357 :     for (i=l-1; i >= 0; i--)
     835             :     {
     836        1357 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     837        1357 :       if (degpol(u))
     838             :       {
     839         233 :         u = FpX_normalize(u, p);
     840         233 :         gel(f, l*j-i) = u;
     841         233 :         e = FpX_div(e, u, p);
     842             :       }
     843        1357 :       if (!degpol(e)) break;
     844             :     }
     845             :   }
     846         623 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
     847         623 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     848         623 :   return f;
     849             : }
     850             : 
     851             : static void
     852           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     853             : {
     854           0 :   long n = degpol(Tp), r = n/d, ct = 0;
     855             :   GEN T, f, ff, p2;
     856           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     857           0 :   p2 = shifti(p,-1);
     858           0 :   T = FpX_get_red(Tp, p);
     859           0 :   XP = FpX_rem(XP, T, p);
     860             :   while (1)
     861           0 :   {
     862           0 :     pari_sp btop = avma;
     863             :     long i;
     864           0 :     GEN g = random_FpX(n, varn(Tp), p);
     865           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     866           0 :     if (signe(t) == 0) continue;
     867           0 :     for(i=1; i<=10; i++)
     868             :     {
     869           0 :       pari_sp btop2 = avma;
     870           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     871           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     872           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     873           0 :       set_avma(btop2);
     874             :     }
     875           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     876           0 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
     877           0 :     set_avma(btop);
     878             :   }
     879           0 :   f = FpX_normalize(f, p);
     880           0 :   ff = FpX_div(Tp, f ,p);
     881           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     882           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     883             : }
     884             : 
     885             : static void
     886         462 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     887             : {
     888             :   pari_sp av;
     889         462 :   GEN Tp = get_FpX_mod(T);
     890         462 :   long n = degpol(hp), vT = varn(Tp), ct = 0;
     891             :   GEN u1, u2, f1, f2, R, h;
     892         462 :   h = FpX_get_red(hp, p);
     893         462 :   t = FpX_rem(t, T, p);
     894         462 :   av = avma;
     895             :   do
     896             :   {
     897         761 :     set_avma(av);
     898         761 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     899         761 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     900         761 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
     901         761 :   } while (degpol(u1)==0 || degpol(u1)==n);
     902         462 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     903         462 :   f1 = FpX_normalize(f1, p);
     904         462 :   u2 = FpX_div(hp, u1, p);
     905         462 :   f2 = FpX_div(Tp, f1, p);
     906         462 :   if (degpol(u1)==1)
     907         336 :     gel(V, idx) = f1;
     908             :   else
     909         126 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     910         462 :   idx += degpol(f1)/d;
     911         462 :   if (degpol(u2)==1)
     912         350 :     gel(V, idx) = f2;
     913             :   else
     914         112 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     915         462 : }
     916             : 
     917             : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
     918             : static void
     919         224 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     920             : {
     921         224 :   long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
     922             :   GEN T, h, t;
     923             :   pari_timer ti;
     924             : 
     925         224 :   T = FpX_get_red(Tp, p);
     926         224 :   XP = FpX_rem(XP, T, p);
     927         224 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     928             :   do
     929             :   {
     930         224 :     GEN g = random_FpX(n, vT, p);
     931         224 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     932         224 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
     933         224 :     h = FpXQ_minpoly(t, T, p);
     934         224 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
     935         224 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
     936         224 :   } while (degpol(h) != r);
     937         224 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
     938         224 : }
     939             : 
     940             : static GEN
     941         655 : FpX_factor_Shoup(GEN T, GEN p)
     942             : {
     943         655 :   long i, n, s = 0;
     944             :   GEN XP, D, V;
     945         655 :   long e = expi(p);
     946             :   pari_timer ti;
     947         655 :   n = get_FpX_degree(T);
     948         655 :   T = FpX_get_red(T, p);
     949         655 :   if (DEBUGLEVEL>=6) timer_start(&ti);
     950         655 :   XP = FpX_Frobenius(T, p);
     951         655 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
     952         655 :   D = FpX_ddf_Shoup(T, XP, p);
     953         655 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
     954         655 :   s = ddf_to_nbfact(D);
     955         655 :   V = cgetg(s+1, t_COL);
     956        5507 :   for (i = 1, s = 1; i <= n; i++)
     957             :   {
     958        4852 :     GEN Di = gel(D,i);
     959        4852 :     long ni = degpol(Di), ri = ni/i;
     960        4852 :     if (ni == 0) continue;
     961         664 :     Di = FpX_normalize(Di, p);
     962         664 :     if (ni == i) { gel(V, s++) = Di; continue; }
     963         224 :     if (ri <= e*expu(e))
     964         224 :       FpX_edf(Di, XP, i, p, V, s);
     965             :     else
     966           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
     967         224 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
     968         224 :     s += ri;
     969             :   }
     970         655 :   return V;
     971             : }
     972             : 
     973             : long
     974      551559 : ddf_to_nbfact(GEN D)
     975             : {
     976      551559 :   long l = lg(D), i, s = 0;
     977      551559 :   for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
     978      551556 :   return s;
     979             : }
     980             : 
     981             : /* Yun algorithm: Assume p > degpol(T) */
     982             : static GEN
     983         672 : FpX_factor_Yun(GEN T, GEN p)
     984             : {
     985         672 :   long n = degpol(T), i = 1;
     986         672 :   GEN a, b, c, d = FpX_deriv(T, p);
     987         672 :   GEN V = cgetg(n+1,t_VEC);
     988         672 :   a = FpX_gcd(T, d, p);
     989         672 :   if (degpol(a) == 0) return mkvec(T);
     990         488 :   b = FpX_div(T, a, p);
     991             :   do
     992             :   {
     993        2432 :     c = FpX_div(d, a, p);
     994        2432 :     d = FpX_sub(c, FpX_deriv(b, p), p);
     995        2432 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
     996        2432 :     gel(V, i++) = a;
     997        2432 :     b = FpX_div(b, a, p);
     998        2432 :   } while (degpol(b));
     999         488 :   setlg(V, i); return V;
    1000             : }
    1001             : GEN
    1002         777 : FpX_factor_squarefree(GEN T, GEN p)
    1003             : {
    1004         777 :   if (lgefint(p)==3)
    1005             :   {
    1006         777 :     ulong pp = (ulong)p[2];
    1007         777 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1008         777 :     return FlxV_to_ZXV(u);
    1009             :   }
    1010           0 :   return FpX_factor_Yun(T, p);
    1011             : }
    1012             : 
    1013             : long
    1014         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1015             : {
    1016         168 :   pari_sp av = avma;
    1017             :   GEN lc, F;
    1018         168 :   long i, l, n = degpol(f), v = varn(f);
    1019         168 :   if (n % k) return 0;
    1020         168 :   if (lgefint(p)==3)
    1021             :   {
    1022         126 :     ulong pp = p[2];
    1023         126 :     GEN fp = ZX_to_Flx(f, pp);
    1024         126 :     if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
    1025         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
    1026         105 :     return 1;
    1027             :   }
    1028          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1029          42 :   if (!lc) { av = avma; return 0; }
    1030          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1031        1491 :   for(i=1; i <= l; i++)
    1032        1456 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1033          35 :   if (pt_r)
    1034             :   {
    1035          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1036        1484 :     for (i=l; i>=1; i--)
    1037             :     {
    1038        1449 :       if (i%k) continue;
    1039         294 :       s = FpX_mul(s, gel(F,i), p);
    1040         294 :       r = FpX_mul(r, s, p);
    1041             :     }
    1042          35 :     *pt_r = gerepileupto(av, r);
    1043           0 :   } else av = avma;
    1044          35 :   return 1;
    1045             : }
    1046             : 
    1047             : static GEN
    1048         616 : FpX_factor_Cantor(GEN T, GEN p)
    1049             : {
    1050         616 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1051         616 :   long i, j, l = lg(V);
    1052         616 :   F = cgetg(l, t_VEC);
    1053         616 :   E = cgetg(l, t_VEC);
    1054        1748 :   for (i=1, j=1; i < l; i++)
    1055        1132 :     if (degpol(gel(V,i)))
    1056             :     {
    1057         655 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1058         655 :       gel(F, j) = Fj;
    1059         655 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1060         655 :       j++;
    1061             :     }
    1062         616 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1063             : }
    1064             : 
    1065             : static GEN
    1066           0 : FpX_ddf_i(GEN T, GEN p)
    1067             : {
    1068             :   GEN XP;
    1069           0 :   T = FpX_get_red(T, p);
    1070           0 :   XP = FpX_Frobenius(T, p);
    1071           0 :   return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
    1072             : }
    1073             : 
    1074             : GEN
    1075           7 : FpX_ddf(GEN f, GEN p)
    1076             : {
    1077           7 :   pari_sp av = avma;
    1078             :   GEN F;
    1079           7 :   switch(ZX_factmod_init(&f, p))
    1080             :   {
    1081           7 :     case 0:  F = F2x_ddf(f);
    1082           7 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    1083           0 :     case 1:  F = Flx_ddf(f,p[2]);
    1084           0 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    1085           0 :     default: F = FpX_ddf_i(f,p); break;
    1086             :   }
    1087           7 :   return gerepilecopy(av, F);
    1088             : }
    1089             : 
    1090             : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
    1091             : static GEN
    1092          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1093             : {
    1094             :   GEN V;
    1095             :   long i, l;
    1096          14 :   if (lgefint(p) == 3)
    1097             :   {
    1098           0 :     ulong pp = p[2];
    1099           0 :     return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
    1100             :   }
    1101          14 :   T = FpX_get_red(T, p);
    1102          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
    1103          28 :   for (i=1; i < l; i++)
    1104          14 :     gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
    1105          14 :   return vddf_to_simplefact(V, get_FpX_degree(T));
    1106             : }
    1107             : 
    1108             : static int
    1109           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1110             : {
    1111           0 :   pari_sp av = avma;
    1112             :   pari_timer ti;
    1113             :   long n;
    1114           0 :   GEN T = get_FpX_mod(Tp);
    1115           0 :   GEN dT = FpX_deriv(T, p);
    1116             :   GEN XP, D;
    1117           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    1118           0 :   n = get_FpX_degree(T);
    1119           0 :   T = FpX_get_red(Tp, p);
    1120           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1121           0 :   XP = FpX_Frobenius(T, p);
    1122           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1123           0 :   D = FpX_ddf_Shoup(T, XP, p);
    1124           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1125           0 :   return gc_bool(av, degpol(gel(D,n)) == n);
    1126             : }
    1127             : 
    1128             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1129             : 
    1130             : /*Assume that p is large and odd*/
    1131             : static GEN
    1132        1385 : FpX_factor_i(GEN f, GEN pp, long flag)
    1133             : {
    1134        1385 :   long d = degpol(f);
    1135        1385 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1136         630 :   switch(flag)
    1137             :   {
    1138         616 :     default: return FpX_factor_Cantor(f, pp);
    1139          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1140           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1141             :   }
    1142             : }
    1143             : 
    1144             : long
    1145           0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1146             : {
    1147           0 :   pari_sp av = avma;
    1148           0 :   long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
    1149           0 :   return gc_long(av,s);
    1150             : }
    1151             : 
    1152             : long
    1153           0 : FpX_nbfact(GEN T, GEN p)
    1154             : {
    1155           0 :   pari_sp av = avma;
    1156           0 :   GEN XP = FpX_Frobenius(T, p);
    1157           0 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1158           0 :   return gc_long(av,n);
    1159             : }
    1160             : 
    1161             : /* p > 2 */
    1162             : static GEN
    1163           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1164             : {
    1165           7 :   switch(d)
    1166             :   {
    1167             :     case -1:
    1168           0 :     case 0: return NULL;
    1169           0 :     case 1: return gen_1;
    1170             :   }
    1171           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1172             : }
    1173             : /* p > 2 */
    1174             : static GEN
    1175          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1176             : {
    1177          14 :   switch(d)
    1178             :   {
    1179           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1180           0 :     case 0: return trivial_fact();
    1181           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1182             :   }
    1183          14 :   switch(FpX_quad_factortype(f, p)) {
    1184           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1185           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1186           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1187             :   }
    1188             : }
    1189             : 
    1190             : GEN
    1191          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1192             : GEN
    1193       82903 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1194             : 
    1195             : /* not gerepile safe */
    1196             : static GEN
    1197         734 : FpX_factor_2(GEN f, GEN p, long d)
    1198             : {
    1199             :   GEN r, s, R, S;
    1200             :   long v;
    1201             :   int sgn;
    1202         734 :   switch(d)
    1203             :   {
    1204           7 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1205          30 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1206          45 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1207             :   }
    1208         652 :   r = FpX_quad_root(f, p, 1);
    1209         652 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1210         138 :   v = varn(f);
    1211         138 :   s = FpX_otherroot(f, r, p);
    1212         138 :   if (signe(r)) r = subii(p, r);
    1213         138 :   if (signe(s)) s = subii(p, s);
    1214         138 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1215         138 :   R = deg1pol_shallow(gen_1, r, v);
    1216         138 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1217          47 :   S = deg1pol_shallow(gen_1, s, v);
    1218          47 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1219             : }
    1220             : static GEN
    1221         755 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1222             : {
    1223         755 :   switch(flag) {
    1224           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1225          14 :     case 1: return FpX_degfact_2(f, p, d);
    1226         734 :     default: return FpX_factor_2(f, p, d);
    1227             :   }
    1228             : }
    1229             : 
    1230             : static int
    1231       73055 : F2x_quad_factortype(GEN x)
    1232       73055 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1233             : 
    1234             : static GEN
    1235           7 : F2x_is_irred_2(GEN f, long d)
    1236           7 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1237             : 
    1238             : static GEN
    1239       10059 : F2x_degfact_2(GEN f, long d)
    1240             : {
    1241       10059 :   if (!d) return trivial_fact();
    1242       10059 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1243        9856 :   switch(F2x_quad_factortype(f)) {
    1244        3213 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1245        3234 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1246        3409 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1247             :   }
    1248             : }
    1249             : 
    1250             : static GEN
    1251      117894 : F2x_factor_2(GEN f, long d)
    1252             : {
    1253      117894 :   long v = f[1];
    1254      117894 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1255      109598 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1256       61843 :   switch(F2x_quad_factortype(f))
    1257             :   {
    1258       14091 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1259       27844 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1260       19908 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1261             :   }
    1262             : }
    1263             : static GEN
    1264      127960 : F2x_factor_deg2(GEN f, long d, long flag)
    1265             : {
    1266      127960 :   switch(flag) {
    1267           7 :     case 2: return F2x_is_irred_2(f, d);
    1268       10059 :     case 1: return F2x_degfact_2(f, d);
    1269      117894 :     default: return F2x_factor_2(f, d);
    1270             :   }
    1271             : }
    1272             : 
    1273             : /* xt = NULL or x^(p-1)/2 mod g */
    1274             : static void
    1275        4951 : split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
    1276             : {
    1277        4951 :   ulong q = p >> 1;
    1278        4951 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1279        4951 :   long d = degpol(a);
    1280        4951 :   if (d < 0)
    1281             :   {
    1282             :     ulong i;
    1283          42 :     split_add_done(S, (GEN)1);
    1284          42 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1285             :   } else {
    1286        4909 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1287        4909 :     if (d)
    1288             :     {
    1289        4909 :       if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
    1290        4909 :       a = Flx_gcd(a, xt, p);
    1291        4909 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1292             :     }
    1293             :   }
    1294        4951 : }
    1295             : static void
    1296        4951 : split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
    1297             : {
    1298        4951 :   ulong q = p >> 1;
    1299        4951 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1300        4951 :   long d = degpol(a);
    1301        4951 :   if (d < 0)
    1302             :   {
    1303          28 :     ulong i, z = nonsquare_Fl(p);
    1304          28 :     split_add_done(S, (GEN)z);
    1305          28 :     for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1306             :   } else {
    1307        4923 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1308        4923 :     if (d)
    1309             :     {
    1310        4923 :       if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
    1311        4923 :       a = Flx_gcd(a, xt, p);
    1312        4923 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1313             :     }
    1314             :   }
    1315        4951 : }
    1316             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1317             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1318             : static int
    1319     4928950 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
    1320             : {
    1321     4928950 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1322     4928878 :   long d = degpol(g);
    1323     4928698 :   if (d < 0) return 0;
    1324     4928698 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1325     4928716 :   if (!d) return 1;
    1326     4913092 :   if ((p >> 4) <= (ulong)d)
    1327             :   { /* small p; split directly using x^((p-1)/2) +/- 1 */
    1328       14685 :     GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
    1329        9734 :                                 : NULL;
    1330        4951 :     split_squares(S, g, p, xt);
    1331        4951 :     split_nonsquares(S, g, p, xt);
    1332             :   } else { /* large p; use x^(p-1) - 1 directly */
    1333     4908141 :     a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
    1334     4900971 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1335     4900971 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1336     4901585 :     g = Flx_gcd(g,a, p);
    1337     4902897 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1338             :   }
    1339     4912466 :   return 1;
    1340             : }
    1341             : 
    1342             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1343             : static GEN
    1344    21999404 : Flx_roots_i(GEN f, ulong p)
    1345             : {
    1346             :   GEN pol, g;
    1347    21999404 :   long v = Flx_valrem(f, &g);
    1348             :   ulong q;
    1349             :   struct split_t S;
    1350             : 
    1351             :   /* optimization: test for small degree first */
    1352    22007114 :   switch(degpol(g))
    1353             :   {
    1354             :     case 1: {
    1355       25515 :       ulong r = p - g[2];
    1356       25515 :       return v? mkvecsmall2(0, r): mkvecsmall(r);
    1357             :     }
    1358             :     case 2: {
    1359    17041190 :       ulong r = Flx_quad_root(g, p, 1), s;
    1360    17072433 :       if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1361    11668762 :       s = Flx_otherroot(g,r, p);
    1362    11705429 :       if (r < s)
    1363     2922562 :         return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1364     8782867 :       else if (r > s)
    1365     8782706 :         return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1366             :       else
    1367         161 :         return v? mkvecsmall2(0, s): mkvecsmall(s);
    1368             :     }
    1369             :   }
    1370     4929614 :   q = p >> 1;
    1371     4929614 :   split_init(&S, lg(f)-1);
    1372     4928962 :   settyp(S.done, t_VECSMALL);
    1373     4928962 :   if (v) split_add_done(&S, (GEN)0);
    1374     4928962 :   if (! split_Flx_cut_out_roots(&S, g, p))
    1375          42 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1376     4927894 :   pol = polx_Flx(f[1]);
    1377    10185904 :   for (pol[2]=1; ; pol[2]++)
    1378     5257674 :   {
    1379    10185904 :     long j, l = lg(S.todo);
    1380    10185904 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1381     5257469 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1382    10871713 :     for (j = 1; j < l; j++)
    1383             :     {
    1384     5614039 :       GEN b, c = gel(S.todo,j);
    1385             :       ulong r, s;
    1386     5614039 :       switch(degpol(c))
    1387             :       {
    1388             :         case 1:
    1389     4804149 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1390     4804625 :           j--; l--; break;
    1391             :         case 2:
    1392      371329 :           r = Flx_quad_root(c, p, 0);
    1393      371342 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1394      371335 :           s = Flx_otherroot(c,r, p);
    1395      371338 :           split_done(&S, j, (GEN)r, (GEN)s);
    1396      371334 :           j--; l--; break;
    1397             :         default:
    1398      438571 :           b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
    1399      438548 :           if (degpol(b) <= 0) continue;
    1400      337395 :           b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
    1401      337393 :           if (!degpol(b)) continue;
    1402      337268 :           b = Flx_normalize(b, p);
    1403      337281 :           c = Flx_div(c,b, p);
    1404      337274 :           split_todo(&S, j, b, c);
    1405             :       }
    1406             :     }
    1407             :   }
    1408             : }
    1409             : 
    1410             : GEN
    1411    21940609 : Flx_roots(GEN f, ulong p)
    1412             : {
    1413    21940609 :   pari_sp av = avma;
    1414    21940609 :   switch(lg(f))
    1415             :   {
    1416           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1417           0 :     case 3: set_avma(av); return cgetg(1, t_VECSMALL);
    1418             :   }
    1419    21942915 :   if (p == 2) return Flx_root_mod_2(f);
    1420    21942915 :   return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
    1421             : }
    1422             : 
    1423             : /* assume x reduced mod p, monic. */
    1424             : static int
    1425     1082977 : Flx_quad_factortype(GEN x, ulong p)
    1426             : {
    1427     1082977 :   ulong b = x[3], c = x[2];
    1428     1082977 :   return krouu(Fl_disc_bc(b, c, p), p);
    1429             : }
    1430             : static GEN
    1431          56 : Flx_is_irred_2(GEN f, ulong p, long d)
    1432             : {
    1433          56 :   if (!d) return NULL;
    1434          56 :   if (d == 1) return gen_1;
    1435          56 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1436             : }
    1437             : static GEN
    1438     1100064 : Flx_degfact_2(GEN f, ulong p, long d)
    1439             : {
    1440     1100064 :   if (!d) return trivial_fact();
    1441     1100064 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1442     1082921 :   switch(Flx_quad_factortype(f, p)) {
    1443      516950 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1444      553084 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1445       12887 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1446             :   }
    1447             : }
    1448             : /* p > 2 */
    1449             : static GEN
    1450      460530 : Flx_factor_2(GEN f, ulong p, long d)
    1451             : {
    1452             :   ulong r, s;
    1453             :   GEN R,S;
    1454      460530 :   long v = f[1];
    1455      460530 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1456      445872 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1457      345336 :   r = Flx_quad_root(f, p, 1);
    1458      345357 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1459      233034 :   s = Flx_otherroot(f, r, p);
    1460      233033 :   r = Fl_neg(r, p);
    1461      233033 :   s = Fl_neg(s, p);
    1462      233033 :   if (s < r) lswap(s,r);
    1463      233033 :   R = mkvecsmall3(v,r,1);
    1464      233031 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1465      200279 :   S = mkvecsmall3(v,s,1);
    1466      200280 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1467             : }
    1468             : static GEN
    1469     1560651 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1470             : {
    1471     1560651 :   switch(flag) {
    1472          56 :     case 2: return Flx_is_irred_2(f, p, d);
    1473     1100064 :     case 1: return Flx_degfact_2(f, p, d);
    1474      460531 :     default: return Flx_factor_2(f, p, d);
    1475             :   }
    1476             : }
    1477             : 
    1478             : void
    1479       20062 : F2xV_to_FlxV_inplace(GEN v)
    1480             : {
    1481             :   long i;
    1482       20062 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_Flx(gel(v,i));
    1483       20062 : }
    1484             : void
    1485      764894 : FlxV_to_ZXV_inplace(GEN v)
    1486             : {
    1487             :   long i;
    1488      764894 :   for(i=1;i<lg(v);i++) gel(v,i)= Flx_to_ZX(gel(v,i));
    1489      764818 : }
    1490             : void
    1491      158099 : F2xV_to_ZXV_inplace(GEN v)
    1492             : {
    1493             :   long i;
    1494      158099 :   for(i=1;i<lg(v);i++) gel(v,i)= F2x_to_ZX(gel(v,i));
    1495      158099 : }
    1496             : 
    1497             : static GEN
    1498       10646 : F2x_Berlekamp_ker(GEN u)
    1499             : {
    1500       10646 :   pari_sp ltop=avma;
    1501       10646 :   long j,N = F2x_degree(u);
    1502             :   GEN Q;
    1503             :   pari_timer T;
    1504       10645 :   timer_start(&T);
    1505       10647 :   Q = F2x_matFrobenius(u);
    1506      249400 :   for (j=1; j<=N; j++)
    1507      238753 :     F2m_flip(Q,j,j);
    1508       10647 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
    1509       10647 :   Q = F2m_ker_sp(Q,0);
    1510       10647 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
    1511       10647 :   return gerepileupto(ltop,Q);
    1512             : }
    1513             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    1514             : static long
    1515       14622 : F2x_split_Berlekamp(GEN *t)
    1516             : {
    1517       14622 :   GEN u = *t, a, b, vker;
    1518       14622 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    1519             : 
    1520       14623 :   if (du == 1) return 1;
    1521       11262 :   if (du == 2)
    1522             :   {
    1523         616 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    1524             :     {
    1525           0 :       t[0] = mkvecsmall2(sv, 2);
    1526           0 :       t[1] = mkvecsmall2(sv, 3);
    1527           0 :       return 2;
    1528             :     }
    1529         616 :     return 1;
    1530             :   }
    1531             : 
    1532       10646 :   vker = F2x_Berlekamp_ker(u);
    1533       10647 :   lb = lgcols(vker);
    1534       10653 :   d = lg(vker)-1;
    1535       10653 :   ir = 0;
    1536             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    1537       53005 :   for (L=1; L<d; )
    1538             :   {
    1539             :     GEN pol;
    1540       31706 :     if (d == 2)
    1541        1742 :       pol = F2v_to_F2x(gel(vker,2), sv);
    1542             :     else
    1543             :     {
    1544       29964 :       GEN v = zero_zv(lb);
    1545       29971 :       v[1] = du;
    1546       29971 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    1547      112093 :       for (i=2; i<=d; i++)
    1548       82122 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    1549       29971 :       pol = F2v_to_F2x(v, sv);
    1550             :     }
    1551       98549 :     for (i=ir; i<L && L<d; i++)
    1552             :     {
    1553       66850 :       a = t[i]; la = F2x_degree(a);
    1554       66847 :       if (la == 1) { set_irred(i); }
    1555       66675 :       else if (la == 2)
    1556             :       {
    1557         733 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    1558             :         {
    1559           0 :           t[i] = mkvecsmall2(sv, 2);
    1560           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    1561             :         }
    1562         733 :         set_irred(i);
    1563             :       }
    1564             :       else
    1565             :       {
    1566       65942 :         pari_sp av = avma;
    1567             :         long lb;
    1568       65942 :         b = F2x_rem(pol, a);
    1569       65937 :         if (F2x_degree(b) <= 0) { set_avma(av); continue; }
    1570       21704 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    1571       21701 :         if (lb && lb < la)
    1572             :         {
    1573       21701 :           t[L] = F2x_div(a,b);
    1574       21698 :           t[i]= b; L++;
    1575             :         }
    1576           0 :         else set_avma(av);
    1577             :       }
    1578             :     }
    1579             :   }
    1580       10646 :   return d;
    1581             : }
    1582             : /* assume deg f > 2 */
    1583             : static GEN
    1584       11758 : F2x_Berlekamp_i(GEN f, long flag)
    1585             : {
    1586       11758 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    1587             :   GEN y, E, t, V;
    1588             : 
    1589       11757 :   val = F2x_valrem(f, &f);
    1590       11759 :   if (flag == 2 && val) return NULL;
    1591       11752 :   V = F2x_factor_squarefree(f); lV = lg(V);
    1592       11751 :   if (flag == 2 && lV > 2) return NULL;
    1593             : 
    1594             :   /* to hold factors and exponents */
    1595       11681 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    1596       11683 :   E = cgetg(d+1,t_VECSMALL);
    1597       11683 :   lfact = 1;
    1598       11683 :   if (val) {
    1599        3480 :     if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
    1600        3480 :     E[1] = val; lfact++;
    1601             :   }
    1602             : 
    1603       49357 :   for (k=1; k<lV; k++)
    1604             :   {
    1605       37827 :     if (F2x_degree(gel(V, k))==0) continue;
    1606       14622 :     gel(t,lfact) = gel(V, k);
    1607       14622 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    1608       14624 :     if (flag == 2 && d != 1) return NULL;
    1609       14470 :     if (flag == 1)
    1610        1589 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    1611       14470 :     for (j=0; j<d; j++) E[lfact+j] = k;
    1612       14470 :     lfact += d;
    1613             :   }
    1614       11530 :   if (flag == 2) return gen_1; /* irreducible */
    1615       11516 :   setlg(t, lfact);
    1616       11516 :   setlg(E, lfact); y = mkvec2(t,E);
    1617             :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    1618       11516 :               : sort_factor_pol(y, cmpGuGu);
    1619             : }
    1620             : 
    1621             : /* Adapted from Shoup NTL */
    1622             : GEN
    1623       77531 : F2x_factor_squarefree(GEN f)
    1624             : {
    1625             :   GEN r, t, v, tv;
    1626       77531 :   long i, q, n = F2x_degree(f);
    1627       77529 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1628      137281 :   for(q = 1;;q *= 2)
    1629             :   {
    1630      197030 :     r = F2x_gcd(f, F2x_deriv(f));
    1631      137278 :     if (F2x_degree(r) == 0)
    1632             :     {
    1633       63342 :       gel(u, q) = f;
    1634       63342 :       break;
    1635             :     }
    1636       73935 :     t = F2x_div(f, r);
    1637       73935 :     if (F2x_degree(t) > 0)
    1638             :     {
    1639             :       long j;
    1640       70589 :       for(j = 1;;j++)
    1641             :       {
    1642      109395 :         v = F2x_gcd(r, t);
    1643       70589 :         tv = F2x_div(t, v);
    1644       70589 :         if (F2x_degree(tv) > 0)
    1645       33142 :           gel(u, j*q) = tv;
    1646       70587 :         if (F2x_degree(v) <= 0) break;
    1647       38806 :         r = F2x_div(r, v);
    1648       38806 :         t = v;
    1649             :       }
    1650       31783 :       if (F2x_degree(r) == 0) break;
    1651             :     }
    1652       59749 :     f = F2x_sqrt(r);
    1653             :   }
    1654      565403 :   for (i = n; i; i--)
    1655      565078 :     if (F2x_degree(gel(u,i))) break;
    1656       77588 :   setlg(u,i+1); return u;
    1657             : }
    1658             : 
    1659             : static GEN
    1660       81471 : F2x_ddf_simple(GEN T, GEN XP)
    1661             : {
    1662       81471 :   pari_sp av = avma, av2;
    1663             :   GEN f, z, Tr, X;
    1664       81471 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1665       81472 :   if (n == 0) return cgetg(1, t_VEC);
    1666       81472 :   if (n == 1) return mkvec(T);
    1667       38396 :   z = XP; Tr = T; X = polx_F2x(v);
    1668       38397 :   f = const_vec(n, pol1_F2x(v));
    1669       38397 :   av2 = avma;
    1670      135570 :   for (j = 1; j <= B; j++)
    1671             :   {
    1672      101743 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1673      101755 :     if (F2x_degree(u))
    1674             :     {
    1675       22213 :       gel(f, j) = u;
    1676       22213 :       Tr = F2x_div(Tr, u);
    1677       22208 :       av2 = avma;
    1678       79537 :     } else z = gerepileuptoleaf(av2, z);
    1679      101771 :     if (!F2x_degree(Tr)) break;
    1680       97174 :     z = F2xq_sqr(z, Tr);
    1681             :   }
    1682       38392 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1683       38395 :   return gerepilecopy(av, f);
    1684             : }
    1685             : 
    1686             : GEN
    1687           7 : F2x_ddf(GEN T)
    1688             : {
    1689             :   GEN XP;
    1690           7 :   T = F2x_get_red(T);
    1691           7 :   XP = F2x_Frobenius(T);
    1692           7 :   return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
    1693             : }
    1694             : 
    1695             : static GEN
    1696        6623 : F2xq_frobtrace(GEN a, long d, GEN T)
    1697             : {
    1698        6623 :   pari_sp av = avma;
    1699             :   long i;
    1700        6623 :   GEN x = a;
    1701       27171 :   for(i=1; i<d; i++)
    1702             :   {
    1703       20548 :     x = F2x_add(a, F2xq_sqr(x,T));
    1704       20548 :     if (gc_needed(av, 2))
    1705           0 :       x = gerepileuptoleaf(av, x);
    1706             :   }
    1707        6623 :   return x;
    1708             : }
    1709             : 
    1710             : static void
    1711        9889 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1712             : {
    1713        9889 :   long n = F2x_degree(Tp), r = n/d;
    1714             :   GEN T, f, ff;
    1715        9889 :   if (r==1) { gel(V, idx) = Tp; return; }
    1716        3369 :   T = Tp;
    1717        3369 :   XP = F2x_rem(XP, T);
    1718             :   while (1)
    1719        3254 :   {
    1720        6623 :     pari_sp btop = avma;
    1721             :     long df;
    1722        6623 :     GEN g = random_F2x(n, Tp[1]);
    1723        6623 :     GEN t = F2xq_frobtrace(g, d, T);
    1724        6623 :     if (lgpol(t) == 0) continue;
    1725        4934 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1726        4934 :     if (df > 0 && df < n) break;
    1727        1565 :     set_avma(btop);
    1728             :   }
    1729        3369 :   ff = F2x_div(Tp, f);
    1730        3369 :   F2x_edf_simple(f, XP, d, V, idx);
    1731        3369 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1732             : }
    1733             : 
    1734             : static GEN
    1735       81465 : F2x_factor_Shoup(GEN T)
    1736             : {
    1737       81465 :   long i, n, s = 0;
    1738             :   GEN XP, D, V;
    1739             :   pari_timer ti;
    1740       81465 :   n = F2x_degree(T);
    1741       81464 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1742       81464 :   XP = F2x_Frobenius(T);
    1743       81464 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1744       81464 :   D = F2x_ddf_simple(T, XP);
    1745       81466 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1746      352905 :   for (i = 1; i <= n; i++)
    1747      271442 :     s += F2x_degree(gel(D,i))/i;
    1748       81463 :   V = cgetg(s+1, t_COL);
    1749      352881 :   for (i = 1, s = 1; i <= n; i++)
    1750             :   {
    1751      271416 :     GEN Di = gel(D,i);
    1752      271416 :     long ni = F2x_degree(Di), ri = ni/i;
    1753      271402 :     if (ni == 0) continue;
    1754       99099 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1755        3151 :     F2x_edf_simple(Di, XP, i, V, s);
    1756        3151 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1757        3151 :     s += ri;
    1758             :   }
    1759       81465 :   return V;
    1760             : }
    1761             : 
    1762             : static GEN
    1763       65778 : F2x_factor_Cantor(GEN T)
    1764             : {
    1765       65778 :   GEN E, F, V = F2x_factor_squarefree(T);
    1766       65775 :   long i, j, l = lg(V);
    1767       65775 :   E = cgetg(l, t_VEC);
    1768       65776 :   F = cgetg(l, t_VEC);
    1769      260975 :   for (i=1, j=1; i < l; i++)
    1770      195197 :     if (F2x_degree(gel(V,i)))
    1771             :     {
    1772       81465 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1773       81465 :       gel(F, j) = Fj;
    1774       81465 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1775       81466 :       j++;
    1776             :     }
    1777       65778 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1778             : }
    1779             : 
    1780             : #if 0
    1781             : static GEN
    1782             : F2x_simplefact_Shoup(GEN T)
    1783             : {
    1784             :   long i, n, s = 0, j = 1, k;
    1785             :   GEN XP, D, V;
    1786             :   pari_timer ti;
    1787             :   n = F2x_degree(T);
    1788             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1789             :   XP = F2x_Frobenius(T);
    1790             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1791             :   D = F2x_ddf_simple(T, XP);
    1792             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1793             :   for (i = 1; i <= n; i++)
    1794             :     s += F2x_degree(gel(D,i))/i;
    1795             :   V = cgetg(s+1, t_VECSMALL);
    1796             :   for (i = 1; i <= n; i++)
    1797             :   {
    1798             :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1799             :     if (ni == 0) continue;
    1800             :     for (k = 1; k <= ri; k++)
    1801             :       V[j++] = i;
    1802             :   }
    1803             :   return V;
    1804             : }
    1805             : static GEN
    1806             : F2x_simplefact_Cantor(GEN T)
    1807             : {
    1808             :   GEN E, F, V = F2x_factor_squarefree(T);
    1809             :   long i, j, l = lg(V);
    1810             :   F = cgetg(l, t_VEC);
    1811             :   E = cgetg(l, t_VEC);
    1812             :   for (i=1, j=1; i < l; i++)
    1813             :     if (F2x_degree(gel(V,i)))
    1814             :     {
    1815             :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1816             :       gel(F, j) = Fj;
    1817             :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1818             :       j++;
    1819             :     }
    1820             :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1821             : }
    1822             : static int
    1823             : F2x_isirred_Cantor(GEN T)
    1824             : {
    1825             :   pari_sp av = avma;
    1826             :   pari_timer ti;
    1827             :   long n;
    1828             :   GEN dT = F2x_deriv(T);
    1829             :   GEN XP, D;
    1830             :   if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
    1831             :   n = F2x_degree(T);
    1832             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1833             :   XP = F2x_Frobenius(T);
    1834             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1835             :   D = F2x_ddf_simple(T, XP);
    1836             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1837             :   return gc_bool(av, F2x_degree(gel(D,n)) == n);
    1838             : }
    1839             : #endif
    1840             : 
    1841             : /* driver for Cantor factorization, assume deg f > 2; not competitive for
    1842             :  * flag != 0, or as deg f increases */
    1843             : static GEN
    1844       65778 : F2x_Cantor_i(GEN f, long flag)
    1845             : {
    1846             :   switch(flag)
    1847             :   {
    1848       65778 :     default: return F2x_factor_Cantor(f);
    1849             : #if 0
    1850             :     case 1: return F2x_simplefact_Cantor(f);
    1851             :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1852             : #endif
    1853             :   }
    1854             : }
    1855             : static GEN
    1856      205496 : F2x_factor_i(GEN f, long flag)
    1857             : {
    1858      205496 :   long d = F2x_degree(f);
    1859      205495 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1860       75463 :   return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
    1861      143313 :                                : F2x_Berlekamp_i(f, flag);
    1862             : }
    1863             : 
    1864             : GEN
    1865           0 : F2x_degfact(GEN f)
    1866             : {
    1867           0 :   pari_sp av = avma;
    1868           0 :   GEN z = F2x_factor_i(f, 1);
    1869           0 :   return gerepilecopy(av, z);
    1870             : }
    1871             : 
    1872             : int
    1873         238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
    1874             : 
    1875             : /* Adapted from Shoup NTL */
    1876             : GEN
    1877      791248 : Flx_factor_squarefree(GEN f, ulong p)
    1878             : {
    1879      791248 :   long i, q, n = degpol(f);
    1880      791246 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1881      845653 :   for(q = 1;;q *= p)
    1882       54402 :   {
    1883      845653 :     GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
    1884      845626 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1885      100793 :     t = Flx_div(f, r, p);
    1886      100793 :     if (degpol(t) > 0)
    1887             :     {
    1888             :       long j;
    1889      149252 :       for(j = 1;;j++)
    1890             :       {
    1891      251324 :         v = Flx_gcd(r, t, p);
    1892      149252 :         tv = Flx_div(t, v, p);
    1893      149252 :         if (degpol(tv) > 0)
    1894       71629 :           gel(u, j*q) = Flx_normalize(tv, p);
    1895      149251 :         if (degpol(v) <= 0) break;
    1896      102072 :         r = Flx_div(r, v, p);
    1897      102072 :         t = v;
    1898             :       }
    1899       47180 :       if (degpol(r) == 0) break;
    1900             :     }
    1901       54402 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1902             :   }
    1903     4299981 :   for (i = n; i; i--)
    1904     4299971 :     if (degpol(gel(u,i))) break;
    1905      791256 :   setlg(u,i+1); return u;
    1906             : }
    1907             : 
    1908             : long
    1909         126 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1910             : {
    1911         126 :   pari_sp av = avma;
    1912             :   ulong lc;
    1913             :   GEN F;
    1914         126 :   long i, n = degpol(f), v = f[1], l;
    1915         126 :   if (n % k) return 0;
    1916         126 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    1917         126 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    1918         126 :   F = Flx_factor_squarefree(f, p); l = lg(F)-1;
    1919        6405 :   for (i = 1; i <= l; i++)
    1920        6300 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1921         105 :   if (pt_r)
    1922             :   {
    1923         105 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    1924        6384 :     for(i = l; i >= 1; i--)
    1925             :     {
    1926        6279 :       if (i%k) continue;
    1927        1274 :       s = Flx_mul(s, gel(F,i), p);
    1928        1274 :       r = Flx_mul(r, s, p);
    1929             :     }
    1930         105 :     *pt_r = gerepileuptoleaf(av, r);
    1931           0 :   } else set_avma(av);
    1932         105 :   return 1;
    1933             : }
    1934             : 
    1935             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    1936             : static GEN
    1937     1172293 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
    1938             : {
    1939     1172293 :   pari_sp av = avma;
    1940             :   GEN b, g, h, F, f, Tr, xq;
    1941             :   long i, j, n, v, bo, ro;
    1942             :   long B, l, m;
    1943             :   pari_timer ti;
    1944     1172293 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    1945     1172295 :   if (n == 0) return cgetg(1, t_VEC);
    1946     1170112 :   if (n == 1) return mkvec(get_Flx_mod(T));
    1947     1071556 :   B = n/2;
    1948     1071556 :   l = usqrt(B);
    1949     1071558 :   m = (B+l-1)/l;
    1950     1071558 :   T = Flx_get_red(T, p);
    1951     1071534 :   b = cgetg(l+2, t_VEC);
    1952     1071543 :   gel(b, 1) = polx_Flx(v);
    1953     1071557 :   gel(b, 2) = XP;
    1954     1071557 :   bo = brent_kung_optpow(n, l-1, 1);
    1955     1071566 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    1956     1071566 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1957     1071566 :   if (expu(p) <= ro)
    1958      239516 :     for (i = 3; i <= l+1; i++)
    1959      135077 :       gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
    1960             :   else
    1961             :   {
    1962      967130 :     xq = Flxq_powers(gel(b, 2), bo,  T, p);
    1963      967099 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
    1964     1112259 :     for (i = 3; i <= l+1; i++)
    1965      145150 :       gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
    1966             :   }
    1967     1071548 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
    1968     1071548 :   xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
    1969     1071543 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
    1970     1071543 :   g = cgetg(m+1, t_VEC);
    1971     1071563 :   gel(g, 1) = gel(xq, 2);
    1972     1953347 :   for(i = 2; i <= m; i++)
    1973      881796 :     gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
    1974     1071551 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
    1975     1071551 :   h = cgetg(m+1, t_VEC);
    1976     3024913 :   for (j = 1; j <= m; j++)
    1977             :   {
    1978     1953353 :     pari_sp av = avma;
    1979     1953353 :     GEN gj = gel(g, j);
    1980     1953353 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    1981     2996254 :     for (i = 2; i <= l; i++)
    1982     1042917 :       e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
    1983     1953337 :     gel(h, j) = gerepileupto(av, e);
    1984             :   }
    1985     1071560 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
    1986     1071560 :   Tr = get_Flx_mod(T);
    1987     1071555 :   F = cgetg(m+1, t_VEC);
    1988     3024850 :   for (j = 1; j <= m; j++)
    1989             :   {
    1990     1953305 :     GEN u = Flx_gcd(Tr, gel(h, j), p);
    1991     1953279 :     if (degpol(u))
    1992             :     {
    1993      811125 :       u = Flx_normalize(u, p);
    1994      811126 :       Tr = Flx_div(Tr, u, p);
    1995             :     }
    1996     1953263 :     gel(F, j) = u;
    1997             :   }
    1998     1071545 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
    1999     1071545 :   f = const_vec(n, pol1_Flx(v));
    2000     3024863 :   for (j = 1; j <= m; j++)
    2001             :   {
    2002     1953311 :     GEN e = gel(F, j);
    2003     2145143 :     for (i=l-1; i >= 0; i--)
    2004             :     {
    2005     2145127 :       GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
    2006     2145051 :       if (degpol(u))
    2007             :       {
    2008      832595 :         gel(f, l*j-i) = u;
    2009      832595 :         e = Flx_div(e, u, p);
    2010             :       }
    2011     2145042 :       if (!degpol(e)) break;
    2012             :     }
    2013             :   }
    2014     1071552 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
    2015     1071552 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    2016     1071554 :   return gerepilecopy(av, f);
    2017             : }
    2018             : 
    2019             : static void
    2020      131670 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2021             : {
    2022      131670 :   long n = degpol(Tp), r = n/d;
    2023             :   GEN T, f, ff;
    2024             :   ulong p2;
    2025      131670 :   if (r==1) { gel(V, idx) = Tp; return; }
    2026       56864 :   p2 = p>>1;
    2027       56864 :   T = Flx_get_red(Tp, p);
    2028       56863 :   XP = Flx_rem(XP, T, p);
    2029             :   while (1)
    2030        6563 :   {
    2031       63427 :     pari_sp btop = avma;
    2032             :     long i;
    2033       63427 :     GEN g = random_Flx(n, Tp[1], p);
    2034       63427 :     GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2035       63427 :     if (lgpol(t) == 0) continue;
    2036      139503 :     for(i=1; i<=10; i++)
    2037             :     {
    2038      134610 :       pari_sp btop2 = avma;
    2039      134610 :       GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
    2040      134611 :       f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
    2041      134609 :       if (degpol(f) > 0 && degpol(f) < n) break;
    2042       77745 :       set_avma(btop2);
    2043             :     }
    2044       61757 :     if (degpol(f) > 0 && degpol(f) < n) break;
    2045        4893 :     set_avma(btop);
    2046             :   }
    2047       56864 :   f = Flx_normalize(f, p);
    2048       56863 :   ff = Flx_div(Tp, f ,p);
    2049       56862 :   Flx_edf_simple(f, XP, d, p, V, idx);
    2050       56862 :   Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
    2051             : }
    2052             : static void
    2053             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
    2054             : 
    2055             : static void
    2056      354886 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
    2057             : {
    2058             :   pari_sp av;
    2059      354886 :   GEN Tp = get_Flx_mod(T);
    2060      354886 :   long n = degpol(hp), vT = Tp[1];
    2061             :   GEN u1, u2, f1, f2;
    2062      354887 :   ulong p2 = p>>1;
    2063             :   GEN R, h;
    2064      354887 :   h = Flx_get_red(hp, p);
    2065      354869 :   t = Flx_rem(t, T, p);
    2066      354813 :   av = avma;
    2067             :   do
    2068             :   {
    2069      582621 :     set_avma(av);
    2070      582625 :     R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
    2071      582558 :     u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
    2072      582669 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2073      354879 :   f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
    2074      354879 :   f1 = Flx_normalize(f1, p);
    2075      354907 :   u2 = Flx_div(hp, u1, p);
    2076      354895 :   f2 = Flx_div(Tp, f1, p);
    2077      354871 :   if (degpol(u1)==1)
    2078             :   {
    2079      257254 :     if (degpol(f1)==d)
    2080      252836 :       gel(V, idx) = f1;
    2081             :     else
    2082        4409 :       Flx_edf(f1, XP, d, p, V, idx);
    2083             :   }
    2084             :   else
    2085       97657 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
    2086      354896 :   idx += degpol(f1)/d;
    2087      354878 :   if (degpol(u2)==1)
    2088             :   {
    2089      255427 :     if (degpol(f2)==d)
    2090      251437 :       gel(V, idx) = f2;
    2091             :     else
    2092        3989 :       Flx_edf(f2, XP, d, p, V, idx);
    2093             :   }
    2094             :   else
    2095       99491 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
    2096      354927 : }
    2097             : 
    2098             : static void
    2099      157802 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
    2100             : {
    2101      157802 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2102             :   GEN T, h, t;
    2103             :   pari_timer ti;
    2104      157801 :   if (r==1) { gel(V, idx) = Tp; return; }
    2105      157801 :   T = Flx_get_red(Tp, p);
    2106      157791 :   XP = Flx_rem(XP, T, p);
    2107      157789 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2108             :   do
    2109             :   {
    2110      162869 :     GEN g = random_Flx(n, vT, p);
    2111      162877 :     t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
    2112      162884 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2113      162884 :     h = Flxq_minpoly(t, T, p);
    2114      162878 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2115      162878 :   } while (degpol(h) <= 1);
    2116      157798 :   Flx_edf_rec(T, XP, h, t, d, p, V, idx);
    2117             : }
    2118             : 
    2119             : static GEN
    2120      447406 : Flx_factor_Shoup(GEN T, ulong p)
    2121             : {
    2122      447406 :   long i, n, s = 0;
    2123             :   GEN XP, D, V;
    2124      447406 :   long e = expu(p);
    2125             :   pari_timer ti;
    2126      447410 :   n = get_Flx_degree(T);
    2127      447406 :   T = Flx_get_red(T, p);
    2128      447371 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2129      447371 :   XP = Flx_Frobenius(T, p);
    2130      447404 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2131      447404 :   D = Flx_ddf_Shoup(T, XP, p);
    2132      447426 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2133      447426 :   s = ddf_to_nbfact(D);
    2134      447427 :   V = cgetg(s+1, t_COL);
    2135     2518191 :   for (i = 1, s = 1; i <= n; i++)
    2136             :   {
    2137     2070776 :     GEN Di = gel(D,i);
    2138     2070776 :     long ni = degpol(Di), ri = ni/i;
    2139     2070777 :     if (ni == 0) continue;
    2140      563024 :     Di = Flx_normalize(Di, p);
    2141      563025 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2142      167351 :     if (ri <= e*expu(e))
    2143      149405 :       Flx_edf(Di, XP, i, p, V, s);
    2144             :     else
    2145       17946 :       Flx_edf_simple(Di, XP, i, p, V, s);
    2146      167336 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2147      167336 :     s += ri;
    2148             :   }
    2149      447415 :   return V;
    2150             : }
    2151             : 
    2152             : static GEN
    2153      423683 : Flx_factor_Cantor(GEN T, ulong p)
    2154             : {
    2155      423683 :   GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
    2156      423672 :   long i, j, l = lg(V);
    2157      423672 :   F = cgetg(l, t_VEC);
    2158      423675 :   E = cgetg(l, t_VEC);
    2159     1066107 :   for (i=1, j=1; i < l; i++)
    2160      642420 :     if (degpol(gel(V,i)))
    2161             :     {
    2162      447407 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p);
    2163      447414 :       gel(F, j) = Fj;
    2164      447414 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2165      447418 :       j++;
    2166             :     }
    2167      423687 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2168             : }
    2169             : 
    2170             : GEN
    2171           0 : Flx_ddf(GEN T, ulong p)
    2172             : {
    2173             :   GEN XP;
    2174           0 :   T = Flx_get_red(T, p);
    2175           0 :   XP = Flx_Frobenius(T, p);
    2176           0 :   return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
    2177             : }
    2178             : 
    2179             : static GEN
    2180      366665 : Flx_simplefact_Cantor(GEN T, ulong p)
    2181             : {
    2182             :   GEN V;
    2183             :   long i, l;
    2184      366665 :   T = Flx_get_red(T, p);
    2185      366665 :   V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
    2186      736509 :   for (i=1; i < l; i++)
    2187      369844 :     gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
    2188      366665 :   return vddf_to_simplefact(V, get_Flx_degree(T));
    2189             : }
    2190             : 
    2191             : static int
    2192        1057 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2193             : {
    2194        1057 :   pari_sp av = avma;
    2195             :   pari_timer ti;
    2196             :   long n;
    2197        1057 :   GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
    2198        1057 :   if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    2199         770 :   n = get_Flx_degree(T);
    2200         770 :   T = Flx_get_red(Tp, p);
    2201         770 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2202         770 :   XP = Flx_Frobenius(T, p);
    2203         770 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2204         770 :   D = Flx_ddf_Shoup(T, XP, p);
    2205         770 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2206         770 :   return gc_bool(av, degpol(gel(D,n)) == n);
    2207             : }
    2208             : 
    2209             : /* f monic */
    2210             : static GEN
    2211     2383996 : Flx_factor_i(GEN f, ulong pp, long flag)
    2212             : {
    2213             :   long d;
    2214     2383996 :   if (pp==2) { /*We need to handle 2 specially */
    2215       31941 :     GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
    2216       31941 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2217       31941 :     return F;
    2218             :   }
    2219     2352055 :   d = degpol(f);
    2220     2352052 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2221      791404 :   switch(flag)
    2222             :   {
    2223      423682 :     default: return Flx_factor_Cantor(f, pp);
    2224      366665 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2225        1057 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2226             :   }
    2227             : }
    2228             : 
    2229             : GEN
    2230     1478601 : Flx_degfact(GEN f, ulong p)
    2231             : {
    2232     1478601 :   pari_sp av = avma;
    2233     1478601 :   GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
    2234     1478601 :   return gerepilecopy(av, z);
    2235             : }
    2236             : 
    2237             : /* T must be squarefree mod p*/
    2238             : GEN
    2239      251338 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2240             : {
    2241             :   GEN XP, D;
    2242             :   pari_timer ti;
    2243      251338 :   long i, s, n = get_Flx_degree(T);
    2244      251338 :   GEN V = const_vecsmall(n, 0);
    2245      251338 :   pari_sp av = avma;
    2246      251338 :   T = Flx_get_red(T, p);
    2247      251338 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2248      251338 :   XP = Flx_Frobenius(T, p);
    2249      251338 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2250      251338 :   D = Flx_ddf_Shoup(T, XP, p);
    2251      251338 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2252     1415376 :   for (i = 1, s = 0; i <= n; i++)
    2253             :   {
    2254     1164038 :     V[i] = degpol(gel(D,i))/i;
    2255     1164038 :     s += V[i];
    2256             :   }
    2257      251338 :   *nb = s;
    2258      251338 :   set_avma(av); return V;
    2259             : }
    2260             : 
    2261             : long
    2262      102935 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2263             : {
    2264      102935 :   pari_sp av = avma;
    2265      102935 :   long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
    2266      102935 :   return gc_long(av,s);
    2267             : }
    2268             : 
    2269             : /* T must be squarefree mod p*/
    2270             : long
    2271      102935 : Flx_nbfact(GEN T, ulong p)
    2272             : {
    2273      102935 :   pari_sp av = avma;
    2274      102935 :   GEN XP = Flx_Frobenius(T, p);
    2275      102935 :   long n = Flx_nbfact_Frobenius(T, XP, p);
    2276      102935 :   return gc_long(av,n);
    2277             : }
    2278             : 
    2279             : int
    2280        1057 : Flx_is_irred(GEN f, ulong p)
    2281             : {
    2282        1057 :   pari_sp av = avma;
    2283        1057 :   f = Flx_normalize(f,p);
    2284        1057 :   return gc_bool(av, !!Flx_factor_i(f,p,2));
    2285             : }
    2286             : 
    2287             : /* Use this function when you think f is reducible, and that there are lots of
    2288             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2289             : int
    2290          77 : FpX_is_irred(GEN f, GEN p)
    2291             : {
    2292          77 :   pari_sp av = avma;
    2293             :   int z;
    2294          77 :   switch(ZX_factmod_init(&f,p))
    2295             :   {
    2296          14 :     case 0:  z = !!F2x_factor_i(f,2); break;
    2297          56 :     case 1:  z = !!Flx_factor_i(f,p[2],2); break;
    2298           7 :     default: z = !!FpX_factor_i(f,p,2); break;
    2299             :   }
    2300          77 :   return gc_bool(av,z);
    2301             : }
    2302             : GEN
    2303          42 : FpX_degfact(GEN f, GEN p) {
    2304          42 :   pari_sp av = avma;
    2305             :   GEN F;
    2306          42 :   switch(ZX_factmod_init(&f,p))
    2307             :   {
    2308           7 :     case 0:  F = F2x_factor_i(f,1); break;
    2309           7 :     case 1:  F = Flx_factor_i(f,p[2],1); break;
    2310          28 :     default: F = FpX_factor_i(f,p,1); break;
    2311             :   }
    2312          42 :   return gerepilecopy(av, F);
    2313             : }
    2314             : 
    2315             : #if 0
    2316             : /* set x <-- x + c*y mod p */
    2317             : /* x is not required to be normalized.*/
    2318             : static void
    2319             : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2320             : {
    2321             :   long i, lx, ly;
    2322             :   ulong *x=(ulong *)gx;
    2323             :   ulong *y=(ulong *)gy;
    2324             :   if (!c) return;
    2325             :   lx = lg(gx);
    2326             :   ly = lg(gy);
    2327             :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2328             :   if (SMALL_ULONG(p))
    2329             :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2330             :   else
    2331             :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2332             : }
    2333             : #endif
    2334             : 
    2335             : GEN
    2336      921749 : FpX_factor(GEN f, GEN p)
    2337             : {
    2338      921749 :   pari_sp av = avma;
    2339             :   GEN F;
    2340      921749 :   switch(ZX_factmod_init(&f, p))
    2341             :   {
    2342      158092 :     case 0:  F = F2x_factor_i(f,0);
    2343      158092 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    2344      762313 :     case 1:  F = Flx_factor_i(f,p[2],0);
    2345      762322 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    2346        1350 :     default: F = FpX_factor_i(f,p,0); break;
    2347             :   }
    2348      921753 :   return gerepilecopy(av, F);
    2349             : }
    2350             : 
    2351             : GEN
    2352      141965 : Flx_factor(GEN f, ulong p)
    2353             : {
    2354      141965 :   pari_sp av = avma;
    2355      141965 :   return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
    2356             : }
    2357             : GEN
    2358       15202 : F2x_factor(GEN f)
    2359             : {
    2360       15202 :   pari_sp av = avma;
    2361       15202 :   return gerepilecopy(av, F2x_factor_i(f,0));
    2362             : }

Generated by: LCOV version 1.13