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.18.0 lcov report (development 29792-e7296d7d08) Lines: 1308 1426 91.7 %
Date: 2024-12-05 09:09:40 Functions: 114 124 91.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_factormod
      19             : 
      20             : /***********************************************************************/
      21             : /**                                                                   **/
      22             : /**               Factorisation over finite field                     **/
      23             : /**                                                                   **/
      24             : /***********************************************************************/
      25             : 
      26             : /*******************************************************************/
      27             : /*                                                                 */
      28             : /*           ROOTS MODULO a prime p (no multiplicities)            */
      29             : /*                                                                 */
      30             : /*******************************************************************/
      31             : /* Replace F by a monic normalized FpX having the same factors;
      32             :  * assume p prime and *F a ZX */
      33             : static int
      34     4343406 : ZX_factmod_init(GEN *F, GEN p)
      35             : {
      36     4343406 :   if (lgefint(p) == 3)
      37             :   {
      38     4340899 :     ulong pp = p[2];
      39     4340899 :     if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
      40     2910590 :     *F = ZX_to_Flx(*F, pp);
      41     2910650 :     if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
      42     2910652 :     return 1;
      43             :   }
      44        2507 :   *F = FpX_red(*F, p);
      45        2508 :   if (lg(*F) > 3) *F = FpX_normalize(*F, p);
      46        2508 :   return 2;
      47             : }
      48             : static GEN
      49      659524 : ZX_rootmod_init(GEN F, GEN p)
      50      659524 : { return lgefint(p) == 3? ZX_to_Flx(F, p[2]): FpX_red(F, p); }
      51             : 
      52             : /* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
      53             : static GEN
      54         600 : all_roots_mod_p(ulong p, int not_0)
      55             : {
      56             :   GEN r;
      57             :   ulong i;
      58         600 :   if (not_0) {
      59         412 :     r = cgetg(p, t_VECSMALL);
      60        1276 :     for (i = 1; i < p; i++) r[i] = i;
      61             :   } else {
      62         188 :     r = cgetg(p+1, t_VECSMALL);
      63         780 :     for (i = 0; i < p; i++) r[i+1] = i;
      64             :   }
      65         600 :   return r;
      66             : }
      67             : 
      68             : /* X^n - 1 */
      69             : static GEN
      70        2186 : Flx_Xnm1(long sv, long n, ulong p)
      71             : {
      72        2186 :   GEN t = cgetg(n+3, t_VECSMALL);
      73             :   long i;
      74        2186 :   t[1] = sv;
      75        2186 :   t[2] = p - 1;
      76        6927 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      77        2186 :   t[i] = 1; return t;
      78             : }
      79             : /* X^n + 1 */
      80             : static GEN
      81        2076 : Flx_Xn1(long sv, long n, ulong p)
      82             : {
      83        2076 :   GEN t = cgetg(n+3, t_VECSMALL);
      84             :   long i;
      85             :   (void) p;
      86        2076 :   t[1] = sv;
      87        2076 :   t[2] = 1;
      88        6729 :   for (i = 3; i <= n+1; i++) t[i] = 0;
      89        2076 :   t[i] = 1; return t;
      90             : }
      91             : 
      92             : /* assume lg(f) > 3 */
      93             : static GEN
      94       92718 : Flx_root_mod_2(GEN f)
      95             : {
      96       92718 :   long i, n = lg(f)-1, c = f[2];
      97       92718 :   int z0 = !c;
      98       92718 :   c ^= 1; /* c = f[2] + f[n] mod 2, we know f[n] is odd */
      99      166801 :   for (i=3; i < n; i++) c ^= f[i];
     100             :   /* c = 0 iff f(1) = 0 (mod 2) */
     101       92718 :   if (z0) return c? mkvecsmall(0): mkvecsmall2(0, 1);
     102       13931 :   return c? cgetg(1, t_VECSMALL): mkvecsmall(1);
     103             : }
     104             : /* assume lg(f) > 3 */
     105             : static ulong
     106          91 : Flx_oneroot_mod_2(GEN f)
     107             : {
     108          91 :   long i, n, c = f[2];
     109          91 :   if (!c) return 0;
     110          91 :   n = lg(f)-1; c = 0; /* = f[2] + f[n] (mod 2); both are odd */
     111         182 :   for (i=3; i < n; i++) c ^= f[i];
     112          91 :   return c? 2: 1;
     113             : }
     114             : 
     115             : static GEN FpX_roots_i(GEN f, GEN p);
     116             : 
     117             : static int
     118    16300265 : cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
     119             : 
     120             : /* assume that f is a ZX and p a prime */
     121             : GEN
     122      425640 : FpX_roots(GEN f, GEN p)
     123             : {
     124      425640 :   pari_sp av = avma;
     125      425640 :   GEN y; f = ZX_rootmod_init(f, p);
     126      425638 :   switch(lg(f))
     127             :   {
     128          14 :     case 2: pari_err_ROOTS0("FpX_roots");
     129       47438 :     case 3: return cgetg(1,t_COL);
     130             :   }
     131      378186 :   if (typ(f) == t_VECSMALL)
     132             :   {
     133      369946 :     ulong pp = p[2];
     134      369946 :     if (pp == 2)
     135       92711 :       y = Flx_root_mod_2(f);
     136             :     else
     137             :     {
     138      277235 :       if (!odd(pp)) pari_err_PRIME("FpX_roots", p);
     139      277236 :       y = Flx_roots_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
     140             :     }
     141      369940 :     y = Flc_to_ZC(y);
     142             :   }
     143             :   else
     144        8240 :     y = FpX_roots_i(f, p);
     145      378173 :   return gerepileupto(av, y);
     146             : }
     147             : 
     148             : /* assume x reduced mod p > 2, monic. */
     149             : static int
     150          21 : FpX_quad_factortype(GEN x, GEN p)
     151             : {
     152          21 :   GEN b = gel(x,3), c = gel(x,2);
     153          21 :   GEN D = subii(sqri(b), shifti(c,2));
     154          21 :   return kronecker(D,p);
     155             : }
     156             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     157             : static GEN
     158       14457 : FpX_quad_root(GEN x, GEN p, int unknown)
     159             : {
     160       14457 :   GEN s, D, b = gel(x,3), c = gel(x,2);
     161             : 
     162       14457 :   if (absequaliu(p, 2)) {
     163           0 :     if (!signe(b)) return c;
     164           0 :     return signe(c)? NULL: gen_1;
     165             :   }
     166       14457 :   D = subii(sqri(b), shifti(c,2));
     167       14457 :   D = remii(D,p);
     168       14457 :   if (unknown && kronecker(D,p) == -1) return NULL;
     169             : 
     170       13867 :   s = Fp_sqrt(D,p);
     171             :   /* p is not prime, go on and give e.g. maxord a chance to recover */
     172       13867 :   if (!s) return NULL;
     173       13859 :   return Fp_halve(Fp_sub(s,b, p), p);
     174             : }
     175             : static GEN
     176        9994 : FpX_otherroot(GEN x, GEN r, GEN p)
     177        9994 : { return Fp_neg(Fp_add(gel(x,3), r, p), p); }
     178             : 
     179             : /* disc(x^2+bx+c) = b^2 - 4c */
     180             : static ulong
     181    27414882 : Fl_disc_bc(ulong b, ulong c, ulong p)
     182    27414882 : { return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
     183             : /* p > 2; allow pi = 0 */
     184             : static ulong
     185    24927880 : Flx_quad_root(GEN x, ulong p, ulong pi, int unknown)
     186             : {
     187    24927880 :   ulong s, b = x[3], c = x[2];
     188    24927880 :   ulong D = Fl_disc_bc(b, c, p);
     189    24900721 :   if (unknown && krouu(D,p) == -1) return p;
     190    16455829 :   s = Fl_sqrt_pre(D, p, pi);
     191    16559880 :   if (s==~0UL) return p;
     192    16559867 :   return Fl_halve(Fl_sub(s,b, p), p);
     193             : }
     194             : static ulong
     195    14751824 : Flx_otherroot(GEN x, ulong r, ulong p)
     196    14751824 : { return Fl_neg(Fl_add(x[3], r, p), p); }
     197             : 
     198             : /* 'todo' contains the list of factors to be split.
     199             :  * 'done' the list of finished factors, no longer touched */
     200             : struct split_t { GEN todo, done; };
     201             : static void
     202     5075861 : split_init(struct split_t *S, long max)
     203             : {
     204     5075861 :   S->todo = vectrunc_init(max);
     205     5075476 :   S->done = vectrunc_init(max);
     206     5075014 : }
     207             : #if 0
     208             : /* move todo[i] to done */
     209             : static void
     210             : split_convert(struct split_t *S, long i)
     211             : {
     212             :   long n = lg(S->todo)-1;
     213             :   vectrunc_append(S->done, gel(S->todo,i));
     214             :   if (n) gel(S->todo,i) = gel(S->todo, n);
     215             :   setlg(S->todo, n);
     216             : }
     217             : #endif
     218             : /* append t to todo */
     219             : static void
     220     5443271 : split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
     221             : /* delete todo[i], add t to done */
     222             : static void
     223     5443911 : split_moveto_done(struct split_t *S, long i, GEN t)
     224             : {
     225     5443911 :   long n = lg(S->todo)-1;
     226     5443911 :   vectrunc_append(S->done, t);
     227     5444043 :   if (n) gel(S->todo,i) = gel(S->todo, n);
     228     5444043 :   setlg(S->todo, n);
     229             : 
     230     5443906 : }
     231             : /* append t to done */
     232             : static void
     233      506225 : split_add_done(struct split_t *S, GEN t)
     234      506225 : { vectrunc_append(S->done, t); }
     235             : /* split todo[i] into a and b */
     236             : static void
     237      416617 : split_todo(struct split_t *S, long i, GEN a, GEN b)
     238             : {
     239      416617 :   gel(S->todo, i) = a;
     240      416617 :   split_add(S, b);
     241      416618 : }
     242             : /* split todo[i] into a and b, moved to done */
     243             : static void
     244      466352 : split_done(struct split_t *S, long i, GEN a, GEN b)
     245             : {
     246      466352 :   split_moveto_done(S, i, a);
     247      466358 :   split_add_done(S, b);
     248      466361 : }
     249             : 
     250             : /* by splitting, assume p > 2 prime, deg(f) > 0 */
     251             : static GEN
     252        8240 : FpX_roots_i(GEN f, GEN p)
     253             : {
     254             :   GEN pol, pol0, a, q;
     255             :   struct split_t S;
     256             : 
     257        8240 :   f = FpX_normalize(f, p);
     258        8240 :   split_init(&S, lg(f)-1);
     259        8240 :   settyp(S.done, t_COL);
     260        8240 :   if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
     261        8240 :   switch(degpol(f))
     262             :   {
     263           7 :     case 0: return ZC_copy(S.done);
     264          14 :     case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
     265        3524 :     case 2: {
     266        3524 :       GEN s, r = FpX_quad_root(f, p, 1);
     267        3524 :       if (r) {
     268        3524 :         split_add_done(&S, r);
     269        3524 :         s = FpX_otherroot(f,r, p);
     270             :         /* f not known to be square free yet */
     271        3524 :         if (!equalii(r, s)) split_add_done(&S, s);
     272             :       }
     273        3524 :       return sort(S.done);
     274             :     }
     275             :   }
     276             : 
     277        4695 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     278        4695 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     279        4695 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     280        4695 :   a = FpX_gcd(f,a, p);
     281        4695 :   if (!degpol(a)) return ZC_copy(S.done);
     282        4315 :   split_add(&S, FpX_normalize(a,p));
     283             : 
     284        4315 :   q = shifti(p,-1);
     285        4315 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     286        4315 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     287        4315 :   for (pol0[2] = 1;; pol0[2]++)
     288       10749 :   {
     289       15064 :     long j, l = lg(S.todo);
     290       15064 :     if (l == 1) return sort(S.done);
     291       10756 :     if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
     292       28738 :     for (j = 1; j < l; j++)
     293             :     {
     294       17989 :       GEN b, r, s, c = gel(S.todo,j);
     295       17989 :       switch(degpol(c))
     296             :       { /* convert linear and quadratics to roots, try to split the rest */
     297        4330 :         case 1:
     298        4330 :           split_moveto_done(&S, j, subii(p, gel(c,2)));
     299        4330 :           j--; l--; break;
     300        6219 :         case 2:
     301        6219 :           r = FpX_quad_root(c, p, 0);
     302        6219 :           if (!r) pari_err_PRIME("polrootsmod",p);
     303        6212 :           s = FpX_otherroot(c,r, p);
     304        6212 :           split_done(&S, j, r, s);
     305        6212 :           j--; l--; break;
     306        7440 :         default:
     307        7440 :           b = FpXQ_pow(pol,q, c,p);
     308        7440 :           if (degpol(b) <= 0) continue;
     309        6234 :           b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
     310        6234 :           if (!degpol(b)) continue;
     311        6234 :           b = FpX_normalize(b, p);
     312        6234 :           c = FpX_div(c,b, p);
     313        6234 :           split_todo(&S, j, b, c);
     314             :       }
     315             :     }
     316             :   }
     317             : }
     318             : 
     319             : /* Assume f is normalized; allow pi = 0 */
     320             : static ulong
     321      416410 : Flx_cubic_root(GEN ff, ulong p, ulong pi)
     322             : {
     323      416410 :   GEN f = Flx_normalize(ff,p);
     324      416408 :   ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
     325             :   ulong t, t2, A, B2, B, A3, A33, S, P, D;
     326      416408 :   if (pi)
     327             :   {
     328      416409 :     t = Fl_mul_pre(a, p3, p, pi);
     329      416412 :     t2 = Fl_sqr_pre(t, p, pi);
     330      416409 :     A = Fl_sub(b, Fl_triple(t2, p), p);
     331      416411 :     B = Fl_sub(c, Fl_mul_pre(t, Fl_add(A, t2, p), p, pi), p);
     332      416416 :     A3 =  Fl_mul_pre(A, p3, p, pi);
     333      416416 :     B2 = Fl_sqr_pre(B, p, pi);
     334             :   }
     335             :   else
     336             :   {
     337           0 :     t = Fl_mul(a, p3, p);
     338           0 :     t2 = Fl_sqr(t, p);
     339           0 :     A = Fl_sub(b, Fl_triple(t2, p), p);
     340           0 :     B = Fl_sub(c, Fl_mul(t, Fl_add(A, t2, p), p), p);
     341           0 :     A3 =  Fl_mul(A, p3, p);
     342           0 :     B2 = Fl_sqr(B, p);
     343             :   }
     344      416414 :   A33 = Fl_powu_pre(A3, 3, p, pi);
     345      416405 :   D = Fl_add(B2, Fl_double(Fl_double(A33, p), p), p);
     346      416403 :   S = Fl_neg(B,p);
     347      416403 :   P = Fl_neg(A3,p);
     348      416401 :   if (krouu(D,p) >= 0)
     349             :   {
     350      341676 :     ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
     351      341679 :     ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
     352      341682 :     if (p%3==2) /* 1 solutions */
     353      138581 :       vS1 = Fl_powu_pre(S1, p - p3, p, pi);
     354             :     else
     355             :     {
     356      203101 :       vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
     357      203103 :       if (vS1==~0UL) return p; /*0 solutions*/
     358             :       /*3 solutions*/
     359             :     }
     360      222622 :     if (!P) return Fl_sub(vS1, t, p);
     361       94781 :     vS2 = pi? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): Fl_div(P, vS1, p);
     362       94782 :     return Fl_sub(Fl_add(vS1,vS2, p), t, p);
     363             :   }
     364             :   else
     365             :   {
     366       74735 :     pari_sp av = avma;
     367       74735 :     GEN S1 = mkvecsmall2(Fl_halve(S, p), (p + 1UL) >> 1);
     368       74736 :     GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
     369             :     ulong Sa;
     370       74736 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     371       74736 :     Sa = vS1[1];
     372       74736 :     if (p%3==1) /*1 solutions*/
     373             :     {
     374       29300 :       ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
     375       29300 :       if (Fa!=P) Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
     376             :     }
     377       74736 :     set_avma(av);
     378       74736 :     return Fl_sub(Fl_double(Sa,p),t,p);
     379             :   }
     380             : }
     381             : 
     382             : /* Assume f is normalized */
     383             : static GEN
     384         119 : FpX_cubic_root(GEN ff, GEN p)
     385             : {
     386         119 :   GEN f = FpX_normalize(ff,p);
     387         119 :   GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
     388         119 :   ulong pm3 = umodiu(p,3);
     389          28 :   GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
     390         119 :                   : diviuexact(addiu(p,1),3);
     391         119 :   GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
     392         119 :   GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
     393         119 :   GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
     394         119 :   GEN A3 =  Fp_mul(A, p3, p), A33 = Fp_powu(A3, 3, p);
     395         119 :   GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
     396         119 :   GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
     397         119 :   if (kronecker(D,p) >= 0)
     398             :   {
     399          28 :     GEN s = Fp_sqrt(D, p), vS1, vS2;
     400          28 :     GEN S1 = S==s ? S: Fp_halve(Fp_sub(S, s, p), p);
     401          28 :     if (pm3 == 2) /* 1 solutions */
     402           0 :       vS1 = Fp_pow(S1, diviuexact(addiu(shifti(p, 1), 1), 3), p);
     403             :     else
     404             :     {
     405          28 :       vS1 = Fp_sqrtn(S1, utoi(3), p, NULL);
     406          28 :       if (!vS1) return p; /*0 solutions*/
     407             :       /*3 solutions*/
     408             :     }
     409          28 :     vS2 = P? Fp_mul(P, Fp_inv(vS1, p), p): 0;
     410          28 :     return Fp_sub(Fp_add(vS1,vS2, p), t, p);
     411             :   }
     412             :   else
     413             :   {
     414          91 :     pari_sp av = avma;
     415          91 :     GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
     416          91 :     GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
     417          91 :     GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
     418             :     GEN Sa;
     419          91 :     if (!vS1) return p; /*0 solutions, p%3==2*/
     420          91 :     Sa = gel(vS1,2);
     421          91 :     if (pm3 == 1) /*1 solutions*/
     422             :     {
     423           0 :       GEN Fa = FpXQ_norm(vS1, T, p);
     424           0 :       if (!equalii(Fa,P))
     425           0 :         Sa = Fp_mul(Sa, Fp_div(Fa, P, p),p);
     426             :     }
     427          91 :     set_avma(av);
     428          91 :     return Fp_sub(shifti(Sa,1),t,p);
     429             :   }
     430             : }
     431             : 
     432             : /* assume p > 2 prime; if fl is set, assume that f splits mod p */
     433             : static ulong
     434     4149855 : Flx_oneroot_pre_i(GEN f, ulong p, ulong pi, long fl)
     435             : {
     436             :   GEN pol, a;
     437             :   ulong q, PI;
     438             :   long da;
     439             : 
     440     4149855 :   if (Flx_val(f)) return 0;
     441     4148194 :   da = degpol(f); f = Flx_normalize(f, p);
     442     4147885 :   if (da == 1) return Fl_neg(f[2], p);
     443     4134292 :   PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
     444     4135156 :   switch(da)
     445             :   {
     446     3391839 :     case 2: return Flx_quad_root(f, p, PI, 1);
     447      402384 :     case 3: if (p>3) return Flx_cubic_root(f, p, PI); /*FALL THROUGH*/
     448             :   }
     449      347611 :   if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
     450      347611 :   if (!fl)
     451             :   {
     452      310295 :     a = Flxq_powu_pre(polx_Flx(f[1]), p - 1, f,p,pi);
     453      310189 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
     454      310189 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
     455      310199 :     a = Flx_gcd_pre(f,a, p, pi);
     456       37316 :   } else a = f;
     457      347574 :   da = degpol(a);
     458      347571 :   if (!da) return p;
     459      247926 :   a = Flx_normalize(a,p);
     460             : 
     461      247955 :   q = p >> 1;
     462      247955 :   pol = polx_Flx(f[1]);
     463      364509 :   for(pol[2] = 1;; pol[2]++)
     464             :   {
     465      364509 :     if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
     466      364523 :     switch(da)
     467             :     {
     468      155862 :       case 1: return Fl_neg(a[2], p);
     469       71405 :       case 2: return Flx_quad_root(a, p, PI, 0);
     470       20704 :       case 3: if (p>3) return Flx_cubic_root(a, p, PI); /*FALL THROUGH*/
     471             :       default: {
     472      116552 :         GEN b = Flxq_powu_pre(pol,q, a,p,pi);
     473             :         long db;
     474      116544 :         if (degpol(b) <= 0) continue;
     475      110785 :         b = Flx_gcd_pre(a,Flx_Fl_add(b,p-1,p), p, pi);
     476      110784 :         db = degpol(b); if (!db) continue;
     477      110784 :         b = Flx_normalize(b, p);
     478      110790 :         if (db <= (da >> 1)) {
     479       68298 :           a = b;
     480       68298 :           da = db;
     481             :         } else {
     482       42492 :           a = Flx_div_pre(a,b, p, pi);
     483       42496 :           da -= db;
     484             :         }
     485             :       }
     486             :     }
     487             :   }
     488             : }
     489             : ulong
     490     4111019 : Flx_oneroot_pre(GEN f, ulong p, ulong pi)
     491     4111019 : { return Flx_oneroot_pre_i(f, p, pi, 0); }
     492             : ulong
     493       38779 : Flx_oneroot_split_pre(GEN f, ulong p, ulong pi)
     494       38779 : { return Flx_oneroot_pre_i(f, p, pi, 1); }
     495             : 
     496             : /* assume p > 3 prime */
     497             : static GEN
     498        5169 : FpX_oneroot_i(GEN f, GEN p)
     499             : {
     500             :   GEN pol, pol0, a, q;
     501             :   long da;
     502             : 
     503        5169 :   if (ZX_val(f)) return gen_0;
     504        4847 :   f = FpX_normalize(f, p);
     505        4847 :   switch(degpol(f))
     506             :   {
     507         821 :     case 1: return subii(p, gel(f,2));
     508        3837 :     case 2: return FpX_quad_root(f, p, 1);
     509         119 :     case 3: return FpX_cubic_root(f, p);
     510             :   }
     511             : 
     512          70 :   a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
     513          70 :   if (lg(a) < 3) pari_err_PRIME("rootmod",p);
     514          70 :   a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
     515          70 :   a = FpX_gcd(f,a, p);
     516          70 :   da = degpol(a);
     517          70 :   if (!da) return NULL;
     518          70 :   a = FpX_normalize(a,p);
     519             : 
     520          70 :   q = shifti(p,-1);
     521          70 :   pol0 = icopy(gen_1); /* constant term, will vary in place */
     522          70 :   pol = deg1pol_shallow(gen_1, pol0, varn(f));
     523         224 :   for (pol0[2]=1; ; pol0[2]++)
     524             :   {
     525         224 :     if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
     526         224 :     switch(da)
     527             :     {
     528          42 :       case 1: return subii(p, gel(a,2));
     529          28 :       case 2: return FpX_quad_root(a, p, 0);
     530         154 :       default: {
     531         154 :         GEN b = FpXQ_pow(pol,q, a,p);
     532             :         long db;
     533         154 :         if (degpol(b) <= 0) continue;
     534         147 :         b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
     535         147 :         db = degpol(b); if (!db) continue;
     536         147 :         b = FpX_normalize(b, p);
     537         147 :         if (db <= (da >> 1)) {
     538         105 :           a = b;
     539         105 :           da = db;
     540             :         } else {
     541          42 :           a = FpX_div(a,b, p);
     542          42 :           da -= db;
     543             :         }
     544             :       }
     545             :     }
     546             :   }
     547             : }
     548             : 
     549             : ulong
     550        2492 : Flx_oneroot(GEN f, ulong p)
     551             : {
     552        2492 :   pari_sp av = avma;
     553        2492 :   switch(lg(f))
     554             :   {
     555           0 :     case 2: return 0;
     556           0 :     case 3: return p;
     557             :   }
     558        2492 :   if (p == 2) return Flx_oneroot_mod_2(f);
     559        2492 :   return gc_ulong(av, Flx_oneroot_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)));
     560             : }
     561             : 
     562             : ulong
     563          14 : Flx_oneroot_split(GEN f, ulong p)
     564             : {
     565          14 :   pari_sp av = avma;
     566          14 :   switch(lg(f))
     567             :   {
     568           0 :     case 2: return 0;
     569           0 :     case 3: return p;
     570             :   }
     571          14 :   if (p == 2) return Flx_oneroot_mod_2(f);
     572          14 :   return gc_ulong(av, Flx_oneroot_split_pre(f, p, 0));
     573             : }
     574             : 
     575             : /* assume that p is prime */
     576             : GEN
     577      233884 : FpX_oneroot(GEN f, GEN p)
     578             : {
     579      233884 :   pari_sp av = avma;
     580      233884 :   f = ZX_rootmod_init(f, p);
     581      233884 :   switch(lg(f))
     582             :   {
     583           0 :     case 2: set_avma(av); return gen_0;
     584           0 :     case 3: return gc_NULL(av);
     585             :   }
     586      233884 :   if (typ(f) == t_VECSMALL)
     587             :   {
     588      228715 :     ulong r, pp = p[2];
     589      228715 :     if (pp == 2)
     590          91 :       r = Flx_oneroot_mod_2(f);
     591             :     else
     592      228624 :       r = Flx_oneroot_pre(f, pp, SMALL_ULONG(pp)? 0: get_Fl_red(pp));
     593      228715 :     set_avma(av);
     594      228715 :     return (r == pp)? NULL: utoi(r);
     595             :   }
     596        5169 :   f = FpX_oneroot_i(f, p);
     597        5169 :   if (!f) return gc_NULL(av);
     598        5169 :   return gerepileuptoint(av, f);
     599             : }
     600             : 
     601             : /* returns a root of unity in F_p that is suitable for finding a factor   */
     602             : /* of degree deg_factor of a polynomial of degree deg; the order is       */
     603             : /* returned in n                                                          */
     604             : /* A good choice seems to be n close to deg/deg_factor; we choose n       */
     605             : /* twice as big and decrement until it divides p-1.                       */
     606             : static GEN
     607         154 : good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
     608             : {
     609         154 :    pari_sp ltop = avma;
     610             :    GEN pm, factn, power, base, zeta;
     611             :    long n;
     612             : 
     613         154 :    pm = subis (p, 1ul);
     614         336 :    for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
     615         154 :    factn = Z_factor(stoi(n));
     616         154 :    power = diviuexact (pm, n);
     617         154 :    base = gen_1;
     618             :    do {
     619         259 :       base = addis (base, 1l);
     620         259 :       zeta = Fp_pow (base, power, p);
     621             :    }
     622         259 :    while (!equaliu (Fp_order (zeta, factn, p), n));
     623         154 :    *pt_n = n;
     624         154 :    return gerepileuptoint (ltop, zeta);
     625             : }
     626             : 
     627             : GEN
     628        1092 : FpX_oneroot_split(GEN fact, GEN p)
     629             : {
     630        1092 :   pari_sp av = avma;
     631             :   long n, deg_f, i, dmin;
     632             :   GEN prim, expo, minfactor, xplusa, zeta, xpow;
     633        1092 :   fact = FpX_normalize(fact, p);
     634        1092 :   deg_f = degpol(fact);
     635        1092 :   if (deg_f <= 3) return FpX_oneroot(fact, p);
     636         133 :   minfactor = fact; /* factor of minimal degree found so far */
     637         133 :   dmin = degpol(minfactor);
     638         133 :   xplusa = pol_x(varn(fact));
     639         287 :   while (dmin > 3)
     640             :   {
     641             :     /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
     642             :     /* zeta varies over the roots of unity in F_p                         */
     643         154 :     fact = minfactor; deg_f = dmin;
     644         154 :     zeta = gen_1;
     645         154 :     prim = good_root_of_unity(p, deg_f, 1, &n);
     646         154 :     expo = diviuexact(subiu(p, 1), n);
     647             :     /* update X+a, avoid a=0 */
     648         154 :     gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
     649         154 :     xpow = FpXQ_pow (xplusa, expo, fact, p);
     650         301 :     for (i = 0; i < n; i++)
     651             :     {
     652         231 :       GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
     653         231 :       long dtmp = degpol(tmp);
     654         231 :       if (dtmp > 0 && dtmp < deg_f)
     655             :       {
     656         154 :         fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
     657         154 :         if (dtmp < dmin)
     658             :         {
     659         154 :           minfactor = FpX_normalize (tmp, p);
     660         154 :           dmin = dtmp;
     661         154 :           if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
     662             :             /* stop early to avoid too many gcds */
     663             :             break;
     664             :         }
     665             :       }
     666         147 :       zeta = Fp_mul (zeta, prim, p);
     667             :     }
     668             :   }
     669         133 :   return gerepileuptoint(av, FpX_oneroot(minfactor, p));
     670             : }
     671             : 
     672             : /*******************************************************************/
     673             : /*                                                                 */
     674             : /*                     FACTORISATION MODULO p                      */
     675             : /*                                                                 */
     676             : /*******************************************************************/
     677             : 
     678             : /* F / E  a vector of vectors of factors / exponents of virtual length l
     679             :  * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
     680             : static GEN
     681     1882598 : FE_concat(GEN F, GEN E, long l)
     682             : {
     683     1882598 :   setlg(E,l); E = shallowconcat1(E);
     684     1882598 :   setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
     685             : }
     686             : 
     687             : static GEN
     688          14 : ddf_to_ddf2_i(GEN V, long fl)
     689             : {
     690             :   GEN F, D;
     691          14 :   long i, j, l = lg(V);
     692          14 :   F = cgetg(l, t_VEC);
     693          14 :   D = cgetg(l, t_VECSMALL);
     694         112 :   for (i = j = 1; i < l; i++)
     695             :   {
     696          98 :     GEN Vi = gel(V,i);
     697          98 :     if ((fl==2 && F2x_degree(Vi) == 0)
     698          98 :       ||(fl==0 && degpol(Vi) == 0)) continue;
     699          35 :     gel(F,j) = Vi;
     700          35 :     uel(D,j) = i; j++;
     701             :   }
     702          14 :   setlg(F,j);
     703          14 :   setlg(D,j); return mkvec2(F,D);
     704             : }
     705             : 
     706             : GEN
     707           7 : ddf_to_ddf2(GEN V)
     708           7 : { return ddf_to_ddf2_i(V, 0); }
     709             : 
     710             : static GEN
     711           7 : F2x_ddf_to_ddf2(GEN V)
     712           7 : { return ddf_to_ddf2_i(V, 2); }
     713             : 
     714             : GEN
     715     5318021 : vddf_to_simplefact(GEN V, long d)
     716             : {
     717             :   GEN E, F;
     718     5318021 :   long i, j, c, l = lg(V);
     719     5318021 :   F = cgetg(d+1, t_VECSMALL);
     720     5317174 :   E = cgetg(d+1, t_VECSMALL);
     721    10716859 :   for (i = c = 1; i < l; i++)
     722             :   {
     723     5399916 :     GEN Vi = gel(V,i);
     724     5399916 :     long l = lg(Vi);
     725    27225185 :     for (j = 1; j < l; j++)
     726             :     {
     727    21825339 :       long k, n = degpol(gel(Vi,j)) / j;
     728    33049752 :       for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
     729             :     }
     730             :   }
     731     5316943 :   setlg(F,c);
     732     5316965 :   setlg(E,c);
     733     5317100 :   return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
     734             : }
     735             : 
     736             : /* product of terms of degree 1 in factorization of f */
     737             : GEN
     738      234191 : FpX_split_part(GEN f, GEN p)
     739             : {
     740      234191 :   long n = degpol(f);
     741      234191 :   GEN z, X = pol_x(varn(f));
     742      234191 :   if (n <= 1) return f;
     743      232223 :   f = FpX_red(f, p);
     744      232224 :   z = FpX_sub(FpX_Frobenius(f, p), X, p);
     745      232224 :   return FpX_gcd(z,f,p);
     746             : }
     747             : 
     748             : /* Compute the number of roots in Fp without counting multiplicity
     749             :  * return -1 for 0 polynomial. lc(f) must be prime to p. */
     750             : long
     751      138146 : FpX_nbroots(GEN f, GEN p)
     752             : {
     753      138146 :   pari_sp av = avma;
     754      138146 :   GEN z = FpX_split_part(f, p);
     755      138146 :   return gc_long(av, degpol(z));
     756             : }
     757             : 
     758             : /* 1 < deg(f) <= p */
     759             : static int
     760       81347 : Flx_is_totally_split_i(GEN f, ulong p)
     761             : {
     762       81347 :   GEN F = Flx_Frobenius(f, p);
     763       81347 :   return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
     764             : }
     765             : int
     766       81354 : Flx_is_totally_split(GEN f, ulong p)
     767             : {
     768       81354 :   pari_sp av = avma;
     769       81354 :   ulong n = degpol(f);
     770       81354 :   if (n <= 1) return 1;
     771       81347 :   if (n > p) return 0; /* includes n < 0 */
     772       81347 :   return gc_bool(av, Flx_is_totally_split_i(f,p));
     773             : }
     774             : int
     775           0 : FpX_is_totally_split(GEN f, GEN p)
     776             : {
     777           0 :   pari_sp av = avma;
     778           0 :   ulong n = degpol(f);
     779             :   int u;
     780           0 :   if (n <= 1) return 1;
     781           0 :   if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
     782           0 :   if (lgefint(p) != 3)
     783           0 :     u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
     784             :   else
     785             :   {
     786           0 :     ulong pp = (ulong)p[2];
     787           0 :     u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
     788             :   }
     789           0 :   return gc_bool(av, u);
     790             : }
     791             : 
     792             : long
     793     4388585 : Flx_nbroots(GEN f, ulong p)
     794             : {
     795     4388585 :   long n = degpol(f);
     796             :   ulong pi;
     797     4388585 :   pari_sp av = avma;
     798             :   GEN z;
     799     4388585 :   if (n <= 1) return n;
     800     4374858 :   if (n == 2)
     801             :   {
     802             :     ulong D;
     803      578478 :     if (p==2) return (f[2]==0) + (f[2]!=f[3]);
     804      433737 :     D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
     805      433742 :     return 1 + krouu(D,p);
     806             :   }
     807     3796380 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     808     3796380 :   z = Flx_sub(Flx_Frobenius_pre(f, p, pi), polx_Flx(f[1]), p);
     809     3796372 :   z = Flx_gcd_pre(z, f, p, pi);
     810     3796383 :   return gc_long(av, degpol(z));
     811             : }
     812             : 
     813             : long
     814        4256 : FpX_ddf_degree(GEN T, GEN XP, GEN p)
     815             : {
     816        4256 :   pari_sp av = avma;
     817             :   GEN X, b, g, xq;
     818             :   long i, j, n, v, B, l, m;
     819             :   pari_timer ti;
     820             :   hashtable h;
     821             : 
     822        4256 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     823        4256 :   X = pol_x(v);
     824        4256 :   if (ZX_equal(X,XP)) return 1;
     825        4256 :   B = n/2;
     826        4256 :   l = usqrt(B);
     827        4256 :   m = (B+l-1)/l;
     828        4256 :   T = FpX_get_red(T, p);
     829        4256 :   hash_init_GEN(&h, l+2, ZX_equal, 1);
     830        4256 :   hash_insert_long(&h, X,  0);
     831        4256 :   hash_insert_long(&h, XP, 1);
     832        4256 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     833        4256 :   b = XP;
     834        4256 :   xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
     835        4256 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
     836       10178 :   for (i = 3; i <= l+1; i++)
     837             :   {
     838        6601 :     b = FpX_FpXQV_eval(b, xq, T, p);
     839        6601 :     if (gequalX(b)) return gc_long(av,i-1);
     840        5922 :     hash_insert_long(&h, b, i-1);
     841             :   }
     842        3577 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
     843        3577 :   g = b;
     844        3577 :   xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
     845        3577 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
     846       12208 :   for(i = 2; i <= m+1; i++)
     847             :   {
     848       10619 :     g = FpX_FpXQV_eval(g, xq, T, p);
     849       10619 :     if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
     850             :   }
     851        1589 :   return gc_long(av,n);
     852             : }
     853             : 
     854             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
     855             : static GEN
     856        1021 : FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
     857             : {
     858             :   GEN b, g, h, F, f, Tr, xq;
     859             :   long i, j, n, v, B, l, m;
     860             :   pari_timer ti;
     861             : 
     862        1021 :   n = get_FpX_degree(T); v = get_FpX_var(T);
     863        1021 :   if (n == 0) return cgetg(1, t_VEC);
     864        1021 :   if (n == 1) return mkvec(get_FpX_mod(T));
     865         848 :   B = n/2;
     866         848 :   l = usqrt(B);
     867         848 :   m = (B+l-1)/l;
     868         848 :   T = FpX_get_red(T, p);
     869         848 :   b = cgetg(l+2, t_VEC);
     870         848 :   gel(b, 1) = pol_x(v);
     871         848 :   gel(b, 2) = XP;
     872         848 :   if (DEBUGLEVEL>=7) timer_start(&ti);
     873         848 :   xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
     874         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
     875        1077 :   for (i = 3; i <= l+1; i++)
     876         229 :     gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
     877         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
     878         848 :   xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
     879         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
     880         848 :   g = cgetg(m+1, t_VEC);
     881         848 :   gel(g, 1) = gel(xq, 2);
     882        1843 :   for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
     883         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
     884         848 :   h = cgetg(m+1, t_VEC);
     885        2691 :   for (j = 1; j <= m; j++)
     886             :   {
     887        1843 :     pari_sp av = avma;
     888        1843 :     GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
     889        3076 :     for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
     890        1843 :     gel(h,j) = gerepileupto(av, e);
     891             :   }
     892         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
     893         848 :   Tr = get_FpX_mod(T);
     894         848 :   F = cgetg(m+1, t_VEC);
     895        2691 :   for (j = 1; j <= m; j++)
     896             :   {
     897        1843 :     GEN u = FpX_gcd(Tr, gel(h,j), p);
     898        1843 :     if (degpol(u))
     899             :     {
     900         481 :       u = FpX_normalize(u, p);
     901         481 :       Tr = FpX_div(Tr, u, p);
     902             :     }
     903        1843 :     gel(F,j) = u;
     904             :   }
     905         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
     906         848 :   f = const_vec(n, pol_1(v));
     907        2691 :   for (j = 1; j <= m; j++)
     908             :   {
     909        1843 :     GEN e = gel(F, j);
     910        1915 :     for (i=l-1; i >= 0; i--)
     911             :     {
     912        1915 :       GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
     913        1915 :       if (degpol(u))
     914             :       {
     915         511 :         u = FpX_normalize(u, p);
     916         511 :         gel(f, l*j-i) = u;
     917         511 :         e = FpX_div(e, u, p);
     918             :       }
     919        1915 :       if (!degpol(e)) break;
     920             :     }
     921             :   }
     922         848 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
     923         848 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
     924         848 :   return f;
     925             : }
     926             : 
     927             : static void
     928           0 : FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     929             : {
     930           0 :   long n = degpol(Tp), r = n/d, ct = 0;
     931             :   GEN T, f, ff, p2;
     932           0 :   if (r==1) { gel(V, idx) = Tp; return; }
     933           0 :   p2 = shifti(p,-1);
     934           0 :   T = FpX_get_red(Tp, p);
     935           0 :   XP = FpX_rem(XP, T, p);
     936             :   while (1)
     937           0 :   {
     938           0 :     pari_sp btop = avma;
     939             :     long i;
     940           0 :     GEN g = random_FpX(n, varn(Tp), p);
     941           0 :     GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
     942           0 :     if (signe(t) == 0) continue;
     943           0 :     for(i=1; i<=10; i++)
     944             :     {
     945           0 :       pari_sp btop2 = avma;
     946           0 :       GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
     947           0 :       f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
     948           0 :       if (degpol(f) > 0 && degpol(f) < n) break;
     949           0 :       set_avma(btop2);
     950             :     }
     951           0 :     if (degpol(f) > 0 && degpol(f) < n) break;
     952           0 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
     953           0 :     set_avma(btop);
     954             :   }
     955           0 :   f = FpX_normalize(f, p);
     956           0 :   ff = FpX_div(Tp, f ,p);
     957           0 :   FpX_edf_simple(f, XP, d, p, V, idx);
     958           0 :   FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
     959             : }
     960             : 
     961             : static void
     962        1075 : FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
     963             : {
     964             :   pari_sp av;
     965        1075 :   GEN Tp = get_FpX_mod(T);
     966        1075 :   long n = degpol(hp), vT = varn(Tp), ct = 0;
     967             :   GEN u1, u2, f1, f2, R, h;
     968        1075 :   h = FpX_get_red(hp, p);
     969        1075 :   t = FpX_rem(t, T, p);
     970        1075 :   av = avma;
     971             :   do
     972             :   {
     973        1660 :     set_avma(av);
     974        1660 :     R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
     975        1660 :     u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
     976        1660 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
     977        1660 :   } while (degpol(u1)==0 || degpol(u1)==n);
     978        1075 :   f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
     979        1075 :   f1 = FpX_normalize(f1, p);
     980        1075 :   u2 = FpX_div(hp, u1, p);
     981        1075 :   f2 = FpX_div(Tp, f1, p);
     982        1075 :   if (degpol(u1)==1)
     983         784 :     gel(V, idx) = f1;
     984             :   else
     985         291 :     FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
     986        1075 :   idx += degpol(f1)/d;
     987        1075 :   if (degpol(u2)==1)
     988         728 :     gel(V, idx) = f2;
     989             :   else
     990         347 :     FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
     991        1075 : }
     992             : 
     993             : /* assume Tp a squarefree product of r > 1 irred. factors of degree d */
     994             : static void
     995         437 : FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
     996             : {
     997         437 :   long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
     998             :   GEN T, h, t;
     999             :   pari_timer ti;
    1000             : 
    1001         437 :   T = FpX_get_red(Tp, p);
    1002         437 :   XP = FpX_rem(XP, T, p);
    1003         437 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    1004             :   do
    1005             :   {
    1006         437 :     GEN g = random_FpX(n, vT, p);
    1007         437 :     t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
    1008         437 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
    1009         437 :     h = FpXQ_minpoly(t, T, p);
    1010         437 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
    1011         437 :     if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
    1012         437 :   } while (degpol(h) != r);
    1013         437 :   FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
    1014         437 : }
    1015             : 
    1016             : static GEN
    1017        1007 : FpX_factor_Shoup(GEN T, GEN p)
    1018             : {
    1019        1007 :   long i, n, s = 0;
    1020             :   GEN XP, D, V;
    1021        1007 :   long e = expi(p);
    1022             :   pari_timer ti;
    1023        1007 :   n = get_FpX_degree(T);
    1024        1007 :   T = FpX_get_red(T, p);
    1025        1007 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1026        1007 :   XP = FpX_Frobenius(T, p);
    1027        1007 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1028        1007 :   D = FpX_ddf_Shoup(T, XP, p);
    1029        1007 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1030        1007 :   s = ddf_to_nbfact(D);
    1031        1007 :   V = cgetg(s+1, t_COL);
    1032        7292 :   for (i = 1, s = 1; i <= n; i++)
    1033             :   {
    1034        6285 :     GEN Di = gel(D,i);
    1035        6285 :     long ni = degpol(Di), ri = ni/i;
    1036        6285 :     if (ni == 0) continue;
    1037        1125 :     Di = FpX_normalize(Di, p);
    1038        1125 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1039         437 :     if (ri <= e*expu(e))
    1040         437 :       FpX_edf(Di, XP, i, p, V, s);
    1041             :     else
    1042           0 :       FpX_edf_simple(Di, XP, i, p, V, s);
    1043         437 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
    1044         437 :     s += ri;
    1045             :   }
    1046        1007 :   return V;
    1047             : }
    1048             : 
    1049             : long
    1050     2315108 : ddf_to_nbfact(GEN D)
    1051             : {
    1052     2315108 :   long l = lg(D), i, s = 0;
    1053    14646570 :   for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
    1054     2315102 :   return s;
    1055             : }
    1056             : 
    1057             : /* Yun algorithm: Assume p > degpol(T) */
    1058             : static GEN
    1059        1696 : FpX_factor_Yun(GEN T, GEN p)
    1060             : {
    1061        1696 :   long n = degpol(T), i = 1;
    1062        1696 :   GEN a, b, c, d = FpX_deriv(T, p);
    1063        1696 :   GEN V = cgetg(n+1,t_VEC);
    1064        1696 :   a = FpX_gcd(T, d, p);
    1065        1696 :   if (degpol(a) == 0) return mkvec(T);
    1066         584 :   b = FpX_div(T, a, p);
    1067             :   do
    1068             :   {
    1069        2638 :     c = FpX_div(d, a, p);
    1070        2638 :     d = FpX_sub(c, FpX_deriv(b, p), p);
    1071        2638 :     a = FpX_normalize(FpX_gcd(b, d, p), p);
    1072        2638 :     gel(V, i++) = a;
    1073        2638 :     b = FpX_div(b, a, p);
    1074        2638 :   } while (degpol(b));
    1075         584 :   setlg(V, i); return V;
    1076             : }
    1077             : GEN
    1078      341852 : FpX_factor_squarefree(GEN T, GEN p)
    1079             : {
    1080      341852 :   if (lgefint(p)==3)
    1081             :   {
    1082      341118 :     ulong pp = (ulong)p[2];
    1083      341118 :     GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
    1084      341118 :     return FlxV_to_ZXV(u);
    1085             :   }
    1086         734 :   return FpX_factor_Yun(T, p);
    1087             : }
    1088             : 
    1089             : GEN
    1090      221697 : FpX_roots_mult(GEN T, long n, GEN p)
    1091             : {
    1092      221697 :   pari_sp av = avma;
    1093      221697 :   GEN V = FpX_factor_squarefree(T, p), W;
    1094      221697 :   long l = lg(V), i;
    1095      221697 :   if (l <= n) { set_avma(av); return cgetg(1,t_COL); }
    1096       36776 :   W = cgetg(l-n+1,t_VEC);
    1097      120325 :   for (i = n; i < l; i++)
    1098       83549 :     gel(W,i-n+1) = FpX_roots(gel(V,i), p);
    1099       36776 :   return gerepileupto(av, sort(shallowconcat1(W)));
    1100             : }
    1101             : 
    1102             : long
    1103         168 : FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
    1104             : {
    1105         168 :   pari_sp av = avma;
    1106             :   GEN lc, F;
    1107         168 :   long i, l, n = degpol(f), v = varn(f);
    1108         168 :   if (n % k) return 0;
    1109         168 :   if (lgefint(p)==3)
    1110             :   {
    1111         126 :     ulong pp = p[2];
    1112         126 :     GEN fp = ZX_to_Flx(f, pp);
    1113         126 :     if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
    1114         105 :     if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
    1115         105 :     return 1;
    1116             :   }
    1117          42 :   lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
    1118          42 :   if (!lc) { av = avma; return 0; }
    1119          42 :   F = FpX_factor_Yun(f, p); l = lg(F)-1;
    1120        1491 :   for(i=1; i <= l; i++)
    1121        1456 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    1122          35 :   if (pt_r)
    1123             :   {
    1124          35 :     GEN r = scalarpol(lc, v), s = pol_1(v);
    1125        1484 :     for (i=l; i>=1; i--)
    1126             :     {
    1127        1449 :       if (i%k) continue;
    1128         294 :       s = FpX_mul(s, gel(F,i), p);
    1129         294 :       r = FpX_mul(r, s, p);
    1130             :     }
    1131          35 :     *pt_r = gerepileupto(av, r);
    1132           0 :   } else av = avma;
    1133          35 :   return 1;
    1134             : }
    1135             : 
    1136             : static GEN
    1137         906 : FpX_factor_Cantor(GEN T, GEN p)
    1138             : {
    1139         906 :   GEN E, F, V = FpX_factor_Yun(T, p);
    1140         906 :   long i, j, l = lg(V);
    1141         906 :   F = cgetg(l, t_VEC);
    1142         906 :   E = cgetg(l, t_VEC);
    1143        2426 :   for (i=1, j=1; i < l; i++)
    1144        1520 :     if (degpol(gel(V,i)))
    1145             :     {
    1146        1007 :       GEN Fj = FpX_factor_Shoup(gel(V,i), p);
    1147        1007 :       gel(F, j) = Fj;
    1148        1007 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1149        1007 :       j++;
    1150             :     }
    1151         906 :   return sort_factor_pol(FE_concat(F,E,j), cmpii);
    1152             : }
    1153             : 
    1154             : static GEN
    1155           0 : FpX_ddf_i(GEN T, GEN p)
    1156             : {
    1157             :   GEN XP;
    1158           0 :   T = FpX_get_red(T, p);
    1159           0 :   XP = FpX_Frobenius(T, p);
    1160           0 :   return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
    1161             : }
    1162             : 
    1163             : GEN
    1164           7 : FpX_ddf(GEN f, GEN p)
    1165             : {
    1166           7 :   pari_sp av = avma;
    1167             :   GEN F;
    1168           7 :   switch(ZX_factmod_init(&f, p))
    1169             :   {
    1170           7 :     case 0:  F = F2x_ddf(f);
    1171           7 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    1172           0 :     case 1:  F = Flx_ddf(f,p[2]);
    1173           0 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    1174           0 :     default: F = FpX_ddf_i(f,p); break;
    1175             :   }
    1176           7 :   return gerepilecopy(av, F);
    1177             : }
    1178             : 
    1179             : static GEN Flx_simplefact_Cantor(GEN T, ulong p);
    1180             : static GEN
    1181          14 : FpX_simplefact_Cantor(GEN T, GEN p)
    1182             : {
    1183             :   GEN V;
    1184             :   long i, l;
    1185          14 :   if (lgefint(p) == 3)
    1186             :   {
    1187           0 :     ulong pp = p[2];
    1188           0 :     return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
    1189             :   }
    1190          14 :   T = FpX_get_red(T, p);
    1191          14 :   V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
    1192          28 :   for (i=1; i < l; i++)
    1193          14 :     gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
    1194          14 :   return vddf_to_simplefact(V, get_FpX_degree(T));
    1195             : }
    1196             : 
    1197             : static int
    1198           0 : FpX_isirred_Cantor(GEN Tp, GEN p)
    1199             : {
    1200           0 :   pari_sp av = avma;
    1201             :   pari_timer ti;
    1202             :   long n;
    1203           0 :   GEN T = get_FpX_mod(Tp);
    1204           0 :   GEN dT = FpX_deriv(T, p);
    1205             :   GEN XP, D;
    1206           0 :   if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
    1207           0 :   n = get_FpX_degree(T);
    1208           0 :   T = FpX_get_red(Tp, p);
    1209           0 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1210           0 :   XP = FpX_Frobenius(T, p);
    1211           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
    1212           0 :   D = FpX_ddf_Shoup(T, XP, p);
    1213           0 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
    1214           0 :   return gc_bool(av, degpol(gel(D,n)) == n);
    1215             : }
    1216             : 
    1217             : static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
    1218             : 
    1219             : /*Assume that p is large and odd*/
    1220             : static GEN
    1221        2508 : FpX_factor_i(GEN f, GEN pp, long flag)
    1222             : {
    1223        2508 :   long d = degpol(f);
    1224        2508 :   if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
    1225         920 :   switch(flag)
    1226             :   {
    1227         906 :     default: return FpX_factor_Cantor(f, pp);
    1228          14 :     case 1: return FpX_simplefact_Cantor(f, pp);
    1229           0 :     case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
    1230             :   }
    1231             : }
    1232             : 
    1233             : long
    1234           0 : FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
    1235             : {
    1236           0 :   pari_sp av = avma;
    1237           0 :   long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
    1238           0 :   return gc_long(av,s);
    1239             : }
    1240             : 
    1241             : long
    1242           0 : FpX_nbfact(GEN T, GEN p)
    1243             : {
    1244           0 :   pari_sp av = avma;
    1245           0 :   GEN XP = FpX_Frobenius(T, p);
    1246           0 :   long n = FpX_nbfact_Frobenius(T, XP, p);
    1247           0 :   return gc_long(av,n);
    1248             : }
    1249             : 
    1250             : /* p > 2 */
    1251             : static GEN
    1252           7 : FpX_is_irred_2(GEN f, GEN p, long d)
    1253             : {
    1254           7 :   switch(d)
    1255             :   {
    1256           0 :     case -1:
    1257           0 :     case 0: return NULL;
    1258           0 :     case 1: return gen_1;
    1259             :   }
    1260           7 :   return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
    1261             : }
    1262             : /* p > 2 */
    1263             : static GEN
    1264          14 : FpX_degfact_2(GEN f, GEN p, long d)
    1265             : {
    1266          14 :   switch(d)
    1267             :   {
    1268           0 :     case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
    1269           0 :     case 0: return trivial_fact();
    1270           0 :     case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
    1271             :   }
    1272          14 :   switch(FpX_quad_factortype(f, p)) {
    1273           7 :     case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1274           7 :     case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
    1275           0 :     default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
    1276             :   }
    1277             : }
    1278             : 
    1279             : GEN
    1280          70 : prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
    1281             : GEN
    1282      349933 : trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
    1283             : 
    1284             : /* not gerepile safe */
    1285             : static GEN
    1286        1567 : FpX_factor_2(GEN f, GEN p, long d)
    1287             : {
    1288             :   GEN r, s, R, S;
    1289             :   long v;
    1290             :   int sgn;
    1291        1567 :   switch(d)
    1292             :   {
    1293           7 :     case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
    1294          37 :     case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1295         674 :     case  1: retmkvec2(mkcol(f), mkvecsmall(1));
    1296             :   }
    1297         849 :   r = FpX_quad_root(f, p, 1);
    1298         849 :   if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
    1299         258 :   v = varn(f);
    1300         258 :   s = FpX_otherroot(f, r, p);
    1301         258 :   if (signe(r)) r = subii(p, r);
    1302         258 :   if (signe(s)) s = subii(p, s);
    1303         258 :   sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
    1304         258 :   R = deg1pol_shallow(gen_1, r, v);
    1305         258 :   if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
    1306         156 :   S = deg1pol_shallow(gen_1, s, v);
    1307         156 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1308             : }
    1309             : static GEN
    1310        1588 : FpX_factor_deg2(GEN f, GEN p, long d, long flag)
    1311             : {
    1312        1588 :   switch(flag) {
    1313           7 :     case 2: return FpX_is_irred_2(f, p, d);
    1314          14 :     case 1: return FpX_degfact_2(f, p, d);
    1315        1567 :     default: return FpX_factor_2(f, p, d);
    1316             :   }
    1317             : }
    1318             : 
    1319             : static int
    1320      451461 : F2x_quad_factortype(GEN x)
    1321      451461 : { return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
    1322             : 
    1323             : static GEN
    1324          14 : F2x_is_irred_2(GEN f, long d)
    1325          14 : { return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
    1326             : 
    1327             : static GEN
    1328       19676 : F2x_degfact_2(GEN f, long d)
    1329             : {
    1330       19676 :   if (!d) return trivial_fact();
    1331       19676 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1332       19480 :   switch(F2x_quad_factortype(f)) {
    1333        5964 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1334        6244 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1335        7272 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1336             :   }
    1337             : }
    1338             : 
    1339             : static GEN
    1340     1164901 : F2x_factor_2(GEN f, long d)
    1341             : {
    1342     1164901 :   long v = f[1];
    1343     1164901 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1344      909922 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1345      427181 :   switch(F2x_quad_factortype(f))
    1346             :   {
    1347       71705 :   case -1: return mkvec2(mkcol(f), mkvecsmall(1));
    1348      327983 :   case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
    1349       27516 :   default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
    1350             :   }
    1351             : }
    1352             : static GEN
    1353     1184591 : F2x_factor_deg2(GEN f, long d, long flag)
    1354             : {
    1355     1184591 :   switch(flag) {
    1356          14 :     case 2: return F2x_is_irred_2(f, d);
    1357       19676 :     case 1: return F2x_degfact_2(f, d);
    1358     1164901 :     default: return F2x_factor_2(f, d);
    1359             :   }
    1360             : }
    1361             : 
    1362             : /* xt = NULL or x^(p-1)/2 mod g */
    1363             : static void
    1364       30282 : split_squares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
    1365             : {
    1366       30282 :   ulong q = p >> 1;
    1367       30282 :   GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
    1368       30282 :   long d = degpol(a);
    1369       30282 :   if (d < 0)
    1370             :   {
    1371             :     ulong i;
    1372         407 :     split_add_done(S, (GEN)1);
    1373         407 :     if (!pi)
    1374        1555 :       for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
    1375             :     else
    1376           0 :       for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr_pre(i,p,pi));
    1377             :   } else {
    1378       29875 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1379       29875 :     if (d)
    1380             :     {
    1381       29786 :       if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
    1382       29786 :       a = Flx_gcd_pre(a, xt, p, pi);
    1383       29786 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1384             :     }
    1385             :   }
    1386       30282 : }
    1387             : static void
    1388       30282 : split_nonsquares(struct split_t *S, GEN g, ulong p, ulong pi, GEN xt)
    1389             : {
    1390       30282 :   ulong q = p >> 1;
    1391       30282 :   GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
    1392       30282 :   long d = degpol(a);
    1393       30282 :   if (d < 0)
    1394             :   {
    1395         249 :     ulong i, z = nonsquare_Fl(p);
    1396         249 :     split_add_done(S, (GEN)z);
    1397         249 :     if (!pi)
    1398         522 :       for (i = 2; i <= q; i++)
    1399         273 :         split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
    1400             :     else
    1401           0 :       for (i = 2; i <= q; i++)
    1402           0 :         split_add_done(S, (GEN)Fl_mul_pre(z, Fl_sqr_pre(i,p,pi), p,pi));
    1403             :   } else {
    1404       30033 :     if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
    1405       30033 :     if (d)
    1406             :     {
    1407       29676 :       if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
    1408       29676 :       a = Flx_gcd_pre(a, xt, p, pi);
    1409       29676 :       if (degpol(a)) split_add(S, Flx_normalize(a, p));
    1410             :     }
    1411             :   }
    1412       30282 : }
    1413             : /* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
    1414             :  * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
    1415             : static int
    1416     5066869 : split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p, ulong pi)
    1417             : {
    1418     5066869 :   GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
    1419     5066776 :   long d = degpol(g);
    1420     5066782 :   if (d < 0) return 0;
    1421     5066182 :   if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
    1422     5066173 :   if (!d) return 1;
    1423     5036201 :   if ((p >> 4) <= (ulong)d)
    1424             :   { /* small p; split directly using x^((p-1)/2) +/- 1 */
    1425       27600 :     GEN xt = ((ulong)d < (p>>1))? Flx_rem_pre(monomial_Flx(1, p>>1, g[1]), g, p, pi)
    1426       30282 :                                 : NULL;
    1427       30282 :     split_squares(S, g, p, pi, xt);
    1428       30282 :     split_nonsquares(S, g, p, pi, xt);
    1429             :   } else { /* large p; use x^(p-1) - 1 directly */
    1430     5005919 :     a = Flxq_powu_pre(polx_Flx(f[1]), p-1, g, p, pi);
    1431     4993066 :     if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
    1432     4993066 :     a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
    1433     4995597 :     g = Flx_gcd_pre(g,a, p,pi);
    1434     5000355 :     if (degpol(g)) split_add(S, Flx_normalize(g,p));
    1435             :   }
    1436     5035598 :   return 1;
    1437             : }
    1438             : 
    1439             : /* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
    1440             : GEN
    1441    24668670 : Flx_roots_pre(GEN f, ulong p, ulong pi)
    1442             : {
    1443             :   GEN pol;
    1444    24668670 :   long v = Flx_valrem(f, &f), n = degpol(f);
    1445             :   ulong q, PI;
    1446             :   struct split_t S;
    1447             : 
    1448    24658984 :   f = Flx_normalize(f, p);
    1449             :   /* optimization: test for small degree first */
    1450    24667203 :   if (n == 1)
    1451             :   {
    1452      114918 :     q = p - f[2];
    1453      114918 :     return v? mkvecsmall2(0, q): mkvecsmall(q);
    1454             :   }
    1455    24552285 :   PI = pi? pi: get_Fl_red(p); /* PI for Fp, pi for Fp[x] */
    1456    24553845 :   if (n == 2)
    1457             :   {
    1458    19485624 :     ulong r = Flx_quad_root(f, p, PI, 1), s;
    1459    19520586 :     if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
    1460    13312220 :     s = Flx_otherroot(f,r, p);
    1461    13356719 :     if (r < s)
    1462     3356277 :       return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
    1463    10000442 :     else if (r > s)
    1464    10035924 :       return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
    1465             :     else
    1466        4832 :       return v? mkvecsmall2(0, s): mkvecsmall(s);
    1467             :   }
    1468     5068221 :   if (SMALL_ULONG(p)) pi = 0; /* bilinear ops faster without Fl_*_pre */
    1469     5068221 :   q = p >> 1;
    1470     5068221 :   split_init(&S, lg(f)-1);
    1471     5066832 :   settyp(S.done, t_VECSMALL);
    1472     5066832 :   if (v) split_add_done(&S, (GEN)0);
    1473     5066832 :   if (! split_Flx_cut_out_roots(&S, f, p, pi))
    1474         600 :     return all_roots_mod_p(p, lg(S.done) == 1);
    1475     5065458 :   pol = polx_Flx(f[1]);
    1476     5066370 :   for (pol[2]=1; ; pol[2]++)
    1477     5467395 :   {
    1478    10533765 :     long j, l = lg(S.todo);
    1479    10533765 :     if (l == 1) { vecsmall_sort(S.done); return S.done; }
    1480     5467515 :     if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
    1481    11424914 :     for (j = 1; j < l; j++)
    1482             :     {
    1483     5957519 :       GEN b, c = gel(S.todo,j);
    1484             :       ulong r, s;
    1485     5957519 :       switch(degpol(c))
    1486             :       {
    1487     4973190 :         case 1:
    1488     4973190 :           split_moveto_done(&S, j, (GEN)(p - c[2]));
    1489     4973190 :           j--; l--; break;
    1490      460140 :         case 2:
    1491      460140 :           r = Flx_quad_root(c, p, PI, 0);
    1492      460153 :           if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
    1493      460146 :           s = Flx_otherroot(c,r, p);
    1494      460141 :           split_done(&S, j, (GEN)r, (GEN)s);
    1495      460148 :           j--; l--; break;
    1496      524070 :         default:
    1497      524070 :           b = Flxq_powu_pre(pol,q, c,p,pi); /* pol^(p-1)/2 */
    1498      524175 :           if (degpol(b) <= 0) continue;
    1499      410780 :           b = Flx_gcd_pre(c,Flx_Fl_add(b,p-1,p), p, pi);
    1500      410867 :           if (!degpol(b)) continue;
    1501      410351 :           b = Flx_normalize(b, p);
    1502      410385 :           c = Flx_div_pre(c,b, p,pi);
    1503      410378 :           split_todo(&S, j, b, c);
    1504             :       }
    1505             :     }
    1506             :   }
    1507             : }
    1508             : 
    1509             : GEN
    1510        1372 : Flx_roots(GEN f, ulong p)
    1511             : {
    1512        1372 :   pari_sp av = avma;
    1513             :   ulong pi;
    1514        1372 :   switch(lg(f))
    1515             :   {
    1516           0 :     case 2: pari_err_ROOTS0("Flx_roots");
    1517           0 :     case 3: set_avma(av); return cgetg(1, t_VECSMALL);
    1518             :   }
    1519        1372 :   if (p == 2) return Flx_root_mod_2(f);
    1520        1365 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    1521        1365 :   return gerepileuptoleaf(av, Flx_roots_pre(f, p, pi));
    1522             : }
    1523             : 
    1524             : /* assume x reduced mod p, monic. */
    1525             : static int
    1526     2490163 : Flx_quad_factortype(GEN x, ulong p)
    1527             : {
    1528     2490163 :   ulong b = x[3], c = x[2];
    1529     2490163 :   return krouu(Fl_disc_bc(b, c, p), p);
    1530             : }
    1531             : static GEN
    1532          56 : Flx_is_irred_2(GEN f, ulong p, long d)
    1533             : {
    1534          56 :   if (!d) return NULL;
    1535          56 :   if (d == 1) return gen_1;
    1536          56 :   return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
    1537             : }
    1538             : static GEN
    1539     2516204 : Flx_degfact_2(GEN f, ulong p, long d)
    1540             : {
    1541     2516204 :   if (!d) return trivial_fact();
    1542     2516204 :   if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
    1543     2490164 :   switch(Flx_quad_factortype(f, p)) {
    1544     1180492 :     case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
    1545     1279228 :     case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
    1546       30380 :     default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
    1547             :   }
    1548             : }
    1549             : /* p > 2 */
    1550             : static GEN
    1551     1932445 : Flx_factor_2(GEN f, ulong p, long d)
    1552             : {
    1553             :   ulong r, s;
    1554             :   GEN R,S;
    1555     1932445 :   long v = f[1];
    1556     1932445 :   if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
    1557     1855848 :   if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
    1558     1536419 :   r = Flx_quad_root(f, p, get_Fl_red(p), 1);
    1559     1536464 :   if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
    1560      923194 :   s = Flx_otherroot(f, r, p);
    1561      923191 :   r = Fl_neg(r, p);
    1562      923190 :   s = Fl_neg(s, p);
    1563      923191 :   if (s < r) lswap(s,r);
    1564      923191 :   R = mkvecsmall3(v,r,1);
    1565      923191 :   if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
    1566      812372 :   S = mkvecsmall3(v,s,1);
    1567      812375 :   return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
    1568             : }
    1569             : static GEN
    1570     4448620 : Flx_factor_deg2(GEN f, ulong p, long d, long flag)
    1571             : {
    1572     4448620 :   switch(flag) {
    1573          56 :     case 2: return Flx_is_irred_2(f, p, d);
    1574     2516229 :     case 1: return Flx_degfact_2(f, p, d);
    1575     1932335 :     default: return Flx_factor_2(f, p, d);
    1576             :   }
    1577             : }
    1578             : 
    1579             : static GEN
    1580       23324 : F2x_Berlekamp_ker(GEN u)
    1581             : {
    1582       23324 :   pari_sp ltop=avma;
    1583       23324 :   long j,N = F2x_degree(u);
    1584             :   GEN Q;
    1585             :   pari_timer T;
    1586       23323 :   timer_start(&T);
    1587       23325 :   Q = F2x_matFrobenius(u);
    1588      318577 :   for (j=1; j<=N; j++)
    1589      295251 :     F2m_flip(Q,j,j);
    1590       23326 :   if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
    1591       23326 :   Q = F2m_ker_sp(Q,0);
    1592       23326 :   if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
    1593       23326 :   return gerepileupto(ltop,Q);
    1594             : }
    1595             : #define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
    1596             : static long
    1597       35825 : F2x_split_Berlekamp(GEN *t)
    1598             : {
    1599       35825 :   GEN u = *t, a, b, vker;
    1600       35825 :   long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
    1601             : 
    1602       35822 :   if (du == 1) return 1;
    1603       27367 :   if (du == 2)
    1604             :   {
    1605        4045 :     if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
    1606             :     {
    1607           0 :       t[0] = mkvecsmall2(sv, 2);
    1608           0 :       t[1] = mkvecsmall2(sv, 3);
    1609           0 :       return 2;
    1610             :     }
    1611        4045 :     return 1;
    1612             :   }
    1613             : 
    1614       23322 :   vker = F2x_Berlekamp_ker(u);
    1615       23326 :   lb = lgcols(vker);
    1616       23343 :   d = lg(vker)-1;
    1617       23343 :   ir = 0;
    1618             :   /* t[i] irreducible for i < ir, still to be treated for i < L */
    1619       60437 :   for (L=1; L<d; )
    1620             :   {
    1621             :     GEN pol;
    1622       37113 :     if (d == 2)
    1623        5552 :       pol = F2v_to_F2x(gel(vker,2), sv);
    1624             :     else
    1625             :     {
    1626       31561 :       GEN v = zero_zv(lb);
    1627       31564 :       v[1] = du;
    1628       31564 :       v[2] = random_Fl(2); /*Assume vker[1]=1*/
    1629      118829 :       for (i=2; i<=d; i++)
    1630       87265 :         if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
    1631       31564 :       pol = F2v_to_F2x(v, sv);
    1632             :     }
    1633      111756 :     for (i=ir; i<L && L<d; i++)
    1634             :     {
    1635       74662 :       a = t[i]; la = F2x_degree(a);
    1636       74654 :       if (la == 1) { set_irred(i); }
    1637       74439 :       else if (la == 2)
    1638             :       {
    1639         725 :         if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
    1640             :         {
    1641           0 :           t[i] = mkvecsmall2(sv, 2);
    1642           0 :           t[L] = mkvecsmall2(sv, 3); L++;
    1643             :         }
    1644         725 :         set_irred(i);
    1645             :       }
    1646             :       else
    1647             :       {
    1648       73714 :         pari_sp av = avma;
    1649             :         long lb;
    1650       73714 :         b = F2x_rem(pol, a);
    1651       73680 :         if (F2x_degree(b) <= 0) { set_avma(av); continue; }
    1652       26613 :         b = F2x_gcd(a,b); lb = F2x_degree(b);
    1653       26611 :         if (lb && lb < la)
    1654             :         {
    1655       26611 :           t[L] = F2x_div(a,b);
    1656       26612 :           t[i]= b; L++;
    1657             :         }
    1658           0 :         else set_avma(av);
    1659             :       }
    1660             :     }
    1661             :   }
    1662       23324 :   return d;
    1663             : }
    1664             : /* assume deg f > 2 */
    1665             : static GEN
    1666       35390 : F2x_Berlekamp_i(GEN f, long flag)
    1667             : {
    1668       35390 :   long lfact, val, d = F2x_degree(f), j, k, lV;
    1669             :   GEN y, E, t, V;
    1670             : 
    1671       35389 :   val = F2x_valrem(f, &f);
    1672       35390 :   if (flag == 2 && val) return NULL;
    1673       35376 :   V = F2x_factor_squarefree(f); lV = lg(V);
    1674       35373 :   if (flag == 2 && lV > 2) return NULL;
    1675             : 
    1676             :   /* to hold factors and exponents */
    1677       35303 :   t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
    1678       35303 :   E = cgetg(d+1,t_VECSMALL);
    1679       35306 :   lfact = 1;
    1680       35306 :   if (val) {
    1681       11275 :     if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
    1682       11275 :     E[1] = val; lfact++;
    1683             :   }
    1684             : 
    1685      117224 :   for (k=1; k<lV; k++)
    1686             :   {
    1687       82072 :     if (F2x_degree(gel(V, k))==0) continue;
    1688       35825 :     gel(t,lfact) = gel(V, k);
    1689       35825 :     d = F2x_split_Berlekamp(&gel(t,lfact));
    1690       35824 :     if (flag == 2 && d != 1) return NULL;
    1691       35670 :     if (flag == 1)
    1692       47783 :       for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
    1693       97541 :     for (j=0; j<d; j++) E[lfact+j] = k;
    1694       35670 :     lfact += d;
    1695             :   }
    1696       35152 :   if (flag == 2) return gen_1; /* irreducible */
    1697       35138 :   setlg(t, lfact);
    1698       35136 :   setlg(E, lfact); y = mkvec2(t,E);
    1699       24343 :   return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
    1700       59481 :               : sort_factor_pol(y, cmpGuGu);
    1701             : }
    1702             : 
    1703             : /* Adapted from Shoup NTL */
    1704             : GEN
    1705      335087 : F2x_factor_squarefree(GEN f)
    1706             : {
    1707             :   GEN r, t, v, tv;
    1708      335087 :   long i, q, n = F2x_degree(f);
    1709      335087 :   GEN u = const_vec(n+1, pol1_F2x(f[1]));
    1710      335091 :   for(q = 1;;q *= 2)
    1711             :   {
    1712      543428 :     r = F2x_gcd(f, F2x_deriv(f));
    1713      543423 :     if (F2x_degree(r) == 0)
    1714             :     {
    1715      241000 :       gel(u, q) = f;
    1716      241000 :       break;
    1717             :     }
    1718      302427 :     t = F2x_div(f, r);
    1719      302424 :     if (F2x_degree(t) > 0)
    1720             :     {
    1721             :       long j;
    1722      146481 :       for(j = 1;;j++)
    1723             :       {
    1724      351458 :         v = F2x_gcd(r, t);
    1725      351455 :         tv = F2x_div(t, v);
    1726      351457 :         if (F2x_degree(tv) > 0)
    1727      148452 :           gel(u, j*q) = tv;
    1728      351457 :         if (F2x_degree(v) <= 0) break;
    1729      204975 :         r = F2x_div(r, v);
    1730      204977 :         t = v;
    1731             :       }
    1732      146480 :       if (F2x_degree(r) == 0) break;
    1733             :     }
    1734      208333 :     f = F2x_sqrt(r);
    1735             :   }
    1736     1194243 :   for (i = n; i; i--)
    1737     1191567 :     if (F2x_degree(gel(u,i))) break;
    1738      335066 :   setlg(u,i+1); return u;
    1739             : }
    1740             : 
    1741             : static GEN
    1742      350795 : F2x_ddf_simple(GEN T, GEN XP)
    1743             : {
    1744      350795 :   pari_sp av = avma, av2;
    1745             :   GEN f, z, Tr, X;
    1746      350795 :   long j, n = F2x_degree(T), v = T[1], B = n/2;
    1747      350796 :   if (n == 0) return cgetg(1, t_VEC);
    1748      350798 :   if (n == 1) return mkvec(T);
    1749      136951 :   z = XP; Tr = T; X = polx_F2x(v);
    1750      136950 :   f = const_vec(n, pol1_F2x(v));
    1751      136951 :   av2 = avma;
    1752      356190 :   for (j = 1; j <= B; j++)
    1753             :   {
    1754      230481 :     GEN u = F2x_gcd(Tr, F2x_add(z, X));
    1755      230515 :     if (F2x_degree(u))
    1756             :     {
    1757       50715 :       gel(f, j) = u;
    1758       50715 :       Tr = F2x_div(Tr, u);
    1759       50712 :       av2 = avma;
    1760      179823 :     } else z = gerepileuptoleaf(av2, z);
    1761      230538 :     if (!F2x_degree(Tr)) break;
    1762      219276 :     z = F2xq_sqr(z, Tr);
    1763             :   }
    1764      136928 :   if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
    1765      136946 :   return gerepilecopy(av, f);
    1766             : }
    1767             : 
    1768             : GEN
    1769           7 : F2x_ddf(GEN T)
    1770             : {
    1771             :   GEN XP;
    1772           7 :   T = F2x_get_red(T);
    1773           7 :   XP = F2x_Frobenius(T);
    1774           7 :   return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
    1775             : }
    1776             : 
    1777             : static GEN
    1778       20080 : F2xq_frobtrace(GEN a, long d, GEN T)
    1779             : {
    1780       20080 :   pari_sp av = avma;
    1781             :   long i;
    1782       20080 :   GEN x = a;
    1783       59188 :   for(i=1; i<d; i++)
    1784             :   {
    1785       39109 :     x = F2x_add(a, F2xq_sqr(x,T));
    1786       39108 :     if (gc_needed(av, 2))
    1787           0 :       x = gerepileuptoleaf(av, x);
    1788             :   }
    1789       20079 :   return x;
    1790             : }
    1791             : 
    1792             : static void
    1793       29853 : F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
    1794             : {
    1795       29853 :   long n = F2x_degree(Tp), r = n/d;
    1796             :   GEN T, f, ff;
    1797       29852 :   if (r==1) { gel(V, idx) = Tp; return; }
    1798       10045 :   T = Tp;
    1799       10045 :   XP = F2x_rem(XP, T);
    1800             :   while (1)
    1801       10033 :   {
    1802       20079 :     pari_sp btop = avma;
    1803             :     long df;
    1804       20079 :     GEN g = random_F2x(n, Tp[1]);
    1805       20080 :     GEN t = F2xq_frobtrace(g, d, T);
    1806       20080 :     if (lgpol(t) == 0) continue;
    1807       15132 :     f = F2x_gcd(t, Tp); df = F2x_degree(f);
    1808       15132 :     if (df > 0 && df < n) break;
    1809        5085 :     set_avma(btop);
    1810             :   }
    1811       10047 :   ff = F2x_div(Tp, f);
    1812       10047 :   F2x_edf_simple(f, XP, d, V, idx);
    1813       10047 :   F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
    1814             : }
    1815             : 
    1816             : static GEN
    1817      350793 : F2x_factor_Shoup(GEN T)
    1818             : {
    1819      350793 :   long i, n, s = 0;
    1820             :   GEN XP, D, V;
    1821             :   pari_timer ti;
    1822      350793 :   n = F2x_degree(T);
    1823      350791 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1824      350791 :   XP = F2x_Frobenius(T);
    1825      350789 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1826      350789 :   D = F2x_ddf_simple(T, XP);
    1827      350794 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1828     1094127 :   for (i = 1; i <= n; i++)
    1829      743337 :     s += F2x_degree(gel(D,i))/i;
    1830      350790 :   V = cgetg(s+1, t_COL);
    1831     1094144 :   for (i = 1, s = 1; i <= n; i++)
    1832             :   {
    1833      743333 :     GEN Di = gel(D,i);
    1834      743333 :     long ni = F2x_degree(Di), ri = ni/i;
    1835      743329 :     if (ni == 0) continue;
    1836      390256 :     if (ni == i) { gel(V, s++) = Di; continue; }
    1837        9736 :     F2x_edf_simple(Di, XP, i, V, s);
    1838        9760 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
    1839        9760 :     s += ri;
    1840             :   }
    1841      350811 :   return V;
    1842             : }
    1843             : 
    1844             : static GEN
    1845      299714 : F2x_factor_Cantor(GEN T)
    1846             : {
    1847      299714 :   GEN E, F, V = F2x_factor_squarefree(T);
    1848      299715 :   long i, j, l = lg(V);
    1849      299715 :   E = cgetg(l, t_VEC);
    1850      299715 :   F = cgetg(l, t_VEC);
    1851     1076131 :   for (i=1, j=1; i < l; i++)
    1852      776415 :     if (F2x_degree(gel(V,i)))
    1853             :     {
    1854      350793 :       GEN Fj = F2x_factor_Shoup(gel(V,i));
    1855      350791 :       gel(F, j) = Fj;
    1856      350791 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1857      350793 :       j++;
    1858             :     }
    1859      299716 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    1860             : }
    1861             : 
    1862             : #if 0
    1863             : static GEN
    1864             : F2x_simplefact_Shoup(GEN T)
    1865             : {
    1866             :   long i, n, s = 0, j = 1, k;
    1867             :   GEN XP, D, V;
    1868             :   pari_timer ti;
    1869             :   n = F2x_degree(T);
    1870             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1871             :   XP = F2x_Frobenius(T);
    1872             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1873             :   D = F2x_ddf_simple(T, XP);
    1874             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1875             :   for (i = 1; i <= n; i++)
    1876             :     s += F2x_degree(gel(D,i))/i;
    1877             :   V = cgetg(s+1, t_VECSMALL);
    1878             :   for (i = 1; i <= n; i++)
    1879             :   {
    1880             :     long ni = F2x_degree(gel(D,i)), ri = ni/i;
    1881             :     if (ni == 0) continue;
    1882             :     for (k = 1; k <= ri; k++)
    1883             :       V[j++] = i;
    1884             :   }
    1885             :   return V;
    1886             : }
    1887             : static GEN
    1888             : F2x_simplefact_Cantor(GEN T)
    1889             : {
    1890             :   GEN E, F, V = F2x_factor_squarefree(T);
    1891             :   long i, j, l = lg(V);
    1892             :   F = cgetg(l, t_VEC);
    1893             :   E = cgetg(l, t_VEC);
    1894             :   for (i=1, j=1; i < l; i++)
    1895             :     if (F2x_degree(gel(V,i)))
    1896             :     {
    1897             :       GEN Fj = F2x_simplefact_Shoup(gel(V,i));
    1898             :       gel(F, j) = Fj;
    1899             :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    1900             :       j++;
    1901             :     }
    1902             :   return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
    1903             : }
    1904             : static int
    1905             : F2x_isirred_Cantor(GEN T)
    1906             : {
    1907             :   pari_sp av = avma;
    1908             :   pari_timer ti;
    1909             :   long n;
    1910             :   GEN dT = F2x_deriv(T);
    1911             :   GEN XP, D;
    1912             :   if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
    1913             :   n = F2x_degree(T);
    1914             :   if (DEBUGLEVEL>=6) timer_start(&ti);
    1915             :   XP = F2x_Frobenius(T);
    1916             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
    1917             :   D = F2x_ddf_simple(T, XP);
    1918             :   if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
    1919             :   return gc_bool(av, F2x_degree(gel(D,n)) == n);
    1920             : }
    1921             : #endif
    1922             : 
    1923             : /* driver for Cantor factorization, assume deg f > 2; not competitive for
    1924             :  * flag != 0, or as deg f increases */
    1925             : static GEN
    1926      299714 : F2x_Cantor_i(GEN f, long flag)
    1927             : {
    1928             :   switch(flag)
    1929             :   {
    1930      299714 :     default: return F2x_factor_Cantor(f);
    1931             : #if 0
    1932             :     case 1: return F2x_simplefact_Cantor(f);
    1933             :     case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
    1934             : #endif
    1935             :   }
    1936             : }
    1937             : static GEN
    1938     1519692 : F2x_factor_i(GEN f, long flag)
    1939             : {
    1940     1519692 :   long d = F2x_degree(f);
    1941     1519684 :   if (d <= 2) return F2x_factor_deg2(f,d,flag);
    1942      310508 :   return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
    1943      645604 :                                : F2x_Berlekamp_i(f, flag);
    1944             : }
    1945             : 
    1946             : GEN
    1947           0 : F2x_degfact(GEN f)
    1948             : {
    1949           0 :   pari_sp av = avma;
    1950           0 :   GEN z = F2x_factor_i(f, 1);
    1951           0 :   return gerepilecopy(av, z);
    1952             : }
    1953             : 
    1954             : int
    1955         238 : F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
    1956             : 
    1957             : /* Adapted from Shoup NTL */
    1958             : GEN
    1959     7242311 : Flx_factor_squarefree_pre(GEN f, ulong p, ulong pi)
    1960             : {
    1961     7242311 :   long i, q, n = degpol(f);
    1962     7242089 :   GEN u = const_vec(n+1, pol1_Flx(f[1]));
    1963     7242639 :   for(q = 1;;q *= p)
    1964      160862 :   {
    1965     7403501 :     GEN t, v, tv, r = Flx_gcd_pre(f, Flx_deriv(f, p), p, pi);
    1966     7400857 :     if (degpol(r) == 0) { gel(u, q) = f; break; }
    1967      653378 :     t = Flx_div_pre(f, r, p, pi);
    1968      653379 :     if (degpol(t) > 0)
    1969             :     {
    1970             :       long j;
    1971      499240 :       for(j = 1;;j++)
    1972             :       {
    1973     1234512 :         v = Flx_gcd_pre(r, t, p, pi);
    1974     1234498 :         tv = Flx_div_pre(t, v, p, pi);
    1975     1234505 :         if (degpol(tv) > 0)
    1976      744649 :           gel(u, j*q) = Flx_normalize(tv, p);
    1977     1234505 :         if (degpol(v) <= 0) break;
    1978      735263 :         r = Flx_div_pre(r, v, p, pi);
    1979      735272 :         t = v;
    1980             :       }
    1981      499242 :       if (degpol(r) == 0) break;
    1982             :     }
    1983      160862 :     f = Flx_normalize(Flx_deflate(r, p), p);
    1984             :   }
    1985    30449340 :   for (i = n; i; i--)
    1986    30447046 :     if (degpol(gel(u,i))) break;
    1987     7241631 :   setlg(u,i+1); return u;
    1988             : }
    1989             : GEN
    1990      341118 : Flx_factor_squarefree(GEN f, ulong p)
    1991      341118 : { return Flx_factor_squarefree_pre(f, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    1992             : 
    1993             : long
    1994        3451 : Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
    1995             : {
    1996        3451 :   pari_sp av = avma;
    1997             :   ulong lc, pi;
    1998             :   GEN F;
    1999        3451 :   long i, n = degpol(f), v = f[1], l;
    2000        3451 :   if (n % k) return 0;
    2001        3451 :   lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
    2002        3451 :   if (lc == ULONG_MAX) { av = avma; return 0; }
    2003        3451 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2004        3451 :   F = Flx_factor_squarefree_pre(f, p, pi); l = lg(F)-1;
    2005       38325 :   for (i = 1; i <= l; i++)
    2006       34895 :     if (i%k && degpol(gel(F,i))) return gc_long(av,0);
    2007        3430 :   if (pt_r)
    2008             :   {
    2009        3430 :     GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
    2010       38304 :     for(i = l; i >= 1; i--)
    2011             :     {
    2012       34874 :       if (i%k) continue;
    2013       12453 :       s = Flx_mul_pre(s, gel(F,i), p, pi);
    2014       12453 :       r = Flx_mul_pre(r, s, p, pi);
    2015             :     }
    2016        3430 :     *pt_r = gerepileuptoleaf(av, r);
    2017           0 :   } else set_avma(av);
    2018        3430 :   return 1;
    2019             : }
    2020             : 
    2021             : /* See <http://www.shoup.net/papers/factorimpl.pdf> */
    2022             : static GEN
    2023     9162785 : Flx_ddf_Shoup(GEN T, GEN XP, ulong p, ulong pi)
    2024             : {
    2025     9162785 :   pari_sp av = avma;
    2026             :   GEN b, g, h, F, f, Tr, xq;
    2027             :   long i, j, n, v, bo, ro;
    2028             :   long B, l, m;
    2029             :   pari_timer ti;
    2030     9162785 :   n = get_Flx_degree(T); v = get_Flx_var(T);
    2031     9163895 :   if (n == 0) return cgetg(1, t_VEC);
    2032     9120232 :   if (n == 1) return mkvec(get_Flx_mod(T));
    2033     8736041 :   B = n/2;
    2034     8736041 :   l = usqrt(B);
    2035     8736365 :   m = (B+l-1)/l;
    2036     8736365 :   T = Flx_get_red(T, p);
    2037     8735959 :   b = cgetg(l+2, t_VEC);
    2038     8735607 :   gel(b, 1) = polx_Flx(v);
    2039     8736140 :   gel(b, 2) = XP;
    2040     8736140 :   bo = brent_kung_optpow(n, l-1, 1);
    2041     8736701 :   ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
    2042     8736701 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2043     8736701 :   if (expu(p) <= ro)
    2044      875051 :     for (i = 3; i <= l+1; i++)
    2045      485094 :       gel(b, i) = Flxq_powu_pre(gel(b, i-1), p, T, p, pi);
    2046             :   else
    2047             :   {
    2048     8346391 :     xq = Flxq_powers_pre(gel(b, 2), bo,  T, p, pi);
    2049     8345179 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
    2050     9066088 :     for (i = 3; i <= l+1; i++)
    2051      720746 :       gel(b, i) = Flx_FlxqV_eval_pre(gel(b, i-1), xq, T, p, pi);
    2052             :   }
    2053     8735299 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
    2054     8735299 :   xq = Flxq_powers_pre(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p, pi);
    2055     8735514 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
    2056     8735514 :   g = cgetg(m+1, t_VEC);
    2057     8736282 :   gel(g, 1) = gel(xq, 2);
    2058    14139211 :   for(i = 2; i <= m; i++)
    2059     5402612 :     gel(g, i) = Flx_FlxqV_eval_pre(gel(g, i-1), xq, T, p, pi);
    2060     8736599 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
    2061     8736599 :   h = cgetg(m+1, t_VEC);
    2062    22874370 :   for (j = 1; j <= m; j++)
    2063             :   {
    2064    14138564 :     pari_sp av = avma;
    2065    14138564 :     GEN gj = gel(g, j);
    2066    14138564 :     GEN e = Flx_sub(gj, gel(b, 1), p);
    2067    18009656 :     for (i = 2; i <= l; i++)
    2068     3872975 :       e = Flxq_mul_pre(e, Flx_sub(gj, gel(b, i), p), T, p, pi);
    2069    14136681 :     gel(h, j) = gerepileupto(av, e);
    2070             :   }
    2071     8735806 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
    2072     8735806 :   Tr = get_Flx_mod(T);
    2073     8736373 :   F = cgetg(m+1, t_VEC);
    2074    22875348 :   for (j = 1; j <= m; j++)
    2075             :   {
    2076    14138357 :     GEN u = Flx_gcd_pre(Tr, gel(h, j), p, pi);
    2077    14137568 :     if (degpol(u))
    2078             :     {
    2079     6538508 :       u = Flx_normalize(u, p);
    2080     6539045 :       Tr = Flx_div_pre(Tr, u, p, pi);
    2081             :     }
    2082    14137803 :     gel(F, j) = u;
    2083             :   }
    2084     8736991 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
    2085     8736991 :   f = const_vec(n, pol1_Flx(v));
    2086    22876441 :   for (j = 1; j <= m; j++)
    2087             :   {
    2088    14139679 :     GEN e = gel(F, j);
    2089    15017762 :     for (i=l-1; i >= 0; i--)
    2090             :     {
    2091    15017845 :       GEN u = Flx_gcd_pre(e, Flx_sub(gel(g, j), gel(b, i+1), p), p, pi);
    2092    15015967 :       if (degpol(u))
    2093             :       {
    2094     6638704 :         gel(f, l*j-i) = u;
    2095     6638704 :         e = Flx_div_pre(e, u, p, pi);
    2096             :       }
    2097    15015638 :       if (!degpol(e)) break;
    2098             :     }
    2099             :   }
    2100     8736762 :   if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
    2101     8736762 :   if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
    2102     8736615 :   return gerepilecopy(av, f);
    2103             : }
    2104             : 
    2105             : static void
    2106      408873 : Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
    2107             : {
    2108      408873 :   long n = degpol(Tp), r = n/d;
    2109             :   GEN T, f, ff;
    2110             :   ulong p2;
    2111      408872 :   if (r==1) { gel(V, idx) = Tp; return; }
    2112      174698 :   p2 = p>>1;
    2113      174698 :   T = Flx_get_red_pre(Tp, p, pi);
    2114      174698 :   XP = Flx_rem_pre(XP, T, p, pi);
    2115             :   while (1)
    2116       20665 :   {
    2117      195363 :     pari_sp btop = avma;
    2118             :     long i;
    2119      195363 :     GEN g = random_Flx(n, Tp[1], p);
    2120      195363 :     GEN t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
    2121      195364 :     if (lgpol(t) == 0) continue;
    2122      434324 :     for(i=1; i<=10; i++)
    2123             :     {
    2124      418729 :       pari_sp btop2 = avma;
    2125      418729 :       GEN R = Flxq_powu_pre(Flx_Fl_add(t, random_Fl(p), p), p2, T, p, pi);
    2126      418722 :       f = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), Tp, p, pi);
    2127      418723 :       if (degpol(f) > 0 && degpol(f) < n) break;
    2128      244024 :       set_avma(btop2);
    2129             :     }
    2130      190293 :     if (degpol(f) > 0 && degpol(f) < n) break;
    2131       15596 :     set_avma(btop);
    2132             :   }
    2133      174699 :   f = Flx_normalize(f, p);
    2134      174699 :   ff = Flx_div_pre(Tp, f, p, pi);
    2135      174699 :   Flx_edf_simple(f, XP, d, p, pi, V, idx);
    2136      174697 :   Flx_edf_simple(ff, XP, d, p, pi, V, idx+degpol(f)/d);
    2137             : }
    2138             : static void
    2139             : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx);
    2140             : 
    2141             : static void
    2142     1357041 : Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, ulong pi,
    2143             :   GEN V, long idx)
    2144             : {
    2145             :   pari_sp av;
    2146     1357041 :   GEN Tp = get_Flx_mod(T);
    2147     1357042 :   long n = degpol(hp), vT = Tp[1];
    2148             :   GEN u1, u2, f1, f2;
    2149     1357048 :   ulong p2 = p>>1;
    2150             :   GEN R, h;
    2151     1357048 :   h = Flx_get_red_pre(hp, p, pi);
    2152     1357038 :   t = Flx_rem_pre(t, T, p, pi);
    2153     1356968 :   av = avma;
    2154             :   do
    2155             :   {
    2156     2284471 :     set_avma(av);
    2157     2284481 :     R = Flxq_powu_pre(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p, pi);
    2158     2284229 :     u1 = Flx_gcd_pre(Flx_Fl_add(R, p-1, p), hp, p, pi);
    2159     2284453 :   } while (degpol(u1)==0 || degpol(u1)==n);
    2160     1356953 :   f1 = Flx_gcd_pre(Flx_Flxq_eval_pre(u1, t, T, p, pi), Tp, p, pi);
    2161     1356975 :   f1 = Flx_normalize(f1, p);
    2162     1357035 :   u2 = Flx_div_pre(hp, u1, p, pi);
    2163     1357051 :   f2 = Flx_div_pre(Tp, f1, p, pi);
    2164     1357038 :   if (degpol(u1)==1)
    2165             :   {
    2166     1040111 :     if (degpol(f1)==d)
    2167     1025695 :       gel(V, idx) = f1;
    2168             :     else
    2169       14413 :       Flx_edf(f1, XP, d, p, pi, V, idx);
    2170             :   }
    2171             :   else
    2172      316970 :     Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, pi, V, idx);
    2173     1357082 :   idx += degpol(f1)/d;
    2174     1357063 :   if (degpol(u2)==1)
    2175             :   {
    2176     1034850 :     if (degpol(f2)==d)
    2177     1020869 :       gel(V, idx) = f2;
    2178             :     else
    2179       13981 :       Flx_edf(f2, XP, d, p, pi, V, idx);
    2180             :   }
    2181             :   else
    2182      322239 :     Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, pi, V, idx);
    2183     1357103 : }
    2184             : 
    2185             : static void
    2186      717899 : Flx_edf(GEN Tp, GEN XP, long d, ulong p, ulong pi, GEN V, long idx)
    2187             : {
    2188      717899 :   long n = degpol(Tp), r = n/d, vT = Tp[1];
    2189             :   GEN T, h, t;
    2190             :   pari_timer ti;
    2191      717897 :   if (r==1) { gel(V, idx) = Tp; return; }
    2192      717897 :   T = Flx_get_red_pre(Tp, p, pi);
    2193      717895 :   XP = Flx_rem_pre(XP, T, p, pi);
    2194      717889 :   if (DEBUGLEVEL>=7) timer_start(&ti);
    2195             :   do
    2196             :   {
    2197      739462 :     GEN g = random_Flx(n, vT, p);
    2198      739467 :     t = gel(Flxq_auttrace_pre(mkvec2(XP, g), d, T, p, pi), 2);
    2199      739473 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
    2200      739473 :     h = Flxq_minpoly_pre(t, T, p, pi);
    2201      739461 :     if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
    2202      739461 :   } while (degpol(h) <= 1);
    2203      717888 :   Flx_edf_rec(T, XP, h, t, d, p, pi, V, idx);
    2204             : }
    2205             : 
    2206             : static GEN
    2207     1653917 : Flx_factor_Shoup(GEN T, ulong p, ulong pi)
    2208             : {
    2209     1653917 :   long i, n, s = 0, e = expu(p);
    2210             :   GEN XP, D, V;
    2211             :   pari_timer ti;
    2212     1653917 :   n = get_Flx_degree(T);
    2213     1653916 :   T = Flx_get_red_pre(T, p, pi);
    2214     1653913 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2215     1653913 :   XP = Flx_Frobenius_pre(T, p, pi);
    2216     1653887 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2217     1653887 :   D = Flx_ddf_Shoup(T, XP, p, pi);
    2218     1653939 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2219     1653939 :   s = ddf_to_nbfact(D);
    2220     1653933 :   V = cgetg(s+1, t_COL);
    2221     9050962 :   for (i = 1, s = 1; i <= n; i++)
    2222             :   {
    2223     7397032 :     GEN Di = gel(D,i);
    2224     7397032 :     long ni = degpol(Di), ri = ni/i;
    2225     7397032 :     if (ni == 0) continue;
    2226     2149739 :     Di = Flx_normalize(Di, p);
    2227     2149771 :     if (ni == i) { gel(V, s++) = Di; continue; }
    2228      748974 :     if (ri <= e*expu(e))
    2229      689503 :       Flx_edf(Di, XP, i, p, pi, V, s);
    2230             :     else
    2231       59476 :       Flx_edf_simple(Di, XP, i, p, pi, V, s);
    2232      748977 :     if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
    2233      748976 :     s += ri;
    2234             :   }
    2235     1653930 :   return V;
    2236             : }
    2237             : 
    2238             : static GEN
    2239     1581966 : Flx_factor_Cantor(GEN T, ulong p)
    2240             : {
    2241     1581966 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2242     1581966 :   GEN E, F, V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi);
    2243     1581958 :   long i, j, l = lg(V);
    2244     1581958 :   F = cgetg(l, t_VEC);
    2245     1581960 :   E = cgetg(l, t_VEC);
    2246     3784095 :   for (i=1, j=1; i < l; i++)
    2247     2202119 :     if (degpol(gel(V,i)))
    2248             :     {
    2249     1653917 :       GEN Fj = Flx_factor_Shoup(gel(V,i), p, pi);
    2250     1653930 :       gel(F, j) = Fj;
    2251     1653930 :       gel(E, j) = const_vecsmall(lg(Fj)-1, i);
    2252     1653931 :       j++;
    2253             :     }
    2254     1581976 :   return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
    2255             : }
    2256             : 
    2257             : GEN
    2258           0 : Flx_ddf_pre(GEN T, ulong p, ulong pi)
    2259             : {
    2260             :   GEN XP;
    2261           0 :   T = Flx_get_red_pre(T, p, pi);
    2262           0 :   XP = Flx_Frobenius_pre(T, p, pi);
    2263           0 :   return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p, pi));
    2264             : }
    2265             : GEN
    2266           0 : Flx_ddf(GEN T, ulong p)
    2267           0 : { return Flx_ddf_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2268             : 
    2269             : static GEN
    2270     5316623 : Flx_simplefact_Cantor(GEN T, ulong p)
    2271             : {
    2272     5316623 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2273             :   long i, l;
    2274             :   GEN V;
    2275     5316623 :   T = Flx_get_red_pre(T, p, pi);
    2276     5316116 :   V = Flx_factor_squarefree_pre(get_Flx_mod(T), p, pi); l = lg(V);
    2277    10716283 :   for (i=1; i < l; i++)
    2278     5397940 :     gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius_pre(gel(V,i), p,pi), p,pi);
    2279     5318343 :   return vddf_to_simplefact(V, get_Flx_degree(T));
    2280             : }
    2281             : 
    2282             : static int
    2283        1078 : Flx_isirred_Cantor(GEN Tp, ulong p)
    2284             : {
    2285        1078 :   pari_sp av = avma;
    2286             :   pari_timer ti;
    2287        1078 :   GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
    2288        1078 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2289             :   long n;
    2290        1078 :   if (degpol(Flx_gcd_pre(T, dT, p, pi)) != 0) return gc_bool(av,0);
    2291         791 :   n = get_Flx_degree(T);
    2292         791 :   T = Flx_get_red_pre(Tp, p, pi);
    2293         791 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2294         791 :   XP = Flx_Frobenius_pre(T, p, pi);
    2295         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2296         791 :   D = Flx_ddf_Shoup(T, XP, p, pi);
    2297         791 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2298         791 :   return gc_bool(av, degpol(gel(D,n)) == n);
    2299             : }
    2300             : 
    2301             : /* f monic */
    2302             : static GEN
    2303    11422339 : Flx_factor_i(GEN f, ulong pp, long flag)
    2304             : {
    2305             :   long d;
    2306    11422339 :   if (pp==2) { /*We need to handle 2 specially */
    2307       74053 :     GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
    2308       74052 :     if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
    2309       74052 :     return F;
    2310             :   }
    2311    11348286 :   d = degpol(f);
    2312    11348178 :   if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
    2313     6899720 :   switch(flag)
    2314             :   {
    2315     1581966 :     default: return Flx_factor_Cantor(f, pp);
    2316     5316676 :     case 1: return Flx_simplefact_Cantor(f, pp);
    2317        1078 :     case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
    2318             :   }
    2319             : }
    2320             : 
    2321             : GEN
    2322     7829318 : Flx_degfact(GEN f, ulong p)
    2323             : {
    2324     7829318 :   pari_sp av = avma;
    2325     7829318 :   GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
    2326     7828682 :   return gerepilecopy(av, z);
    2327             : }
    2328             : 
    2329             : /* T must be squarefree mod p*/
    2330             : GEN
    2331     1452542 : Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
    2332             : {
    2333             :   GEN XP, D;
    2334             :   pari_timer ti;
    2335     1452542 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    2336     1452542 :   long i, s, n = get_Flx_degree(T);
    2337     1452520 :   GEN V = const_vecsmall(n, 0);
    2338     1452509 :   pari_sp av = avma;
    2339     1452509 :   T = Flx_get_red_pre(T, p, pi);
    2340     1452515 :   if (DEBUGLEVEL>=6) timer_start(&ti);
    2341     1452515 :   XP = Flx_Frobenius_pre(T, p, pi);
    2342     1452470 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
    2343     1452470 :   D = Flx_ddf_Shoup(T, XP, p, pi);
    2344     1452648 :   if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
    2345     7477182 :   for (i = 1, s = 0; i <= n; i++) { V[i] = degpol(gel(D,i))/i; s += V[i]; }
    2346     1452609 :   *nb = s; set_avma(av); return V;
    2347             : }
    2348             : 
    2349             : long
    2350      659472 : Flx_nbfact_Frobenius_pre(GEN T, GEN XP, ulong p, ulong pi)
    2351             : {
    2352      659472 :   pari_sp av = avma;
    2353      659472 :   long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p, pi));
    2354      659474 :   return gc_long(av,s);
    2355             : }
    2356             : long
    2357           0 : Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
    2358           0 : { return Flx_nbfact_Frobenius_pre(T, XP, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2359             : 
    2360             : /* T must be squarefree mod p*/
    2361             : long
    2362      659468 : Flx_nbfact_pre(GEN T, ulong p, ulong pi)
    2363             : {
    2364      659468 :   pari_sp av = avma;
    2365      659468 :   GEN XP = Flx_Frobenius_pre(T, p, pi);
    2366      659472 :   long n = Flx_nbfact_Frobenius_pre(T, XP, p, pi);
    2367      659474 :   return gc_long(av,n);
    2368             : }
    2369             : long
    2370      659468 : Flx_nbfact(GEN T, ulong p)
    2371      659468 : { return Flx_nbfact_pre(T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
    2372             : 
    2373             : int
    2374        1057 : Flx_is_irred(GEN f, ulong p)
    2375             : {
    2376        1057 :   pari_sp av = avma;
    2377        1057 :   f = Flx_normalize(f,p);
    2378        1057 :   return gc_bool(av, !!Flx_factor_i(f,p,2));
    2379             : }
    2380             : 
    2381             : /* Use this function when you think f is reducible, and that there are lots of
    2382             :  * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
    2383             : int
    2384         112 : FpX_is_irred(GEN f, GEN p)
    2385             : {
    2386         112 :   pari_sp av = avma;
    2387             :   int z;
    2388         112 :   switch(ZX_factmod_init(&f,p))
    2389             :   {
    2390          28 :     case 0:  z = !!F2x_factor_i(f,2); break;
    2391          77 :     case 1:  z = !!Flx_factor_i(f,p[2],2); break;
    2392           7 :     default: z = !!FpX_factor_i(f,p,2); break;
    2393             :   }
    2394         112 :   return gc_bool(av,z);
    2395             : }
    2396             : GEN
    2397       47947 : FpX_degfact(GEN f, GEN p) {
    2398       47947 :   pari_sp av = avma;
    2399             :   GEN F;
    2400       47947 :   switch(ZX_factmod_init(&f,p))
    2401             :   {
    2402         476 :     case 0:  F = F2x_factor_i(f,1); break;
    2403       47443 :     case 1:  F = Flx_factor_i(f,p[2],1); break;
    2404          28 :     default: F = FpX_factor_i(f,p,1); break;
    2405             :   }
    2406       47947 :   return gerepilecopy(av, F);
    2407             : }
    2408             : 
    2409             : #if 0
    2410             : /* set x <-- x + c*y mod p */
    2411             : /* x is not required to be normalized.*/
    2412             : static void
    2413             : Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
    2414             : {
    2415             :   long i, lx, ly;
    2416             :   ulong *x=(ulong *)gx;
    2417             :   ulong *y=(ulong *)gy;
    2418             :   if (!c) return;
    2419             :   lx = lg(gx);
    2420             :   ly = lg(gy);
    2421             :   if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
    2422             :   if (SMALL_ULONG(p))
    2423             :     for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
    2424             :   else
    2425             :     for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
    2426             : }
    2427             : #endif
    2428             : 
    2429             : GEN
    2430     4295341 : FpX_factor(GEN f, GEN p)
    2431             : {
    2432     4295341 :   pari_sp av = avma;
    2433             :   GEN F;
    2434     4295341 :   switch(ZX_factmod_init(&f, p))
    2435             :   {
    2436     1429824 :     case 0:  F = F2x_factor_i(f,0);
    2437     1429830 :              F2xV_to_ZXV_inplace(gel(F,1)); break;
    2438     2863130 :     case 1:  F = Flx_factor_i(f,p[2],0);
    2439     2863149 :              FlxV_to_ZXV_inplace(gel(F,1)); break;
    2440        2432 :     default: F = FpX_factor_i(f,p,0); break;
    2441             :   }
    2442     4295462 :   return gerepilecopy(av, F);
    2443             : }
    2444             : 
    2445             : GEN
    2446      681802 : Flx_factor(GEN f, ulong p)
    2447             : {
    2448      681802 :   pari_sp av = avma;
    2449      681802 :   return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
    2450             : }
    2451             : GEN
    2452       15077 : F2x_factor(GEN f)
    2453             : {
    2454       15077 :   pari_sp av = avma;
    2455       15077 :   return gerepilecopy(av, F2x_factor_i(f,0));
    2456             : }

Generated by: LCOV version 1.16