CoCalc -- FpX_factor.c

Testing latest pari + WASM + node.js... and it works?! Wow.
pari / 2021-07-01 / pari-2.14.0.alpha / src / basemath / FpX_factor.c
²⁸⁴⁷⁸ views
License: GPL3
ubuntu2004
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
ZX_factmod_init(GEN *F, GEN p)
35
{
36
  if (lgefint(p) == 3)
37
  {
38
    ulong pp = p[2];
39
    if (pp == 2) { *F = ZX_to_F2x(*F); return 0; }
40
    *F = ZX_to_Flx(*F, pp);
41
    if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
42
    return 1;
43
  }
44
  *F = FpX_red(*F, p);
45
  if (lg(*F) > 3) *F = FpX_normalize(*F, p);
46
  return 2;
47
}
48
static void
49
ZX_rootmod_init(GEN *F, GEN p)
50
{
51
  if (lgefint(p) == 3)
52
  {
53
    ulong pp = p[2];
54
    *F = ZX_to_Flx(*F, pp);
55
    if (lg(*F) > 3) *F = Flx_normalize(*F, pp);
56
  }
57
  else
58
  {
59
    *F = FpX_red(*F, p);
60
    if (lg(*F) > 3) *F = FpX_normalize(*F, p);
61
  }
62
}
63

64
/* return 1,...,p-1 [not_0 = 1] or 0,...,p [not_0 = 0] */
65
static GEN
66
all_roots_mod_p(ulong p, int not_0)
67
{
68
  GEN r;
69
  ulong i;
70
  if (not_0) {
71
    r = cgetg(p, t_VECSMALL);
72
    for (i = 1; i < p; i++) r[i] = i;
73
  } else {
74
    r = cgetg(p+1, t_VECSMALL);
75
    for (i = 0; i < p; i++) r[i+1] = i;
76
  }
77
  return r;
78
}
79

80
/* X^n - 1 */
81
static GEN
82
Flx_Xnm1(long sv, long n, ulong p)
83
{
84
  GEN t = cgetg(n+3, t_VECSMALL);
85
  long i;
86
  t[1] = sv;
87
  t[2] = p - 1;
88
  for (i = 3; i <= n+1; i++) t[i] = 0;
89
  t[i] = 1; return t;
90
}
91
/* X^n + 1 */
92
static GEN
93
Flx_Xn1(long sv, long n, ulong p)
94
{
95
  GEN t = cgetg(n+3, t_VECSMALL);
96
  long i;
97
  (void) p;
98
  t[1] = sv;
99
  t[2] = 1;
100
  for (i = 3; i <= n+1; i++) t[i] = 0;
101
  t[i] = 1; return t;
102
}
103

104
static GEN
105
Flx_root_mod_2(GEN f)
106
{
107
  int z1, z0 = !(f[2] & 1);
108
  long i,n;
109
  GEN y;
110

111
  for (i=2, n=1; i < lg(f); i++) n += f[i];
112
  z1 = n & 1;
113
  y = cgetg(z0+z1+1, t_VECSMALL); i = 1;
114
  if (z0) y[i++] = 0;
115
  if (z1) y[i  ] = 1;
116
  return y;
117
}
118
static ulong
119
Flx_oneroot_mod_2(GEN f)
120
{
121
  long i,n;
122
  if (!(f[2] & 1)) return 0;
123
  for (i=2, n=1; i < lg(f); i++) n += f[i];
124
  if (n & 1) return 1;
125
  return 2;
126
}
127

128
static GEN FpX_roots_i(GEN f, GEN p);
129
static GEN Flx_roots_i(GEN f, ulong p);
130

131
static int
132
cmpGuGu(GEN a, GEN b) { return (ulong)a < (ulong)b? -1: (a == b? 0: 1); }
133

134
/* Generic driver to computes the roots of f modulo pp, using 'Roots' when
135
 * pp is a small prime.
136
 * if (gpwrap), check types thoroughly and return t_INTMODs, otherwise
137
 * assume that f is an FpX, pp a prime and return t_INTs */
138
static GEN
139
rootmod_aux(GEN f, GEN pp)
140
{
141
  GEN y;
142
  switch(lg(f))
143
  {
144
    case 2: pari_err_ROOTS0("rootmod");
145
    case 3: return cgetg(1,t_COL);
146
  }
147
  if (typ(f) == t_VECSMALL)
148
  {
149
    ulong p = pp[2];
150
    if (p == 2)
151
      y = Flx_root_mod_2(f);
152
    else
153
    {
154
      if (!odd(p)) pari_err_PRIME("rootmod",utoi(p));
155
      y = Flx_roots_i(f, p);
156
    }
157
    y = Flc_to_ZC(y);
158
  }
159
  else
160
    y = FpX_roots_i(f, pp);
161
  return y;
162
}
163
/* assume that f is a ZX and p a prime */
164
GEN
165
FpX_roots(GEN f, GEN p)
166
{
167
  pari_sp av = avma;
168
  GEN y; ZX_rootmod_init(&f, p); y = rootmod_aux(f, p);
169
  return gerepileupto(av, y);
170
}
171

172
/* assume x reduced mod p > 2, monic. */
173
static int
174
FpX_quad_factortype(GEN x, GEN p)
175
{
176
  GEN b = gel(x,3), c = gel(x,2);
177
  GEN D = subii(sqri(b), shifti(c,2));
178
  return kronecker(D,p);
179
}
180
/* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
181
static GEN
182
FpX_quad_root(GEN x, GEN p, int unknown)
183
{
184
  GEN s, D, b = gel(x,3), c = gel(x,2);
185

186
  if (absequaliu(p, 2)) {
187
    if (!signe(b)) return c;
188
    return signe(c)? NULL: gen_1;
189
  }
190
  D = subii(sqri(b), shifti(c,2));
191
  D = remii(D,p);
192
  if (unknown && kronecker(D,p) == -1) return NULL;
193

194
  s = Fp_sqrt(D,p);
195
  /* p is not prime, go on and give e.g. maxord a chance to recover */
196
  if (!s) return NULL;
197
  return Fp_halve(Fp_sub(s,b, p), p);
198
}
199
static GEN
200
FpX_otherroot(GEN x, GEN r, GEN p)
201
{ return Fp_neg(Fp_add(gel(x,3), r, p), p); }
202

203
/* disc(x^2+bx+c) = b^2 - 4c */
204
static ulong
205
Fl_disc_bc(ulong b, ulong c, ulong p)
206
{ return Fl_sub(Fl_sqr(b,p), Fl_double(Fl_double(c,p),p), p); }
207
/* p > 2 */
208
static ulong
209
Flx_quad_root(GEN x, ulong p, int unknown)
210
{
211
  ulong s, b = x[3], c = x[2];
212
  ulong D = Fl_disc_bc(b, c, p);
213
  if (unknown && krouu(D,p) == -1) return p;
214
  s = Fl_sqrt(D,p);
215
  if (s==~0UL) return p;
216
  return Fl_halve(Fl_sub(s,b, p), p);
217
}
218
static ulong
219
Flx_otherroot(GEN x, ulong r, ulong p)
220
{ return Fl_neg(Fl_add(x[3], r, p), p); }
221

222
/* 'todo' contains the list of factors to be split.
223
 * 'done' the list of finished factors, no longer touched */
224
struct split_t { GEN todo, done; };
225
static void
226
split_init(struct split_t *S, long max)
227
{
228
  S->todo = vectrunc_init(max);
229
  S->done = vectrunc_init(max);
230
}
231
#if 0
232
/* move todo[i] to done */
233
static void
234
split_convert(struct split_t *S, long i)
235
{
236
  long n = lg(S->todo)-1;
237
  vectrunc_append(S->done, gel(S->todo,i));
238
  if (n) gel(S->todo,i) = gel(S->todo, n);
239
  setlg(S->todo, n);
240
}
241
#endif
242
/* append t to todo */
243
static void
244
split_add(struct split_t *S, GEN t) { vectrunc_append(S->todo, t); }
245
/* delete todo[i], add t to done */
246
static void
247
split_moveto_done(struct split_t *S, long i, GEN t)
248
{
249
  long n = lg(S->todo)-1;
250
  vectrunc_append(S->done, t);
251
  if (n) gel(S->todo,i) = gel(S->todo, n);
252
  setlg(S->todo, n);
253

254
}
255
/* append t to done */
256
static void
257
split_add_done(struct split_t *S, GEN t)
258
{ vectrunc_append(S->done, t); }
259
/* split todo[i] into a and b */
260
static void
261
split_todo(struct split_t *S, long i, GEN a, GEN b)
262
{
263
  gel(S->todo, i) = a;
264
  split_add(S, b);
265
}
266
/* split todo[i] into a and b, moved to done */
267
static void
268
split_done(struct split_t *S, long i, GEN a, GEN b)
269
{
270
  split_moveto_done(S, i, a);
271
  split_add_done(S, b);
272
}
273

274
/* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
275
static GEN
276
FpX_roots_i(GEN f, GEN p)
277
{
278
  GEN pol, pol0, a, q;
279
  struct split_t S;
280

281
  split_init(&S, lg(f)-1);
282
  settyp(S.done, t_COL);
283
  if (ZX_valrem(f, &f)) split_add_done(&S, gen_0);
284
  switch(degpol(f))
285
  {
286
    case 0: return ZC_copy(S.done);
287
    case 1: split_add_done(&S, subii(p, gel(f,2))); return ZC_copy(S.done);
288
    case 2: {
289
      GEN s, r = FpX_quad_root(f, p, 1);
290
      if (r) {
291
        split_add_done(&S, r);
292
        s = FpX_otherroot(f,r, p);
293
        /* f not known to be square free yet */
294
        if (!equalii(r, s)) split_add_done(&S, s);
295
      }
296
      return sort(S.done);
297
    }
298
  }
299

300
  a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
301
  if (lg(a) < 3) pari_err_PRIME("rootmod",p);
302
  a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
303
  a = FpX_gcd(f,a, p);
304
  if (!degpol(a)) return ZC_copy(S.done);
305
  split_add(&S, FpX_normalize(a,p));
306

307
  q = shifti(p,-1);
308
  pol0 = icopy(gen_1); /* constant term, will vary in place */
309
  pol = deg1pol_shallow(gen_1, pol0, varn(f));
310
  for (pol0[2] = 1;; pol0[2]++)
311
  {
312
    long j, l = lg(S.todo);
313
    if (l == 1) return sort(S.done);
314
    if (pol0[2] == 100 && !BPSW_psp(p)) pari_err_PRIME("polrootsmod",p);
315
    for (j = 1; j < l; j++)
316
    {
317
      GEN b, r, s, c = gel(S.todo,j);
318
      switch(degpol(c))
319
      { /* convert linear and quadratics to roots, try to split the rest */
320
        case 1:
321
          split_moveto_done(&S, j, subii(p, gel(c,2)));
322
          j--; l--; break;
323
        case 2:
324
          r = FpX_quad_root(c, p, 0);
325
          if (!r) pari_err_PRIME("polrootsmod",p);
326
          s = FpX_otherroot(c,r, p);
327
          split_done(&S, j, r, s);
328
          j--; l--; break;
329
        default:
330
          b = FpXQ_pow(pol,q, c,p);
331
          if (degpol(b) <= 0) continue;
332
          b = FpX_gcd(c,FpX_Fp_sub_shallow(b,gen_1,p), p);
333
          if (!degpol(b)) continue;
334
          b = FpX_normalize(b, p);
335
          c = FpX_div(c,b, p);
336
          split_todo(&S, j, b, c);
337
      }
338
    }
339
  }
340
}
341

342
/* Assume f is normalized */
343
static ulong
344
Flx_cubic_root(GEN ff, ulong p)
345
{
346
  GEN f = Flx_normalize(ff,p);
347
  ulong pi = get_Fl_red(p);
348
  ulong a = f[4], b=f[3], c=f[2], p3 = p%3==1 ? (2*p+1)/3 :(p+1)/3;
349
  ulong t = Fl_mul_pre(a, p3, p, pi), t2 = Fl_sqr_pre(t, p, pi);
350
  ulong A = Fl_sub(b, Fl_triple(t2, p), p);
351
  ulong B = Fl_addmul_pre(c, t, Fl_sub(Fl_double(t2, p), b, p), p, pi);
352
  ulong A3 =  Fl_mul_pre(A, p3, p, pi);
353
  ulong A32 = Fl_sqr_pre(A3, p, pi), A33 = Fl_mul_pre(A3, A32, p, pi);
354
  ulong S = Fl_neg(B,p), P = Fl_neg(A3,p);
355
  ulong D = Fl_add(Fl_sqr_pre(S, p, pi), Fl_double(Fl_double(A33, p), p), p);
356
  if (krouu(D,p) >= 0)
357
  {
358
    ulong s = Fl_sqrt_pre(D, p, pi), vS1, vS2;
359
    ulong S1 = S==s ? S: Fl_halve(Fl_sub(S, s, p), p);
360
    if (p%3==2) /* 1 solutions */
361
      vS1 = Fl_powu_pre(S1, (2*p-1)/3, p, pi);
362
    else
363
    {
364
      vS1 = Fl_sqrtl_pre(S1, 3, p, pi);
365
      if (vS1==~0UL) return p; /*0 solutions*/
366
      /*3 solutions*/
367
    }
368
    vS2 = P? Fl_mul_pre(P, Fl_inv(vS1, p), p, pi): 0;
369
    return Fl_sub(Fl_add(vS1,vS2, p), t, p);
370
  }
371
  else
372
  {
373
    pari_sp av = avma;
374
    GEN S1 = mkvecsmall2(Fl_halve(S, p), Fl_halve(1UL, p));
375
    GEN vS1 = Fl2_sqrtn_pre(S1, utoi(3), D, p, pi, NULL);
376
    ulong Sa;
377
    if (!vS1) return p; /*0 solutions, p%3==2*/
378
    Sa = vS1[1];
379
    if (p%3==1) /*1 solutions*/
380
    {
381
      ulong Fa = Fl2_norm_pre(vS1, D, p, pi);
382
      if (Fa!=P)
383
        Sa = Fl_mul(Sa, Fl_div(Fa, P, p),p);
384
    }
385
    set_avma(av);
386
    return Fl_sub(Fl_double(Sa,p),t,p);
387
  }
388
}
389

390
/* Assume f is normalized */
391
static GEN
392
FpX_cubic_root(GEN ff, GEN p)
393
{
394
  GEN f = FpX_normalize(ff,p);
395
  GEN a = gel(f,4), b = gel(f,3), c = gel(f,2);
396
  ulong pm3 = umodiu(p,3);
397
  GEN p3 = pm3==1 ? diviuexact(addiu(shifti(p,1),1),3)
398
                  : diviuexact(addiu(p,1),3);
399
  GEN t = Fp_mul(a, p3, p), t2 = Fp_sqr(t, p);
400
  GEN A = Fp_sub(b, Fp_mulu(t2, 3, p), p);
401
  GEN B = Fp_addmul(c, t, Fp_sub(shifti(t2, 1), b, p), p);
402
  GEN A3 =  Fp_mul(A, p3, p);
403
  GEN A32 = Fp_sqr(A3, p), A33 = Fp_mul(A3, A32, p);
404
  GEN S = Fp_neg(B,p), P = Fp_neg(A3,p);
405
  GEN D = Fp_add(Fp_sqr(S, p), shifti(A33, 2), p);
406
  if (kronecker(D,p) >= 0)
407
  {
408
    GEN s = Fp_sqrt(D, p), vS1, vS2;
409
    GEN S1 = S==s ? S: Fp_halve(Fp_sub(S, s, p), p);
410
    if (pm3 == 2) /* 1 solutions */
411
      vS1 = Fp_pow(S1, diviuexact(addiu(shifti(p, 1), 1), 3), p);
412
    else
413
    {
414
      vS1 = Fp_sqrtn(S1, utoi(3), p, NULL);
415
      if (!vS1) return p; /*0 solutions*/
416
      /*3 solutions*/
417
    }
418
    vS2 = P? Fp_mul(P, Fp_inv(vS1, p), p): 0;
419
    return Fp_sub(Fp_add(vS1,vS2, p), t, p);
420
  }
421
  else
422
  {
423
    pari_sp av = avma;
424
    GEN T = deg2pol_shallow(gen_1, gen_0, negi(D), 0);
425
    GEN S1 = deg1pol_shallow(Fp_halve(gen_1, p), Fp_halve(S, p), 0);
426
    GEN vS1 = FpXQ_sqrtn(S1, utoi(3), T, p, NULL);
427
    GEN Sa;
428
    if (!vS1) return p; /*0 solutions, p%3==2*/
429
    Sa = gel(vS1,2);
430
    if (pm3 == 1) /*1 solutions*/
431
    {
432
      GEN Fa = FpXQ_norm(vS1, T, p);
433
      if (!equalii(Fa,P))
434
        Sa = Fp_mul(Sa, Fp_div(Fa, P, p),p);
435
    }
436
    set_avma(av);
437
    return Fp_sub(shifti(Sa,1),t,p);
438
  }
439
}
440

441
/* assume p > 2 prime */
442
static ulong
443
Flx_oneroot_i(GEN f, ulong p, long fl)
444
{
445
  GEN pol, a;
446
  ulong q;
447
  long da;
448

449
  if (Flx_val(f)) return 0;
450
  switch(degpol(f))
451
  {
452
    case 1: return Fl_neg(f[2], p);
453
    case 2: return Flx_quad_root(f, p, 1);
454
    case 3: if (p>3) return Flx_cubic_root(f, p); /*FALL THROUGH*/
455
  }
456

457
  if (!fl)
458
  {
459
    a = Flxq_powu(polx_Flx(f[1]), p - 1, f,p);
460
    if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
461
    a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod f */
462
    a = Flx_gcd(f,a, p);
463
  } else a = f;
464
  da = degpol(a);
465
  if (!da) return p;
466
  a = Flx_normalize(a,p);
467

468
  q = p >> 1;
469
  pol = polx_Flx(f[1]);
470
  for(pol[2] = 1;; pol[2]++)
471
  {
472
    if (pol[2] == 1000 && !uisprime(p)) pari_err_PRIME("Flx_oneroot",utoipos(p));
473
    switch(da)
474
    {
475
      case 1: return Fl_neg(a[2], p);
476
      case 2: return Flx_quad_root(a, p, 0);
477
      case 3: if (p>3) return Flx_cubic_root(a, p); /*FALL THROUGH*/
478
      default: {
479
        GEN b = Flxq_powu(pol,q, a,p);
480
        long db;
481
        if (degpol(b) <= 0) continue;
482
        b = Flx_gcd(a,Flx_Fl_add(b,p-1,p), p);
483
        db = degpol(b); if (!db) continue;
484
        b = Flx_normalize(b, p);
485
        if (db <= (da >> 1)) {
486
          a = b;
487
          da = db;
488
        } else {
489
          a = Flx_div(a,b, p);
490
          da -= db;
491
        }
492
      }
493
    }
494
  }
495
}
496

497
/* assume p > 3 prime */
498
static GEN
499
FpX_oneroot_i(GEN f, GEN p)
500
{
501
  GEN pol, pol0, a, q;
502
  long da;
503

504
  if (ZX_val(f)) return gen_0;
505
  switch(degpol(f))
506
  {
507
    case 1: return subii(p, gel(f,2));
508
    case 2: return FpX_quad_root(f, p, 1);
509
    case 3: return FpX_cubic_root(f, p);
510
  }
511

512
  a = FpXQ_pow(pol_x(varn(f)), subiu(p,1), f,p);
513
  if (lg(a) < 3) pari_err_PRIME("rootmod",p);
514
  a = FpX_Fp_sub_shallow(a, gen_1, p); /* a = x^(p-1) - 1 mod f */
515
  a = FpX_gcd(f,a, p);
516
  da = degpol(a);
517
  if (!da) return NULL;
518
  a = FpX_normalize(a,p);
519

520
  q = shifti(p,-1);
521
  pol0 = icopy(gen_1); /* constant term, will vary in place */
522
  pol = deg1pol_shallow(gen_1, pol0, varn(f));
523
  for (pol0[2]=1; ; pol0[2]++)
524
  {
525
    if (pol0[2] == 1000 && !BPSW_psp(p)) pari_err_PRIME("FpX_oneroot",p);
526
    switch(da)
527
    {
528
      case 1: return subii(p, gel(a,2));
529
      case 2: return FpX_quad_root(a, p, 0);
530
      default: {
531
        GEN b = FpXQ_pow(pol,q, a,p);
532
        long db;
533
        if (degpol(b) <= 0) continue;
534
        b = FpX_gcd(a,FpX_Fp_sub_shallow(b,gen_1,p), p);
535
        db = degpol(b); if (!db) continue;
536
        b = FpX_normalize(b, p);
537
        if (db <= (da >> 1)) {
538
          a = b;
539
          da = db;
540
        } else {
541
          a = FpX_div(a,b, p);
542
          da -= db;
543
        }
544
      }
545
    }
546
  }
547
}
548

549
ulong
550
Flx_oneroot(GEN f, ulong p)
551
{
552
  pari_sp av = avma;
553
  ulong r;
554
  switch(lg(f))
555
  {
556
    case 2: return 0;
557
    case 3: return p;
558
  }
559
  if (p == 2) return Flx_oneroot_mod_2(f);
560
  r = Flx_oneroot_i(Flx_normalize(f, p), p, 0);
561
  return gc_ulong(av,r);
562
}
563

564
ulong
565
Flx_oneroot_split(GEN f, ulong p)
566
{
567
  pari_sp av = avma;
568
  ulong r;
569
  switch(lg(f))
570
  {
571
    case 2: return 0;
572
    case 3: return p;
573
  }
574
  if (p == 2) return Flx_oneroot_mod_2(f);
575
  r = Flx_oneroot_i(Flx_normalize(f, p), p, 1);
576
  return gc_ulong(av,r);
577
}
578

579
/* assume that p is prime */
580
GEN
581
FpX_oneroot(GEN f, GEN pp) {
582
  pari_sp av = avma;
583
  ZX_rootmod_init(&f, pp);
584
  switch(lg(f))
585
  {
586
    case 2: set_avma(av); return gen_0;
587
    case 3: return gc_NULL(av);
588
  }
589
  if (typ(f) == t_VECSMALL)
590
  {
591
    ulong r, p = pp[2];
592
    if (p == 2)
593
      r = Flx_oneroot_mod_2(f);
594
    else
595
      r = Flx_oneroot_i(f, p, 0);
596
    set_avma(av);
597
    return (r == p)? NULL: utoi(r);
598
  }
599
  f = FpX_oneroot_i(f, pp);
600
  if (!f) return gc_NULL(av);
601
  return gerepileuptoint(av, f);
602
}
603

604
/* returns a root of unity in F_p that is suitable for finding a factor   */
605
/* of degree deg_factor of a polynomial of degree deg; the order is       */
606
/* returned in n                                                          */
607
/* A good choice seems to be n close to deg/deg_factor; we choose n       */
608
/* twice as big and decrement until it divides p-1.                       */
609
static GEN
610
good_root_of_unity(GEN p, long deg, long deg_factor, long *pt_n)
611
{
612
   pari_sp ltop = avma;
613
   GEN pm, factn, power, base, zeta;
614
   long n;
615

616
   pm = subis (p, 1ul);
617
   for (n = deg / 2 / deg_factor + 1; !dvdiu (pm, n); n--);
618
   factn = Z_factor(stoi(n));
619
   power = diviuexact (pm, n);
620
   base = gen_1;
621
   do {
622
      base = addis (base, 1l);
623
      zeta = Fp_pow (base, power, p);
624
   }
625
   while (!equaliu (Fp_order (zeta, factn, p), n));
626
   *pt_n = n;
627
   return gerepileuptoint (ltop, zeta);
628
}
629

630
GEN
631
FpX_oneroot_split(GEN fact, GEN p)
632
{
633
  pari_sp av = avma;
634
  long n, deg_f, i, dmin;
635
  GEN prim, expo, minfactor, xplusa, zeta, xpow;
636
  fact = FpX_normalize(fact, p);
637
  deg_f = degpol(fact);
638
  if (deg_f <= 3) return FpX_oneroot(fact, p);
639
  minfactor = fact; /* factor of minimal degree found so far */
640
  dmin = degpol(minfactor);
641
  xplusa = pol_x(varn(fact));
642
  while (dmin > 3)
643
  {
644
    /* split minfactor by computing its gcd with (X+a)^exp-zeta, where    */
645
    /* zeta varies over the roots of unity in F_p                         */
646
    fact = minfactor; deg_f = dmin;
647
    zeta = gen_1;
648
    prim = good_root_of_unity(p, deg_f, 1, &n);
649
    expo = diviuexact(subiu(p, 1), n);
650
    /* update X+a, avoid a=0 */
651
    gel (xplusa, 2) = addis (gel (xplusa, 2), 1);
652
    xpow = FpXQ_pow (xplusa, expo, fact, p);
653
    for (i = 0; i < n; i++)
654
    {
655
      GEN tmp = FpX_gcd(FpX_Fp_sub(xpow, zeta, p), fact, p);
656
      long dtmp = degpol(tmp);
657
      if (dtmp > 0 && dtmp < deg_f)
658
      {
659
        fact = FpX_div(fact, tmp, p); deg_f = degpol(fact);
660
        if (dtmp < dmin)
661
        {
662
          minfactor = FpX_normalize (tmp, p);
663
          dmin = dtmp;
664
          if (dmin == 1 || dmin <= (2 * deg_f) / n - 1)
665
            /* stop early to avoid too many gcds */
666
            break;
667
        }
668
      }
669
      zeta = Fp_mul (zeta, prim, p);
670
    }
671
  }
672
  return gerepileuptoint(av, FpX_oneroot(minfactor, p));
673
}
674

675
/*******************************************************************/
676
/*                                                                 */
677
/*                     FACTORISATION MODULO p                      */
678
/*                                                                 */
679
/*******************************************************************/
680

681
/* F / E  a vector of vectors of factors / exponents of virtual length l
682
 * (their real lg may be larger). Set their lg to j, concat and return [F,E] */
683
static GEN
684
FE_concat(GEN F, GEN E, long l)
685
{
686
  setlg(E,l); E = shallowconcat1(E);
687
  setlg(F,l); F = shallowconcat1(F); return mkvec2(F,E);
688
}
689

690
static GEN
691
ddf_to_ddf2_i(GEN V, long fl)
692
{
693
  GEN F, D;
694
  long i, j, l = lg(V);
695
  F = cgetg(l, t_VEC);
696
  D = cgetg(l, t_VECSMALL);
697
  for (i = j = 1; i < l; i++)
698
  {
699
    GEN Vi = gel(V,i);
700
    if ((fl==2 && F2x_degree(Vi) == 0)
701
      ||(fl==0 && degpol(Vi) == 0)) continue;
702
    gel(F,j) = Vi;
703
    uel(D,j) = i; j++;
704
  }
705
  setlg(F,j);
706
  setlg(D,j); return mkvec2(F,D);
707
}
708

709
GEN
710
ddf_to_ddf2(GEN V)
711
{ return ddf_to_ddf2_i(V, 0); }
712

713
static GEN
714
F2x_ddf_to_ddf2(GEN V)
715
{ return ddf_to_ddf2_i(V, 2); }
716

717
GEN
718
vddf_to_simplefact(GEN V, long d)
719
{
720
  GEN E, F;
721
  long i, j, c, l = lg(V);
722
  F = cgetg(d+1, t_VECSMALL);
723
  E = cgetg(d+1, t_VECSMALL);
724
  for (i = c = 1; i < l; i++)
725
  {
726
    GEN Vi = gel(V,i);
727
    long l = lg(Vi);
728
    for (j = 1; j < l; j++)
729
    {
730
      long k, n = degpol(gel(Vi,j)) / j;
731
      for (k = 1; k <= n; k++) { uel(F,c) = j; uel(E,c) = i; c++; }
732
    }
733
  }
734
  setlg(F,c);
735
  setlg(E,c);
736
  return sort_factor(mkvec2(F,E), (void*)&cmpGuGu, cmp_nodata);
737
}
738

739
/* product of terms of degree 1 in factorization of f */
740
GEN
741
FpX_split_part(GEN f, GEN p)
742
{
743
  long n = degpol(f);
744
  GEN z, X = pol_x(varn(f));
745
  if (n <= 1) return f;
746
  f = FpX_red(f, p);
747
  z = FpX_sub(FpX_Frobenius(f, p), X, p);
748
  return FpX_gcd(z,f,p);
749
}
750

751
/* Compute the number of roots in Fp without counting multiplicity
752
 * return -1 for 0 polynomial. lc(f) must be prime to p. */
753
long
754
FpX_nbroots(GEN f, GEN p)
755
{
756
  pari_sp av = avma;
757
  GEN z = FpX_split_part(f, p);
758
  return gc_long(av, degpol(z));
759
}
760

761
/* 1 < deg(f) <= p */
762
static int
763
Flx_is_totally_split_i(GEN f, ulong p)
764
{
765
  GEN F = Flx_Frobenius(f, p);
766
  return degpol(F)==1 && uel(F,2)==0UL && uel(F,3)==1UL;
767
}
768
int
769
Flx_is_totally_split(GEN f, ulong p)
770
{
771
  pari_sp av = avma;
772
  ulong n = degpol(f);
773
  if (n <= 1) return 1;
774
  if (n > p) return 0; /* includes n < 0 */
775
  return gc_bool(av, Flx_is_totally_split_i(f,p));
776
}
777
int
778
FpX_is_totally_split(GEN f, GEN p)
779
{
780
  pari_sp av = avma;
781
  ulong n = degpol(f);
782
  int u;
783
  if (n <= 1) return 1;
784
  if (abscmpui(n, p) > 0) return 0; /* includes n < 0 */
785
  if (lgefint(p) != 3)
786
    u = gequalX(FpX_Frobenius(FpX_red(f,p), p));
787
  else
788
  {
789
    ulong pp = (ulong)p[2];
790
    u = Flx_is_totally_split_i(ZX_to_Flx(f,pp), pp);
791
  }
792
  return gc_bool(av, u);
793
}
794

795
long
796
Flx_nbroots(GEN f, ulong p)
797
{
798
  long n = degpol(f);
799
  pari_sp av = avma;
800
  GEN z;
801
  if (n <= 1) return n;
802
  if (n == 2)
803
  {
804
    ulong D;
805
    if (p==2) return (f[2]==0) + (f[2]!=f[3]);
806
    D = Fl_sub(Fl_sqr(f[3], p), Fl_mul(Fl_mul(f[4], f[2], p), 4%p, p), p);
807
    return 1 + krouu(D,p);
808
  }
809
  z = Flx_sub(Flx_Frobenius(f, p), polx_Flx(f[1]), p);
810
  z = Flx_gcd(z, f, p);
811
  return gc_long(av, degpol(z));
812
}
813

814
long
815
FpX_ddf_degree(GEN T, GEN XP, GEN p)
816
{
817
  pari_sp av = avma;
818
  GEN X, b, g, xq;
819
  long i, j, n, v, B, l, m;
820
  pari_timer ti;
821
  hashtable h;
822

823
  n = get_FpX_degree(T); v = get_FpX_var(T);
824
  X = pol_x(v);
825
  if (ZX_equal(X,XP)) return 1;
826
  B = n/2;
827
  l = usqrt(B);
828
  m = (B+l-1)/l;
829
  T = FpX_get_red(T, p);
830
  hash_init_GEN(&h, l+2, ZX_equal, 1);
831
  hash_insert_long(&h, X,  0);
832
  hash_insert_long(&h, XP, 1);
833
  if (DEBUGLEVEL>=7) timer_start(&ti);
834
  b = XP;
835
  xq = FpXQ_powers(b, brent_kung_optpow(n, l-1, 1),  T, p);
836
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq baby");
837
  for (i = 3; i <= l+1; i++)
838
  {
839
    b = FpX_FpXQV_eval(b, xq, T, p);
840
    if (gequalX(b)) return gc_long(av,i-1);
841
    hash_insert_long(&h, b, i-1);
842
  }
843
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: baby");
844
  g = b;
845
  xq = FpXQ_powers(g, brent_kung_optpow(n, m, 1),  T, p);
846
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_degree: xq giant");
847
  for(i = 2; i <= m+1; i++)
848
  {
849
    g = FpX_FpXQV_eval(g, xq, T, p);
850
    if (hash_haskey_long(&h, g, &j)) return gc_long(av, l*i-j);
851
  }
852
  return gc_long(av,n);
853
}
854

855
/* See <http://www.shoup.net/papers/factorimpl.pdf> */
856
static GEN
857
FpX_ddf_Shoup(GEN T, GEN XP, GEN p)
858
{
859
  GEN b, g, h, F, f, Tr, xq;
860
  long i, j, n, v, B, l, m;
861
  pari_timer ti;
862

863
  n = get_FpX_degree(T); v = get_FpX_var(T);
864
  if (n == 0) return cgetg(1, t_VEC);
865
  if (n == 1) return mkvec(get_FpX_mod(T));
866
  B = n/2;
867
  l = usqrt(B);
868
  m = (B+l-1)/l;
869
  T = FpX_get_red(T, p);
870
  b = cgetg(l+2, t_VEC);
871
  gel(b, 1) = pol_x(v);
872
  gel(b, 2) = XP;
873
  if (DEBUGLEVEL>=7) timer_start(&ti);
874
  xq = FpXQ_powers(gel(b, 2), brent_kung_optpow(n, l-1, 1),  T, p);
875
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq baby");
876
  for (i = 3; i <= l+1; i++)
877
    gel(b, i) = FpX_FpXQV_eval(gel(b, i-1), xq, T, p);
878
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: baby");
879
  xq = FpXQ_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
880
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: xq giant");
881
  g = cgetg(m+1, t_VEC);
882
  gel(g, 1) = gel(xq, 2);
883
  for(i = 2; i <= m; i++) gel(g, i) = FpX_FpXQV_eval(gel(g, i-1), xq, T, p);
884
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: giant");
885
  h = cgetg(m+1, t_VEC);
886
  for (j = 1; j <= m; j++)
887
  {
888
    pari_sp av = avma;
889
    GEN gj = gel(g,j), e = FpX_sub(gj, gel(b,1), p);
890
    for (i = 2; i <= l; i++) e = FpXQ_mul(e, FpX_sub(gj, gel(b,i), p), T, p);
891
    gel(h,j) = gerepileupto(av, e);
892
  }
893
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: diff");
894
  Tr = get_FpX_mod(T);
895
  F = cgetg(m+1, t_VEC);
896
  for (j = 1; j <= m; j++)
897
  {
898
    GEN u = FpX_gcd(Tr, gel(h,j), p);
899
    if (degpol(u))
900
    {
901
      u = FpX_normalize(u, p);
902
      Tr = FpX_div(Tr, u, p);
903
    }
904
    gel(F,j) = u;
905
  }
906
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: F");
907
  f = const_vec(n, pol_1(v));
908
  for (j = 1; j <= m; j++)
909
  {
910
    GEN e = gel(F, j);
911
    for (i=l-1; i >= 0; i--)
912
    {
913
      GEN u = FpX_gcd(e, FpX_sub(gel(g, j), gel(b, i+1), p), p);
914
      if (degpol(u))
915
      {
916
        u = FpX_normalize(u, p);
917
        gel(f, l*j-i) = u;
918
        e = FpX_div(e, u, p);
919
      }
920
      if (!degpol(e)) break;
921
    }
922
  }
923
  if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_ddf_Shoup: f");
924
  if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
925
  return f;
926
}
927

928
static void
929
FpX_edf_simple(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
930
{
931
  long n = degpol(Tp), r = n/d, ct = 0;
932
  GEN T, f, ff, p2;
933
  if (r==1) { gel(V, idx) = Tp; return; }
934
  p2 = shifti(p,-1);
935
  T = FpX_get_red(Tp, p);
936
  XP = FpX_rem(XP, T, p);
937
  while (1)
938
  {
939
    pari_sp btop = avma;
940
    long i;
941
    GEN g = random_FpX(n, varn(Tp), p);
942
    GEN t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
943
    if (signe(t) == 0) continue;
944
    for(i=1; i<=10; i++)
945
    {
946
      pari_sp btop2 = avma;
947
      GEN R = FpXQ_pow(FpX_Fp_add(t, randomi(p), p), p2, T, p);
948
      f = FpX_gcd(FpX_Fp_sub(R, gen_1, p), Tp, p);
949
      if (degpol(f) > 0 && degpol(f) < n) break;
950
      set_avma(btop2);
951
    }
952
    if (degpol(f) > 0 && degpol(f) < n) break;
953
    if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_simple",p);
954
    set_avma(btop);
955
  }
956
  f = FpX_normalize(f, p);
957
  ff = FpX_div(Tp, f ,p);
958
  FpX_edf_simple(f, XP, d, p, V, idx);
959
  FpX_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
960
}
961

962
static void
963
FpX_edf_rec(GEN T, GEN hp, GEN t, long d, GEN p2, GEN p, GEN V, long idx)
964
{
965
  pari_sp av;
966
  GEN Tp = get_FpX_mod(T);
967
  long n = degpol(hp), vT = varn(Tp), ct = 0;
968
  GEN u1, u2, f1, f2, R, h;
969
  h = FpX_get_red(hp, p);
970
  t = FpX_rem(t, T, p);
971
  av = avma;
972
  do
973
  {
974
    set_avma(av);
975
    R = FpXQ_pow(deg1pol(gen_1, randomi(p), vT), p2, h, p);
976
    u1 = FpX_gcd(FpX_Fp_sub(R, gen_1, p), hp, p);
977
    if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf_rec",p);
978
  } while (degpol(u1)==0 || degpol(u1)==n);
979
  f1 = FpX_gcd(FpX_FpXQ_eval(u1, t, T, p), Tp, p);
980
  f1 = FpX_normalize(f1, p);
981
  u2 = FpX_div(hp, u1, p);
982
  f2 = FpX_div(Tp, f1, p);
983
  if (degpol(u1)==1)
984
    gel(V, idx) = f1;
985
  else
986
    FpX_edf_rec(FpX_get_red(f1, p), u1, t, d, p2, p, V, idx);
987
  idx += degpol(f1)/d;
988
  if (degpol(u2)==1)
989
    gel(V, idx) = f2;
990
  else
991
    FpX_edf_rec(FpX_get_red(f2, p), u2, t, d, p2, p, V, idx);
992
}
993

994
/* assume Tp a squarefree product of r > 1 irred. factors of degree d */
995
static void
996
FpX_edf(GEN Tp, GEN XP, long d, GEN p, GEN V, long idx)
997
{
998
  long n = degpol(Tp), r = n/d, vT = varn(Tp), ct = 0;
999
  GEN T, h, t;
1000
  pari_timer ti;
1001

1002
  T = FpX_get_red(Tp, p);
1003
  XP = FpX_rem(XP, T, p);
1004
  if (DEBUGLEVEL>=7) timer_start(&ti);
1005
  do
1006
  {
1007
    GEN g = random_FpX(n, vT, p);
1008
    t = gel(FpXQ_auttrace(mkvec2(XP, g), d, T, p), 2);
1009
    if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_auttrace");
1010
    h = FpXQ_minpoly(t, T, p);
1011
    if (DEBUGLEVEL>=7) timer_printf(&ti,"FpX_edf: FpXQ_minpoly");
1012
    if (++ct == 10 && !BPSW_psp(p)) pari_err_PRIME("FpX_edf",p);
1013
  } while (degpol(h) != r);
1014
  FpX_edf_rec(T, h, t, d, shifti(p, -1), p, V, idx);
1015
}
1016

1017
static GEN
1018
FpX_factor_Shoup(GEN T, GEN p)
1019
{
1020
  long i, n, s = 0;
1021
  GEN XP, D, V;
1022
  long e = expi(p);
1023
  pari_timer ti;
1024
  n = get_FpX_degree(T);
1025
  T = FpX_get_red(T, p);
1026
  if (DEBUGLEVEL>=6) timer_start(&ti);
1027
  XP = FpX_Frobenius(T, p);
1028
  if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1029
  D = FpX_ddf_Shoup(T, XP, p);
1030
  if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1031
  s = ddf_to_nbfact(D);
1032
  V = cgetg(s+1, t_COL);
1033
  for (i = 1, s = 1; i <= n; i++)
1034
  {
1035
    GEN Di = gel(D,i);
1036
    long ni = degpol(Di), ri = ni/i;
1037
    if (ni == 0) continue;
1038
    Di = FpX_normalize(Di, p);
1039
    if (ni == i) { gel(V, s++) = Di; continue; }
1040
    if (ri <= e*expu(e))
1041
      FpX_edf(Di, XP, i, p, V, s);
1042
    else
1043
      FpX_edf_simple(Di, XP, i, p, V, s);
1044
    if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_edf(%ld)",i);
1045
    s += ri;
1046
  }
1047
  return V;
1048
}
1049

1050
long
1051
ddf_to_nbfact(GEN D)
1052
{
1053
  long l = lg(D), i, s = 0;
1054
  for(i = 1; i < l; i++) s += degpol(gel(D,i))/i;
1055
  return s;
1056
}
1057

1058
/* Yun algorithm: Assume p > degpol(T) */
1059
static GEN
1060
FpX_factor_Yun(GEN T, GEN p)
1061
{
1062
  long n = degpol(T), i = 1;
1063
  GEN a, b, c, d = FpX_deriv(T, p);
1064
  GEN V = cgetg(n+1,t_VEC);
1065
  a = FpX_gcd(T, d, p);
1066
  if (degpol(a) == 0) return mkvec(T);
1067
  b = FpX_div(T, a, p);
1068
  do
1069
  {
1070
    c = FpX_div(d, a, p);
1071
    d = FpX_sub(c, FpX_deriv(b, p), p);
1072
    a = FpX_normalize(FpX_gcd(b, d, p), p);
1073
    gel(V, i++) = a;
1074
    b = FpX_div(b, a, p);
1075
  } while (degpol(b));
1076
  setlg(V, i); return V;
1077
}
1078
GEN
1079
FpX_factor_squarefree(GEN T, GEN p)
1080
{
1081
  if (lgefint(p)==3)
1082
  {
1083
    ulong pp = (ulong)p[2];
1084
    GEN u = Flx_factor_squarefree(ZX_to_Flx(T,pp), pp);
1085
    return FlxV_to_ZXV(u);
1086
  }
1087
  return FpX_factor_Yun(T, p);
1088
}
1089

1090
long
1091
FpX_ispower(GEN f, ulong k, GEN p, GEN *pt_r)
1092
{
1093
  pari_sp av = avma;
1094
  GEN lc, F;
1095
  long i, l, n = degpol(f), v = varn(f);
1096
  if (n % k) return 0;
1097
  if (lgefint(p)==3)
1098
  {
1099
    ulong pp = p[2];
1100
    GEN fp = ZX_to_Flx(f, pp);
1101
    if (!Flx_ispower(fp, k, pp, pt_r)) return gc_long(av,0);
1102
    if (pt_r) *pt_r = gerepileupto(av, Flx_to_ZX(*pt_r)); else set_avma(av);
1103
    return 1;
1104
  }
1105
  lc = Fp_sqrtn(leading_coeff(f), stoi(k), p, NULL);
1106
  if (!lc) { av = avma; return 0; }
1107
  F = FpX_factor_Yun(f, p); l = lg(F)-1;
1108
  for(i=1; i <= l; i++)
1109
    if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1110
  if (pt_r)
1111
  {
1112
    GEN r = scalarpol(lc, v), s = pol_1(v);
1113
    for (i=l; i>=1; i--)
1114
    {
1115
      if (i%k) continue;
1116
      s = FpX_mul(s, gel(F,i), p);
1117
      r = FpX_mul(r, s, p);
1118
    }
1119
    *pt_r = gerepileupto(av, r);
1120
  } else av = avma;
1121
  return 1;
1122
}
1123

1124
static GEN
1125
FpX_factor_Cantor(GEN T, GEN p)
1126
{
1127
  GEN E, F, V = FpX_factor_Yun(T, p);
1128
  long i, j, l = lg(V);
1129
  F = cgetg(l, t_VEC);
1130
  E = cgetg(l, t_VEC);
1131
  for (i=1, j=1; i < l; i++)
1132
    if (degpol(gel(V,i)))
1133
    {
1134
      GEN Fj = FpX_factor_Shoup(gel(V,i), p);
1135
      gel(F, j) = Fj;
1136
      gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1137
      j++;
1138
    }
1139
  return sort_factor_pol(FE_concat(F,E,j), cmpii);
1140
}
1141

1142
static GEN
1143
FpX_ddf_i(GEN T, GEN p)
1144
{
1145
  GEN XP;
1146
  T = FpX_get_red(T, p);
1147
  XP = FpX_Frobenius(T, p);
1148
  return ddf_to_ddf2(FpX_ddf_Shoup(T, XP, p));
1149
}
1150

1151
GEN
1152
FpX_ddf(GEN f, GEN p)
1153
{
1154
  pari_sp av = avma;
1155
  GEN F;
1156
  switch(ZX_factmod_init(&f, p))
1157
  {
1158
    case 0:  F = F2x_ddf(f);
1159
             F2xV_to_ZXV_inplace(gel(F,1)); break;
1160
    case 1:  F = Flx_ddf(f,p[2]);
1161
             FlxV_to_ZXV_inplace(gel(F,1)); break;
1162
    default: F = FpX_ddf_i(f,p); break;
1163
  }
1164
  return gerepilecopy(av, F);
1165
}
1166

1167
static GEN Flx_simplefact_Cantor(GEN T, ulong p);
1168
static GEN
1169
FpX_simplefact_Cantor(GEN T, GEN p)
1170
{
1171
  GEN V;
1172
  long i, l;
1173
  if (lgefint(p) == 3)
1174
  {
1175
    ulong pp = p[2];
1176
    return Flx_simplefact_Cantor(ZX_to_Flx(T,pp), pp);
1177
  }
1178
  T = FpX_get_red(T, p);
1179
  V = FpX_factor_Yun(get_FpX_mod(T), p); l = lg(V);
1180
  for (i=1; i < l; i++)
1181
    gel(V,i) = FpX_ddf_Shoup(gel(V,i), FpX_Frobenius(gel(V,i), p), p);
1182
  return vddf_to_simplefact(V, get_FpX_degree(T));
1183
}
1184

1185
static int
1186
FpX_isirred_Cantor(GEN Tp, GEN p)
1187
{
1188
  pari_sp av = avma;
1189
  pari_timer ti;
1190
  long n;
1191
  GEN T = get_FpX_mod(Tp);
1192
  GEN dT = FpX_deriv(T, p);
1193
  GEN XP, D;
1194
  if (degpol(FpX_gcd(T, dT, p)) != 0) return gc_bool(av,0);
1195
  n = get_FpX_degree(T);
1196
  T = FpX_get_red(Tp, p);
1197
  if (DEBUGLEVEL>=6) timer_start(&ti);
1198
  XP = FpX_Frobenius(T, p);
1199
  if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_Frobenius");
1200
  D = FpX_ddf_Shoup(T, XP, p);
1201
  if (DEBUGLEVEL>=6) timer_printf(&ti,"FpX_ddf_Shoup");
1202
  return gc_bool(av, degpol(gel(D,n)) == n);
1203
}
1204

1205
static GEN FpX_factor_deg2(GEN f, GEN p, long d, long flag);
1206

1207
/*Assume that p is large and odd*/
1208
static GEN
1209
FpX_factor_i(GEN f, GEN pp, long flag)
1210
{
1211
  long d = degpol(f);
1212
  if (d <= 2) return FpX_factor_deg2(f,pp,d,flag);
1213
  switch(flag)
1214
  {
1215
    default: return FpX_factor_Cantor(f, pp);
1216
    case 1: return FpX_simplefact_Cantor(f, pp);
1217
    case 2: return FpX_isirred_Cantor(f, pp)? gen_1: NULL;
1218
  }
1219
}
1220

1221
long
1222
FpX_nbfact_Frobenius(GEN T, GEN XP, GEN p)
1223
{
1224
  pari_sp av = avma;
1225
  long s = ddf_to_nbfact(FpX_ddf_Shoup(T, XP, p));
1226
  return gc_long(av,s);
1227
}
1228

1229
long
1230
FpX_nbfact(GEN T, GEN p)
1231
{
1232
  pari_sp av = avma;
1233
  GEN XP = FpX_Frobenius(T, p);
1234
  long n = FpX_nbfact_Frobenius(T, XP, p);
1235
  return gc_long(av,n);
1236
}
1237

1238
/* p > 2 */
1239
static GEN
1240
FpX_is_irred_2(GEN f, GEN p, long d)
1241
{
1242
  switch(d)
1243
  {
1244
    case -1:
1245
    case 0: return NULL;
1246
    case 1: return gen_1;
1247
  }
1248
  return FpX_quad_factortype(f, p) == -1? gen_1: NULL;
1249
}
1250
/* p > 2 */
1251
static GEN
1252
FpX_degfact_2(GEN f, GEN p, long d)
1253
{
1254
  switch(d)
1255
  {
1256
    case -1:retmkvec2(mkvecsmall(-1),mkvecsmall(1));
1257
    case 0: return trivial_fact();
1258
    case 1: retmkvec2(mkvecsmall(1), mkvecsmall(1));
1259
  }
1260
  switch(FpX_quad_factortype(f, p)) {
1261
    case  1: retmkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1262
    case -1: retmkvec2(mkvecsmall(2), mkvecsmall(1));
1263
    default: retmkvec2(mkvecsmall(1), mkvecsmall(2));
1264
  }
1265
}
1266

1267
GEN
1268
prime_fact(GEN x) { retmkmat2(mkcolcopy(x), mkcol(gen_1)); }
1269
GEN
1270
trivial_fact(void) { retmkmat2(cgetg(1,t_COL), cgetg(1,t_COL)); }
1271

1272
/* not gerepile safe */
1273
static GEN
1274
FpX_factor_2(GEN f, GEN p, long d)
1275
{
1276
  GEN r, s, R, S;
1277
  long v;
1278
  int sgn;
1279
  switch(d)
1280
  {
1281
    case -1: retmkvec2(mkcol(pol_0(varn(f))), mkvecsmall(1));
1282
    case  0: retmkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1283
    case  1: retmkvec2(mkcol(f), mkvecsmall(1));
1284
  }
1285
  r = FpX_quad_root(f, p, 1);
1286
  if (!r) return mkvec2(mkcol(f), mkvecsmall(1));
1287
  v = varn(f);
1288
  s = FpX_otherroot(f, r, p);
1289
  if (signe(r)) r = subii(p, r);
1290
  if (signe(s)) s = subii(p, s);
1291
  sgn = cmpii(s, r); if (sgn < 0) swap(s,r);
1292
  R = deg1pol_shallow(gen_1, r, v);
1293
  if (!sgn) return mkvec2(mkcol(R), mkvecsmall(2));
1294
  S = deg1pol_shallow(gen_1, s, v);
1295
  return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1296
}
1297
static GEN
1298
FpX_factor_deg2(GEN f, GEN p, long d, long flag)
1299
{
1300
  switch(flag) {
1301
    case 2: return FpX_is_irred_2(f, p, d);
1302
    case 1: return FpX_degfact_2(f, p, d);
1303
    default: return FpX_factor_2(f, p, d);
1304
  }
1305
}
1306

1307
static int
1308
F2x_quad_factortype(GEN x)
1309
{ return x[2] == 7 ? -1: x[2] == 6 ? 1 :0; }
1310

1311
static GEN
1312
F2x_is_irred_2(GEN f, long d)
1313
{ return d == 1 || (d==2 && F2x_quad_factortype(f) == -1)? gen_1: NULL; }
1314

1315
static GEN
1316
F2x_degfact_2(GEN f, long d)
1317
{
1318
  if (!d) return trivial_fact();
1319
  if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1320
  switch(F2x_quad_factortype(f)) {
1321
    case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1322
    case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1323
    default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1324
  }
1325
}
1326

1327
static GEN
1328
F2x_factor_2(GEN f, long d)
1329
{
1330
  long v = f[1];
1331
  if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1332
  if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1333
  switch(F2x_quad_factortype(f))
1334
  {
1335
  case -1: return mkvec2(mkcol(f), mkvecsmall(1));
1336
  case 0:  return mkvec2(mkcol(mkvecsmall2(v,2+F2x_coeff(f,0))), mkvecsmall(2));
1337
  default: return mkvec2(mkcol2(mkvecsmall2(v,2),mkvecsmall2(v,3)), mkvecsmall2(1,1));
1338
  }
1339
}
1340
static GEN
1341
F2x_factor_deg2(GEN f, long d, long flag)
1342
{
1343
  switch(flag) {
1344
    case 2: return F2x_is_irred_2(f, d);
1345
    case 1: return F2x_degfact_2(f, d);
1346
    default: return F2x_factor_2(f, d);
1347
  }
1348
}
1349

1350
/* xt = NULL or x^(p-1)/2 mod g */
1351
static void
1352
split_squares(struct split_t *S, GEN g, ulong p, GEN xt)
1353
{
1354
  ulong q = p >> 1;
1355
  GEN a = Flx_mod_Xnm1(g, q, p); /* mod x^(p-1)/2 - 1 */
1356
  long d = degpol(a);
1357
  if (d < 0)
1358
  {
1359
    ulong i;
1360
    split_add_done(S, (GEN)1);
1361
    for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_sqr(i,p));
1362
  } else {
1363
    if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1364
    if (d)
1365
    {
1366
      if (xt) xt = Flx_Fl_add(xt, p-1, p); else xt = Flx_Xnm1(g[1], q, p);
1367
      a = Flx_gcd(a, xt, p);
1368
      if (degpol(a)) split_add(S, Flx_normalize(a, p));
1369
    }
1370
  }
1371
}
1372
static void
1373
split_nonsquares(struct split_t *S, GEN g, ulong p, GEN xt)
1374
{
1375
  ulong q = p >> 1;
1376
  GEN a = Flx_mod_Xn1(g, q, p); /* mod x^(p-1)/2 + 1 */
1377
  long d = degpol(a);
1378
  if (d < 0)
1379
  {
1380
    ulong i, z = nonsquare_Fl(p);
1381
    split_add_done(S, (GEN)z);
1382
    for (i = 2; i <= q; i++) split_add_done(S, (GEN)Fl_mul(z, Fl_sqr(i,p), p));
1383
  } else {
1384
    if (a != g) { (void)Flx_valrem(a, &a); d = degpol(a); }
1385
    if (d)
1386
    {
1387
      if (xt) xt = Flx_Fl_add(xt, 1, p); else xt = Flx_Xn1(g[1], q, p);
1388
      a = Flx_gcd(a, xt, p);
1389
      if (degpol(a)) split_add(S, Flx_normalize(a, p));
1390
    }
1391
  }
1392
}
1393
/* p > 2. f monic Flx, f(0) != 0. Add to split_t structs coprime factors
1394
 * of g = \prod_{f(a) = 0} (X - a). Return 0 when f(x) = 0 for all x in Fp* */
1395
static int
1396
split_Flx_cut_out_roots(struct split_t *S, GEN f, ulong p)
1397
{
1398
  GEN a, g = Flx_mod_Xnm1(f, p-1, p); /* f mod x^(p-1) - 1 */
1399
  long d = degpol(g);
1400
  if (d < 0) return 0;
1401
  if (g != f) { (void)Flx_valrem(g, &g); d = degpol(g); } /*kill powers of x*/
1402
  if (!d) return 1;
1403
  if ((p >> 4) <= (ulong)d)
1404
  { /* small p; split directly using x^((p-1)/2) +/- 1 */
1405
    GEN xt = ((ulong)d < (p>>1))? Flx_rem(monomial_Flx(1, p>>1, g[1]), g, p)
1406
                                : NULL;
1407
    split_squares(S, g, p, xt);
1408
    split_nonsquares(S, g, p, xt);
1409
  } else { /* large p; use x^(p-1) - 1 directly */
1410
    a = Flxq_powu(polx_Flx(f[1]), p-1, g,p);
1411
    if (lg(a) < 3) pari_err_PRIME("rootmod",utoipos(p));
1412
    a = Flx_Fl_add(a, p-1, p); /* a = x^(p-1) - 1 mod g */
1413
    g = Flx_gcd(g,a, p);
1414
    if (degpol(g)) split_add(S, Flx_normalize(g,p));
1415
  }
1416
  return 1;
1417
}
1418

1419
/* by splitting, assume p > 2 prime, deg(f) > 0, and f monic */
1420
static GEN
1421
Flx_roots_i(GEN f, ulong p)
1422
{
1423
  GEN pol, g;
1424
  long v = Flx_valrem(f, &g);
1425
  ulong q;
1426
  struct split_t S;
1427

1428
  /* optimization: test for small degree first */
1429
  switch(degpol(g))
1430
  {
1431
    case 1: {
1432
      ulong r = p - g[2];
1433
      return v? mkvecsmall2(0, r): mkvecsmall(r);
1434
    }
1435
    case 2: {
1436
      ulong r = Flx_quad_root(g, p, 1), s;
1437
      if (r == p) return v? mkvecsmall(0): cgetg(1,t_VECSMALL);
1438
      s = Flx_otherroot(g,r, p);
1439
      if (r < s)
1440
        return v? mkvecsmall3(0, r, s): mkvecsmall2(r, s);
1441
      else if (r > s)
1442
        return v? mkvecsmall3(0, s, r): mkvecsmall2(s, r);
1443
      else
1444
        return v? mkvecsmall2(0, s): mkvecsmall(s);
1445
    }
1446
  }
1447
  q = p >> 1;
1448
  split_init(&S, lg(f)-1);
1449
  settyp(S.done, t_VECSMALL);
1450
  if (v) split_add_done(&S, (GEN)0);
1451
  if (! split_Flx_cut_out_roots(&S, g, p))
1452
    return all_roots_mod_p(p, lg(S.done) == 1);
1453
  pol = polx_Flx(f[1]);
1454
  for (pol[2]=1; ; pol[2]++)
1455
  {
1456
    long j, l = lg(S.todo);
1457
    if (l == 1) { vecsmall_sort(S.done); return S.done; }
1458
    if (pol[2] == 100 && !uisprime(p)) pari_err_PRIME("polrootsmod",utoipos(p));
1459
    for (j = 1; j < l; j++)
1460
    {
1461
      GEN b, c = gel(S.todo,j);
1462
      ulong r, s;
1463
      switch(degpol(c))
1464
      {
1465
        case 1:
1466
          split_moveto_done(&S, j, (GEN)(p - c[2]));
1467
          j--; l--; break;
1468
        case 2:
1469
          r = Flx_quad_root(c, p, 0);
1470
          if (r == p) pari_err_PRIME("polrootsmod",utoipos(p));
1471
          s = Flx_otherroot(c,r, p);
1472
          split_done(&S, j, (GEN)r, (GEN)s);
1473
          j--; l--; break;
1474
        default:
1475
          b = Flxq_powu(pol,q, c,p); /* pol^(p-1)/2 */
1476
          if (degpol(b) <= 0) continue;
1477
          b = Flx_gcd(c,Flx_Fl_add(b,p-1,p), p);
1478
          if (!degpol(b)) continue;
1479
          b = Flx_normalize(b, p);
1480
          c = Flx_div(c,b, p);
1481
          split_todo(&S, j, b, c);
1482
      }
1483
    }
1484
  }
1485
}
1486

1487
GEN
1488
Flx_roots(GEN f, ulong p)
1489
{
1490
  pari_sp av = avma;
1491
  switch(lg(f))
1492
  {
1493
    case 2: pari_err_ROOTS0("Flx_roots");
1494
    case 3: set_avma(av); return cgetg(1, t_VECSMALL);
1495
  }
1496
  if (p == 2) return Flx_root_mod_2(f);
1497
  return gerepileuptoleaf(av, Flx_roots_i(Flx_normalize(f, p), p));
1498
}
1499

1500
/* assume x reduced mod p, monic. */
1501
static int
1502
Flx_quad_factortype(GEN x, ulong p)
1503
{
1504
  ulong b = x[3], c = x[2];
1505
  return krouu(Fl_disc_bc(b, c, p), p);
1506
}
1507
static GEN
1508
Flx_is_irred_2(GEN f, ulong p, long d)
1509
{
1510
  if (!d) return NULL;
1511
  if (d == 1) return gen_1;
1512
  return Flx_quad_factortype(f, p) == -1? gen_1: NULL;
1513
}
1514
static GEN
1515
Flx_degfact_2(GEN f, ulong p, long d)
1516
{
1517
  if (!d) return trivial_fact();
1518
  if (d == 1) return mkvec2(mkvecsmall(1), mkvecsmall(1));
1519
  switch(Flx_quad_factortype(f, p)) {
1520
    case 1: return mkvec2(mkvecsmall2(1,1), mkvecsmall2(1,1));
1521
    case -1:return mkvec2(mkvecsmall(2), mkvecsmall(1));
1522
    default: return mkvec2(mkvecsmall(1), mkvecsmall(2));
1523
  }
1524
}
1525
/* p > 2 */
1526
static GEN
1527
Flx_factor_2(GEN f, ulong p, long d)
1528
{
1529
  ulong r, s;
1530
  GEN R,S;
1531
  long v = f[1];
1532
  if (!d) return mkvec2(cgetg(1,t_COL), cgetg(1,t_VECSMALL));
1533
  if (labs(d) == 1) return mkvec2(mkcol(f), mkvecsmall(1));
1534
  r = Flx_quad_root(f, p, 1);
1535
  if (r==p) return mkvec2(mkcol(f), mkvecsmall(1));
1536
  s = Flx_otherroot(f, r, p);
1537
  r = Fl_neg(r, p);
1538
  s = Fl_neg(s, p);
1539
  if (s < r) lswap(s,r);
1540
  R = mkvecsmall3(v,r,1);
1541
  if (s == r) return mkvec2(mkcol(R), mkvecsmall(2));
1542
  S = mkvecsmall3(v,s,1);
1543
  return mkvec2(mkcol2(R,S), mkvecsmall2(1,1));
1544
}
1545
static GEN
1546
Flx_factor_deg2(GEN f, ulong p, long d, long flag)
1547
{
1548
  switch(flag) {
1549
    case 2: return Flx_is_irred_2(f, p, d);
1550
    case 1: return Flx_degfact_2(f, p, d);
1551
    default: return Flx_factor_2(f, p, d);
1552
  }
1553
}
1554

1555
static GEN
1556
F2x_Berlekamp_ker(GEN u)
1557
{
1558
  pari_sp ltop=avma;
1559
  long j,N = F2x_degree(u);
1560
  GEN Q;
1561
  pari_timer T;
1562
  timer_start(&T);
1563
  Q = F2x_matFrobenius(u);
1564
  for (j=1; j<=N; j++)
1565
    F2m_flip(Q,j,j);
1566
  if(DEBUGLEVEL>=9) timer_printf(&T,"Berlekamp matrix");
1567
  Q = F2m_ker_sp(Q,0);
1568
  if(DEBUGLEVEL>=9) timer_printf(&T,"kernel");
1569
  return gerepileupto(ltop,Q);
1570
}
1571
#define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}
1572
static long
1573
F2x_split_Berlekamp(GEN *t)
1574
{
1575
  GEN u = *t, a, b, vker;
1576
  long lb, d, i, ir, L, la, sv = u[1], du = F2x_degree(u);
1577

1578
  if (du == 1) return 1;
1579
  if (du == 2)
1580
  {
1581
    if (F2x_quad_factortype(u) == 1) /* 0 is a root: shouldn't occur */
1582
    {
1583
      t[0] = mkvecsmall2(sv, 2);
1584
      t[1] = mkvecsmall2(sv, 3);
1585
      return 2;
1586
    }
1587
    return 1;
1588
  }
1589

1590
  vker = F2x_Berlekamp_ker(u);
1591
  lb = lgcols(vker);
1592
  d = lg(vker)-1;
1593
  ir = 0;
1594
  /* t[i] irreducible for i < ir, still to be treated for i < L */
1595
  for (L=1; L<d; )
1596
  {
1597
    GEN pol;
1598
    if (d == 2)
1599
      pol = F2v_to_F2x(gel(vker,2), sv);
1600
    else
1601
    {
1602
      GEN v = zero_zv(lb);
1603
      v[1] = du;
1604
      v[2] = random_Fl(2); /*Assume vker[1]=1*/
1605
      for (i=2; i<=d; i++)
1606
        if (random_Fl(2)) F2v_add_inplace(v, gel(vker,i));
1607
      pol = F2v_to_F2x(v, sv);
1608
    }
1609
    for (i=ir; i<L && L<d; i++)
1610
    {
1611
      a = t[i]; la = F2x_degree(a);
1612
      if (la == 1) { set_irred(i); }
1613
      else if (la == 2)
1614
      {
1615
        if (F2x_quad_factortype(a) == 1) /* 0 is a root: shouldn't occur */
1616
        {
1617
          t[i] = mkvecsmall2(sv, 2);
1618
          t[L] = mkvecsmall2(sv, 3); L++;
1619
        }
1620
        set_irred(i);
1621
      }
1622
      else
1623
      {
1624
        pari_sp av = avma;
1625
        long lb;
1626
        b = F2x_rem(pol, a);
1627
        if (F2x_degree(b) <= 0) { set_avma(av); continue; }
1628
        b = F2x_gcd(a,b); lb = F2x_degree(b);
1629
        if (lb && lb < la)
1630
        {
1631
          t[L] = F2x_div(a,b);
1632
          t[i]= b; L++;
1633
        }
1634
        else set_avma(av);
1635
      }
1636
    }
1637
  }
1638
  return d;
1639
}
1640
/* assume deg f > 2 */
1641
static GEN
1642
F2x_Berlekamp_i(GEN f, long flag)
1643
{
1644
  long lfact, val, d = F2x_degree(f), j, k, lV;
1645
  GEN y, E, t, V;
1646

1647
  val = F2x_valrem(f, &f);
1648
  if (flag == 2 && val) return NULL;
1649
  V = F2x_factor_squarefree(f); lV = lg(V);
1650
  if (flag == 2 && lV > 2) return NULL;
1651

1652
  /* to hold factors and exponents */
1653
  t = cgetg(d+1, flag? t_VECSMALL: t_VEC);
1654
  E = cgetg(d+1,t_VECSMALL);
1655
  lfact = 1;
1656
  if (val) {
1657
    if (flag == 1) t[1] = 1; else gel(t,1) = polx_F2x(f[1]);
1658
    E[1] = val; lfact++;
1659
  }
1660

1661
  for (k=1; k<lV; k++)
1662
  {
1663
    if (F2x_degree(gel(V, k))==0) continue;
1664
    gel(t,lfact) = gel(V, k);
1665
    d = F2x_split_Berlekamp(&gel(t,lfact));
1666
    if (flag == 2 && d != 1) return NULL;
1667
    if (flag == 1)
1668
      for (j=0; j<d; j++) t[lfact+j] = F2x_degree(gel(t,lfact+j));
1669
    for (j=0; j<d; j++) E[lfact+j] = k;
1670
    lfact += d;
1671
  }
1672
  if (flag == 2) return gen_1; /* irreducible */
1673
  setlg(t, lfact);
1674
  setlg(E, lfact); y = mkvec2(t,E);
1675
  return flag ? sort_factor(y, (void*)&cmpGuGu, cmp_nodata)
1676
              : sort_factor_pol(y, cmpGuGu);
1677
}
1678

1679
/* Adapted from Shoup NTL */
1680
GEN
1681
F2x_factor_squarefree(GEN f)
1682
{
1683
  GEN r, t, v, tv;
1684
  long i, q, n = F2x_degree(f);
1685
  GEN u = const_vec(n+1, pol1_F2x(f[1]));
1686
  for(q = 1;;q *= 2)
1687
  {
1688
    r = F2x_gcd(f, F2x_deriv(f));
1689
    if (F2x_degree(r) == 0)
1690
    {
1691
      gel(u, q) = f;
1692
      break;
1693
    }
1694
    t = F2x_div(f, r);
1695
    if (F2x_degree(t) > 0)
1696
    {
1697
      long j;
1698
      for(j = 1;;j++)
1699
      {
1700
        v = F2x_gcd(r, t);
1701
        tv = F2x_div(t, v);
1702
        if (F2x_degree(tv) > 0)
1703
          gel(u, j*q) = tv;
1704
        if (F2x_degree(v) <= 0) break;
1705
        r = F2x_div(r, v);
1706
        t = v;
1707
      }
1708
      if (F2x_degree(r) == 0) break;
1709
    }
1710
    f = F2x_sqrt(r);
1711
  }
1712
  for (i = n; i; i--)
1713
    if (F2x_degree(gel(u,i))) break;
1714
  setlg(u,i+1); return u;
1715
}
1716

1717
static GEN
1718
F2x_ddf_simple(GEN T, GEN XP)
1719
{
1720
  pari_sp av = avma, av2;
1721
  GEN f, z, Tr, X;
1722
  long j, n = F2x_degree(T), v = T[1], B = n/2;
1723
  if (n == 0) return cgetg(1, t_VEC);
1724
  if (n == 1) return mkvec(T);
1725
  z = XP; Tr = T; X = polx_F2x(v);
1726
  f = const_vec(n, pol1_F2x(v));
1727
  av2 = avma;
1728
  for (j = 1; j <= B; j++)
1729
  {
1730
    GEN u = F2x_gcd(Tr, F2x_add(z, X));
1731
    if (F2x_degree(u))
1732
    {
1733
      gel(f, j) = u;
1734
      Tr = F2x_div(Tr, u);
1735
      av2 = avma;
1736
    } else z = gerepileuptoleaf(av2, z);
1737
    if (!F2x_degree(Tr)) break;
1738
    z = F2xq_sqr(z, Tr);
1739
  }
1740
  if (F2x_degree(Tr)) gel(f, F2x_degree(Tr)) = Tr;
1741
  return gerepilecopy(av, f);
1742
}
1743

1744
GEN
1745
F2x_ddf(GEN T)
1746
{
1747
  GEN XP;
1748
  T = F2x_get_red(T);
1749
  XP = F2x_Frobenius(T);
1750
  return F2x_ddf_to_ddf2(F2x_ddf_simple(T, XP));
1751
}
1752

1753
static GEN
1754
F2xq_frobtrace(GEN a, long d, GEN T)
1755
{
1756
  pari_sp av = avma;
1757
  long i;
1758
  GEN x = a;
1759
  for(i=1; i<d; i++)
1760
  {
1761
    x = F2x_add(a, F2xq_sqr(x,T));
1762
    if (gc_needed(av, 2))
1763
      x = gerepileuptoleaf(av, x);
1764
  }
1765
  return x;
1766
}
1767

1768
static void
1769
F2x_edf_simple(GEN Tp, GEN XP, long d, GEN V, long idx)
1770
{
1771
  long n = F2x_degree(Tp), r = n/d;
1772
  GEN T, f, ff;
1773
  if (r==1) { gel(V, idx) = Tp; return; }
1774
  T = Tp;
1775
  XP = F2x_rem(XP, T);
1776
  while (1)
1777
  {
1778
    pari_sp btop = avma;
1779
    long df;
1780
    GEN g = random_F2x(n, Tp[1]);
1781
    GEN t = F2xq_frobtrace(g, d, T);
1782
    if (lgpol(t) == 0) continue;
1783
    f = F2x_gcd(t, Tp); df = F2x_degree(f);
1784
    if (df > 0 && df < n) break;
1785
    set_avma(btop);
1786
  }
1787
  ff = F2x_div(Tp, f);
1788
  F2x_edf_simple(f, XP, d, V, idx);
1789
  F2x_edf_simple(ff, XP, d, V, idx+F2x_degree(f)/d);
1790
}
1791

1792
static GEN
1793
F2x_factor_Shoup(GEN T)
1794
{
1795
  long i, n, s = 0;
1796
  GEN XP, D, V;
1797
  pari_timer ti;
1798
  n = F2x_degree(T);
1799
  if (DEBUGLEVEL>=6) timer_start(&ti);
1800
  XP = F2x_Frobenius(T);
1801
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1802
  D = F2x_ddf_simple(T, XP);
1803
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1804
  for (i = 1; i <= n; i++)
1805
    s += F2x_degree(gel(D,i))/i;
1806
  V = cgetg(s+1, t_COL);
1807
  for (i = 1, s = 1; i <= n; i++)
1808
  {
1809
    GEN Di = gel(D,i);
1810
    long ni = F2x_degree(Di), ri = ni/i;
1811
    if (ni == 0) continue;
1812
    if (ni == i) { gel(V, s++) = Di; continue; }
1813
    F2x_edf_simple(Di, XP, i, V, s);
1814
    if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_edf(%ld)",i);
1815
    s += ri;
1816
  }
1817
  return V;
1818
}
1819

1820
static GEN
1821
F2x_factor_Cantor(GEN T)
1822
{
1823
  GEN E, F, V = F2x_factor_squarefree(T);
1824
  long i, j, l = lg(V);
1825
  E = cgetg(l, t_VEC);
1826
  F = cgetg(l, t_VEC);
1827
  for (i=1, j=1; i < l; i++)
1828
    if (F2x_degree(gel(V,i)))
1829
    {
1830
      GEN Fj = F2x_factor_Shoup(gel(V,i));
1831
      gel(F, j) = Fj;
1832
      gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1833
      j++;
1834
    }
1835
  return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
1836
}
1837

1838
#if 0
1839
static GEN
1840
F2x_simplefact_Shoup(GEN T)
1841
{
1842
  long i, n, s = 0, j = 1, k;
1843
  GEN XP, D, V;
1844
  pari_timer ti;
1845
  n = F2x_degree(T);
1846
  if (DEBUGLEVEL>=6) timer_start(&ti);
1847
  XP = F2x_Frobenius(T);
1848
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1849
  D = F2x_ddf_simple(T, XP);
1850
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1851
  for (i = 1; i <= n; i++)
1852
    s += F2x_degree(gel(D,i))/i;
1853
  V = cgetg(s+1, t_VECSMALL);
1854
  for (i = 1; i <= n; i++)
1855
  {
1856
    long ni = F2x_degree(gel(D,i)), ri = ni/i;
1857
    if (ni == 0) continue;
1858
    for (k = 1; k <= ri; k++)
1859
      V[j++] = i;
1860
  }
1861
  return V;
1862
}
1863
static GEN
1864
F2x_simplefact_Cantor(GEN T)
1865
{
1866
  GEN E, F, V = F2x_factor_squarefree(T);
1867
  long i, j, l = lg(V);
1868
  F = cgetg(l, t_VEC);
1869
  E = cgetg(l, t_VEC);
1870
  for (i=1, j=1; i < l; i++)
1871
    if (F2x_degree(gel(V,i)))
1872
    {
1873
      GEN Fj = F2x_simplefact_Shoup(gel(V,i));
1874
      gel(F, j) = Fj;
1875
      gel(E, j) = const_vecsmall(lg(Fj)-1, i);
1876
      j++;
1877
    }
1878
  return sort_factor(FE_concat(F,E,j), (void*)&cmpGuGu, cmp_nodata);
1879
}
1880
static int
1881
F2x_isirred_Cantor(GEN T)
1882
{
1883
  pari_sp av = avma;
1884
  pari_timer ti;
1885
  long n;
1886
  GEN dT = F2x_deriv(T);
1887
  GEN XP, D;
1888
  if (F2x_degree(F2x_gcd(T, dT)) != 0) return gc_bool(av,0);
1889
  n = F2x_degree(T);
1890
  if (DEBUGLEVEL>=6) timer_start(&ti);
1891
  XP = F2x_Frobenius(T);
1892
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_Frobenius");
1893
  D = F2x_ddf_simple(T, XP);
1894
  if (DEBUGLEVEL>=6) timer_printf(&ti,"F2x_ddf_simple");
1895
  return gc_bool(av, F2x_degree(gel(D,n)) == n);
1896
}
1897
#endif
1898

1899
/* driver for Cantor factorization, assume deg f > 2; not competitive for
1900
 * flag != 0, or as deg f increases */
1901
static GEN
1902
F2x_Cantor_i(GEN f, long flag)
1903
{
1904
  switch(flag)
1905
  {
1906
    default: return F2x_factor_Cantor(f);
1907
#if 0
1908
    case 1: return F2x_simplefact_Cantor(f);
1909
    case 2: return F2x_isirred_Cantor(f)? gen_1: NULL;
1910
#endif
1911
  }
1912
}
1913
static GEN
1914
F2x_factor_i(GEN f, long flag)
1915
{
1916
  long d = F2x_degree(f);
1917
  if (d <= 2) return F2x_factor_deg2(f,d,flag);
1918
  return (flag == 0 && d <= 20)? F2x_Cantor_i(f, flag)
1919
                               : F2x_Berlekamp_i(f, flag);
1920
}
1921

1922
GEN
1923
F2x_degfact(GEN f)
1924
{
1925
  pari_sp av = avma;
1926
  GEN z = F2x_factor_i(f, 1);
1927
  return gerepilecopy(av, z);
1928
}
1929

1930
int
1931
F2x_is_irred(GEN f) { return !!F2x_factor_i(f, 2); }
1932

1933
/* Adapted from Shoup NTL */
1934
GEN
1935
Flx_factor_squarefree(GEN f, ulong p)
1936
{
1937
  long i, q, n = degpol(f);
1938
  GEN u = const_vec(n+1, pol1_Flx(f[1]));
1939
  for(q = 1;;q *= p)
1940
  {
1941
    GEN t, v, tv, r = Flx_gcd(f, Flx_deriv(f, p), p);
1942
    if (degpol(r) == 0) { gel(u, q) = f; break; }
1943
    t = Flx_div(f, r, p);
1944
    if (degpol(t) > 0)
1945
    {
1946
      long j;
1947
      for(j = 1;;j++)
1948
      {
1949
        v = Flx_gcd(r, t, p);
1950
        tv = Flx_div(t, v, p);
1951
        if (degpol(tv) > 0)
1952
          gel(u, j*q) = Flx_normalize(tv, p);
1953
        if (degpol(v) <= 0) break;
1954
        r = Flx_div(r, v, p);
1955
        t = v;
1956
      }
1957
      if (degpol(r) == 0) break;
1958
    }
1959
    f = Flx_normalize(Flx_deflate(r, p), p);
1960
  }
1961
  for (i = n; i; i--)
1962
    if (degpol(gel(u,i))) break;
1963
  setlg(u,i+1); return u;
1964
}
1965

1966
long
1967
Flx_ispower(GEN f, ulong k, ulong p, GEN *pt_r)
1968
{
1969
  pari_sp av = avma;
1970
  ulong lc;
1971
  GEN F;
1972
  long i, n = degpol(f), v = f[1], l;
1973
  if (n % k) return 0;
1974
  lc = Fl_sqrtn(Flx_lead(f), k, p, NULL);
1975
  if (lc == ULONG_MAX) { av = avma; return 0; }
1976
  F = Flx_factor_squarefree(f, p); l = lg(F)-1;
1977
  for (i = 1; i <= l; i++)
1978
    if (i%k && degpol(gel(F,i))) return gc_long(av,0);
1979
  if (pt_r)
1980
  {
1981
    GEN r = Fl_to_Flx(lc, v), s = pol1_Flx(v);
1982
    for(i = l; i >= 1; i--)
1983
    {
1984
      if (i%k) continue;
1985
      s = Flx_mul(s, gel(F,i), p);
1986
      r = Flx_mul(r, s, p);
1987
    }
1988
    *pt_r = gerepileuptoleaf(av, r);
1989
  } else set_avma(av);
1990
  return 1;
1991
}
1992

1993
/* See <http://www.shoup.net/papers/factorimpl.pdf> */
1994
static GEN
1995
Flx_ddf_Shoup(GEN T, GEN XP, ulong p)
1996
{
1997
  pari_sp av = avma;
1998
  GEN b, g, h, F, f, Tr, xq;
1999
  long i, j, n, v, bo, ro;
2000
  long B, l, m;
2001
  pari_timer ti;
2002
  n = get_Flx_degree(T); v = get_Flx_var(T);
2003
  if (n == 0) return cgetg(1, t_VEC);
2004
  if (n == 1) return mkvec(get_Flx_mod(T));
2005
  B = n/2;
2006
  l = usqrt(B);
2007
  m = (B+l-1)/l;
2008
  T = Flx_get_red(T, p);
2009
  b = cgetg(l+2, t_VEC);
2010
  gel(b, 1) = polx_Flx(v);
2011
  gel(b, 2) = XP;
2012
  bo = brent_kung_optpow(n, l-1, 1);
2013
  ro = l<=1 ? 0:(bo-1)/(l-1) + ((n-1)/bo);
2014
  if (DEBUGLEVEL>=7) timer_start(&ti);
2015
  if (expu(p) <= ro)
2016
    for (i = 3; i <= l+1; i++)
2017
      gel(b, i) = Flxq_powu(gel(b, i-1), p, T, p);
2018
  else
2019
  {
2020
    xq = Flxq_powers(gel(b, 2), bo,  T, p);
2021
    if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq baby");
2022
    for (i = 3; i <= l+1; i++)
2023
      gel(b, i) = Flx_FlxqV_eval(gel(b, i-1), xq, T, p);
2024
  }
2025
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: baby");
2026
  xq = Flxq_powers(gel(b, l+1), brent_kung_optpow(n, m-1, 1),  T, p);
2027
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: xq giant");
2028
  g = cgetg(m+1, t_VEC);
2029
  gel(g, 1) = gel(xq, 2);
2030
  for(i = 2; i <= m; i++)
2031
    gel(g, i) = Flx_FlxqV_eval(gel(g, i-1), xq, T, p);
2032
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: giant");
2033
  h = cgetg(m+1, t_VEC);
2034
  for (j = 1; j <= m; j++)
2035
  {
2036
    pari_sp av = avma;
2037
    GEN gj = gel(g, j);
2038
    GEN e = Flx_sub(gj, gel(b, 1), p);
2039
    for (i = 2; i <= l; i++)
2040
      e = Flxq_mul(e, Flx_sub(gj, gel(b, i), p), T, p);
2041
    gel(h, j) = gerepileupto(av, e);
2042
  }
2043
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: diff");
2044
  Tr = get_Flx_mod(T);
2045
  F = cgetg(m+1, t_VEC);
2046
  for (j = 1; j <= m; j++)
2047
  {
2048
    GEN u = Flx_gcd(Tr, gel(h, j), p);
2049
    if (degpol(u))
2050
    {
2051
      u = Flx_normalize(u, p);
2052
      Tr = Flx_div(Tr, u, p);
2053
    }
2054
    gel(F, j) = u;
2055
  }
2056
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: F");
2057
  f = const_vec(n, pol1_Flx(v));
2058
  for (j = 1; j <= m; j++)
2059
  {
2060
    GEN e = gel(F, j);
2061
    for (i=l-1; i >= 0; i--)
2062
    {
2063
      GEN u = Flx_gcd(e, Flx_sub(gel(g, j), gel(b, i+1), p), p);
2064
      if (degpol(u))
2065
      {
2066
        gel(f, l*j-i) = u;
2067
        e = Flx_div(e, u, p);
2068
      }
2069
      if (!degpol(e)) break;
2070
    }
2071
  }
2072
  if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_ddf_Shoup: f");
2073
  if (degpol(Tr)) gel(f, degpol(Tr)) = Tr;
2074
  return gerepilecopy(av, f);
2075
}
2076

2077
static void
2078
Flx_edf_simple(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2079
{
2080
  long n = degpol(Tp), r = n/d;
2081
  GEN T, f, ff;
2082
  ulong p2;
2083
  if (r==1) { gel(V, idx) = Tp; return; }
2084
  p2 = p>>1;
2085
  T = Flx_get_red(Tp, p);
2086
  XP = Flx_rem(XP, T, p);
2087
  while (1)
2088
  {
2089
    pari_sp btop = avma;
2090
    long i;
2091
    GEN g = random_Flx(n, Tp[1], p);
2092
    GEN t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2093
    if (lgpol(t) == 0) continue;
2094
    for(i=1; i<=10; i++)
2095
    {
2096
      pari_sp btop2 = avma;
2097
      GEN R = Flxq_powu(Flx_Fl_add(t, random_Fl(p), p), p2, T, p);
2098
      f = Flx_gcd(Flx_Fl_add(R, p-1, p), Tp, p);
2099
      if (degpol(f) > 0 && degpol(f) < n) break;
2100
      set_avma(btop2);
2101
    }
2102
    if (degpol(f) > 0 && degpol(f) < n) break;
2103
    set_avma(btop);
2104
  }
2105
  f = Flx_normalize(f, p);
2106
  ff = Flx_div(Tp, f ,p);
2107
  Flx_edf_simple(f, XP, d, p, V, idx);
2108
  Flx_edf_simple(ff, XP, d, p, V, idx+degpol(f)/d);
2109
}
2110
static void
2111
Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx);
2112

2113
static void
2114
Flx_edf_rec(GEN T, GEN XP, GEN hp, GEN t, long d, ulong p, GEN V, long idx)
2115
{
2116
  pari_sp av;
2117
  GEN Tp = get_Flx_mod(T);
2118
  long n = degpol(hp), vT = Tp[1];
2119
  GEN u1, u2, f1, f2;
2120
  ulong p2 = p>>1;
2121
  GEN R, h;
2122
  h = Flx_get_red(hp, p);
2123
  t = Flx_rem(t, T, p);
2124
  av = avma;
2125
  do
2126
  {
2127
    set_avma(av);
2128
    R = Flxq_powu(mkvecsmall3(vT, random_Fl(p), 1), p2, h, p);
2129
    u1 = Flx_gcd(Flx_Fl_add(R, p-1, p), hp, p);
2130
  } while (degpol(u1)==0 || degpol(u1)==n);
2131
  f1 = Flx_gcd(Flx_Flxq_eval(u1, t, T, p), Tp, p);
2132
  f1 = Flx_normalize(f1, p);
2133
  u2 = Flx_div(hp, u1, p);
2134
  f2 = Flx_div(Tp, f1, p);
2135
  if (degpol(u1)==1)
2136
  {
2137
    if (degpol(f1)==d)
2138
      gel(V, idx) = f1;
2139
    else
2140
      Flx_edf(f1, XP, d, p, V, idx);
2141
  }
2142
  else
2143
    Flx_edf_rec(Flx_get_red(f1, p), XP, u1, t, d, p, V, idx);
2144
  idx += degpol(f1)/d;
2145
  if (degpol(u2)==1)
2146
  {
2147
    if (degpol(f2)==d)
2148
      gel(V, idx) = f2;
2149
    else
2150
      Flx_edf(f2, XP, d, p, V, idx);
2151
  }
2152
  else
2153
    Flx_edf_rec(Flx_get_red(f2, p), XP, u2, t, d, p, V, idx);
2154
}
2155

2156
static void
2157
Flx_edf(GEN Tp, GEN XP, long d, ulong p, GEN V, long idx)
2158
{
2159
  long n = degpol(Tp), r = n/d, vT = Tp[1];
2160
  GEN T, h, t;
2161
  pari_timer ti;
2162
  if (r==1) { gel(V, idx) = Tp; return; }
2163
  T = Flx_get_red(Tp, p);
2164
  XP = Flx_rem(XP, T, p);
2165
  if (DEBUGLEVEL>=7) timer_start(&ti);
2166
  do
2167
  {
2168
    GEN g = random_Flx(n, vT, p);
2169
    t = gel(Flxq_auttrace(mkvec2(XP, g), d, T, p), 2);
2170
    if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_auttrace");
2171
    h = Flxq_minpoly(t, T, p);
2172
    if (DEBUGLEVEL>=7) timer_printf(&ti,"Flx_edf: Flxq_minpoly");
2173
  } while (degpol(h) <= 1);
2174
  Flx_edf_rec(T, XP, h, t, d, p, V, idx);
2175
}
2176

2177
static GEN
2178
Flx_factor_Shoup(GEN T, ulong p)
2179
{
2180
  long i, n, s = 0;
2181
  GEN XP, D, V;
2182
  long e = expu(p);
2183
  pari_timer ti;
2184
  n = get_Flx_degree(T);
2185
  T = Flx_get_red(T, p);
2186
  if (DEBUGLEVEL>=6) timer_start(&ti);
2187
  XP = Flx_Frobenius(T, p);
2188
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2189
  D = Flx_ddf_Shoup(T, XP, p);
2190
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2191
  s = ddf_to_nbfact(D);
2192
  V = cgetg(s+1, t_COL);
2193
  for (i = 1, s = 1; i <= n; i++)
2194
  {
2195
    GEN Di = gel(D,i);
2196
    long ni = degpol(Di), ri = ni/i;
2197
    if (ni == 0) continue;
2198
    Di = Flx_normalize(Di, p);
2199
    if (ni == i) { gel(V, s++) = Di; continue; }
2200
    if (ri <= e*expu(e))
2201
      Flx_edf(Di, XP, i, p, V, s);
2202
    else
2203
      Flx_edf_simple(Di, XP, i, p, V, s);
2204
    if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_edf(%ld)",i);
2205
    s += ri;
2206
  }
2207
  return V;
2208
}
2209

2210
static GEN
2211
Flx_factor_Cantor(GEN T, ulong p)
2212
{
2213
  GEN E, F, V = Flx_factor_squarefree(get_Flx_mod(T), p);
2214
  long i, j, l = lg(V);
2215
  F = cgetg(l, t_VEC);
2216
  E = cgetg(l, t_VEC);
2217
  for (i=1, j=1; i < l; i++)
2218
    if (degpol(gel(V,i)))
2219
    {
2220
      GEN Fj = Flx_factor_Shoup(gel(V,i), p);
2221
      gel(F, j) = Fj;
2222
      gel(E, j) = const_vecsmall(lg(Fj)-1, i);
2223
      j++;
2224
    }
2225
  return sort_factor_pol(FE_concat(F,E,j), cmpGuGu);
2226
}
2227

2228
GEN
2229
Flx_ddf(GEN T, ulong p)
2230
{
2231
  GEN XP;
2232
  T = Flx_get_red(T, p);
2233
  XP = Flx_Frobenius(T, p);
2234
  return ddf_to_ddf2(Flx_ddf_Shoup(T, XP, p));
2235
}
2236

2237
static GEN
2238
Flx_simplefact_Cantor(GEN T, ulong p)
2239
{
2240
  GEN V;
2241
  long i, l;
2242
  T = Flx_get_red(T, p);
2243
  V = Flx_factor_squarefree(get_Flx_mod(T), p); l = lg(V);
2244
  for (i=1; i < l; i++)
2245
    gel(V,i) = Flx_ddf_Shoup(gel(V,i), Flx_Frobenius(gel(V,i), p), p);
2246
  return vddf_to_simplefact(V, get_Flx_degree(T));
2247
}
2248

2249
static int
2250
Flx_isirred_Cantor(GEN Tp, ulong p)
2251
{
2252
  pari_sp av = avma;
2253
  pari_timer ti;
2254
  long n;
2255
  GEN T = get_Flx_mod(Tp), dT = Flx_deriv(T, p), XP, D;
2256
  if (degpol(Flx_gcd(T, dT, p)) != 0) return gc_bool(av,0);
2257
  n = get_Flx_degree(T);
2258
  T = Flx_get_red(Tp, p);
2259
  if (DEBUGLEVEL>=6) timer_start(&ti);
2260
  XP = Flx_Frobenius(T, p);
2261
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2262
  D = Flx_ddf_Shoup(T, XP, p);
2263
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2264
  return gc_bool(av, degpol(gel(D,n)) == n);
2265
}
2266

2267
/* f monic */
2268
static GEN
2269
Flx_factor_i(GEN f, ulong pp, long flag)
2270
{
2271
  long d;
2272
  if (pp==2) { /*We need to handle 2 specially */
2273
    GEN F = F2x_factor_i(Flx_to_F2x(f),flag);
2274
    if (flag==0) F2xV_to_FlxV_inplace(gel(F,1));
2275
    return F;
2276
  }
2277
  d = degpol(f);
2278
  if (d <= 2) return Flx_factor_deg2(f,pp,d,flag);
2279
  switch(flag)
2280
  {
2281
    default: return Flx_factor_Cantor(f, pp);
2282
    case 1: return Flx_simplefact_Cantor(f, pp);
2283
    case 2: return Flx_isirred_Cantor(f, pp)? gen_1: NULL;
2284
  }
2285
}
2286

2287
GEN
2288
Flx_degfact(GEN f, ulong p)
2289
{
2290
  pari_sp av = avma;
2291
  GEN z = Flx_factor_i(Flx_normalize(f,p),p,1);
2292
  return gerepilecopy(av, z);
2293
}
2294

2295
/* T must be squarefree mod p*/
2296
GEN
2297
Flx_nbfact_by_degree(GEN T, long *nb, ulong p)
2298
{
2299
  GEN XP, D;
2300
  pari_timer ti;
2301
  long i, s, n = get_Flx_degree(T);
2302
  GEN V = const_vecsmall(n, 0);
2303
  pari_sp av = avma;
2304
  T = Flx_get_red(T, p);
2305
  if (DEBUGLEVEL>=6) timer_start(&ti);
2306
  XP = Flx_Frobenius(T, p);
2307
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_Frobenius");
2308
  D = Flx_ddf_Shoup(T, XP, p);
2309
  if (DEBUGLEVEL>=6) timer_printf(&ti,"Flx_ddf_Shoup");
2310
  for (i = 1, s = 0; i <= n; i++)
2311
  {
2312
    V[i] = degpol(gel(D,i))/i;
2313
    s += V[i];
2314
  }
2315
  *nb = s;
2316
  set_avma(av); return V;
2317
}
2318

2319
long
2320
Flx_nbfact_Frobenius(GEN T, GEN XP, ulong p)
2321
{
2322
  pari_sp av = avma;
2323
  long s = ddf_to_nbfact(Flx_ddf_Shoup(T, XP, p));
2324
  return gc_long(av,s);
2325
}
2326

2327
/* T must be squarefree mod p*/
2328
long
2329
Flx_nbfact(GEN T, ulong p)
2330
{
2331
  pari_sp av = avma;
2332
  GEN XP = Flx_Frobenius(T, p);
2333
  long n = Flx_nbfact_Frobenius(T, XP, p);
2334
  return gc_long(av,n);
2335
}
2336

2337
int
2338
Flx_is_irred(GEN f, ulong p)
2339
{
2340
  pari_sp av = avma;
2341
  f = Flx_normalize(f,p);
2342
  return gc_bool(av, !!Flx_factor_i(f,p,2));
2343
}
2344

2345
/* Use this function when you think f is reducible, and that there are lots of
2346
 * factors. If you believe f has few factors, use FpX_nbfact(f,p)==1 instead */
2347
int
2348
FpX_is_irred(GEN f, GEN p)
2349
{
2350
  pari_sp av = avma;
2351
  int z;
2352
  switch(ZX_factmod_init(&f,p))
2353
  {
2354
    case 0:  z = !!F2x_factor_i(f,2); break;
2355
    case 1:  z = !!Flx_factor_i(f,p[2],2); break;
2356
    default: z = !!FpX_factor_i(f,p,2); break;
2357
  }
2358
  return gc_bool(av,z);
2359
}
2360
GEN
2361
FpX_degfact(GEN f, GEN p) {
2362
  pari_sp av = avma;
2363
  GEN F;
2364
  switch(ZX_factmod_init(&f,p))
2365
  {
2366
    case 0:  F = F2x_factor_i(f,1); break;
2367
    case 1:  F = Flx_factor_i(f,p[2],1); break;
2368
    default: F = FpX_factor_i(f,p,1); break;
2369
  }
2370
  return gerepilecopy(av, F);
2371
}
2372

2373
#if 0
2374
/* set x <-- x + c*y mod p */
2375
/* x is not required to be normalized.*/
2376
static void
2377
Flx_addmul_inplace(GEN gx, GEN gy, ulong c, ulong p)
2378
{
2379
  long i, lx, ly;
2380
  ulong *x=(ulong *)gx;
2381
  ulong *y=(ulong *)gy;
2382
  if (!c) return;
2383
  lx = lg(gx);
2384
  ly = lg(gy);
2385
  if (lx<ly) pari_err_BUG("lx<ly in Flx_addmul_inplace");
2386
  if (SMALL_ULONG(p))
2387
    for (i=2; i<ly;  i++) x[i] = (x[i] + c*y[i]) % p;
2388
  else
2389
    for (i=2; i<ly;  i++) x[i] = Fl_add(x[i], Fl_mul(c,y[i],p),p);
2390
}
2391
#endif
2392

2393
GEN
2394
FpX_factor(GEN f, GEN p)
2395
{
2396
  pari_sp av = avma;
2397
  GEN F;
2398
  switch(ZX_factmod_init(&f, p))
2399
  {
2400
    case 0:  F = F2x_factor_i(f,0);
2401
             F2xV_to_ZXV_inplace(gel(F,1)); break;
2402
    case 1:  F = Flx_factor_i(f,p[2],0);
2403
             FlxV_to_ZXV_inplace(gel(F,1)); break;
2404
    default: F = FpX_factor_i(f,p,0); break;
2405
  }
2406
  return gerepilecopy(av, F);
2407
}
2408

2409
GEN
2410
Flx_factor(GEN f, ulong p)
2411
{
2412
  pari_sp av = avma;
2413
  return gerepilecopy(av, Flx_factor_i(Flx_normalize(f,p),p,0));
2414
}
2415
GEN
2416
F2x_factor(GEN f)
2417
{
2418
  pari_sp av = avma;
2419
  return gerepilecopy(av, F2x_factor_i(f,0));
2420
}
2421

2422
Product

Resources

Company
