Line data Source code
1 : /* Copyright (C) 2000-2004 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 : /*******************************************************************/
16 : /* */
17 : /* SUBFIELDS OF A NUMBER FIELD */
18 : /* J. Klueners and M. Pohst, J. Symb. Comp. (1996), vol. 11 */
19 : /* */
20 : /*******************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_nfsubfields
25 :
26 : typedef struct _poldata {
27 : GEN pol;
28 : GEN dis; /* |disc(pol)| */
29 : GEN roo; /* roots(pol) */
30 : GEN den; /* multiple of index(pol) */
31 : } poldata;
32 : typedef struct _primedata {
33 : GEN p; /* prime */
34 : GEN pol; /* pol mod p, squarefree */
35 : GEN ff; /* factorization of pol mod p */
36 : GEN Z; /* cycle structure of the above [ Frobenius orbits ] */
37 : long lcm; /* lcm of the above */
38 : GEN T; /* ffinit(p, lcm) */
39 :
40 : GEN fk; /* factorization of pol over F_(p^lcm) */
41 : GEN firstroot; /* *[i] = index of first root of fk[i] */
42 : GEN interp; /* *[i] = interpolation polynomial for fk[i]
43 : * [= 1 on the first root firstroot[i], 0 on the others] */
44 : GEN bezoutC; /* Bezout coefficients attached to the ff[i] */
45 : GEN Trk; /* used to compute traces (cf poltrace) */
46 : } primedata;
47 : typedef struct _blockdata {
48 : poldata *PD; /* data depending from pol */
49 : primedata *S;/* data depending from pol, p */
50 : GEN DATA; /* data depending from pol, p, degree, # translations [updated] */
51 : long N; /* deg(PD.pol) */
52 : long d; /* subfield degree */
53 : long size;/* block degree = N/d */
54 : long fl;
55 : } blockdata;
56 :
57 : static GEN print_block_system(blockdata *B, GEN Y, GEN BS);
58 : static GEN test_block(blockdata *B, GEN L, GEN D);
59 :
60 : /* COMBINATORIAL PART: generate potential block systems */
61 :
62 : #define BIL 32 /* for 64bit machines also */
63 : /* Computation of potential block systems of given size d attached to a
64 : * rational prime p: give a row vector of row vectors containing the
65 : * potential block systems of imprimitivity; a potential block system is a
66 : * vector of row vectors (enumeration of the roots). */
67 : static GEN
68 33234 : calc_block(blockdata *B, GEN Z, GEN Y, GEN SB)
69 : {
70 33234 : long r = lg(Z), lK, i, j, t, tp, T, u, nn, lnon, lY;
71 : GEN K, n, non, pn, pnon, e, Yp, Zp, Zpp;
72 33234 : pari_sp av0 = avma;
73 :
74 33234 : if (DEBUGLEVEL>3)
75 : {
76 0 : err_printf("lg(Z) = %ld, lg(Y) = %ld\n", r,lg(Y));
77 0 : if (DEBUGLEVEL > 5)
78 : {
79 0 : err_printf("Z = %Ps\n",Z);
80 0 : err_printf("Y = %Ps\n",Y);
81 : }
82 : }
83 33234 : lnon = minss(BIL, r);
84 33235 : e = new_chunk(BIL);
85 33236 : n = new_chunk(r);
86 33236 : non = new_chunk(lnon);
87 33236 : pnon = new_chunk(lnon);
88 33236 : pn = new_chunk(lnon);
89 :
90 33236 : Zp = cgetg(lnon,t_VEC);
91 33236 : Zpp = cgetg(lnon,t_VEC); nn = 0;
92 68530 : for (i=1; i<r; i++) { n[i] = lg(gel(Z,i))-1; nn += n[i]; }
93 33236 : lY = lg(Y); Yp = cgetg(lY+1,t_VEC);
94 35525 : for (j=1; j<lY; j++) gel(Yp,j) = gel(Y,j);
95 :
96 : {
97 33236 : pari_sp av = avma;
98 33236 : long k = nn / B->size;
99 67655 : for (j = 1; j < r; j++)
100 34769 : if (n[j] % k) break;
101 33236 : if (j == r)
102 : {
103 32886 : gel(Yp,lY) = Z;
104 32886 : SB = print_block_system(B, Yp, SB);
105 32885 : set_avma(av);
106 : }
107 : }
108 33235 : gel(Yp,lY) = Zp;
109 :
110 33235 : K = divisorsu(n[1]); lK = lg(K);
111 141375 : for (i=1; i<lK; i++)
112 : {
113 108140 : long ngcd = n[1], k = K[i], dk = B->size*k, lpn = 0;
114 115497 : for (j=2; j<r; j++)
115 7357 : if (n[j]%k == 0)
116 : {
117 7301 : if (++lpn >= BIL) pari_err_OVERFLOW("calc_block");
118 7301 : pn[lpn] = n[j]; pnon[lpn] = j;
119 7301 : ngcd = ugcd(ngcd, n[j]);
120 : }
121 108140 : if (dk % ngcd) continue;
122 66946 : T = 1L<<lpn;
123 66946 : if (lpn == r-2)
124 : {
125 66890 : T--; /* done already above --> print_block_system */
126 66890 : if (!T) continue;
127 : }
128 :
129 3766 : if (dk == n[1])
130 : { /* empty subset, t = 0. Split out for clarity */
131 1680 : Zp[1] = Z[1]; setlg(Zp, 2);
132 3563 : for (u=1,j=2; j<r; j++) Zpp[u++] = Z[j];
133 1680 : setlg(Zpp, u);
134 1680 : SB = calc_block(B, Zpp, Yp, SB);
135 : }
136 :
137 4690 : for (t = 1; t < T; t++)
138 : { /* loop through all nonempty subsets of [1..lpn] */
139 2898 : for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
140 : {
141 1974 : if (tp&1) { nn += pn[u]; e[u] = 1; } else e[u] = 0;
142 : }
143 924 : if (dk != nn) continue;
144 :
145 1505 : for (j=1; j<r; j++) non[j]=0;
146 371 : Zp[1] = Z[1];
147 1134 : for (u=2,j=1; j<=lpn; j++)
148 763 : if (e[j]) { Zp[u] = Z[pnon[j]]; non[pnon[j]] = 1; u++; }
149 371 : setlg(Zp, u);
150 1134 : for (u=1,j=2; j<r; j++)
151 763 : if (!non[j]) Zpp[u++] = Z[j];
152 371 : setlg(Zpp, u);
153 371 : SB = calc_block(B, Zpp, Yp, SB);
154 : }
155 : }
156 33235 : return gc_const(av0, SB);
157 : }
158 :
159 : /* product of permutations. Put the result in perm1. */
160 : static void
161 230048 : perm_mul_i(GEN perm1, GEN perm2)
162 : {
163 230048 : long i, N = lg(perm1);
164 230048 : pari_sp av = avma;
165 230048 : GEN perm = new_chunk(N);
166 10451504 : for (i=1; i<N; i++) perm[i] = perm1[perm2[i]];
167 10451504 : for (i=1; i<N; i++) perm1[i]= perm[i];
168 230048 : set_avma(av);
169 230048 : }
170 :
171 : /* cy is a cycle; compute cy^l as a permutation */
172 : static GEN
173 47789 : cycle_power_to_perm(GEN perm,GEN cy,long l)
174 : {
175 47789 : long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
176 :
177 47789 : lp = l % lcy;
178 1989575 : for (i=1; i<N; i++) perm[i] = i;
179 47789 : if (lp)
180 : {
181 42511 : pari_sp av = avma;
182 42511 : GEN p1 = new_chunk(N);
183 42511 : b = cy[1];
184 459158 : for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
185 42511 : perm[b] = cy[1];
186 1802003 : for (i=1; i<N; i++) p1[i] = perm[i];
187 :
188 224770 : for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
189 42511 : set_avma(av);
190 : }
191 47789 : return perm;
192 : }
193 :
194 : /* image du block system D par la permutation perm */
195 : static GEN
196 22491 : im_block_by_perm(GEN D,GEN perm)
197 : {
198 22491 : long i, lb = lg(D);
199 22491 : GEN Dn = cgetg(lb,t_VEC);
200 233786 : for (i=1; i<lb; i++) gel(Dn,i) = vecsmallpermute(perm, gel(D,i));
201 22491 : return Dn;
202 : }
203 :
204 : static void
205 35154 : append(GEN D, GEN a)
206 : {
207 35154 : long i,l = lg(D), m = lg(a);
208 35154 : GEN x = D + (l-1);
209 121198 : for (i=1; i<m; i++) gel(x,i) = gel(a,i);
210 35154 : setlg(D, l+m-1);
211 35154 : }
212 :
213 : static GEN
214 32886 : print_block_system(blockdata *B, GEN Y, GEN SB)
215 : {
216 32886 : long i, j, l, ll, lp, u, v, ns, r = lg(Y), N = B->N;
217 : long *k, *n, **e, *t;
218 32886 : GEN D, De, Z, cyperm, perm, VOID = cgetg(1, t_VECSMALL);
219 :
220 32886 : if (DEBUGLEVEL>5) err_printf("Y = %Ps\n",Y);
221 32886 : n = new_chunk(N+1);
222 32886 : D = vectrunc_init(N+1);
223 32886 : t = new_chunk(r+1);
224 32886 : k = new_chunk(r+1);
225 32886 : Z = cgetg(r+1, t_VEC);
226 68040 : for (ns=0,i=1; i<r; i++)
227 : {
228 35154 : GEN Yi = gel(Y,i);
229 35154 : long ki = 0, si = lg(Yi)-1;
230 :
231 72233 : for (j=1; j<=si; j++) { n[j] = lg(gel(Yi,j))-1; ki += n[j]; }
232 35154 : ki /= B->size;
233 35154 : De = cgetg(ki+1,t_VEC);
234 121198 : for (j=1; j<=ki; j++) gel(De,j) = VOID;
235 72233 : for (j=1; j<=si; j++)
236 : {
237 37079 : GEN cy = gel(Yi,j);
238 214963 : for (l=1,lp=0; l<=n[j]; l++)
239 : {
240 177884 : lp++; if (lp > ki) lp = 1;
241 177884 : gel(De,lp) = vecsmall_append(gel(De,lp), cy[l]);
242 : }
243 : }
244 35154 : append(D, De);
245 35154 : if (si>1 && ki>1)
246 : {
247 1883 : GEN p1 = cgetg(si,t_VEC);
248 3808 : for (j=2; j<=si; j++) p1[j-1] = Yi[j];
249 1883 : ns++;
250 1883 : t[ns] = si-1;
251 1883 : k[ns] = ki-1;
252 1883 : gel(Z,ns) = p1;
253 : }
254 : }
255 32886 : if (DEBUGLEVEL>2) err_printf("\nns = %ld\n",ns);
256 32886 : if (!ns) return test_block(B, SB, D);
257 :
258 1862 : setlg(Z, ns+1);
259 1862 : e = (long**)new_chunk(ns+1);
260 3745 : for (i=1; i<=ns; i++)
261 : {
262 1883 : e[i] = new_chunk(t[i]+1);
263 3808 : for (j=1; j<=t[i]; j++) e[i][j] = 0;
264 : }
265 1862 : cyperm= cgetg(N+1,t_VECSMALL);
266 1862 : perm = cgetg(N+1,t_VECSMALL); i = ns;
267 : do
268 : {
269 22491 : pari_sp av = avma;
270 762951 : for (u=1; u<=N; u++) perm[u] = u;
271 45738 : for (u=1; u<=ns; u++)
272 71036 : for (v=1; v<=t[u]; v++)
273 47789 : perm_mul_i(perm, cycle_power_to_perm(cyperm, gmael(Z,u,v), e[u][v]));
274 22491 : SB = test_block(B, SB, im_block_by_perm(D,perm));
275 22491 : set_avma(av);
276 :
277 : /* i = 1..ns, j = 1..t[i], e[i][j] loop through 0..k[i].
278 : * TODO: flatten to 1-dimensional loop */
279 22491 : if (++e[ns][t[ns]] > k[ns])
280 : {
281 3038 : j = t[ns]-1;
282 3157 : while (j>=1 && e[ns][j] == k[ns]) j--;
283 4186 : if (j >= 1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l] = 0; }
284 : else
285 : {
286 1967 : i = ns-1;
287 1988 : while (i>=1)
288 : {
289 126 : j = t[i];
290 147 : while (j>=1 && e[i][j] == k[i]) j--;
291 126 : if (j<1) i--;
292 : else
293 : {
294 105 : e[i][j]++;
295 105 : for (l=j+1; l<=t[i]; l++) e[i][l] = 0;
296 210 : for (ll=i+1; ll<=ns; ll++)
297 210 : for (l=1; l<=t[ll]; l++) e[ll][l] = 0;
298 105 : break;
299 : }
300 : }
301 : }
302 : }
303 : }
304 22491 : while (i > 0);
305 1862 : return SB;
306 : }
307 :
308 : /* ALGEBRAIC PART: test potential block systems */
309 :
310 : static GEN
311 40411 : polsimplify(GEN x)
312 : {
313 40411 : long i,lx = lg(x);
314 238742 : for (i=2; i<lx; i++)
315 198331 : if (typ(gel(x,i)) == t_POL) gel(x,i) = constant_coeff(gel(x,i));
316 40411 : return x;
317 : }
318 :
319 : /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
320 : static long
321 40411 : ok_coeffs(GEN g,GEN M)
322 : {
323 40411 : long i, lg = lg(g)-1; /* g is monic, and cst term is ok */
324 98123 : for (i=3; i<lg; i++)
325 64852 : if (abscmpii(gel(g,i), gel(M,i)) > 0) return 0;
326 33271 : return 1;
327 : }
328 :
329 : /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
330 : static GEN
331 174090 : trace(GEN x, GEN Trq, GEN p)
332 : {
333 : long i, l;
334 : GEN s;
335 174090 : if (typ(x) == t_INT) return Fp_mul(x, gel(Trq,1), p);
336 174090 : l = lg(x)-1; if (l == 1) return gen_0;
337 174090 : x++; s = mulii(gel(x,1), gel(Trq,1));
338 883795 : for (i=2; i<l; i++)
339 709705 : s = addii(s, mulii(gel(x,i), gel(Trq,i)));
340 174090 : return modii(s, p);
341 : }
342 :
343 : /* assume x in Fq[X], return Tr_{Fq[X]/Fp[X]}(x), varn(X) = 0 */
344 : static GEN
345 36581 : poltrace(GEN x, GEN Trq, GEN p)
346 : {
347 36581 : if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
348 210669 : pari_APPLY_pol(trace(gel(x,i),Trq,p));
349 : }
350 :
351 : /* Find h in Fp[X] such that h(a[i]) = listdelta[i] for all modular factors
352 : * ff[i], where a[i] is a fixed root of ff[i] in Fq = Z[Y]/(p,T) [namely the
353 : * first one in FpX_factorff_irred output]. Let f = ff[i], A the given root,
354 : * then h mod f is Tr_Fq/Fp ( h(A) f(X)/(X-A)f'(A) ), most of the expression
355 : * being precomputed. The complete h is recovered via chinese remaindering */
356 : static GEN
357 32506 : chinese_retrieve_pol(GEN DATA, primedata *S, GEN listdelta)
358 : {
359 32506 : GEN interp, bezoutC, h, p = S->p, pol = FpX_red(gel(DATA,1), p);
360 : long i, l;
361 32507 : interp = gel(DATA,9);
362 32507 : bezoutC= gel(DATA,6);
363 :
364 32507 : h = NULL; l = lg(interp);
365 69088 : for (i=1; i<l; i++)
366 : { /* h(firstroot[i]) = listdelta[i] */
367 36581 : GEN t = FqX_Fq_mul(gel(interp,i), gel(listdelta,i), S->T, p);
368 36581 : t = poltrace(t, gel(S->Trk,i), p);
369 36580 : t = FpX_mul(t, gel(bezoutC,i), p);
370 36581 : h = h? FpX_add(h,t,p): t;
371 : }
372 32507 : return FpX_rem(h, pol, p);
373 : }
374 :
375 : /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
376 : * (cf subfield) was a block system; then there
377 : * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
378 : * in Fp[X] (cf chinese_retrieve_pol). Try to lift it; den is a
379 : * multiplicative bound for denominator of lift. */
380 : static GEN
381 32508 : embedding(GEN g, GEN DATA, primedata *S, GEN den, GEN listdelta)
382 : {
383 32508 : GEN TR, w0_Q, w0, w1_Q, w1, wpow, h0, gp, T, q2, q, maxp, a, p = S->p;
384 : long rt;
385 : pari_sp av;
386 :
387 32508 : T = gel(DATA,1); rt = brent_kung_optpow(degpol(T), 4, 3);
388 32508 : maxp= gel(DATA,7);
389 32508 : gp = RgX_deriv(g); av = avma;
390 32505 : w0 = chinese_retrieve_pol(DATA, S, listdelta);
391 32506 : w0_Q = centermod(gmul(w0,den), p);
392 32508 : h0 = FpXQ_inv(FpX_FpXQ_eval(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
393 32506 : wpow = NULL; q = sqri(p);
394 : for(;;)
395 : {/* Given g,w0,h0 in Z[x], s.t. h0.g'(w0) = 1 and g(w0) = 0 mod (T,p), find
396 : * [w1,h1] satisfying the same conditions mod p^2, [w1,h1] = [w0,h0] (mod p)
397 : * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
398 115849 : if (DEBUGLEVEL>1)
399 0 : err_printf("lifting embedding mod p^k = %Ps^%ld\n",S->p, Z_pval(q,S->p));
400 :
401 : /* w1 := w0 - h0 g(w0) mod (T,q) */
402 115849 : if (wpow) a = FpX_FpXQV_eval(g,wpow, T,q);
403 32504 : else a = FpX_FpXQ_eval(g,w0, T,q); /* first time */
404 : /* now, a = 0 (p) */
405 115861 : a = FpXQ_mul(ZX_neg(h0), ZX_Z_divexact(a, p), T,p);
406 115855 : w1 = ZX_add(w0, ZX_Z_mul(a, p));
407 :
408 115847 : w1_Q = centermod(ZX_Z_mul(w1, remii(den,q)), q);
409 115859 : if (ZX_equal(w1_Q, w0_Q))
410 : {
411 21734 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
412 21734 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
413 : }
414 94128 : else if (cmpii(q,maxp) > 0)
415 : {
416 14971 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
417 14973 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
418 497 : if (DEBUGLEVEL) err_printf("coeff too big for embedding\n");
419 497 : return NULL;
420 : }
421 83356 : gerepileall(av, 5, &w1,&h0,&w1_Q,&q,&p);
422 83362 : q2 = sqri(q);
423 83347 : wpow = FpXQ_powers(w1, rt, T, q2);
424 : /* h0 := h0 * (2 - h0 g'(w1)) mod (T,q)
425 : * = h0 + h0 * (1 - h0 g'(w1)) */
426 83349 : a = FpXQ_mul(ZX_neg(h0), FpX_FpXQV_eval(gp, FpXV_red(wpow,q),T,q), T,q);
427 83351 : a = ZX_Z_add_shallow(a, gen_1); /* 1 - h0 g'(w1) = 0 (p) */
428 83351 : a = FpXQ_mul(h0, ZX_Z_divexact(a, p), T,p);
429 83355 : h0 = ZX_add(h0, ZX_Z_mul(a, p));
430 83345 : w0 = w1; w0_Q = w1_Q; p = q; q = q2;
431 : }
432 32011 : TR = gel(DATA,5);
433 32011 : if (!gequal0(TR)) w1_Q = RgX_translate(w1_Q, TR);
434 32011 : return gdiv(w1_Q,den);
435 : }
436 :
437 : /* return U list of polynomials s.t U[i] = 1 mod fk[i] and 0 mod fk[j] for all
438 : * other j */
439 : static GEN
440 30840 : get_bezout(GEN pol, GEN fk, GEN p)
441 : {
442 30840 : long i, l = lg(fk);
443 30840 : GEN A, B, d, u, v, U = cgetg(l, t_VEC);
444 63335 : for (i=1; i<l; i++)
445 : {
446 32493 : A = gel(fk,i);
447 32493 : B = FpX_div(pol, A, p);
448 32493 : d = FpX_extgcd(A,B,p, &u, &v);
449 32494 : if (degpol(d) > 0) pari_err_COPRIME("get_bezout",A,B);
450 32494 : d = constant_coeff(d);
451 32494 : if (!gequal1(d)) v = FpX_Fp_div(v, d, p);
452 32493 : gel(U,i) = FpX_mul(B,v, p);
453 : }
454 30842 : return U;
455 : }
456 :
457 : static GEN
458 30841 : init_traces(GEN ff, GEN T, GEN p)
459 : {
460 30841 : long N = degpol(T),i,j,k, r = lg(ff);
461 30841 : GEN Frob = FpX_matFrobenius(T,p);
462 : GEN y,p1,p2,Trk,pow,pow1;
463 :
464 30842 : k = degpol(gel(ff,r-1)); /* largest degree in modular factorization */
465 30842 : pow = cgetg(k+1, t_VEC);
466 30842 : gel(pow,1) = gen_0; /* dummy */
467 30842 : gel(pow,2) = Frob;
468 30842 : pow1= cgetg(k+1, t_VEC); /* 1st line */
469 110614 : for (i=3; i<=k; i++)
470 79772 : gel(pow,i) = FpM_mul(gel(pow,i-1), Frob, p);
471 30842 : gel(pow1,1) = gen_0; /* dummy */
472 141456 : for (i=2; i<=k; i++)
473 : {
474 110614 : p1 = cgetg(N+1, t_VEC);
475 110614 : gel(pow1,i) = p1; p2 = gel(pow,i);
476 778372 : for (j=1; j<=N; j++) gel(p1,j) = gcoeff(p2,1,j);
477 : }
478 :
479 : /* Trk[i] = line 1 of x -> x + x^p + ... + x^{p^(i-1)} */
480 30842 : Trk = pow; /* re-use (destroy) pow */
481 30842 : gel(Trk,1) = vec_ei(N,1);
482 141446 : for (i=2; i<=k; i++)
483 110608 : gel(Trk,i) = gadd(gel(Trk,i-1), gel(pow1,i));
484 30838 : y = cgetg(r, t_VEC);
485 63332 : for (i=1; i<r; i++) y[i] = Trk[degpol(gel(ff,i))];
486 30840 : return y;
487 : }
488 :
489 : static void
490 30838 : init_primedata(primedata *S)
491 : {
492 30838 : long i, j, l, lff = lg(S->ff), v = fetch_var(), N = degpol(S->pol);
493 30840 : GEN T, p = S->p;
494 :
495 30840 : if (S->lcm == degpol(gel(S->ff,lff-1)))
496 : {
497 30826 : T = leafcopy(gel(S->ff,lff-1));
498 30828 : setvarn(T, v);
499 : }
500 : else
501 14 : T = init_Fq(p, S->lcm, v);
502 30842 : S->T = T;
503 30842 : S->firstroot = cgetg(lff, t_VECSMALL);
504 30842 : S->interp = cgetg(lff, t_VEC);
505 30842 : S->fk = cgetg(N+1, t_VEC);
506 63334 : for (l=1,j=1; j<lff; j++)
507 : { /* compute roots and fix ordering (Frobenius cycles) */
508 32494 : GEN F = gel(S->ff, j), deg1 = FpX_factorff_irred(F, T,p);
509 32494 : GEN H = gel(deg1,1), a = Fq_neg(constant_coeff(H), T,p);
510 32493 : GEN Q = FqX_div(F, H, T,p);
511 32493 : GEN q = Fq_inv(FqX_eval(Q, a, T,p), T,p);
512 32494 : gel(S->interp,j) = FqX_Fq_mul(Q, q, T,p); /* = 1 at a, 0 at other roots */
513 32492 : S->firstroot[j] = l;
514 182914 : for (i=1; i<lg(deg1); i++,l++) gel(S->fk, l) = gel(deg1, i);
515 : }
516 30840 : S->Trk = init_traces(S->ff, T,p);
517 30840 : S->bezoutC = get_bezout(S->pol, S->ff, p);
518 30842 : }
519 :
520 : static int
521 30855 : choose_prime(primedata *S, GEN pol)
522 : {
523 30855 : long i, j, k, r, lcm, oldr, K, N = degpol(pol);
524 : ulong p, pp;
525 : GEN Z, ff, n, oldn;
526 : pari_sp av;
527 : forprime_t T;
528 :
529 30855 : u_forprime_init(&T, (N*N) >> 2, ULONG_MAX);
530 30856 : oldr = S->lcm = LONG_MAX;
531 30856 : S->ff = oldn = NULL; pp = 0; /* gcc -Wall */
532 30856 : av = avma; K = N + 10;
533 99710 : for(k = 1; k < K || !S->ff; k++,set_avma(av))
534 : {
535 : GEN Tp;
536 98253 : if (k > 5 * N) return 0;
537 : do
538 : {
539 116650 : p = u_forprime_next(&T);
540 116650 : Tp = ZX_to_Flx(pol, p);
541 : }
542 116650 : while (!Flx_is_squarefree(Tp, p));
543 98254 : ff = gel(Flx_factor(Tp, p), 1);
544 98255 : r = lg(ff)-1;
545 98255 : if (r == N || r >= BIL) continue;
546 :
547 95805 : n = cgetg(r+1, t_VECSMALL); lcm = n[1] = degpol(gel(ff,1));
548 256044 : for (j=2; j<=r; j++) { n[j] = degpol(gel(ff,j)); lcm = ulcm(lcm, n[j]); }
549 95804 : if (r > oldr || (r == oldr && (lcm <= S->lcm || S->lcm > 2*N)))
550 45454 : continue;
551 50350 : if (DEBUGLEVEL) err_printf("p = %lu,\tlcm = %ld,\torbits: %Ps\n",p,lcm,n);
552 :
553 50350 : pp = p;
554 50350 : oldr = r;
555 50350 : oldn = n;
556 50350 : S->ff = ff;
557 50350 : S->lcm = lcm; if (r == 1) break;
558 20950 : av = avma;
559 : }
560 30856 : if (oldr > 6) return 0;
561 30842 : if (DEBUGLEVEL) err_printf("Chosen prime: p = %ld\n", pp);
562 30842 : FlxV_to_ZXV_inplace(S->ff);
563 30842 : S->p = utoipos(pp);
564 30842 : S->pol = FpX_red(pol, S->p); init_primedata(S);
565 30842 : n = oldn; r = lg(n); S->Z = Z = cgetg(r,t_VEC);
566 63336 : for (k=0,i=1; i<r; i++)
567 : {
568 32494 : GEN t = cgetg(n[i]+1, t_VECSMALL); gel(Z,i) = t;
569 182924 : for (j=1; j<=n[i]; j++) t[j] = ++k;
570 : }
571 30842 : return 1;
572 : }
573 :
574 : /* maxroot t_REAL */
575 : static GEN
576 31947 : bound_for_coeff(long m, GEN R, GEN *maxroot)
577 : {
578 31947 : GEN b1, b2, M, v, C = vecbinomial(m-1);
579 31948 : long i, r1, l = lg(R);
580 :
581 48594 : for (r1 = 1; r1 < l; r1++)
582 47145 : if (typ(gel(R,r1)) != t_REAL) break;
583 31948 : r1--;
584 31948 : R = gabs(R,0); *maxroot = vecmax(R);
585 48594 : for (b1 = gen_1, i = 1; i <= r1; i++)
586 16646 : if (gcmpgs(gel(R,i), 1) > 0) b1 = gmul(b1, gel(R,i));
587 179566 : for (b2 = gen_1 ; i < l; i++)
588 147621 : if (gcmpgs(gel(R,i), 1) > 0) b2 = gmul(b2, gel(R,i));
589 31945 : M = gmul(b1, gsqr(b2)); /* Mahler measure */
590 31948 : v = cgetg(m+2, t_VEC); gel(v,1) = gel(v,2) = gen_0; /* unused */
591 79783 : for (i = 1; i < m; i++) /* binom(m-1, i) * M + binom(m-1, i-1) */
592 47838 : gel(v, i+2) = ceil_safe(gadd(gmul(gel(C, i+1), M), gel(C, i)));
593 31945 : return v;
594 : }
595 :
596 : static GEN
597 9527 : RgV_translate(GEN x, GEN t) { pari_APPLY_same(gadd(t, gel(x,i))); }
598 : static GEN
599 9527 : RgV_negtranslate(GEN x, GEN t) { pari_APPLY_same(gsub(t, gel(x,i))); }
600 :
601 : /* d = requested degree for subfield. Return DATA, valid for given pol, S and d
602 : * If DATA != NULL, translate pol [ --> pol(X+1) ] and update DATA
603 : * 1: polynomial pol
604 : * 2: p^e (for Hensel lifts) such that p^e > max(M),
605 : * 3: Hensel lift to precision p^e of DATA[4]
606 : * 4: roots of pol in F_(p^S->lcm),
607 : * 5: number of polynomial changes (translations)
608 : * 6: Bezout coefficients attached to the S->ff[i]
609 : * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod pol.
610 : * 8: bound M for polynomials defining subfields x PD->den
611 : * 9: *[i] = interpolation polynomial for S->ff[i] [= 1 on the first root
612 : S->firstroot[i], 0 on the others] */
613 : static void
614 31947 : compute_data(blockdata *B)
615 : {
616 : GEN ffL, roo, pe, p1, p2, fk, MM, maxroot, pol;
617 31947 : primedata *S = B->S;
618 31947 : GEN p = S->p, T = S->T, ff = S->ff, DATA = B->DATA;
619 31947 : long i, l, e, N, lff = lg(ff);
620 :
621 31947 : if (DEBUGLEVEL>1) err_printf("Entering compute_data()\n\n");
622 31947 : pol = B->PD->pol; N = degpol(pol);
623 31947 : roo = B->PD->roo;
624 31947 : if (DATA)
625 : {
626 763 : GEN TR = addiu(gel(DATA,5), 1), mTR = negi(TR), interp, bezoutC;
627 :
628 763 : if (DEBUGLEVEL>1) err_printf("... update (translate) an existing DATA\n\n");
629 763 : gel(DATA,5) = TR;
630 763 : pol = RgX_translate(gel(DATA,1), gen_m1);
631 763 : roo = RgV_translate(roo, TR);
632 763 : fk = RgV_negtranslate(gel(DATA,4),
633 763 : deg1pol_shallow(gen_1, gen_m1, varn(pol)));
634 763 : bezoutC = gel(DATA,6); l = lg(bezoutC);
635 763 : interp = gel(DATA,9);
636 2212 : for (i=1; i<l; i++)
637 : {
638 1449 : if (degpol(gel(interp,i)) > 0) /* do not turn pol_1(0) into gen_1 */
639 1449 : gel(interp,i) = FpXX_red(RgX_translate(gel(interp,i), gen_m1), p);
640 1449 : if (degpol(gel(bezoutC,i)) > 0)
641 1351 : gel(bezoutC,i)= FpXX_red(RgX_translate(gel(bezoutC,i), gen_m1), p);
642 : }
643 763 : ff = cgetg(lff, t_VEC); /* copy, do not overwrite! */
644 2212 : for (i=1; i<lff; i++)
645 1449 : gel(ff,i) = FpX_red(RgX_translate(gel(S->ff,i), mTR), p);
646 : }
647 : else
648 : {
649 31184 : DATA = cgetg(10,t_VEC);
650 31184 : fk = S->fk;
651 31184 : gel(DATA,5) = gen_0;
652 31184 : gel(DATA,6) = leafcopy(S->bezoutC);
653 31184 : gel(DATA,9) = leafcopy(S->interp);
654 : }
655 31947 : gel(DATA,1) = pol;
656 31947 : MM = gmul2n(bound_for_coeff(B->d, roo, &maxroot), 1);
657 31945 : gel(DATA,8) = MM;
658 31945 : e = logintall(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [d-1 test] */
659 31947 : gel(DATA,2) = pe;
660 31947 : gel(DATA,4) = roots_from_deg1(fk);
661 :
662 : /* compute fhk = ZpX_liftfact(pol,fk,T,p,e,pe) in 2 steps
663 : * 1) lift in Zp to precision p^e */
664 31947 : ffL = ZpX_liftfact(pol, ff, pe, p, e);
665 66415 : for (l=i=1; i<lff; i++)
666 : { /* 2) lift factorization of ff[i] in Qp[X] / T */
667 34467 : long l2 = l + degpol(gel(ffL,i));
668 34467 : gel(ffL,i) = ZqX_liftfact(gel(ffL,i), vecslice(fk, l, l2-1), T, pe, p, e);
669 34468 : l = l2;
670 : }
671 31948 : gel(DATA,3) = roots_from_deg1(shallowconcat1(ffL));
672 :
673 31947 : p1 = mulur(N, powruhalf(utor(N-1,DEFAULTPREC), N-1));
674 31945 : p2 = powru(maxroot, B->size + N*(N-1)/2);
675 31944 : p1 = divrr(mulrr(p1,p2), gsqrt(B->PD->dis,DEFAULTPREC));
676 31948 : gel(DATA,7) = mulii(shifti(ceil_safe(p1), 1), B->PD->den);
677 :
678 31946 : if (DEBUGLEVEL>1) {
679 0 : err_printf("f = %Ps\n",DATA[1]);
680 0 : err_printf("p = %Ps, lift to p^%ld\n", p, e);
681 0 : err_printf("2 * Hadamard bound * ind = %Ps\n",DATA[7]);
682 0 : err_printf("2 * M = %Ps\n",DATA[8]);
683 : }
684 31946 : if (B->DATA) { DATA = gclone(DATA); if (isclone(B->DATA)) gunclone(B->DATA); }
685 31946 : B->DATA = DATA;
686 31946 : }
687 :
688 : /* g = polynomial, h = embedding. Return [[g,h]] */
689 : static GEN
690 1463 : _subfield(GEN g, GEN h) { return mkvec(mkvec2(g,h)); }
691 :
692 : /* Return a subfield, gen_0 [ change p ] or NULL [ not a subfield ] */
693 : static GEN
694 54278 : subfield(GEN A, blockdata *B)
695 : {
696 54278 : long N, i, j, d, lf, m = lg(A)-1;
697 : GEN M, pe, pol, fhk, g, e, d_1_term, delta, listdelta, whichdelta;
698 54278 : GEN T = B->S->T, p = B->S->p, firstroot = B->S->firstroot;
699 :
700 54278 : pol= gel(B->DATA,1); N = degpol(pol); d = N/m; /* m | N */
701 54278 : pe = gel(B->DATA,2);
702 54278 : fhk= gel(B->DATA,3);
703 54278 : M = gel(B->DATA,8);
704 :
705 54278 : delta = cgetg(m+1,t_VEC);
706 54278 : whichdelta = cgetg(N+1, t_VECSMALL);
707 54284 : d_1_term = gen_0;
708 344627 : for (i=1; i<=m; i++)
709 : {
710 290351 : GEN Ai = gel(A,i), p1 = gel(fhk,Ai[1]);
711 902516 : for (j=2; j<=d; j++)
712 612169 : p1 = Fq_mul(p1, gel(fhk,Ai[j]), T, pe);
713 290347 : gel(delta,i) = p1;
714 290347 : if (DEBUGLEVEL>5) err_printf("delta[%ld] = %Ps\n",i,p1);
715 : /* g = prod (X - delta[i])
716 : * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
717 : /* fk[k] belongs to block number whichdelta[k] */
718 1192859 : for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
719 290347 : if (typ(p1) == t_POL) p1 = constant_coeff(p1);
720 290347 : d_1_term = addii(d_1_term, p1);
721 : }
722 54276 : d_1_term = centermod(d_1_term, pe); /* Tr(g) */
723 54276 : if (abscmpii(d_1_term, gel(M,3)) > 0) {
724 13867 : if (DEBUGLEVEL>1) err_printf("d-1 test failed\n");
725 13867 : return NULL;
726 : }
727 40409 : g = FqV_roots_to_pol(delta, T, pe, 0);
728 40411 : g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
729 40411 : if (!ok_coeffs(g,M)) {
730 7140 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
731 7140 : if (DEBUGLEVEL>1) err_printf("coeff too big for pol g(x)\n");
732 7140 : return NULL;
733 : }
734 33271 : if (!FpX_is_squarefree(g, p)) {
735 763 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
736 763 : if (DEBUGLEVEL>1) err_printf("changing f(x): p divides disc(g)\n");
737 763 : compute_data(B);
738 763 : return subfield(A, B);
739 : }
740 :
741 32508 : lf = lg(firstroot); listdelta = cgetg(lf, t_VEC);
742 69090 : for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
743 32508 : if (DEBUGLEVEL) err_printf("candidate = %Ps\n", g);
744 32508 : e = embedding(g, B->DATA, B->S, B->PD->den, listdelta);
745 32508 : if (!e) return NULL;
746 32011 : if (DEBUGLEVEL) err_printf("... OK!\n");
747 32011 : return B->fl==1? mkvec(g):_subfield(g, e);
748 : }
749 :
750 : /* L list of current subfields, test whether potential block D is a block,
751 : * if so, append corresponding subfield */
752 : static GEN
753 53515 : test_block(blockdata *B, GEN L, GEN D)
754 : {
755 53515 : pari_sp av = avma;
756 53515 : GEN sub = subfield(D, B);
757 53515 : if (sub) {
758 32011 : GEN old = L;
759 32011 : L = gclone( L? shallowconcat(L, sub): sub );
760 32010 : guncloneNULL(old);
761 : }
762 53514 : return gc_const(av,L);
763 : }
764 :
765 : /* subfields of degree d */
766 : static GEN
767 31184 : subfields_of_given_degree(blockdata *B)
768 : {
769 31184 : pari_sp av = avma;
770 : GEN L;
771 :
772 31184 : if (DEBUGLEVEL) err_printf("\n* Look for subfields of degree %ld\n\n", B->d);
773 31184 : B->DATA = NULL; compute_data(B);
774 31183 : L = calc_block(B, B->S->Z, cgetg(1,t_VEC), NULL);
775 31184 : if (DEBUGLEVEL>9)
776 0 : err_printf("\nSubfields of degree %ld: %Ps\n", B->d, L? L: cgetg(1,t_VEC));
777 31184 : if (isclone(B->DATA)) gunclone(B->DATA);
778 31184 : return gc_const(av,L);
779 : }
780 :
781 : static void
782 35 : setvarn2(GEN t, long v) { setvarn(gel(t,1),v); setvarn(gel(t,2),v); }
783 : static GEN
784 29820 : fix_var(GEN x, long v, long fl)
785 : {
786 29820 : long i, l = lg(x);
787 29820 : if (!v) return x;
788 28 : if (fl)
789 42 : for (i = 1; i < l; i++) setvarn(gel(x,i), v);
790 : else
791 49 : for (i = 1; i < l; i++) setvarn2(gel(x,i), v);
792 28 : return x;
793 : }
794 :
795 : static void
796 30842 : subfields_poldata(GEN nf, GEN T, poldata *PD)
797 : {
798 : GEN L, dis;
799 :
800 30842 : PD->pol = T;
801 30842 : if (nf)
802 : {
803 168 : PD->den = nf_get_zkden(nf);
804 168 : PD->roo = nf_get_roots(nf);
805 168 : PD->dis = mulii(absi_shallow(nf_get_disc(nf)), sqri(nf_get_index(nf)));
806 : }
807 : else
808 : {
809 30674 : PD->den = initgaloisborne(T,NULL,nbits2prec(bit_accuracy(ZX_max_lg(T))), &L,NULL,&dis);
810 30673 : PD->roo = L;
811 30673 : PD->dis = absi_shallow(dis);
812 : }
813 30841 : }
814 :
815 : static GEN nfsubfields_fa(GEN nf, long d, long fl);
816 : static GEN
817 595 : subfieldsall(GEN nf0, long fl)
818 : {
819 595 : pari_sp av = avma;
820 : long N, ld, i, v;
821 : GEN nf, G, T, dg, LSB, NLSB;
822 : poldata PD;
823 : primedata S;
824 : blockdata B;
825 :
826 : /* much easier if nf is Galois (WSS) */
827 595 : G = galoisinit(nf0, NULL);
828 595 : T = get_nfpol(nf0, &nf);
829 595 : if (G != gen_0)
830 : {
831 : GEN L, S;
832 : long l;
833 56 : L = lift_shallow( galoissubfields(G, fl, varn(T)) );
834 56 : l = lg(L); S = cgetg(l, t_VECSMALL);
835 924 : for (i=1; i<l; i++) S[i] = lg(fl==1? gel(L,i): gmael(L,i,1));
836 56 : return gerepilecopy(av, vecpermute(L, vecsmall_indexsort(S)));
837 : }
838 539 : v = varn(T); N = degpol(T);
839 539 : dg = divisorsu(N); ld = lg(dg)-1;
840 539 : LSB = fl==1 ? mkvec(pol_x(v)): _subfield(pol_x(v), pol_0(v));
841 539 : if (ld <= 2)
842 : {
843 168 : if (ld == 2)
844 168 : LSB = shallowconcat(LSB, fl==1? mkvec(T): _subfield(T, pol_x(v)));
845 168 : return gerepilecopy(av, LSB);
846 : }
847 371 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
848 371 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, 0, fl); }
849 357 : subfields_poldata(nf, T, &PD);
850 :
851 357 : if (DEBUGLEVEL) err_printf("\n***** Entering subfields\n\npol = %Ps\n",T);
852 357 : B.PD = &PD;
853 357 : B.S = &S;
854 357 : B.N = N;
855 357 : B.fl = fl;
856 1057 : for (i=ld-1; i>1; i--)
857 : {
858 700 : B.size = dg[i];
859 700 : B.d = N / B.size;
860 700 : NLSB = subfields_of_given_degree(&B);
861 700 : if (NLSB) { LSB = gconcat(LSB, NLSB); gunclone(NLSB); }
862 : }
863 357 : (void)delete_var(); /* from init_primedata() */
864 357 : LSB = shallowconcat(LSB, fl==1? mkvec(T):_subfield(T, pol_x(0)));
865 357 : if (DEBUGLEVEL) err_printf("\n***** Leaving subfields\n\n");
866 357 : return fix_var(gerepilecopy(av, LSB), v, fl);
867 : }
868 :
869 : GEN
870 32751 : nfsubfields0(GEN nf0, long d, long fl)
871 : {
872 32751 : pari_sp av = avma;
873 : long N, v0;
874 : GEN nf, LSB, T, G;
875 : poldata PD;
876 : primedata S;
877 : blockdata B;
878 32751 : if (fl<0 || fl>1) pari_err_FLAG("nfsubfields");
879 32752 : if (typ(nf0)==t_VEC && lg(nf0)==3) return nfsubfields_fa(nf0, d, fl);
880 32507 : if (!d) return subfieldsall(nf0, fl);
881 :
882 : /* treat trivial cases */
883 31912 : T = get_nfpol(nf0, &nf); v0 = varn(T); N = degpol(T);
884 31911 : RgX_check_ZX(T,"nfsubfields");
885 31911 : if (d == N)
886 28 : return gerepilecopy(av, fl==1 ? mkvec(T) : _subfield(T, pol_x(v0)));
887 31883 : if (d == 1)
888 28 : return gerepilecopy(av, fl==1 ? mkvec(pol_x(v0)) : _subfield(pol_x(v0), zeropol(v0)));
889 31855 : if (d < 1 || d > N || N % d) return cgetg(1,t_VEC);
890 :
891 : /* much easier if nf is Galois (WSS) */
892 31813 : G = galoisinit(nf0, NULL);
893 31814 : if (G != gen_0)
894 : { /* Bingo */
895 1330 : GEN L = galoissubgroups(G), F;
896 1330 : long k,i, l = lg(L), o = N/d;
897 1330 : F = cgetg(l, t_VEC);
898 1330 : k = 1;
899 5936 : for (i=1; i<l; i++)
900 : {
901 4606 : GEN H = gel(L,i);
902 4606 : if (group_order(H) == o)
903 1540 : gel(F,k++) = lift_shallow(galoisfixedfield(G, gel(H,1), fl, v0));
904 : }
905 1330 : setlg(F, k);
906 1330 : return gerepilecopy(av, F);
907 : }
908 30484 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
909 30484 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, d, fl); }
910 30485 : subfields_poldata(nf, T, &PD);
911 30484 : B.PD = &PD;
912 30484 : B.S = &S;
913 30484 : B.N = N;
914 30484 : B.d = d;
915 30484 : B.size = N/d;
916 30484 : B.fl = fl;
917 30484 : LSB = subfields_of_given_degree(&B);
918 30484 : (void)delete_var(); /* from init_primedata */
919 30484 : set_avma(av);
920 30484 : if (!LSB) return cgetg(1, t_VEC);
921 29462 : G = gcopy(LSB); gunclone(LSB);
922 29463 : return fix_var(G, v0, fl);
923 : }
924 :
925 : GEN
926 329 : nfsubfields(GEN nf0, long d)
927 329 : { return nfsubfields0(nf0, d, 0); }
928 :
929 : /******************************/
930 : /* */
931 : /* Maximal CM subfield */
932 : /* Aurel Page (2019) */
933 : /* */
934 : /******************************/
935 :
936 : /* ero: maximum exponent+1 of roots of pol */
937 : static GEN
938 3325 : try_subfield_generator(GEN pol, GEN v, long e, long p, long ero, long fl)
939 : {
940 : GEN a, P, Q;
941 : long d, bound, i, B, bi, ed;
942 :
943 3325 : a = gtopolyrev(v, varn(pol));
944 3325 : P = Flxq_charpoly(ZX_to_Flx(a,p), ZX_to_Flx(pol,p), p);
945 3325 : Flx_ispower(P, e, p, &Q);
946 3325 : if (!Flx_is_squarefree(Q,p)) return NULL;
947 1701 : d = degpol(pol)/e;
948 1701 : B = 0;
949 26607 : for (i=1; i<lg(v); i++)
950 : {
951 24906 : bi = (i-1)*ero + expi(gel(v,i));
952 24906 : if (bi > B) B = bi;
953 : }
954 1701 : ed = expu(d);
955 1701 : B += ed+1;
956 1701 : bound = 0;
957 8253 : for (i=0; 2*i<=d; i++)
958 : {
959 6552 : if (!i) bi = d*B;
960 4851 : else bi = (d-i)*B + i*(3+ed-expu(i));
961 6552 : if (bi > bound) bound = bi;
962 : }
963 1701 : Q = ZXQ_minpoly(a,pol,d,bound);
964 1701 : return fl==1? Q: mkvec2(Q, a);
965 : }
966 :
967 : /* subfield sub of nf of degree d assuming:
968 : - V is contained in sub
969 : - V is not contained in a proper subfield of sub
970 : ero: maximum exponent+1 of roots of pol
971 : output as nfsubfields:
972 : - pair [g,h] where g absolute equation for the subfield and h expresses
973 : - one of the roots of g in terms of the generator of nf
974 : */
975 : static GEN
976 1876 : subfield_generator(GEN pol, GEN V, long d, long ero, long fl)
977 : {
978 1876 : long p, i, e, vp = varn(pol);
979 1876 : GEN a = NULL, v = cgetg(lg(V),t_COL), B;
980 :
981 1876 : if (d==1) return fl ? pol_x(vp): mkvec2(pol_x(vp), pol_0(vp));
982 1701 : e = degpol(pol)/d;
983 1701 : p = 1009;
984 3318 : for (i=1; i<lg(V); i++)
985 : {
986 3297 : a = try_subfield_generator(pol, gel(V,i), e, p, ero, fl);
987 3297 : if (a) return a;
988 1617 : p = unextprime(p+1);
989 : }
990 21 : B = stoi(10);
991 : while(1)
992 : {
993 126 : for (i=1; i<lg(v); i++) gel(v,i) = randomi(B);
994 28 : a = try_subfield_generator(pol, QM_QC_mul(V,v), e, p, ero, fl);
995 28 : if (a) return a;
996 7 : p = unextprime(p+1);
997 : }
998 : return NULL;/*LCOV_EXCL_LINE*/
999 : }
1000 :
1001 : static GEN
1002 37961 : RgXY_to_RgC(GEN P, long dx, long dy)
1003 : {
1004 : GEN res, c;
1005 37961 : long i, j, k, d = degpol(P);
1006 37961 : if (d > dy) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1007 37961 : res = cgetg((dx+1)*(dy+1)+1, t_COL);
1008 37961 : k = 1;
1009 93618 : for (i=0; i<=d; i++)
1010 : {
1011 55657 : c = gel(P,i+2);
1012 55657 : if (typ(c)==t_POL)
1013 : {
1014 51408 : long dc = degpol(c);
1015 51408 : if (dc > dx) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1016 1021790 : for (j=0; j<=dc; j++)
1017 970382 : gel(res,k++) = gel(c,j+2);
1018 : } else
1019 : {
1020 4249 : gel(res,k++) = c; j=1;
1021 : }
1022 344365 : for ( ; j<=dx; j++)
1023 288708 : gel(res,k++) = gen_0;
1024 : }
1025 87304 : for( ; i<=dy; i++)
1026 1075242 : for (j=0; j<=dx; j++)
1027 1025899 : gel(res,k++) = gen_0;
1028 37961 : return res;
1029 : }
1030 :
1031 : /* lambda: t_VEC of t_INT; 0 means ignore this factor */
1032 : static GEN
1033 2583 : twoembequation(GEN pol, GEN fa, GEN lambda)
1034 : {
1035 : GEN m, vpolx, poly;
1036 2583 : long i,j, lfa = lg(fa), dx = degpol(pol);
1037 2583 : long vx = varn(pol), vy = varn(gel(fa,1)); /* vx < vy ! */
1038 :
1039 2583 : if (varncmp(vx,vy) <= 0) pari_err_BUG("twoembequation [incorrect variable priorities]");
1040 :
1041 2583 : lambda = shallowcopy(lambda);
1042 2583 : fa = shallowcopy(fa);
1043 2583 : j = 1;
1044 27839 : for (i=1; i<lfa; i++)
1045 25256 : if (signe(gel(lambda,i)))
1046 : {
1047 2618 : gel(lambda,j) = negi(gel(lambda,i));
1048 2618 : gel(fa,j) = gel(fa,i);
1049 2618 : j++;
1050 : }
1051 2583 : setlg(lambda, j);
1052 2583 : setlg(fa, j); lfa = j;
1053 :
1054 2583 : vpolx = ZXQ_powers(pol_x(vx),dx-1,pol);
1055 2583 : m = cgetg(dx+1, t_MAT);
1056 40292 : for (j=1; j <= dx; j++)
1057 37709 : gel(m,j) = cgetg(lfa, t_COL);
1058 5201 : for(i=1; i<lfa; i++)
1059 : {
1060 2618 : long dy = degpol(gel(fa,i));
1061 2618 : poly = pol_1(vy);
1062 40579 : for (j=1; j <= dx; j++)
1063 : {
1064 37961 : gcoeff(m,i,j) = RgXY_to_RgC(gadd(ZX_Z_mul(gel(vpolx,j),gel(lambda,i)),poly), dx, dy);
1065 37961 : poly = RgXQX_rem(RgX_shift(poly,1), gel(fa,i), pol);
1066 : }
1067 : }
1068 40292 : for(j=1; j<=dx; j++) gel(m,j) = shallowconcat1(gel(m,j));
1069 2583 : return QM_ker(m);
1070 : }
1071 :
1072 : static void
1073 1561 : subfields_cleanup(GEN* nf, GEN* pol, long* n, GEN* fa)
1074 : {
1075 1561 : *fa = NULL;
1076 1561 : if (typ(*nf) != t_VEC && typ(*nf) != t_POL) pari_err_TYPE("subfields_cleanup", *nf);
1077 1554 : if (typ(*nf) == t_VEC && lg(*nf) == 3)
1078 : {
1079 301 : *fa = gel(*nf,2);
1080 301 : *nf = gel(*nf,1);
1081 301 : if (typ(*fa)!=t_MAT || lg(*fa)!=3)
1082 14 : pari_err_TYPE("subfields_cleanup [fa should be a factorisation matrix]", *fa);
1083 : }
1084 1540 : if (typ(*nf) == t_POL)
1085 : {
1086 784 : *pol = *nf;
1087 784 : *nf = NULL;
1088 784 : if (!RgX_is_ZX(*pol)) pari_err_TYPE("subfields_cleanup [not integral]", *pol);
1089 777 : if (!equali1(leading_coeff(*pol))) pari_err_TYPE("subfields_cleanup [not monic]", *pol);
1090 770 : *n = degpol(*pol);
1091 770 : if (*n<=0) pari_err_TYPE("subfields_cleanup [constant polynomial]", *pol);
1092 : }
1093 : else
1094 : {
1095 756 : *nf = checknf(*nf);
1096 735 : *pol = nf_get_pol(*nf);
1097 735 : *n = degpol(*pol);
1098 : }
1099 1498 : if(*fa)
1100 : {
1101 273 : long v = varn(*pol);
1102 273 : GEN o = gcoeff(*fa,1,1);
1103 273 : if (varncmp(varn(o),v) >= 0) pari_err_PRIORITY("nfsubfields_fa", o, "<=", v);
1104 : }
1105 1477 : }
1106 :
1107 : static GEN
1108 280 : rootsuptoconj(GEN pol, long prec)
1109 : {
1110 : GEN ro;
1111 : long n, i;
1112 280 : ro = roots(pol,prec);
1113 280 : n = lg(ro)-1;
1114 1498 : for (i=1; i<=n/2; i++)
1115 1218 : gel(ro,i) = gel(ro,2*i-1);
1116 280 : setlg(ro,n/2+1);
1117 280 : return ro;
1118 : }
1119 : static GEN
1120 1085 : cmsubfield_get_roots(GEN pol, GEN nf, long n, long* r2, long *prec)
1121 : {
1122 : GEN ro;
1123 1085 : if (nf)
1124 : {
1125 735 : if (nf_get_r1(nf)) return NULL;
1126 371 : *r2 = nf_get_r2(nf);
1127 371 : *prec = nf_get_prec(nf);
1128 371 : ro = nf_get_roots(nf);
1129 : }
1130 : else
1131 : {
1132 350 : if (n%2 || sturm(pol)) return NULL;
1133 280 : *r2 = n/2;
1134 280 : *prec = MEDDEFAULTPREC;
1135 280 : ro = rootsuptoconj(pol, *prec);
1136 : }
1137 651 : return ro;
1138 : }
1139 :
1140 : static GEN
1141 595 : subfields_get_fa(GEN pol, GEN nf, GEN fa)
1142 : {
1143 595 : if (!fa)
1144 : {
1145 392 : GEN poly = shallowcopy(pol);
1146 392 : setvarn(poly, fetch_var_higher());
1147 392 : fa = nffactor(nf? nf: pol, poly);
1148 : }
1149 595 : return liftpol_shallow(gel(fa,1));
1150 : }
1151 :
1152 : static long
1153 287 : subfields_get_ero(GEN pol, GEN nf)
1154 : {
1155 574 : return 1 + gexpo(nf? nf_get_roots(nf):
1156 287 : QX_complex_roots(pol, LOWDEFAULTPREC));
1157 : }
1158 :
1159 : static GEN
1160 280 : try_imag(GEN x, GEN c, GEN pol, long v, ulong p, GEN emb, GEN galpol, long fl)
1161 : {
1162 280 : GEN a = Q_primpart(RgX_sub(RgX_RgXQ_eval(x,c,pol),x));
1163 280 : if (Flx_is_squarefree(Flxq_charpoly(ZX_to_Flx(a,p),ZX_to_Flx(pol,p),p),p))
1164 : {
1165 168 : pol = ZXQ_charpoly(a, pol, v);
1166 168 : return fl ? pol : mkvec2(pol, RgX_RgXQ_eval(a, emb, galpol));
1167 : }
1168 112 : return NULL;
1169 : }
1170 :
1171 : static GEN
1172 210 : galoissubfieldcm(GEN G, long fl)
1173 : {
1174 210 : pari_sp av = avma;
1175 : GEN c, H, elts, g, Hset, c2, gene, sub, pol, emb, a, galpol, B, b;
1176 : long n, i, j, nH, ind, v, d;
1177 210 : ulong p = 1009;
1178 :
1179 210 : galpol = gal_get_pol(G);
1180 210 : n = degpol(galpol);
1181 210 : v = varn(galpol);
1182 210 : c = galois_get_conj(G);
1183 : /* compute the list of c*g*c*g^(-1) : product of all pairs of conjugations
1184 : * maximal CM subfield is the field fixed by those elements, if c does not
1185 : * belong to the group they generate */
1186 210 : checkgroup(G, &elts);
1187 210 : elts = gen_sort_shallow(elts,(void*)vecsmall_lexcmp,cmp_nodata);
1188 210 : H = vecsmall_ei(n,1); /* indices of elements of H */
1189 210 : Hset = zero_F2v(n);
1190 210 : F2v_set(Hset,1);
1191 210 : nH = 1;
1192 1456 : for (i=2; i<=n; i++)
1193 : {
1194 1246 : g = gel(elts,i);
1195 1246 : c2 = perm_mul(c,perm_conj(g,c));
1196 1246 : if (!F2v_coeff(Hset,c2[1]))
1197 : {
1198 182 : nH++;
1199 182 : H[nH] = c2[1];
1200 182 : F2v_set(Hset,c2[1]);
1201 : }
1202 : }
1203 : /* group generated */
1204 210 : gene = gcopy(H);
1205 210 : setlg(gene,nH+1);
1206 210 : i = 1; /* last element that has been multiplied by the generators */
1207 392 : while (i < nH)
1208 : {
1209 1218 : for (j=1; j<lg(gene); j++)
1210 : {
1211 1036 : g = gel(elts,gene[j]);
1212 1036 : ind = g[H[i]]; /* index of the product */
1213 1036 : if (!F2v_coeff(Hset,ind))
1214 : {
1215 0 : nH++;
1216 0 : if (ind==c[1] || 2*nH>n) return gc_const(av, gen_0);
1217 0 : H[nH] = ind;
1218 0 : F2v_set(Hset,ind);
1219 : }
1220 : }
1221 182 : i++;
1222 : }
1223 210 : H = cgetg(lg(gene), t_VEC);
1224 602 : for (i=1; i<lg(H); i++)
1225 392 : gel(H,i) = gel(elts,gene[i]);
1226 210 : sub = galoisfixedfield(G, H, 0, -1);
1227 :
1228 : /* compute a totally imaginary generator */
1229 210 : pol = gel(sub,1);
1230 210 : emb = liftpol_shallow(gel(sub,2));
1231 210 : d = degpol(pol);
1232 210 : if (!(ZX_deflate_order(pol)%2) && sturm(RgX_deflate(pol,2))==d/2)
1233 : {
1234 42 : setvarn(pol,v);
1235 42 : return fl==1 ? pol: mkvec2(pol,emb);
1236 : }
1237 :
1238 : /* compute action of c on the subfield from that on the large field */
1239 168 : c = galoispermtopol(G,c);
1240 168 : if (d<n)
1241 : {
1242 35 : GEN M = cgetg(d+1,t_MAT), contc, contM;
1243 35 : gel(M,1) = col_ei(n,1); a = pol_1(v);
1244 98 : for (i=2; i<=d; i++)
1245 : {
1246 63 : a = RgX_rem(QX_mul(a,emb), galpol);
1247 63 : gel(M,i) = RgX_to_RgC(a,n);
1248 : }
1249 35 : c = RgX_RgXQ_eval(emb,c,galpol);
1250 35 : c = Q_primitive_part(c,&contc);
1251 35 : c = RgX_to_RgC(c,n);
1252 35 : M = Q_primitive_part(M,&contM);
1253 35 : c = RgM_RgC_invimage(M,c);
1254 35 : if (contc)
1255 : {
1256 21 : if (contM) contc = gdiv(contc,contM);
1257 21 : c = RgV_Rg_mul(c, contc);
1258 : }
1259 14 : else if (contM) c = RgV_Rg_mul(c, ginv(contM));
1260 35 : c = RgV_to_RgX(c, v);
1261 : }
1262 :
1263 : /* search for a generator of the form c(b)-b */
1264 273 : for (i=1; i<d; i++)
1265 : {
1266 238 : a = try_imag(pol_xn(i,v),c,pol,v,p,emb,galpol,fl);
1267 238 : if (a) return a;
1268 105 : p = unextprime(p+1);
1269 : }
1270 35 : B = stoi(10);
1271 35 : b = pol_xn(d-1,v);
1272 : while(1)
1273 : {
1274 210 : for (i=2; i<lg(b); i++) gel(b,i) = randomi(B);
1275 42 : a = try_imag(b,c,pol,v,p,emb,galpol,fl);
1276 42 : if (a) return a;
1277 7 : p = unextprime(p+1);
1278 : }
1279 : return NULL;/*LCOV_EXCL_LINE*/
1280 : }
1281 :
1282 : static GEN
1283 133 : quadsubfieldcm(GEN pol, long fl)
1284 : {
1285 133 : GEN a = gel(pol,3), b = gel(pol,2), d, P;
1286 133 : long v = varn(pol);
1287 133 : if (mpodd(a))
1288 35 : { b = mului(4, b); d = gen_2; }
1289 : else
1290 98 : { a = divis(a,2); d = gen_1; }
1291 133 : P = deg2pol_shallow(gen_1, gen_0, subii(b, sqri(a)), v);
1292 133 : return fl==1 ? P: mkvec2(P, deg1pol_shallow(d,a,v));
1293 : }
1294 :
1295 : GEN
1296 1148 : nfsubfieldscm(GEN nf, long fl)
1297 : {
1298 1148 : pari_sp av = avma;
1299 : GEN fa, lambda, V, res, ro, a, aa, ev, minev, pol, G;
1300 1148 : long i, j, n, r2, minj=0, prec, emax, emin, e, precbound, ero;
1301 :
1302 1148 : subfields_cleanup(&nf, &pol, &n, &fa);
1303 1085 : ro = cmsubfield_get_roots(pol, nf, n, &r2, &prec);
1304 1085 : if (!ro) return gc_const(av, gen_0);
1305 : /* now r2 == 2*n */
1306 :
1307 651 : if (n==2) return gerepilecopy(av, quadsubfieldcm(pol, fl));
1308 518 : G = galoisinit(nf? nf: pol, NULL);
1309 518 : if (G != gen_0) return gerepilecopy(av, galoissubfieldcm(G, fl));
1310 :
1311 308 : ero = 0;
1312 1624 : for (i=1; i<lg(ro); i++)
1313 : {
1314 1316 : e = 1+gexpo(gel(ro,i));
1315 1316 : if (e > ero) ero = e;
1316 : }
1317 308 : ero++;
1318 308 : fa = subfields_get_fa(pol, nf, fa);
1319 :
1320 308 : emax = 1;
1321 308 : emin = -1;
1322 1624 : for (i=1; i<lg(ro); i++)
1323 4963 : for (j=i+1; j<lg(ro); j++)
1324 : {
1325 3647 : e = gexpo(gsub(gel(ro,i),gel(ro,j)));
1326 3647 : if (e > emax) emax = e;
1327 3647 : if (e < emin) emin = e;
1328 : }
1329 308 : precbound = n*(emax-emin) + gexpo(fa) + n*n + 5;
1330 308 : precbound = 3 + precbound/BITS_IN_LONG;
1331 308 : if (prec < precbound)
1332 : {
1333 0 : prec = precbound;
1334 0 : ro = rootsuptoconj(pol, prec);
1335 : }
1336 :
1337 308 : lambda = zerovec(lg(fa)-1);
1338 1624 : for (i=1; i<=r2; i++)
1339 : {
1340 1316 : a = gel(ro,i);
1341 1316 : aa = conj_i(a);
1342 9422 : for (j=1; j<lg(fa); j++)
1343 : {
1344 8106 : ev = cxnorm(poleval(poleval(gel(fa,j),aa),a));
1345 8106 : if (j==1 || cmprr(minev,ev)>0) { minj = j; minev = ev; }
1346 : }
1347 1316 : gel(lambda,minj) = gen_m1;
1348 : }
1349 :
1350 308 : V = twoembequation(pol, fa, lambda);
1351 308 : if (lg(V)==1) { delete_var(); return gc_const(av, gen_0); }
1352 259 : res = subfield_generator(pol, V, 2*(lg(V)-1), ero, fl);
1353 259 : delete_var();
1354 259 : return gerepilecopy(av, res);
1355 : }
1356 :
1357 : static int
1358 19474 : field_is_contained(GEN V, GEN W, int strict)
1359 : {
1360 : GEN VW;
1361 19474 : ulong p = 1073741827;
1362 : /* distinct overfield must have different dimension */
1363 19474 : if (strict && lg(V) == lg(W)) return 0;
1364 : /* dimension of overfield must be multiple */
1365 14553 : if ((lg(W)-1) % (lg(V)-1)) return 0;
1366 10402 : VW = shallowconcat(V,W);
1367 10402 : if (Flm_rank(ZM_to_Flm(VW,p),p) > lg(W)-1) return 0;
1368 4235 : return ZM_rank(VW) == lg(W)-1;
1369 : }
1370 :
1371 : /***********************************************/
1372 : /* */
1373 : /* Maximal, generating, all subfields */
1374 : /* Aurel Page (2019) */
1375 : /* after van Hoeij, Klueners, Novocin */
1376 : /* Journal of Symbolic Computation 52 (2013) */
1377 : /* */
1378 : /***********************************************/
1379 :
1380 : const long subf_MAXIMAL = 1; /* return the maximal subfields */
1381 : const long subf_GENERATING = 2; /* return the generating subfields */
1382 : static GEN
1383 287 : maxgen_subfields(GEN pol, GEN fa, long flag)
1384 : {
1385 287 : pari_sp av = avma;
1386 287 : GEN principal, ismax, isgene, Lmax = NULL, Lgene, res, V, W, W1;
1387 287 : long i, i2, j, flmax, flgene, nbmax = 0, nbgene = 0;
1388 :
1389 287 : if (!flag) return cgetg(1,t_VEC);
1390 287 : flmax = (flag & subf_MAXIMAL)!=0;
1391 287 : flgene = (flag & subf_GENERATING)!=0;
1392 :
1393 : /* compute principal subfields */
1394 287 : principal = cgetg(lg(fa),t_VEC);
1395 2562 : for (i=1; i<lg(fa); i++)
1396 2275 : gel(principal,i) = twoembequation(pol, fa, vec_ei(lg(fa)-1,i));
1397 287 : principal = gen_sort_uniq(principal, (void*)&cmp_universal, &cmp_nodata);
1398 : /* remove nf and duplicates (sort_uniq possibly not enough) */
1399 287 : i2 = 1;
1400 1694 : for (i=1; i<lg(principal)-1; i++)
1401 : {
1402 1407 : long dup = 0;
1403 1407 : V = gel(principal,i);
1404 1407 : j = i2-1;
1405 2877 : while (j > 0 && lg(gel(principal,j)) == lg(V))
1406 : {
1407 1470 : if (field_is_contained(gel(principal,j),V,0)) { dup=1; break; }
1408 1470 : j--;
1409 : }
1410 1407 : if (!dup) gel(principal,i2++) = V;
1411 : }
1412 287 : setlg(principal, i2);
1413 :
1414 : /* a subfield is generating iff all overfields contain the first overfield */
1415 287 : ismax = cgetg(lg(principal),t_VECSMALL);
1416 287 : isgene = cgetg(lg(principal),t_VECSMALL);
1417 1694 : for (i=1; i<lg(principal); i++)
1418 : {
1419 1407 : V = gel(principal,i);
1420 1407 : ismax[i] = flmax;
1421 1407 : isgene[i] = flgene;
1422 1407 : W1 = NULL; /* intersection of strict overfields */
1423 4858 : for (j=i+1; j<lg(principal); j++)
1424 : {
1425 3696 : W = gel(principal,j);
1426 3696 : if (!field_is_contained(V,W,1)) continue;
1427 693 : ismax[i] = 0;
1428 693 : if (!flgene) break;
1429 483 : if (!W1) { W1 = W; continue; }
1430 189 : if (!field_is_contained(W1,W,1))
1431 : {
1432 63 : W1 = intersect(W1,W);
1433 63 : if (lg(W1)==lg(V)) { isgene[i]=0; break; }
1434 : }
1435 : }
1436 : }
1437 :
1438 1694 : for (i=1; i<lg(principal); i++)
1439 : {
1440 1407 : nbmax += ismax[i];
1441 1407 : nbgene += isgene[i];
1442 : }
1443 :
1444 287 : if (flmax)
1445 : {
1446 98 : Lmax = cgetg(nbmax+1, t_VEC);
1447 98 : j=1;
1448 518 : for (i=1; i<lg(principal); i++)
1449 420 : if (ismax[i]) gel(Lmax,j++) = gel(principal,i);
1450 : }
1451 :
1452 287 : if (flgene)
1453 : {
1454 189 : Lgene = cgetg(nbgene+1, t_VEC);
1455 189 : j=1;
1456 1176 : for (i=1; i<lg(principal); i++)
1457 987 : if (isgene[i]) gel(Lgene,j++) = gel(principal,i);
1458 : }
1459 :
1460 287 : if (!flgene) res = Lmax;
1461 189 : else if (!flmax) res = Lgene;
1462 0 : else res = mkvec2(Lmax,Lgene);
1463 287 : return gerepilecopy(av, res);
1464 : }
1465 :
1466 : GEN
1467 154 : nfsubfieldsmax(GEN nf, long fl)
1468 : {
1469 154 : pari_sp av = avma;
1470 : GEN pol, fa, Lmax, V;
1471 : long n, i, ero;
1472 :
1473 154 : subfields_cleanup(&nf, &pol, &n, &fa);
1474 154 : if (n==1) { set_avma(av); return cgetg(1,t_VEC); }
1475 140 : if (uisprime(n))
1476 63 : return gerepilecopy(av, fl==1 ? mkvec(pol_x(varn(pol)))
1477 21 : : mkvec(mkvec2(pol_x(varn(pol)),gen_0)));
1478 98 : ero = subfields_get_ero(pol, nf);
1479 98 : fa = subfields_get_fa(pol, nf, fa);
1480 98 : Lmax = maxgen_subfields(pol, fa, subf_MAXIMAL);
1481 308 : for (i=1; i<lg(Lmax); i++)
1482 : {
1483 210 : V = gel(Lmax,i);
1484 210 : gel(Lmax,i) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1485 : }
1486 98 : delete_var();
1487 98 : return gerepilecopy(av, Lmax);
1488 : }
1489 :
1490 : static void
1491 1764 : heap_climb(GEN* H, long i)
1492 : {
1493 : long j;
1494 1764 : if (i==1) return;
1495 1302 : j = i/2;
1496 1302 : if (cmp_universal(gel(*H,i),gel(*H,j)) > 0)
1497 : {
1498 532 : swap(gel(*H,i), gel(*H,j));
1499 532 : return heap_climb(H,j);
1500 : }
1501 : }
1502 :
1503 : static void
1504 1232 : heap_push(GEN* H, long *len, GEN x)
1505 : {
1506 1232 : if (*len+1 == lg(*H))
1507 : {
1508 14 : GEN H2 = zerovec(2*(*len));
1509 : long i;
1510 154 : for(i=1; i<lg(*H); i++)
1511 140 : gel(H2,i) = gel(*H,i);
1512 14 : *H = H2;
1513 : }
1514 1232 : (*len)++;
1515 1232 : gel(*H,*len) = x;
1516 1232 : return heap_climb(H,*len);
1517 : }
1518 :
1519 : static void
1520 2569 : heap_descend(GEN* H, long len, long i)
1521 : {
1522 2569 : long maxi = i, j = 2*i;
1523 2569 : if (j > len) return;
1524 1337 : if (cmp_universal(gel(*H,j),gel(*H,i)) > 0) maxi = j;
1525 1337 : j++;
1526 1337 : if (j<=len && cmp_universal(gel(*H,j),gel(*H,maxi))>0) maxi = j;
1527 1337 : if (maxi == i) return;
1528 1148 : swap(gel(*H,i), gel(*H,maxi));
1529 1148 : return heap_descend(H,len,maxi);
1530 : }
1531 :
1532 : static void
1533 1421 : heap_pop(GEN *H, long *len, GEN* top)
1534 : {
1535 1421 : *top = gel(*H,1);
1536 1421 : gel(*H,1) = gel(*H,*len);
1537 1421 : (*len)--;
1538 1421 : return heap_descend(H,*len,1);
1539 : };
1540 :
1541 : static GEN
1542 259 : nfsubfields_fa(GEN nf, long d, long fl)
1543 : {
1544 259 : pari_sp av = avma;
1545 : GEN pol, fa, gene, res, res2, H, V, v, W, w, data;
1546 : long n, r, i, j, nres, len, s, newfield, ero, vp;
1547 :
1548 259 : subfields_cleanup(&nf, &pol, &n, &fa); vp = varn(pol);
1549 238 : if (d && (d<1 || d>n || n%d)) return gerepilecopy(av, cgetg(1,t_VEC));
1550 245 : if (!d && uisprime(n)) return gerepilecopy(av,
1551 0 : fl==1 ? mkvec2( pol_x(varn(pol)), pol)
1552 14 : : mkvec2( mkvec2(pol_x(vp),pol_0(vp)), mkvec2(pol,pol_x(vp))));
1553 231 : if (n==1 || d==1) return gerepilecopy(av,
1554 14 : fl==1 ? mkvec(pol_x(varn(pol))): _subfield(pol_x(vp),pol_0(vp)));
1555 217 : if (d==n) return gerepilecopy(av,
1556 14 : fl==1 ? mkvec(pol): _subfield(pol,pol_x(vp)));
1557 189 : ero = subfields_get_ero(pol, nf);
1558 189 : fa = subfields_get_fa(pol, nf, fa);
1559 189 : gene = maxgen_subfields(pol, fa, subf_GENERATING);
1560 :
1561 189 : if (d)
1562 : {
1563 : /* keep only generating subfields of degree a multiple of d */
1564 14 : j=1;
1565 147 : for (i=1; i<lg(gene); i++)
1566 133 : if ((lg(gel(gene,i))-1) % d == 0)
1567 : {
1568 98 : gel(gene,j) = gel(gene,i);
1569 98 : j++;
1570 : }
1571 14 : setlg(gene,j);
1572 : }
1573 189 : r = lg(gene)-1;
1574 :
1575 189 : res = zerovec(10);
1576 189 : nres = 0;
1577 189 : H = zerovec(10);
1578 189 : gel(H,1) = mkvec3(matid(n),zero_F2v(r),mkvecsmall(0));
1579 189 : len = 1;
1580 :
1581 1610 : while (len>0)
1582 : {
1583 1421 : heap_pop(&H, &len, &data);
1584 1421 : V = gel(data,1);
1585 1421 : v = gel(data,2);
1586 1421 : s = gel(data,3)[1];
1587 6153 : for (i=s+1; i<=r; i++)
1588 4732 : if (!F2v_coeff(v,i))
1589 : {
1590 3675 : W = vec_Q_primpart(intersect(V, gel(gene,i)));
1591 3675 : w = F2v_copy(v);
1592 3675 : F2v_set(w, i);
1593 3675 : newfield = 1;
1594 18130 : for (j=1; j<=r; j++)
1595 16800 : if (!F2v_coeff(w,j) && field_is_contained(W,gel(gene,j),1))
1596 : {
1597 3416 : if (j<i) { newfield = 0; break; }
1598 1071 : F2v_set(w,j);
1599 : }
1600 3675 : if (newfield && (!d || (lg(W)-1)%d==0)) heap_push(&H, &len, mkvec3(W,w,mkvecsmall(i)));
1601 : }
1602 :
1603 1421 : if (!d || lg(V)-1==d)
1604 : {
1605 1407 : nres++;
1606 1407 : if (nres == lg(res))
1607 : {
1608 28 : res2 = zerovec(2*lg(res));
1609 392 : for(j=1; j<lg(res); j++) gel(res2,j) = gel(res,j);
1610 28 : res = res2;
1611 : }
1612 1407 : gel(res,nres) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1613 : }
1614 : }
1615 189 : setlg(res,nres+1);
1616 189 : vecreverse_inplace(res);
1617 :
1618 189 : delete_var();
1619 189 : return gerepilecopy(av, res);
1620 : }
|
