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 32883 : set_avma(av);
106 : }
107 : }
108 33233 : gel(Yp,lY) = Zp;
109 :
110 33233 : K = divisorsu(n[1]); lK = lg(K);
111 141377 : for (i=1; i<lK; i++)
112 : {
113 108141 : long ngcd = n[1], k = K[i], dk = B->size*k, lpn = 0;
114 115498 : 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 108141 : 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 33236 : 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 214959 : for (l=1,lp=0; l<=n[j]; l++)
239 : {
240 177880 : lp++; if (lp > ki) lp = 1;
241 177880 : 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 40410 : ok_coeffs(GEN g,GEN M)
322 : {
323 40410 : long i, lg = lg(g)-1; /* g is monic, and cst term is ok */
324 98119 : for (i=3; i<lg; i++)
325 64849 : if (abscmpii(gel(g,i), gel(M,i)) > 0) return 0;
326 33270 : return 1;
327 : }
328 :
329 : /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
330 : static GEN
331 174095 : trace(GEN x, GEN Trq, GEN p)
332 : {
333 : long i, l;
334 : GEN s;
335 174095 : if (typ(x) == t_INT) return Fp_mul(x, gel(Trq,1), p);
336 174095 : l = lg(x)-1; if (l == 1) return gen_0;
337 174095 : x++; s = mulii(gel(x,1), gel(Trq,1));
338 883924 : for (i=2; i<l; i++)
339 709834 : 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 : long i,l;
348 : GEN y;
349 36581 : if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
350 36581 : y = cgetg_copy(x, &l); y[1] = x[1];
351 210664 : for (i=2; i<l; i++) gel(y,i) = trace(gel(x,i),Trq,p);
352 36570 : return normalizepol_lg(y, l);
353 : }
354 :
355 : /* Find h in Fp[X] such that h(a[i]) = listdelta[i] for all modular factors
356 : * ff[i], where a[i] is a fixed root of ff[i] in Fq = Z[Y]/(p,T) [namely the
357 : * first one in FpX_factorff_irred output]. Let f = ff[i], A the given root,
358 : * then h mod f is Tr_Fq/Fp ( h(A) f(X)/(X-A)f'(A) ), most of the expression
359 : * being precomputed. The complete h is recovered via chinese remaindering */
360 : static GEN
361 32508 : chinese_retrieve_pol(GEN DATA, primedata *S, GEN listdelta)
362 : {
363 32508 : GEN interp, bezoutC, h, p = S->p, pol = FpX_red(gel(DATA,1), p);
364 : long i, l;
365 32506 : interp = gel(DATA,9);
366 32506 : bezoutC= gel(DATA,6);
367 :
368 32506 : h = NULL; l = lg(interp);
369 69083 : for (i=1; i<l; i++)
370 : { /* h(firstroot[i]) = listdelta[i] */
371 36580 : GEN t = FqX_Fq_mul(gel(interp,i), gel(listdelta,i), S->T, p);
372 36581 : t = poltrace(t, gel(S->Trk,i), p);
373 36580 : t = FpX_mul(t, gel(bezoutC,i), p);
374 36578 : h = h? FpX_add(h,t,p): t;
375 : }
376 32503 : return FpX_rem(h, pol, p);
377 : }
378 :
379 : /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
380 : * (cf subfield) was a block system; then there
381 : * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
382 : * in Fp[X] (cf chinese_retrieve_pol). Try to lift it; den is a
383 : * multiplicative bound for denominator of lift. */
384 : static GEN
385 32508 : embedding(GEN g, GEN DATA, primedata *S, GEN den, GEN listdelta)
386 : {
387 32508 : GEN TR, w0_Q, w0, w1_Q, w1, wpow, h0, gp, T, q2, q, maxp, a, p = S->p;
388 : long rt;
389 : pari_sp av;
390 :
391 32508 : T = gel(DATA,1); rt = brent_kung_optpow(degpol(T), 4, 3);
392 32508 : maxp= gel(DATA,7);
393 32508 : gp = RgX_deriv(g); av = avma;
394 32507 : w0 = chinese_retrieve_pol(DATA, S, listdelta);
395 32505 : w0_Q = centermod(gmul(w0,den), p);
396 32506 : h0 = FpXQ_inv(FpX_FpXQ_eval(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
397 32507 : wpow = NULL; q = sqri(p);
398 : for(;;)
399 : {/* Given g,w0,h0 in Z[x], s.t. h0.g'(w0) = 1 and g(w0) = 0 mod (T,p), find
400 : * [w1,h1] satisfying the same conditions mod p^2, [w1,h1] = [w0,h0] (mod p)
401 : * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
402 115854 : if (DEBUGLEVEL>1)
403 0 : err_printf("lifting embedding mod p^k = %Ps^%ld\n",S->p, Z_pval(q,S->p));
404 :
405 : /* w1 := w0 - h0 g(w0) mod (T,q) */
406 115854 : if (wpow) a = FpX_FpXQV_eval(g,wpow, T,q);
407 32505 : else a = FpX_FpXQ_eval(g,w0, T,q); /* first time */
408 : /* now, a = 0 (p) */
409 115859 : a = FpXQ_mul(ZX_neg(h0), ZX_Z_divexact(a, p), T,p);
410 115861 : w1 = ZX_add(w0, ZX_Z_mul(a, p));
411 :
412 115835 : w1_Q = centermod(ZX_Z_mul(w1, remii(den,q)), q);
413 115857 : if (ZX_equal(w1_Q, w0_Q))
414 : {
415 21731 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
416 21729 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
417 : }
418 94130 : else if (cmpii(q,maxp) > 0)
419 : {
420 14974 : GEN G = is_pm1(den)? g: RgX_rescale(g,den);
421 14974 : if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
422 497 : if (DEBUGLEVEL) err_printf("coeff too big for embedding\n");
423 497 : return NULL;
424 : }
425 83353 : gerepileall(av, 5, &w1,&h0,&w1_Q,&q,&p);
426 83359 : q2 = sqri(q);
427 83347 : wpow = FpXQ_powers(w1, rt, T, q2);
428 : /* h0 := h0 * (2 - h0 g'(w1)) mod (T,q)
429 : * = h0 + h0 * (1 - h0 g'(w1)) */
430 83350 : a = FpXQ_mul(ZX_neg(h0), FpX_FpXQV_eval(gp, FpXV_red(wpow,q),T,q), T,q);
431 83351 : a = ZX_Z_add_shallow(a, gen_1); /* 1 - h0 g'(w1) = 0 (p) */
432 83346 : a = FpXQ_mul(h0, ZX_Z_divexact(a, p), T,p);
433 83352 : h0 = ZX_add(h0, ZX_Z_mul(a, p));
434 83349 : w0 = w1; w0_Q = w1_Q; p = q; q = q2;
435 : }
436 32010 : TR = gel(DATA,5);
437 32010 : if (!gequal0(TR)) w1_Q = RgX_translate(w1_Q, TR);
438 32009 : return gdiv(w1_Q,den);
439 : }
440 :
441 : /* return U list of polynomials s.t U[i] = 1 mod fk[i] and 0 mod fk[j] for all
442 : * other j */
443 : static GEN
444 30840 : get_bezout(GEN pol, GEN fk, GEN p)
445 : {
446 30840 : long i, l = lg(fk);
447 30840 : GEN A, B, d, u, v, U = cgetg(l, t_VEC);
448 63333 : for (i=1; i<l; i++)
449 : {
450 32493 : A = gel(fk,i);
451 32493 : B = FpX_div(pol, A, p);
452 32492 : d = FpX_extgcd(A,B,p, &u, &v);
453 32493 : if (degpol(d) > 0) pari_err_COPRIME("get_bezout",A,B);
454 32493 : d = constant_coeff(d);
455 32493 : if (!gequal1(d)) v = FpX_Fp_div(v, d, p);
456 32493 : gel(U,i) = FpX_mul(B,v, p);
457 : }
458 30840 : return U;
459 : }
460 :
461 : static GEN
462 30841 : init_traces(GEN ff, GEN T, GEN p)
463 : {
464 30841 : long N = degpol(T),i,j,k, r = lg(ff);
465 30841 : GEN Frob = FpX_matFrobenius(T,p);
466 : GEN y,p1,p2,Trk,pow,pow1;
467 :
468 30842 : k = degpol(gel(ff,r-1)); /* largest degree in modular factorization */
469 30842 : pow = cgetg(k+1, t_VEC);
470 30842 : gel(pow,1) = gen_0; /* dummy */
471 30842 : gel(pow,2) = Frob;
472 30842 : pow1= cgetg(k+1, t_VEC); /* 1st line */
473 110614 : for (i=3; i<=k; i++)
474 79772 : gel(pow,i) = FpM_mul(gel(pow,i-1), Frob, p);
475 30842 : gel(pow1,1) = gen_0; /* dummy */
476 141455 : for (i=2; i<=k; i++)
477 : {
478 110612 : p1 = cgetg(N+1, t_VEC);
479 110613 : gel(pow1,i) = p1; p2 = gel(pow,i);
480 778367 : for (j=1; j<=N; j++) gel(p1,j) = gcoeff(p2,1,j);
481 : }
482 :
483 : /* Trk[i] = line 1 of x -> x + x^p + ... + x^{p^(i-1)} */
484 30843 : Trk = pow; /* re-use (destroy) pow */
485 30843 : gel(Trk,1) = vec_ei(N,1);
486 141441 : for (i=2; i<=k; i++)
487 110604 : gel(Trk,i) = gadd(gel(Trk,i-1), gel(pow1,i));
488 30837 : y = cgetg(r, t_VEC);
489 63332 : for (i=1; i<r; i++) y[i] = Trk[degpol(gel(ff,i))];
490 30840 : return y;
491 : }
492 :
493 : static void
494 30842 : init_primedata(primedata *S)
495 : {
496 30842 : long i, j, l, lff = lg(S->ff), v = fetch_var(), N = degpol(S->pol);
497 30842 : GEN T, p = S->p;
498 :
499 30842 : if (S->lcm == degpol(gel(S->ff,lff-1)))
500 : {
501 30828 : T = leafcopy(gel(S->ff,lff-1));
502 30828 : setvarn(T, v);
503 : }
504 : else
505 14 : T = init_Fq(p, S->lcm, v);
506 30842 : S->T = T;
507 30842 : S->firstroot = cgetg(lff, t_VECSMALL);
508 30842 : S->interp = cgetg(lff, t_VEC);
509 30842 : S->fk = cgetg(N+1, t_VEC);
510 63334 : for (l=1,j=1; j<lff; j++)
511 : { /* compute roots and fix ordering (Frobenius cycles) */
512 32494 : GEN F = gel(S->ff, j), deg1 = FpX_factorff_irred(F, T,p);
513 32493 : GEN H = gel(deg1,1), a = Fq_neg(constant_coeff(H), T,p);
514 32491 : GEN Q = FqX_div(F, H, T,p);
515 32494 : GEN q = Fq_inv(FqX_eval(Q, a, T,p), T,p);
516 32494 : gel(S->interp,j) = FqX_Fq_mul(Q, q, T,p); /* = 1 at a, 0 at other roots */
517 32492 : S->firstroot[j] = l;
518 182910 : for (i=1; i<lg(deg1); i++,l++) gel(S->fk, l) = gel(deg1, i);
519 : }
520 30840 : S->Trk = init_traces(S->ff, T,p);
521 30840 : S->bezoutC = get_bezout(S->pol, S->ff, p);
522 30840 : }
523 :
524 : static int
525 30856 : choose_prime(primedata *S, GEN pol)
526 : {
527 30856 : long i, j, k, r, lcm, oldr, K, N = degpol(pol);
528 : ulong p, pp;
529 : GEN Z, ff, n, oldn;
530 : pari_sp av;
531 : forprime_t T;
532 :
533 30856 : u_forprime_init(&T, (N*N) >> 2, ULONG_MAX);
534 30856 : oldr = S->lcm = LONG_MAX;
535 30856 : S->ff = oldn = NULL; pp = 0; /* gcc -Wall */
536 30856 : av = avma; K = N + 10;
537 99710 : for(k = 1; k < K || !S->ff; k++,set_avma(av))
538 : {
539 : GEN Tp;
540 98254 : if (k > 5 * N) return 0;
541 : do
542 : {
543 116654 : p = u_forprime_next(&T);
544 116653 : Tp = ZX_to_Flx(pol, p);
545 : }
546 116651 : while (!Flx_is_squarefree(Tp, p));
547 98251 : ff = gel(Flx_factor(Tp, p), 1);
548 98254 : r = lg(ff)-1;
549 98254 : if (r == N || r >= BIL) continue;
550 :
551 95804 : n = cgetg(r+1, t_VECSMALL); lcm = n[1] = degpol(gel(ff,1));
552 256044 : for (j=2; j<=r; j++) { n[j] = degpol(gel(ff,j)); lcm = ulcm(lcm, n[j]); }
553 95804 : if (r > oldr || (r == oldr && (lcm <= S->lcm || S->lcm > 2*N)))
554 45453 : continue;
555 50351 : if (DEBUGLEVEL) err_printf("p = %lu,\tlcm = %ld,\torbits: %Ps\n",p,lcm,n);
556 :
557 50351 : pp = p;
558 50351 : oldr = r;
559 50351 : oldn = n;
560 50351 : S->ff = ff;
561 50351 : S->lcm = lcm; if (r == 1) break;
562 20951 : av = avma;
563 : }
564 30856 : if (oldr > 6) return 0;
565 30842 : if (DEBUGLEVEL) err_printf("Chosen prime: p = %ld\n", pp);
566 30842 : FlxV_to_ZXV_inplace(S->ff);
567 30840 : S->p = utoipos(pp);
568 30840 : S->pol = FpX_red(pol, S->p); init_primedata(S);
569 30840 : n = oldn; r = lg(n); S->Z = Z = cgetg(r,t_VEC);
570 63334 : for (k=0,i=1; i<r; i++)
571 : {
572 32493 : GEN t = cgetg(n[i]+1, t_VECSMALL); gel(Z,i) = t;
573 182919 : for (j=1; j<=n[i]; j++) t[j] = ++k;
574 : }
575 30841 : return 1;
576 : }
577 :
578 : /* maxroot t_REAL */
579 : static GEN
580 31948 : bound_for_coeff(long m, GEN rr, GEN *maxroot)
581 : {
582 31948 : long i,r1, lrr=lg(rr);
583 31948 : GEN p1,b1,b2,B,M, C = matpascal(m-1);
584 :
585 48594 : for (r1=1; r1 < lrr; r1++)
586 47145 : if (typ(gel(rr,r1)) != t_REAL) break;
587 31948 : r1--;
588 :
589 31948 : rr = gabs(rr,0); *maxroot = vecmax(rr);
590 196219 : for (i=1; i<lrr; i++)
591 164272 : if (gcmp(gel(rr,i), gen_1) < 0) gel(rr,i) = gen_1;
592 48591 : for (b1=gen_1,i=1; i<=r1; i++) b1 = gmul(b1, gel(rr,i));
593 179564 : for (b2=gen_1 ; i<lrr; i++) b2 = gmul(b2, gel(rr,i));
594 31945 : B = gmul(b1, gsqr(b2)); /* Mahler measure */
595 31948 : M = cgetg(m+2, t_VEC); gel(M,1) = gel(M,2) = gen_0; /* unused */
596 79783 : for (i=1; i<m; i++)
597 : {
598 47836 : p1 = gadd(gmul(gcoeff(C, m, i+1), B),/* binom(m-1, i) */
599 47837 : gcoeff(C, m, i)); /* binom(m-1, i-1) */
600 47834 : gel(M,i+2) = ceil_safe(p1);
601 : }
602 31946 : return M;
603 : }
604 :
605 : /* d = requested degree for subfield. Return DATA, valid for given pol, S and d
606 : * If DATA != NULL, translate pol [ --> pol(X+1) ] and update DATA
607 : * 1: polynomial pol
608 : * 2: p^e (for Hensel lifts) such that p^e > max(M),
609 : * 3: Hensel lift to precision p^e of DATA[4]
610 : * 4: roots of pol in F_(p^S->lcm),
611 : * 5: number of polynomial changes (translations)
612 : * 6: Bezout coefficients attached to the S->ff[i]
613 : * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod pol.
614 : * 8: bound M for polynomials defining subfields x PD->den
615 : * 9: *[i] = interpolation polynomial for S->ff[i] [= 1 on the first root
616 : S->firstroot[i], 0 on the others] */
617 : static void
618 31948 : compute_data(blockdata *B)
619 : {
620 : GEN ffL, roo, pe, p1, p2, fk, fhk, MM, maxroot, pol;
621 31948 : primedata *S = B->S;
622 31948 : GEN p = S->p, T = S->T, ff = S->ff, DATA = B->DATA;
623 31948 : long i, j, l, e, N, lff = lg(ff);
624 :
625 31948 : if (DEBUGLEVEL>1) err_printf("Entering compute_data()\n\n");
626 31948 : pol = B->PD->pol; N = degpol(pol);
627 31948 : roo = B->PD->roo;
628 31948 : if (DATA)
629 : {
630 763 : GEN Xm1 = gsub(pol_x(varn(pol)), gen_1);
631 763 : GEN TR = addiu(gel(DATA,5), 1);
632 763 : GEN mTR = negi(TR), interp, bezoutC;
633 :
634 763 : if (DEBUGLEVEL>1) err_printf("... update (translate) an existing DATA\n\n");
635 :
636 763 : gel(DATA,5) = TR;
637 763 : pol = RgX_translate(gel(DATA,1), gen_m1);
638 763 : p1 = cgetg_copy(roo, &l);
639 9527 : for (i=1; i<l; i++) gel(p1,i) = gadd(TR, gel(roo,i));
640 763 : roo = p1;
641 :
642 763 : fk = gel(DATA,4); l = lg(fk);
643 9527 : for (i=1; i<l; i++) gel(fk,i) = gsub(Xm1, gel(fk,i));
644 :
645 763 : bezoutC = gel(DATA,6); l = lg(bezoutC);
646 763 : interp = gel(DATA,9);
647 2212 : for (i=1; i<l; i++)
648 : {
649 1449 : if (degpol(gel(interp,i)) > 0) /* do not turn pol_1(0) into gen_1 */
650 : {
651 1449 : p1 = RgX_translate(gel(interp,i), gen_m1);
652 1449 : gel(interp,i) = FpXX_red(p1, p);
653 : }
654 1449 : if (degpol(gel(bezoutC,i)) > 0)
655 : {
656 1351 : p1 = RgX_translate(gel(bezoutC,i), gen_m1);
657 1351 : gel(bezoutC,i) = FpXX_red(p1, p);
658 : }
659 : }
660 763 : ff = cgetg(lff, t_VEC); /* copy, do not overwrite! */
661 2212 : for (i=1; i<lff; i++)
662 1449 : gel(ff,i) = FpX_red(RgX_translate(gel(S->ff,i), mTR), p);
663 : }
664 : else
665 : {
666 31185 : DATA = cgetg(10,t_VEC);
667 31185 : fk = S->fk;
668 31185 : gel(DATA,5) = gen_0;
669 31185 : gel(DATA,6) = leafcopy(S->bezoutC);
670 31185 : gel(DATA,9) = leafcopy(S->interp);
671 : }
672 31948 : gel(DATA,1) = pol;
673 31948 : MM = gmul2n(bound_for_coeff(B->d, roo, &maxroot), 1);
674 31946 : gel(DATA,8) = MM;
675 31946 : e = logintall(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [d-1 test] */
676 31948 : gel(DATA,2) = pe;
677 31948 : gel(DATA,4) = roots_from_deg1(fk);
678 :
679 : /* compute fhk = ZpX_liftfact(pol,fk,T,p,e,pe) in 2 steps
680 : * 1) lift in Zp to precision p^e */
681 31948 : ffL = ZpX_liftfact(pol, ff, pe, p, e);
682 31948 : fhk = NULL;
683 66416 : for (l=i=1; i<lff; i++)
684 : { /* 2) lift factorization of ff[i] in Qp[X] / T */
685 34468 : GEN F, L = gel(ffL,i);
686 34468 : long di = degpol(L);
687 34468 : F = cgetg(di+1, t_VEC);
688 199010 : for (j=1; j<=di; j++) F[j] = fk[l++];
689 34468 : L = ZqX_liftfact(L, F, T, pe, p, e);
690 34468 : fhk = fhk? shallowconcat(fhk, L): L;
691 : }
692 31948 : gel(DATA,3) = roots_from_deg1(fhk);
693 :
694 31948 : p1 = mulur(N, powruhalf(utor(N-1,DEFAULTPREC), N-1));
695 31948 : p2 = powru(maxroot, B->size + N*(N-1)/2);
696 31947 : p1 = divrr(mulrr(p1,p2), gsqrt(B->PD->dis,DEFAULTPREC));
697 31948 : gel(DATA,7) = mulii(shifti(ceil_safe(p1), 1), B->PD->den);
698 :
699 31946 : if (DEBUGLEVEL>1) {
700 0 : err_printf("f = %Ps\n",DATA[1]);
701 0 : err_printf("p = %Ps, lift to p^%ld\n", p, e);
702 0 : err_printf("2 * Hadamard bound * ind = %Ps\n",DATA[7]);
703 0 : err_printf("2 * M = %Ps\n",DATA[8]);
704 : }
705 31946 : if (B->DATA) { DATA = gclone(DATA); if (isclone(B->DATA)) gunclone(B->DATA); }
706 31946 : B->DATA = DATA;
707 31946 : }
708 :
709 : /* g = polynomial, h = embedding. Return [[g,h]] */
710 : static GEN
711 1463 : _subfield(GEN g, GEN h) { return mkvec(mkvec2(g,h)); }
712 :
713 : /* Return a subfield, gen_0 [ change p ] or NULL [ not a subfield ] */
714 : static GEN
715 54278 : subfield(GEN A, blockdata *B)
716 : {
717 54278 : long N, i, j, d, lf, m = lg(A)-1;
718 : GEN M, pe, pol, fhk, g, e, d_1_term, delta, listdelta, whichdelta;
719 54278 : GEN T = B->S->T, p = B->S->p, firstroot = B->S->firstroot;
720 :
721 54278 : pol= gel(B->DATA,1); N = degpol(pol); d = N/m; /* m | N */
722 54278 : pe = gel(B->DATA,2);
723 54278 : fhk= gel(B->DATA,3);
724 54278 : M = gel(B->DATA,8);
725 :
726 54278 : delta = cgetg(m+1,t_VEC);
727 54278 : whichdelta = cgetg(N+1, t_VECSMALL);
728 54284 : d_1_term = gen_0;
729 344622 : for (i=1; i<=m; i++)
730 : {
731 290348 : GEN Ai = gel(A,i), p1 = gel(fhk,Ai[1]);
732 902512 : for (j=2; j<=d; j++)
733 612166 : p1 = Fq_mul(p1, gel(fhk,Ai[j]), T, pe);
734 290346 : gel(delta,i) = p1;
735 290346 : if (DEBUGLEVEL>5) err_printf("delta[%ld] = %Ps\n",i,p1);
736 : /* g = prod (X - delta[i])
737 : * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
738 : /* fk[k] belongs to block number whichdelta[k] */
739 1192856 : for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
740 290346 : if (typ(p1) == t_POL) p1 = constant_coeff(p1);
741 290345 : d_1_term = addii(d_1_term, p1);
742 : }
743 54274 : d_1_term = centermod(d_1_term, pe); /* Tr(g) */
744 54273 : if (abscmpii(d_1_term, gel(M,3)) > 0) {
745 13867 : if (DEBUGLEVEL>1) err_printf("d-1 test failed\n");
746 13867 : return NULL;
747 : }
748 40407 : g = FqV_roots_to_pol(delta, T, pe, 0);
749 40411 : g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
750 40410 : if (!ok_coeffs(g,M)) {
751 7140 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
752 7140 : if (DEBUGLEVEL>1) err_printf("coeff too big for pol g(x)\n");
753 7140 : return NULL;
754 : }
755 33270 : if (!FpX_is_squarefree(g, p)) {
756 763 : if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
757 763 : if (DEBUGLEVEL>1) err_printf("changing f(x): p divides disc(g)\n");
758 763 : compute_data(B);
759 763 : return subfield(A, B);
760 : }
761 :
762 32508 : lf = lg(firstroot); listdelta = cgetg(lf, t_VEC);
763 69090 : for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
764 32508 : if (DEBUGLEVEL) err_printf("candidate = %Ps\n", g);
765 32508 : e = embedding(g, B->DATA, B->S, B->PD->den, listdelta);
766 32505 : if (!e) return NULL;
767 32008 : if (DEBUGLEVEL) err_printf("... OK!\n");
768 32008 : return B->fl==1? mkvec(g):_subfield(g, e);
769 : }
770 :
771 : /* L list of current subfields, test whether potential block D is a block,
772 : * if so, append corresponding subfield */
773 : static GEN
774 53515 : test_block(blockdata *B, GEN L, GEN D)
775 : {
776 53515 : pari_sp av = avma;
777 53515 : GEN sub = subfield(D, B);
778 53512 : if (sub) {
779 32008 : GEN old = L;
780 32008 : L = gclone( L? shallowconcat(L, sub): sub );
781 32008 : guncloneNULL(old);
782 : }
783 53512 : return gc_const(av,L);
784 : }
785 :
786 : /* subfields of degree d */
787 : static GEN
788 31185 : subfields_of_given_degree(blockdata *B)
789 : {
790 31185 : pari_sp av = avma;
791 : GEN L;
792 :
793 31185 : if (DEBUGLEVEL) err_printf("\n* Look for subfields of degree %ld\n\n", B->d);
794 31185 : B->DATA = NULL; compute_data(B);
795 31183 : L = calc_block(B, B->S->Z, cgetg(1,t_VEC), NULL);
796 31184 : if (DEBUGLEVEL>9)
797 0 : err_printf("\nSubfields of degree %ld: %Ps\n", B->d, L? L: cgetg(1,t_VEC));
798 31184 : if (isclone(B->DATA)) gunclone(B->DATA);
799 31184 : return gc_const(av,L);
800 : }
801 :
802 : static void
803 35 : setvarn2(GEN t, long v) { setvarn(gel(t,1),v); setvarn(gel(t,2),v); }
804 : static GEN
805 29820 : fix_var(GEN x, long v, long fl)
806 : {
807 29820 : long i, l = lg(x);
808 29820 : if (!v) return x;
809 28 : if (fl)
810 42 : for (i = 1; i < l; i++) setvarn(gel(x,i), v);
811 : else
812 49 : for (i = 1; i < l; i++) setvarn2(gel(x,i), v);
813 28 : return x;
814 : }
815 :
816 : static void
817 30841 : subfields_poldata(GEN nf, GEN T, poldata *PD)
818 : {
819 : GEN L, dis;
820 :
821 30841 : PD->pol = T;
822 30841 : if (nf)
823 : {
824 168 : PD->den = nf_get_zkden(nf);
825 168 : PD->roo = nf_get_roots(nf);
826 168 : PD->dis = mulii(absi_shallow(nf_get_disc(nf)), sqri(nf_get_index(nf)));
827 : }
828 : else
829 : {
830 30673 : PD->den = initgaloisborne(T,NULL,nbits2prec(bit_accuracy(ZX_max_lg(T))), &L,NULL,&dis);
831 30674 : PD->roo = L;
832 30674 : PD->dis = absi_shallow(dis);
833 : }
834 30842 : }
835 :
836 : static GEN nfsubfields_fa(GEN nf, long d, long fl);
837 : static GEN
838 595 : subfieldsall(GEN nf0, long fl)
839 : {
840 595 : pari_sp av = avma;
841 : long N, ld, i, v;
842 : GEN nf, G, T, dg, LSB, NLSB;
843 : poldata PD;
844 : primedata S;
845 : blockdata B;
846 :
847 : /* much easier if nf is Galois (WSS) */
848 595 : G = galoisinit(nf0, NULL);
849 595 : T = get_nfpol(nf0, &nf);
850 595 : if (G != gen_0)
851 : {
852 : GEN L, S;
853 : long l;
854 56 : L = lift_shallow( galoissubfields(G, fl, varn(T)) );
855 56 : l = lg(L); S = cgetg(l, t_VECSMALL);
856 924 : for (i=1; i<l; i++) S[i] = lg(fl==1? gel(L,i): gmael(L,i,1));
857 56 : return gerepilecopy(av, vecpermute(L, vecsmall_indexsort(S)));
858 : }
859 539 : v = varn(T); N = degpol(T);
860 539 : dg = divisorsu(N); ld = lg(dg)-1;
861 539 : LSB = fl==1 ? mkvec(pol_x(v)): _subfield(pol_x(v), pol_0(v));
862 539 : if (ld <= 2)
863 : {
864 168 : if (ld == 2)
865 168 : LSB = shallowconcat(LSB, fl==1? mkvec(T): _subfield(T, pol_x(v)));
866 168 : return gerepilecopy(av, LSB);
867 : }
868 371 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
869 371 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, 0, fl); }
870 357 : subfields_poldata(nf, T, &PD);
871 :
872 357 : if (DEBUGLEVEL) err_printf("\n***** Entering subfields\n\npol = %Ps\n",T);
873 357 : B.PD = &PD;
874 357 : B.S = &S;
875 357 : B.N = N;
876 357 : B.fl = fl;
877 1057 : for (i=ld-1; i>1; i--)
878 : {
879 700 : B.size = dg[i];
880 700 : B.d = N / B.size;
881 700 : NLSB = subfields_of_given_degree(&B);
882 700 : if (NLSB) { LSB = gconcat(LSB, NLSB); gunclone(NLSB); }
883 : }
884 357 : (void)delete_var(); /* from init_primedata() */
885 357 : LSB = shallowconcat(LSB, fl==1? mkvec(T):_subfield(T, pol_x(0)));
886 357 : if (DEBUGLEVEL) err_printf("\n***** Leaving subfields\n\n");
887 357 : return fix_var(gerepilecopy(av, LSB), v, fl);
888 : }
889 :
890 : GEN
891 32753 : nfsubfields0(GEN nf0, long d, long fl)
892 : {
893 32753 : pari_sp av = avma;
894 : long N, v0;
895 : GEN nf, LSB, T, G;
896 : poldata PD;
897 : primedata S;
898 : blockdata B;
899 32753 : if (fl<0 || fl>1) pari_err_FLAG("nfsubfields");
900 32753 : if (typ(nf0)==t_VEC && lg(nf0)==3) return nfsubfields_fa(nf0, d, fl);
901 32508 : if (!d) return subfieldsall(nf0, fl);
902 :
903 : /* treat trivial cases */
904 31913 : T = get_nfpol(nf0, &nf); v0 = varn(T); N = degpol(T);
905 31913 : RgX_check_ZX(T,"nfsubfields");
906 31913 : if (d == N)
907 28 : return gerepilecopy(av, fl==1 ? mkvec(T) : _subfield(T, pol_x(v0)));
908 31885 : if (d == 1)
909 28 : return gerepilecopy(av, fl==1 ? mkvec(pol_x(v0)) : _subfield(pol_x(v0), zeropol(v0)));
910 31857 : if (d < 1 || d > N || N % d) return cgetg(1,t_VEC);
911 :
912 : /* much easier if nf is Galois (WSS) */
913 31815 : G = galoisinit(nf0, NULL);
914 31815 : if (G != gen_0)
915 : { /* Bingo */
916 1330 : GEN L = galoissubgroups(G), F;
917 1330 : long k,i, l = lg(L), o = N/d;
918 1330 : F = cgetg(l, t_VEC);
919 1330 : k = 1;
920 5936 : for (i=1; i<l; i++)
921 : {
922 4606 : GEN H = gel(L,i);
923 4606 : if (group_order(H) == o)
924 1540 : gel(F,k++) = lift_shallow(galoisfixedfield(G, gel(H,1), fl, v0));
925 : }
926 1330 : setlg(F, k);
927 1330 : return gerepilecopy(av, F);
928 : }
929 30485 : if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
930 30485 : if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, d, fl); }
931 30484 : subfields_poldata(nf, T, &PD);
932 30485 : B.PD = &PD;
933 30485 : B.S = &S;
934 30485 : B.N = N;
935 30485 : B.d = d;
936 30485 : B.size = N/d;
937 30485 : B.fl = fl;
938 30485 : LSB = subfields_of_given_degree(&B);
939 30484 : (void)delete_var(); /* from init_primedata */
940 30484 : set_avma(av);
941 30484 : if (!LSB) return cgetg(1, t_VEC);
942 29462 : G = gcopy(LSB); gunclone(LSB);
943 29463 : return fix_var(G, v0, fl);
944 : }
945 :
946 : GEN
947 329 : nfsubfields(GEN nf0, long d)
948 329 : { return nfsubfields0(nf0, d, 0); }
949 :
950 : /******************************/
951 : /* */
952 : /* Maximal CM subfield */
953 : /* Aurel Page (2019) */
954 : /* */
955 : /******************************/
956 :
957 : /* ero: maximum exponent+1 of roots of pol */
958 : static GEN
959 3325 : try_subfield_generator(GEN pol, GEN v, long e, long p, long ero, long fl)
960 : {
961 : GEN a, P, Q;
962 : long d, bound, i, B, bi, ed;
963 :
964 3325 : a = gtopolyrev(v, varn(pol));
965 3325 : P = Flxq_charpoly(ZX_to_Flx(a,p), ZX_to_Flx(pol,p), p);
966 3325 : Flx_ispower(P, e, p, &Q);
967 3325 : if (!Flx_is_squarefree(Q,p)) return NULL;
968 1701 : d = degpol(pol)/e;
969 1701 : B = 0;
970 26607 : for (i=1; i<lg(v); i++)
971 : {
972 24906 : bi = (i-1)*ero + expi(gel(v,i));
973 24906 : if (bi > B) B = bi;
974 : }
975 1701 : ed = expu(d);
976 1701 : B += ed+1;
977 1701 : bound = 0;
978 8253 : for (i=0; 2*i<=d; i++)
979 : {
980 6552 : if (!i) bi = d*B;
981 4851 : else bi = (d-i)*B + i*(3+ed-expu(i));
982 6552 : if (bi > bound) bound = bi;
983 : }
984 1701 : Q = ZXQ_minpoly(a,pol,d,bound);
985 1701 : return fl==1? Q: mkvec2(Q, a);
986 : }
987 :
988 : /* subfield sub of nf of degree d assuming:
989 : - V is contained in sub
990 : - V is not contained in a proper subfield of sub
991 : ero: maximum exponent+1 of roots of pol
992 : output as nfsubfields:
993 : - pair [g,h] where g absolute equation for the subfield and h expresses
994 : - one of the roots of g in terms of the generator of nf
995 : */
996 : static GEN
997 1876 : subfield_generator(GEN pol, GEN V, long d, long ero, long fl)
998 : {
999 1876 : long p, i, e, vp = varn(pol);
1000 1876 : GEN a = NULL, v = cgetg(lg(V),t_COL), B;
1001 :
1002 1876 : if (d==1) return fl ? pol_x(vp): mkvec2(pol_x(vp), pol_0(vp));
1003 1701 : e = degpol(pol)/d;
1004 1701 : p = 1009;
1005 3318 : for (i=1; i<lg(V); i++)
1006 : {
1007 3297 : a = try_subfield_generator(pol, gel(V,i), e, p, ero, fl);
1008 3297 : if (a) return a;
1009 1617 : p = unextprime(p+1);
1010 : }
1011 21 : B = stoi(10);
1012 : while(1)
1013 : {
1014 126 : for (i=1; i<lg(v); i++) gel(v,i) = randomi(B);
1015 28 : a = try_subfield_generator(pol, QM_QC_mul(V,v), e, p, ero, fl);
1016 28 : if (a) return a;
1017 7 : p = unextprime(p+1);
1018 : }
1019 : return NULL;/*LCOV_EXCL_LINE*/
1020 : }
1021 :
1022 : static GEN
1023 37961 : RgXY_to_RgC(GEN P, long dx, long dy)
1024 : {
1025 : GEN res, c;
1026 37961 : long i, j, k, d = degpol(P);
1027 37961 : if (d > dy) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1028 37961 : res = cgetg((dx+1)*(dy+1)+1, t_COL);
1029 37961 : k = 1;
1030 93618 : for (i=0; i<=d; i++)
1031 : {
1032 55657 : c = gel(P,i+2);
1033 55657 : if (typ(c)==t_POL)
1034 : {
1035 51408 : long dc = degpol(c);
1036 51408 : if (dc > dx) pari_err_BUG("RgXY_to_RgC [incorrect degree]");
1037 1021790 : for (j=0; j<=dc; j++)
1038 970382 : gel(res,k++) = gel(c,j+2);
1039 : } else
1040 : {
1041 4249 : gel(res,k++) = c; j=1;
1042 : }
1043 344365 : for ( ; j<=dx; j++)
1044 288708 : gel(res,k++) = gen_0;
1045 : }
1046 87304 : for( ; i<=dy; i++)
1047 1075242 : for (j=0; j<=dx; j++)
1048 1025899 : gel(res,k++) = gen_0;
1049 37961 : return res;
1050 : }
1051 :
1052 : /* lambda: t_VEC of t_INT; 0 means ignore this factor */
1053 : static GEN
1054 2583 : twoembequation(GEN pol, GEN fa, GEN lambda)
1055 : {
1056 : GEN m, vpolx, poly;
1057 2583 : long i,j, lfa = lg(fa), dx = degpol(pol);
1058 2583 : long vx = varn(pol), vy = varn(gel(fa,1)); /* vx < vy ! */
1059 :
1060 2583 : if (varncmp(vx,vy) <= 0) pari_err_BUG("twoembequation [incorrect variable priorities]");
1061 :
1062 2583 : lambda = shallowcopy(lambda);
1063 2583 : fa = shallowcopy(fa);
1064 2583 : j = 1;
1065 27839 : for (i=1; i<lfa; i++)
1066 25256 : if (signe(gel(lambda,i)))
1067 : {
1068 2618 : gel(lambda,j) = negi(gel(lambda,i));
1069 2618 : gel(fa,j) = gel(fa,i);
1070 2618 : j++;
1071 : }
1072 2583 : setlg(lambda, j);
1073 2583 : setlg(fa, j); lfa = j;
1074 :
1075 2583 : vpolx = ZXQ_powers(pol_x(vx),dx-1,pol);
1076 2583 : m = cgetg(dx+1, t_MAT);
1077 40292 : for (j=1; j <= dx; j++)
1078 37709 : gel(m,j) = cgetg(lfa, t_COL);
1079 5201 : for(i=1; i<lfa; i++)
1080 : {
1081 2618 : long dy = degpol(gel(fa,i));
1082 2618 : poly = pol_1(vy);
1083 40579 : for (j=1; j <= dx; j++)
1084 : {
1085 37961 : gcoeff(m,i,j) = RgXY_to_RgC(gadd(ZX_Z_mul(gel(vpolx,j),gel(lambda,i)),poly), dx, dy);
1086 37961 : poly = RgXQX_rem(RgX_shift(poly,1), gel(fa,i), pol);
1087 : }
1088 : }
1089 40292 : for(j=1; j<=dx; j++) gel(m,j) = shallowconcat1(gel(m,j));
1090 2583 : return QM_ker(m);
1091 : }
1092 :
1093 : static void
1094 1561 : subfields_cleanup(GEN* nf, GEN* pol, long* n, GEN* fa)
1095 : {
1096 1561 : *fa = NULL;
1097 1561 : if (typ(*nf) != t_VEC && typ(*nf) != t_POL) pari_err_TYPE("subfields_cleanup", *nf);
1098 1554 : if (typ(*nf) == t_VEC && lg(*nf) == 3)
1099 : {
1100 301 : *fa = gel(*nf,2);
1101 301 : *nf = gel(*nf,1);
1102 301 : if (typ(*fa)!=t_MAT || lg(*fa)!=3)
1103 14 : pari_err_TYPE("subfields_cleanup [fa should be a factorisation matrix]", *fa);
1104 : }
1105 1540 : if (typ(*nf) == t_POL)
1106 : {
1107 784 : *pol = *nf;
1108 784 : *nf = NULL;
1109 784 : if (!RgX_is_ZX(*pol)) pari_err_TYPE("subfields_cleanup [not integral]", *pol);
1110 777 : if (!equali1(leading_coeff(*pol))) pari_err_TYPE("subfields_cleanup [not monic]", *pol);
1111 770 : *n = degpol(*pol);
1112 770 : if (*n<=0) pari_err_TYPE("subfields_cleanup [constant polynomial]", *pol);
1113 : }
1114 : else
1115 : {
1116 756 : *nf = checknf(*nf);
1117 735 : *pol = nf_get_pol(*nf);
1118 735 : *n = degpol(*pol);
1119 : }
1120 1498 : if(*fa)
1121 : {
1122 273 : long v = varn(*pol);
1123 273 : GEN o = gcoeff(*fa,1,1);
1124 273 : if (varncmp(varn(o),v) >= 0) pari_err_PRIORITY("nfsubfields_fa", o, "<=", v);
1125 : }
1126 1477 : }
1127 :
1128 : static GEN
1129 280 : rootsuptoconj(GEN pol, long prec)
1130 : {
1131 : GEN ro;
1132 : long n, i;
1133 280 : ro = roots(pol,prec);
1134 280 : n = lg(ro)-1;
1135 1498 : for (i=1; i<=n/2; i++)
1136 1218 : gel(ro,i) = gel(ro,2*i-1);
1137 280 : setlg(ro,n/2+1);
1138 280 : return ro;
1139 : }
1140 : static GEN
1141 1085 : cmsubfield_get_roots(GEN pol, GEN nf, long n, long* r2, long *prec)
1142 : {
1143 : GEN ro;
1144 1085 : if (nf)
1145 : {
1146 735 : if (nf_get_r1(nf)) return NULL;
1147 371 : *r2 = nf_get_r2(nf);
1148 371 : *prec = nf_get_prec(nf);
1149 371 : ro = nf_get_roots(nf);
1150 : }
1151 : else
1152 : {
1153 350 : if (n%2 || sturm(pol)) return NULL;
1154 280 : *r2 = n/2;
1155 280 : *prec = MEDDEFAULTPREC;
1156 280 : ro = rootsuptoconj(pol, *prec);
1157 : }
1158 651 : return ro;
1159 : }
1160 :
1161 : static GEN
1162 595 : subfields_get_fa(GEN pol, GEN nf, GEN fa)
1163 : {
1164 595 : if (!fa)
1165 : {
1166 392 : GEN poly = shallowcopy(pol);
1167 392 : setvarn(poly, fetch_var_higher());
1168 392 : fa = nffactor(nf? nf: pol, poly);
1169 : }
1170 595 : return liftpol_shallow(gel(fa,1));
1171 : }
1172 :
1173 : static long
1174 287 : subfields_get_ero(GEN pol, GEN nf)
1175 : {
1176 574 : return 1 + gexpo(nf? nf_get_roots(nf):
1177 287 : QX_complex_roots(pol, LOWDEFAULTPREC));
1178 : }
1179 :
1180 : static GEN
1181 280 : try_imag(GEN x, GEN c, GEN pol, long v, ulong p, GEN emb, GEN galpol, long fl)
1182 : {
1183 280 : GEN a = Q_primpart(RgX_sub(RgX_RgXQ_eval(x,c,pol),x));
1184 280 : if (Flx_is_squarefree(Flxq_charpoly(ZX_to_Flx(a,p),ZX_to_Flx(pol,p),p),p))
1185 : {
1186 168 : pol = ZXQ_charpoly(a, pol, v);
1187 168 : return fl ? pol : mkvec2(pol, RgX_RgXQ_eval(a, emb, galpol));
1188 : }
1189 112 : return NULL;
1190 : }
1191 :
1192 : static GEN
1193 210 : galoissubfieldcm(GEN G, long fl)
1194 : {
1195 210 : pari_sp av = avma;
1196 : GEN c, H, elts, g, Hset, c2, gene, sub, pol, emb, a, galpol, B, b;
1197 : long n, i, j, nH, ind, v, d;
1198 210 : ulong p = 1009;
1199 :
1200 210 : galpol = gal_get_pol(G);
1201 210 : n = degpol(galpol);
1202 210 : v = varn(galpol);
1203 210 : c = galois_get_conj(G);
1204 : /* compute the list of c*g*c*g^(-1) : product of all pairs of conjugations
1205 : * maximal CM subfield is the field fixed by those elements, if c does not
1206 : * belong to the group they generate */
1207 210 : checkgroup(G, &elts);
1208 210 : elts = gen_sort_shallow(elts,(void*)vecsmall_lexcmp,cmp_nodata);
1209 210 : H = vecsmall_ei(n,1); /* indices of elements of H */
1210 210 : Hset = zero_F2v(n);
1211 210 : F2v_set(Hset,1);
1212 210 : nH = 1;
1213 1456 : for (i=2; i<=n; i++)
1214 : {
1215 1246 : g = gel(elts,i);
1216 1246 : c2 = perm_mul(c,perm_conj(g,c));
1217 1246 : if (!F2v_coeff(Hset,c2[1]))
1218 : {
1219 182 : nH++;
1220 182 : H[nH] = c2[1];
1221 182 : F2v_set(Hset,c2[1]);
1222 : }
1223 : }
1224 : /* group generated */
1225 210 : gene = gcopy(H);
1226 210 : setlg(gene,nH+1);
1227 210 : i = 1; /* last element that has been multiplied by the generators */
1228 392 : while (i < nH)
1229 : {
1230 1218 : for (j=1; j<lg(gene); j++)
1231 : {
1232 1036 : g = gel(elts,gene[j]);
1233 1036 : ind = g[H[i]]; /* index of the product */
1234 1036 : if (!F2v_coeff(Hset,ind))
1235 : {
1236 0 : nH++;
1237 0 : if (ind==c[1] || 2*nH>n) return gc_const(av, gen_0);
1238 0 : H[nH] = ind;
1239 0 : F2v_set(Hset,ind);
1240 : }
1241 : }
1242 182 : i++;
1243 : }
1244 210 : H = cgetg(lg(gene), t_VEC);
1245 602 : for (i=1; i<lg(H); i++)
1246 392 : gel(H,i) = gel(elts,gene[i]);
1247 210 : sub = galoisfixedfield(G, H, 0, -1);
1248 :
1249 : /* compute a totally imaginary generator */
1250 210 : pol = gel(sub,1);
1251 210 : emb = liftpol_shallow(gel(sub,2));
1252 210 : d = degpol(pol);
1253 210 : if (!(ZX_deflate_order(pol)%2) && sturm(RgX_deflate(pol,2))==d/2)
1254 : {
1255 42 : setvarn(pol,v);
1256 42 : return fl==1 ? pol: mkvec2(pol,emb);
1257 : }
1258 :
1259 : /* compute action of c on the subfield from that on the large field */
1260 168 : c = galoispermtopol(G,c);
1261 168 : if (d<n)
1262 : {
1263 35 : GEN M = cgetg(d+1,t_MAT), contc, contM;
1264 35 : gel(M,1) = col_ei(n,1); a = pol_1(v);
1265 98 : for (i=2; i<=d; i++)
1266 : {
1267 63 : a = RgX_rem(QX_mul(a,emb), galpol);
1268 63 : gel(M,i) = RgX_to_RgC(a,n);
1269 : }
1270 35 : c = RgX_RgXQ_eval(emb,c,galpol);
1271 35 : c = Q_primitive_part(c,&contc);
1272 35 : c = RgX_to_RgC(c,n);
1273 35 : M = Q_primitive_part(M,&contM);
1274 35 : c = RgM_RgC_invimage(M,c);
1275 35 : if (contc)
1276 : {
1277 21 : if (contM) contc = gdiv(contc,contM);
1278 21 : c = RgV_Rg_mul(c, contc);
1279 : }
1280 14 : else if (contM) c = RgV_Rg_mul(c, ginv(contM));
1281 35 : c = RgV_to_RgX(c, v);
1282 : }
1283 :
1284 : /* search for a generator of the form c(b)-b */
1285 273 : for (i=1; i<d; i++)
1286 : {
1287 238 : a = try_imag(pol_xn(i,v),c,pol,v,p,emb,galpol,fl);
1288 238 : if (a) return a;
1289 105 : p = unextprime(p+1);
1290 : }
1291 35 : B = stoi(10);
1292 35 : b = pol_xn(d-1,v);
1293 : while(1)
1294 : {
1295 210 : for (i=2; i<lg(b); i++) gel(b,i) = randomi(B);
1296 42 : a = try_imag(b,c,pol,v,p,emb,galpol,fl);
1297 42 : if (a) return a;
1298 7 : p = unextprime(p+1);
1299 : }
1300 : return NULL;/*LCOV_EXCL_LINE*/
1301 : }
1302 :
1303 : static GEN
1304 133 : quadsubfieldcm(GEN pol, long fl)
1305 : {
1306 133 : GEN a = gel(pol,3), b = gel(pol,2), d, P;
1307 133 : long v = varn(pol);
1308 133 : if (mpodd(a))
1309 35 : { b = mului(4, b); d = gen_2; }
1310 : else
1311 98 : { a = divis(a,2); d = gen_1; }
1312 133 : P = deg2pol_shallow(gen_1, gen_0, subii(b, sqri(a)), v);
1313 133 : return fl==1 ? P: mkvec2(P, deg1pol_shallow(d,a,v));
1314 : }
1315 :
1316 : GEN
1317 1148 : nfsubfieldscm(GEN nf, long fl)
1318 : {
1319 1148 : pari_sp av = avma;
1320 : GEN fa, lambda, V, res, ro, a, aa, ev, minev, pol, G;
1321 1148 : long i, j, n, r2, minj=0, prec, emax, emin, e, precbound, ero;
1322 :
1323 1148 : subfields_cleanup(&nf, &pol, &n, &fa);
1324 1085 : ro = cmsubfield_get_roots(pol, nf, n, &r2, &prec);
1325 1085 : if (!ro) return gc_const(av, gen_0);
1326 : /* now r2 == 2*n */
1327 :
1328 651 : if (n==2) return gerepilecopy(av, quadsubfieldcm(pol, fl));
1329 518 : G = galoisinit(nf? nf: pol, NULL);
1330 518 : if (G != gen_0) return gerepilecopy(av, galoissubfieldcm(G, fl));
1331 :
1332 308 : ero = 0;
1333 1624 : for (i=1; i<lg(ro); i++)
1334 : {
1335 1316 : e = 1+gexpo(gel(ro,i));
1336 1316 : if (e > ero) ero = e;
1337 : }
1338 308 : ero++;
1339 308 : fa = subfields_get_fa(pol, nf, fa);
1340 :
1341 308 : emax = 1;
1342 308 : emin = -1;
1343 1624 : for (i=1; i<lg(ro); i++)
1344 4963 : for (j=i+1; j<lg(ro); j++)
1345 : {
1346 3647 : e = gexpo(gsub(gel(ro,i),gel(ro,j)));
1347 3647 : if (e > emax) emax = e;
1348 3647 : if (e < emin) emin = e;
1349 : }
1350 308 : precbound = n*(emax-emin) + gexpo(fa) + n*n + 5;
1351 308 : precbound = 3 + precbound/BITS_IN_LONG;
1352 308 : if (prec < precbound)
1353 : {
1354 0 : prec = precbound;
1355 0 : ro = rootsuptoconj(pol, prec);
1356 : }
1357 :
1358 308 : lambda = zerovec(lg(fa)-1);
1359 1624 : for (i=1; i<=r2; i++)
1360 : {
1361 1316 : a = gel(ro,i);
1362 1316 : aa = conj_i(a);
1363 9422 : for (j=1; j<lg(fa); j++)
1364 : {
1365 8106 : ev = cxnorm(poleval(poleval(gel(fa,j),aa),a));
1366 8106 : if (j==1 || cmprr(minev,ev)>0) { minj = j; minev = ev; }
1367 : }
1368 1316 : gel(lambda,minj) = gen_m1;
1369 : }
1370 :
1371 308 : V = twoembequation(pol, fa, lambda);
1372 308 : if (lg(V)==1) { delete_var(); return gc_const(av, gen_0); }
1373 259 : res = subfield_generator(pol, V, 2*(lg(V)-1), ero, fl);
1374 259 : delete_var();
1375 259 : return gerepilecopy(av, res);
1376 : }
1377 :
1378 : static int
1379 19474 : field_is_contained(GEN V, GEN W, int strict)
1380 : {
1381 : GEN VW;
1382 19474 : ulong p = 1073741827;
1383 : /* distinct overfield must have different dimension */
1384 19474 : if (strict && lg(V) == lg(W)) return 0;
1385 : /* dimension of overfield must be multiple */
1386 14553 : if ((lg(W)-1) % (lg(V)-1)) return 0;
1387 10402 : VW = shallowconcat(V,W);
1388 10402 : if (Flm_rank(ZM_to_Flm(VW,p),p) > lg(W)-1) return 0;
1389 4235 : return ZM_rank(VW) == lg(W)-1;
1390 : }
1391 :
1392 : /***********************************************/
1393 : /* */
1394 : /* Maximal, generating, all subfields */
1395 : /* Aurel Page (2019) */
1396 : /* after van Hoeij, Klueners, Novocin */
1397 : /* Journal of Symbolic Computation 52 (2013) */
1398 : /* */
1399 : /***********************************************/
1400 :
1401 : const long subf_MAXIMAL = 1; /* return the maximal subfields */
1402 : const long subf_GENERATING = 2; /* return the generating subfields */
1403 : static GEN
1404 287 : maxgen_subfields(GEN pol, GEN fa, long flag)
1405 : {
1406 287 : pari_sp av = avma;
1407 287 : GEN principal, ismax, isgene, Lmax = NULL, Lgene, res, V, W, W1;
1408 287 : long i, i2, j, flmax, flgene, nbmax = 0, nbgene = 0;
1409 :
1410 287 : if (!flag) return cgetg(1,t_VEC);
1411 287 : flmax = (flag & subf_MAXIMAL)!=0;
1412 287 : flgene = (flag & subf_GENERATING)!=0;
1413 :
1414 : /* compute principal subfields */
1415 287 : principal = cgetg(lg(fa),t_VEC);
1416 2562 : for (i=1; i<lg(fa); i++)
1417 2275 : gel(principal,i) = twoembequation(pol, fa, vec_ei(lg(fa)-1,i));
1418 287 : principal = gen_sort_uniq(principal, (void*)&cmp_universal, &cmp_nodata);
1419 : /* remove nf and duplicates (sort_uniq possibly not enough) */
1420 287 : i2 = 1;
1421 1694 : for (i=1; i<lg(principal)-1; i++)
1422 : {
1423 1407 : long dup = 0;
1424 1407 : V = gel(principal,i);
1425 1407 : j = i2-1;
1426 2877 : while (j > 0 && lg(gel(principal,j)) == lg(V))
1427 : {
1428 1470 : if (field_is_contained(gel(principal,j),V,0)) { dup=1; break; }
1429 1470 : j--;
1430 : }
1431 1407 : if (!dup) gel(principal,i2++) = V;
1432 : }
1433 287 : setlg(principal, i2);
1434 :
1435 : /* a subfield is generating iff all overfields contain the first overfield */
1436 287 : ismax = cgetg(lg(principal),t_VECSMALL);
1437 287 : isgene = cgetg(lg(principal),t_VECSMALL);
1438 1694 : for (i=1; i<lg(principal); i++)
1439 : {
1440 1407 : V = gel(principal,i);
1441 1407 : ismax[i] = flmax;
1442 1407 : isgene[i] = flgene;
1443 1407 : W1 = NULL; /* intersection of strict overfields */
1444 4858 : for (j=i+1; j<lg(principal); j++)
1445 : {
1446 3696 : W = gel(principal,j);
1447 3696 : if (!field_is_contained(V,W,1)) continue;
1448 693 : ismax[i] = 0;
1449 693 : if (!flgene) break;
1450 483 : if (!W1) { W1 = W; continue; }
1451 189 : if (!field_is_contained(W1,W,1))
1452 : {
1453 63 : W1 = intersect(W1,W);
1454 63 : if (lg(W1)==lg(V)) { isgene[i]=0; break; }
1455 : }
1456 : }
1457 : }
1458 :
1459 1694 : for (i=1; i<lg(principal); i++)
1460 : {
1461 1407 : nbmax += ismax[i];
1462 1407 : nbgene += isgene[i];
1463 : }
1464 :
1465 287 : if (flmax)
1466 : {
1467 98 : Lmax = cgetg(nbmax+1, t_VEC);
1468 98 : j=1;
1469 518 : for (i=1; i<lg(principal); i++)
1470 420 : if (ismax[i]) gel(Lmax,j++) = gel(principal,i);
1471 : }
1472 :
1473 287 : if (flgene)
1474 : {
1475 189 : Lgene = cgetg(nbgene+1, t_VEC);
1476 189 : j=1;
1477 1176 : for (i=1; i<lg(principal); i++)
1478 987 : if (isgene[i]) gel(Lgene,j++) = gel(principal,i);
1479 : }
1480 :
1481 287 : if (!flgene) res = Lmax;
1482 189 : else if (!flmax) res = Lgene;
1483 0 : else res = mkvec2(Lmax,Lgene);
1484 287 : return gerepilecopy(av, res);
1485 : }
1486 :
1487 : GEN
1488 154 : nfsubfieldsmax(GEN nf, long fl)
1489 : {
1490 154 : pari_sp av = avma;
1491 : GEN pol, fa, Lmax, V;
1492 : long n, i, ero;
1493 :
1494 154 : subfields_cleanup(&nf, &pol, &n, &fa);
1495 154 : if (n==1) { set_avma(av); return cgetg(1,t_VEC); }
1496 140 : if (uisprime(n))
1497 63 : return gerepilecopy(av, fl==1 ? mkvec(pol_x(varn(pol)))
1498 21 : : mkvec(mkvec2(pol_x(varn(pol)),gen_0)));
1499 98 : ero = subfields_get_ero(pol, nf);
1500 98 : fa = subfields_get_fa(pol, nf, fa);
1501 98 : Lmax = maxgen_subfields(pol, fa, subf_MAXIMAL);
1502 308 : for (i=1; i<lg(Lmax); i++)
1503 : {
1504 210 : V = gel(Lmax,i);
1505 210 : gel(Lmax,i) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1506 : }
1507 98 : delete_var();
1508 98 : return gerepilecopy(av, Lmax);
1509 : }
1510 :
1511 : static void
1512 1764 : heap_climb(GEN* H, long i)
1513 : {
1514 : long j;
1515 1764 : if (i==1) return;
1516 1302 : j = i/2;
1517 1302 : if (cmp_universal(gel(*H,i),gel(*H,j)) > 0)
1518 : {
1519 532 : swap(gel(*H,i), gel(*H,j));
1520 532 : return heap_climb(H,j);
1521 : }
1522 : }
1523 :
1524 : static void
1525 1232 : heap_push(GEN* H, long *len, GEN x)
1526 : {
1527 1232 : if (*len+1 == lg(*H))
1528 : {
1529 14 : GEN H2 = zerovec(2*(*len));
1530 : long i;
1531 154 : for(i=1; i<lg(*H); i++)
1532 140 : gel(H2,i) = gel(*H,i);
1533 14 : *H = H2;
1534 : }
1535 1232 : (*len)++;
1536 1232 : gel(*H,*len) = x;
1537 1232 : return heap_climb(H,*len);
1538 : }
1539 :
1540 : static void
1541 2569 : heap_descend(GEN* H, long len, long i)
1542 : {
1543 2569 : long maxi = i, j = 2*i;
1544 2569 : if (j > len) return;
1545 1337 : if (cmp_universal(gel(*H,j),gel(*H,i)) > 0) maxi = j;
1546 1337 : j++;
1547 1337 : if (j<=len && cmp_universal(gel(*H,j),gel(*H,maxi))>0) maxi = j;
1548 1337 : if (maxi == i) return;
1549 1148 : swap(gel(*H,i), gel(*H,maxi));
1550 1148 : return heap_descend(H,len,maxi);
1551 : }
1552 :
1553 : static void
1554 1421 : heap_pop(GEN *H, long *len, GEN* top)
1555 : {
1556 1421 : *top = gel(*H,1);
1557 1421 : gel(*H,1) = gel(*H,*len);
1558 1421 : (*len)--;
1559 1421 : return heap_descend(H,*len,1);
1560 : };
1561 :
1562 : static GEN
1563 259 : nfsubfields_fa(GEN nf, long d, long fl)
1564 : {
1565 259 : pari_sp av = avma;
1566 : GEN pol, fa, gene, res, res2, H, V, v, W, w, data;
1567 : long n, r, i, j, nres, len, s, newfield, ero, vp;
1568 :
1569 259 : subfields_cleanup(&nf, &pol, &n, &fa); vp = varn(pol);
1570 238 : if (d && (d<1 || d>n || n%d)) return gerepilecopy(av, cgetg(1,t_VEC));
1571 245 : if (!d && uisprime(n)) return gerepilecopy(av,
1572 0 : fl==1 ? mkvec2( pol_x(varn(pol)), pol)
1573 14 : : mkvec2( mkvec2(pol_x(vp),pol_0(vp)), mkvec2(pol,pol_x(vp))));
1574 231 : if (n==1 || d==1) return gerepilecopy(av,
1575 14 : fl==1 ? mkvec(pol_x(varn(pol))): _subfield(pol_x(vp),pol_0(vp)));
1576 217 : if (d==n) return gerepilecopy(av,
1577 14 : fl==1 ? mkvec(pol): _subfield(pol,pol_x(vp)));
1578 189 : ero = subfields_get_ero(pol, nf);
1579 189 : fa = subfields_get_fa(pol, nf, fa);
1580 189 : gene = maxgen_subfields(pol, fa, subf_GENERATING);
1581 :
1582 189 : if (d)
1583 : {
1584 : /* keep only generating subfields of degree a multiple of d */
1585 14 : j=1;
1586 147 : for (i=1; i<lg(gene); i++)
1587 133 : if ((lg(gel(gene,i))-1) % d == 0)
1588 : {
1589 98 : gel(gene,j) = gel(gene,i);
1590 98 : j++;
1591 : }
1592 14 : setlg(gene,j);
1593 : }
1594 189 : r = lg(gene)-1;
1595 :
1596 189 : res = zerovec(10);
1597 189 : nres = 0;
1598 189 : H = zerovec(10);
1599 189 : gel(H,1) = mkvec3(matid(n),zero_F2v(r),mkvecsmall(0));
1600 189 : len = 1;
1601 :
1602 1610 : while (len>0)
1603 : {
1604 1421 : heap_pop(&H, &len, &data);
1605 1421 : V = gel(data,1);
1606 1421 : v = gel(data,2);
1607 1421 : s = gel(data,3)[1];
1608 6153 : for (i=s+1; i<=r; i++)
1609 4732 : if (!F2v_coeff(v,i))
1610 : {
1611 3675 : W = vec_Q_primpart(intersect(V, gel(gene,i)));
1612 3675 : w = F2v_copy(v);
1613 3675 : F2v_set(w, i);
1614 3675 : newfield = 1;
1615 18130 : for (j=1; j<=r; j++)
1616 16800 : if (!F2v_coeff(w,j) && field_is_contained(W,gel(gene,j),1))
1617 : {
1618 3416 : if (j<i) { newfield = 0; break; }
1619 1071 : F2v_set(w,j);
1620 : }
1621 3675 : if (newfield && (!d || (lg(W)-1)%d==0)) heap_push(&H, &len, mkvec3(W,w,mkvecsmall(i)));
1622 : }
1623 :
1624 1421 : if (!d || lg(V)-1==d)
1625 : {
1626 1407 : nres++;
1627 1407 : if (nres == lg(res))
1628 : {
1629 28 : res2 = zerovec(2*lg(res));
1630 392 : for(j=1; j<lg(res); j++) gel(res2,j) = gel(res,j);
1631 28 : res = res2;
1632 : }
1633 1407 : gel(res,nres) = subfield_generator(pol, V, lg(V)-1, ero, fl);
1634 : }
1635 : }
1636 189 : setlg(res,nres+1);
1637 189 : vecreverse_inplace(res);
1638 :
1639 189 : delete_var();
1640 189 : return gerepilecopy(av, res);
1641 : }
|
