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 28358 : calc_block(blockdata *B, GEN Z, GEN Y, GEN SB)
69 : {
70 28358 : 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 28358 : pari_sp av0 = avma;
73 :
74 28358 : 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 28358 : lnon = minss(BIL, r);
84 28358 : e = new_chunk(BIL);
85 28358 : n = new_chunk(r);
86 28358 : non = new_chunk(lnon);
87 28358 : pnon = new_chunk(lnon);
88 28358 : pn = new_chunk(lnon);
89 :
90 28358 : Zp = cgetg(lnon,t_VEC);
91 28358 : Zpp = cgetg(lnon,t_VEC); nn = 0;
92 58420 : for (i=1; i<r; i++) { n[i] = lg(gel(Z,i))-1; nn += n[i]; }
93 28358 : lY = lg(Y); Yp = cgetg(lY+1,t_VEC);
94 30267 : for (j=1; j<lY; j++) gel(Yp,j) = gel(Y,j);
95 :
96 : {
97 28358 : pari_sp av = avma;
98 28358 : long k = nn / B->size;
99 57715 : for (j = 1; j < r; j++)
100 29635 : if (n[j] % k) break;
101 28358 : if (j == r)
102 : {
103 28080 : gel(Yp,lY) = Z;
104 28080 : SB = print_block_system(B, Yp, SB);
105 28080 : set_avma(av);
106 : }
107 : }
108 28358 : gel(Yp,lY) = Zp;
109 :
110 28358 : K = divisorsu(n[1]); lK = lg(K);
111 120524 : for (i=1; i<lK; i++)
112 : {
113 92166 : long ngcd = n[1], k = K[i], dk = B->size*k, lpn = 0;
114 98269 : for (j=2; j<r; j++)
115 6103 : if (n[j]%k == 0)
116 : {
117 6063 : if (++lpn >= BIL) pari_err_OVERFLOW("calc_block");
118 6063 : pn[lpn] = n[j]; pnon[lpn] = j;
119 6063 : ngcd = ugcd(ngcd, n[j]);
120 : }
121 92166 : if (dk % ngcd) continue;
122 57051 : T = 1L<<lpn;
123 57051 : if (lpn == r-2)
124 : {
125 57011 : T--; /* done already above --> print_block_system */
126 57011 : if (!T) continue;
127 : }
128 :
129 3100 : if (dk == n[1])
130 : { /* empty subset, t = 0. Split out for clarity */
131 1402 : Zp[1] = Z[1]; setlg(Zp, 2);
132 2973 : for (u=1,j=2; j<r; j++) Zpp[u++] = Z[j];
133 1402 : setlg(Zpp, u);
134 1402 : SB = calc_block(B, Zpp, Yp, SB);
135 : }
136 :
137 3860 : for (t = 1; t < T; t++)
138 : { /* loop through all nonempty subsets of [1..lpn] */
139 2370 : for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
140 : {
141 1610 : if (tp&1) { nn += pn[u]; e[u] = 1; } else e[u] = 0;
142 : }
143 760 : if (dk != nn) continue;
144 :
145 1267 : for (j=1; j<r; j++) non[j]=0;
146 313 : Zp[1] = Z[1];
147 954 : for (u=2,j=1; j<=lpn; j++)
148 641 : if (e[j]) { Zp[u] = Z[pnon[j]]; non[pnon[j]] = 1; u++; }
149 313 : setlg(Zp, u);
150 954 : for (u=1,j=2; j<r; j++)
151 641 : if (!non[j]) Zpp[u++] = Z[j];
152 313 : setlg(Zpp, u);
153 313 : SB = calc_block(B, Zpp, Yp, SB);
154 : }
155 : }
156 28358 : return gc_const(av0, SB);
157 : }
158 :
159 : /* product of permutations. Put the result in perm1. */
160 : static void
161 167352 : perm_mul_i(GEN perm1, GEN perm2)
162 : {
163 167352 : long i, N = lg(perm1);
164 167352 : pari_sp av = avma;
165 167352 : GEN perm = new_chunk(N);
166 7504680 : for (i=1; i<N; i++) perm[i] = perm1[perm2[i]];
167 7504680 : for (i=1; i<N; i++) perm1[i]= perm[i];
168 167352 : set_avma(av);
169 167352 : }
170 :
171 : /* cy is a cycle; compute cy^l as a permutation */
172 : static GEN
173 35399 : cycle_power_to_perm(GEN perm,GEN cy,long l)
174 : {
175 35399 : long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
176 :
177 35399 : lp = l % lcy;
178 1437485 : for (i=1; i<N; i++) perm[i] = i;
179 35399 : if (lp)
180 : {
181 31387 : pari_sp av = avma;
182 31387 : GEN p1 = new_chunk(N);
183 31387 : b = cy[1];
184 333550 : for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
185 31387 : perm[b] = cy[1];
186 1300383 : for (i=1; i<N; i++) p1[i] = perm[i];
187 :
188 163340 : for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
189 31387 : set_avma(av);
190 : }
191 35399 : return perm;
192 : }
193 :
194 : /* image du block system D par la permutation perm */
195 : static GEN
196 17297 : im_block_by_perm(GEN D,GEN perm)
197 : {
198 17297 : long i, lb = lg(D);
199 17297 : GEN Dn = cgetg(lb,t_VEC);
200 175514 : for (i=1; i<lb; i++) gel(Dn,i) = vecsmallpermute(perm, gel(D,i));
201 17297 : return Dn;
202 : }
203 :
204 : static void
205 29974 : append(GEN D, GEN a)
206 : {
207 29974 : long i,l = lg(D), m = lg(a);
208 29974 : GEN x = D + (l-1);
209 103235 : for (i=1; i<m; i++) gel(x,i) = gel(a,i);
210 29974 : setlg(D, l+m-1);
211 29974 : }
212 :
213 : static GEN
214 28080 : print_block_system(blockdata *B, GEN Y, GEN SB)
215 : {
216 28080 : long i, j, l, ll, lp, u, v, ns, r = lg(Y), N = B->N;
217 : long *k, *n, **e, *t;
218 28080 : GEN D, De, Z, cyperm, perm, VOID = cgetg(1, t_VECSMALL);
219 :
220 28080 : if (DEBUGLEVEL>5) err_printf("Y = %Ps\n",Y);
221 28080 : n = new_chunk(N+1);
222 28080 : D = vectrunc_init(N+1);
223 28080 : t = new_chunk(r+1);
224 28080 : k = new_chunk(r+1);
225 28080 : Z = cgetg(r+1, t_VEC);
226 58054 : for (ns=0,i=1; i<r; i++)
227 : {
228 29974 : GEN Yi = gel(Y,i);
229 29974 : long ki = 0, si = lg(Yi)-1;
230 :
231 61553 : for (j=1; j<=si; j++) { n[j] = lg(gel(Yi,j))-1; ki += n[j]; }
232 29974 : ki /= B->size;
233 29974 : De = cgetg(ki+1,t_VEC);
234 103235 : for (j=1; j<=ki; j++) gel(De,j) = VOID;
235 61553 : for (j=1; j<=si; j++)
236 : {
237 31579 : GEN cy = gel(Yi,j);
238 182261 : for (l=1,lp=0; l<=n[j]; l++)
239 : {
240 150682 : lp++; if (lp > ki) lp = 1;
241 150682 : gel(De,lp) = vecsmall_append(gel(De,lp), cy[l]);
242 : }
243 : }
244 29974 : append(D, De);
245 29974 : if (si>1 && ki>1)
246 : {
247 1573 : GEN p1 = cgetg(si,t_VEC);
248 3178 : for (j=2; j<=si; j++) p1[j-1] = Yi[j];
249 1573 : ns++;
250 1573 : t[ns] = si-1;
251 1573 : k[ns] = ki-1;
252 1573 : gel(Z,ns) = p1;
253 : }
254 : }
255 28080 : if (DEBUGLEVEL>2) err_printf("\nns = %ld\n",ns);
256 28080 : if (!ns) return test_block(B, SB, D);
257 :
258 1558 : setlg(Z, ns+1);
259 1558 : e = (long**)new_chunk(ns+1);
260 3131 : for (i=1; i<=ns; i++)
261 : {
262 1573 : e[i] = new_chunk(t[i]+1);
263 3178 : for (j=1; j<=t[i]; j++) e[i][j] = 0;
264 : }
265 1558 : cyperm= cgetg(N+1,t_VECSMALL);
266 1558 : perm = cgetg(N+1,t_VECSMALL); i = ns;
267 : do
268 : {
269 17297 : pari_sp av = avma;
270 560909 : for (u=1; u<=N; u++) perm[u] = u;
271 35134 : for (u=1; u<=ns; u++)
272 53236 : for (v=1; v<=t[u]; v++)
273 35399 : perm_mul_i(perm, cycle_power_to_perm(cyperm, gmael(Z,u,v), e[u][v]));
274 17297 : SB = test_block(B, SB, im_block_by_perm(D,perm));
275 17297 : 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 17297 : if (++e[ns][t[ns]] > k[ns])
280 : {
281 2404 : j = t[ns]-1;
282 2491 : while (j>=1 && e[ns][j] == k[ns]) j--;
283 3230 : if (j >= 1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l] = 0; }
284 : else
285 : {
286 1633 : i = ns-1;
287 1648 : while (i>=1)
288 : {
289 90 : j = t[i];
290 105 : while (j>=1 && e[i][j] == k[i]) j--;
291 90 : if (j<1) i--;
292 : else
293 : {
294 75 : e[i][j]++;
295 75 : for (l=j+1; l<=t[i]; l++) e[i][l] = 0;
296 150 : for (ll=i+1; ll<=ns; ll++)
297 150 : for (l=1; l<=t[ll]; l++) e[ll][l] = 0;
298 75 : break;
299 : }
300 : }
301 : }
302 : }
303 : }
304 17297 : while (i > 0);
305 1558 : return SB;
306 : }
307 :
308 : /* ALGEBRAIC PART: test potential block systems */
309 :
310 : static GEN
311 33874 : polsimplify(GEN x)
312 : {
313 33874 : long i,lx = lg(x);
314 194844 : for (i=2; i<lx; i++)
315 160970 : if (typ(gel(x,i)) == t_POL) gel(x,i) = constant_coeff(gel(x,i));
316 33874 : return x;
317 : }
318 :
319 : /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
320 : static long
321 33874 : ok_coeffs(GEN g,GEN M)
322 : {
323 33874 : long i, lg = lg(g)-1; /* g is monic, and cst term is ok */
324 82832 : for (i=3; i<lg; i++)
325 54441 : if (abscmpii(gel(g,i), gel(M,i)) > 0) return 0;
326 28391 : return 1;
327 : }
328 :
329 : /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
330 : static GEN
331 147404 : trace(GEN x, GEN Trq, GEN p)
332 : {
333 : long i, l;
334 : GEN s;
335 147404 : if (typ(x) == t_INT) return Fp_mul(x, gel(Trq,1), p);
336 147404 : l = lg(x)-1; if (l == 1) return gen_0;
337 147404 : x++; s = mulii(gel(x,1), gel(Trq,1));
338 738797 : for (i=2; i<l; i++)
339 591393 : s = addii(s, mulii(gel(x,i), gel(Trq,i)));
340 147404 : 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 31142 : poltrace(GEN x, GEN Trq, GEN p)
346 : {
347 31142 : if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
348 178546 : 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 27760 : chinese_retrieve_pol(GEN DATA, primedata *S, GEN listdelta)
358 : {
359 27760 : GEN interp, bezoutC, h, p = S->p, pol = FpX_red(gel(DATA,1), p);
360 : long i, l;
361 27760 : interp = gel(DATA,9);
362 27760 : bezoutC= gel(DATA,6);
363 :
364 27760 : h = NULL; l = lg(interp);
365 58902 : for (i=1; i<l; i++)
366 : { /* h(firstroot[i]) = listdelta[i] */
367 31142 : GEN t = FqX_Fq_mul(gel(interp,i), gel(listdelta,i), S->T, p);
368 31142 : t = poltrace(t, gel(S->Trk,i), p);
369 31142 : t = FpX_mul(t, gel(bezoutC,i), p);
370 31142 : h = h? FpX_add(h,t,p): t;
371 : }
372 27760 : 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 27760 : embedding(GEN g, GEN DATA, primedata *S, GEN den, GEN listdelta)
382 : {
383 27760 : 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 27760 : T = gel(DATA,1); rt = brent_kung_optpow(degpol(T), 4, 3);
388 27760 : maxp= gel(DATA,7);
389 27760 : gp = RgX_deriv(g); av = avma;
390 27760 : w0 = chinese_retrieve_pol(DATA, S, listdelta);
391 27760 : w0_Q = centermod(gmul(w0,den), p);
392 27760 : h0 = FpXQ_inv(FpX_FpXQ_eval(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
393 27760 : 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 98878 : 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 98878 : if (wpow) a = FpX_FpXQV_eval(g,wpow, T,q);
403 27760 : else a = FpX_FpXQ_eval(g,w0, T,q); /* first time */
404 : /* now, a = 0 (p) */
405 98878 : a = FpXQ_mul(ZX_neg(h0), ZX_Z_divexact(a, p), T,p);
406 98878 : w1 = ZX_add(w0, ZX_Z_mul(a, p));
407 :
408 98878 : w1_Q = centermod(ZX_Z_mul(w1, remii(den,q)), q);
409 98878 : if (ZX_equal(w1_Q, w0_Q))
410 : {
411 18531 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
412 18531 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
413 : }
414 80347 : else if (cmpii(q,maxp) > 0)
415 : {
416 12784 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
417 12784 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
418 426 : if (DEBUGLEVEL) err_printf("coeff too big for embedding\n");
419 426 : return NULL;
420 : }
421 71118 : (void)gc_all(av, 5, &w1,&h0,&w1_Q,&q,&p);
422 71118 : q2 = sqri(q);
423 71118 : 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 71118 : a = FpXQ_mul(ZX_neg(h0), FpX_FpXQV_eval(gp, FpXV_red(wpow,q),T,q), T,q);
427 71118 : a = ZX_Z_add_shallow(a, gen_1); /* 1 - h0 g'(w1) = 0 (p) */
428 71118 : a = FpXQ_mul(h0, ZX_Z_divexact(a, p), T,p);
429 71118 : h0 = ZX_add(h0, ZX_Z_mul(a, p));
430 71118 : w0 = w1; w0_Q = w1_Q; p = q; q = q2;
431 : }
432 27334 : TR = gel(DATA,5);
433 27334 : if (!gequal0(TR)) w1_Q = RgX_Rg_translate(w1_Q, TR);
434 27334 : 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 26397 : get_bezout(GEN pol, GEN fk, GEN p)
441 : {
442 26397 : long i, l = lg(fk);
443 26397 : GEN A, B, d, u, v, U = cgetg(l, t_VEC);
444 54184 : for (i=1; i<l; i++)
445 : {
446 27787 : A = gel(fk,i);
447 27787 : B = FpX_div(pol, A, p);
448 27787 : d = FpX_extgcd(A,B,p, &u, &v);
449 27787 : if (degpol(d) > 0) pari_err_COPRIME("get_bezout",A,B);
450 27787 : d = constant_coeff(d);
451 27787 : if (!gequal1(d)) v = FpX_Fp_div(v, d, p);
452 27787 : gel(U,i) = FpX_mul(B,v, p);
453 : }
454 26397 : return U;
455 : }
456 :
457 : static GEN
458 26397 : init_traces(GEN ff, GEN T, GEN p)
459 : {
460 26397 : long N = degpol(T),i,j,k, r = lg(ff);
461 26397 : GEN Frob = FpX_matFrobenius(T,p);
462 : GEN y,p1,p2,Trk,pow,pow1;
463 :
464 26397 : k = degpol(gel(ff,r-1)); /* largest degree in modular factorization */
465 26397 : pow = cgetg(k+1, t_VEC);
466 26397 : gel(pow,1) = gen_0; /* dummy */
467 26397 : gel(pow,2) = Frob;
468 26397 : pow1= cgetg(k+1, t_VEC); /* 1st line */
469 94457 : for (i=3; i<=k; i++)
470 68060 : gel(pow,i) = FpM_mul(gel(pow,i-1), Frob, p);
471 26397 : gel(pow1,1) = gen_0; /* dummy */
472 120854 : for (i=2; i<=k; i++)
473 : {
474 94457 : p1 = cgetg(N+1, t_VEC);
475 94457 : gel(pow1,i) = p1; p2 = gel(pow,i);
476 657157 : 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 26397 : Trk = pow; /* re-use (destroy) pow */
481 26397 : gel(Trk,1) = vec_ei(N,1);
482 120854 : for (i=2; i<=k; i++)
483 94457 : gel(Trk,i) = gadd(gel(Trk,i-1), gel(pow1,i));
484 26397 : y = cgetg(r, t_VEC);
485 54184 : for (i=1; i<r; i++) y[i] = Trk[degpol(gel(ff,i))];
486 26397 : return y;
487 : }
488 :
489 : static void
490 26397 : init_primedata(primedata *S)
491 : {
492 26397 : long i, j, l, lff = lg(S->ff), v = fetch_var(), N = degpol(S->pol);
493 26397 : GEN T, p = S->p;
494 :
495 26397 : if (S->lcm == degpol(gel(S->ff,lff-1)))
496 : {
497 26387 : T = leafcopy(gel(S->ff,lff-1));
498 26387 : setvarn(T, v);
499 : }
500 : else
501 10 : T = init_Fq(p, S->lcm, v);
502 26397 : S->T = T;
503 26397 : S->firstroot = cgetg(lff, t_VECSMALL);
504 26397 : S->interp = cgetg(lff, t_VEC);
505 26397 : S->fk = cgetg(N+1, t_VEC);
506 54184 : for (l=1,j=1; j<lff; j++)
507 : { /* compute roots and fix ordering (Frobenius cycles) */
508 27787 : GEN F = gel(S->ff, j), deg1 = FpX_factorff_irred(F, T,p);
509 27787 : GEN H = gel(deg1,1), a = Fq_neg(constant_coeff(H), T,p);
510 27787 : GEN Q = FqX_div(F, H, T,p);
511 27787 : GEN q = Fq_inv(FqX_eval(Q, a, T,p), T,p);
512 27787 : gel(S->interp,j) = FqX_Fq_mul(Q, q, T,p); /* = 1 at a, 0 at other roots */
513 27787 : S->firstroot[j] = l;
514 156187 : for (i=1; i<lg(deg1); i++,l++) gel(S->fk, l) = gel(deg1, i);
515 : }
516 26397 : S->Trk = init_traces(S->ff, T,p);
517 26397 : S->bezoutC = get_bezout(S->pol, S->ff, p);
518 26397 : }
519 :
520 : static int
521 26407 : choose_prime(primedata *S, GEN pol)
522 : {
523 26407 : 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 26407 : u_forprime_init(&T, (N*N) >> 2, ULONG_MAX);
530 26407 : oldr = S->lcm = LONG_MAX;
531 26407 : S->ff = oldn = NULL; pp = 0; /* gcc -Wall */
532 26407 : av = avma; K = N + 10;
533 84886 : for(k = 1; k < K || !S->ff; k++,set_avma(av))
534 : {
535 : GEN Tp;
536 83662 : if (k > 5 * N) return 0;
537 : do
538 : {
539 99328 : p = u_forprime_next(&T);
540 99328 : Tp = ZX_to_Flx(pol, p);
541 : }
542 99328 : while (!Flx_is_squarefree(Tp, p));
543 83662 : ff = gel(Flx_factor(Tp, p), 1);
544 83662 : r = lg(ff)-1;
545 83662 : if (r == N || r >= BIL) continue;
546 :
547 81579 : n = cgetg(r+1, t_VECSMALL); lcm = n[1] = degpol(gel(ff,1));
548 216028 : for (j=2; j<=r; j++) { n[j] = degpol(gel(ff,j)); lcm = ulcm(lcm, n[j]); }
549 81579 : if (r > oldr || (r == oldr && (lcm <= S->lcm || S->lcm > 2*N)))
550 38475 : continue;
551 43104 : if (DEBUGLEVEL) err_printf("p = %lu,\tlcm = %ld,\torbits: %Ps\n",p,lcm,n);
552 :
553 43104 : pp = p;
554 43104 : oldr = r;
555 43104 : oldn = n;
556 43104 : S->ff = ff;
557 43104 : S->lcm = lcm; if (r == 1) break;
558 17921 : av = avma;
559 : }
560 26407 : if (oldr > 6) return 0;
561 26397 : if (DEBUGLEVEL) err_printf("Chosen prime: p = %ld\n", pp);
562 26397 : FlxV_to_ZXV_inplace(S->ff);
563 26397 : S->p = utoipos(pp);
564 26397 : S->pol = FpX_red(pol, S->p); init_primedata(S);
565 26397 : n = oldn; r = lg(n); S->Z = Z = cgetg(r,t_VEC);
566 54184 : for (k=0,i=1; i<r; i++)
567 : {
568 27787 : GEN t = cgetg(n[i]+1, t_VECSMALL); gel(Z,i) = t;
569 156187 : for (j=1; j<=n[i]; j++) t[j] = ++k;
570 : }
571 26397 : return 1;
572 : }
573 :
574 : /* maxroot t_REAL */
575 : static GEN
576 27274 : bound_for_coeff(long m, GEN R, GEN *maxroot)
577 : {
578 27274 : GEN b1, b2, M, v, C = vecbinomial(m-1);
579 27274 : long i, r1, l = lg(R);
580 :
581 41280 : for (r1 = 1; r1 < l; r1++)
582 40038 : if (typ(gel(R,r1)) != t_REAL) break;
583 27274 : r1--;
584 27274 : R = gabs(R,0); *maxroot = vecmax(R);
585 41280 : for (b1 = gen_1, i = 1; i <= r1; i++)
586 14006 : if (gcmpgs(gel(R,i), 1) > 0) b1 = gmul(b1, gel(R,i));
587 152530 : for (b2 = gen_1 ; i < l; i++)
588 125256 : if (gcmpgs(gel(R,i), 1) > 0) b2 = gmul(b2, gel(R,i));
589 27274 : M = gmul(b1, gsqr(b2)); /* Mahler measure */
590 27274 : v = cgetg(m+2, t_VEC); gel(v,1) = gel(v,2) = gen_0; /* unused */
591 67956 : for (i = 1; i < m; i++) /* binom(m-1, i) * M + binom(m-1, i-1) */
592 40682 : gel(v, i+2) = ceil_safe(gadd(gmul(gel(C, i+1), M), gel(C, i)));
593 27274 : return v;
594 : }
595 :
596 : static GEN
597 7893 : RgV_translate(GEN x, GEN t) { pari_APPLY_same(gadd(t, gel(x,i))); }
598 : static GEN
599 7893 : 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 27274 : compute_data(blockdata *B)
615 : {
616 : GEN ffL, roo, pe, p1, p2, fk, MM, maxroot, pol;
617 27274 : primedata *S = B->S;
618 27274 : GEN p = S->p, T = S->T, ff = S->ff, DATA = B->DATA;
619 27274 : long i, l, e, N, lff = lg(ff);
620 :
621 27274 : if (DEBUGLEVEL>1) err_printf("Entering compute_data()\n\n");
622 27274 : pol = B->PD->pol; N = degpol(pol);
623 27274 : roo = B->PD->roo;
624 27274 : if (DATA)
625 : {
626 631 : GEN TR = addiu(gel(DATA,5), 1), mTR = negi(TR), interp, bezoutC;
627 :
628 631 : if (DEBUGLEVEL>1) err_printf("... update (translate) an existing DATA\n\n");
629 631 : gel(DATA,5) = TR;
630 631 : pol = RgX_Rg_translate(gel(DATA,1), gen_m1);
631 631 : roo = RgV_translate(roo, TR);
632 631 : fk = RgV_negtranslate(gel(DATA,4),
633 631 : deg1pol_shallow(gen_1, gen_m1, varn(pol)));
634 631 : bezoutC = gel(DATA,6); l = lg(bezoutC);
635 631 : interp = gel(DATA,9);
636 1836 : for (i=1; i<l; i++)
637 : {
638 1205 : if (degpol(gel(interp,i)) > 0) /* do not turn pol_1(0) into gen_1 */
639 1205 : gel(interp,i) = FpXX_red(RgX_Rg_translate(gel(interp,i), gen_m1), p);
640 1205 : if (degpol(gel(bezoutC,i)) > 0)
641 1132 : gel(bezoutC,i)= FpXX_red(RgX_Rg_translate(gel(bezoutC,i), gen_m1), p);
642 : }
643 631 : ff = cgetg(lff, t_VEC); /* copy, do not overwrite! */
644 1836 : for (i=1; i<lff; i++)
645 1205 : gel(ff,i) = FpX_red(RgX_Rg_translate(gel(S->ff,i), mTR), p);
646 : }
647 : else
648 : {
649 26643 : DATA = cgetg(10,t_VEC);
650 26643 : fk = S->fk;
651 26643 : gel(DATA,5) = gen_0;
652 26643 : gel(DATA,6) = leafcopy(S->bezoutC);
653 26643 : gel(DATA,9) = leafcopy(S->interp);
654 : }
655 27274 : gel(DATA,1) = pol;
656 27274 : MM = gmul2n(bound_for_coeff(B->d, roo, &maxroot), 1);
657 27274 : gel(DATA,8) = MM;
658 27274 : e = logintall(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [d-1 test] */
659 27274 : gel(DATA,2) = pe;
660 27274 : 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 27274 : ffL = ZpX_liftfact(pol, ff, pe, p, e);
665 56642 : for (l=i=1; i<lff; i++)
666 : { /* 2) lift factorization of ff[i] in Qp[X] / T */
667 29368 : long l2 = l + degpol(gel(ffL,i));
668 29368 : gel(ffL,i) = ZqX_liftfact(gel(ffL,i), vecslice(fk, l, l2-1), T, pe, p, e);
669 29368 : l = l2;
670 : }
671 27274 : gel(DATA,3) = roots_from_deg1(shallowconcat1(ffL));
672 :
673 27274 : p1 = mulur(N, powruhalf(utor(N-1,DEFAULTPREC), N-1));
674 27274 : p2 = powru(maxroot, B->size + N*(N-1)/2);
675 27274 : p1 = divrr(mulrr(p1,p2), gsqrt(B->PD->dis,DEFAULTPREC));
676 27274 : gel(DATA,7) = mulii(shifti(ceil_safe(p1), 1), B->PD->den);
677 :
678 27274 : 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 27274 : if (B->DATA) { DATA = gclone(DATA); if (isclone(B->DATA)) gunclone(B->DATA); }
685 27274 : B->DATA = DATA;
686 27274 : }
687 :
688 : /* g = polynomial, h = embedding. Return [[g,h]] */
689 : static GEN
690 1164 : _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 44450 : subfield(GEN A, blockdata *B)
695 : {
696 44450 : 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 44450 : GEN T = B->S->T, p = B->S->p, firstroot = B->S->firstroot;
699 :
700 44450 : pol= gel(B->DATA,1); N = degpol(pol); d = N/m; /* m | N */
701 44450 : pe = gel(B->DATA,2);
702 44450 : fhk= gel(B->DATA,3);
703 44450 : M = gel(B->DATA,8);
704 :
705 44450 : delta = cgetg(m+1,t_VEC);
706 44450 : whichdelta = cgetg(N+1, t_VECSMALL);
707 44450 : d_1_term = gen_0;
708 270083 : for (i=1; i<=m; i++)
709 : {
710 225633 : GEN Ai = gel(A,i), p1 = gel(fhk,Ai[1]);
711 681284 : for (j=2; j<=d; j++)
712 455651 : p1 = Fq_mul(p1, gel(fhk,Ai[j]), T, pe);
713 225633 : gel(delta,i) = p1;
714 225633 : 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 906917 : for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
719 225633 : if (typ(p1) == t_POL) p1 = constant_coeff(p1);
720 225633 : d_1_term = addii(d_1_term, p1);
721 : }
722 44450 : d_1_term = centermod(d_1_term, pe); /* Tr(g) */
723 44450 : if (abscmpii(d_1_term, gel(M,3)) > 0) {
724 10576 : if (DEBUGLEVEL>1) err_printf("d-1 test failed\n");
725 10576 : return NULL;
726 : }
727 33874 : g = FqV_roots_to_pol(delta, T, pe, 0);
728 33874 : g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
729 33874 : if (!ok_coeffs(g,M)) {
730 5483 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
731 5483 : if (DEBUGLEVEL>1) err_printf("coeff too big for pol g(x)\n");
732 5483 : return NULL;
733 : }
734 28391 : if (!FpX_is_squarefree(g, p)) {
735 631 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
736 631 : if (DEBUGLEVEL>1) err_printf("changing f(x): p divides disc(g)\n");
737 631 : compute_data(B);
738 631 : return subfield(A, B);
739 : }
740 :
741 27760 : lf = lg(firstroot); listdelta = cgetg(lf, t_VEC);
742 58902 : for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
743 27760 : if (DEBUGLEVEL) err_printf("candidate = %Ps\n", g);
744 27760 : e = embedding(g, B->DATA, B->S, B->PD->den, listdelta);
745 27760 : if (!e) return NULL;
746 27334 : if (DEBUGLEVEL) err_printf("... OK!\n");
747 27334 : 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 43819 : test_block(blockdata *B, GEN L, GEN D)
754 : {
755 43819 : pari_sp av = avma;
756 43819 : GEN sub = subfield(D, B);
757 43819 : if (sub) {
758 27334 : GEN old = L;
759 27334 : L = gclone( L? shallowconcat(L, sub): sub );
760 27334 : guncloneNULL(old);
761 : }
762 43819 : return gc_const(av,L);
763 : }
764 :
765 : /* subfields of degree d */
766 : static GEN
767 26643 : subfields_of_given_degree(blockdata *B)
768 : {
769 26643 : pari_sp av = avma;
770 : GEN L;
771 :
772 26643 : if (DEBUGLEVEL) err_printf("\n* Look for subfields of degree %ld\n\n", B->d);
773 26643 : B->DATA = NULL; compute_data(B);
774 26643 : L = calc_block(B, B->S->Z, cgetg(1,t_VEC), NULL);
775 26643 : if (DEBUGLEVEL>9)
776 0 : err_printf("\nSubfields of degree %ld: %Ps\n", B->d, L? L: cgetg(1,t_VEC));
777 26643 : if (isclone(B->DATA)) gunclone(B->DATA);
778 26643 : return gc_const(av,L);
779 : }
780 :
781 : static void
782 30 : setvarn2(GEN t, long v) { setvarn(gel(t,1),v); setvarn(gel(t,2),v); }
783 : static GEN
784 25524 : fix_var(GEN x, long v, long fl)
785 : {
786 25524 : long i, l = lg(x);
787 25524 : if (!v) return x;
788 22 : if (fl)
789 30 : for (i = 1; i < l; i++) setvarn(gel(x,i), v);
790 : else
791 42 : for (i = 1; i < l; i++) setvarn2(gel(x,i), v);
792 22 : return x;
793 : }
794 :
795 : static void
796 26397 : subfields_poldata(GEN nf, GEN T, poldata *PD)
797 : {
798 : GEN L, dis;
799 :
800 26397 : PD->pol = T;
801 26397 : if (nf)
802 : {
803 144 : PD->den = nf_get_zkden(nf);
804 144 : PD->roo = nf_get_roots(nf);
805 144 : PD->dis = mulii(absi_shallow(nf_get_disc(nf)), sqri(nf_get_index(nf)));
806 : }
807 : else
808 : {
809 26253 : PD->den = initgaloisborne(T,NULL,nbits2prec(bit_accuracy(ZX_max_lg(T))), &L,NULL,&dis);
810 26253 : PD->roo = L;
811 26253 : PD->dis = absi_shallow(dis);
812 : }
813 26397 : }
814 :
815 : static GEN nfsubfields_fa(GEN nf, long d, long fl);
816 : static GEN
817 473 : subfieldsall(GEN nf0, long fl)
818 : {
819 473 : 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 473 : G = galoisinit(nf0, NULL);
828 473 : T = get_nfpol(nf0, &nf);
829 473 : if (G != gen_0)
830 : {
831 : GEN L, S;
832 : long l;
833 40 : L = lift_shallow( galoissubfields(G, fl, varn(T)) );
834 40 : l = lg(L); S = cgetg(l, t_VECSMALL);
835 660 : for (i=1; i<l; i++) S[i] = lg(fl==1? gel(L,i): gmael(L,i,1));
836 40 : return gc_GEN(av, vecpermute(L, vecsmall_indexsort(S)));
837 : }
838 433 : v = varn(T); N = degpol(T);
839 433 : dg = divisorsu(N); ld = lg(dg)-1;
840 433 : LSB = fl==1 ? mkvec(pol_x(v)): _subfield(pol_x(v), pol_0(v));
841 433 : if (ld <= 2)
842 : {
843 144 : if (ld == 2)
844 144 : LSB = shallowconcat(LSB, fl==1? mkvec(T): _subfield(T, pol_x(v)));
845 144 : return gc_GEN(av, LSB);
846 : }
847 289 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
848 289 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, 0, fl); }
849 279 : subfields_poldata(nf, T, &PD);
850 :
851 279 : if (DEBUGLEVEL) err_printf("\n***** Entering subfields\n\npol = %Ps\n",T);
852 279 : B.PD = &PD;
853 279 : B.S = &S;
854 279 : B.N = N;
855 279 : B.fl = fl;
856 804 : for (i=ld-1; i>1; i--)
857 : {
858 525 : B.size = dg[i];
859 525 : B.d = N / B.size;
860 525 : NLSB = subfields_of_given_degree(&B);
861 525 : if (NLSB) { LSB = gconcat(LSB, NLSB); gunclone(NLSB); }
862 : }
863 279 : (void)delete_var(); /* from init_primedata() */
864 279 : LSB = shallowconcat(LSB, fl==1? mkvec(T):_subfield(T, pol_x(0)));
865 279 : if (DEBUGLEVEL) err_printf("\n***** Leaving subfields\n\n");
866 279 : return fix_var(gc_GEN(av, LSB), v, fl);
867 : }
868 :
869 : GEN
870 27970 : nfsubfields0(GEN nf0, long d, long fl)
871 : {
872 27970 : pari_sp av = avma;
873 : long N, v0;
874 : GEN nf, LSB, T, G;
875 : poldata PD;
876 : primedata S;
877 : blockdata B;
878 27970 : if (fl<0 || fl>1) pari_err_FLAG("nfsubfields");
879 27970 : if (typ(nf0)==t_VEC && lg(nf0)==3) return nfsubfields_fa(nf0, d, fl);
880 27795 : if (!d) return subfieldsall(nf0, fl);
881 :
882 : /* treat trivial cases */
883 27322 : T = get_nfpol(nf0, &nf); v0 = varn(T); N = degpol(T);
884 27322 : RgX_check_ZX(T,"nfsubfields");
885 27322 : if (d == N)
886 20 : return gc_GEN(av, fl==1 ? mkvec(T) : _subfield(T, pol_x(v0)));
887 27302 : if (d == 1)
888 20 : return gc_GEN(av, fl==1 ? mkvec(pol_x(v0)) : _subfield(pol_x(v0), zeropol(v0)));
889 27282 : if (d < 1 || d > N || N % d) return cgetg(1,t_VEC);
890 :
891 : /* much easier if nf is Galois (WSS) */
892 27252 : G = galoisinit(nf0, NULL);
893 27252 : if (G != gen_0)
894 : { /* Bingo */
895 1134 : GEN L = galoissubgroups(G), F;
896 1134 : long k,i, l = lg(L), o = N/d;
897 1134 : F = cgetg(l, t_VEC);
898 1134 : k = 1;
899 5008 : for (i=1; i<l; i++)
900 : {
901 3874 : GEN H = gel(L,i);
902 3874 : if (group_order(H) == o)
903 1305 : gel(F,k++) = lift_shallow(galoisfixedfield(G, gel(H,1), fl, v0));
904 : }
905 1134 : setlg(F, k);
906 1134 : return gc_GEN(av, F);
907 : }
908 26118 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
909 26118 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, d, fl); }
910 26118 : subfields_poldata(nf, T, &PD);
911 26118 : B.PD = &PD;
912 26118 : B.S = &S;
913 26118 : B.N = N;
914 26118 : B.d = d;
915 26118 : B.size = N/d;
916 26118 : B.fl = fl;
917 26118 : LSB = subfields_of_given_degree(&B);
918 26118 : (void)delete_var(); /* from init_primedata */
919 26118 : set_avma(av);
920 26118 : if (!LSB) return cgetg(1, t_VEC);
921 25245 : G = gcopy(LSB); gunclone(LSB);
922 25245 : return fix_var(G, v0, fl);
923 : }
924 :
925 : GEN
926 282 : nfsubfields(GEN nf0, long d)
927 282 : { 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 2375 : 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 2375 : a = gtopolyrev(v, varn(pol));
944 2375 : P = Flxq_charpoly(ZX_to_Flx(a,p), ZX_to_Flx(pol,p), p);
945 2375 : Flx_ispower(P, e, p, &Q);
946 2375 : if (!Flx_is_squarefree(Q,p)) return NULL;
947 1220 : d = degpol(pol)/e;
948 1220 : B = 0;
949 19050 : for (i=1; i<lg(v); i++)
950 : {
951 17830 : bi = (i-1)*ero + expi(gel(v,i));
952 17830 : if (bi > B) B = bi;
953 : }
954 1220 : ed = expu(d);
955 1220 : B += ed+1;
956 1220 : bound = 0;
957 5913 : for (i=0; 2*i<=d; i++)
958 : {
959 4693 : if (!i) bi = d*B;
960 3473 : else bi = (d-i)*B + i*(3+ed-expu(i));
961 4693 : if (bi > bound) bound = bi;
962 : }
963 1220 : Q = ZXQ_minpoly(a,pol,d,bound);
964 1220 : 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 1345 : subfield_generator(GEN pol, GEN V, long d, long ero, long fl)
977 : {
978 1345 : long p, i, e, vp = varn(pol);
979 1345 : GEN a = NULL, v = cgetg(lg(V),t_COL), B;
980 :
981 1345 : if (d==1) return fl ? pol_x(vp): mkvec2(pol_x(vp), pol_0(vp));
982 1220 : e = degpol(pol)/d;
983 1220 : p = 1009;
984 2375 : for (i=1; i<lg(V); i++)
985 : {
986 2360 : a = try_subfield_generator(pol, gel(V,i), e, p, ero, fl);
987 2360 : if (a) return a;
988 1155 : p = unextprime(p+1);
989 : }
990 15 : B = stoi(10);
991 : while(1)
992 : {
993 65 : for (i=1; i<lg(v); i++) gel(v,i) = randomi(B);
994 15 : a = try_subfield_generator(pol, QM_QC_mul(V,v), e, p, ero, fl);
995 15 : if (a) return a;
996 0 : p = unextprime(p+1);
997 : }
998 : return NULL;/*LCOV_EXCL_LINE*/
999 : }
1000 :
1001 : static GEN
1002 27183 : RgXY_to_RgC(GEN P, long dx, long dy)
1003 : {
1004 : GEN res, c;
1005 27183 : long i, j, k, d = degpol(P);
1006 27183 : if (d > dy) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1007 27183 : res = cgetg((dx+1)*(dy+1)+1, t_COL);
1008 27183 : k = 1;
1009 67074 : for (i=0; i<=d; i++)
1010 : {
1011 39891 : c = gel(P,i+2);
1012 39891 : if (typ(c)==t_POL)
1013 : {
1014 36825 : long dc = degpol(c);
1015 36825 : if (dc > dx) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1016 730499 : for (j=0; j<=dc; j++)
1017 693674 : gel(res,k++) = gel(c,j+2);
1018 : } else
1019 : {
1020 3066 : gel(res,k++) = c; j=1;
1021 : }
1022 246580 : for ( ; j<=dx; j++)
1023 206689 : gel(res,k++) = gen_0;
1024 : }
1025 62532 : for( ; i<=dy; i++)
1026 768930 : for (j=0; j<=dx; j++)
1027 733581 : gel(res,k++) = gen_0;
1028 27183 : return res;
1029 : }
1030 :
1031 : /* lambda: t_VEC of t_INT; 0 means ignore this factor */
1032 : static GEN
1033 1855 : twoembequation(GEN pol, GEN fa, GEN lambda)
1034 : {
1035 : GEN m, vpolx, poly;
1036 1855 : long i,j, lfa = lg(fa), dx = degpol(pol);
1037 1855 : long vx = varn(pol), vy = varn(gel(fa,1)); /* vx < vy ! */
1038 :
1039 1855 : if (varncmp(vx,vy) <= 0) pari_err_BUG("twoembequation [incorrect variable priorities]");
1040 :
1041 1855 : lambda = shallowcopy(lambda);
1042 1855 : fa = shallowcopy(fa);
1043 1855 : j = 1;
1044 19931 : for (i=1; i<lfa; i++)
1045 18076 : if (signe(gel(lambda,i)))
1046 : {
1047 1881 : gel(lambda,j) = negi(gel(lambda,i));
1048 1881 : gel(fa,j) = gel(fa,i);
1049 1881 : j++;
1050 : }
1051 1855 : setlg(lambda, j);
1052 1855 : setlg(fa, j); lfa = j;
1053 :
1054 1855 : vpolx = ZXQ_powers(pol_x(vx),dx-1,pol);
1055 1855 : m = cgetg(dx+1, t_MAT);
1056 28850 : for (j=1; j <= dx; j++)
1057 26995 : gel(m,j) = cgetg(lfa, t_COL);
1058 3736 : for(i=1; i<lfa; i++)
1059 : {
1060 1881 : long dy = degpol(gel(fa,i));
1061 1881 : poly = pol_1(vy);
1062 29064 : for (j=1; j <= dx; j++)
1063 : {
1064 27183 : gcoeff(m,i,j) = RgXY_to_RgC(gadd(ZX_Z_mul(gel(vpolx,j),gel(lambda,i)),poly), dx, dy);
1065 27183 : poly = RgXQX_rem(RgX_shift(poly,1), gel(fa,i), pol);
1066 : }
1067 : }
1068 28850 : for(j=1; j<=dx; j++) gel(m,j) = shallowconcat1(gel(m,j));
1069 1855 : return QM_ker(m);
1070 : }
1071 :
1072 : static void
1073 1210 : subfields_cleanup(GEN* nf, GEN* pol, long* n, GEN* fa)
1074 : {
1075 1210 : *fa = NULL;
1076 1210 : if (typ(*nf) != t_VEC && typ(*nf) != t_POL) pari_err_TYPE("subfields_cleanup", *nf);
1077 1205 : if (typ(*nf) == t_VEC && lg(*nf) == 3)
1078 : {
1079 215 : *fa = gel(*nf,2);
1080 215 : *nf = gel(*nf,1);
1081 215 : if (typ(*fa)!=t_MAT || lg(*fa)!=3)
1082 10 : pari_err_TYPE("subfields_cleanup [fa should be a factorisation matrix]", *fa);
1083 : }
1084 1195 : if (typ(*nf) == t_POL)
1085 : {
1086 560 : *pol = *nf;
1087 560 : *nf = NULL;
1088 560 : if (!RgX_is_ZX(*pol)) pari_err_TYPE("subfields_cleanup [not integral]", *pol);
1089 555 : if (!equali1(leading_coeff(*pol))) pari_err_TYPE("subfields_cleanup [not monic]", *pol);
1090 550 : *n = degpol(*pol);
1091 550 : if (*n<=0) pari_err_TYPE("subfields_cleanup [constant polynomial]", *pol);
1092 : }
1093 : else
1094 : {
1095 635 : *nf = checknf(*nf);
1096 620 : *pol = nf_get_pol(*nf);
1097 620 : *n = degpol(*pol);
1098 : }
1099 1165 : if(*fa)
1100 : {
1101 195 : long v = varn(*pol);
1102 195 : GEN o = gcoeff(*fa,1,1);
1103 195 : if (varncmp(varn(o),v) >= 0) pari_err_PRIORITY("nfsubfields_fa", o, "<=", v);
1104 : }
1105 1150 : }
1106 :
1107 : static GEN
1108 200 : rootsuptoconj(GEN pol, long prec)
1109 : {
1110 : GEN ro;
1111 : long n, i;
1112 200 : ro = roots(pol,prec);
1113 200 : n = lg(ro)-1;
1114 1070 : for (i=1; i<=n/2; i++)
1115 870 : gel(ro,i) = gel(ro,2*i-1);
1116 200 : setlg(ro,n/2+1);
1117 200 : return ro;
1118 : }
1119 : static GEN
1120 870 : cmsubfield_get_roots(GEN pol, GEN nf, long n, long* r2, long *prec)
1121 : {
1122 : GEN ro;
1123 870 : if (nf)
1124 : {
1125 620 : if (nf_get_r1(nf)) return NULL;
1126 308 : *r2 = nf_get_r2(nf);
1127 308 : *prec = nf_get_prec(nf);
1128 308 : ro = nf_get_roots(nf);
1129 : }
1130 : else
1131 : {
1132 250 : if (n%2 || sturm(pol)) return NULL;
1133 200 : *r2 = n/2;
1134 200 : *prec = MEDDEFAULTPREC;
1135 200 : ro = rootsuptoconj(pol, *prec);
1136 : }
1137 508 : return ro;
1138 : }
1139 :
1140 : static GEN
1141 435 : subfields_get_fa(GEN pol, GEN nf, GEN fa)
1142 : {
1143 435 : if (!fa)
1144 : {
1145 290 : GEN poly = shallowcopy(pol);
1146 290 : setvarn(poly, fetch_var_higher());
1147 290 : fa = nffactor(nf? nf: pol, poly);
1148 : }
1149 435 : return liftpol_shallow(gel(fa,1));
1150 : }
1151 :
1152 : static long
1153 205 : subfields_get_ero(GEN pol, GEN nf)
1154 : {
1155 410 : return 1 + gexpo(nf? nf_get_roots(nf):
1156 205 : QX_complex_roots(pol, LOWDEFAULTPREC));
1157 : }
1158 :
1159 : static GEN
1160 224 : try_imag(GEN x, GEN c, GEN pol, long v, ulong p, GEN emb, GEN galpol, long fl)
1161 : {
1162 224 : GEN a = Q_primpart(RgX_sub(RgX_RgXQ_eval(x,c,pol),x));
1163 224 : if (Flx_is_squarefree(Flxq_charpoly(ZX_to_Flx(a,p),ZX_to_Flx(pol,p),p),p))
1164 : {
1165 134 : pol = ZXQ_charpoly(a, pol, v);
1166 134 : return fl ? pol : mkvec2(pol, RgX_RgXQ_eval(a, emb, galpol));
1167 : }
1168 90 : return NULL;
1169 : }
1170 :
1171 : static GEN
1172 166 : galoissubfieldcm(GEN G, long fl)
1173 : {
1174 166 : 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 166 : ulong p = 1009;
1178 :
1179 166 : galpol = gal_get_pol(G);
1180 166 : n = degpol(galpol);
1181 166 : v = varn(galpol);
1182 166 : 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 166 : checkgroup(G, &elts);
1187 166 : elts = gen_sort_shallow(elts,(void*)vecsmall_lexcmp,cmp_nodata);
1188 166 : H = vecsmall_ei(n,1); /* indices of elements of H */
1189 166 : Hset = zero_F2v(n);
1190 166 : F2v_set(Hset,1);
1191 166 : nH = 1;
1192 1124 : for (i=2; i<=n; i++)
1193 : {
1194 958 : g = gel(elts,i);
1195 958 : c2 = perm_mul(c,perm_conj(g,c));
1196 958 : if (!F2v_coeff(Hset,c2[1]))
1197 : {
1198 134 : nH++;
1199 134 : H[nH] = c2[1];
1200 134 : F2v_set(Hset,c2[1]);
1201 : }
1202 : }
1203 : /* group generated */
1204 166 : gene = gcopy(H);
1205 166 : setlg(gene,nH+1);
1206 166 : i = 1; /* last element that has been multiplied by the generators */
1207 300 : while (i < nH)
1208 : {
1209 886 : for (j=1; j<lg(gene); j++)
1210 : {
1211 752 : g = gel(elts,gene[j]);
1212 752 : ind = g[H[i]]; /* index of the product */
1213 752 : 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 134 : i++;
1222 : }
1223 166 : H = cgetg(lg(gene), t_VEC);
1224 466 : for (i=1; i<lg(H); i++)
1225 300 : gel(H,i) = gel(elts,gene[i]);
1226 166 : sub = galoisfixedfield(G, H, 0, -1);
1227 :
1228 : /* compute a totally imaginary generator */
1229 166 : pol = gel(sub,1);
1230 166 : emb = liftpol_shallow(gel(sub,2));
1231 166 : d = degpol(pol);
1232 166 : if (!(ZX_deflate_order(pol)%2) && sturm(RgX_deflate(pol,2))==d/2)
1233 : {
1234 32 : setvarn(pol,v);
1235 32 : return fl==1 ? pol: mkvec2(pol,emb);
1236 : }
1237 :
1238 : /* compute action of c on the subfield from that on the large field */
1239 134 : c = galoispermtopol(G,c);
1240 134 : if (d<n)
1241 : {
1242 26 : GEN M = cgetg(d+1,t_MAT), contc, contM;
1243 26 : gel(M,1) = col_ei(n,1); a = pol_1(v);
1244 72 : for (i=2; i<=d; i++)
1245 : {
1246 46 : a = RgX_rem(QX_mul(a,emb), galpol);
1247 46 : gel(M,i) = RgX_to_RgC(a,n);
1248 : }
1249 26 : c = RgX_RgXQ_eval(emb,c,galpol);
1250 26 : c = Q_primitive_part(c,&contc);
1251 26 : c = RgX_to_RgC(c,n);
1252 26 : M = Q_primitive_part(M,&contM);
1253 26 : c = RgM_RgC_invimage(M,c);
1254 26 : if (contc)
1255 : {
1256 16 : if (contM) contc = gdiv(contc,contM);
1257 16 : c = RgV_Rg_mul(c, contc);
1258 : }
1259 10 : else if (contM) c = RgV_Rg_mul(c, ginv(contM));
1260 26 : c = RgV_to_RgX(c, v);
1261 : }
1262 :
1263 : /* search for a generator of the form c(b)-b */
1264 218 : for (i=1; i<d; i++)
1265 : {
1266 190 : a = try_imag(pol_xn(i,v),c,pol,v,p,emb,galpol,fl);
1267 190 : if (a) return a;
1268 84 : p = unextprime(p+1);
1269 : }
1270 28 : B = stoi(10);
1271 28 : b = pol_xn(d-1,v);
1272 : while(1)
1273 : {
1274 170 : for (i=2; i<lg(b); i++) gel(b,i) = randomi(B);
1275 34 : a = try_imag(b,c,pol,v,p,emb,galpol,fl);
1276 34 : if (a) return a;
1277 6 : p = unextprime(p+1);
1278 : }
1279 : return NULL;/*LCOV_EXCL_LINE*/
1280 : }
1281 :
1282 : static GEN
1283 112 : quadsubfieldcm(GEN pol, long fl)
1284 : {
1285 112 : GEN a = gel(pol,3), b = gel(pol,2), d, P;
1286 112 : long v = varn(pol);
1287 112 : if (mpodd(a))
1288 26 : { b = mului(4, b); d = gen_2; }
1289 : else
1290 86 : { a = divis(a,2); d = gen_1; }
1291 112 : P = deg2pol_shallow(gen_1, gen_0, subii(b, sqri(a)), v);
1292 112 : return fl==1 ? P: mkvec2(P, deg1pol_shallow(d,a,v));
1293 : }
1294 :
1295 : GEN
1296 915 : nfsubfieldscm(GEN nf, long fl)
1297 : {
1298 915 : pari_sp av = avma;
1299 : GEN fa, lambda, V, res, ro, a, aa, ev, minev, pol, G;
1300 915 : long i, j, n, r2, minj=0, prec, emax, emin, e, precbound, ero;
1301 :
1302 915 : subfields_cleanup(&nf, &pol, &n, &fa);
1303 870 : ro = cmsubfield_get_roots(pol, nf, n, &r2, &prec);
1304 870 : if (!ro) return gc_const(av, gen_0);
1305 : /* now r2 == 2*n */
1306 :
1307 508 : if (n==2) return gc_GEN(av, quadsubfieldcm(pol, fl));
1308 396 : G = galoisinit(nf? nf: pol, NULL);
1309 396 : if (G != gen_0) return gc_GEN(av, galoissubfieldcm(G, fl));
1310 :
1311 230 : ero = 0;
1312 1200 : for (i=1; i<lg(ro); i++)
1313 : {
1314 970 : e = 1+gexpo(gel(ro,i));
1315 970 : if (e > ero) ero = e;
1316 : }
1317 230 : ero++;
1318 230 : fa = subfields_get_fa(pol, nf, fa);
1319 :
1320 230 : emax = 1;
1321 230 : emin = -1;
1322 1200 : for (i=1; i<lg(ro); i++)
1323 3610 : for (j=i+1; j<lg(ro); j++)
1324 : {
1325 2640 : e = gexpo(gsub(gel(ro,i),gel(ro,j)));
1326 2640 : if (e > emax) emax = e;
1327 2640 : if (e < emin) emin = e;
1328 : }
1329 230 : precbound = n*(emax-emin) + gexpo(fa) + n*n + 5;
1330 230 : precbound = 3 + precbound/BITS_IN_LONG;
1331 230 : if (prec < precbound)
1332 : {
1333 0 : prec = precbound;
1334 0 : ro = rootsuptoconj(pol, prec);
1335 : }
1336 :
1337 230 : lambda = zerovec(lg(fa)-1);
1338 1200 : for (i=1; i<=r2; i++)
1339 : {
1340 970 : a = gel(ro,i);
1341 970 : aa = conj_i(a);
1342 6884 : for (j=1; j<lg(fa); j++)
1343 : {
1344 5914 : ev = cxnorm(poleval(poleval(gel(fa,j),aa),a));
1345 5914 : if (j==1 || cmprr(minev,ev)>0) { minj = j; minev = ev; }
1346 : }
1347 970 : gel(lambda,minj) = gen_m1;
1348 : }
1349 :
1350 230 : V = twoembequation(pol, fa, lambda);
1351 230 : if (lg(V)==1) { delete_var(); return gc_const(av, gen_0); }
1352 190 : res = subfield_generator(pol, V, 2*(lg(V)-1), ero, fl);
1353 190 : delete_var();
1354 190 : return gc_GEN(av, res);
1355 : }
1356 :
1357 : static int
1358 13910 : field_is_contained(GEN V, GEN W, int strict)
1359 : {
1360 : GEN VW;
1361 13910 : ulong p = 1073741827;
1362 : /* distinct overfield must have different dimension */
1363 13910 : if (strict && lg(V) == lg(W)) return 0;
1364 : /* dimension of overfield must be multiple */
1365 10395 : if ((lg(W)-1) % (lg(V)-1)) return 0;
1366 7430 : VW = shallowconcat(V,W);
1367 7430 : if (Flm_rank(ZM_to_Flm(VW,p),p) > lg(W)-1) return 0;
1368 3025 : 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 205 : maxgen_subfields(GEN pol, GEN fa, long flag)
1384 : {
1385 205 : pari_sp av = avma;
1386 205 : GEN principal, ismax, isgene, Lmax = NULL, Lgene, res, V, W, W1;
1387 205 : long i, i2, j, flmax, flgene, nbmax = 0, nbgene = 0;
1388 :
1389 205 : if (!flag) return cgetg(1,t_VEC);
1390 205 : flmax = (flag & subf_MAXIMAL)!=0;
1391 205 : flgene = (flag & subf_GENERATING)!=0;
1392 :
1393 : /* compute principal subfields */
1394 205 : principal = cgetg(lg(fa),t_VEC);
1395 1830 : for (i=1; i<lg(fa); i++)
1396 1625 : gel(principal,i) = twoembequation(pol, fa, vec_ei(lg(fa)-1,i));
1397 205 : principal = gen_sort_uniq(principal, (void*)&cmp_universal, &cmp_nodata);
1398 : /* remove nf and duplicates (sort_uniq possibly not enough) */
1399 205 : i2 = 1;
1400 1210 : for (i=1; i<lg(principal)-1; i++)
1401 : {
1402 1005 : long dup = 0;
1403 1005 : V = gel(principal,i);
1404 1005 : j = i2-1;
1405 2055 : while (j > 0 && lg(gel(principal,j)) == lg(V))
1406 : {
1407 1050 : if (field_is_contained(gel(principal,j),V,0)) { dup=1; break; }
1408 1050 : j--;
1409 : }
1410 1005 : if (!dup) gel(principal,i2++) = V;
1411 : }
1412 205 : setlg(principal, i2);
1413 :
1414 : /* a subfield is generating iff all overfields contain the first overfield */
1415 205 : ismax = cgetg(lg(principal),t_VECSMALL);
1416 205 : isgene = cgetg(lg(principal),t_VECSMALL);
1417 1210 : for (i=1; i<lg(principal); i++)
1418 : {
1419 1005 : V = gel(principal,i);
1420 1005 : ismax[i] = flmax;
1421 1005 : isgene[i] = flgene;
1422 1005 : W1 = NULL; /* intersection of strict overfields */
1423 3470 : for (j=i+1; j<lg(principal); j++)
1424 : {
1425 2640 : W = gel(principal,j);
1426 2640 : if (!field_is_contained(V,W,1)) continue;
1427 495 : ismax[i] = 0;
1428 495 : if (!flgene) break;
1429 345 : if (!W1) { W1 = W; continue; }
1430 135 : if (!field_is_contained(W1,W,1))
1431 : {
1432 45 : W1 = intersect(W1,W);
1433 45 : if (lg(W1)==lg(V)) { isgene[i]=0; break; }
1434 : }
1435 : }
1436 : }
1437 :
1438 1210 : for (i=1; i<lg(principal); i++)
1439 : {
1440 1005 : nbmax += ismax[i];
1441 1005 : nbgene += isgene[i];
1442 : }
1443 :
1444 205 : if (flmax)
1445 : {
1446 70 : Lmax = cgetg(nbmax+1, t_VEC);
1447 70 : j=1;
1448 370 : for (i=1; i<lg(principal); i++)
1449 300 : if (ismax[i]) gel(Lmax,j++) = gel(principal,i);
1450 : }
1451 :
1452 205 : if (flgene)
1453 : {
1454 135 : Lgene = cgetg(nbgene+1, t_VEC);
1455 135 : j=1;
1456 840 : for (i=1; i<lg(principal); i++)
1457 705 : if (isgene[i]) gel(Lgene,j++) = gel(principal,i);
1458 : }
1459 :
1460 205 : if (!flgene) res = Lmax;
1461 135 : else if (!flmax) res = Lgene;
1462 0 : else res = mkvec2(Lmax,Lgene);
1463 205 : return gc_GEN(av, res);
1464 : }
1465 :
1466 : GEN
1467 110 : nfsubfieldsmax(GEN nf, long fl)
1468 : {
1469 110 : pari_sp av = avma;
1470 : GEN pol, fa, Lmax, V;
1471 : long n, i, ero;
1472 :
1473 110 : subfields_cleanup(&nf, &pol, &n, &fa);
1474 110 : if (n==1) retgc_const(av, cgetg(1, t_VEC));
1475 100 : if (uisprime(n))
1476 45 : return gc_GEN(av, fl==1 ? mkvec(pol_x(varn(pol)))
1477 15 : : mkvec(mkvec2(pol_x(varn(pol)),gen_0)));
1478 70 : ero = subfields_get_ero(pol, nf);
1479 70 : fa = subfields_get_fa(pol, nf, fa);
1480 70 : Lmax = maxgen_subfields(pol, fa, subf_MAXIMAL);
1481 220 : for (i=1; i<lg(Lmax); i++)
1482 : {
1483 150 : V = gel(Lmax,i);
1484 150 : gel(Lmax,i) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1485 : }
1486 70 : delete_var();
1487 70 : return gc_GEN(av, Lmax);
1488 : }
1489 :
1490 : static void
1491 1260 : heap_climb(GEN* H, long i)
1492 : {
1493 : long j;
1494 1260 : if (i==1) return;
1495 930 : j = i/2;
1496 930 : if (cmp_universal(gel(*H,i),gel(*H,j)) > 0)
1497 : {
1498 380 : swap(gel(*H,i), gel(*H,j));
1499 380 : return heap_climb(H,j);
1500 : }
1501 : }
1502 :
1503 : static void
1504 880 : heap_push(GEN* H, long *len, GEN x)
1505 : {
1506 880 : if (*len+1 == lg(*H))
1507 : {
1508 10 : GEN H2 = zerovec(2*(*len));
1509 : long i;
1510 110 : for(i=1; i<lg(*H); i++)
1511 100 : gel(H2,i) = gel(*H,i);
1512 10 : *H = H2;
1513 : }
1514 880 : (*len)++;
1515 880 : gel(*H,*len) = x;
1516 880 : return heap_climb(H,*len);
1517 : }
1518 :
1519 : static void
1520 1835 : heap_descend(GEN* H, long len, long i)
1521 : {
1522 1835 : long maxi = i, j = 2*i;
1523 1835 : if (j > len) return;
1524 955 : if (cmp_universal(gel(*H,j),gel(*H,i)) > 0) maxi = j;
1525 955 : j++;
1526 955 : if (j<=len && cmp_universal(gel(*H,j),gel(*H,maxi))>0) maxi = j;
1527 955 : if (maxi == i) return;
1528 820 : swap(gel(*H,i), gel(*H,maxi));
1529 820 : return heap_descend(H,len,maxi);
1530 : }
1531 :
1532 : static void
1533 1015 : heap_pop(GEN *H, long *len, GEN* top)
1534 : {
1535 1015 : *top = gel(*H,1);
1536 1015 : gel(*H,1) = gel(*H,*len);
1537 1015 : (*len)--;
1538 1015 : return heap_descend(H,*len,1);
1539 : };
1540 :
1541 : static GEN
1542 185 : nfsubfields_fa(GEN nf, long d, long fl)
1543 : {
1544 185 : 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 185 : subfields_cleanup(&nf, &pol, &n, &fa); vp = varn(pol);
1549 170 : if (d && (d<1 || d>n || n%d)) return gc_GEN(av, cgetg(1,t_VEC));
1550 175 : if (!d && uisprime(n)) return gc_GEN(av,
1551 0 : fl==1 ? mkvec2( pol_x(varn(pol)), pol)
1552 10 : : mkvec2( mkvec2(pol_x(vp),pol_0(vp)), mkvec2(pol,pol_x(vp))));
1553 165 : if (n==1 || d==1) return gc_GEN(av,
1554 10 : fl==1 ? mkvec(pol_x(varn(pol))): _subfield(pol_x(vp),pol_0(vp)));
1555 155 : if (d==n) return gc_GEN(av,
1556 10 : fl==1 ? mkvec(pol): _subfield(pol,pol_x(vp)));
1557 135 : ero = subfields_get_ero(pol, nf);
1558 135 : fa = subfields_get_fa(pol, nf, fa);
1559 135 : gene = maxgen_subfields(pol, fa, subf_GENERATING);
1560 :
1561 135 : if (d)
1562 : {
1563 : /* keep only generating subfields of degree a multiple of d */
1564 10 : j=1;
1565 105 : for (i=1; i<lg(gene); i++)
1566 95 : if ((lg(gel(gene,i))-1) % d == 0)
1567 : {
1568 70 : gel(gene,j) = gel(gene,i);
1569 70 : j++;
1570 : }
1571 10 : setlg(gene,j);
1572 : }
1573 135 : r = lg(gene)-1;
1574 :
1575 135 : res = zerovec(10);
1576 135 : nres = 0;
1577 135 : H = zerovec(10);
1578 135 : gel(H,1) = mkvec3(matid(n),zero_F2v(r),mkvecsmall(0));
1579 135 : len = 1;
1580 :
1581 1150 : while (len>0)
1582 : {
1583 1015 : heap_pop(&H, &len, &data);
1584 1015 : V = gel(data,1);
1585 1015 : v = gel(data,2);
1586 1015 : s = gel(data,3)[1];
1587 4395 : for (i=s+1; i<=r; i++)
1588 3380 : if (!F2v_coeff(v,i))
1589 : {
1590 2625 : W = vec_Q_primpart(intersect(V, gel(gene,i)));
1591 2625 : w = F2v_copy(v);
1592 2625 : F2v_set(w, i);
1593 2625 : newfield = 1;
1594 12950 : for (j=1; j<=r; j++)
1595 12000 : if (!F2v_coeff(w,j) && field_is_contained(W,gel(gene,j),1))
1596 : {
1597 2440 : if (j<i) { newfield = 0; break; }
1598 765 : F2v_set(w,j);
1599 : }
1600 2625 : if (newfield && (!d || (lg(W)-1)%d==0)) heap_push(&H, &len, mkvec3(W,w,mkvecsmall(i)));
1601 : }
1602 :
1603 1015 : if (!d || lg(V)-1==d)
1604 : {
1605 1005 : nres++;
1606 1005 : if (nres == lg(res))
1607 : {
1608 20 : res2 = zerovec(2*lg(res));
1609 280 : for(j=1; j<lg(res); j++) gel(res2,j) = gel(res,j);
1610 20 : res = res2;
1611 : }
1612 1005 : gel(res,nres) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1613 : }
1614 : }
1615 135 : setlg(res,nres+1);
1616 135 : vecreverse_inplace(res);
1617 :
1618 135 : delete_var();
1619 135 : return gc_GEN(av, res);
1620 : }
|
