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 426232 : v13checkrnf(GEN v)
30 426232 : { return typ(gel(v,6)) == t_VEC; }
31 : static int
32 50666 : rawcheckbnf(GEN v) { return typ(v)==t_VEC && lg(v)==11; }
33 : static int
34 49791 : 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 41482 : 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 924 : 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 217203 : 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 665 : v9checkgal(GEN v)
47 665 : { GEN x = gel(v,2); return typ(x) == t_VEC && lg(x) == 4; }
48 :
49 : int
50 437246 : checkrnf_i(GEN rnf)
51 437246 : { return (typ(rnf)==t_VEC && lg(rnf)==13 && v13checkrnf(rnf)); }
52 :
53 : void
54 425525 : checkrnf(GEN rnf)
55 425525 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
56 :
57 : GEN
58 5426990 : checkbnf_i(GEN X)
59 : {
60 5426990 : if (typ(X) == t_VEC)
61 5426815 : switch (lg(X))
62 : {
63 5410297 : case 11:
64 5410297 : if (typ(gel(X,6)) != t_INT) return NULL; /* pre-2.2.4 format */
65 5410297 : if (lg(gel(X,10)) != 4) return NULL; /* pre-2.8.1 format */
66 5410297 : return X;
67 15995 : case 7: return checkbnf_i(bnr_get_bnf(X));
68 : }
69 698 : return NULL;
70 : }
71 :
72 : GEN
73 92293424 : checknf_i(GEN X)
74 : {
75 92293424 : if (typ(X)==t_VEC)
76 91643861 : switch(lg(X))
77 : {
78 91143526 : case 10: return X;
79 484491 : case 11: return checknf_i(bnf_get_nf(X));
80 2891 : case 7: return checknf_i(bnr_get_bnf(X));
81 8883 : case 3: if (typ(gel(X,2)) == t_POLMOD) return checknf_i(gel(X,1));
82 : }
83 662509 : return NULL;
84 : }
85 :
86 : GEN
87 3386568 : checkbnf(GEN x)
88 : {
89 3386568 : GEN bnf = checkbnf_i(x);
90 3386565 : if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
91 3386545 : return bnf;
92 : }
93 :
94 : GEN
95 89115130 : checknf(GEN x)
96 : {
97 89115130 : GEN nf = checknf_i(x);
98 89115258 : if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
99 89115213 : return nf;
100 : }
101 :
102 : GEN
103 2022650 : checkbnr_i(GEN bnr)
104 : {
105 2022650 : if (typ(bnr)!=t_VEC || lg(bnr)!=7 || !checkbnf_i(bnr_get_bnf(bnr)))
106 84 : return NULL;
107 2022564 : return bnr;
108 : }
109 : void
110 1839398 : checkbnr(GEN bnr)
111 : {
112 1839398 : if (!checkbnr_i(bnr))
113 14 : pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
114 1839378 : }
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 1791348 : checkbid_i(GEN bid)
125 : {
126 : GEN f;
127 1791348 : if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(bid_get_U(bid)) != t_VEC)
128 269909 : return NULL;
129 1521439 : f = bid_get_mod(bid);
130 1521439 : if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
131 1521439 : return bid;
132 : }
133 : void
134 1521439 : checkbid(GEN bid)
135 : {
136 1521439 : if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
137 1521432 : }
138 : void
139 24115 : checkabgrp(GEN v)
140 : {
141 24115 : if (typ(v) == t_VEC) switch(lg(v))
142 : {
143 21714 : case 4: if (typ(gel(v,3)) != t_VEC) break;
144 24115 : case 3: if (typ(gel(v,2)) != t_VEC) break;
145 24087 : if (typ(gel(v,1)) != t_INT) break;
146 24087 : return;/*OK*/
147 0 : default: break;
148 : }
149 28 : pari_err_TYPE("checkabgrp",v);
150 : }
151 :
152 : GEN
153 619021 : checknfelt_mod(GEN nf, GEN x, const char *s)
154 : {
155 619021 : GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
156 619018 : if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
157 618883 : return a;
158 : }
159 :
160 : int
161 6793 : check_ZKmodule_i(GEN M)
162 : {
163 6793 : return (typ(M) ==t_VEC && lg(M) >= 3
164 6793 : && typ(gel(M,1)) == t_MAT
165 6793 : && typ(gel(M,2)) == t_VEC
166 13586 : && lgcols(M) == lg(gel(M,2)));
167 : }
168 : void
169 6737 : check_ZKmodule(GEN M, const char *s)
170 6737 : { if (!check_ZKmodule_i(M)) pari_err_TYPE(s, M); }
171 :
172 : static long
173 126476 : typv6(GEN x)
174 : {
175 126476 : if (typ(gel(x,1)) == t_VEC && lg(gel(x,3)) == 3)
176 : {
177 23597 : GEN t = gel(x,3);
178 23597 : if (typ(t) != t_VEC) return typ_NULL;
179 23597 : t = gel(x,5);
180 23597 : switch(typ(gel(x,5)))
181 : {
182 455 : case t_VEC: return typ_BID;
183 23142 : case t_MAT: return typ_BIDZ;
184 0 : default: return typ_NULL;
185 : }
186 : }
187 102879 : 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 28504 : get_bnf(GEN x, long *t)
194 : {
195 28504 : 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 27881 : case t_VEC:
200 27881 : switch(lg(x))
201 : {
202 2247 : case 5: if (typ(gel(x,1)) != t_INT) break;
203 2191 : *t = typ_QUA; return NULL;
204 14721 : case 6: *t = typv6(x); return NULL;
205 2380 : case 7: *t = typ_BNR;
206 2380 : x = bnr_get_bnf(x);
207 2380 : if (!rawcheckbnf(x)) break;
208 2324 : return x;
209 84 : case 9:
210 84 : if (!v9checkgal(x)) break;
211 84 : *t = typ_GAL; return NULL;
212 602 : case 10:
213 602 : if (!v10checknf(x)) break;
214 602 : *t = typ_NF; return NULL;
215 476 : case 11:
216 476 : if (!v11checkbnf(x)) break;
217 476 : *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 116564 : get_nf(GEN x, long *t)
234 : {
235 116564 : 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 113575 : case t_VEC:
240 113575 : 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 99554 : 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 574 : case 9:
254 574 : if (!v9checkgal(x)) break;
255 574 : *t = typ_GAL; return NULL;
256 1617 : case 10:
257 1617 : if (!v10checknf(x)) break;
258 1617 : *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 420 : case 13:
263 420 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_nf(x); }
264 392 : if (!v13checkrnf(x)) break;
265 392 : *t = typ_RNF; return NULL;
266 3486 : case 17: *t = typ_ELL; return NULL;
267 : }
268 399 : break;
269 147 : 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 2842 : *t = typ_NULL; return NULL;
275 : }
276 :
277 : long
278 268590 : nftyp(GEN x)
279 : {
280 268590 : switch(typ(x))
281 : {
282 28 : case t_POL : return typ_POL;
283 7 : case t_QUAD: return typ_Q;
284 268548 : case t_VEC:
285 268548 : switch(lg(x))
286 : {
287 259 : case 13:
288 259 : if (v13checkgchar(x)) return typ_GCHAR;
289 252 : if (!v13checkrnf(x)) break;
290 245 : return typ_RNF;
291 214984 : case 10:
292 214984 : if (!v10checknf(x)) break;
293 214977 : return typ_NF;
294 259 : case 11:
295 259 : if (!v11checkbnf(x)) break;
296 259 : return typ_BNF;
297 40761 : case 7:
298 40761 : x = bnr_get_bnf(x);
299 40761 : if (!rawcheckbnf(x) || !v11checkbnf(x)) break;
300 40747 : return typ_BNR;
301 12201 : case 6:
302 12201 : 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 : }
309 112 : return typ_NULL;
310 : }
311 :
312 : /*************************************************************************/
313 : /** **/
314 : /** GALOIS GROUP **/
315 : /** **/
316 : /*************************************************************************/
317 :
318 : GEN
319 7627 : tschirnhaus(GEN x)
320 : {
321 7627 : pari_sp av = avma, av2;
322 7627 : long a, v = varn(x);
323 7627 : GEN u, y = cgetg(5,t_POL);
324 :
325 7627 : if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
326 7627 : if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
327 7627 : if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
328 7627 : y[1] = evalsigne(1)|evalvarn(0);
329 : do
330 : {
331 8287 : a = random_bits(2); if (a==0) a = 1; gel(y,4) = stoi(a);
332 8287 : a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
333 8287 : a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
334 8287 : u = RgXQ_charpoly(y,x,v); av2 = avma;
335 : }
336 8287 : while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
337 7627 : if (DEBUGLEVEL>1)
338 0 : err_printf("Tschirnhaus transform. New pol: %Ps",u);
339 7627 : set_avma(av2); return gc_upto(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 845876 : ZX_Z_normalize(GEN pol, GEN *ptk)
347 : {
348 845876 : long i,j, sk, n = degpol(pol); /* > 0 */
349 : GEN k, fa, P, E, a, POL;
350 :
351 845869 : if (ptk) *ptk = gen_1;
352 845869 : if (!n) return pol;
353 845862 : a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
354 1627888 : for (i = n-2; i >= 0; i--)
355 : {
356 1371734 : k = gcdii(k, gel(a,i));
357 1371606 : if (is_pm1(k)) return pol;
358 : }
359 256154 : sk = signe(k);
360 256154 : if (!sk) return pol; /* monomial! */
361 254068 : fa = absZ_factor_limit(k, 0); k = gen_1;
362 254073 : P = gel(fa,1);
363 254073 : E = gel(fa,2);
364 254073 : POL = leafcopy(pol); a = POL+2;
365 547088 : for (i = lg(P)-1; i > 0; i--)
366 : {
367 293010 : GEN p = gel(P,i), pv, pvj;
368 293010 : long vmin = itos(gel(E,i));
369 : /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
370 1210789 : for (j=n-1; j>=0; j--)
371 : {
372 : long v;
373 917780 : if (!signe(gel(a,j))) continue;
374 789263 : v = Z_pval(gel(a,j), p) / (n - j);
375 789261 : if (v < vmin) vmin = v;
376 : }
377 293009 : if (!vmin) continue;
378 19717 : pvj = pv = powiu(p,vmin); k = mulii(k, pv);
379 : /* a[j] /= p^(v*(n-j)) */
380 79433 : for (j=n-1; j>=0; j--)
381 : {
382 59716 : if (j < n-1) pvj = mulii(pvj, pv);
383 59716 : gel(a,j) = diviiexact(gel(a,j), pvj);
384 : }
385 : }
386 254078 : if (ptk) *ptk = k;
387 254078 : 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 843185 : ZX_primitive_to_monic(GEN pol, GEN *pL)
395 : {
396 843185 : long i,j, n = degpol(pol);
397 843183 : GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
398 :
399 843178 : if (is_pm1(lc))
400 : {
401 840289 : if (pL) *pL = gen_1;
402 840289 : return signe(lc) < 0? ZX_neg(pol): pol;
403 : }
404 2891 : if (signe(lc) < 0)
405 35 : POL = ZX_neg(pol);
406 : else
407 2856 : POL = leafcopy(pol);
408 2891 : a = POL+2; lc = gel(a,n);
409 2891 : fa = absZ_factor_limit(lc,0); L = gen_1;
410 2891 : P = gel(fa,1);
411 2891 : E = gel(fa,2);
412 5971 : for (i = lg(P)-1; i > 0; i--)
413 : {
414 3080 : GEN p = gel(P,i), pk, pku;
415 3080 : long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
416 :
417 3080 : 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 17311 : for (j=n-1; j>0; j--)
420 : {
421 14231 : if (!signe(gel(a,j))) continue;
422 8834 : v = Z_pval(gel(a,j), p);
423 8932 : while (v + d < k * j) { k++; d += n; }
424 : }
425 3080 : pk = powiu(p,k); j0 = d/k;
426 3080 : L = mulii(L, pk);
427 :
428 3080 : pku = powiu(p,d - k*j0);
429 : /* a[j] *= p^(d - kj) */
430 7329 : for (j=j0; j>=0; j--)
431 : {
432 4249 : if (j < j0) pku = mulii(pku, pk);
433 4249 : gel(a,j) = mulii(gel(a,j), pku);
434 : }
435 3080 : j0++;
436 3080 : pku = powiu(p,k*j0 - d);
437 : /* a[j] /= p^(kj - d) */
438 19222 : for (j=j0; j<=n; j++)
439 : {
440 16142 : if (j > j0) pku = mulii(pku, pk);
441 16142 : gel(a,j) = diviiexact(gel(a,j), pku);
442 : }
443 : }
444 2891 : if (pL) *pL = L;
445 2891 : 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 842874 : ZX_Q_normalize(GEN pol, GEN *pL)
453 : {
454 842874 : GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
455 842875 : POL = ZX_Z_normalize(POL, pL);
456 842815 : if (pL) *pL = gdiv(lc, *pL);
457 842890 : return POL;
458 : }
459 :
460 : GEN
461 6677213 : ZX_Q_mul(GEN A, GEN z)
462 : {
463 6677213 : pari_sp av = avma;
464 6677213 : long i, l = lg(A);
465 : GEN d, n, Ad, B, u;
466 6677213 : if (typ(z)==t_INT) return ZX_Z_mul(A,z);
467 5621545 : n = gel(z, 1); d = gel(z, 2);
468 5621545 : Ad = RgX_to_RgC(FpX_red(A, d), l-2);
469 5621543 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
470 5621543 : B = cgetg(l, t_POL);
471 5621542 : B[1] = A[1];
472 5621542 : if (equali1(u))
473 : {
474 6082499 : for(i=2; i<l; i++)
475 3942411 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
476 : } else
477 : {
478 20532803 : for(i=2; i<l; i++)
479 : {
480 17051363 : GEN di = gcdii(gel(Ad, i-1), u);
481 17051358 : GEN ni = mulii(n, diviiexact(gel(A,i), di));
482 17051333 : if (equalii(d, di))
483 3876585 : gel(B, i) = ni;
484 : else
485 13174768 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
486 : }
487 : }
488 5621528 : return gc_GEN(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 gc_upto(av, res);
539 : }
540 :
541 : static GEN
542 96156 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
543 : {
544 96156 : pari_sp av = avma;
545 96156 : long i = lg(pol)-1;
546 : GEN res, den;
547 96156 : if (i==1) return gen_0;
548 96156 : pol = Q_remove_denom(pol, &den);
549 96156 : res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
550 399578 : for (i-- ; i>=2; i--)
551 303415 : res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
552 96163 : if (den) res = RgC_Rg_div(res, den);
553 96163 : return gc_upto(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 gc_upto(av, res);
571 : }
572 :
573 : /* compute s(x), not stack clean */
574 : static GEN
575 96381 : ZC_galoisapply(GEN nf, GEN s, GEN x)
576 : {
577 96381 : x = nf_to_scalar_or_alg(nf, x);
578 96379 : if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
579 96162 : 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 6276 : zk_galoisapplymod(GEN nf, GEN z, GEN S, GEN p)
585 : {
586 : GEN den, pe, pe1, denpe, R;
587 :
588 6276 : z = nf_to_scalar_or_alg(nf, z);
589 6276 : if (typ(z) != t_POL) return z;
590 6276 : 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 gc_upto(av,y);
648 : }
649 7 : case t_COL:
650 7 : return gc_upto(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 127593 : galoisapply(GEN nf, GEN aut, GEN x)
674 : {
675 127593 : pari_sp av = avma;
676 : long lx;
677 : GEN y;
678 :
679 127593 : nf = checknf(nf);
680 127593 : 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 1442 : case t_POL:
687 1442 : aut = algtobasis(nf, aut);
688 1442 : y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
689 1442 : return gc_upto(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 gc_GEN(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 gc_upto(av, y);
702 : }
703 0 : break;
704 :
705 89632 : case t_COL:
706 89632 : aut = algtobasis(nf, aut);
707 89633 : return gc_upto(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 gc_upto(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 13916 : nfgaloismatrixapply(GEN nf, GEN M, GEN x)
722 : {
723 13916 : pari_sp av = avma;
724 : long lx;
725 : GEN y;
726 :
727 13916 : nf = checknf(nf);
728 13916 : 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 gc_upto(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 gc_GEN(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 gc_upto(av, y);
748 : }
749 0 : break;
750 :
751 6027 : case t_COL: return RgM_RgC_mul(M, x);
752 :
753 7154 : case t_MAT: /* ideal */
754 7154 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
755 7154 : if (nbrows(x) != nf_get_degree(nf)) break;
756 7154 : return gc_upto(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 96488 : nfgaloismatrix(GEN nf, GEN s)
765 : {
766 96488 : pari_sp av2, av = avma;
767 : GEN zk, D, M, H, m;
768 : long k, n;
769 :
770 96488 : nf = checknf(nf);
771 96487 : zk = nf_get_zkprimpart(nf); n = lg(zk)-1;
772 96487 : M = cgetg(n+1, t_MAT);
773 96487 : gel(M,1) = col_ei(n, 1); /* s(1) = 1 */
774 96487 : if (n == 1) return M;
775 96487 : av2 = avma;
776 96487 : if (typ(s) != t_COL) s = algtobasis(nf, s);
777 96487 : D = nf_get_zkden(nf);
778 96487 : H = RgV_to_RgM(zk, n);
779 96487 : if (n == 2)
780 : {
781 63433 : GEN t = gel(H,2); /* D * s(w_2) */
782 63433 : t = ZC_Z_add(ZC_Z_mul(s, gel(t,2)), gel(t,1));
783 63433 : gel(M,2) = gc_upto(av2, gdiv(t, D));
784 63434 : return M;
785 : }
786 33054 : m = zk_multable(nf, s);
787 33052 : gel(M,2) = s; /* M[,k] = s(x^(k-1)) */
788 143552 : for (k = 3; k <= n; k++) gel(M,k) = ZM_ZC_mul(m, gel(M,k-1));
789 33053 : M = ZM_mul(M, H);
790 33054 : if (!equali1(D)) M = ZM_Z_divexact(M, D);
791 33054 : return gc_upto(av, M);
792 : }
793 :
794 : static GEN
795 8190 : get_aut(GEN nf, GEN gal, GEN aut, GEN g)
796 : {
797 8190 : return aut ? gel(aut, g[1]): poltobasis(nf, galoispermtopol(gal, g));
798 : }
799 :
800 : static GEN
801 1491 : idealquasifrob(GEN nf, GEN gal, GEN grp, GEN pr, GEN subg, GEN *S, GEN aut)
802 : {
803 1491 : pari_sp av = avma;
804 1491 : long i, n = nf_get_degree(nf), f = pr_get_f(pr);
805 1491 : GEN pi = pr_get_gen(pr);
806 5680 : for (i=1; i<=n; i++)
807 : {
808 5680 : GEN g = gel(grp,i);
809 5680 : if ((!subg && perm_orderu(g) == (ulong)f)
810 5194 : || (subg && perm_relorder(g, subg)==f))
811 : {
812 1956 : *S = get_aut(nf, gal, aut, g);
813 1956 : 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]) = gc_upto(av, vec);
833 : }
834 1379 : return aut;
835 : }
836 :
837 : static void
838 2513 : gal_check_pol(const char *f, GEN x, GEN y)
839 2513 : { if (!RgX_equal_var(x,y)) pari_err_MODULUS(f,x,y); }
840 :
841 : /* true nf */
842 : GEN
843 98 : idealfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN aut)
844 : {
845 98 : pari_sp av = avma;
846 98 : GEN S=NULL, g=NULL; /*-Wall*/
847 : GEN T, p, a, b, modpr;
848 : long f, n, s;
849 98 : f = pr_get_f(pr); n = nf_get_degree(nf);
850 98 : if (f==1) { set_avma(av); return identity_perm(n); }
851 98 : g = idealquasifrob(nf, gal, gal_get_group(gal), pr, NULL, &S, aut);
852 98 : if (f==2) return gc_leaf(av, g);
853 42 : modpr = zk_to_Fq_init(nf,&pr,&T,&p);
854 42 : a = pol_x(nf_get_varn(nf));
855 42 : b = nf_to_Fq(nf, zk_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
856 84 : for (s = 1; s < f-1; s++)
857 : {
858 42 : a = Fq_pow(a, p, T, p);
859 42 : if (ZX_equal(a, b)) break;
860 : }
861 42 : g = perm_powu(g, Fl_inv(s, f));
862 42 : return gc_upto(av, g);
863 : }
864 :
865 : GEN
866 105 : idealfrobenius(GEN nf, GEN gal, GEN pr)
867 : {
868 105 : nf = checknf(nf);
869 105 : checkgal(gal);
870 105 : checkprid(pr);
871 105 : gal_check_pol("idealfrobenius",nf_get_pol(nf),gal_get_pol(gal));
872 105 : if (pr_get_e(pr)>1) pari_err_DOMAIN("idealfrobenius","pr.e", ">", gen_1,pr);
873 98 : 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 gc_upto(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 gc_leaf(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 = galoissubgroups(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 gc_GEN(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 679094 : get_nfpol(GEN x, GEN *nf)
1092 : {
1093 679094 : if (typ(x) == t_POL) { *nf = NULL; return x; }
1094 391101 : *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 gc_GEN(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 gc_GEN(av,y);
1221 : }
1222 :
1223 : static GEN
1224 7350 : partmap_reverse(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1225 : {
1226 7350 : pari_sp av = avma;
1227 7350 : GEN rnf = rnfequation2(a, t), z;
1228 7350 : if (!RgX_equal(b, gel(rnf,1)))
1229 7 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1230 7343 : z = liftpol_shallow(gel(rnf, 2));
1231 7343 : setvarn(z, v);
1232 7343 : if (!isint1(lb)) z = RgX_unscale(z, lb);
1233 7343 : if (!isint1(la)) z = RgX_Rg_div(z, la);
1234 7343 : return gc_GEN(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 8078 : if (degpol(G)) { N = RgX_div(N,G); D = RgX_div(D,G); }
1252 8078 : if (!isint1(lb)) { N = RgX_unscale(N, lb); D = RgX_unscale(D, lb); }
1253 8078 : if (!isint1(la)) D = RgX_Rg_mul(D, la);
1254 8078 : return gc_GEN(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 7231 : nfisincl_from_fact(GEN a, long da, GEN b, GEN la, GEN lb, long vb, GEN y,
1263 : long flag)
1264 : {
1265 7231 : long i, k, l = lg(y), db = degpol(b), d = db / da;
1266 7231 : GEN x = cgetg(l, t_VEC);
1267 7630 : for (i= k = 1; i < l; i++)
1268 : {
1269 7448 : GEN t = gel(y,i);
1270 7448 : if (degpol(t) != d) continue;
1271 7350 : gel(x, k++) = partmap_reverse(a, b, t, la, lb, vb);
1272 7343 : 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 7343 : nfisincl0(GEN fa, GEN fb, long flag)
1299 : {
1300 7343 : pari_sp av = avma;
1301 : long vb;
1302 : GEN a, b, nfa, nfb, x, y, la, lb;
1303 : int newvar;
1304 7343 : if (flag < 0 || flag > 3) pari_err_FLAG("nfisincl");
1305 :
1306 7343 : a = get_nfpol(fa, &nfa);
1307 7343 : b = get_nfpol(fb, &nfb);
1308 7343 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nsisincl"); }
1309 7343 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nsisincl"); }
1310 7343 : 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 gc_GEN(av,x);
1319 : }
1320 7322 : if (flag==0 && !tests_OK(a, nfa, b, nfb, 0)) { set_avma(av); return gen_0; }
1321 :
1322 7203 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1323 7203 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1324 7203 : vb = varn(b); newvar = (varncmp(varn(a),vb) <= 0);
1325 7203 : if (newvar) { b = leafcopy(b); setvarn(b, fetch_var_higher()); }
1326 7203 : y = lift_shallow(gel(nffactor(nfa,b),1));
1327 7203 : if (flag==2)
1328 0 : x = nfisincl_from_fact_frac(a, b, la, lb, vb, y);
1329 : else
1330 7203 : x = nfisincl_from_fact(nfa, degpol(a), b, la, lb, vb, y, flag);
1331 7196 : if (newvar) (void)delete_var();
1332 7196 : return gc_GEN(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 80640 : 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, vecsmallpermute(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 gc_GEN(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 236337 : get_roots(GEN x, long r1, long prec)
1492 : {
1493 : long i, ru;
1494 : GEN z;
1495 236337 : if (typ(x) != t_POL)
1496 : {
1497 0 : z = leafcopy(x);
1498 0 : ru = (lg(z)-1 + r1) >> 1;
1499 : }
1500 : else
1501 : {
1502 236337 : long n = degpol(x);
1503 236337 : z = (r1 == n)? ZX_realroots_irred(x, prec): QX_complex_roots(x,prec);
1504 236343 : ru = (n+r1)>>1;
1505 : }
1506 512821 : for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
1507 236343 : 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 394244 : quicktrace(GEN x, GEN sym)
1519 : {
1520 394244 : GEN p1 = gen_0;
1521 : long i;
1522 :
1523 394244 : if (typ(x) != t_POL) return gmul(x, gel(sym,1));
1524 394244 : if (signe(x))
1525 : {
1526 394245 : sym--;
1527 2896146 : for (i=lg(x)-1; i>1; i--)
1528 2501930 : p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
1529 : }
1530 394215 : return p1;
1531 : }
1532 :
1533 : static GEN
1534 85466 : get_Tr(GEN mul, GEN x, GEN basden)
1535 : {
1536 85466 : GEN t, bas = gel(basden,1), den = gel(basden,2);
1537 85466 : long i, j, n = lg(bas)-1;
1538 85466 : GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
1539 :
1540 85466 : gel(TW,1) = utoipos(n);
1541 269938 : for (i=2; i<=n; i++)
1542 : {
1543 184474 : t = quicktrace(gel(bas,i), sym);
1544 184471 : if (den && gel(den,i)) t = diviiexact(t,gel(den,i));
1545 184471 : gel(TW,i) = t; /* tr(w[i]) */
1546 : }
1547 85464 : gel(T,1) = TW;
1548 269940 : for (i=2; i<=n; i++)
1549 : {
1550 184473 : gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
1551 742547 : for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
1552 558071 : gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
1553 : }
1554 85467 : return T;
1555 : }
1556 :
1557 : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
1558 : static GEN
1559 206655 : get_bas_den(GEN bas)
1560 : {
1561 206655 : GEN b,d,den, dbas = leafcopy(bas);
1562 206655 : long i, l = lg(bas);
1563 206655 : int power = 1;
1564 206655 : den = cgetg(l,t_VEC);
1565 959214 : for (i=1; i<l; i++)
1566 : {
1567 752562 : b = Q_remove_denom(gel(bas,i), &d);
1568 752554 : gel(dbas,i) = b;
1569 752554 : gel(den,i) = d; if (d) power = 0;
1570 : }
1571 206652 : if (power) den = NULL; /* power basis */
1572 206652 : return mkvec2(dbas, den);
1573 : }
1574 :
1575 : /* return multiplication table for S->basis */
1576 : static GEN
1577 85465 : nf_multable(nfmaxord_t *S, GEN invbas)
1578 : {
1579 85465 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1580 85465 : long i,j, n = degpol(T);
1581 85465 : GEN mul = cgetg(n*n+1,t_MAT);
1582 :
1583 : /* i = 1 split for efficiency, assume w[1] = 1 */
1584 355405 : for (j=1; j<=n; j++)
1585 269940 : gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
1586 269939 : for (i=2; i<=n; i++)
1587 742552 : for (j=i; j<=n; j++)
1588 : {
1589 558078 : pari_sp av = avma;
1590 558078 : GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1591 558081 : z = ZM_ZX_mul(invbas, z); /* integral column */
1592 558053 : if (den)
1593 : {
1594 353333 : GEN d = mul_denom(gel(den,i), gel(den,j));
1595 353331 : if (d) z = ZC_Z_divexact(z, d);
1596 : }
1597 558048 : gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gc_upto(av,z);
1598 : }
1599 85464 : return mul;
1600 : }
1601 :
1602 : /* as get_Tr, mul_table not precomputed */
1603 : static GEN
1604 28897 : make_Tr(nfmaxord_t *S)
1605 : {
1606 28897 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1607 28897 : long i,j, n = degpol(T);
1608 28897 : 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 28896 : c = cgetg(n+1,t_COL); gel(M,1) = c;
1612 28896 : gel(c, 1) = utoipos(n);
1613 86303 : for (j=2; j<=n; j++)
1614 : {
1615 57409 : pari_sp av = avma;
1616 57409 : t = quicktrace(gel(w,j), sym);
1617 57406 : if (den)
1618 : {
1619 36721 : d = gel(den,j);
1620 36721 : if (d) t = diviiexact(t, d);
1621 : }
1622 57405 : gel(c,j) = gc_INT(av, t);
1623 : }
1624 86296 : for (i=2; i<=n; i++)
1625 : {
1626 57406 : c = cgetg(n+1,t_COL); gel(M,i) = c;
1627 209772 : for (j=1; j<i ; j++) gel(c,j) = gcoeff(M,i,j);
1628 209769 : for ( ; j<=n; j++)
1629 : {
1630 152367 : pari_sp av = avma;
1631 152367 : t = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1632 152366 : t = quicktrace(t, sym);
1633 152354 : if (den)
1634 : {
1635 119235 : d = mul_denom(gel(den,i),gel(den,j));
1636 119237 : if (d) t = diviiexact(t, d);
1637 : }
1638 152356 : gel(c,j) = gc_INT(av, t); /* Tr (W[i]W[j]) */
1639 : }
1640 : }
1641 28890 : return M;
1642 : }
1643 :
1644 : /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
1645 : static void
1646 269601 : make_M(nffp_t *F, int trunc)
1647 : {
1648 269601 : GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
1649 : GEN m, d, M;
1650 269601 : long i, j, l = lg(ro), n = lg(bas);
1651 269601 : M = cgetg(n,t_MAT);
1652 269599 : gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
1653 941215 : for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
1654 872348 : for (i=1; i<l; i++)
1655 : {
1656 602779 : GEN r = gel(ro,i), ri;
1657 602779 : ri = (gexpo(r) > 1)? ginv(r): NULL;
1658 2926302 : for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
1659 : }
1660 269569 : if (den)
1661 534443 : for (j=2; j<n; j++)
1662 : {
1663 414873 : d = gel(den,j); if (!d) continue;
1664 350747 : m = gel(M,j);
1665 1754389 : for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
1666 : }
1667 :
1668 269570 : if (trunc && gprecision(M) > F->prec)
1669 : {
1670 18361 : M = gprec_w(M, F->prec);
1671 18361 : F->ro = gprec_w(ro,F->prec);
1672 : }
1673 269570 : F->M = M;
1674 269570 : }
1675 :
1676 : /* return G real such that G~ * G = T_2 */
1677 : static void
1678 269601 : make_G(nffp_t *F)
1679 : {
1680 269601 : GEN G, M = F->M;
1681 269601 : long i, j, k, r1 = F->r1, l = lg(M);
1682 :
1683 269601 : if (r1 == l-1) { F->G = M; return; }
1684 221101 : G = cgetg(l, t_MAT);
1685 1024345 : for (j = 1; j < l; j++)
1686 : {
1687 803299 : GEN g, m = gel(M,j);
1688 803299 : gel(G,j) = g = cgetg(l, t_COL);
1689 1335604 : for (k = i = 1; i <= r1; i++) gel(g,k++) = gel(m,i);
1690 2657110 : for ( ; k < l; i++)
1691 : {
1692 1853866 : GEN r = gel(m,i);
1693 1853866 : if (typ(r) == t_COMPLEX)
1694 : {
1695 1515429 : GEN a = gel(r,1), b = gel(r,2);
1696 1515429 : gel(g,k++) = mpadd(a, b);
1697 1515401 : gel(g,k++) = mpsub(a, b);
1698 : }
1699 : else
1700 : {
1701 338437 : gel(g,k++) = r;
1702 338437 : gel(g,k++) = r;
1703 : }
1704 : }
1705 : }
1706 221046 : F->G = G;
1707 : }
1708 :
1709 : static long
1710 302230 : prec_fix(long prec)
1711 : {
1712 : #ifndef LONG_IS_64BIT
1713 : /* make sure that default accuracy is the same on 32/64bit */
1714 43246 : if (odd(prec2lg(prec))) prec+=EXTRAPRECWORD;
1715 : #endif
1716 302230 : return prec;
1717 : }
1718 : static void
1719 269600 : make_M_G(nffp_t *F, int trunc)
1720 : {
1721 : long n, eBD, prec;
1722 269600 : if (F->extraprec < 0)
1723 : { /* not initialized yet; compute roots so that absolute accuracy
1724 : * of M & G >= prec */
1725 : double er;
1726 269579 : n = degpol(F->T);
1727 269579 : eBD = 1 + gexpo(gel(F->basden,1));
1728 269579 : er = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->T);
1729 269576 : if (er < 0) er = 0;
1730 269576 : F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
1731 : }
1732 269596 : prec = prec_fix(F->prec + F->extraprec);
1733 269596 : if (!F->ro || gprecision(gel(F->ro,1)) < prec)
1734 236337 : F->ro = get_roots(F->T, F->r1, prec);
1735 :
1736 269602 : make_M(F, trunc);
1737 269601 : make_G(F);
1738 269593 : }
1739 :
1740 : static void
1741 177822 : nffp_init(nffp_t *F, nfmaxord_t *S, long prec)
1742 : {
1743 177822 : F->T = S->T;
1744 177822 : F->r1 = S->r1;
1745 177822 : F->basden = S->basden;
1746 177822 : F->ro = NULL;
1747 177822 : F->extraprec = -1;
1748 177822 : F->prec = prec;
1749 177822 : }
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 2597 : get_nfindex(GEN bas)
1756 : {
1757 2597 : pari_sp av = avma;
1758 2597 : long n = lg(bas)-1, i;
1759 : GEN D, d, mat;
1760 :
1761 : /* assume bas[1] = 1 */
1762 2597 : D = gel(bas,1);
1763 2597 : if (! is_pm1(simplify_shallow(D))) pari_err_TYPE("get_nfindex", D);
1764 2597 : D = gen_1;
1765 12418 : for (i = 2; i <= n; i++)
1766 : { /* after nfbasis, basis is upper triangular! */
1767 9828 : GEN B = gel(bas,i), lc;
1768 9828 : if (degpol(B) != i-1) break;
1769 9821 : lc = gel(B, i+1);
1770 9821 : switch (typ(lc))
1771 : {
1772 4013 : case t_INT: continue;
1773 5808 : 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 2597 : 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 2597 : return gc_INT(av, D);
1787 : }
1788 : /* make sure all components of S are initialized */
1789 : static void
1790 170927 : nfmaxord_complete(nfmaxord_t *S)
1791 : {
1792 170927 : if (!S->dT) S->dT = ZX_disc(S->T);
1793 170927 : if (!S->index)
1794 : {
1795 2597 : if (S->dK) /* fast */
1796 0 : S->index = sqrti( diviiexact(S->dT, S->dK) );
1797 : else
1798 2597 : S->index = get_nfindex(S->basis);
1799 : }
1800 170927 : if (!S->dK) S->dK = diviiexact(S->dT, sqri(S->index));
1801 170927 : if (S->r1 < 0) S->r1 = ZX_sturm_irred(S->T);
1802 170927 : if (!S->basden) S->basden = get_bas_den(S->basis);
1803 170928 : }
1804 :
1805 : GEN
1806 85463 : nfmaxord_to_nf(nfmaxord_t *S, GEN ro, long prec)
1807 : {
1808 85463 : GEN nf = cgetg(10,t_VEC);
1809 85465 : GEN T = S->T, Tr, D, w, A, dA, MDI, mat = cgetg(9,t_VEC);
1810 85466 : long n = degpol(T);
1811 : nffp_t F;
1812 85466 : nfmaxord_complete(S);
1813 85466 : nffp_init(&F,S,prec);
1814 85466 : F.ro = ro;
1815 85466 : make_M_G(&F, 0);
1816 :
1817 85466 : gel(nf,1) = S->T;
1818 85466 : gel(nf,2) = mkvec2s(S->r1, (n - S->r1)>>1);
1819 85465 : gel(nf,3) = S->dK;
1820 85465 : gel(nf,4) = S->index;
1821 85465 : gel(nf,5) = mat;
1822 85465 : if (gprecision(gel(F.ro,1)) > prec) F.ro = gprec_wtrunc(F.ro, prec);
1823 85464 : gel(nf,6) = F.ro;
1824 85464 : w = S->basis;
1825 85464 : if (!is_pm1(S->index)) w = Q_remove_denom(w, NULL);
1826 85464 : gel(nf,7) = w;
1827 85464 : gel(nf,8) = ZM_inv(RgV_to_RgM(w,n), NULL);
1828 85465 : gel(nf,9) = nf_multable(S, nf_get_invzk(nf));
1829 85466 : gel(mat,1) = F.M;
1830 85466 : gel(mat,2) = F.G;
1831 :
1832 85466 : Tr = get_Tr(gel(nf,9), T, S->basden);
1833 85464 : gel(mat,6) = A = ZM_inv(Tr, &dA); /* dA T^-1, primitive */
1834 85465 : 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 85465 : MDI = idealtwoelt(nf, A);
1838 85466 : gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
1839 85465 : gel(mat,7) = MDI;
1840 85465 : if (is_pm1(S->index))
1841 : { /* principal ideal (T'), whose norm is |dK| */
1842 49660 : D = zk_scalar_or_multable(nf, ZX_deriv(T));
1843 49661 : if (typ(D) == t_MAT) D = ZM_hnfmod(D, absi_shallow(S->dK));
1844 : }
1845 : else
1846 : {
1847 35805 : GEN c = diviiexact(dA, gcoeff(A,1,1));
1848 35805 : D = idealHNF_inv_Z(nf, A); /* (A\cap Z) / A */
1849 35805 : if (!is_pm1(c)) D = ZM_Z_mul(D, c);
1850 : }
1851 85466 : gel(mat,3) = RM_round_maxrank(F.G);
1852 85462 : gel(mat,4) = Tr;
1853 85462 : gel(mat,5) = D;
1854 85462 : w = S->dKP; if (!w) { w = gel(absZ_factor(S->dK), 1); settyp(w, t_VEC); }
1855 85462 : gel(mat,8) = w; return nf;
1856 : }
1857 :
1858 : static GEN
1859 3374 : primes_certify(GEN dK, GEN dKP)
1860 : {
1861 3374 : long i, l = lg(dKP);
1862 3374 : GEN v, w, D = dK;
1863 3374 : v = vectrunc_init(l);
1864 3374 : w = vectrunc_init(l);
1865 10479 : for (i = 1; i < l; i++)
1866 : {
1867 7105 : GEN p = gel(dKP,i);
1868 7105 : vectrunc_append(isprime(p)? w: v, p);
1869 7105 : (void)Z_pvalrem(D, p, &D);
1870 : }
1871 3374 : if (!is_pm1(D))
1872 : {
1873 0 : D = absi_shallow(D);
1874 0 : vectrunc_append(isprime(D)? w: v, D);
1875 : }
1876 3374 : return mkvec2(v,w);
1877 : }
1878 : GEN
1879 3374 : nfcertify(GEN nf)
1880 : {
1881 3374 : pari_sp av = avma;
1882 : GEN vw;
1883 3374 : nf = checknf(nf);
1884 3374 : vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
1885 3374 : return gc_GEN(av, gel(vw,1));
1886 : }
1887 :
1888 : /* set *pro to roots of S->T */
1889 : static GEN
1890 73981 : get_red_G(nfmaxord_t *S, GEN *pro)
1891 : {
1892 73981 : pari_sp av = avma;
1893 73981 : GEN G, u, u0 = NULL;
1894 73981 : long prec, n = degpol(S->T);
1895 : nffp_t F;
1896 :
1897 73981 : prec = nbits2prec(n+32);
1898 73981 : nffp_init(&F, S, prec);
1899 : for (;;)
1900 : {
1901 73981 : F.prec = prec; make_M_G(&F, 0); G = F.G;
1902 73983 : if (u0) G = RgM_mul(G, u0);
1903 73983 : 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 73983 : if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST)))
1907 : {
1908 73983 : if (lg(u)-1 == n) break;
1909 : /* singular ==> loss of accuracy */
1910 0 : if (u0) u0 = gc_upto(av, RgM_mul(u0,u));
1911 0 : else u0 = gc_GEN(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 73983 : if (u0) u = RgM_mul(u0,u);
1918 73983 : *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 103565 : set_LLL_basis(nfmaxord_t *S, GEN *pro, long flag, double DELTA)
1925 : {
1926 103565 : GEN B = S->basis;
1927 103565 : long N = degpol(S->T);
1928 103565 : if (S->r1 < 0)
1929 : {
1930 18200 : S->r1 = ZX_sturm_irred(S->T);
1931 18200 : if (odd(N - S->r1)) pari_err_IRREDPOL("set_LLL_basis", S->T);
1932 : }
1933 103558 : if (!S->basden) S->basden = get_bas_den(B);
1934 103558 : *pro = NULL; if (flag & nf_NOLLL) return;
1935 102879 : if (S->r1 == N) {
1936 28897 : pari_sp av = avma;
1937 28897 : GEN u = ZM_lll(make_Tr(S), DELTA, LLL_GRAM|LLL_KEEP_FIRST|LLL_IM);
1938 28897 : B = gc_upto(av, RgV_RgM_mul(B, u));
1939 : }
1940 : else
1941 73982 : B = RgV_RgM_mul(B, get_red_G(S, pro));
1942 102881 : S->basis = B;
1943 102881 : S->basden = get_bas_den(B);
1944 : }
1945 :
1946 : /* = 1 iff |a| > |b| or equality and a > 0 */
1947 : static int
1948 128702 : cmpii_polred(GEN a, GEN b)
1949 : {
1950 128702 : int fl = abscmpii(a, b);
1951 : long sa, sb;
1952 128702 : if (fl) return fl;
1953 106980 : sa = signe(a);
1954 106980 : sb = signe(b);
1955 106980 : if (sa == sb) return 0;
1956 936 : return sa == 1? 1: -1;
1957 : }
1958 : static int
1959 30428 : ZX_cmp(GEN x, GEN y)
1960 30428 : { 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 52159 : ZX_is_better(GEN y, GEN x, GEN *dx)
1965 : {
1966 : pari_sp av;
1967 : int cmp;
1968 : GEN d;
1969 52159 : if (!*dx) *dx = ZX_disc(x);
1970 52159 : av = avma; d = ZX_disc(y);
1971 52159 : cmp = abscmpii(d, *dx);
1972 52159 : if (cmp < 0) { *dx = d; return 1; }
1973 47648 : 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 87344 : nf_input_type(GEN x)
2014 : {
2015 87344 : GEN T, V, DKP = NULL;
2016 : long i, d, v;
2017 87344 : switch(typ(x))
2018 : {
2019 80694 : case t_POL: return t_POL;
2020 6650 : case t_VEC:
2021 6650 : switch(lg(x))
2022 : {
2023 2849 : case 4: DKP = gel(x,3);
2024 6622 : case 3: break;
2025 28 : default: return t_VEC; /* nf or incorrect */
2026 : }
2027 6622 : T = gel(x,1); V = gel(x,2);
2028 6622 : if (typ(T) != t_POL) return -1;
2029 6622 : switch(typ(V))
2030 : {
2031 210 : case t_INT: case t_MAT: return t_POL;
2032 6412 : case t_VEC: case t_COL:
2033 6412 : if (RgV_is_ZV(V)) return t_POL;
2034 3514 : break;
2035 0 : default: return -1;
2036 : }
2037 3514 : d = degpol(T); v = varn(T);
2038 3514 : if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
2039 21504 : for (i = 1; i <= d; i++)
2040 : { /* check integer basis */
2041 18011 : GEN c = gel(V,i);
2042 18011 : switch(typ(c))
2043 : {
2044 35 : case t_INT: break;
2045 17976 : 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 3493 : if (DKP && (typ(DKP) != t_VEC || !RgV_is_ZV(DKP))) return -1;
2051 3493 : return 0;
2052 : }
2053 0 : return t_VEC; /* nf or incorrect */
2054 : }
2055 :
2056 : /* cater for obsolete nf_PARTIALFACT flag */
2057 : static void
2058 3941 : nfinit_basic_partial(nfmaxord_t *S, GEN T)
2059 : {
2060 3941 : if (typ(T) == t_POL) { nfmaxord(S, mkvec2(T,utoipos(500000)), 0); }
2061 14 : else nfinit_basic(S, T);
2062 3941 : }
2063 : static void
2064 14343 : nfinit_basic_flag(nfmaxord_t *S, GEN x, long flag)
2065 : {
2066 14343 : if (flag & nf_PARTIALFACT)
2067 35 : nfinit_basic_partial(S, x);
2068 : else
2069 14308 : nfinit_basic(S, x);
2070 14336 : }
2071 :
2072 : /* true nf */
2073 : static GEN
2074 91785 : nf_basden(GEN nf)
2075 : {
2076 91785 : GEN zkD = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
2077 91785 : D = equali1(D)? NULL: const_vec(lg(zkD)-1, D);
2078 91785 : return mkvec2(zkD, D);
2079 : }
2080 : void
2081 87344 : nfinit_basic(nfmaxord_t *S, GEN T)
2082 : {
2083 87344 : switch (nf_input_type(T))
2084 : {
2085 83802 : 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 3493 : case 0: /* monic integral polynomial + integer basis (+ ramified primes)*/
2099 3493 : S->T = S->T0 = gel(T,1);
2100 3493 : S->basis = gel(T,2);
2101 3493 : S->basden = NULL;
2102 3493 : S->index = NULL;
2103 3493 : S->dK = NULL;
2104 3493 : S->dKP = NULL;
2105 3493 : if (lg(T) == 4)
2106 : {
2107 2849 : GEN v = gel(T,3); if (typ(v) == t_COL) v = shallowtrans(v);
2108 2849 : S->dKP = v;
2109 : }
2110 3493 : S->dT = NULL;
2111 3493 : S->r1 = -1; break;
2112 20 : default: /* -1 */
2113 20 : pari_err_TYPE("nfinit_basic", T);
2114 : }
2115 3521 : S->dTP = S->dTE = S->dKE = NULL;
2116 3521 : S->unscale = gen_1;
2117 : }
2118 :
2119 : GEN
2120 85466 : nfinit_complete(nfmaxord_t *S, long flag, long prec)
2121 : {
2122 85466 : GEN nf, unscale = S->unscale, rev = NULL;
2123 :
2124 85466 : if (!ZX_is_irred(S->T)) pari_err_IRREDPOL("nfinit",S->T);
2125 85462 : 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 85462 : if (!(flag & nf_RED) && !isint1(unscale))
2131 : { /* implies lc(x0) = 1 and L := 1/unscale is integral */
2132 4585 : long d = degpol(S->T0);
2133 4585 : GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
2134 4585 : GEN f = powiu(L, (d*(d-1)) >> 1);
2135 4585 : S->T = S->T0; /* restore original user-supplied x0, unscale data */
2136 4585 : S->unscale = gen_1;
2137 4585 : S->dT = gmul(S->dT, sqri(f));
2138 4585 : S->basis = RgXV_unscale(S->basis, unscale);
2139 4585 : S->index = gmul(S->index, f);
2140 : }
2141 85462 : nfmaxord_complete(S); /* more expensive after set_LLL_basis */
2142 85463 : 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 85183 : GEN ro; set_LLL_basis(S, &ro, flag, 0.99);
2155 85182 : nf = nfmaxord_to_nf(S, ro, prec);
2156 : }
2157 85462 : 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 85462 : 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 11721 : nfinit0(GEN x, long flag,long prec)
2175 : {
2176 11721 : const pari_sp av = avma;
2177 : nfmaxord_t S;
2178 11721 : if (flag < 0 || flag > 7) pari_err_FLAG("nfinit");
2179 11721 : if (checkrnf_i(x)) return rnf_build_nfabs(x, prec);
2180 11700 : nfinit_basic(&S, x);
2181 11679 : return gc_GEN(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 7318 : nfinit(GEN x, long prec) { return nfinit0(x, 0, prec); }
2189 :
2190 : /* assume x a bnr/bnf/nf */
2191 : long
2192 1273646 : nf_get_prec(GEN x)
2193 : {
2194 1273646 : GEN nf = checknf(x), ro = nf_get_roots(nf);
2195 1273643 : return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
2196 : }
2197 :
2198 : /* true nf */
2199 : GEN
2200 91757 : nfnewprec_shallow(GEN nf, long prec)
2201 : {
2202 91757 : GEN m, NF = leafcopy(nf);
2203 : nffp_t F;
2204 :
2205 91757 : F.T = nf_get_pol(nf);
2206 91757 : F.ro = NULL;
2207 91757 : F.r1 = nf_get_r1(nf);
2208 91757 : F.basden = nf_basden(nf);
2209 91757 : F.extraprec = -1;
2210 91757 : F.prec = prec; make_M_G(&F, 0);
2211 91756 : gel(NF,5) = m = leafcopy(gel(NF,5));
2212 91757 : gel(m,1) = F.M;
2213 91757 : gel(m,2) = F.G;
2214 91757 : gel(NF,6) = F.ro; return NF;
2215 : }
2216 :
2217 : /* true rnf */
2218 : GEN
2219 56 : rnfnewprec_shallow(GEN rnf, long prec)
2220 : {
2221 56 : GEN RNF = leafcopy(rnf);
2222 56 : gel(RNF,10) = nfnewprec_shallow(gel(RNF,10), prec);
2223 56 : if (obj_check(RNF, rnf_NFABS)) rnfcomplete(RNF);
2224 56 : return RNF;
2225 : }
2226 : GEN
2227 0 : rnfnewprec(GEN rnf, long prec)
2228 : {
2229 0 : pari_sp av = avma;
2230 0 : checkrnf(rnf);
2231 0 : return gc_GEN(av, rnfnewprec_shallow(rnf, prec));
2232 : }
2233 :
2234 : GEN
2235 161 : nfnewprec(GEN nf, long prec)
2236 : {
2237 161 : pari_sp av = avma;
2238 : GEN z;
2239 161 : switch(nftyp(nf))
2240 : {
2241 56 : default: pari_err_TYPE("nfnewprec", nf);
2242 0 : case typ_RNF: z = rnfnewprec_shallow(nf,prec); break;
2243 14 : case typ_BNF: z = bnfnewprec_shallow(nf,prec); break;
2244 84 : case typ_NF: z = nfnewprec_shallow(nf, prec); break;
2245 7 : case typ_BNR: return bnrnewprec(nf,prec);
2246 : }
2247 98 : return gc_GEN(av, z);
2248 : }
2249 :
2250 : /********************************************************************/
2251 : /** **/
2252 : /** POLRED **/
2253 : /** **/
2254 : /********************************************************************/
2255 : GEN
2256 0 : embednorm_T2(GEN x, long r1)
2257 : {
2258 0 : pari_sp av = avma;
2259 0 : GEN p = RgV_sumpart(x, r1);
2260 0 : GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
2261 0 : if (q != gen_0) p = gadd(p, gmul2n(q,1));
2262 0 : return avma == av? gcopy(p): gc_upto(av, p);
2263 : }
2264 :
2265 : /* simplified version of gnorm for scalar, noncomplex inputs, without GC */
2266 : static GEN
2267 167488 : real_norm(GEN x)
2268 : {
2269 167488 : switch(typ(x))
2270 : {
2271 0 : case t_INT: return sqri(x);
2272 167488 : case t_REAL: return sqrr(x);
2273 0 : case t_FRAC: return sqrfrac(x);
2274 : }
2275 0 : pari_err_TYPE("real_norm", x);
2276 : return NULL;/*LCOV_EXCL_LINE*/
2277 : }
2278 : /* simplified version of gnorm, without GC */
2279 : static GEN
2280 33685460 : complex_norm(GEN x)
2281 : {
2282 33685460 : return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
2283 : }
2284 : /* return T2(x), argument r1 needed in case x has components whose type
2285 : * is unexpected, e.g. all of them t_INT for embed(gen_1) */
2286 : GEN
2287 72148 : embed_T2(GEN x, long r1)
2288 : {
2289 72148 : pari_sp av = avma;
2290 72148 : long i, l = lg(x);
2291 72148 : GEN c, s = NULL, t = NULL;
2292 72148 : if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
2293 239636 : for (i = 1; i <= r1; i++)
2294 : {
2295 167488 : c = real_norm(gel(x,i));
2296 167490 : s = s? gadd(s, c): c;
2297 : }
2298 206154 : for (; i < l; i++)
2299 : {
2300 134008 : c = complex_norm(gel(x,i));
2301 134008 : t = t? gadd(t, c): c;
2302 : }
2303 72146 : if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
2304 72146 : return gc_upto(av, s);
2305 : }
2306 : /* return N(x) */
2307 : GEN
2308 17286837 : embed_norm(GEN x, long r1)
2309 : {
2310 17286837 : pari_sp av = avma;
2311 17286837 : long i, l = lg(x);
2312 17286837 : GEN c, s = NULL, t = NULL;
2313 17286837 : if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
2314 44729793 : for (i = 1; i <= r1; i++)
2315 : {
2316 27570222 : c = gel(x,i);
2317 27570222 : s = s? gmul(s, c): c;
2318 : }
2319 50710696 : for (; i < l; i++)
2320 : {
2321 33551440 : c = complex_norm(gel(x,i));
2322 33550541 : t = t? gmul(t, c): c;
2323 : }
2324 17159256 : if (t) s = s? gmul(s,t): t;
2325 17159448 : return gc_upto(av, s);
2326 : }
2327 :
2328 : typedef struct {
2329 : long r1, v, prec;
2330 : GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
2331 : GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
2332 : GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
2333 : GEN bound; /* T2 norm of the polynomial defining nf */
2334 : long expo_best_disc; /* expo(disc(x)), best generator so far */
2335 : } CG_data;
2336 :
2337 : /* characteristic pol of x (given by embeddings) */
2338 : static GEN
2339 264603 : get_pol(CG_data *d, GEN x)
2340 : {
2341 : long e;
2342 264603 : GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
2343 264603 : return (e > -5)? NULL: g;
2344 : }
2345 :
2346 : /* characteristic pol of x (given as vector on (w_i)) */
2347 : static GEN
2348 100461 : get_polchar(CG_data *d, GEN x)
2349 100461 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
2350 :
2351 : /* Choose a canonical polynomial in the pair { Pmin_a, Pmin_{-a} }, i.e.
2352 : * { z(X), (-1)^(deg z) z(-Z) } and keeping the smallest wrt cmpii_polred
2353 : * Either leave z alone (return 1) or set z <- (-1)^n z(-X). In place. */
2354 : int
2355 89881 : ZX_canon_neg(GEN z)
2356 : {
2357 : long i, s;
2358 132114 : for (i = lg(z)-2; i >= 2; i -= 2)
2359 : { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
2360 127970 : s = signe(gel(z,i));
2361 127970 : if (!s) continue;
2362 : /* non trivial */
2363 85737 : if (s < 0) break; /* z(X) < (-1)^n z(-X) */
2364 :
2365 210635 : for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
2366 36425 : return 1;
2367 : }
2368 53456 : return 0;
2369 : }
2370 : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
2371 : * Return NULL on failure: discriminant too large or non primitive */
2372 : static GEN
2373 136699 : try_polmin(CG_data *d, nfmaxord_t *S, GEN v, long flag, GEN *ai)
2374 : {
2375 136699 : const long best = flag & nf_ABSOLUTE;
2376 : long ed;
2377 136699 : pari_sp av = avma;
2378 : GEN g;
2379 136699 : if (best)
2380 : {
2381 136132 : ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
2382 136131 : set_avma(av); if (d->expo_best_disc < ed) return NULL;
2383 : }
2384 : else
2385 567 : ed = 0;
2386 82535 : g = get_pol(d, v);
2387 : /* accuracy too low, compute algebraically */
2388 82536 : if (!g) { set_avma(av); g = ZXQ_charpoly(*ai, S->T, varn(S->T)); }
2389 82536 : g = ZX_radical(g);
2390 82536 : if (best && degpol(g) != degpol(S->T)) return gc_NULL(av);
2391 30606 : g = gc_GEN(av, g);
2392 30606 : d->expo_best_disc = ed;
2393 30606 : if (flag & nf_ORIG)
2394 : {
2395 27239 : if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
2396 27239 : if (!isint1(S->unscale)) *ai = RgX_unscale(*ai, S->unscale);
2397 : }
2398 : else
2399 3367 : (void)ZX_canon_neg(g);
2400 30606 : if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
2401 30606 : return g;
2402 : }
2403 :
2404 : /* does x generate the correct field ? */
2405 : static GEN
2406 100461 : chk_gen(void *data, GEN x)
2407 : {
2408 100461 : pari_sp av = avma, av1;
2409 100461 : GEN h, g = get_polchar((CG_data*)data,x);
2410 100461 : if (!g) pari_err_PREC("chk_gen");
2411 100461 : av1 = avma;
2412 100461 : h = ZX_gcd(g, ZX_deriv(g));
2413 100461 : if (degpol(h)) return gc_NULL(av);
2414 56092 : if (DEBUGLEVEL>3) err_printf(" generator: %Ps\n",g);
2415 56092 : set_avma(av1); return gc_upto(av, g);
2416 : }
2417 :
2418 : static long
2419 32634 : chk_gen_prec(long N, long bit)
2420 32634 : { return prec_fix(nbits2prec(10 + (long)log2((double)N) + bit)); }
2421 :
2422 : /* v = [P,A] two vectors (of ZX and ZV resp.) of same length; remove duplicate
2423 : * polynomials in P, updating A, in place. Among elements having the same
2424 : * characteristic pol, choose the smallest according to ZV_abscmp */
2425 : static void
2426 14000 : remove_duplicates(GEN v)
2427 : {
2428 14000 : GEN x, a, P = gel(v,1), A = gel(v,2);
2429 14000 : long k, i, l = lg(P);
2430 14000 : pari_sp av = avma;
2431 :
2432 14000 : if (l < 2) return;
2433 14000 : (void)sort_factor_pol(mkvec2(P, A), cmpii);
2434 14000 : x = gel(P,1); a = gel(A,1);
2435 54103 : for (k=1,i=2; i<l; i++)
2436 40103 : if (ZX_equal(gel(P,i), x))
2437 : {
2438 19747 : if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
2439 : }
2440 : else
2441 : {
2442 20356 : gel(A,k) = a;
2443 20356 : gel(P,k) = x;
2444 20356 : k++;
2445 20356 : x = gel(P,i); a = gel(A,i);
2446 : }
2447 14000 : l = k+1;
2448 14000 : gel(A,k) = a; setlg(A,l);
2449 14000 : gel(P,k) = x; setlg(P,l); set_avma(av);
2450 : }
2451 :
2452 : static void
2453 18382 : polred_init(nfmaxord_t *S, nffp_t *F, CG_data *d)
2454 : {
2455 18382 : long e, prec, n = degpol(S->T);
2456 : double log2rho;
2457 : GEN ro;
2458 18382 : set_LLL_basis(S, &ro, 0, 0.9999);
2459 : /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
2460 18375 : log2rho = ro ? (double)gexpo(ro): fujiwara_bound(S->T);
2461 18375 : e = n * (long)(log2rho + log2((double)n)) + 1;
2462 18375 : if (e < 0) e = 0; /* can occur if n = 1 */
2463 18375 : prec = chk_gen_prec(n, e);
2464 18375 : nffp_init(F,S,prec);
2465 18375 : F->ro = ro;
2466 18375 : make_M_G(F, 1);
2467 :
2468 18375 : d->v = varn(S->T);
2469 18375 : d->expo_best_disc = -1;
2470 18375 : d->ZKembed = NULL;
2471 18375 : d->M = NULL;
2472 18375 : d->u = NULL;
2473 18375 : d->r1= S->r1;
2474 18375 : }
2475 : static GEN
2476 14273 : findmindisc(GEN y)
2477 : {
2478 14273 : GEN x = gel(y,1), dx = NULL;
2479 14273 : long i, l = lg(y);
2480 35980 : for (i = 2; i < l; i++)
2481 : {
2482 21707 : GEN yi = gel(y,i);
2483 21707 : if (ZX_is_better(yi,x,&dx)) x = yi;
2484 : }
2485 14273 : return x;
2486 : }
2487 : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
2488 : * [ the best wrt discriminant ordering ], but keep all imprimitive
2489 : * polynomials */
2490 : static void
2491 4116 : filter(GEN y, GEN b, long n)
2492 : {
2493 : GEN x, a, dx;
2494 4116 : long i, k = 1, l = lg(y);
2495 4116 : a = x = dx = NULL;
2496 34771 : for (i = 1; i < l; i++)
2497 : {
2498 30655 : GEN yi = gel(y,i), ai = gel(b,i);
2499 30655 : if (degpol(yi) == n)
2500 : {
2501 30501 : if (!dx) dx = ZX_disc(yi); else if (!ZX_is_better(yi,x,&dx)) continue;
2502 5696 : x = yi; a = ai; continue;
2503 : }
2504 154 : gel(y,k) = yi;
2505 154 : gel(b,k) = ai; k++;
2506 : }
2507 4116 : if (dx)
2508 : {
2509 4019 : gel(y,k) = x;
2510 4019 : gel(b,k) = a; k++;
2511 : }
2512 4116 : setlg(y, k);
2513 4116 : setlg(b, k);
2514 4116 : }
2515 :
2516 : static GEN
2517 4151 : polred_aux(nfmaxord_t *S, GEN *pro, long flag)
2518 : { /* only keep polynomials of max degree and best discriminant */
2519 4151 : const long best = flag & nf_ABSOLUTE;
2520 4151 : const long orig = flag & nf_ORIG;
2521 4151 : GEN M, b, y, x = S->T;
2522 4151 : long maxi, i, j, k, v = varn(x), n = lg(S->basis)-1;
2523 : nffp_t F;
2524 : CG_data d;
2525 :
2526 4151 : if (n == 1)
2527 : {
2528 28 : if (!best)
2529 : {
2530 14 : GEN X = pol_x(v);
2531 14 : return orig? mkmat2(mkcol(X),mkcol(gen_1)): mkvec(X);
2532 : }
2533 : else
2534 14 : return orig? trivial_fact(): cgetg(1,t_VEC);
2535 : }
2536 :
2537 4123 : polred_init(S, &F, &d);
2538 4116 : if (pro) *pro = F.ro;
2539 4116 : M = F.M;
2540 4116 : if (best)
2541 : {
2542 4067 : if (!S->dT) S->dT = ZX_disc(S->T);
2543 4067 : d.expo_best_disc = expi(S->dT);
2544 : }
2545 :
2546 : /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
2547 4116 : y = cgetg(n*n + 1, t_VEC);
2548 4116 : b = cgetg(n*n + 1, t_COL);
2549 4116 : k = 1;
2550 4116 : if (!best) { gel(y,1) = pol_x(v); gel(b,1) = gen_0; k++; }
2551 27055 : for (i = 2; i <= n; i++)
2552 : {
2553 : GEN ch, ai;
2554 22939 : ai = gel(S->basis,i);
2555 22939 : ch = try_polmin(&d, S, gel(M,i), flag, &ai);
2556 22939 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2557 : }
2558 4116 : maxi = minss(n, 3);
2559 16051 : for (i = 1; i <= maxi; i++)
2560 68817 : for (j = i+1; j <= n; j++)
2561 : {
2562 : GEN ch, ai, v;
2563 56882 : ai = gadd(gel(S->basis,i), gel(S->basis,j));
2564 56882 : v = RgV_add(gel(M,i), gel(M,j));
2565 : /* defining polynomial for Q(w_i+w_j) */
2566 56881 : ch = try_polmin(&d, S, v, flag, &ai);
2567 56882 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2568 :
2569 56882 : ai = gsub(gel(S->basis,i), gel(S->basis,j));
2570 56882 : v = RgV_sub(gel(M,i), gel(M,j));
2571 : /* defining polynomial for Q(w_i-w_j) */
2572 56882 : ch = try_polmin(&d, S, v, flag, &ai);
2573 56882 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2574 : }
2575 4116 : setlg(y, k);
2576 4116 : setlg(b, k); filter(y, b, n);
2577 4116 : if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
2578 3472 : settyp(y, t_COL);
2579 3472 : (void)sort_factor_pol(mkmat2(y, b), cmpii);
2580 3472 : return mkmat2(b, y);
2581 : }
2582 :
2583 : /* FIXME: obsolete */
2584 : static GEN
2585 70 : Polred(GEN x, long flag, GEN fa)
2586 : {
2587 70 : pari_sp av = avma;
2588 : nfmaxord_t S;
2589 70 : if (fa)
2590 0 : nfinit_basic(&S, mkvec2(x,fa));
2591 : else
2592 70 : nfinit_basic_flag(&S, x, flag);
2593 63 : return gc_GEN(av, polred_aux(&S, NULL, flag));
2594 : }
2595 :
2596 : /* finds "best" polynomial in polred_aux list, defaulting to S->T if none of
2597 : * them is primitive. *px a ZX, characteristic polynomial of Mod(*pb,S->T),
2598 : * *pdx its discriminant if pdx != NULL. Set *pro = polroots(S->T) */
2599 : static void
2600 4088 : polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
2601 : {
2602 4088 : GEN y, dx, x = S->T; /* default value */
2603 : long i, l;
2604 4088 : y = polred_aux(S, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
2605 4081 : dx = S->dT;
2606 4081 : if (pb)
2607 : {
2608 3458 : GEN a, b = deg1pol_shallow(S->unscale, gen_0, varn(x));
2609 3458 : a = gel(y,1); l = lg(a);
2610 3458 : y = gel(y,2);
2611 6812 : for (i=1; i<l; i++)
2612 : {
2613 3354 : GEN yi = gel(y,i);
2614 3354 : if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
2615 : }
2616 3458 : *pb = b;
2617 : }
2618 : else
2619 : {
2620 623 : l = lg(y);
2621 1239 : for (i=1; i<l; i++)
2622 : {
2623 616 : GEN yi = gel(y,i);
2624 616 : if (ZX_is_better(yi,x,&dx)) x = yi;
2625 : }
2626 : }
2627 4081 : if (pdx) { if (!dx) dx = ZX_disc(x); *pdx = dx; }
2628 4081 : *px = x;
2629 4081 : }
2630 : static GEN
2631 3906 : polredbest_i(GEN T, long flag)
2632 : {
2633 : nfmaxord_t S;
2634 : GEN a;
2635 3906 : nfinit_basic_partial(&S, T);
2636 3906 : polredbest_aux(&S, NULL, &T, NULL, flag? &a: NULL);
2637 3899 : if (flag == 2)
2638 350 : T = mkvec2(T, a);
2639 3549 : else if (flag == 1)
2640 : {
2641 2926 : GEN b = (S.T0 == T)? pol_x(varn(T)): QXQ_reverse(a, S.T0);
2642 : /* charpoly(Mod(a,T0)) = T; charpoly(Mod(b,T)) = S.x */
2643 2926 : if (degpol(T) == 1) b = grem(b,T);
2644 2926 : T = mkvec2(T, mkpolmod(b,T));
2645 : }
2646 3899 : return T;
2647 : }
2648 : GEN
2649 3549 : polredbest(GEN T, long flag)
2650 : {
2651 3549 : pari_sp av = avma;
2652 3549 : if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
2653 3549 : return gc_GEN(av, polredbest_i(T, flag));
2654 : }
2655 : GEN
2656 21 : polredord(GEN x)
2657 : {
2658 21 : pari_sp av = avma;
2659 : GEN v, lt;
2660 : long i, n, vx;
2661 :
2662 21 : if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
2663 21 : x = Q_primpart(x); RgX_check_ZX(x,"polredord");
2664 21 : n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
2665 21 : if (n == 1) return gc_GEN(av, mkvec(x));
2666 14 : lt = leading_coeff(x); vx = varn(x);
2667 14 : if (is_pm1(lt))
2668 : {
2669 7 : if (signe(lt) < 0) x = ZX_neg(x);
2670 7 : v = pol_x_powers(n, vx);
2671 : }
2672 : else
2673 : { GEN L;
2674 : /* basis for Dedekind order */
2675 7 : v = cgetg(n+1, t_VEC);
2676 7 : gel(v,1) = scalarpol_shallow(lt, vx);
2677 14 : for (i = 2; i <= n; i++)
2678 7 : gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
2679 7 : gel(v,1) = pol_1(vx);
2680 7 : x = ZX_Q_normalize(x, &L);
2681 7 : v = gsubst(v, vx, monomial(ginv(L),1,vx));
2682 14 : for (i=2; i <= n; i++)
2683 7 : if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = pol_xn(i-1, vx);
2684 : }
2685 14 : return gc_upto(av, polred(mkvec2(x, v)));
2686 : }
2687 :
2688 : /* DEPRECATED: backward compatibility */
2689 : GEN
2690 56 : polred0(GEN x, long flag)
2691 : {
2692 56 : long fl = 0;
2693 56 : if (flag & 1) fl |= nf_PARTIALFACT;
2694 56 : if (flag & 2) fl |= nf_ORIG;
2695 56 : return Polred(x, fl, NULL);
2696 : }
2697 :
2698 : GEN
2699 14 : polred(GEN x) { return Polred(x, 0, NULL); }
2700 : GEN
2701 0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
2702 : GEN
2703 0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
2704 : GEN
2705 0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
2706 : GEN
2707 0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
2708 : GEN
2709 0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
2710 :
2711 : /********************************************************************/
2712 : /** **/
2713 : /** POLREDABS **/
2714 : /** **/
2715 : /********************************************************************/
2716 : /* set V[k] := matrix of multiplication by nk.zk[k] */
2717 : static GEN
2718 22708 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
2719 : {
2720 22708 : GEN v, Mk = cgetg(N+1, t_MAT);
2721 : long i, e;
2722 35374 : for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
2723 147930 : for ( ; i <=N; i++)
2724 : {
2725 125222 : v = vecmul(gel(M,k), gel(M,i));
2726 125222 : v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
2727 125223 : gel(Mk,i) = grndtoi(v, &e);
2728 125222 : if (e > -5) return NULL;
2729 : }
2730 22708 : gel(V,k) = Mk; return Mk;
2731 : }
2732 :
2733 : static GEN
2734 13083 : ZM_image_shallow(GEN M, long *pr)
2735 : {
2736 : long j, k, r;
2737 13083 : GEN y, d = ZM_pivots(M, &k);
2738 13083 : r = lg(M)-1 - k;
2739 13083 : y = cgetg(r+1,t_MAT);
2740 58944 : for (j=k=1; j<=r; k++)
2741 45861 : if (d[k]) gel(y,j++) = gel(M,k);
2742 13083 : *pr = r; return y;
2743 : }
2744 :
2745 : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
2746 : * Q from algo 2.7.6] */
2747 : static GEN
2748 14266 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
2749 : {
2750 14266 : CG_data *d = (CG_data*)chk->data;
2751 : GEN P, V, D, inv, bound, S, M;
2752 14266 : long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
2753 14266 : long i, j, prec, firstprim = 0, skipfirst = 0;
2754 : pari_sp av;
2755 :
2756 14266 : d->u = U;
2757 14266 : d->ZKembed = M = RgM_mul(d->M, U);
2758 :
2759 14266 : av = avma; bound = d->bound;
2760 14266 : D = cgetg(N+1, t_VECSMALL);
2761 95851 : for (i = 1; i <= N; i++)
2762 : {
2763 81592 : pari_sp av2 = avma;
2764 81592 : P = get_pol(d, gel(M,i));
2765 81592 : if (!P) pari_err_PREC("chk_gen_init");
2766 81585 : P = gc_GEN(av2, ZX_radical(P));
2767 81585 : D[i] = degpol(P);
2768 81585 : if (D[i] == N)
2769 : { /* primitive element */
2770 57609 : GEN B = embed_T2(gel(M,i), r1);
2771 57609 : if (!firstprim) firstprim = i; /* index of first primitive element */
2772 57609 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2773 57609 : if (gcmp(B,bound) < 0) bound = gc_leaf(av2, B);
2774 : }
2775 : else
2776 : {
2777 23976 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
2778 23976 : if (firstprim)
2779 : { /* cycle basis vectors so that primitive elements come last */
2780 2879 : GEN u = d->u, e = M;
2781 2879 : GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
2782 2879 : long tS = D[i];
2783 8875 : for (j = i; j > firstprim; j--)
2784 : {
2785 5996 : u[j] = u[j-1];
2786 5996 : e[j] = e[j-1];
2787 5996 : R[j] = R[j-1];
2788 5996 : D[j] = D[j-1];
2789 : }
2790 2879 : gel(u,firstprim) = tu;
2791 2879 : gel(e,firstprim) = te;
2792 2879 : gel(R,firstprim) = tR;
2793 2879 : D[firstprim] = tS; firstprim++;
2794 : }
2795 : }
2796 : }
2797 14259 : if (!firstprim)
2798 : { /* try (a little) to find primitive elements to improve bound */
2799 28 : GEN x = cgetg(N+1, t_VECSMALL);
2800 28 : if (DEBUGLEVEL>1)
2801 0 : err_printf("chk_gen_init: difficult field, trying random elements\n");
2802 308 : for (i = 0; i < 10; i++)
2803 : {
2804 : GEN e, B;
2805 3010 : for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
2806 280 : e = RgM_zc_mul(M, x);
2807 280 : B = embed_T2(e, r1);
2808 280 : if (gcmp(B,bound) >= 0) continue;
2809 16 : P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
2810 16 : if (!ZX_is_squarefree(P)) continue;
2811 16 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2812 16 : bound = B ;
2813 : }
2814 : }
2815 :
2816 14259 : if (firstprim != 1)
2817 : {
2818 13984 : inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
2819 13984 : V = gel(inv,1);
2820 54615 : for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
2821 : /* V corresponds to 1_Z */
2822 13984 : V = grndtoi(V, &j);
2823 13984 : if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
2824 13984 : S = mkmat(V); /* 1 */
2825 :
2826 13984 : V = cgetg(N+1, t_VEC);
2827 35864 : for (i = 1; i <= N; i++,skipfirst++)
2828 : { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
2829 : GEN Mx, M2;
2830 35864 : long j, k, h, rkM, dP = D[i];
2831 :
2832 35864 : if (dP == N) break; /* primitive */
2833 22708 : Mx = set_mulid(V, M, inv, r1, r2, N, i);
2834 22708 : if (!Mx) break; /* prec. problem. Stop */
2835 22708 : if (dP == 1) continue;
2836 9461 : rkM = lg(S)-1;
2837 9461 : M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
2838 9461 : gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
2839 9461 : k = 2;
2840 19072 : for (h = 1; h < dP; h++)
2841 : {
2842 : long r; /* add to M2 the elts of S * nf.zk[i] */
2843 39864 : for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
2844 13083 : setlg(M2, k); k = 1;
2845 13083 : S = ZM_image_shallow(shallowconcat(S,M2), &r);
2846 13911 : if (r == rkM) break;
2847 10439 : if (r > rkM)
2848 : {
2849 10439 : rkM = r;
2850 10439 : if (rkM == N) break;
2851 : }
2852 : }
2853 9461 : if (rkM == N) break;
2854 : /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
2855 : }
2856 : }
2857 : /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
2858 14259 : chk->skipfirst = skipfirst;
2859 14259 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
2860 :
2861 : /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
2862 14259 : bound = gc_leaf(av, bound);
2863 14259 : prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
2864 14259 : if (DEBUGLEVEL)
2865 0 : err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
2866 14259 : if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
2867 14259 : if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
2868 14259 : return bound;
2869 : }
2870 :
2871 : static GEN
2872 14273 : polredabs_i(GEN x, nfmaxord_t *S, GEN *u, long flag)
2873 : {
2874 14273 : FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
2875 : nffp_t F;
2876 : CG_data d;
2877 : GEN v, y, a;
2878 : long i, l;
2879 :
2880 14273 : nfinit_basic_flag(S, x, flag);
2881 14273 : x = S->T0;
2882 14273 : if (degpol(x) == 1)
2883 : {
2884 14 : long vx = varn(x);
2885 14 : *u = NULL;
2886 14 : return mkvec2(mkvec( pol_x(vx) ),
2887 14 : mkvec( deg1pol_shallow(gen_1, negi(gel(S->T,2)), vx) ));
2888 : }
2889 14259 : chk.data = (void*)&d;
2890 14259 : polred_init(S, &F, &d);
2891 14259 : d.bound = embed_T2(F.ro, d.r1);
2892 14259 : if (realprec(d.bound) > F.prec) d.bound = rtor(d.bound, F.prec);
2893 : for (;;)
2894 21 : {
2895 14280 : GEN R = R_from_QR(F.G, F.prec);
2896 14280 : if (R)
2897 : {
2898 14266 : d.prec = F.prec;
2899 14266 : d.M = F.M;
2900 14266 : v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
2901 14266 : if (v) break;
2902 : }
2903 21 : F.prec = precdbl(F.prec);
2904 21 : F.ro = NULL;
2905 21 : make_M_G(&F, 1);
2906 21 : if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",F.prec);
2907 : }
2908 14259 : y = gel(v,1);
2909 14259 : a = gel(v,2); l = lg(a);
2910 69972 : for (i = 1; i < l; i++) /* normalize wrt z -> -z */
2911 55713 : if (ZX_canon_neg(gel(y,i)) && (flag & (nf_ORIG|nf_RAW)))
2912 525 : gel(a,i) = ZC_neg(gel(a,i));
2913 14259 : *u = d.u; return v;
2914 : }
2915 :
2916 : GEN
2917 14000 : polredabs0(GEN x, long flag)
2918 : {
2919 14000 : pari_sp av = avma;
2920 : GEN Y, A, u, v;
2921 : nfmaxord_t S;
2922 : long i, l;
2923 :
2924 14000 : v = polredabs_i(x, &S, &u, flag);
2925 14000 : remove_duplicates(v);
2926 14000 : Y = gel(v,1);
2927 14000 : A = gel(v,2);
2928 14000 : l = lg(A); if (l == 1) pari_err_BUG("polredabs (missing vector)");
2929 14000 : if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
2930 14000 : if (!(flag & nf_ALL))
2931 : {
2932 14000 : GEN y = findmindisc(Y);
2933 28259 : for (i = 1; i < l; i++)
2934 28259 : if (ZX_equal(gel(Y,i), y)) break;
2935 14000 : Y = mkvec(gel(Y,i));
2936 14000 : A = mkvec(gel(A,i)); l = 2;
2937 : }
2938 14273 : if (flag & (nf_RAW|nf_ORIG)) for (i = 1; i < l; i++)
2939 : {
2940 273 : GEN y = gel(Y,i), a = gel(A,i);
2941 273 : if (u) a = RgV_RgC_mul(S.basis, ZM_ZC_mul(u, a));
2942 273 : if (flag & nf_ORIG)
2943 : {
2944 266 : a = QXQ_reverse(a, S.T);
2945 266 : if (!isint1(S.unscale)) a = gdiv(a, S.unscale); /* not RgX_Rg_div */
2946 266 : a = mkpolmod(a,y);
2947 : }
2948 273 : gel(Y,i) = mkvec2(y, a);
2949 : }
2950 14000 : return gc_GEN(av, (flag & nf_ALL)? Y: gel(Y,1));
2951 : }
2952 :
2953 : GEN
2954 0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
2955 : GEN
2956 13426 : polredabs(GEN x) { return polredabs0(x,0); }
2957 : GEN
2958 0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
2959 :
2960 : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
2961 : * flag & nf_ORIG: + element (base change)
2962 : * flag & nf_ABSOLUTE: absolute polynomial */
2963 : static GEN
2964 637 : rnfpolred_i(GEN nf, GEN R, long flag, long best)
2965 : {
2966 637 : const char *f = best? "rnfpolredbest": "rnfpolredabs";
2967 637 : const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
2968 637 : GEN listP = NULL, red, pol, A, P, T, rnfeq;
2969 637 : pari_sp av = avma;
2970 637 : long even = 0;
2971 :
2972 637 : if (typ(R) == t_VEC) {
2973 14 : if (lg(R) != 3) pari_err_TYPE(f,R);
2974 14 : listP = gel(R,2);
2975 14 : R = gel(R,1);
2976 : }
2977 637 : if (typ(R) != t_POL) pari_err_TYPE(f,R);
2978 637 : nf = checknf(nf);
2979 637 : T = nf_get_pol(nf);
2980 637 : R = RgX_nffix(f, T, R, 0);
2981 637 : if (best || (flag & nf_PARTIALFACT))
2982 : {
2983 364 : rnfeq = abs? nf_rnfeq(nf, R): nf_rnfeqsimple(nf, R);
2984 364 : pol = gel(rnfeq,1);
2985 364 : if (listP) pol = mkvec2(pol, listP);
2986 357 : red = best? polredbest_i(pol, abs? 1: 2)
2987 364 : : polredabs0(pol, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
2988 364 : P = gel(red,1);
2989 364 : A = gel(red,2);
2990 : }
2991 : else
2992 : {
2993 : nfmaxord_t S;
2994 : GEN rnf, u, v, y, a;
2995 : long i, j, l;
2996 : pari_timer ti;
2997 273 : if (DEBUGLEVEL>1) timer_start(&ti);
2998 273 : rnf = rnfinit(nf, R);
2999 273 : rnfeq = rnf_get_map(rnf);
3000 273 : pol = rnf_zkabs(rnf);
3001 273 : if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
3002 273 : v = polredabs_i(pol, &S, &u, nf_ORIG);
3003 273 : pol = gel(pol,1);
3004 273 : y = gel(v,1); P = findmindisc(y);
3005 273 : a = gel(v,2);
3006 273 : l = lg(y); A = cgetg(l, t_VEC);
3007 1897 : for (i = j = 1; i < l; i++)
3008 1624 : if (ZX_equal(gel(y,i),P)) gel(A,j++) = gel(a,i);
3009 273 : setlg(A,j); /* mod(A[i], pol) are all roots of P in Q[X]/(pol) */
3010 273 : even = j > 2 && !odd(degpol(P)) && !odd(ZX_deflate_order(P));
3011 1365 : for (i = 1; i < j; i++)
3012 : {
3013 1092 : GEN t = gel(A,i);
3014 1092 : if (u) t = ZM_ZC_mul(u,t);
3015 1092 : gel(A,i) = RgV_RgC_mul(S.basis, t);
3016 : }
3017 : }
3018 637 : if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
3019 637 : if (flag & nf_ABSOLUTE)
3020 : {
3021 14 : if (flag & nf_ORIG)
3022 : {
3023 7 : GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
3024 7 : GEN k = gel(rnfeq,3); /* Mod(variable(R),R) + k*a root of pol */
3025 7 : if (typ(A) == t_VEC) A = gel(A,1); /* any root will do */
3026 7 : a = RgX_RgXQ_eval(a, lift_shallow(A), P); /* Mod(a, P) root of T */
3027 7 : P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
3028 : }
3029 14 : return gc_GEN(av, P);
3030 : }
3031 623 : if (typ(A) != t_VEC)
3032 : {
3033 357 : A = eltabstorel_lift(rnfeq, A);
3034 357 : P = lift_if_rational( RgXQ_charpoly(A, R, varn(R)) );
3035 : }
3036 : else
3037 : { /* canonical factor */
3038 266 : long i, l = lg(A), v = varn(R);
3039 266 : GEN besta = NULL;
3040 1316 : for (i = 1; i < l; i++)
3041 : {
3042 1050 : GEN a = eltabstorel_lift(rnfeq, gel(A,i));
3043 1050 : GEN p = lift_if_rational( RgXQ_charpoly(a, R, v) );
3044 1050 : if (even)
3045 : {
3046 532 : GEN q = RgX_unscale(p, gen_m1);
3047 532 : if (cmp_universal(q, p) < 0) { p = q; a = gneg(a); }
3048 : }
3049 1050 : if (i == 1 || cmp_universal(p, P) < 0) { P = p; besta = a; }
3050 : }
3051 266 : A = besta;
3052 : }
3053 623 : if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,R),P));
3054 623 : return gc_GEN(av, P);
3055 : }
3056 : GEN
3057 280 : rnfpolredabs(GEN nf, GEN R, long flag)
3058 280 : { return rnfpolred_i(nf,R,flag, 0); }
3059 : GEN
3060 357 : rnfpolredbest(GEN nf, GEN R, long flag)
3061 : {
3062 357 : if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
3063 357 : return rnfpolred_i(nf,R,flag, 1);
3064 : }
|
