Line data Source code
1 : /* Copyright (C) 2000 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 : /* NUMBER FIELDS */
18 : /* */
19 : /**************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_nf
24 :
25 : int new_galois_format = 0;
26 :
27 : /* v a t_VEC, lg(v) = 13, sanity check for true rnf */
28 : static int
29 298609 : v13checkrnf(GEN v)
30 298609 : { return typ(gel(v,6)) == t_VEC; }
31 : static int
32 46417 : rawcheckbnf(GEN v) { return typ(v)==t_VEC && lg(v)==11; }
33 : static int
34 45458 : rawchecknf(GEN v) { return typ(v)==t_VEC && lg(v)==10; }
35 : /* v a t_VEC, lg(v) = 11, sanity check for true bnf */
36 : static int
37 37254 : v11checkbnf(GEN v) { return rawchecknf(bnf_get_nf(v)); }
38 : /* v a t_VEC, lg(v) = 13, sanity check for true nf and true bnf */
39 : static int
40 819 : v13checkgchar(GEN v) { return rawchecknf(gchar_get_nf(v)) && rawcheckbnf(gchar_get_bnf(v)); }
41 : /* v a t_VEC, lg(v) = 10, sanity check for true nf */
42 : static int
43 75894 : v10checknf(GEN v) { return typ(gel(v,1))==t_POL; }
44 : /* v a t_VEC, lg(v) = 9, sanity check for true gal */
45 : static int
46 658 : v9checkgal(GEN v)
47 658 : { GEN x = gel(v,2); return typ(x) == t_VEC && lg(x) == 4; }
48 :
49 : int
50 308076 : checkrnf_i(GEN rnf)
51 308076 : { return (typ(rnf)==t_VEC && lg(rnf)==13 && v13checkrnf(rnf)); }
52 :
53 : void
54 298007 : checkrnf(GEN rnf)
55 298007 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
56 :
57 : GEN
58 5239470 : checkbnf_i(GEN X)
59 : {
60 5239470 : if (typ(X) == t_VEC)
61 5239296 : switch (lg(X))
62 : {
63 5235562 : case 11:
64 5235562 : if (typ(gel(X,6)) != t_INT) return NULL; /* pre-2.2.4 format */
65 5235562 : if (lg(gel(X,10)) != 4) return NULL; /* pre-2.8.1 format */
66 5235562 : return X;
67 3262 : case 7: return checkbnf_i(bnr_get_bnf(X));
68 : }
69 646 : return NULL;
70 : }
71 :
72 : GEN
73 78187444 : checknf_i(GEN X)
74 : {
75 78187444 : if (typ(X)==t_VEC)
76 77538404 : switch(lg(X))
77 : {
78 77056249 : case 10: return X;
79 475321 : case 11: return checknf_i(bnf_get_nf(X));
80 2884 : case 7: return checknf_i(bnr_get_bnf(X));
81 2212 : case 3: if (typ(gel(X,2)) == t_POLMOD) return checknf_i(gel(X,1));
82 : }
83 652983 : return NULL;
84 : }
85 :
86 : GEN
87 3231402 : checkbnf(GEN x)
88 : {
89 3231402 : GEN bnf = checkbnf_i(x);
90 3231405 : if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
91 3231385 : return bnf;
92 : }
93 :
94 : GEN
95 76214530 : checknf(GEN x)
96 : {
97 76214530 : GEN nf = checknf_i(x);
98 76214276 : if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
99 76214248 : return nf;
100 : }
101 :
102 : GEN
103 2003161 : checkbnr_i(GEN bnr)
104 : {
105 2003161 : if (typ(bnr)!=t_VEC || lg(bnr)!=7 || !checkbnf_i(bnr_get_bnf(bnr)))
106 84 : return NULL;
107 2003067 : return bnr;
108 : }
109 : void
110 1820081 : checkbnr(GEN bnr)
111 : {
112 1820081 : if (!checkbnr_i(bnr))
113 14 : pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
114 1820053 : }
115 :
116 : void
117 0 : checksqmat(GEN x, long N)
118 : {
119 0 : if (typ(x)!=t_MAT) pari_err_TYPE("checksqmat",x);
120 0 : if (lg(x) == 1 || lgcols(x) != N+1) pari_err_DIM("checksqmat");
121 0 : }
122 :
123 : GEN
124 1636656 : checkbid_i(GEN bid)
125 : {
126 : GEN f;
127 1636656 : if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(bid_get_U(bid)) != t_VEC)
128 254202 : return NULL;
129 1382453 : f = bid_get_mod(bid);
130 1382457 : if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
131 1382457 : return bid;
132 : }
133 : void
134 1379896 : checkbid(GEN bid)
135 : {
136 1379896 : if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
137 1379888 : }
138 : void
139 19810 : checkabgrp(GEN v)
140 : {
141 19810 : if (typ(v) == t_VEC) switch(lg(v))
142 : {
143 17423 : case 4: if (typ(gel(v,3)) != t_VEC) break;
144 19810 : case 3: if (typ(gel(v,2)) != t_VEC) break;
145 19782 : if (typ(gel(v,1)) != t_INT) break;
146 19782 : return;/*OK*/
147 0 : default: break;
148 : }
149 28 : pari_err_TYPE("checkabgrp",v);
150 : }
151 :
152 : GEN
153 270635 : checknfelt_mod(GEN nf, GEN x, const char *s)
154 : {
155 270635 : GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
156 270638 : if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
157 270526 : return a;
158 : }
159 :
160 : int
161 10575 : check_ZKmodule_i(GEN M)
162 : {
163 10575 : return (typ(M) ==t_VEC && lg(M) >= 3
164 10575 : && typ(gel(M,1)) == t_MAT
165 10575 : && typ(gel(M,2)) == t_VEC
166 21150 : && lgcols(M) == lg(gel(M,2)));
167 : }
168 : void
169 10519 : check_ZKmodule(GEN M, const char *s)
170 10519 : { if (!check_ZKmodule_i(M)) pari_err_TYPE(s, M); }
171 :
172 : static long
173 113141 : typv6(GEN x)
174 : {
175 113141 : if (typ(gel(x,1)) == t_VEC && lg(gel(x,3)) == 3)
176 : {
177 12705 : GEN t = gel(x,3);
178 12705 : if (typ(t) != t_VEC) return typ_NULL;
179 12705 : t = gel(x,5);
180 12705 : switch(typ(gel(x,5)))
181 : {
182 448 : case t_VEC: return typ_BID;
183 12257 : case t_MAT: return typ_BIDZ;
184 0 : default: return typ_NULL;
185 : }
186 : }
187 100436 : if (typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT) return typ_PRID;
188 2527 : if (typ(gel(x,1)) == t_INT && typ(gel(x,2)) == t_VEC) return typ_QUA;
189 196 : return typ_NULL;
190 : }
191 :
192 : GEN
193 22092 : get_bnf(GEN x, long *t)
194 : {
195 22092 : switch(typ(x))
196 : {
197 56 : case t_POL: *t = typ_POL; return NULL;
198 56 : case t_QUAD: *t = typ_Q ; return NULL;
199 21469 : case t_VEC:
200 21469 : switch(lg(x))
201 : {
202 2247 : case 5: if (typ(gel(x,1)) != t_INT) break;
203 2191 : *t = typ_QUA; return NULL;
204 8463 : case 6: *t = typv6(x); return NULL;
205 2373 : case 7: *t = typ_BNR;
206 2373 : x = bnr_get_bnf(x);
207 2373 : if (!rawcheckbnf(x)) break;
208 2317 : return x;
209 84 : case 9:
210 84 : if (!v9checkgal(x)) break;
211 84 : *t = typ_GAL; return NULL;
212 462 : case 10:
213 462 : if (!v10checknf(x)) break;
214 462 : *t = typ_NF; return NULL;
215 469 : case 11:
216 469 : if (!v11checkbnf(x)) break;
217 469 : *t = typ_BNF; return x;
218 245 : case 13:
219 245 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_bnf(x); }
220 56 : if (!v13checkrnf(x)) break;
221 56 : *t = typ_RNF; return NULL;
222 266 : case 17: *t = typ_ELL; return NULL;
223 : }
224 6972 : break;
225 112 : case t_COL:
226 112 : if (get_prid(x)) { *t = typ_MODPR; return NULL; }
227 56 : break;
228 : }
229 7427 : *t = typ_NULL; return NULL;
230 : }
231 :
232 : GEN
233 115535 : get_nf(GEN x, long *t)
234 : {
235 115535 : switch(typ(x))
236 : {
237 133 : case t_POL : *t = typ_POL; return NULL;
238 133 : case t_QUAD: *t = typ_Q ; return NULL;
239 112588 : case t_VEC:
240 112588 : switch(lg(x))
241 : {
242 140 : case 3:
243 140 : if (typ(gel(x,2)) != t_POLMOD) break;
244 140 : return get_nf(gel(x,1),t);
245 133 : case 5:
246 133 : if (typ(gel(x,1)) != t_INT) break;
247 0 : *t = typ_QUA; return NULL;
248 99078 : case 6: *t = typv6(x); return NULL;
249 7301 : case 7:
250 7301 : x = bnr_get_bnf(x);
251 7301 : if (!rawcheckbnf(x) || !rawchecknf(x = bnf_get_nf(x))) break;
252 7168 : *t = typ_BNR; return x;
253 567 : case 9:
254 567 : if (!v9checkgal(x)) break;
255 567 : *t = typ_GAL; return NULL;
256 1141 : case 10:
257 1141 : if (!v10checknf(x)) break;
258 1141 : *t = typ_NF; return x;
259 217 : case 11:
260 217 : if (!rawchecknf(x = bnf_get_nf(x))) break;
261 217 : *t = typ_BNF; return x;
262 392 : case 13:
263 392 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_nf(x); }
264 364 : if (!v13checkrnf(x)) break;
265 364 : *t = typ_RNF; return NULL;
266 3486 : case 17: *t = typ_ELL; return NULL;
267 : }
268 399 : break;
269 133 : case t_QFB: *t = typ_QFB; return NULL;
270 266 : case t_COL:
271 266 : if (get_prid(x)) { *t = typ_MODPR; return NULL; }
272 133 : break;
273 : }
274 2814 : *t = typ_NULL; return NULL;
275 : }
276 :
277 : long
278 116998 : nftyp(GEN x)
279 : {
280 116998 : switch(typ(x))
281 : {
282 28 : case t_POL : return typ_POL;
283 7 : case t_QUAD: return typ_Q;
284 116956 : case t_VEC:
285 116956 : switch(lg(x))
286 : {
287 182 : case 13:
288 182 : if (v13checkgchar(x)) return typ_GCHAR;
289 175 : if (!v13checkrnf(x)) break;
290 175 : return typ_RNF;
291 74291 : case 10:
292 74291 : if (!v10checknf(x)) break;
293 74284 : return typ_NF;
294 280 : case 11:
295 280 : if (!v11checkbnf(x)) break;
296 280 : return typ_BNF;
297 36519 : case 7:
298 36519 : x = bnr_get_bnf(x);
299 36519 : if (!rawcheckbnf(x) || !v11checkbnf(x)) break;
300 36505 : return typ_BNR;
301 5600 : case 6:
302 5600 : return typv6(x);
303 7 : case 9:
304 7 : if (!v9checkgal(x)) break;
305 0 : return typ_GAL;
306 7 : case 17: return typ_ELL;
307 : }
308 14 : }
309 105 : return typ_NULL;
310 : }
311 :
312 : /*************************************************************************/
313 : /** **/
314 : /** GALOIS GROUP **/
315 : /** **/
316 : /*************************************************************************/
317 :
318 : GEN
319 7654 : tschirnhaus(GEN x)
320 : {
321 7654 : pari_sp av = avma, av2;
322 7654 : long a, v = varn(x);
323 7654 : GEN u, y = cgetg(5,t_POL);
324 :
325 7654 : if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
326 7654 : if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
327 7654 : if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
328 7654 : y[1] = evalsigne(1)|evalvarn(0);
329 : do
330 : {
331 8326 : a = random_bits(2); if (a==0) a = 1; gel(y,4) = stoi(a);
332 8326 : a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
333 8326 : a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
334 8326 : u = RgXQ_charpoly(y,x,v); av2 = avma;
335 : }
336 8326 : while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
337 7654 : if (DEBUGLEVEL>1)
338 0 : err_printf("Tschirnhaus transform. New pol: %Ps",u);
339 7654 : set_avma(av2); return gerepileupto(av,u);
340 : }
341 :
342 : /* Assume pol in Z[X], monic of degree n. Find L in Z such that
343 : * POL = L^(-n) pol(L x) is monic in Z[X]. Return POL and set *ptk = L.
344 : * No GC. */
345 : GEN
346 843817 : ZX_Z_normalize(GEN pol, GEN *ptk)
347 : {
348 843817 : long i,j, sk, n = degpol(pol); /* > 0 */
349 : GEN k, fa, P, E, a, POL;
350 :
351 843817 : if (ptk) *ptk = gen_1;
352 843817 : if (!n) return pol;
353 843810 : a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
354 1622949 : for (i = n-2; i >= 0; i--)
355 : {
356 1367318 : k = gcdii(k, gel(a,i));
357 1367119 : if (is_pm1(k)) return pol;
358 : }
359 255631 : sk = signe(k);
360 255631 : if (!sk) return pol; /* monomial! */
361 253699 : fa = absZ_factor_limit(k, 0); k = gen_1;
362 253718 : P = gel(fa,1);
363 253718 : E = gel(fa,2);
364 253718 : POL = leafcopy(pol); a = POL+2;
365 544906 : for (i = lg(P)-1; i > 0; i--)
366 : {
367 291180 : GEN p = gel(P,i), pv, pvj;
368 291180 : long vmin = itos(gel(E,i));
369 : /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
370 1203870 : for (j=n-1; j>=0; j--)
371 : {
372 : long v;
373 912686 : if (!signe(gel(a,j))) continue;
374 784603 : v = Z_pval(gel(a,j), p) / (n - j);
375 784605 : if (v < vmin) vmin = v;
376 : }
377 291184 : if (!vmin) continue;
378 18263 : pvj = pv = powiu(p,vmin); k = mulii(k, pv);
379 : /* a[j] /= p^(v*(n-j)) */
380 74164 : for (j=n-1; j>=0; j--)
381 : {
382 55903 : if (j < n-1) pvj = mulii(pvj, pv);
383 55903 : gel(a,j) = diviiexact(gel(a,j), pvj);
384 : }
385 : }
386 253726 : if (ptk) *ptk = k;
387 253726 : return POL;
388 : }
389 :
390 : /* Assume pol != 0 in Z[X]. Find C in Q, L in Z such that POL = C pol(x/L) monic
391 : * in Z[X]. Return POL and set *pL = L. Wasteful (but correct) if pol is not
392 : * primitive: better if caller used Q_primpart already. No GC. */
393 : GEN
394 841188 : ZX_primitive_to_monic(GEN pol, GEN *pL)
395 : {
396 841188 : long i,j, n = degpol(pol);
397 841188 : GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
398 :
399 841193 : if (is_pm1(lc))
400 : {
401 838309 : if (pL) *pL = gen_1;
402 838309 : return signe(lc) < 0? ZX_neg(pol): pol;
403 : }
404 2884 : if (signe(lc) < 0)
405 35 : POL = ZX_neg(pol);
406 : else
407 2849 : POL = leafcopy(pol);
408 2884 : a = POL+2; lc = gel(a,n);
409 2884 : fa = absZ_factor_limit(lc,0); L = gen_1;
410 2884 : P = gel(fa,1);
411 2884 : E = gel(fa,2);
412 5936 : for (i = lg(P)-1; i > 0; i--)
413 : {
414 3052 : GEN p = gel(P,i), pk, pku;
415 3052 : long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
416 :
417 3052 : if (d < 0) { k++; d += n; }
418 : /* k = ceil(e[i] / n); find d, k such that p^d pol(x / p^k) monic */
419 17171 : for (j=n-1; j>0; j--)
420 : {
421 14119 : if (!signe(gel(a,j))) continue;
422 8722 : v = Z_pval(gel(a,j), p);
423 8799 : while (v + d < k * j) { k++; d += n; }
424 : }
425 3052 : pk = powiu(p,k); j0 = d/k;
426 3052 : L = mulii(L, pk);
427 :
428 3052 : pku = powiu(p,d - k*j0);
429 : /* a[j] *= p^(d - kj) */
430 7175 : for (j=j0; j>=0; j--)
431 : {
432 4123 : if (j < j0) pku = mulii(pku, pk);
433 4123 : gel(a,j) = mulii(gel(a,j), pku);
434 : }
435 3052 : j0++;
436 3052 : pku = powiu(p,k*j0 - d);
437 : /* a[j] /= p^(kj - d) */
438 19152 : for (j=j0; j<=n; j++)
439 : {
440 16100 : if (j > j0) pku = mulii(pku, pk);
441 16100 : gel(a,j) = diviiexact(gel(a,j), pku);
442 : }
443 : }
444 2884 : if (pL) *pL = L;
445 2884 : return POL;
446 : }
447 : /* Assume pol != 0 in Z[X]. Find C,L in Q such that POL = C pol(x/L)
448 : * monic in Z[X]. Return POL and set *pL = L.
449 : * Wasteful (but correct) if pol is not primitive: better if caller used
450 : * Q_primpart already. No GC. */
451 : GEN
452 840871 : ZX_Q_normalize(GEN pol, GEN *pL)
453 : {
454 840871 : GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
455 840877 : POL = ZX_Z_normalize(POL, pL);
456 840809 : if (pL) *pL = gdiv(lc, *pL);
457 840898 : return POL;
458 : }
459 :
460 : GEN
461 5803677 : ZX_Q_mul(GEN A, GEN z)
462 : {
463 5803677 : pari_sp av = avma;
464 5803677 : long i, l = lg(A);
465 : GEN d, n, Ad, B, u;
466 5803677 : if (typ(z)==t_INT) return ZX_Z_mul(A,z);
467 4872176 : n = gel(z, 1); d = gel(z, 2);
468 4872176 : Ad = RgX_to_RgC(FpX_red(A, d), l-2);
469 4872176 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
470 4872170 : B = cgetg(l, t_POL);
471 4872177 : B[1] = A[1];
472 4872177 : if (equali1(u))
473 : {
474 4689610 : for(i=2; i<l; i++)
475 2900492 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
476 : } else
477 : {
478 18590288 : for(i=2; i<l; i++)
479 : {
480 15507233 : GEN di = gcdii(gel(Ad, i-1), u);
481 15507216 : GEN ni = mulii(n, diviiexact(gel(A,i), di));
482 15507199 : if (equalii(d, di))
483 3655857 : gel(B, i) = ni;
484 : else
485 11851368 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
486 : }
487 : }
488 4872173 : return gerepilecopy(av, B);
489 : }
490 :
491 : /* T != 0 in Z[x], returns a monic polynomial U in Z[x] generating the
492 : * same field: there exist C in Q, L in Z such that U(x) = C T(x/L).
493 : * Set *L = NULL if L = 1, and to L otherwise. No garbage collecting. */
494 : GEN
495 0 : ZX_to_monic(GEN T, GEN *L)
496 : {
497 0 : GEN lc = leading_coeff(T);
498 0 : if (is_pm1(lc)) { *L = gen_1; return signe(lc) > 0? T: ZX_neg(T); }
499 0 : return ZX_primitive_to_monic(Q_primpart(T), L);
500 : }
501 :
502 : GEN
503 42 : poltomonic(GEN T, GEN *L)
504 : {
505 42 : pari_sp av = avma;
506 42 : if (typ(T) != t_POL || !RgX_is_QX(T)) pari_err_TYPE("poltomonic", T);
507 42 : if (degpol(T) < 0) pari_err_ROOTS0("poltomonic");
508 35 : T = ZX_Q_normalize(Q_primpart(T), L); return gc_all(av, L? 2: 1, &T, L);
509 : }
510 :
511 : GEN
512 8897 : ZXX_Q_mul(GEN A, GEN z)
513 : {
514 : long i, l;
515 : GEN B;
516 8897 : if (typ(z)==t_INT) return ZXX_Z_mul(A,z);
517 8365 : B = cgetg_copy(A, &l);
518 8365 : B[1] = A[1];
519 116767 : for (i=2; i<l; i++)
520 : {
521 108402 : GEN Ai = gel(A,i);
522 108402 : gel(B,i) = typ(Ai)==t_POL ? ZX_Q_mul(Ai, z): gmul(Ai, z);
523 : }
524 8365 : return B;
525 : }
526 :
527 : /* Evaluate pol in s using nfelt arithmetic and Horner rule */
528 : GEN
529 11151 : nfpoleval(GEN nf, GEN pol, GEN s)
530 : {
531 11151 : pari_sp av=avma;
532 11151 : long i=lg(pol)-1;
533 : GEN res;
534 11151 : if (i==1) return gen_0;
535 11151 : res = nf_to_scalar_or_basis(nf, gel(pol,i));
536 27944 : for (i-- ; i>=2; i--)
537 16793 : res = nfadd(nf, nfmul(nf, s, res), gel(pol,i));
538 11151 : return gerepileupto(av, res);
539 : }
540 :
541 : static GEN
542 94942 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
543 : {
544 94942 : pari_sp av = avma;
545 94942 : long i = lg(pol)-1;
546 : GEN res, den;
547 94942 : if (i==1) return gen_0;
548 94942 : pol = Q_remove_denom(pol, &den);
549 94944 : res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
550 394989 : for (i-- ; i>=2; i--)
551 300045 : res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
552 94944 : if (den) res = RgC_Rg_div(res, den);
553 94944 : return gerepileupto(av, res);
554 : }
555 :
556 : GEN
557 8824 : FpX_FpC_nfpoleval(GEN nf, GEN pol, GEN a, GEN p)
558 : {
559 8824 : pari_sp av=avma;
560 8824 : long i=lg(pol)-1, n=nf_get_degree(nf);
561 : GEN res, Ma;
562 8824 : if (i==1) return zerocol(n);
563 8824 : Ma = FpM_red(zk_multable(nf, a), p);
564 8824 : res = scalarcol(gel(pol,i),n);
565 70246 : for (i-- ; i>=2; i--)
566 : {
567 61422 : res = FpM_FpC_mul(Ma, res, p);
568 61422 : gel(res,1) = Fp_add(gel(res,1), gel(pol,i), p);
569 : }
570 8824 : return gerepileupto(av, res);
571 : }
572 :
573 : /* compute s(x), not stack clean */
574 : static GEN
575 95131 : ZC_galoisapply(GEN nf, GEN s, GEN x)
576 : {
577 95131 : x = nf_to_scalar_or_alg(nf, x);
578 95132 : if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
579 94943 : return QX_table_nfpoleval(nf, x, zk_multable(nf, s));
580 : }
581 :
582 : /* true nf; S = automorphism in basis form, return an FpC = S(z) mod p */
583 : GEN
584 6255 : zk_galoisapplymod(GEN nf, GEN z, GEN S, GEN p)
585 : {
586 : GEN den, pe, pe1, denpe, R;
587 :
588 6255 : z = nf_to_scalar_or_alg(nf, z);
589 6255 : if (typ(z) != t_POL) return z;
590 6255 : if (gequalX(z)) return FpC_red(S, p); /* common, e.g. modpr_genFq */
591 5877 : z = Q_remove_denom(z,&den);
592 5877 : denpe = pe = NULL;
593 5877 : pe1 = p;
594 5877 : if (den)
595 : {
596 5450 : ulong e = Z_pvalrem(den, p, &den);
597 5450 : if (e) { pe = powiu(p, e); pe1 = mulii(pe, p); }
598 5450 : denpe = Zp_inv(den, p, e+1);
599 : }
600 5877 : R = FpX_FpC_nfpoleval(nf, FpX_red(z, pe1), FpC_red(S, pe1), pe1);
601 5877 : if (denpe) R = FpC_Fp_mul(R, denpe, pe1);
602 5877 : if (pe) R = gdivexact(R, pe);
603 5877 : return R;
604 : }
605 :
606 : /* true nf */
607 : static GEN
608 7 : pr_make(GEN nf, GEN p, GEN u, GEN e, GEN f)
609 : {
610 7 : GEN t = FpM_deplin(zk_multable(nf, u), p);
611 7 : t = zk_scalar_or_multable(nf, t);
612 7 : return mkvec5(p, u, e, f, t);
613 : }
614 : static GEN
615 7 : pr_galoisapply(GEN nf, GEN pr, GEN aut)
616 : {
617 7 : GEN p = pr_get_p(pr), u = zk_galoisapplymod(nf, pr_get_gen(pr), aut, p);
618 7 : return pr_make(nf, p, u, gel(pr,3), gel(pr,4));
619 : }
620 : static GEN
621 0 : pr_galoismatrixapply(GEN nf, GEN pr, GEN M)
622 : {
623 0 : GEN p = pr_get_p(pr), u = FpC_red(ZM_ZC_mul(M, pr_get_gen(pr)), p);
624 0 : return pr_make(nf, p, u, gel(pr,3), gel(pr,4));
625 : }
626 :
627 : static GEN
628 7 : vecgaloisapply(GEN nf, GEN aut, GEN x)
629 21 : { pari_APPLY_same(galoisapply(nf, aut, gel(x,i))); }
630 : static GEN
631 0 : vecgaloismatrixapply(GEN nf, GEN aut, GEN x)
632 0 : { pari_APPLY_same(nfgaloismatrixapply(nf, aut, gel(x,i))); }
633 :
634 : /* x: famat or standard algebraic number, aut automorphism in ZC form
635 : * simplified from general galoisapply */
636 : static GEN
637 49 : elt_galoisapply(GEN nf, GEN aut, GEN x)
638 : {
639 49 : pari_sp av = avma;
640 49 : switch(typ(x))
641 : {
642 7 : case t_INT: return icopy(x);
643 7 : case t_FRAC: return gcopy(x);
644 7 : case t_POLMOD: x = gel(x,2); /* fall through */
645 14 : case t_POL: {
646 14 : GEN y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
647 14 : return gerepileupto(av,y);
648 : }
649 7 : case t_COL:
650 7 : return gerepileupto(av, ZC_galoisapply(nf, aut, x));
651 14 : case t_MAT:
652 14 : switch(lg(x)) {
653 7 : case 1: return cgetg(1, t_MAT);
654 7 : case 3: retmkmat2(vecgaloisapply(nf,aut,gel(x,1)), ZC_copy(gel(x,2)));
655 : }
656 : }
657 0 : pari_err_TYPE("galoisapply",x);
658 : return NULL; /* LCOV_EXCL_LINE */
659 : }
660 : /* M automorphism in matrix form */
661 : static GEN
662 0 : elt_galoismatrixapply(GEN nf, GEN M, GEN x)
663 : {
664 0 : if (typ(x) == t_MAT)
665 0 : switch(lg(x)) {
666 0 : case 1: return cgetg(1, t_MAT);
667 0 : case 3: retmkmat2(vecgaloismatrixapply(nf,M,gel(x,1)), ZC_copy(gel(x,2)));
668 : }
669 0 : return nfgaloismatrixapply(nf, M, x);
670 : }
671 :
672 : GEN
673 126372 : galoisapply(GEN nf, GEN aut, GEN x)
674 : {
675 126372 : pari_sp av = avma;
676 : long lx;
677 : GEN y;
678 :
679 126372 : nf = checknf(nf);
680 126372 : switch(typ(x))
681 : {
682 378 : case t_INT: return icopy(x);
683 7 : case t_FRAC: return gcopy(x);
684 :
685 35 : case t_POLMOD: x = gel(x,2); /* fall through */
686 1295 : case t_POL:
687 1295 : aut = algtobasis(nf, aut);
688 1295 : y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
689 1295 : return gerepileupto(av,y);
690 :
691 56 : case t_VEC:
692 56 : aut = algtobasis(nf, aut);
693 56 : switch(lg(x))
694 : {
695 7 : case 6:
696 7 : if (pr_is_inert(x)) { set_avma(av); return gcopy(x); }
697 7 : return gerepilecopy(av, pr_galoisapply(nf, x, aut));
698 49 : case 3: y = cgetg(3,t_VEC);
699 49 : gel(y,1) = galoisapply(nf, aut, gel(x,1));
700 49 : gel(y,2) = elt_galoisapply(nf, aut, gel(x,2));
701 49 : return gerepileupto(av, y);
702 : }
703 0 : break;
704 :
705 88558 : case t_COL:
706 88558 : aut = algtobasis(nf, aut);
707 88559 : return gerepileupto(av, ZC_galoisapply(nf, aut, x));
708 :
709 36078 : case t_MAT: /* ideal */
710 36078 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
711 36078 : if (nbrows(x) != nf_get_degree(nf)) break;
712 36078 : y = RgM_mul(nfgaloismatrix(nf,aut), x);
713 36078 : return gerepileupto(av, idealhnf_shallow(nf,y));
714 : }
715 0 : pari_err_TYPE("galoisapply",x);
716 : return NULL; /* LCOV_EXCL_LINE */
717 : }
718 :
719 : /* M automorphism in galoismatrix form */
720 : GEN
721 13909 : nfgaloismatrixapply(GEN nf, GEN M, GEN x)
722 : {
723 13909 : pari_sp av = avma;
724 : long lx;
725 : GEN y;
726 :
727 13909 : nf = checknf(nf);
728 13909 : switch(typ(x))
729 : {
730 735 : case t_INT: return icopy(x);
731 0 : case t_FRAC: return gcopy(x);
732 :
733 0 : case t_POLMOD: x = gel(x,2); /* fall through */
734 0 : case t_POL:
735 0 : x = algtobasis(nf, x);
736 0 : return gerepileupto(av, basistoalg(nf, RgM_RgC_mul(M, x)));
737 :
738 0 : case t_VEC:
739 0 : switch(lg(x))
740 : {
741 0 : case 6:
742 0 : if (pr_is_inert(x)) { set_avma(av); return gcopy(x); }
743 0 : return gerepilecopy(av, pr_galoismatrixapply(nf, x, M));
744 0 : case 3: y = cgetg(3,t_VEC);
745 0 : gel(y,1) = nfgaloismatrixapply(nf, M, gel(x,1));
746 0 : gel(y,2) = elt_galoismatrixapply(nf, M, gel(x,2));
747 0 : return gerepileupto(av, y);
748 : }
749 0 : break;
750 :
751 6027 : case t_COL: return RgM_RgC_mul(M, x);
752 :
753 7147 : case t_MAT: /* ideal */
754 7147 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
755 7147 : if (nbrows(x) != nf_get_degree(nf)) break;
756 7147 : return gerepileupto(av, idealhnf_shallow(nf,RgM_mul(M, x)));
757 : }
758 0 : pari_err_TYPE("galoisapply",x);
759 : return NULL; /* LCOV_EXCL_LINE */
760 : }
761 :
762 : /* compute action of automorphism s on nf.zk */
763 : GEN
764 94525 : nfgaloismatrix(GEN nf, GEN s)
765 : {
766 94525 : pari_sp av2, av = avma;
767 : GEN zk, D, M, H, m;
768 : long k, n;
769 :
770 94525 : nf = checknf(nf);
771 94525 : zk = nf_get_zkprimpart(nf); n = lg(zk)-1;
772 94525 : M = cgetg(n+1, t_MAT);
773 94524 : gel(M,1) = col_ei(n, 1); /* s(1) = 1 */
774 94525 : if (n == 1) return M;
775 94525 : av2 = avma;
776 94525 : if (typ(s) != t_COL) s = algtobasis(nf, s);
777 94525 : D = nf_get_zkden(nf);
778 94524 : H = RgV_to_RgM(zk, n);
779 94526 : if (n == 2)
780 : {
781 62970 : GEN t = gel(H,2); /* D * s(w_2) */
782 62970 : t = ZC_Z_add(ZC_Z_mul(s, gel(t,2)), gel(t,1));
783 62970 : gel(M,2) = gerepileupto(av2, gdiv(t, D));
784 62971 : return M;
785 : }
786 31556 : m = zk_multable(nf, s);
787 31556 : gel(M,2) = s; /* M[,k] = s(x^(k-1)) */
788 138208 : for (k = 3; k <= n; k++) gel(M,k) = ZM_ZC_mul(m, gel(M,k-1));
789 31556 : M = ZM_mul(M, H);
790 31555 : if (!equali1(D)) M = ZM_Z_divexact(M, D);
791 31554 : return gerepileupto(av, M);
792 : }
793 :
794 : static GEN
795 8162 : get_aut(GEN nf, GEN gal, GEN aut, GEN g)
796 : {
797 8162 : return aut ? gel(aut, g[1]): poltobasis(nf, galoispermtopol(gal, g));
798 : }
799 :
800 : static GEN
801 1463 : idealquasifrob(GEN nf, GEN gal, GEN grp, GEN pr, GEN subg, GEN *S, GEN aut)
802 : {
803 1463 : pari_sp av = avma;
804 1463 : long i, n = nf_get_degree(nf), f = pr_get_f(pr);
805 1463 : GEN pi = pr_get_gen(pr);
806 5624 : for (i=1; i<=n; i++)
807 : {
808 5624 : GEN g = gel(grp,i);
809 5624 : if ((!subg && perm_orderu(g) == (ulong)f)
810 5166 : || (subg && perm_relorder(g, subg)==f))
811 : {
812 1928 : *S = get_aut(nf, gal, aut, g);
813 1928 : if (ZC_prdvd(ZC_galoisapply(nf, *S, pi), pr)) return g;
814 465 : set_avma(av);
815 : }
816 : }
817 0 : pari_err_BUG("idealquasifrob [Frobenius not found]");
818 : return NULL;/*LCOV_EXCL_LINE*/
819 : }
820 :
821 : GEN
822 1379 : nfgaloispermtobasis(GEN nf, GEN gal)
823 : {
824 1379 : GEN grp = gal_get_group(gal);
825 1379 : long i, n = lg(grp)-1;
826 1379 : GEN aut = cgetg(n+1, t_VEC);
827 15323 : for(i=1; i<=n; i++)
828 : {
829 13944 : pari_sp av = avma;
830 13944 : GEN g = gel(grp, i);
831 13944 : GEN vec = poltobasis(nf, galoispermtopol(gal, g));
832 13944 : gel(aut, g[1]) = gerepileupto(av, vec);
833 : }
834 1379 : return aut;
835 : }
836 :
837 : static void
838 2485 : gal_check_pol(const char *f, GEN x, GEN y)
839 2485 : { if (!RgX_equal_var(x,y)) pari_err_MODULUS(f,x,y); }
840 :
841 : /* true nf */
842 : GEN
843 70 : idealfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN aut)
844 : {
845 70 : pari_sp av = avma;
846 70 : GEN S=NULL, g=NULL; /*-Wall*/
847 : GEN T, p, a, b, modpr;
848 : long f, n, s;
849 70 : f = pr_get_f(pr); n = nf_get_degree(nf);
850 70 : if (f==1) { set_avma(av); return identity_perm(n); }
851 70 : g = idealquasifrob(nf, gal, gal_get_group(gal), pr, NULL, &S, aut);
852 70 : if (f==2) return gerepileuptoleaf(av, g);
853 21 : modpr = zk_to_Fq_init(nf,&pr,&T,&p);
854 21 : a = pol_x(nf_get_varn(nf));
855 21 : b = nf_to_Fq(nf, zk_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
856 42 : for (s = 1; s < f-1; s++)
857 : {
858 21 : a = Fq_pow(a, p, T, p);
859 21 : if (ZX_equal(a, b)) break;
860 : }
861 21 : g = perm_powu(g, Fl_inv(s, f));
862 21 : return gerepileupto(av, g);
863 : }
864 :
865 : GEN
866 77 : idealfrobenius(GEN nf, GEN gal, GEN pr)
867 : {
868 77 : nf = checknf(nf);
869 77 : checkgal(gal);
870 77 : checkprid(pr);
871 77 : gal_check_pol("idealfrobenius",nf_get_pol(nf),gal_get_pol(gal));
872 77 : if (pr_get_e(pr)>1) pari_err_DOMAIN("idealfrobenius","pr.e", ">", gen_1,pr);
873 70 : return idealfrobenius_aut(nf, gal, pr, NULL);
874 : }
875 :
876 : /* true nf */
877 : GEN
878 616 : idealramfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN ram, GEN aut)
879 : {
880 616 : pari_sp av = avma;
881 616 : GEN S=NULL, g=NULL; /*-Wall*/
882 : GEN T, p, a, b, modpr;
883 : GEN isog, deco;
884 : long f, n, s;
885 616 : f = pr_get_f(pr); n = nf_get_degree(nf);
886 616 : if (f==1) { set_avma(av); return identity_perm(n); }
887 399 : modpr = zk_to_Fq_init(nf,&pr,&T,&p);
888 399 : deco = group_elts(gel(ram,1), nf_get_degree(nf));
889 399 : isog = group_set(gel(ram,2), nf_get_degree(nf));
890 399 : g = idealquasifrob(nf, gal, deco, pr, isog, &S, aut);
891 399 : a = pol_x(nf_get_varn(nf));
892 399 : b = nf_to_Fq(nf, zk_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
893 854 : for (s=0; !ZX_equal(a, b); s++)
894 455 : a = Fq_pow(a, p, T, p);
895 399 : g = perm_powu(g, Fl_inv(s, f));
896 399 : return gerepileupto(av, g);
897 : }
898 :
899 : GEN
900 0 : idealramfrobenius(GEN nf, GEN gal, GEN pr, GEN ram)
901 : {
902 0 : return idealramfrobenius_aut(nf, gal, pr, ram, NULL);
903 : }
904 :
905 : static GEN
906 1834 : idealinertiagroup(GEN nf, GEN gal, GEN aut, GEN pr)
907 : {
908 1834 : long i, n = nf_get_degree(nf);
909 1834 : GEN p, T, modpr = zk_to_Fq_init(nf,&pr,&T,&p);
910 1834 : GEN b = modpr_genFq(modpr);
911 1834 : long e = pr_get_e(pr), coprime = ugcd(e, pr_get_f(pr)) == 1;
912 1834 : GEN grp = gal_get_group(gal), pi = pr_get_gen(pr);
913 1834 : pari_sp ltop = avma;
914 7893 : for (i=1; i<=n; i++)
915 : {
916 7893 : GEN iso = gel(grp,i);
917 7893 : if (perm_orderu(iso) == (ulong)e)
918 : {
919 3329 : GEN S = get_aut(nf, gal, aut, iso);
920 3329 : if (ZC_prdvd(ZC_galoisapply(nf, S, pi), pr)
921 2352 : && (coprime || gequalX(nf_to_Fq(nf, galoisapply(nf,S,b), modpr))))
922 1834 : return iso;
923 1495 : set_avma(ltop);
924 : }
925 : }
926 0 : pari_err_BUG("idealinertiagroup [no isotropic element]");
927 : return NULL;/*LCOV_EXCL_LINE*/
928 : }
929 :
930 : static GEN
931 1911 : idealramgroupstame(GEN nf, GEN gal, GEN aut, GEN pr)
932 : {
933 1911 : pari_sp av = avma;
934 : GEN iso, frob, giso, isog, S, res;
935 1911 : long e = pr_get_e(pr), f = pr_get_f(pr);
936 1911 : GEN grp = gal_get_group(gal);
937 1911 : if (e == 1)
938 : {
939 77 : if (f==1)
940 0 : return cgetg(1,t_VEC);
941 77 : frob = idealquasifrob(nf, gal, grp, pr, NULL, &S, aut);
942 77 : set_avma(av);
943 77 : res = cgetg(2, t_VEC);
944 77 : gel(res, 1) = cyclicgroup(frob, f);
945 77 : return res;
946 : }
947 1834 : res = cgetg(3, t_VEC);
948 1834 : av = avma;
949 1834 : iso = idealinertiagroup(nf, gal, aut, pr);
950 1834 : set_avma(av);
951 1834 : giso = cyclicgroup(iso, e);
952 1834 : gel(res, 2) = giso;
953 1834 : if (f==1)
954 : {
955 917 : gel(res, 1) = giso;
956 917 : return res;
957 : }
958 917 : av = avma;
959 917 : isog = group_set(giso, nf_get_degree(nf));
960 917 : frob = idealquasifrob(nf, gal, grp, pr, isog, &S, aut);
961 917 : set_avma(av);
962 917 : gel(res, 1) = dicyclicgroup(iso,frob,e,f);
963 917 : return res;
964 : }
965 :
966 : /* true nf, p | e */
967 : static GEN
968 497 : idealramgroupswild(GEN nf, GEN gal, GEN aut, GEN pr)
969 : {
970 497 : pari_sp av2, av = avma;
971 497 : GEN p, T, idx, g, gbas, pi, pibas, Dpi, modpr = zk_to_Fq_init(nf,&pr,&T,&p);
972 497 : long bound, i, vDpi, vDg, n = nf_get_degree(nf);
973 497 : long e = pr_get_e(pr);
974 497 : long f = pr_get_f(pr);
975 : ulong nt,rorder;
976 497 : GEN pg, ppi, grp = gal_get_group(gal);
977 :
978 : /* G_i = {s: v(s(pi) - pi) > i} trivial for i > bound;
979 : * v_pr(Diff) = sum_{i = 0}^{bound} (#G_i - 1) >= e-1 + bound*(p-1)*/
980 497 : bound = (idealval(nf, nf_get_diff(nf), pr) - (e-1)) / (itou(p)-1);
981 497 : (void) u_pvalrem(n,p,&nt);
982 497 : rorder = e*f*(n/nt);
983 497 : idx = const_vecsmall(n,-1);
984 497 : pg = NULL;
985 497 : vDg = 0;
986 497 : if (f == 1)
987 154 : g = gbas = NULL;
988 : else
989 : {
990 : GEN Dg;
991 343 : g = nf_to_scalar_or_alg(nf, modpr_genFq(modpr));
992 343 : if (!gequalX(g)) /* p | nf.index */
993 : {
994 7 : g = Q_remove_denom(g, &Dg);
995 7 : vDg = Z_pval(Dg,p);
996 7 : pg = powiu(p, vDg + 1);
997 7 : g = FpX_red(g, pg);
998 : }
999 343 : gbas = nf_to_scalar_or_basis(nf, g);
1000 : }
1001 497 : pi = nf_to_scalar_or_alg(nf, pr_get_gen(pr));
1002 497 : pi = Q_remove_denom(pi, &Dpi);
1003 497 : vDpi = Dpi ? Z_pval(Dpi, p): 0;
1004 497 : ppi = powiu(p, vDpi + (bound + e)/e);
1005 497 : pi = FpX_red(pi, ppi);
1006 497 : pibas = nf_to_scalar_or_basis(nf, pi);
1007 497 : av2 = avma;
1008 4704 : for (i = 2; i <= n; i++)
1009 : {
1010 4207 : GEN S, Spi, piso, iso = gel(grp, i);
1011 4207 : long j, o, ix = iso[1];
1012 4207 : if (idx[ix] >= 0 || rorder % (o = (long)perm_orderu(iso))) continue;
1013 :
1014 2905 : piso = iso;
1015 2905 : S = get_aut(nf, gal, aut, iso);
1016 2905 : Spi = FpX_FpC_nfpoleval(nf, pi, FpC_red(S, ppi), ppi);
1017 : /* Computation made knowing that the valuation is <= bound + 1. Correct
1018 : * to maximal value if reduction mod ppi altered this */
1019 2905 : idx[ix] = minss(bound+1, idealval(nf, gsub(Spi,pibas), pr) - e*vDpi);
1020 2905 : if (idx[ix] == 0) idx[ix] = -1;
1021 2457 : else if (g)
1022 : {
1023 1848 : GEN Sg = pg? FpX_FpC_nfpoleval(nf, g, FpC_red(S, pg), pg): S;
1024 1848 : if (vDg)
1025 42 : { if (nfval(nf, gsub(Sg, gbas), pr) - e*vDg <= 0) idx[ix] = 0; }
1026 : else /* same, more efficient */
1027 1806 : { if (!ZC_prdvd(gsub(Sg, gbas), pr)) idx[ix] = 0; }
1028 : }
1029 5488 : for (j = 2; j < o; j++)
1030 : {
1031 2583 : piso = perm_mul(piso,iso);
1032 2583 : if (ugcd(j,o)==1) idx[ piso[1] ] = idx[ix];
1033 : }
1034 2905 : set_avma(av2);
1035 : }
1036 497 : return gerepileuptoleaf(av, idx);
1037 : }
1038 :
1039 : GEN
1040 2408 : idealramgroups_aut(GEN nf, GEN gal, GEN pr, GEN aut)
1041 : {
1042 2408 : pari_sp av = avma;
1043 : GEN tbl, idx, res, set, sub;
1044 : long i, j, e, n, maxm, p;
1045 : ulong et;
1046 2408 : nf = checknf(nf);
1047 2408 : checkgal(gal);
1048 2408 : checkprid(pr);
1049 2408 : gal_check_pol("idealramgroups",nf_get_pol(nf),gal_get_pol(gal));
1050 2408 : e = pr_get_e(pr); n = nf_get_degree(nf);
1051 2408 : p = itos(pr_get_p(pr));
1052 2408 : if (e%p) return idealramgroupstame(nf, gal, aut, pr);
1053 497 : (void) u_lvalrem(e,p,&et);
1054 497 : idx = idealramgroupswild(nf, gal, aut, pr);
1055 497 : sub = group_subgroups(galois_group(gal));
1056 497 : tbl = subgroups_tableset(sub, n);
1057 497 : maxm = vecsmall_max(idx)+1;
1058 497 : res = cgetg(maxm+1,t_VEC);
1059 497 : set = zero_F2v(n); F2v_set(set,1);
1060 2499 : for(i=maxm; i>0; i--)
1061 : {
1062 : long ix;
1063 20468 : for(j=1;j<=n;j++)
1064 18466 : if (idx[j]==i-1)
1065 3521 : F2v_set(set,j);
1066 2002 : ix = tableset_find_index(tbl, set);
1067 2002 : if (ix==0) pari_err_BUG("idealramgroups");
1068 2002 : gel(res,i) = gel(sub, ix);
1069 : }
1070 497 : return gerepilecopy(av, res);
1071 : }
1072 :
1073 : GEN
1074 126 : idealramgroups(GEN nf, GEN gal, GEN pr)
1075 : {
1076 126 : return idealramgroups_aut(nf, gal, pr, NULL);
1077 : }
1078 :
1079 : /* x = relative polynomial nf = absolute nf, bnf = absolute bnf */
1080 : GEN
1081 112 : get_bnfpol(GEN x, GEN *bnf, GEN *nf)
1082 : {
1083 112 : *bnf = checkbnf_i(x);
1084 112 : *nf = checknf_i(x);
1085 112 : if (*nf) x = nf_get_pol(*nf);
1086 112 : if (typ(x) != t_POL) pari_err_TYPE("get_bnfpol",x);
1087 112 : return x;
1088 : }
1089 :
1090 : GEN
1091 667810 : get_nfpol(GEN x, GEN *nf)
1092 : {
1093 667810 : if (typ(x) == t_POL) { *nf = NULL; return x; }
1094 381091 : *nf = checknf(x); return nf_get_pol(*nf);
1095 : }
1096 :
1097 : static GEN
1098 539 : incl_disc(GEN nfa, GEN a, int nolocal)
1099 : {
1100 : GEN d;
1101 539 : if (nfa) return nf_get_disc(nfa);
1102 469 : if (nolocal) return NULL;
1103 462 : d = ZX_disc(a);
1104 462 : if (!signe(d)) pari_err_IRREDPOL("nfisincl",a);
1105 455 : return d;
1106 : }
1107 :
1108 : static int
1109 49 : badp(GEN fa, GEN db, long q)
1110 : {
1111 49 : GEN P = gel(fa,1), E = gel(fa,2);
1112 49 : long i, l = lg(P);
1113 105 : for (i = 1; i < l; i++)
1114 56 : if (mod2(gel(E,i)) && !dvdii(db, powiu(gel(P,i),q))) return 1;
1115 49 : return 0;
1116 : }
1117 :
1118 : /* is isomorphism / inclusion (a \subset b) compatible with what we know about
1119 : * basic invariants ? (degree, signature, discriminant); test for isomorphism
1120 : * if fliso is set and for inclusion otherwise */
1121 : static int
1122 287 : tests_OK(GEN a, GEN nfa, GEN b, GEN nfb, long fliso)
1123 : {
1124 : GEN da2, da, db, fa, P, U;
1125 287 : long i, l, q, m = degpol(a), n = degpol(b);
1126 :
1127 287 : if (m <= 0) pari_err_IRREDPOL("nfisincl",a);
1128 287 : if (n <= 0) pari_err_IRREDPOL("nfisincl",b);
1129 280 : q = n / m; /* relative degree */
1130 280 : if (fliso) { if (n != m) return 0; } else { if (n % m) return 0; }
1131 280 : if (m == 1) return 1;
1132 :
1133 : /*local test expensive if n^2 >> m^4 <=> q = n/m >> m */
1134 273 : db = incl_disc(nfb, b, q > m);
1135 273 : da = db? incl_disc(nfa, a, 0): NULL;
1136 266 : if (nfa && nfb) /* both nf structures available */
1137 : {
1138 7 : long r1a = nf_get_r1(nfa), r1b = nf_get_r1(nfb);
1139 0 : return fliso ? (r1a == r1b && equalii(da, db))
1140 7 : : (r1b <= r1a * q && dvdii(db, powiu(da, q)));
1141 : }
1142 259 : if (!db) return 1;
1143 252 : if (fliso) return issquare(gdiv(da,db));
1144 :
1145 203 : if (nfa)
1146 : {
1147 7 : P = nf_get_ramified_primes(nfa); l = lg(P);
1148 14 : for (i = 1; i < l; i++)
1149 7 : if (Z_pval(db, gel(P,i)) < q * Z_pval(da, gel(P,i))) return 0;
1150 7 : return 1;
1151 : }
1152 196 : else if (nfb)
1153 : {
1154 28 : P = nf_get_ramified_primes(nfb); l = lg(P);
1155 56 : for (i = 1; i < l; i++)
1156 : {
1157 28 : GEN p = gel(P,i);
1158 28 : long va = Z_pval(nfdisc(mkvec2(a, mkvec(p))), p);
1159 28 : if (va && Z_pval(db, gel(P,i)) < va * q) return 0;
1160 : }
1161 28 : return 1;
1162 : }
1163 : /* da = dK A^2, db = dL B^2, dL = dK^q * N(D)
1164 : * da = da1 * da2, da2 maximal s.t. (da2, db) = 1: let p a prime divisor of
1165 : * da2 then p \nmid da1 * dK and p | A => v_p(da) = v_p(da2) is even */
1166 168 : da2 = Z_ppo(da, db);
1167 168 : if (!is_pm1(da2))
1168 : { /* replace da by da1 all of whose prime divisors divide db */
1169 126 : da2 = absi_shallow(da2);
1170 126 : if (!Z_issquare(da2)) return 0;
1171 14 : da = diviiexact(da, da2);
1172 : }
1173 56 : if (is_pm1(da)) return 1;
1174 49 : fa = absZ_factor_limit_strict(da, 0, &U);
1175 49 : if (badp(fa, db, q)) return 0;
1176 49 : if (U && mod2(gel(U,2)) && expi(gel(U,1)) < 150)
1177 : { /* cofactor is small, finish */
1178 0 : fa = absZ_factor(gel(U,1));
1179 0 : if (badp(fa, db, q)) return 0;
1180 : }
1181 49 : return 1;
1182 : }
1183 :
1184 : GEN
1185 77 : nfisisom(GEN a, GEN b)
1186 : {
1187 77 : pari_sp av = avma;
1188 : long i, va, vb, lx;
1189 : GEN nfa, nfb, y, la, lb;
1190 77 : int newvar, sw = 0;
1191 :
1192 77 : a = get_nfpol(a, &nfa);
1193 77 : b = get_nfpol(b, &nfb);
1194 77 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nfisisom"); }
1195 77 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nfisisom"); }
1196 77 : if (ZX_equal(a, b))
1197 : {
1198 21 : y = galoisconj(nfb? nfb: b, NULL); settyp(y, t_VEC);
1199 21 : return gerepilecopy(av,y);
1200 : }
1201 56 : if (nfa && !nfb) { swap(a,b); nfb = nfa; nfa = NULL; sw = 1; }
1202 56 : if (!tests_OK(a, nfa, b, nfb, 1)) { set_avma(av); return gen_0; }
1203 :
1204 49 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1205 49 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1206 49 : va = varn(a); vb = varn(b); newvar = (varncmp(vb,va) <= 0);
1207 49 : if (newvar) { a = leafcopy(a); setvarn(a, fetch_var_higher()); }
1208 49 : y = lift_shallow(nfroots(nfb,a));
1209 49 : if (newvar) (void)delete_var();
1210 49 : lx = lg(y); if (lx==1) { set_avma(av); return gen_0; }
1211 49 : if (sw) { vb = va; b = leafcopy(b); setvarn(b, vb); }
1212 308 : for (i=1; i<lx; i++)
1213 : {
1214 259 : GEN t = gel(y,i);
1215 259 : if (typ(t) == t_POL) setvarn(t, vb); else t = scalarpol(t, vb);
1216 259 : if (lb != gen_1) t = RgX_unscale(t, lb);
1217 259 : if (la != gen_1) t = RgX_Rg_div(t, la);
1218 259 : gel(y,i) = sw? RgXQ_reverse(t, b): t;
1219 : }
1220 49 : settyp(y, t_VEC); return gerepilecopy(av,y);
1221 : }
1222 :
1223 : static GEN
1224 7315 : partmap_reverse(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1225 : {
1226 7315 : pari_sp av = avma;
1227 7315 : GEN rnf = rnfequation2(a, t), z;
1228 7315 : if (!RgX_equal(b, gel(rnf,1)))
1229 7 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1230 7307 : z = liftpol_shallow(gel(rnf, 2));
1231 7307 : setvarn(z, v);
1232 7307 : if (!isint1(lb)) z = RgX_unscale(z, lb);
1233 7307 : if (!isint1(la)) z = RgX_Rg_div(z, la);
1234 7307 : return gerepilecopy(av, z);
1235 : }
1236 :
1237 : static GEN
1238 8085 : partmap_reverse_frac(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1239 : {
1240 8085 : pari_sp av = avma;
1241 8085 : long k = 0;
1242 : GEN N, D, G, L, de;
1243 8085 : GEN C = ZX_ZXY_resultant_all(a, Q_remove_denom(t,&de), &k, &L);
1244 8085 : if (k || degpol(b) != degpol(C))
1245 0 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1246 8085 : L = Q_primpart(L);
1247 8085 : N = gel(L,1); if (!signe(N)) { set_avma(av); return pol_0(v); }
1248 8078 : D = gel(L,2);
1249 8078 : N = RgX_neg(N); setvarn(N, v); setvarn(D, v);
1250 8078 : G = QX_gcd(N,D);
1251 8077 : if (degpol(G)) { N = RgX_div(N,G); D = RgX_div(D,G); }
1252 8077 : if (!isint1(lb)) { N = RgX_unscale(N, lb); D = RgX_unscale(D, lb); }
1253 8077 : if (!isint1(la)) D = RgX_Rg_mul(D, la);
1254 8077 : return gerepilecopy(av, mkrfrac(N,D));
1255 : }
1256 :
1257 : GEN
1258 8085 : partmap_reverse_frac_worker(GEN t, GEN a, GEN b, GEN la, GEN lb, long v)
1259 8085 : { return partmap_reverse_frac(a, b, t, la, lb, v); }
1260 :
1261 : static GEN
1262 7196 : nfisincl_from_fact(GEN a, long da, GEN b, GEN la, GEN lb, long vb, GEN y,
1263 : long flag)
1264 : {
1265 7196 : long i, k, l = lg(y), db = degpol(b), d = db / da;
1266 7196 : GEN x = cgetg(l, t_VEC);
1267 7595 : for (i= k = 1; i < l; i++)
1268 : {
1269 7413 : GEN t = gel(y,i);
1270 7413 : if (degpol(t) != d) continue;
1271 7315 : gel(x, k++) = partmap_reverse(a, b, t, la, lb, vb);
1272 7308 : if (flag) return gel(x,1);
1273 : }
1274 182 : if (k==1) return gen_0;
1275 91 : setlg(x, k);
1276 91 : gen_sort_inplace(x, (void*)&cmp_RgX, &cmp_nodata, NULL);
1277 91 : return x;
1278 : }
1279 :
1280 : static GEN
1281 1232 : nfisincl_from_fact_frac(GEN a, GEN b, GEN la, GEN lb, long vb, GEN y)
1282 : {
1283 1232 : long i, k, l = lg(y), d = degpol(b) / degpol(a);
1284 1232 : GEN worker, x = cgetg(l, t_VEC);
1285 9317 : for (i = k = 1; i < l; i++)
1286 : {
1287 8085 : GEN t = gel(y,i);
1288 8085 : if (degpol(t) != d) continue;
1289 8085 : gel(x, k++) = t;
1290 : }
1291 1232 : if (k==1) return gen_0;
1292 1232 : worker = snm_closure(is_entry("_partmap_reverse_frac_worker"),
1293 : mkvec5(a,b,la,lb,stoi(vb)));
1294 1232 : setlg(x, k); return gen_parapply(worker, x);
1295 : }
1296 :
1297 : GEN
1298 7308 : nfisincl0(GEN fa, GEN fb, long flag)
1299 : {
1300 7308 : pari_sp av = avma;
1301 : long vb;
1302 : GEN a, b, nfa, nfb, x, y, la, lb;
1303 : int newvar;
1304 7308 : if (flag < 0 || flag > 3) pari_err_FLAG("nfisincl");
1305 :
1306 7308 : a = get_nfpol(fa, &nfa);
1307 7308 : b = get_nfpol(fb, &nfb);
1308 7308 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nsisincl"); }
1309 7308 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nsisincl"); }
1310 7308 : if (ZX_equal(a, b) && flag<=1)
1311 : {
1312 21 : if (flag==1)
1313 : {
1314 7 : x = pol_x(varn(b));
1315 7 : return degpol(b) > 1 ? x: RgX_rem(x,b);
1316 : }
1317 14 : x = galoisconj(fb, NULL); settyp(x, t_VEC);
1318 14 : return gerepilecopy(av,x);
1319 : }
1320 7287 : if (flag==0 && !tests_OK(a, nfa, b, nfb, 0)) { set_avma(av); return gen_0; }
1321 :
1322 7168 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1323 7168 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1324 7168 : vb = varn(b); newvar = (varncmp(varn(a),vb) <= 0);
1325 7168 : if (newvar) { b = leafcopy(b); setvarn(b, fetch_var_higher()); }
1326 7168 : y = lift_shallow(gel(nffactor(nfa,b),1));
1327 7168 : if (flag==2)
1328 0 : x = nfisincl_from_fact_frac(a, b, la, lb, vb, y);
1329 : else
1330 7168 : x = nfisincl_from_fact(nfa, degpol(a), b, la, lb, vb, y, flag);
1331 7160 : if (newvar) (void)delete_var();
1332 7160 : return gerepilecopy(av,x);
1333 : }
1334 :
1335 : GEN
1336 112 : nfisincl(GEN fa, GEN fb) { return nfisincl0(fa, fb, 0); }
1337 :
1338 : static GEN
1339 8092 : RgF_to_Flxq(GEN F, GEN T, ulong p)
1340 : {
1341 : GEN N, D, iD;
1342 8092 : if (typ(F)==t_POL) return RgX_to_Flx(F, p);
1343 8085 : N = RgX_to_Flx(gel(F,1), p); D = RgX_to_Flx(gel(F,2), p);
1344 8085 : iD = Flxq_invsafe(D, T, p);
1345 8085 : if (!iD) return NULL;
1346 8078 : return Flxq_mul(N, iD, T, p);
1347 : }
1348 :
1349 : #define pari_APPLY_abort(EXPR)\
1350 : { \
1351 : long i, _l; \
1352 : GEN _y = cgetg_copy(x, &_l);\
1353 : for (i=1; i<_l; i++) \
1354 : { GEN _z = EXPR;\
1355 : if (!_z) return _z;\
1356 : gel(_y,i) = _z;\
1357 : } return _y;\
1358 : }
1359 :
1360 : static GEN
1361 1239 : RgFV_to_FlxqV(GEN x, GEN T, ulong p)
1362 9324 : { pari_APPLY_abort(RgF_to_Flxq(gel(x,i), T, p)) }
1363 :
1364 : static GEN
1365 1232 : nfsplitting_auto(GEN g, GEN R)
1366 : {
1367 : pari_sp av;
1368 : forprime_t T;
1369 1232 : long i, d = degpol(g);
1370 : ulong p;
1371 1232 : GEN P, K, N, G, q, den = Q_denom(R), Rp, Gp;
1372 1232 : u_forprime_init(&T, d*maxss(expu(d)-3, 2), ULONG_MAX);
1373 1232 : av = avma;
1374 79408 : for(;; set_avma(av))
1375 : {
1376 80640 : p = u_forprime_next(&T);
1377 80640 : if (dvdiu(den,p)) continue;
1378 80640 : Gp = ZX_to_Flx(g, p);
1379 80639 : if (!Flx_is_totally_split(Gp, p)) continue;
1380 1239 : P = Flx_roots(Gp, p);
1381 1239 : Rp = RgFV_to_FlxqV(R, Gp, p);
1382 1239 : if (Rp) break;
1383 7 : if (DEBUGLEVEL) err_printf("nfsplitting_auto: bad p : %lu\n",p);
1384 : }
1385 1232 : if (d == 1) return mkvec3(g, mkcol(gel(Rp,1)), utoi(p));
1386 1225 : K = Flm_Flc_invimage(FlxV_to_Flm(Rp, d), vecsmall_ei(d, 2), p);
1387 1225 : if (!K) pari_err_BUG("nfsplitting_auto");
1388 1225 : N = Flm_transpose(FlxV_Flv_multieval(Rp, P, p));
1389 1225 : q = perm_inv(vecsmall_indexsort(gel(N,1)));
1390 1225 : G = cgetg(d+1, t_COL);
1391 35707 : for (i=1; i<=d; i++)
1392 : {
1393 34482 : GEN r = perm_mul(vecsmall_indexsort(gel(N,i)), q);
1394 34482 : gel(G,i) = FlxV_Flc_mul(Rp, vecpermute(K, r), p);
1395 : }
1396 1225 : return mkvec3(g, G, utoi(p));
1397 : }
1398 :
1399 : static GEN
1400 1421 : nfsplitting_composite(GEN P)
1401 : {
1402 1421 : GEN F = gel(ZX_factor(P), 1), Q = NULL;
1403 1421 : long i, n = lg(F)-1;
1404 2849 : for (i = 1; i <= n; i++)
1405 : {
1406 1428 : GEN Fi = gel(F, i);
1407 1428 : if (degpol(Fi) == 1) continue;
1408 1400 : Q = Q ? veclast(compositum(Q, Fi)): Fi;
1409 : }
1410 1421 : return Q ? Q: pol_x(varn(P));
1411 : }
1412 : GEN
1413 1435 : nfsplitting0(GEN T0, GEN D, long flag)
1414 : {
1415 1435 : pari_sp av = avma;
1416 : long d, Ds, v;
1417 1435 : GEN T, F, K, N = NULL, lT = NULL;
1418 1435 : if (flag < 0 || flag > 3) pari_err_FLAG("nfsplitting");
1419 1435 : T = T0 = get_nfpol(T0, &K);
1420 1428 : if (!K)
1421 : {
1422 : GEN c;
1423 1407 : if (typ(T) != t_POL) pari_err_TYPE("nfsplitting",T);
1424 1407 : T = Q_primitive_part(T, &c);
1425 1407 : lT = leading_coeff(T); if (isint1(lT)) lT = NULL;
1426 1407 : if (flag && (c || lT)) pari_err_TYPE("nfsplitting", T0);
1427 1400 : RgX_check_ZX(T,"nfsplitting");
1428 : }
1429 1421 : T = nfsplitting_composite(T);
1430 1421 : if (flag && !ZX_equal(T, T0)) pari_err_IRREDPOL("nfsplitting", T0);
1431 1407 : d = degpol(T); v = varn(T);
1432 1407 : if (d <= 1 && !flag) { set_avma(av); return pol_x(v); }
1433 1379 : if (!K) {
1434 1358 : if (lT) T = polredbest(T,0);
1435 1358 : K = T;
1436 : }
1437 1379 : if (D)
1438 1239 : { if (typ(D) != t_INT || signe(D) < 1) pari_err_TYPE("nfsplitting",D); }
1439 140 : else if (d <= 7 ||
1440 35 : (d <= 11 && pari_is_dir(stack_strcat(pari_datadir, "/galdata"))))
1441 126 : D = gel(polgalois(T,DEFAULTPREC), 1);
1442 : else
1443 14 : D = mpfact(d);
1444 1379 : Ds = itos_or_0(D);
1445 1379 : T = leafcopy(T); setvarn(T, fetch_var_higher());
1446 1379 : for(F = T;;)
1447 1659 : {
1448 3038 : GEN P = gel(nffactor(K, F), 1), Q = veclast(P);
1449 3038 : if (degpol(gel(P,1)) == degpol(Q))
1450 : {
1451 1288 : if (!flag) break;
1452 1260 : P = liftall_shallow(P);
1453 1260 : if (flag==1)
1454 28 : N = nfisincl_from_fact(K, d, F, gen_1, gen_1, v, P, flag);
1455 : else
1456 1232 : N = nfisincl_from_fact_frac(T0, F, gen_1, gen_1, v, P);
1457 1260 : break;
1458 : }
1459 1750 : F = rnfequation(K,Q);
1460 1750 : if (degpol(F) == Ds && !flag) break;
1461 : }
1462 1379 : if (umodiu(D,degpol(F)))
1463 : {
1464 14 : char *sD = itostr(D);
1465 14 : pari_warn(warner,stack_strcat("ignoring incorrect degree bound ",sD));
1466 : }
1467 1379 : setvarn(F, v); (void)delete_var();
1468 1379 : if (flag) F = flag == 3? nfsplitting_auto(F, N): mkvec2(F, N);
1469 1379 : return gerepilecopy(av, F);
1470 : }
1471 :
1472 : GEN
1473 0 : nfsplitting(GEN T, GEN D) { return nfsplitting0(T, D, 0); }
1474 :
1475 : /*************************************************************************/
1476 : /** **/
1477 : /** INITALG **/
1478 : /** **/
1479 : /*************************************************************************/
1480 : typedef struct {
1481 : GEN T;
1482 : GEN ro; /* roots of T */
1483 : long r1;
1484 : GEN basden;
1485 : long prec;
1486 : long extraprec; /* possibly -1 = irrelevant or not computed */
1487 : GEN M, G; /* possibly NULL = irrelevant or not computed */
1488 : } nffp_t;
1489 :
1490 : static GEN
1491 206147 : get_roots(GEN x, long r1, long prec)
1492 : {
1493 : long i, ru;
1494 : GEN z;
1495 206147 : if (typ(x) != t_POL)
1496 : {
1497 0 : z = leafcopy(x);
1498 0 : ru = (lg(z)-1 + r1) >> 1;
1499 : }
1500 : else
1501 : {
1502 206147 : long n = degpol(x);
1503 206146 : z = (r1 == n)? ZX_realroots_irred(x, prec): QX_complex_roots(x,prec);
1504 206151 : ru = (n+r1)>>1;
1505 : }
1506 449641 : for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
1507 206151 : z[0]=evaltyp(t_VEC)|evallg(ru+1); return z;
1508 : }
1509 :
1510 : GEN
1511 0 : nf_get_allroots(GEN nf)
1512 : {
1513 0 : return embed_roots(nf_get_roots(nf), nf_get_r1(nf));
1514 : }
1515 :
1516 : /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
1517 : static GEN
1518 381515 : quicktrace(GEN x, GEN sym)
1519 : {
1520 381515 : GEN p1 = gen_0;
1521 : long i;
1522 :
1523 381515 : if (typ(x) != t_POL) return gmul(x, gel(sym,1));
1524 381515 : if (signe(x))
1525 : {
1526 381516 : sym--;
1527 2798983 : for (i=lg(x)-1; i>1; i--)
1528 2417524 : p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
1529 : }
1530 381458 : return p1;
1531 : }
1532 :
1533 : static GEN
1534 83502 : get_Tr(GEN mul, GEN x, GEN basden)
1535 : {
1536 83502 : GEN t, bas = gel(basden,1), den = gel(basden,2);
1537 83502 : long i, j, n = lg(bas)-1;
1538 83502 : GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
1539 :
1540 83501 : gel(TW,1) = utoipos(n);
1541 262661 : for (i=2; i<=n; i++)
1542 : {
1543 179161 : t = quicktrace(gel(bas,i), sym);
1544 179154 : if (den && gel(den,i)) t = diviiexact(t,gel(den,i));
1545 179152 : gel(TW,i) = t; /* tr(w[i]) */
1546 : }
1547 83500 : gel(T,1) = TW;
1548 262652 : for (i=2; i<=n; i++)
1549 : {
1550 179155 : gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
1551 707405 : for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
1552 528253 : gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
1553 : }
1554 83497 : return T;
1555 : }
1556 :
1557 : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
1558 : static GEN
1559 202350 : get_bas_den(GEN bas)
1560 : {
1561 202350 : GEN b,d,den, dbas = leafcopy(bas);
1562 202351 : long i, l = lg(bas);
1563 202351 : int power = 1;
1564 202351 : den = cgetg(l,t_VEC);
1565 937809 : for (i=1; i<l; i++)
1566 : {
1567 735459 : b = Q_remove_denom(gel(bas,i), &d);
1568 735446 : gel(dbas,i) = b;
1569 735446 : gel(den,i) = d; if (d) power = 0;
1570 : }
1571 202350 : if (power) den = NULL; /* power basis */
1572 202350 : return mkvec2(dbas, den);
1573 : }
1574 :
1575 : /* return multiplication table for S->basis */
1576 : static GEN
1577 83499 : nf_multable(nfmaxord_t *S, GEN invbas)
1578 : {
1579 83499 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1580 83499 : long i,j, n = degpol(T);
1581 83499 : GEN mul = cgetg(n*n+1,t_MAT);
1582 :
1583 : /* i = 1 split for efficiency, assume w[1] = 1 */
1584 346166 : for (j=1; j<=n; j++)
1585 262661 : gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
1586 262667 : for (i=2; i<=n; i++)
1587 707419 : for (j=i; j<=n; j++)
1588 : {
1589 528257 : pari_sp av = avma;
1590 528257 : GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1591 528256 : z = ZM_ZX_mul(invbas, z); /* integral column */
1592 528228 : if (den)
1593 : {
1594 345679 : GEN d = mul_denom(gel(den,i), gel(den,j));
1595 345674 : if (d) z = ZC_Z_divexact(z, d);
1596 : }
1597 528229 : gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
1598 : }
1599 83505 : return mul;
1600 : }
1601 :
1602 : /* as get_Tr, mul_table not precomputed */
1603 : static GEN
1604 27945 : make_Tr(nfmaxord_t *S)
1605 : {
1606 27945 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1607 27945 : long i,j, n = degpol(T);
1608 27945 : GEN c, t, d, M = cgetg(n+1,t_MAT), sym = polsym(T, n-1);
1609 :
1610 : /* W[i] = w[i]/den[i]; assume W[1] = 1, case i = 1 split for efficiency */
1611 27946 : c = cgetg(n+1,t_COL); gel(M,1) = c;
1612 27946 : gel(c, 1) = utoipos(n);
1613 83376 : for (j=2; j<=n; j++)
1614 : {
1615 55430 : pari_sp av = avma;
1616 55430 : t = quicktrace(gel(w,j), sym);
1617 55428 : if (den)
1618 : {
1619 35241 : d = gel(den,j);
1620 35241 : if (d) t = diviiexact(t, d);
1621 : }
1622 55424 : gel(c,j) = gerepileuptoint(av, t);
1623 : }
1624 83370 : for (i=2; i<=n; i++)
1625 : {
1626 55430 : c = cgetg(n+1,t_COL); gel(M,i) = c;
1627 202365 : for (j=1; j<i ; j++) gel(c,j) = gcoeff(M,i,j);
1628 202357 : for ( ; j<=n; j++)
1629 : {
1630 146933 : pari_sp av = avma;
1631 146933 : t = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1632 146936 : t = quicktrace(t, sym);
1633 146917 : if (den)
1634 : {
1635 114648 : d = mul_denom(gel(den,i),gel(den,j));
1636 114645 : if (d) t = diviiexact(t, d);
1637 : }
1638 146911 : gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
1639 : }
1640 : }
1641 27940 : return M;
1642 : }
1643 :
1644 : /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
1645 : static void
1646 237320 : make_M(nffp_t *F, int trunc)
1647 : {
1648 237320 : GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
1649 : GEN m, d, M;
1650 237320 : long i, j, l = lg(ro), n = lg(bas);
1651 237320 : M = cgetg(n,t_MAT);
1652 237323 : gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
1653 839879 : for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
1654 773226 : for (i=1; i<l; i++)
1655 : {
1656 535919 : GEN r = gel(ro,i), ri;
1657 535919 : ri = (gexpo(r) > 1)? ginv(r): NULL;
1658 2669848 : for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
1659 : }
1660 237307 : if (den)
1661 524759 : for (j=2; j<n; j++)
1662 : {
1663 407223 : d = gel(den,j); if (!d) continue;
1664 343979 : m = gel(M,j);
1665 1715419 : for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
1666 : }
1667 :
1668 237291 : if (trunc && gprecision(M) > F->prec)
1669 : {
1670 18116 : M = gprec_w(M, F->prec);
1671 18116 : F->ro = gprec_w(ro,F->prec);
1672 : }
1673 237291 : F->M = M;
1674 237291 : }
1675 :
1676 : /* return G real such that G~ * G = T_2 */
1677 : static void
1678 237322 : make_G(nffp_t *F)
1679 : {
1680 237322 : GEN G, M = F->M;
1681 237322 : long i, j, k, r1 = F->r1, l = lg(M);
1682 :
1683 237322 : if (r1 == l-1) { F->G = M; return; }
1684 190096 : G = cgetg(l, t_MAT);
1685 895234 : for (j = 1; j < l; j++)
1686 : {
1687 705188 : GEN g, m = gel(M,j);
1688 705188 : gel(G,j) = g = cgetg(l, t_COL);
1689 1146723 : for (k = i = 1; i <= r1; i++) gel(g,k++) = gel(m,i);
1690 2406490 : for ( ; k < l; i++)
1691 : {
1692 1701356 : GEN r = gel(m,i);
1693 1701356 : if (typ(r) == t_COMPLEX)
1694 : {
1695 1397392 : GEN a = gel(r,1), b = gel(r,2);
1696 1397392 : gel(g,k++) = mpadd(a, b);
1697 1397351 : gel(g,k++) = mpsub(a, b);
1698 : }
1699 : else
1700 : {
1701 303964 : gel(g,k++) = r;
1702 303964 : gel(g,k++) = r;
1703 : }
1704 : }
1705 : }
1706 190046 : F->G = G;
1707 : }
1708 :
1709 : static long
1710 269483 : prec_fix(long prec)
1711 : {
1712 : #ifndef LONG_IS_64BIT
1713 : /* make sure that default accuracy is the same on 32/64bit */
1714 38558 : if (odd(prec2lg(prec))) prec+=EXTRAPRECWORD;
1715 : #endif
1716 269483 : return prec;
1717 : }
1718 : static void
1719 237320 : make_M_G(nffp_t *F, int trunc)
1720 : {
1721 : long n, eBD, prec;
1722 237320 : if (F->extraprec < 0)
1723 : { /* not initialized yet; compute roots so that absolute accuracy
1724 : * of M & G >= prec */
1725 : double er;
1726 237299 : n = degpol(F->T);
1727 237300 : eBD = 1 + gexpo(gel(F->basden,1));
1728 237302 : er = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->T);
1729 237296 : if (er < 0) er = 0;
1730 237296 : F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
1731 : }
1732 237318 : prec = prec_fix(F->prec + F->extraprec);
1733 237317 : if (!F->ro || gprecision(gel(F->ro,1)) < prec)
1734 206145 : F->ro = get_roots(F->T, F->r1, prec);
1735 :
1736 237320 : make_M(F, trunc);
1737 237322 : make_G(F);
1738 237314 : }
1739 :
1740 : static void
1741 174465 : nffp_init(nffp_t *F, nfmaxord_t *S, long prec)
1742 : {
1743 174465 : F->T = S->T;
1744 174465 : F->r1 = S->r1;
1745 174465 : F->basden = S->basden;
1746 174465 : F->ro = NULL;
1747 174465 : F->extraprec = -1;
1748 174465 : F->prec = prec;
1749 174465 : }
1750 :
1751 : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
1752 : * basis. Assume bas[1] = 1 and that the leading coefficient of elements
1753 : * of bas are of the form 1/b for a t_INT b */
1754 : static GEN
1755 1666 : get_nfindex(GEN bas)
1756 : {
1757 1666 : pari_sp av = avma;
1758 1666 : long n = lg(bas)-1, i;
1759 : GEN D, d, mat;
1760 :
1761 : /* assume bas[1] = 1 */
1762 1666 : D = gel(bas,1);
1763 1666 : if (! is_pm1(simplify_shallow(D))) pari_err_TYPE("get_nfindex", D);
1764 1666 : D = gen_1;
1765 8799 : for (i = 2; i <= n; i++)
1766 : { /* after nfbasis, basis is upper triangular! */
1767 7140 : GEN B = gel(bas,i), lc;
1768 7140 : if (degpol(B) != i-1) break;
1769 7133 : lc = gel(B, i+1);
1770 7133 : switch (typ(lc))
1771 : {
1772 2772 : case t_INT: continue;
1773 4361 : case t_FRAC: if (is_pm1(gel(lc,1)) ) {D = mulii(D, gel(lc,2)); continue;}
1774 0 : default: pari_err_TYPE("get_nfindex", B);
1775 : }
1776 : }
1777 1666 : if (i <= n)
1778 : { /* not triangular after all */
1779 7 : bas = vecslice(bas,i,n);
1780 7 : bas = Q_remove_denom(bas, &d);
1781 7 : if (!d) return D;
1782 7 : mat = RgV_to_RgM(bas, n);
1783 7 : mat = rowslice(mat, i,n);
1784 7 : D = mulii(D, diviiexact(powiu(d, n-i+1), absi_shallow(ZM_det(mat))));
1785 : }
1786 1666 : return gerepileuptoint(av, D);
1787 : }
1788 : /* make sure all components of S are initialized */
1789 : static void
1790 167008 : nfmaxord_complete(nfmaxord_t *S)
1791 : {
1792 167008 : if (!S->dT) S->dT = ZX_disc(S->T);
1793 167008 : if (!S->index)
1794 : {
1795 1666 : if (S->dK) /* fast */
1796 0 : S->index = sqrti( diviiexact(S->dT, S->dK) );
1797 : else
1798 1666 : S->index = get_nfindex(S->basis);
1799 : }
1800 167008 : if (!S->dK) S->dK = diviiexact(S->dT, sqri(S->index));
1801 167008 : if (S->r1 < 0) S->r1 = ZX_sturm_irred(S->T);
1802 167006 : if (!S->basden) S->basden = get_bas_den(S->basis);
1803 167006 : }
1804 :
1805 : GEN
1806 83501 : nfmaxord_to_nf(nfmaxord_t *S, GEN ro, long prec)
1807 : {
1808 83501 : GEN nf = cgetg(10,t_VEC);
1809 83502 : GEN T = S->T, Tr, D, w, A, dA, MDI, mat = cgetg(9,t_VEC);
1810 83506 : long n = degpol(T);
1811 : nffp_t F;
1812 83503 : nfmaxord_complete(S);
1813 83503 : nffp_init(&F,S,prec);
1814 83502 : F.ro = ro;
1815 83502 : make_M_G(&F, 0);
1816 :
1817 83506 : gel(nf,1) = S->T;
1818 83506 : gel(nf,2) = mkvec2s(S->r1, (n - S->r1)>>1);
1819 83506 : gel(nf,3) = S->dK;
1820 83506 : gel(nf,4) = S->index;
1821 83506 : gel(nf,5) = mat;
1822 83506 : if (gprecision(gel(F.ro,1)) > prec) F.ro = gprec_wtrunc(F.ro, prec);
1823 83505 : gel(nf,6) = F.ro;
1824 83505 : w = S->basis;
1825 83505 : if (!is_pm1(S->index)) w = Q_remove_denom(w, NULL);
1826 83503 : gel(nf,7) = w;
1827 83503 : gel(nf,8) = ZM_inv(RgV_to_RgM(w,n), NULL);
1828 83500 : gel(nf,9) = nf_multable(S, nf_get_invzk(nf));
1829 83503 : gel(mat,1) = F.M;
1830 83503 : gel(mat,2) = F.G;
1831 :
1832 83503 : Tr = get_Tr(gel(nf,9), T, S->basden);
1833 83499 : gel(mat,6) = A = ZM_inv(Tr, &dA); /* dA T^-1, primitive */
1834 83504 : A = ZM_hnfmodid(A, dA);
1835 : /* CAVEAT: nf is not complete yet, but the fields needed for
1836 : * idealtwoelt, zk_scalar_or_multable and idealinv are present ! */
1837 83505 : MDI = idealtwoelt(nf, A);
1838 83506 : gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
1839 83501 : gel(mat,7) = MDI;
1840 83501 : if (is_pm1(S->index))
1841 : { /* principal ideal (T'), whose norm is |dK| */
1842 48642 : D = zk_scalar_or_multable(nf, ZX_deriv(T));
1843 48642 : if (typ(D) == t_MAT) D = ZM_hnfmod(D, absi_shallow(S->dK));
1844 : }
1845 : else
1846 : {
1847 34859 : GEN c = diviiexact(dA, gcoeff(A,1,1));
1848 34859 : D = idealHNF_inv_Z(nf, A); /* (A\cap Z) / A */
1849 34860 : if (!is_pm1(c)) D = ZM_Z_mul(D, c);
1850 : }
1851 83505 : gel(mat,3) = RM_round_maxrank(F.G);
1852 83501 : gel(mat,4) = Tr;
1853 83501 : gel(mat,5) = D;
1854 83501 : w = S->dKP; if (!w) { w = gel(absZ_factor(S->dK), 1); settyp(w, t_VEC); }
1855 83501 : gel(mat,8) = w; return nf;
1856 : }
1857 :
1858 : static GEN
1859 3101 : primes_certify(GEN dK, GEN dKP)
1860 : {
1861 3101 : long i, l = lg(dKP);
1862 3101 : GEN v, w, D = dK;
1863 3101 : v = vectrunc_init(l);
1864 3101 : w = vectrunc_init(l);
1865 9695 : for (i = 1; i < l; i++)
1866 : {
1867 6594 : GEN p = gel(dKP,i);
1868 6594 : vectrunc_append(isprime(p)? w: v, p);
1869 6594 : (void)Z_pvalrem(D, p, &D);
1870 : }
1871 3101 : if (!is_pm1(D))
1872 : {
1873 0 : if (signe(D) < 0) D = negi(D);
1874 0 : vectrunc_append(isprime(D)? w: v, D);
1875 : }
1876 3101 : return mkvec2(v,w);
1877 : }
1878 : GEN
1879 3101 : nfcertify(GEN nf)
1880 : {
1881 3101 : pari_sp av = avma;
1882 : GEN vw;
1883 3101 : nf = checknf(nf);
1884 3101 : vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
1885 3101 : return gerepilecopy(av, gel(vw,1));
1886 : }
1887 :
1888 : /* set *pro to roots of S->T */
1889 : static GEN
1890 72835 : get_red_G(nfmaxord_t *S, GEN *pro)
1891 : {
1892 72835 : pari_sp av = avma;
1893 72835 : GEN G, u, u0 = NULL;
1894 72835 : long prec, n = degpol(S->T);
1895 : nffp_t F;
1896 :
1897 72834 : prec = nbits2prec(n+32);
1898 72833 : nffp_init(&F, S, prec);
1899 : for (;;)
1900 : {
1901 72833 : F.prec = prec; make_M_G(&F, 0); G = F.G;
1902 72832 : if (u0) G = RgM_mul(G, u0);
1903 72832 : if (DEBUGLEVEL)
1904 0 : err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
1905 0 : prec + F.extraprec, prec, F.extraprec);
1906 72832 : if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST)))
1907 : {
1908 72833 : if (lg(u)-1 == n) break;
1909 : /* singular ==> loss of accuracy */
1910 0 : if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
1911 0 : else u0 = gerepilecopy(av, u);
1912 : }
1913 0 : prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
1914 0 : F.ro = NULL;
1915 0 : if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
1916 : }
1917 72833 : if (u0) u = RgM_mul(u0,u);
1918 72833 : *pro = F.ro; return u;
1919 : }
1920 :
1921 : /* Compute an LLL-reduced basis for the integer basis of nf(T).
1922 : * set *pro = roots of x if computed [NULL if not computed] */
1923 : static void
1924 101360 : set_LLL_basis(nfmaxord_t *S, GEN *pro, long flag, double DELTA)
1925 : {
1926 101360 : GEN B = S->basis;
1927 101360 : long N = degpol(S->T);
1928 101360 : if (S->r1 < 0)
1929 : {
1930 17955 : S->r1 = ZX_sturm_irred(S->T);
1931 17955 : if (odd(N - S->r1)) pari_err_IRREDPOL("set_LLL_basis", S->T);
1932 : }
1933 101353 : if (!S->basden) S->basden = get_bas_den(B);
1934 101353 : *pro = NULL; if (flag & nf_NOLLL) return;
1935 100779 : if (S->r1 == N) {
1936 27945 : pari_sp av = avma;
1937 27945 : GEN u = ZM_lll(make_Tr(S), DELTA, LLL_GRAM|LLL_KEEP_FIRST|LLL_IM);
1938 27946 : B = gerepileupto(av, RgV_RgM_mul(B, u));
1939 : }
1940 : else
1941 72834 : B = RgV_RgM_mul(B, get_red_G(S, pro));
1942 100780 : S->basis = B;
1943 100780 : S->basden = get_bas_den(B);
1944 : }
1945 :
1946 : /* = 1 iff |a| > |b| or equality and a > 0 */
1947 : static int
1948 125758 : cmpii_polred(GEN a, GEN b)
1949 : {
1950 125758 : int fl = abscmpii(a, b);
1951 : long sa, sb;
1952 125758 : if (fl) return fl;
1953 104116 : sa = signe(a);
1954 104116 : sb = signe(b);
1955 104116 : if (sa == sb) return 0;
1956 746 : return sa == 1? 1: -1;
1957 : }
1958 : static int
1959 29968 : ZX_cmp(GEN x, GEN y)
1960 29968 : { return gen_cmp_RgX((void*)cmpii_polred, x, y); }
1961 : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
1962 : * (if so update *dx); both x and y are monic */
1963 : static int
1964 51612 : ZX_is_better(GEN y, GEN x, GEN *dx)
1965 : {
1966 : pari_sp av;
1967 : int cmp;
1968 : GEN d;
1969 51612 : if (!*dx) *dx = ZX_disc(x);
1970 51612 : av = avma; d = ZX_disc(y);
1971 51611 : cmp = abscmpii(d, *dx);
1972 51610 : if (cmp < 0) { *dx = d; return 1; }
1973 47116 : return gc_bool(av, cmp? 0: (ZX_cmp(y, x) < 0));
1974 : }
1975 :
1976 : static void polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
1977 : /* Seek a simpler, polynomial pol defining the same number field as
1978 : * x (assumed to be monic at this point) */
1979 : static GEN
1980 280 : nfpolred(nfmaxord_t *S, GEN *pro)
1981 : {
1982 280 : GEN x = S->T, dx, b, rev;
1983 280 : long n = degpol(x), v = varn(x);
1984 :
1985 280 : if (n == 1) {
1986 98 : S->T = pol_x(v);
1987 98 : *pro = NULL;
1988 98 : return scalarpol_shallow(negi(gel(x,2)), v);
1989 : }
1990 182 : polredbest_aux(S, pro, &x, &dx, &b);
1991 182 : if (x == S->T) return NULL; /* no improvement */
1992 147 : if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
1993 :
1994 : /* update T */
1995 147 : rev = QXQ_reverse(b, S->T);
1996 147 : S->basis = QXV_QXQ_eval(S->basis, rev, x);
1997 147 : S->index = sqrti( diviiexact(dx,S->dK) );
1998 147 : S->basden = get_bas_den(S->basis);
1999 147 : S->dT = dx;
2000 147 : S->T = x;
2001 147 : *pro = NULL; /* reset */
2002 147 : return rev;
2003 : }
2004 :
2005 : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
2006 : * basis (deprecated), or listP data (nfbasis input format) to specify
2007 : * a set of primes at with the basis order must be maximal.
2008 : * 1) nf type (or unrecognized): return t_VEC
2009 : * 2) ZX or [ZX, listP]: return t_POL
2010 : * 3) [ZX, order basis]: return 0 (deprecated)
2011 : * incorrect: return -1 */
2012 : static long
2013 85339 : nf_input_type(GEN x)
2014 : {
2015 85339 : GEN T, V, DKP = NULL;
2016 : long i, d, v;
2017 85339 : switch(typ(x))
2018 : {
2019 79704 : case t_POL: return t_POL;
2020 5635 : case t_VEC:
2021 5635 : switch(lg(x))
2022 : {
2023 1869 : case 4: DKP = gel(x,3);
2024 5607 : case 3: break;
2025 28 : default: return t_VEC; /* nf or incorrect */
2026 : }
2027 5607 : T = gel(x,1); V = gel(x,2);
2028 5607 : if (typ(T) != t_POL) return -1;
2029 5607 : switch(typ(V))
2030 : {
2031 224 : case t_INT: case t_MAT: return t_POL;
2032 5383 : case t_VEC: case t_COL:
2033 5383 : if (RgV_is_ZV(V)) return t_POL;
2034 2534 : break;
2035 0 : default: return -1;
2036 : }
2037 2534 : d = degpol(T); v = varn(T);
2038 2534 : if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
2039 16289 : for (i = 1; i <= d; i++)
2040 : { /* check integer basis */
2041 13776 : GEN c = gel(V,i);
2042 13776 : switch(typ(c))
2043 : {
2044 35 : case t_INT: break;
2045 13741 : case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
2046 : /* fall through */
2047 14 : default: return -1;
2048 : }
2049 : }
2050 2513 : if (DKP && (typ(DKP) != t_VEC || !RgV_is_ZV(DKP))) return -1;
2051 2513 : return 0;
2052 : }
2053 0 : return t_VEC; /* nf or incorrect */
2054 : }
2055 :
2056 : /* cater for obsolete nf_PARTIALFACT flag */
2057 : static void
2058 3906 : nfinit_basic_partial(nfmaxord_t *S, GEN T)
2059 : {
2060 3906 : if (typ(T) == t_POL) { nfmaxord(S, mkvec2(T,utoipos(500000)), 0); }
2061 14 : else nfinit_basic(S, T);
2062 3906 : }
2063 : static void
2064 14119 : nfinit_basic_flag(nfmaxord_t *S, GEN x, long flag)
2065 : {
2066 14119 : if (flag & nf_PARTIALFACT)
2067 35 : nfinit_basic_partial(S, x);
2068 : else
2069 14084 : nfinit_basic(S, x);
2070 14112 : }
2071 :
2072 : /* true nf */
2073 : static GEN
2074 62863 : nf_basden(GEN nf)
2075 : {
2076 62863 : GEN zkD = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
2077 62863 : D = equali1(D)? NULL: const_vec(lg(zkD)-1, D);
2078 62863 : return mkvec2(zkD, D);
2079 : }
2080 : void
2081 85341 : nfinit_basic(nfmaxord_t *S, GEN T)
2082 : {
2083 85341 : switch (nf_input_type(T))
2084 : {
2085 82777 : case t_POL: nfmaxord(S, T, 0); return;
2086 28 : case t_VEC:
2087 : { /* nf, bnf, bnr */
2088 28 : GEN nf = checknf(T);
2089 28 : S->T = S->T0 = nf_get_pol(nf);
2090 28 : S->basis = nf_get_zk(nf); /* probably useless */
2091 28 : S->basden = nf_basden(nf);
2092 28 : S->index = nf_get_index(nf);
2093 28 : S->dK = nf_get_disc(nf);
2094 28 : S->dKP = nf_get_ramified_primes(nf);
2095 28 : S->dT = mulii(S->dK, sqri(S->index));
2096 28 : S->r1 = nf_get_r1(nf); break;
2097 : }
2098 2513 : case 0: /* monic integral polynomial + integer basis (+ ramified primes)*/
2099 2513 : S->T = S->T0 = gel(T,1);
2100 2513 : S->basis = gel(T,2);
2101 2513 : S->basden = NULL;
2102 2513 : S->index = NULL;
2103 2513 : S->dK = NULL;
2104 2513 : S->dKP = NULL;
2105 2513 : if (lg(T) == 4)
2106 : {
2107 1869 : GEN v = gel(T,3); if (typ(v) == t_COL) v = shallowtrans(v);
2108 1869 : S->dKP = v;
2109 : }
2110 2513 : S->dT = NULL;
2111 2513 : S->r1 = -1; break;
2112 20 : default: /* -1 */
2113 20 : pari_err_TYPE("nfinit_basic", T);
2114 : }
2115 2541 : S->dTP = S->dTE = S->dKE = NULL;
2116 2541 : S->unscale = gen_1;
2117 : }
2118 :
2119 : GEN
2120 83503 : nfinit_complete(nfmaxord_t *S, long flag, long prec)
2121 : {
2122 83503 : GEN nf, unscale = S->unscale, rev = NULL;
2123 :
2124 83503 : if (!ZX_is_irred(S->T)) pari_err_IRREDPOL("nfinit",S->T);
2125 83505 : if (!(flag & nf_RED) && !ZX_is_monic(S->T0))
2126 : {
2127 49 : pari_warn(warner,"nonmonic polynomial. Result of the form [nf,c]");
2128 49 : flag |= nf_RED | nf_ORIG;
2129 : }
2130 83505 : if (!(flag & nf_RED) && !isint1(unscale))
2131 : { /* implies lc(x0) = 1 and L := 1/unscale is integral */
2132 4403 : long d = degpol(S->T0);
2133 4403 : GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
2134 4403 : GEN f = powiu(L, (d*(d-1)) >> 1);
2135 4402 : S->T = S->T0; /* restore original user-supplied x0, unscale data */
2136 4402 : S->unscale = gen_1;
2137 4402 : S->dT = gmul(S->dT, sqri(f));
2138 4402 : S->basis = RgXV_unscale(S->basis, unscale);
2139 4403 : S->index = gmul(S->index, f);
2140 : }
2141 83505 : nfmaxord_complete(S); /* more expensive after set_LLL_basis */
2142 83503 : if (flag & nf_RED)
2143 : {
2144 : GEN ro;
2145 : /* lie to polred: more efficient to update *after* modreverse, than to
2146 : * unscale in the polred subsystem */
2147 280 : S->unscale = gen_1;
2148 280 : rev = nfpolred(S, &ro);
2149 280 : nf = nfmaxord_to_nf(S, ro, prec);
2150 280 : S->unscale = unscale; /* restore */
2151 : }
2152 : else
2153 : {
2154 83223 : GEN ro; set_LLL_basis(S, &ro, flag, 0.99);
2155 83221 : nf = nfmaxord_to_nf(S, ro, prec);
2156 : }
2157 83501 : if (flag & nf_ORIG)
2158 : {
2159 84 : if (!rev)
2160 : { /* no improvement */
2161 28 : long v = varn(S->T);
2162 28 : rev = degpol(S->T) == 1? pol_0(v): pol_x(v);
2163 : }
2164 84 : if (!isint1(unscale)) rev = RgX_Rg_div(rev, unscale);
2165 84 : nf = mkvec2(nf, mkpolmod(rev, S->T));
2166 : }
2167 83501 : return nf;
2168 : }
2169 : /* Initialize the number field defined by the polynomial x (in variable v)
2170 : * flag & nf_RED: try a polred first.
2171 : * flag & nf_ORIG: return [nfinit(x), Mod(a,red)], where
2172 : * Mod(a,red) = Mod(v,x) (i.e return the base change). */
2173 : GEN
2174 10069 : nfinit0(GEN x, long flag,long prec)
2175 : {
2176 10069 : const pari_sp av = avma;
2177 : nfmaxord_t S;
2178 10069 : if (flag < 0 || flag > 7) pari_err_FLAG("nfinit");
2179 10069 : if (checkrnf_i(x)) return rnf_build_nfabs(x, prec);
2180 10048 : nfinit_basic(&S, x);
2181 10027 : return gerepilecopy(av, nfinit_complete(&S, flag, prec));
2182 : }
2183 : GEN
2184 182 : nfinitred(GEN x, long prec) { return nfinit0(x, nf_RED, prec); }
2185 : GEN
2186 0 : nfinitred2(GEN x, long prec) { return nfinit0(x, nf_RED|nf_ORIG, prec); }
2187 : GEN
2188 6226 : nfinit(GEN x, long prec) { return nfinit0(x, 0, prec); }
2189 :
2190 : /* assume x a bnr/bnf/nf */
2191 : long
2192 897156 : nf_get_prec(GEN x)
2193 : {
2194 897156 : GEN nf = checknf(x), ro = nf_get_roots(nf);
2195 897161 : return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
2196 : }
2197 :
2198 : /* true nf */
2199 : GEN
2200 62835 : nfnewprec_shallow(GEN nf, long prec)
2201 : {
2202 62835 : GEN m, NF = leafcopy(nf);
2203 : nffp_t F;
2204 :
2205 62835 : F.T = nf_get_pol(nf);
2206 62835 : F.ro = NULL;
2207 62835 : F.r1 = nf_get_r1(nf);
2208 62835 : F.basden = nf_basden(nf);
2209 62835 : F.extraprec = -1;
2210 62835 : F.prec = prec; make_M_G(&F, 0);
2211 62831 : gel(NF,5) = m = leafcopy(gel(NF,5));
2212 62835 : gel(m,1) = F.M;
2213 62835 : gel(m,2) = F.G;
2214 62835 : gel(NF,6) = F.ro; return NF;
2215 : }
2216 :
2217 : GEN
2218 154 : nfnewprec(GEN nf, long prec)
2219 : {
2220 : GEN z;
2221 154 : switch(nftyp(nf))
2222 : {
2223 49 : default: pari_err_TYPE("nfnewprec", nf);
2224 14 : case typ_BNF: z = bnfnewprec(nf,prec); break;
2225 7 : case typ_BNR: z = bnrnewprec(nf,prec); break;
2226 84 : case typ_NF: {
2227 84 : pari_sp av = avma;
2228 84 : z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
2229 84 : break;
2230 : }
2231 : }
2232 105 : return z;
2233 : }
2234 :
2235 : /********************************************************************/
2236 : /** **/
2237 : /** POLRED **/
2238 : /** **/
2239 : /********************************************************************/
2240 : GEN
2241 0 : embednorm_T2(GEN x, long r1)
2242 : {
2243 0 : pari_sp av = avma;
2244 0 : GEN p = RgV_sumpart(x, r1);
2245 0 : GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
2246 0 : if (q != gen_0) p = gadd(p, gmul2n(q,1));
2247 0 : return avma == av? gcopy(p): gerepileupto(av, p);
2248 : }
2249 :
2250 : /* simplified version of gnorm for scalar, noncomplex inputs, without GC */
2251 : static GEN
2252 160672 : real_norm(GEN x)
2253 : {
2254 160672 : switch(typ(x))
2255 : {
2256 0 : case t_INT: return sqri(x);
2257 160672 : case t_REAL: return sqrr(x);
2258 0 : case t_FRAC: return sqrfrac(x);
2259 : }
2260 0 : pari_err_TYPE("real_norm", x);
2261 : return NULL;/*LCOV_EXCL_LINE*/
2262 : }
2263 : /* simplified version of gnorm, without GC */
2264 : static GEN
2265 34896899 : complex_norm(GEN x)
2266 : {
2267 34896899 : return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
2268 : }
2269 : /* return T2(x), argument r1 needed in case x has components whose type
2270 : * is unexpected, e.g. all of them t_INT for embed(gen_1) */
2271 : GEN
2272 70948 : embed_T2(GEN x, long r1)
2273 : {
2274 70948 : pari_sp av = avma;
2275 70948 : long i, l = lg(x);
2276 70948 : GEN c, s = NULL, t = NULL;
2277 70948 : if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
2278 231619 : for (i = 1; i <= r1; i++)
2279 : {
2280 160672 : c = real_norm(gel(x,i));
2281 160672 : s = s? gadd(s, c): c;
2282 : }
2283 202817 : for (; i < l; i++)
2284 : {
2285 131872 : c = complex_norm(gel(x,i));
2286 131869 : t = t? gadd(t, c): c;
2287 : }
2288 70945 : if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
2289 70945 : return gerepileupto(av, s);
2290 : }
2291 : /* return N(x) */
2292 : GEN
2293 18933642 : embed_norm(GEN x, long r1)
2294 : {
2295 18933642 : pari_sp av = avma;
2296 18933642 : long i, l = lg(x);
2297 18933642 : GEN c, s = NULL, t = NULL;
2298 18933642 : if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
2299 47175155 : for (i = 1; i <= r1; i++)
2300 : {
2301 28370224 : c = gel(x,i);
2302 28370224 : s = s? gmul(s, c): c;
2303 : }
2304 53569080 : for (; i < l; i++)
2305 : {
2306 34765127 : c = complex_norm(gel(x,i));
2307 34763264 : t = t? gmul(t, c): c;
2308 : }
2309 18803953 : if (t) s = s? gmul(s,t): t;
2310 18803962 : return gerepileupto(av, s);
2311 : }
2312 :
2313 : typedef struct {
2314 : long r1, v, prec;
2315 : GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
2316 : GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
2317 : GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
2318 : GEN bound; /* T2 norm of the polynomial defining nf */
2319 : long expo_best_disc; /* expo(disc(x)), best generator so far */
2320 : } CG_data;
2321 :
2322 : /* characteristic pol of x (given by embeddings) */
2323 : static GEN
2324 260240 : get_pol(CG_data *d, GEN x)
2325 : {
2326 : long e;
2327 260240 : GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
2328 260240 : return (e > -5)? NULL: g;
2329 : }
2330 :
2331 : /* characteristic pol of x (given as vector on (w_i)) */
2332 : static GEN
2333 98080 : get_polchar(CG_data *d, GEN x)
2334 98080 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
2335 :
2336 : /* Choose a canonical polynomial in the pair { Pmin_a, Pmin_{-a} }, i.e.
2337 : * { z(X), (-1)^(deg z) z(-Z) } and keeping the smallest wrt cmpii_polred
2338 : * Either leave z alone (return 1) or set z <- (-1)^n z(-X). In place. */
2339 : int
2340 89019 : ZX_canon_neg(GEN z)
2341 : {
2342 : long i, s;
2343 126688 : for (i = lg(z)-2; i >= 2; i -= 2)
2344 : { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
2345 123167 : s = signe(gel(z,i));
2346 123167 : if (!s) continue;
2347 : /* non trivial */
2348 85498 : if (s < 0) break; /* z(X) < (-1)^n z(-X) */
2349 :
2350 211373 : for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
2351 36438 : return 1;
2352 : }
2353 52581 : return 0;
2354 : }
2355 : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
2356 : * Return NULL on failure: discriminant too large or non primitive */
2357 : static GEN
2358 136081 : try_polmin(CG_data *d, nfmaxord_t *S, GEN v, long flag, GEN *ai)
2359 : {
2360 136081 : const long best = flag & nf_ABSOLUTE;
2361 : long ed;
2362 136081 : pari_sp av = avma;
2363 : GEN g;
2364 136081 : if (best)
2365 : {
2366 135206 : ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
2367 135208 : set_avma(av); if (d->expo_best_disc < ed) return NULL;
2368 : }
2369 : else
2370 875 : ed = 0;
2371 82247 : g = get_pol(d, v);
2372 : /* accuracy too low, compute algebraically */
2373 82247 : if (!g) { set_avma(av); g = ZXQ_charpoly(*ai, S->T, varn(S->T)); }
2374 82247 : g = ZX_radical(g);
2375 82247 : if (best && degpol(g) != degpol(S->T)) return gc_NULL(av);
2376 30611 : g = gerepilecopy(av, g);
2377 30612 : d->expo_best_disc = ed;
2378 30612 : if (flag & nf_ORIG)
2379 : {
2380 26930 : if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
2381 26929 : if (!isint1(S->unscale)) *ai = RgX_unscale(*ai, S->unscale);
2382 : }
2383 : else
2384 3682 : (void)ZX_canon_neg(g);
2385 30611 : if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
2386 30611 : return g;
2387 : }
2388 :
2389 : /* does x generate the correct field ? */
2390 : static GEN
2391 98080 : chk_gen(void *data, GEN x)
2392 : {
2393 98080 : pari_sp av = avma, av1;
2394 98080 : GEN h, g = get_polchar((CG_data*)data,x);
2395 98081 : if (!g) pari_err_PREC("chk_gen");
2396 98081 : av1 = avma;
2397 98081 : h = ZX_gcd(g, ZX_deriv(g));
2398 98081 : if (degpol(h)) return gc_NULL(av);
2399 55246 : if (DEBUGLEVEL>3) err_printf(" generator: %Ps\n",g);
2400 55246 : set_avma(av1); return gerepileupto(av, g);
2401 : }
2402 :
2403 : static long
2404 32165 : chk_gen_prec(long N, long bit)
2405 32165 : { return prec_fix(nbits2prec(10 + (long)log2((double)N) + bit)); }
2406 :
2407 : /* v = [P,A] two vectors (of ZX and ZV resp.) of same length; remove duplicate
2408 : * polynomials in P, updating A, in place. Among elements having the same
2409 : * characteristic pol, choose the smallest according to ZV_abscmp */
2410 : static void
2411 13825 : remove_duplicates(GEN v)
2412 : {
2413 13825 : GEN x, a, P = gel(v,1), A = gel(v,2);
2414 13825 : long k, i, l = lg(P);
2415 13825 : pari_sp av = avma;
2416 :
2417 13825 : if (l < 2) return;
2418 13825 : (void)sort_factor_pol(mkvec2(P, A), cmpii);
2419 13825 : x = gel(P,1); a = gel(A,1);
2420 53858 : for (k=1,i=2; i<l; i++)
2421 40033 : if (ZX_equal(gel(P,i), x))
2422 : {
2423 19677 : if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
2424 : }
2425 : else
2426 : {
2427 20356 : gel(A,k) = a;
2428 20356 : gel(P,k) = x;
2429 20356 : k++;
2430 20356 : x = gel(P,i); a = gel(A,i);
2431 : }
2432 13825 : l = k+1;
2433 13825 : gel(A,k) = a; setlg(A,l);
2434 13825 : gel(P,k) = x; setlg(P,l); set_avma(av);
2435 : }
2436 :
2437 : static void
2438 18137 : polred_init(nfmaxord_t *S, nffp_t *F, CG_data *d)
2439 : {
2440 18137 : long e, prec, n = degpol(S->T);
2441 : double log2rho;
2442 : GEN ro;
2443 18137 : set_LLL_basis(S, &ro, 0, 0.9999);
2444 : /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
2445 18130 : log2rho = ro ? (double)gexpo(ro): fujiwara_bound(S->T);
2446 18130 : e = n * (long)(log2rho + log2((double)n)) + 1;
2447 18130 : if (e < 0) e = 0; /* can occur if n = 1 */
2448 18130 : prec = chk_gen_prec(n, e);
2449 18130 : nffp_init(F,S,prec);
2450 18130 : F->ro = ro;
2451 18130 : make_M_G(F, 1);
2452 :
2453 18130 : d->v = varn(S->T);
2454 18130 : d->expo_best_disc = -1;
2455 18130 : d->ZKembed = NULL;
2456 18130 : d->M = NULL;
2457 18130 : d->u = NULL;
2458 18130 : d->r1= S->r1;
2459 18130 : }
2460 : static GEN
2461 14049 : findmindisc(GEN y)
2462 : {
2463 14049 : GEN x = gel(y,1), dx = NULL;
2464 14049 : long i, l = lg(y);
2465 35217 : for (i = 2; i < l; i++)
2466 : {
2467 21168 : GEN yi = gel(y,i);
2468 21168 : if (ZX_is_better(yi,x,&dx)) x = yi;
2469 : }
2470 14049 : return x;
2471 : }
2472 : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
2473 : * [ the best wrt discriminant ordering ], but keep all imprimitive
2474 : * polynomials */
2475 : static void
2476 4095 : filter(GEN y, GEN b, long n)
2477 : {
2478 : GEN x, a, dx;
2479 4095 : long i, k = 1, l = lg(y);
2480 4095 : a = x = dx = NULL;
2481 34770 : for (i = 1; i < l; i++)
2482 : {
2483 30675 : GEN yi = gel(y,i), ai = gel(b,i);
2484 30675 : if (degpol(yi) == n)
2485 : {
2486 30507 : if (!dx) dx = ZX_disc(yi); else if (!ZX_is_better(yi,x,&dx)) continue;
2487 5675 : x = yi; a = ai; continue;
2488 : }
2489 168 : gel(y,k) = yi;
2490 168 : gel(b,k) = ai; k++;
2491 : }
2492 4095 : if (dx)
2493 : {
2494 3991 : gel(y,k) = x;
2495 3991 : gel(b,k) = a; k++;
2496 : }
2497 4095 : setlg(y, k);
2498 4095 : setlg(b, k);
2499 4095 : }
2500 :
2501 : static GEN
2502 4130 : polred_aux(nfmaxord_t *S, GEN *pro, long flag)
2503 : { /* only keep polynomials of max degree and best discriminant */
2504 4130 : const long best = flag & nf_ABSOLUTE;
2505 4130 : const long orig = flag & nf_ORIG;
2506 4130 : GEN M, b, y, x = S->T;
2507 4130 : long maxi, i, j, k, v = varn(x), n = lg(S->basis)-1;
2508 : nffp_t F;
2509 : CG_data d;
2510 :
2511 4130 : if (n == 1)
2512 : {
2513 28 : if (!best)
2514 : {
2515 14 : GEN X = pol_x(v);
2516 14 : return orig? mkmat2(mkcol(X),mkcol(gen_1)): mkvec(X);
2517 : }
2518 : else
2519 14 : return orig? trivial_fact(): cgetg(1,t_VEC);
2520 : }
2521 :
2522 4102 : polred_init(S, &F, &d);
2523 4095 : if (pro) *pro = F.ro;
2524 4095 : M = F.M;
2525 4095 : if (best)
2526 : {
2527 4032 : if (!S->dT) S->dT = ZX_disc(S->T);
2528 4032 : d.expo_best_disc = expi(S->dT);
2529 : }
2530 :
2531 : /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
2532 4095 : y = cgetg(n*n + 1, t_VEC);
2533 4095 : b = cgetg(n*n + 1, t_COL);
2534 4095 : k = 1;
2535 4095 : if (!best) { gel(y,1) = pol_x(v); gel(b,1) = gen_0; k++; }
2536 26935 : for (i = 2; i <= n; i++)
2537 : {
2538 : GEN ch, ai;
2539 22840 : ai = gel(S->basis,i);
2540 22840 : ch = try_polmin(&d, S, gel(M,i), flag, &ai);
2541 22840 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2542 : }
2543 4095 : maxi = minss(n, 3);
2544 15995 : for (i = 1; i <= maxi; i++)
2545 68522 : for (j = i+1; j <= n; j++)
2546 : {
2547 : GEN ch, ai, v;
2548 56622 : ai = gadd(gel(S->basis,i), gel(S->basis,j));
2549 56623 : v = RgV_add(gel(M,i), gel(M,j));
2550 : /* defining polynomial for Q(w_i+w_j) */
2551 56621 : ch = try_polmin(&d, S, v, flag, &ai);
2552 56622 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2553 :
2554 56622 : ai = gsub(gel(S->basis,i), gel(S->basis,j));
2555 56620 : v = RgV_sub(gel(M,i), gel(M,j));
2556 : /* defining polynomial for Q(w_i-w_j) */
2557 56621 : ch = try_polmin(&d, S, v, flag, &ai);
2558 56622 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2559 : }
2560 4095 : setlg(y, k);
2561 4095 : setlg(b, k); filter(y, b, n);
2562 4095 : if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
2563 3416 : settyp(y, t_COL);
2564 3416 : (void)sort_factor_pol(mkmat2(y, b), cmpii);
2565 3416 : return mkmat2(b, y);
2566 : }
2567 :
2568 : /* FIXME: obsolete */
2569 : static GEN
2570 84 : Polred(GEN x, long flag, GEN fa)
2571 : {
2572 84 : pari_sp av = avma;
2573 : nfmaxord_t S;
2574 84 : if (fa)
2575 14 : nfinit_basic(&S, mkvec2(x,fa));
2576 : else
2577 70 : nfinit_basic_flag(&S, x, flag);
2578 77 : return gerepilecopy(av, polred_aux(&S, NULL, flag));
2579 : }
2580 :
2581 : /* finds "best" polynomial in polred_aux list, defaulting to S->T if none of
2582 : * them is primitive. *px a ZX, characteristic polynomial of Mod(*pb,S->T),
2583 : * *pdx its discriminant if pdx != NULL. Set *pro = polroots(S->T) */
2584 : static void
2585 4053 : polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
2586 : {
2587 4053 : GEN y, dx, x = S->T; /* default value */
2588 : long i, l;
2589 4053 : y = polred_aux(S, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
2590 4046 : dx = S->dT;
2591 4046 : if (pb)
2592 : {
2593 3402 : GEN a, b = deg1pol_shallow(S->unscale, gen_0, varn(x));
2594 3402 : a = gel(y,1); l = lg(a);
2595 3402 : y = gel(y,2);
2596 6693 : for (i=1; i<l; i++)
2597 : {
2598 3291 : GEN yi = gel(y,i);
2599 3291 : if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
2600 : }
2601 3402 : *pb = b;
2602 : }
2603 : else
2604 : {
2605 644 : l = lg(y);
2606 1281 : for (i=1; i<l; i++)
2607 : {
2608 637 : GEN yi = gel(y,i);
2609 637 : if (ZX_is_better(yi,x,&dx)) x = yi;
2610 : }
2611 : }
2612 4046 : if (pdx) { if (!dx) dx = ZX_disc(x); *pdx = dx; }
2613 4046 : *px = x;
2614 4046 : }
2615 : static GEN
2616 3871 : polredbest_i(GEN T, long flag)
2617 : {
2618 : nfmaxord_t S;
2619 : GEN a;
2620 3871 : nfinit_basic_partial(&S, T);
2621 3871 : polredbest_aux(&S, NULL, &T, NULL, flag? &a: NULL);
2622 3864 : if (flag == 2)
2623 350 : T = mkvec2(T, a);
2624 3514 : else if (flag == 1)
2625 : {
2626 2870 : GEN b = (S.T0 == T)? pol_x(varn(T)): QXQ_reverse(a, S.T0);
2627 : /* charpoly(Mod(a,T0)) = T; charpoly(Mod(b,T)) = S.x */
2628 2870 : if (degpol(T) == 1) b = grem(b,T);
2629 2870 : T = mkvec2(T, mkpolmod(b,T));
2630 : }
2631 3864 : return T;
2632 : }
2633 : GEN
2634 3514 : polredbest(GEN T, long flag)
2635 : {
2636 3514 : pari_sp av = avma;
2637 3514 : if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
2638 3514 : return gerepilecopy(av, polredbest_i(T, flag));
2639 : }
2640 : /* DEPRECATED: backward compatibility */
2641 : GEN
2642 70 : polred0(GEN x, long flag, GEN fa)
2643 : {
2644 70 : long fl = 0;
2645 70 : if (flag & 1) fl |= nf_PARTIALFACT;
2646 70 : if (flag & 2) fl |= nf_ORIG;
2647 70 : return Polred(x, fl, fa);
2648 : }
2649 :
2650 : GEN
2651 21 : polredord(GEN x)
2652 : {
2653 21 : pari_sp av = avma;
2654 : GEN v, lt;
2655 : long i, n, vx;
2656 :
2657 21 : if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
2658 21 : x = Q_primpart(x); RgX_check_ZX(x,"polredord");
2659 21 : n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
2660 21 : if (n == 1) return gerepilecopy(av, mkvec(x));
2661 14 : lt = leading_coeff(x); vx = varn(x);
2662 14 : if (is_pm1(lt))
2663 : {
2664 7 : if (signe(lt) < 0) x = ZX_neg(x);
2665 7 : v = pol_x_powers(n, vx);
2666 : }
2667 : else
2668 : { GEN L;
2669 : /* basis for Dedekind order */
2670 7 : v = cgetg(n+1, t_VEC);
2671 7 : gel(v,1) = scalarpol_shallow(lt, vx);
2672 14 : for (i = 2; i <= n; i++)
2673 7 : gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
2674 7 : gel(v,1) = pol_1(vx);
2675 7 : x = ZX_Q_normalize(x, &L);
2676 7 : v = gsubst(v, vx, monomial(ginv(L),1,vx));
2677 14 : for (i=2; i <= n; i++)
2678 7 : if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = pol_xn(i-1, vx);
2679 : }
2680 14 : return gerepileupto(av, polred(mkvec2(x, v)));
2681 : }
2682 :
2683 : GEN
2684 14 : polred(GEN x) { return Polred(x, 0, NULL); }
2685 : GEN
2686 0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
2687 : GEN
2688 0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
2689 : GEN
2690 0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
2691 : GEN
2692 0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
2693 : GEN
2694 0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
2695 :
2696 : /********************************************************************/
2697 : /** **/
2698 : /** POLREDABS **/
2699 : /** **/
2700 : /********************************************************************/
2701 : /* set V[k] := matrix of multiplication by nk.zk[k] */
2702 : static GEN
2703 22302 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
2704 : {
2705 22302 : GEN v, Mk = cgetg(N+1, t_MAT);
2706 : long i, e;
2707 34742 : for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
2708 144512 : for ( ; i <=N; i++)
2709 : {
2710 122210 : v = vecmul(gel(M,k), gel(M,i));
2711 122210 : v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
2712 122209 : gel(Mk,i) = grndtoi(v, &e);
2713 122210 : if (e > -5) return NULL;
2714 : }
2715 22302 : gel(V,k) = Mk; return Mk;
2716 : }
2717 :
2718 : static GEN
2719 12689 : ZM_image_shallow(GEN M, long *pr)
2720 : {
2721 : long j, k, r;
2722 12689 : GEN y, d = ZM_pivots(M, &k);
2723 12689 : r = lg(M)-1 - k;
2724 12689 : y = cgetg(r+1,t_MAT);
2725 55729 : for (j=k=1; j<=r; k++)
2726 43040 : if (d[k]) gel(y,j++) = gel(M,k);
2727 12689 : *pr = r; return y;
2728 : }
2729 :
2730 : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
2731 : * Q from algo 2.7.6] */
2732 : static GEN
2733 14042 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
2734 : {
2735 14042 : CG_data *d = (CG_data*)chk->data;
2736 : GEN P, V, D, inv, bound, S, M;
2737 14042 : long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
2738 14042 : long i, j, prec, firstprim = 0, skipfirst = 0;
2739 : pari_sp av;
2740 :
2741 14042 : d->u = U;
2742 14042 : d->ZKembed = M = RgM_mul(d->M, U);
2743 :
2744 14042 : av = avma; bound = d->bound;
2745 14042 : D = cgetg(N+1, t_VECSMALL);
2746 93933 : for (i = 1; i <= N; i++)
2747 : {
2748 79898 : pari_sp av2 = avma;
2749 79898 : P = get_pol(d, gel(M,i));
2750 79898 : if (!P) pari_err_PREC("chk_gen_init");
2751 79891 : P = gerepilecopy(av2, ZX_radical(P));
2752 79891 : D[i] = degpol(P);
2753 79891 : if (D[i] == N)
2754 : { /* primitive element */
2755 56633 : GEN B = embed_T2(gel(M,i), r1);
2756 56632 : if (!firstprim) firstprim = i; /* index of first primitive element */
2757 56632 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2758 56632 : if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
2759 : }
2760 : else
2761 : {
2762 23258 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
2763 23258 : if (firstprim)
2764 : { /* cycle basis vectors so that primitive elements come last */
2765 2660 : GEN u = d->u, e = M;
2766 2660 : GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
2767 2660 : long tS = D[i];
2768 8208 : for (j = i; j > firstprim; j--)
2769 : {
2770 5548 : u[j] = u[j-1];
2771 5548 : e[j] = e[j-1];
2772 5548 : R[j] = R[j-1];
2773 5548 : D[j] = D[j-1];
2774 : }
2775 2660 : gel(u,firstprim) = tu;
2776 2660 : gel(e,firstprim) = te;
2777 2660 : gel(R,firstprim) = tR;
2778 2660 : D[firstprim] = tS; firstprim++;
2779 : }
2780 : }
2781 : }
2782 14035 : if (!firstprim)
2783 : { /* try (a little) to find primitive elements to improve bound */
2784 28 : GEN x = cgetg(N+1, t_VECSMALL);
2785 28 : if (DEBUGLEVEL>1)
2786 0 : err_printf("chk_gen_init: difficult field, trying random elements\n");
2787 308 : for (i = 0; i < 10; i++)
2788 : {
2789 : GEN e, B;
2790 3010 : for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
2791 280 : e = RgM_zc_mul(M, x);
2792 280 : B = embed_T2(e, r1);
2793 280 : if (gcmp(B,bound) >= 0) continue;
2794 14 : P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
2795 14 : if (!ZX_is_squarefree(P)) continue;
2796 14 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2797 14 : bound = B ;
2798 : }
2799 : }
2800 :
2801 14035 : if (firstprim != 1)
2802 : {
2803 13756 : inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
2804 13756 : V = gel(inv,1);
2805 53308 : for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
2806 : /* V corresponds to 1_Z */
2807 13756 : V = grndtoi(V, &j);
2808 13756 : if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
2809 13756 : S = mkmat(V); /* 1 */
2810 :
2811 13756 : V = cgetg(N+1, t_VEC);
2812 35286 : for (i = 1; i <= N; i++,skipfirst++)
2813 : { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
2814 : GEN Mx, M2;
2815 35286 : long j, k, h, rkM, dP = D[i];
2816 :
2817 35286 : if (dP == N) break; /* primitive */
2818 22302 : Mx = set_mulid(V, M, inv, r1, r2, N, i);
2819 22302 : if (!Mx) break; /* prec. problem. Stop */
2820 22302 : if (dP == 1) continue;
2821 9266 : rkM = lg(S)-1;
2822 9266 : M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
2823 9266 : gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
2824 9266 : k = 2;
2825 18535 : for (h = 1; h < dP; h++)
2826 : {
2827 : long r; /* add to M2 the elts of S * nf.zk[i] */
2828 38073 : for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
2829 12689 : setlg(M2, k); k = 1;
2830 12689 : S = ZM_image_shallow(shallowconcat(S,M2), &r);
2831 13461 : if (r == rkM) break;
2832 10041 : if (r > rkM)
2833 : {
2834 10041 : rkM = r;
2835 10041 : if (rkM == N) break;
2836 : }
2837 : }
2838 9266 : if (rkM == N) break;
2839 : /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
2840 : }
2841 : }
2842 : /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
2843 14035 : chk->skipfirst = skipfirst;
2844 14035 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
2845 :
2846 : /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
2847 14035 : bound = gerepileuptoleaf(av, bound);
2848 14035 : prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
2849 14035 : if (DEBUGLEVEL)
2850 0 : err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
2851 14035 : if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
2852 14035 : if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
2853 14035 : return bound;
2854 : }
2855 :
2856 : static GEN
2857 14049 : polredabs_i(GEN x, nfmaxord_t *S, GEN *u, long flag)
2858 : {
2859 14049 : FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
2860 : nffp_t F;
2861 : CG_data d;
2862 : GEN v, y, a;
2863 : long i, l;
2864 :
2865 14049 : nfinit_basic_flag(S, x, flag);
2866 14049 : x = S->T0;
2867 14049 : if (degpol(x) == 1)
2868 : {
2869 14 : long vx = varn(x);
2870 14 : *u = NULL;
2871 14 : return mkvec2(mkvec( pol_x(vx) ),
2872 14 : mkvec( deg1pol_shallow(gen_1, negi(gel(S->T,2)), vx) ));
2873 : }
2874 14035 : chk.data = (void*)&d;
2875 14035 : polred_init(S, &F, &d);
2876 14035 : d.bound = embed_T2(F.ro, d.r1);
2877 14035 : if (realprec(d.bound) > F.prec) d.bound = rtor(d.bound, F.prec);
2878 : for (;;)
2879 21 : {
2880 14056 : GEN R = R_from_QR(F.G, F.prec);
2881 14056 : if (R)
2882 : {
2883 14042 : d.prec = F.prec;
2884 14042 : d.M = F.M;
2885 14042 : v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
2886 14042 : if (v) break;
2887 : }
2888 21 : F.prec = precdbl(F.prec);
2889 21 : F.ro = NULL;
2890 21 : make_M_G(&F, 1);
2891 21 : if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",F.prec);
2892 : }
2893 14035 : y = gel(v,1);
2894 14035 : a = gel(v,2); l = lg(a);
2895 68915 : for (i = 1; i < l; i++) /* normalize wrt z -> -z */
2896 54880 : if (ZX_canon_neg(gel(y,i)) && (flag & (nf_ORIG|nf_RAW)))
2897 546 : gel(a,i) = ZC_neg(gel(a,i));
2898 14035 : *u = d.u; return v;
2899 : }
2900 :
2901 : GEN
2902 13825 : polredabs0(GEN x, long flag)
2903 : {
2904 13825 : pari_sp av = avma;
2905 : GEN Y, A, u, v;
2906 : nfmaxord_t S;
2907 : long i, l;
2908 :
2909 13825 : v = polredabs_i(x, &S, &u, flag);
2910 13825 : remove_duplicates(v);
2911 13825 : Y = gel(v,1);
2912 13825 : A = gel(v,2);
2913 13825 : l = lg(A); if (l == 1) pari_err_BUG("polredabs (missing vector)");
2914 13825 : if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
2915 13825 : if (!(flag & nf_ALL))
2916 : {
2917 13825 : GEN y = findmindisc(Y);
2918 28084 : for (i = 1; i < l; i++)
2919 28084 : if (ZX_equal(gel(Y,i), y)) break;
2920 13825 : Y = mkvec(gel(Y,i));
2921 13825 : A = mkvec(gel(A,i)); l = 2;
2922 : }
2923 14091 : if (flag & (nf_RAW|nf_ORIG)) for (i = 1; i < l; i++)
2924 : {
2925 266 : GEN y = gel(Y,i), a = gel(A,i);
2926 266 : if (u) a = RgV_RgC_mul(S.basis, ZM_ZC_mul(u, a));
2927 266 : if (flag & nf_ORIG)
2928 : {
2929 259 : a = QXQ_reverse(a, S.T);
2930 259 : if (!isint1(S.unscale)) a = gdiv(a, S.unscale); /* not RgX_Rg_div */
2931 259 : a = mkpolmod(a,y);
2932 : }
2933 266 : gel(Y,i) = mkvec2(y, a);
2934 : }
2935 13825 : return gerepilecopy(av, (flag & nf_ALL)? Y: gel(Y,1));
2936 : }
2937 :
2938 : GEN
2939 0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
2940 : GEN
2941 13377 : polredabs(GEN x) { return polredabs0(x,0); }
2942 : GEN
2943 0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
2944 :
2945 : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
2946 : * flag & nf_ORIG: + element (base change)
2947 : * flag & nf_ABSOLUTE: absolute polynomial */
2948 : static GEN
2949 588 : rnfpolred_i(GEN nf, GEN R, long flag, long best)
2950 : {
2951 588 : const char *f = best? "rnfpolredbest": "rnfpolredabs";
2952 588 : const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
2953 588 : GEN listP = NULL, red, pol, A, P, T, rnfeq;
2954 588 : pari_sp av = avma;
2955 :
2956 588 : if (typ(R) == t_VEC) {
2957 14 : if (lg(R) != 3) pari_err_TYPE(f,R);
2958 14 : listP = gel(R,2);
2959 14 : R = gel(R,1);
2960 : }
2961 588 : if (typ(R) != t_POL) pari_err_TYPE(f,R);
2962 588 : nf = checknf(nf);
2963 588 : T = nf_get_pol(nf);
2964 588 : R = RgX_nffix(f, T, R, 0);
2965 588 : if (best || (flag & nf_PARTIALFACT))
2966 : {
2967 364 : rnfeq = abs? nf_rnfeq(nf, R): nf_rnfeqsimple(nf, R);
2968 364 : pol = gel(rnfeq,1);
2969 364 : if (listP) pol = mkvec2(pol, listP);
2970 357 : red = best? polredbest_i(pol, abs? 1: 2)
2971 364 : : polredabs0(pol, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
2972 364 : P = gel(red,1);
2973 364 : A = gel(red,2);
2974 : }
2975 : else
2976 : {
2977 : nfmaxord_t S;
2978 : GEN rnf, u, v, y, a;
2979 : long i, j, l;
2980 : pari_timer ti;
2981 224 : if (DEBUGLEVEL>1) timer_start(&ti);
2982 224 : rnf = rnfinit(nf, R);
2983 224 : rnfeq = rnf_get_map(rnf);
2984 224 : pol = rnf_zkabs(rnf);
2985 224 : if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
2986 224 : v = polredabs_i(pol, &S, &u, nf_ORIG);
2987 224 : pol = gel(pol,1);
2988 224 : y = gel(v,1); P = findmindisc(y);
2989 224 : a = gel(v,2);
2990 224 : l = lg(y); A = cgetg(l, t_VEC);
2991 1260 : for (i = j = 1; i < l; i++)
2992 1036 : if (ZX_equal(gel(y,i),P))
2993 : {
2994 924 : GEN t = gel(a,i);
2995 924 : if (u) t = RgV_RgC_mul(S.basis, ZM_ZC_mul(u,t));
2996 924 : gel(A,j++) = t;
2997 : }
2998 224 : setlg(A,j); /* mod(A[i], pol) are all roots of P in Q[X]/(pol) */
2999 : }
3000 588 : if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
3001 588 : if (flag & nf_ABSOLUTE)
3002 : {
3003 14 : if (flag & nf_ORIG)
3004 : {
3005 7 : GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
3006 7 : GEN k = gel(rnfeq,3); /* Mod(variable(R),R) + k*a root of pol */
3007 7 : if (typ(A) == t_VEC) A = gel(A,1); /* any root will do */
3008 7 : a = RgX_RgXQ_eval(a, lift_shallow(A), P); /* Mod(a, P) root of T */
3009 7 : P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
3010 : }
3011 14 : return gerepilecopy(av, P);
3012 : }
3013 574 : if (typ(A) != t_VEC)
3014 : {
3015 357 : A = eltabstorel_lift(rnfeq, A);
3016 357 : P = lift_if_rational( RgXQ_charpoly(A, R, varn(R)) );
3017 : }
3018 : else
3019 : { /* canonical factor */
3020 217 : long i, l = lg(A), v = varn(R);
3021 217 : GEN besta = NULL;
3022 1099 : for (i = 1; i < l; i++)
3023 : {
3024 882 : GEN a = eltabstorel_lift(rnfeq, gel(A,i));
3025 882 : GEN p = lift_if_rational( RgXQ_charpoly(a, R, v) );
3026 882 : if (i == 1 || cmp_universal(p, P) < 0) { P = p; besta = a; }
3027 : }
3028 217 : A = besta;
3029 : }
3030 574 : if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,R),P));
3031 574 : return gerepilecopy(av, P);
3032 : }
3033 : GEN
3034 231 : rnfpolredabs(GEN nf, GEN R, long flag)
3035 231 : { return rnfpolred_i(nf,R,flag, 0); }
3036 : GEN
3037 357 : rnfpolredbest(GEN nf, GEN R, long flag)
3038 : {
3039 357 : if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
3040 357 : return rnfpolred_i(nf,R,flag, 1);
3041 : }
|
